<jupyter_start><jupyter_text>Prediction   Introduction In labor economics an important question is what determines the wage of workers. This is a causal question, but we could begin to investigate from a predictive perspective.In the following wage example,  𝑌  is the hourly wage of a worker and  𝑋  is a vector of worker's characteristics, e.g., education, experience, gender. Two main questions here are:* How to use job-relevant characteristics, such as education and experience, to best predict wages?* What is the difference in predicted wages between men and women with the same job-relevant characteristics?In this lab, we focus on the prediction question first. DataThe data set we consider is from the March Supplement of the U.S. Current Population Survey, year 2015. We select white non-hispanic individuals, aged 25 to 64 years, and working more than 35 hours per week during at least 50 weeks of the year. We exclude self-employed workers; individuals living in group quarters; individuals in the military, agricultural or private household sectors; individuals with inconsistent reports on earnings and employment status; individuals with allocated or missing information in any of the variables used in the analysis; and individuals with hourly wage below  3 .The variable of interest  𝑌  is the hourly wage rate constructed as the ratio of the annual earnings to the total number of hours worked, which is constructed in turn as the product of number of weeks worked and the usual number of hours worked per week. In our analysis, we also focus on single (never married) workers. The final sample is of size  $n=5150$ . ---<jupyter_code>import numpy as np
import pandas as pd
import statsmodels.api as sm
import sklearn.linear_model as lm
import sklearn.metrics as metrics
import statsmodels.formula.api as smf
from sklearn.model_selection import cross_val_score, train_test_split<jupyter_output><empty_output><jupyter_text>Data AnalysisSet the following file_directory to a place where you downloaded https://github.com/CausalAIBook/MetricsMLNotebooks/blob/main/PM1/wage2015_subsample_inference.rdata<jupyter_code>file_dir = None
df = pd.read_csv(file_dir)
df.head()
df = df.drop("Unnamed: 0")<jupyter_output><empty_output><jupyter_text>Construct variablesWe are constructing the output variable  $Y$  and the matrix  $Z$  which includes the characteristics of workers that are given in the data.<jupyter_code>Y = np.log(df['wage'])
Z = df[df.columns.difference(['wage', 'lwage'])]
Z.shape<jupyter_output><empty_output><jupyter_text>For the outcome variable wage and a subset of the raw regressors, we calculate the empirical mean and other empirical measures to get familiar with the data.<jupyter_code>Z.describe()<jupyter_output><empty_output><jupyter_text>E.g., the share of female workers in our sample is ~44% ( 𝑠𝑒𝑥=1  if female).<jupyter_code># if you want to print this table to latex
Z.describe().to_latex()<jupyter_output><empty_output><jupyter_text>---  Prediction Question Now, we will construct a prediction rule for hourly wage  $Y$, which depends linearly on job-relevant characteristics  $X$:$$𝑌=\beta′𝑋+𝜖.$$Our goals are:* Predict wages using various characteristics of workers.* Assess the predictive performance using the (adjusted) sample MSE, the (adjusted) sample $R^2$  and the out-of-sample MSE and  $R^2$ .We employ two different specifications for prediction:1. **Basic Model**:  $𝑋$  consists of a set of raw regressors (e.g. gender, experience, education indicators, occupation and industry indicators, regional indicators).2. **Flexible Model**:  $𝑋$  consists of all raw regressors from the basic model plus occupation and industry indicators, transformations (e.g.,  $exp^2$  and  $exp^3$ ) and additional two-way interactions of polynomial in experience with other regressors. An example of a regressor created through a two-way interaction is experience times the indicator of having a college degree.Using the **Flexible Model**, enables us to approximate the real relationship by a more complex regression model and therefore to reduce the bias. The **Flexible Model** increases the range of potential shapes of the estimated regression function. In general, flexible models often deliver good prediction accuracy but give models which are harder to interpret.Now, let us fit both models to our data by running ordinary least squares (ols):<jupyter_code># 1. Basic Model
model1 = 'lwage ~ sex + exp1 + shs + hsg+ scl + clg + mw + so + we + C(occ2) + C(ind2)'
results1 = smf.ols(model1, data=df).fit()
print(results1.summary())
# 2. Flexible Model
model2 = 'lwage ~ sex + shs+hsg+scl+clg+C(occ2)+C(ind2)+mw+so+we + (exp1+exp2+exp3+exp4)*(shs+hsg+scl+clg+C(occ2)+C(ind2)+mw+so+we)'
results2 = smf.ols(model2, data=df).fit()
print(results2.summary())<jupyter_output><empty_output><jupyter_text>Next we try the LassoWe use the statsmodels package with the formula api and the sklearn cross-validation package to hypertune the parameter.<jupyter_code># Lasso with cross-validation
model3 = smf.ols(model2, data=df)
X = model3.data.exog
y = model3.data.endog

