<jupyter_start><jupyter_text>Causal Inference of Kang-Schafer simulation.       Run in Google Colab        View source on GitHub   >[Causal Inference of Kang-Schafer simulation.](scrollTo=RYsS3JrpHaJI)>>[Imports](scrollTo=VGDMBD7TRIf4)>>[Correctly Specified Model](scrollTo=X1D4OBiEZhnT)>>[Misspecified Model](scrollTo=0zuCvTGVdMVj)>>[Benchmark Execution Time](scrollTo=ITCnSqDiyjIx) We illustrate empirical calibration to estimate the average treatment effect on the treated (ATT) on Kang-Schafer simulation under both correctly specified and misspecified models, and benchmark the execution time.  For details of simulation setup, please refer to [kang_schafer_population_mean.ipynb](./kang_schafer_population_mean.ipynb).  Imports<jupyter_code>from matplotlib import pyplot as plt
import numpy as np
import pandas as pd
import patsy
import seaborn as sns
import timeit

# install and import ec
!pip install -q git+https://github.com/google/empirical_calibration
import empirical_calibration as ec

sns.set_style('whitegrid')
%config InlineBackend.figure_format='retina'<jupyter_output><empty_output><jupyter_text>Correctly Specified Model We run the simulation 1000 times under correctly specified logistic propensity score and linear outcome regression. For each simulation we estimate the average treatment effect on the treated (ATT) using empirical calibration.<jupyter_code>def estimate_att(formula):
  simulation = ec.data.kang_schafer.Simulation(size=1000)

  t = simulation.treatment
  y = simulation.outcome

  df = pd.DataFrame(
      np.column_stack(
          [simulation.covariates, simulation.transformed_covariates]))
  df.columns = ["z1", "z2", "z3", "z4", "x1", "x2", "x3", "x4"]
  x = patsy.dmatrix(formula, df, return_type="dataframe").values


  weights = ec.maybe_exact_calibrate(covariates=x[t == 0],
                                     target_covariates=x[t == 1])[0]

  return (np.mean(y[t == 1]) - np.mean(y[t == 0]),
          np.mean(y[t == 1]) - np.sum(y[t == 0] * weights))


def show_estimates(estimates, col='weighted'):
  ax = estimates[col].hist(bins=20, alpha=0.8, edgecolor='none')
  plt.axvline(estimates[col].mean(), linestyle='dashed', color='red')
  print('bias of {} is {}'.format(col, estimates[col].mean()))
  print('rmse of {} is {}'.format(col, np.sqrt(np.mean((estimates[col] - 0.) ** 2))))
estimates = pd.DataFrame(
    [estimate_att("-1 + z1 + z2 + z3 + z4") for i in xrange(1000)])
estimates.columns = ['raw', 'weighted']<jupyter_output><empty_output><jupyter_text>The mean of the 1000 ATT estimates after weight correction is very close to the true zero ATT.<jupyter_code>show_estimates(estimates,'raw')
show_estimates(estimates,'weighted')<jupyter_output>bias of weighted is 0.00108065962091
rmse of weighted is 0.0854238555451<jupyter_text>Misspecified Model If the transformed covariates are observed in place of the true covariates, both the propensity score model and outcome regression model become misspecified. We run 1000 simulations and for each simulation estimate the ATT by balancing the transformed covariates. The causal estimate is no longer unbiased.<jupyter_code>estimates_miss = pd.DataFrame([estimate_att("-1 + x1 + x2 + x3 + x4") for i in xrange(1000)])
estimates_miss.columns = ['raw', 'weighted']
show_estimates(estimates_miss)<jupyter_output>bias of weighted is -4.37849137373
rmse of weighted is 4.50105698077<jupyter_text>One reasonable strategy is to expand the set of balancing covariates and hope it will make the model less "misspecified". If we additional balance the two-way interactions and the log transformation, the bias indeed reduces.<jupyter_code>formula = ("-1 + (x1 + x2 + x3 + x4)**2 + I(np.log(x1)) + I(np.log(x2)) + "
           "I(np.log(x3)) + I(np.log(x4))")

estimates_expanded = pd.DataFrame([estimate_att(formula) for i in xrange(1000)])
estimates_expanded.columns = ['raw', 'weighted']
show_estimates(estimates_expanded)<jupyter_output>bias of weighted is -1.88003551222
rmse of weighted is 2.07831395567<jupyter_text>If the model was misspecified in the sense that more covariates are included than necessary, the causal estimate remains unbiased.<jupyter_code>formula = "-1 + z1 + z2 + z3 + z4 + x1 + x2 + x3 + x4"
estimates_redundant = pd.DataFrame([estimate_att(formula) for i in range(1000)])
estimates_redundant.columns = ['raw', 'weighted']
show_estimates(estimates_redundant)<jupyter_output>bias of weighted is -0.00229483538445
rmse of weighted is 0.0891054715229<jupyter_text>Benchmark Execution Time The execution time is generally linear with respect to the sample size. With 1 million control units, it takes around 1 second to find the weights.<jupyter_code>np.random.seed(123)
simulation = ec.data.kang_schafer.Simulation(size=2000)
x1 = simulation.covariates[simulation.treatment == 1]
x0 = simulation.covariates[simulation.treatment == 0]
pd.Series(timeit.repeat(
    'ec.maybe_exact_calibrate(x0, x1)',
    setup='from __main__ import x1, x0, ec',
    repeat=100,
    number=1)).describe()
np.random.seed(123)
simulation = ec.data.kang_schafer.Simulation(size=20000)
x1 = simulation.covariates[simulation.treatment == 1]
x0 = simulation.covariates[simulation.treatment == 0]
pd.Series(timeit.repeat(
    'ec.maybe_exact_calibrate(x0, x1)',
    setup='from __main__ import x1, x0, ec',
    repeat=100,
    number=1)).describe()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.11268129416139828, 'fraction_numerical': 0.04369654146522871, 'mean_word_length': 5.543795620437956, 'pattern_counts': {'":': 0, '<': 25, '<?xml version=': 0, '>': 34, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '7', 'dir': 'jupyter-structured-clean-dedup', 'id': '608990', 'license': 'Apache-2.0', 'n_tokens_mistral': 1991, 'n_tokens_neox': 1755, 'n_words': 480, 'path': 'notebooks/causal_inference_kang_schafer.ipynb', 'repo_name': 'miaojingang/empirical_calibration', 'size': '5378'}