# Estimating the Effect of Schooling on Wages: Intrumental Variables Application

### Summary of Contents:
1. [Introduction](#intro)
2. [NLSYM Dataset](#data)
3. [A Gentle Start: The Naive Approach](#naive)
4. [Using Instrumental Variables: 2SLS](#2sls)
5. [Bonus: Deep Instrumental Variables](#deepiv)

**Important:** This notebook provides some room for experimentation in solving this problem. If you are looking for the end-to-end solution, look for the same file name with the "Complete" suffix added.

# 1. Introduction <a class="anchor" id="intro"></a>

To measure true causal effects of a treatment $T$ on an outcome $Y$ from observational data, we need to record all features $X$ that might influence both $T$ and $Y$. These $X$'s are called confounders. 

When some confounders are not recorded in the data, we might get biased estimates of the treatment effect. Here is an example:
* Children of high-income parents might attain higher levels of education (e.g. college) since they can afford it
* Children of high-income parents might also obtain better paying jobs due to parents' connections and knowledge
* At first sight, it might appear as if education has an effect on income, when in fact this could be fully explained by family background

There are several reasons for not recording all possible confounders, such as incomplete data or a confounder that is difficult to quantify (e.g. parental involvement). However, not all is lost! In cases such as these, we can use instrumental variables $Z$, features that affect the outcome only through their effect on the treatment. 

In this notebook, we use a real-world problem to show how treatment effects can be extracted with the help of instrumental variables. 

# 2. NLSYM Dataset <a class="anchor" id="data"></a>

The **causal impact of schooling on wages** had been studied at length. Although it is generally agreed that there is a positive impact, it is difficult to measure this effect precisely. The core problem is that education levels are not assigned at random in the population and we cannot record all possible confounders. (Think about the value parents assign to education. How would you quantify how valuable parents think their children's education is?). 



To get around this issue, we can use **proximity to a 4-year college** as an instrumental variable. Having a college nearby can allow individuals (especially low-income ones) to complete more years of education. Hence, if there was a positive treatment effect, we would expect these individuals to have higher wages on average. Note that college proximity is a valid IV since it does not affect wages directly.  

We use data from the National Longitudinal Survey of Young Men (NLSYM, 1966) to estimate the average treatment effect (ATE) of education on wages (see also [Card, 1999](https://www.nber.org/papers/w4483)). The NLSYM data contains entries from men ages 14-24 that were interviewed in 1966 and again in 1976. 

The dataset contains the following variables:
* $Y$ (outcome): wages (log)
* $T$ (treatment): years of schooling
* $Z$ (IV): proximity to a 4-year college (binary)
* $X$ (heterogeneity): e.g. parental education
* $W$ (controls): e.g. family composition, location, etc.

The world can then be modelled as:
$$
\begin{align}
Y & = \theta(X) \cdot T + f(W) + \epsilon\\
T & = g(Z, W) + \eta
\end{align}
$$
where $\epsilon, \eta$ are uncorrelated error terms.


```python
# Python imports
import keras
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from scipy.stats import pearsonr
from sklearn.ensemble import RandomForestRegressor
from sklearn.linear_model import LinearRegression
from sklearn.model_selection import KFold
from sklearn.preprocessing import PolynomialFeatures, StandardScaler

# EconML imports
from econml.dml import DMLCateEstimator
from econml.two_stage_least_squares import NonparametricTwoStageLeastSquares
from econml.deepiv import DeepIVEstimator

%matplotlib inline
```

    Using TensorFlow backend.



```python
# Data processing
df = pd.read_csv("data/card.csv", dtype=float)
# Filter out individuals with low education levels (outliers)
data_filter = df['educ'].values >= 6
# Define some variables
T = df['educ'].values[data_filter]
Z = df['nearc4'].values[data_filter]
Y = df['lwage'].values[data_filter]

# Impute missing values with mean, add dummy columns
# Filter outliers (interviewees with less than 6 years of education)
X_df = df[['exper', 'expersq']].copy()
X_df['fatheduc'] = df['fatheduc'].fillna(value=df['fatheduc'].mean())
X_df['fatheduc_nan'] = df['fatheduc'].isnull() * 1
X_df['motheduc'] = df['motheduc'].fillna(value=df['motheduc'].mean())
X_df['motheduc_nan'] = df['motheduc'].isnull() * 1
X_df[['momdad14', 'sinmom14', 'reg661', 'reg662',
        'reg663', 'reg664', 'reg665', 'reg666', 'reg667', 'reg668', 'reg669', 'south66']] = df[['momdad14', 'sinmom14', 
        'reg661', 'reg662','reg663', 'reg664', 'reg665', 'reg666', 'reg667', 'reg668', 'reg669', 'south66']]
X_df[['black', 'smsa', 'south', 'smsa66']] = df[['black', 'smsa', 'south', 'smsa66']]
columns_to_scale = ['fatheduc', 'motheduc', 'exper', 'expersq']
# Scale continuous variables
scaler = StandardScaler()
X_df[columns_to_scale] = scaler.fit_transform(X_df[columns_to_scale])
X = X_df.values[data_filter]
```


```python
# Explore data
X_df.head()
```




<div>
<style scoped>
    .dataframe tbody tr th:only-of-type {
        vertical-align: middle;
    }

    .dataframe tbody tr th {
        vertical-align: top;
    }

    .dataframe thead th {
        text-align: right;
    }
</style>
<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>exper</th>
      <th>expersq</th>
      <th>fatheduc</th>
      <th>fatheduc_nan</th>
      <th>motheduc</th>
      <th>motheduc_nan</th>
      <th>momdad14</th>
      <th>sinmom14</th>
      <th>reg661</th>
      <th>reg662</th>
      <th>...</th>
      <th>reg665</th>
      <th>reg666</th>
      <th>reg667</th>
      <th>reg668</th>
      <th>reg669</th>
      <th>south66</th>
      <th>black</th>
      <th>smsa</th>
      <th>south</th>
      <th>smsa66</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>0</th>
      <td>1.725159</td>
      <td>1.896133</td>
      <td>5.439188e-16</td>
      <td>1</td>
      <td>0.000000</td>
      <td>1</td>
      <td>1.0</td>
      <td>0.0</td>
      <td>1.0</td>
      <td>0.0</td>
      <td>...</td>
      <td>0.0</td>
      <td>0.0</td>
      <td>0.0</td>
      <td>0.0</td>
      <td>0.0</td>
      <td>0.0</td>
      <td>1.0</td>
      <td>1.0</td>
      <td>0.0</td>
      <td>1.0</td>
    </tr>
    <tr>
      <th>1</th>
      <td>0.034739</td>
      <td>-0.172321</td>
      <td>-6.134540e-01</td>
      <td>0</td>
      <td>-0.786159</td>
      <td>0</td>
      <td>1.0</td>
      <td>0.0</td>
      <td>1.0</td>
      <td>0.0</td>
      <td>...</td>
      <td>0.0</td>
      <td>0.0</td>
      <td>0.0</td>
      <td>0.0</td>
      <td>0.0</td>
      <td>0.0</td>
      <td>0.0</td>
      <td>1.0</td>
      <td>0.0</td>
      <td>1.0</td>
    </tr>
    <tr>
      <th>2</th>
      <td>1.725159</td>
      <td>1.896133</td>
      <td>1.223740e+00</td>
      <td>0</td>
      <td>0.553046</td>
      <td>0</td>
      <td>1.0</td>
      <td>0.0</td>
      <td>1.0</td>
      <td>0.0</td>
      <td>...</td>
      <td>0.0</td>
      <td>0.0</td>
      <td>0.0</td>
      <td>0.0</td>
      <td>0.0</td>
      <td>0.0</td>
      <td>0.0</td>
      <td>1.0</td>
      <td>0.0</td>
      <td>1.0</td>
    </tr>
    <tr>
      <th>3</th>
      <td>0.276228</td>
      <td>0.052254</td>
      <td>3.051432e-01</td>
      <td>0</td>
      <td>0.553046</td>
      <td>0</td>
      <td>1.0</td>
      <td>0.0</td>
      <td>0.0</td>
      <td>1.0</td>
      <td>...</td>
      <td>0.0</td>
      <td>0.0</td>
      <td>0.0</td>
      <td>0.0</td>
      <td>0.0</td>
      <td>0.0</td>
      <td>0.0</td>
      <td>1.0</td>
      <td>0.0</td>
      <td>1.0</td>
    </tr>
    <tr>
      <th>4</th>
      <td>1.725159</td>
      <td>1.896133</td>
      <td>-6.134540e-01</td>
      <td>0</td>
      <td>-1.120960</td>
      <td>0</td>
      <td>1.0</td>
      <td>0.0</td>
      <td>0.0</td>
      <td>1.0</td>
      <td>...</td>
      <td>0.0</td>
      <td>0.0</td>
      <td>0.0</td>
      <td>0.0</td>
      <td>0.0</td>
      <td>0.0</td>
      <td>0.0</td>
      <td>1.0</td>
      <td>0.0</td>
      <td>1.0</td>
    </tr>
  </tbody>
</table>
<p>5 rows × 22 columns</p>
</div>



# 3. A Gentle Start: The Naive Approach <a class="anchor" id="naive"></a>

Let's assume we know nothing about instrumental variables and we want to measure the treatment effect of schooling on wages. We can apply an IV-free method like Double Machine Learning (DML) to do this and extract a treatment effect. 


```python
### DIY Section!
# Define a DML Estimator
dml_est = None
```


```python
# Fit the DML estimator and calculate the heterogeneous treatment effect
dml_te = None
# Uncoment after you calculated dml_te
"""
dml_ate = dml_te.mean()
"""
```




    '\ndml_ate = dml_te.mean()\n'




```python
# Uncomment after you calculated the dml_ate
"""
print("Average treatment effect: {0:.3f}".format(dml_ate))
"""
```




    '\nprint("Average treatment effect: {0:.3f}".format(dml_ate))\n'



This treatment effect is smaller than other values obtained in literature via IV. Why could that be? 

Because DML (like all IV-free methods) assumes that the residual errors are uncorrelated (i.e. $Y - \hat{Y}$ is uncorrelated with $T - \hat{T}$). Let's test this assumption:


```python
### DIY Section
# Calculate the residuals Y_res = Y - Y_pred and T_res = T - T_pred
# Bonus: Split data in 2 parts, train on one part and predict on the other to avoid over-fitting 
T_res, Y_res = np.zeros(T.shape[0]), np.zeros(Y.shape[0])
```


```python
# Uncomment this after you calculated T_res and Y_res
"""
plt.scatter(T_res, Y_res)
plt.show()
"""
```




    '\nplt.scatter(T_res, Y_res)\nplt.show()\n'




```python
# Uncomment this after you calculated T_res and Y_res
"""
corr_coeficient = pearsonr(T_res, Y_res)[0]
print("Correlation coefficient between T and Y errors: {0:.2f}".format(corr_coeficient))
"""
```




    '\ncorr_coeficient = pearsonr(T_res, Y_res)[0]\nprint("Correlation coefficient between T and Y errors: {0:.2f}".format(corr_coeficient))\n'



The correlation coefficient between the residuals is quite large, which means that there is some unobserved variables that affect both $T$ and $Y$. To get an accurate estimate in this case, we need to use IVs. 

# 4. Using Intrumental Variables: 2SLS <a class="anchor" id="2sls"></a>

Two stage least square regression procedure (2SLS):
1. Fit a model $T \sim W, Z$
2. Fit a linear model $Y \sim \hat{T}$ where $\hat{T}$ is the prediction of the model in step 1.
The coefficient from 2. above is the average treatment effect.

If interested in heterogeneous treatment effects, fit a model $Y \sim \hat{T}\otimes h(X)$, where $h(X)$ is a chosen featurization of the treatment effect. 

For more information, see the `econml` [documentation](https://econml.azurewebsites.net).


```python
# For average treatment effects, X is a column of 1s
W = X
Z = Z.reshape(-1, 1)
T = T.reshape(-1, 1)
X_ate = np.ones_like(Z)
```


```python
# We apply 2SLS from the EconML library
two_sls_est = NonparametricTwoStageLeastSquares(
    t_featurizer=PolynomialFeatures(degree=1, include_bias=False),
    x_featurizer=PolynomialFeatures(degree=1, include_bias=False),
    z_featurizer=PolynomialFeatures(degree=1, include_bias=False),
    dt_featurizer=None) # dt_featurizer only matters for marginal_effect
```


```python
two_sls_est.fit(Y, T, X_ate, W, Z)
two_sls_ate = two_sls_est.effect(np.ones((1,1)))[0]
```


```python
print("Average treatment effect: {0:.3f}".format(two_sls_ate))
```

    Average treatment effect: 0.134


# 5. Bonus: Deep Instrumental Variables <a class="anchor" id="deepiv"></a>

For very flexible, but fully non-parametric IV methods, you can use neural networks for the two models in 2SLS and a mixture of gaussians for the featurizer $h(X)$. In `econml`, this method is called DeepIV. 

The NLSYM dataset is small (on neural net scale) so applying DeepIV is a bit of a stretch. Nevertheless, we apply DeepIV the NLSYM data as an example. You should not read too much into the results. 


```python
# Define treatment model, T ~ X, Z
treatment_model = keras.Sequential([keras.layers.Dense(64, activation='relu', input_shape=(X.shape[1] + 1,)),
                                    keras.layers.Dropout(rate=0.17),
                                    keras.layers.Dense(32, activation='relu'),
                                    keras.layers.Dropout(rate=0.17),
                                    keras.layers.Dense(1)])
# Define outcome model, Y ~ T_hat, X
response_model = keras.Sequential([keras.layers.Dense(64, activation='relu', input_shape=(X.shape[1] + 1,)),
                                   keras.layers.Dropout(rate=0.17),
                                   keras.layers.Dense(32, activation='relu'),
                                   keras.layers.Dropout(rate=0.17),
                                   keras.layers.Dense(1)])
```


```python
keras_fit_options = { "epochs": 30,
                      "validation_split": 0.3,
                      "callbacks": [keras.callbacks.EarlyStopping(patience=2, restore_best_weights=True)]}

deepIvEst = DeepIVEstimator(n_components = 10, # number of gaussians in our mixture density network
                            m = lambda z, x : treatment_model(keras.layers.concatenate([z,x])), # treatment model
                            h = lambda t, x : response_model(keras.layers.concatenate([t,x])),  # response model
                            n_samples = 1, # number of samples to use to estimate the response
                            use_upper_bound_loss = False, # whether to use an approximation to the true loss
                            n_gradient_samples = 1, # number of samples to use in second estimate of the response (to make loss estimate unbiased)
                            optimizer='adam', # Keras optimizer to use for training - see https://keras.io/optimizers/ 
                            first_stage_options = keras_fit_options, # options for training treatment model
                            second_stage_options = keras_fit_options) # options for training response model
```


```python
deepIvEst.fit(Y, T, X, Z)
deepIv_effect = deepIvEst.effect(X)
```

    WARNING:tensorflow:From C:\Users\moprescu\AppData\Local\Continuum\anaconda3\lib\site-packages\tensorflow\python\ops\math_ops.py:3066: to_int32 (from tensorflow.python.ops.math_ops) is deprecated and will be removed in a future version.
    Instructions for updating:
    Use tf.cast instead.
    Train on 2093 samples, validate on 898 samples
    Epoch 1/30
    2093/2093 [==============================] - 1s 654us/step - loss: 23.8329 - val_loss: 5.9042
    Epoch 2/30
    2093/2093 [==============================] - 0s 95us/step - loss: 5.0409 - val_loss: 3.9612
    Epoch 3/30
    2093/2093 [==============================] - 0s 118us/step - loss: 3.8898 - val_loss: 3.4648
    Epoch 4/30
    2093/2093 [==============================] - 0s 116us/step - loss: 3.5212 - val_loss: 3.2288
    Epoch 5/30
    2093/2093 [==============================] - 0s 118us/step - loss: 3.3116 - val_loss: 3.0616
    Epoch 6/30
    2093/2093 [==============================] - 0s 77us/step - loss: 3.1305 - val_loss: 2.9131
    Epoch 7/30
    2093/2093 [==============================] - 0s 118us/step - loss: 2.9753 - val_loss: 2.7573
    Epoch 8/30
    2093/2093 [==============================] - 0s 116us/step - loss: 2.8158 - val_loss: 2.5689
    Epoch 9/30
    2093/2093 [==============================] - 0s 111us/step - loss: 2.5800 - val_loss: 2.2856
    Epoch 10/30
    2093/2093 [==============================] - 0s 112us/step - loss: 2.2167 - val_loss: 1.7570
    Epoch 11/30
    2093/2093 [==============================] - 0s 114us/step - loss: 1.7293 - val_loss: 1.4546
    Epoch 12/30
    2093/2093 [==============================] - 0s 110us/step - loss: 1.6448 - val_loss: 1.4081
    Epoch 13/30
    2093/2093 [==============================] - 0s 113us/step - loss: 1.6269 - val_loss: 1.4076
    Epoch 14/30
    2093/2093 [==============================] - 0s 117us/step - loss: 1.5911 - val_loss: 1.3996
    Epoch 15/30
    2093/2093 [==============================] - 0s 87us/step - loss: 1.5620 - val_loss: 1.3765
    Epoch 16/30
    2093/2093 [==============================] - 0s 119us/step - loss: 1.5518 - val_loss: 1.3627
    Epoch 17/30
    2093/2093 [==============================] - 0s 110us/step - loss: 1.5245 - val_loss: 1.3463
    Epoch 18/30
    2093/2093 [==============================] - 0s 110us/step - loss: 1.4822 - val_loss: 1.3022
    Epoch 19/30
    2093/2093 [==============================] - 0s 113us/step - loss: 1.5112 - val_loss: 1.3304
    Epoch 20/30
    2093/2093 [==============================] - 0s 108us/step - loss: 1.5124 - val_loss: 1.3415
    Train on 2093 samples, validate on 898 samples
    Epoch 1/30
    2093/2093 [==============================] - 3s 1ms/step - loss: 9.3704 - val_loss: 1.8564
    Epoch 2/30
    2093/2093 [==============================] - 0s 166us/step - loss: 5.7598 - val_loss: 8.5731
    Epoch 3/30
    2093/2093 [==============================] - 0s 164us/step - loss: 5.2696 - val_loss: 1.6304
    Epoch 4/30
    2093/2093 [==============================] - 0s 167us/step - loss: 34.8642 - val_loss: 1.8287
    Epoch 5/30
    2093/2093 [==============================] - 0s 168us/step - loss: 5.5516 - val_loss: 564.5439



```python
print("Average treatment effect: {0:.3f}".format(deepIv_effect.mean()))
```

    Average treatment effect: 0.191



```python
# Heterogeneity of treatment effects
plt.hist(deepIv_effect)
plt.show()
```


```python

```

c4d17cba36e1b5698250600f6f5eb9728a23fee3
ipynb
Jupyter Notebook
notebooks/Hands-on Part II -- Estimating the Effect of Schooling on Wages.ipynb
ashlili0401/econmic
e2137e20af0f90382c2fcf2db53823b8737124d0
MIT
notebooks/Hands-on Part II -- Estimating the Effect of Schooling on Wages.ipynb
ashlili0401/econmic
e2137e20af0f90382c2fcf2db53823b8737124d0
MIT
notebooks/Hands-on Part II -- Estimating the Effect of Schooling on Wages.ipynb
ashlili0401/econmic
e2137e20af0f90382c2fcf2db53823b8737124d0
MIT
Qwen/Qwen-72B
1. YES
2. YES