# train model using CV 
results3 = lm.LassoCV(cv=10, normalize=True, fit_intercept=False,
                      alphas=np.logspace(-6,2, 100)).fit(X,y)
lasso_df = pd.DataFrame()
lasso_df['Var'] = model3.exog_names
lasso_df['Coef.'] = results3.coef_
lasso_df<jupyter_output><empty_output><jupyter_text>Result evaluationNow, we can evaluate the performance of both models based on the (adjusted)  $R^2_{𝑠𝑎𝑚𝑝𝑙𝑒}$  and the (adjusted) $𝑀𝑆𝐸_{𝑠𝑎𝑚𝑝𝑙𝑒}$.<jupyter_code># Print the R^2

# basic model
r2_basic = results1.rsquared
print(f"R-squared for the basic model: {r2_basic}"))

adj_r2_basic = results1.rsquared_adj
print(f"Adjusted R-squared for the basic model: {adj_r2_basic}")

# flexible model
r2_flexible = results2.rsquared
print(f"R-squared for the flexible model: {r2_flexible}")

adj_r2_flexible = results2.rsquared_adj
print(f"Adjusted R-squared for the flexible model: {adj_r2_flexible}")

# Lasso model
lasso_preds = results3.predict(X)
mse_lasso = sum((y-lasso_preds)**2)/y.shape[0]
r2_lasso = 1. - mse_lasso/np.var(y)

print(f"R-squared for the lasso model: {r2_lasso}"))

n = X.shape[0]
p = X.shape[1]
adj_r2_lasso = 1. - (1.-r2_lasso)*((n-1.)/(n-p-1.))
print(f"Adjusted R-squared for the lasso model: {adj_r2_lasso}")
# Print the MSE

# basic model
mse_basic = np.mean(results1.resid**2)
print(f"MSE for the basic model: {mse_basic}")

adj_mse_basic = results1.mse_resid
print(f"Adjusted MSE for the basic model: {adj_mse_basic}")

# flexible model
mse_flexible = np.mean(results2.resid**2)
print(f"MSE for the flexible model: {mse_flexible}")

adj_mse_flexible = results2.mse_resid
print(f"Adjusted MSE for the flexible model: {adj_mse_flexible}")

# Lasso model
#mse_lasso = metrics.mean_squared_error(lasso_preds, y)
print(f"MSE for the lasso model: {mse_lasso}")

adj_mse_lasso = mse_lasso*n/(n-p)
print(f"Adjusted MSE for the lasso model: {adj_mse_lasso}")
# store the results in a table
res_df = pd.DataFrame()

res_df['Model'] = ['Basic reg', 'Flexible reg', 'Flexible Lasso']

res_df['p'] = [results1.params.shape[0],
           results2.params.shape[0],
           results2.params.shape[0]]

res_df['R2'] = [r2_basic, r2_flexible, r2_lasso]
res_df['MSE'] = [mse_basic, mse_flexible, mse_lasso]

res_df['adj_R2'] = [adj_r2_basic, adj_r2_flexible, adj_r2_lasso]
res_df['adj_MSE'] = [adj_mse_basic, adj_mse_flexible, adj_mse_lasso]

# Show results
res_df.head()
# print to Latex
res_df.to_latex()<jupyter_output><empty_output><jupyter_text>Data Splitting Measure the prediction quality of the two models via data splitting:* Randomly split the data into one training sample and one testing sample. Here we just use a simple method (stratified splitting is a more sophiscticated version of splitting that we can consider).* Use the training sample for estimating the parameters of the Basic Model and the Flexible Model.* Use the testing sample for evaluation. Predict the  𝚠𝚊𝚐𝚎  of every observation in the testing sample based on the estimated parameters in the training sample.* Calculate the Mean Squared Prediction Error $MSE_{test}$ based on the testing sample for both prediction models.<jupyter_code># split the data into training and testing sets

# get the indices
n = df.shape[0]
sh = 4./5.
train_idx = np.random.choice(n, int(np.floor(sh*n)), replace=False)
test_idx = np.setdiff1d(np.arange(n), train_idx)
# Basic model

# estimating the parameters in the training sample
mod1 = smf.ols(model1, data=df)
X1 = mod1.data.exog
y1 = mod1.data.endog

