<jupyter_start><jupyter_text>Notebook 4 - g-computation method (teacher notebook) Let's reconsider the LaLonde dataset as we try to implement a g-computation implementation from scratch.<jupyter_code>import matplotlib.pyplot as plt
import pandas as pd
from sklearn.ensemble import RandomForestRegressor
from sklearn.metrics import mean_squared_error
from sklearn.model_selection import train_test_split
from sklearn.linear_model import LinearRegression
from sklearn.neighbors import NearestNeighbors
from sklearn.preprocessing import StandardScaler
from statsmodels.distributions.empirical_distribution import ECDF

%matplotlib inline
plt.rcParams['figure.figsize'] = [10, 5]
# Windows users: replace the following line with `df = pd.read_csv(r"..\..\data\lalonde.csv")`
df = pd.read_csv("../../data/lalonde.csv")<jupyter_output><empty_output><jupyter_text>To remind you, the LaLonde dataset should contain the following 12 columns:age   $\;\;\;\;\;\;$age in years.educ    $\;\;\;\;\;\;$years of schooling.black    $\;\;\;\;\;\;$indicator variable for blacks.hisp    $\;\;\;\;\;\;$indicator variable for Hispanics.married    $\;\;\;\;\;\;$indicator variable for martial status.nodegr    $\;\;\;\;\;\;$indicator variable for high school diploma.re74    $\;\;\;\;\;\;$real earnings in 1974.re75    $\;\;\;\;\;\;$real earnings in 1975.re78    $\;\;\;\;\;\;$real earnings in 1978.u74    $\;\;\;\;\;\;$indicator variable for earnings in 1974 being zero.u75    $\;\;\;\;\;\;$indicator variable for earnings in 1975 being zero.treat    $\;\;\;\;\;\;$an indicator variable for treatment status.<jupyter_code>df.head()
df['ID'] = range(0,445)
df = df[['ID', 'age', 'educ', 'black', 'hisp', 'married', 're78', 'treat']]
# Let's scale our continuous covariates
features = df[['age', 'educ']]

# Use scaler of choice; here Standard scaler is used
scaler = StandardScaler().fit(features.values)
features = scaler.transform(features.values)

df[['scaled_age', 'scaled_educ']] = features<jupyter_output><empty_output><jupyter_text>First, let's train a predictive model with scikit-learn. This model should take in the treatment variable and covariates, and predict the outcome variable (i.e. `re78`) If you're curious about how well the model performs that's probably a good thing to be curious about. Let's split our small dataset into test and training sets. Now, admittedly, we have a tiny dataset, but we're doing this to show that standard ML evaluation methods are still relevant here.<jupyter_code>train_df, test_df = train_test_split(df, test_size = 0.20, random_state = 512)
features = ['scaled_age', 'scaled_educ', 'black', 'hisp', 'married', 'treat']<jupyter_output><empty_output><jupyter_text>Given the tiny data, let's use simple ML like a multiple linear regression model<jupyter_code>temp_model = LinearRegression().fit(X = train_df[features], y = train_df['re78'])
mean_squared_error(temp_model.predict(test_df[features]).reshape(-1,1), test_df['re78'].values.reshape(-1,1))<jupyter_output><empty_output><jupyter_text>You could compare different model types here to see which performs best in a hold-out set. Not surprisingly, given the small dataset, it's going to be hard to find a model that performs well. Ideally, we want a fairly good predicive model so that we can see accurate counterfactual outcomes. Anyways, this is just a demonstration so the model doesn't have to perform outstandingly well. Let's now train on the whole dataset.<jupyter_code>model = LinearRegression().fit(X = df[features], y = df['re78'])<jupyter_output><empty_output><jupyter_text>Now that we have a model we can start making predictions to see counterfactual outcomes     EXERCISE: Take the full dataset and "force" everyone to have received the treatment, what is the mean value of real earnings from 1978 that you predict with your model?<jupyter_code>df_treat = df.copy()
df_treat['treat'] = 1
df_treat['treat'].value_counts()
mean_treat = model.predict(df_treat[features]).mean()
print(f"Mean outcome value is {round(mean_treat, 2)}")<jupyter_output><empty_output><jupyter_text>EXERCISE: Now "force" everyone to have not received the treatment. What is the mean outcome value you get now?<jupyter_code>df_ctrl = df.copy()
df_ctrl['treat'] = 0
df_ctrl['treat'].value_counts()
mean_ctrl = model.predict(df_ctrl[features]).mean()
print(f"Mean outcome value is {round(mean_ctrl, 2)}")<jupyter_output><empty_output><jupyter_text>The difference (6283.43 - 4601.56 = 1681.87) between these two represents the average treat effect (ATE). This is the expected mean difference in real earnings between a population that received the treatment and the same population that didn't receive the treatment. Confounders are adjusted for.   You could obtain confidence intervals around this estimate by taking bootstrap replicates of the original dataset, training a new model on each replicate, and then calculating the ATE for each replicate. The 2.5th and 97.5th percentiles of these values would represent the bounds of your 95% confidence interval.     EXERCISE: Can you calculate the average treatment effect among the treated (ATT)? Hint: just as the metric's name suggests, you would focus exclusively on the originally treated group, and see how a receiving and not receiving the treatment among that group would play out.<jupyter_code># We're only focusing on the originally treated group now
df_treat = df[df['treat'] == 1].copy()
mean_treat = model.predict(df_treat[features]).mean()
# We still only focus on the originally treated group, but we predict on the counterfactual / unobserved event for them 
# (a world where they didn't receive the treatment)
df_ctrl = df[df['treat'] == 1].copy()
df_ctrl['treat'] = 0
mean_ctrl = model.predict(df_ctrl[features]).mean()
print(f"The ATT is {mean_treat - mean_ctrl}")<jupyter_output>The ATT is 1681.8695520414522
{'fraction_non_alphanumeric': 0.11451557685724068, 'fraction_numerical': 0.019342690859294763, 'mean_word_length': 5.506681514476615, 'pattern_counts': {'":': 0, '<': 32, '<?xml version=': 0, '>': 32, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '8', 'dir': 'jupyter-structured-clean-dedup', 'id': '342961', 'license': 'MIT', 'n_tokens_mistral': 1849, 'n_tokens_neox': 1652, 'n_words': 633, 'path': 'notebooks/teacher/4_teacher_g_comp.ipynb', 'repo_name': 'ronikobrosly/causal_inference_intro', 'size': '5842'}