# separate training and testing sets
X1_train = X1[train_idx,:]
X1_test = X1[test_idx,:]
y1_train = y1[train_idx]
y1_test = y1[test_idx]

# estimating the parameters in the training sample
regbasic = sm.OLS(y1_train, X1_train).fit()

# predict out of sample
trainregbasic = regbasic.predict(X1_test)

# calculating out-of-sample MSE
MSE_test1 = sum((y1_test-trainregbasic)**2)/y1_test.shape[0]
R2_test1 = 1. - MSE_test1/np.var(y1_test)

print("Test MSE for the basic model: "+ str(MSE_test1))
print("Test R2 for the basic model: "+ str(R2_test1))<jupyter_output><empty_output><jupyter_text>In the basic model, the  $𝑀𝑆𝐸_{𝑡𝑒𝑠𝑡}$  is quite closed to the  $𝑀𝑆𝐸_{𝑠𝑎𝑚𝑝𝑙𝑒}$.<jupyter_code># Flexible model

# estimating the parameters in the training sample
mod2 = smf.ols(model2, data=df)
X2 = mod2.data.exog
y2 = mod2.data.endog

# separate training and testing sets
X2_train = X2[train_idx,:]
X2_test = X2[test_idx,:]
y2_train = y2[train_idx]
y2_test = y2[test_idx]

# estimating the parameters in the training sample
regflex = sm.OLS(y2_train, X2_train).fit()

# predict out of sample
trainregflex = regflex.predict(X2_test)

# calculating out-of-sample MSE
MSE_test2 = sum((y2_test-trainregflex)**2)/y2_test.shape[0]
R2_test2 = 1. - MSE_test2/np.var(y2_test)

print("Test MSE for the basic model: "+ str(MSE_test2))
print("Test R2 for the basic model: "+ str(R2_test2))<jupyter_output><empty_output><jupyter_text>In the flexible model, the discrepancy between the  $𝑀𝑆𝐸_{𝑡𝑒𝑠𝑡}$  and the  $𝑀𝑆𝐸_{𝑠𝑎𝑚𝑝𝑙𝑒}$  is not large.It is worth to notice that the  $𝑀𝑆𝐸_{𝑡𝑒𝑠𝑡}$  vary across different data splits. Hence, it is a good idea average the out-of-sample MSE over different data splits to get valid results.Nevertheless, we observe that, based on the out-of-sample  $𝑀𝑆𝐸$ , the basic model using ols regression performs is about as well (or slightly better) than the flexible model.Next, let us use lasso regression in the flexible model instead of ols regression. Lasso (least absolute shrinkage and selection operator) is a penalized regression method that can be used to reduce the complexity of a regression model when the number of regressors  $p$  is relatively large in relation to  $n$ .Note that the out-of-sample  𝑀𝑆𝐸  on the test sample can be computed for any other black-box prediction method as well. Thus, let us finally compare the performance of lasso regression in the flexible model to ols regression.<jupyter_code># Flexible model using Lasso

# train model using CV 
reglasso = lm.LassoCV(cv=10, fit_intercept=False,
                      alphas=np.logspace(-6,2, 100)).fit(X2_train,y2_train)

# predict out of sample
trainreglasso = reglasso.predict(X2_test)

# calculating out-of-sample MSE
MSE_test3 = sum((y2_test-trainreglasso)**2)/y2_test.shape[0]
R2_test3 = 1. - MSE_test3/np.var(y2_test)

print("Test MSE for the basic model: "+ str(MSE_test3))
print("Test R2 for the basic model: "+ str(R2_test3))
# store the results in a table
res_df2 = pd.DataFrame()

res_df2['Model'] = ['Basic reg', 'Flexible reg', 'Flexible Lasso']

res_df2['$MSE_{test}$'] = [MSE_test1, MSE_test2, MSE_test3]
res_df2['$R^2_{test}$'] = [R2_test1, R2_test2, R2_test3]

# Show results
res_df2.head()
# print to Latex
res_df2.to_latex()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.0964527027027027, 'fraction_numerical': 0.01875, 'mean_word_length': 5.372981700753498, 'pattern_counts': {'":': 0, '<': 45, '<?xml version=': 0, '>': 45, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '11', 'dir': 'jupyter-structured-clean-dedup', 'id': '361544', 'license': 'MIT', 'n_tokens_mistral': 4250, 'n_tokens_neox': 4115, 'n_words': 1343, 'path': 'PM1/PM1_prediction.ipynb', 'repo_name': 'CausalAIBook/MetricsMLNotebooks', 'size': '11840'}