Treatment effects under conditional independence¶
Author: Josef Perktold
This notebook illustrates the basic usage of the new treatment effect functionality in statsmodels.
The main class is statsmodels.treatment.treatment_effects.TreatmentEffect
.
This class estimates treatment effect and potential outcome using 5 different methods, ipw, ra, aipw, aipw-wls, ipw-ra. The last three methods require both a treatment or selection model and an outcome model. Standard errors and inference are based on the joint GMM representation of selection or treatment model, outcome model and effect functions. The approach for inference follows Stata, however Stata support a wider range of models. Estimation and inference are valid under conditional independence or ignorability.
The outcome model is currently limited to a linear model based on OLS. Treatment is currently restricted to binary treatment which can be either Logit or Probit.
The example follows Cattaneo.
[1]:
import os
import numpy as np
from numpy.testing import assert_allclose
import pandas as pd
from statsmodels.regression.linear_model import OLS
from statsmodels.discrete.discrete_model import Probit
from statsmodels.treatment.treatment_effects import (
TreatmentEffect
)
from statsmodels.treatment.tests.results import results_teffects as res_st
# Load data for example
cur_dir = os.path.abspath(os.path.dirname(res_st.__file__))
file_name = 'cataneo2.csv'
file_path = os.path.join(cur_dir, file_name)
dta_cat = pd.read_csv(file_path)
methods = ['ra', 'ipw', 'aipw', 'aipw_wls', 'ipw_ra']
methods_st = [
("ra", res_st.results_ra),
("ipw", res_st.results_ipw),
("aipw", res_st.results_aipw),
("aipw_wls", res_st.results_aipw_wls),
("ipw_ra", res_st.results_ipwra),
]
# allow wider display of data frames
pd.set_option('display.width', 500)
[2]:
dta_cat.head()
[2]:
bweight | mmarried | mhisp | fhisp | foreign | alcohol | deadkids | mage | medu | fage | ... | prenatal | birthmonth | lbweight | fbaby | prenatal1 | mbsmoke_ | mmarried_ | fbaby_ | prenatal1_ | mage2 | |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
0 | 3459 | married | 0 | 0 | 0 | 0 | 0 | 24 | 14 | 28 | ... | 1 | 12 | 0 | No | Yes | 0 | 1 | 0 | 1 | 576.0 |
1 | 3260 | notmarried | 0 | 0 | 1 | 0 | 0 | 20 | 10 | 0 | ... | 1 | 7 | 0 | No | Yes | 0 | 0 | 0 | 1 | 400.0 |
2 | 3572 | married | 0 | 0 | 1 | 0 | 0 | 22 | 9 | 30 | ... | 1 | 3 | 0 | No | Yes | 0 | 1 | 0 | 1 | 484.0 |
3 | 2948 | married | 0 | 0 | 0 | 0 | 0 | 26 | 12 | 30 | ... | 1 | 1 | 0 | No | Yes | 0 | 1 | 0 | 1 | 676.0 |
4 | 2410 | married | 0 | 0 | 0 | 0 | 0 | 20 | 12 | 21 | ... | 1 | 3 | 1 | Yes | Yes | 0 | 1 | 1 | 1 | 400.0 |
5 rows × 28 columns
Create TreatmentEffect instance and compute ipw¶
The TreatmentEffect class requires - a OLS model instance for the outcome model, - a results instance of the selection model and - a treatment indicator variable.
In the following example we use Probit as the selection model. Using Logit is also supported.
[3]:
# treatment selection model
formula = 'mbsmoke_ ~ mmarried_ + mage + mage2 + fbaby_ + medu'
res_probit = Probit.from_formula(formula, dta_cat).fit()
# outcome model
formula_outcome = 'bweight ~ prenatal1_ + mmarried_ + mage + fbaby_'
mod = OLS.from_formula(formula_outcome, dta_cat)
# treatment indicator variable
tind = np.asarray(dta_cat['mbsmoke_'])
teff = TreatmentEffect(mod, tind, results_select=res_probit)
Optimization terminated successfully.
Current function value: 0.439575
Iterations 6
After creating the TreatmentEffect instance, we can call any of the 5 methods to compute potential outcomes, POM0, POM1, and average treatment effect, ATE. POM0 is the potential outcome for the no treatment group, POM1 is the potential outcome for the treatment group, treatment effect is POM1 - POM0.
For example teff.ipw()
computes POM and ATE using inverse probability weighting. The probability of treatment is also commonly called the propensity score. The summary
of the estimation includes standard errors and confidence interval for POM and ATE.
Standard errors and other inferential statistics are based on the Generalized Method of Moments (GMM) representation of the selection and outcome models and the moment conditions for the results statistic. Method ipw
uses the selection model but not the outcome model. Method ra
uses the outcome model but not the selection model. The doubly robust estimators aipw
, aipw-wls
, ipw-ra
include both selection and outcome models, where at least one of those two has to be correctly
specified to get consistent estimates of the treatment effect. The moment conditions for the target variables, POM0, POM1, and ATE are based on POM0 and ATE. The remaining POM1 is computed as a linear combination of POM0 and ATE.
The internal gmm results are attached to the treatment results as results_gmm
.
By default the treatment effect methods computes average treatment effect, where average is take over the sample observations. Option effect_group
can be used to compute either average treatment effect on the treated, ATT, using effect_group=1
or average treatment effect on the non-treated using effect_group=0
.
[4]:
res = teff.ipw()
res
[4]:
<class 'statsmodels.treatment.treatment_effects.TreatmentEffectResults'>
Test for Constraints
==============================================================================
coef std err z P>|z| [0.025 0.975]
------------------------------------------------------------------------------
ATE -230.6891 25.817 -8.936 0.000 -281.289 -180.089
POM0 3403.4632 9.571 355.586 0.000 3384.704 3422.223
POM1 3172.7741 24.001 132.193 0.000 3125.733 3219.815
==============================================================================
[5]:
res.summary_frame()
[5]:
coef | std err | z | P>|z| | Conf. Int. Low | Conf. Int. Upp. | |
---|---|---|---|---|---|---|
ATE | -230.689070 | 25.816758 | -8.935633 | 4.048542e-19 | -281.288985 | -180.089154 |
POM0 | 3403.463163 | 9.571412 | 355.586324 | 0.000000e+00 | 3384.703540 | 3422.222785 |
POM1 | 3172.774093 | 24.001059 | 132.193085 | 0.000000e+00 | 3125.732881 | 3219.815305 |
[6]:
print(res.results_gmm.summary())
_IPWGMM Results
==============================================================================
Dep. Variable: y Hansen J: 3.988e-09
Model: _IPWGMM Prob (Hansen J): nan
Method: GMM
Date: Sat, 12 Oct 2024
Time: 17:05:40
No. Observations: 4642
==============================================================================
coef std err z P>|z| [0.025 0.975]
------------------------------------------------------------------------------
p 0 -230.6891 25.817 -8.936 0.000 -281.289 -180.089
p 1 3403.4632 9.571 355.586 0.000 3384.704 3422.223
p 2 -1.5583 0.461 -3.380 0.001 -2.462 -0.655
p 3 -0.6485 0.055 -11.711 0.000 -0.757 -0.540
p 4 0.1744 0.036 4.836 0.000 0.104 0.245
p 5 -0.0033 0.001 -4.921 0.000 -0.005 -0.002
p 6 -0.2176 0.050 -4.390 0.000 -0.315 -0.120
p 7 -0.0864 0.010 -8.630 0.000 -0.106 -0.067
==============================================================================
average treatment effect on the treated
see more below
[7]:
teff.ipw(effect_group=1)
[7]:
<class 'statsmodels.treatment.treatment_effects.TreatmentEffectResults'>
Test for Constraints
==============================================================================
coef std err z P>|z| [0.025 0.975]
------------------------------------------------------------------------------
ATE -225.1796 23.658 -9.518 0.000 -271.549 -178.811
POM0 3362.8393 14.198 236.855 0.000 3335.012 3390.667
POM1 3137.6597 19.071 164.526 0.000 3100.281 3175.038
==============================================================================
average treatment effect on the untreated
[8]:
teff.ipw(effect_group=0)
[8]:
<class 'statsmodels.treatment.treatment_effects.TreatmentEffectResults'>
Test for Constraints
==============================================================================
coef std err z P>|z| [0.025 0.975]
------------------------------------------------------------------------------
ATE -231.8782 27.699 -8.371 0.000 -286.168 -177.588
POM0 3412.9116 9.283 367.634 0.000 3394.716 3431.107
POM1 3181.0334 26.120 121.786 0.000 3129.840 3232.227
==============================================================================
Other methods to compute ATE work in the same or similar way as for ipw
for example regression adjustment ra
and double robust ipw_ra
.
[9]:
res_ra = teff.ra()
res_ra
[9]:
<class 'statsmodels.treatment.treatment_effects.TreatmentEffectResults'>
Test for Constraints
==============================================================================
coef std err z P>|z| [0.025 0.975]
------------------------------------------------------------------------------
ATE -239.6392 23.824 -10.059 0.000 -286.333 -192.945
POM0 3403.2423 9.525 357.288 0.000 3384.573 3421.911
POM1 3163.6031 21.864 144.698 0.000 3120.751 3206.455
==============================================================================
[10]:
res_ra.summary_frame()
[10]:
coef | std err | z | P>|z| | Conf. Int. Low | Conf. Int. Upp. | |
---|---|---|---|---|---|---|
ATE | -239.639211 | 23.824021 | -10.058722 | 8.408247e-24 | -286.333435 | -192.944988 |
POM0 | 3403.242272 | 9.525207 | 357.288006 | 0.000000e+00 | 3384.573209 | 3421.911335 |
POM1 | 3163.603060 | 21.863509 | 144.697867 | 0.000000e+00 | 3120.751371 | 3206.454750 |
[11]:
ra2 = teff.ipw_ra(effect_group=1, return_results=True)
ra2.summary_frame()
[11]:
coef | std err | z | P>|z| | Conf. Int. Low | Conf. Int. Upp. | |
---|---|---|---|---|---|---|
ATE | -223.545262 | 23.794008 | -9.395023 | 5.720507e-21 | -270.180660 | -176.909864 |
POM0 | 3361.204984 | 14.465009 | 232.367989 | 0.000000e+00 | 3332.854088 | 3389.555880 |
POM1 | 3137.659722 | 19.070923 | 164.525844 | 0.000000e+00 | 3100.281400 | 3175.038045 |
All methods in TreatmentEffect¶
The following computes and prints ATE and POM for all methods. (We include the call to TreatmentEffect as a reminder.)
[12]:
teff = TreatmentEffect(mod, tind, results_select=res_probit)
for m in methods:
res = getattr(teff, m)()
print("\n", m)
print(res.summary_frame())
ra
coef std err z P>|z| Conf. Int. Low Conf. Int. Upp.
ATE -239.639211 23.824021 -10.058722 8.408247e-24 -286.333435 -192.944988
POM0 3403.242272 9.525207 357.288006 0.000000e+00 3384.573209 3421.911335
POM1 3163.603060 21.863509 144.697867 0.000000e+00 3120.751371 3206.454750
ipw
coef std err z P>|z| Conf. Int. Low Conf. Int. Upp.
ATE -230.689070 25.816758 -8.935633 4.048542e-19 -281.288985 -180.089154
POM0 3403.463163 9.571412 355.586324 0.000000e+00 3384.703540 3422.222785
POM1 3172.774093 24.001059 132.193085 0.000000e+00 3125.732881 3219.815305
aipw
coef std err z P>|z| Conf. Int. Low Conf. Int. Upp.
ATE -230.989648 26.214445 -8.811541 1.234375e-18 -282.369017 -179.610280
POM0 3403.355674 9.568514 355.682783 0.000000e+00 3384.601731 3422.109616
POM1 3172.366025 24.427402 129.869153 0.000000e+00 3124.489197 3220.242854
aipw_wls
coef std err z P>|z| Conf. Int. Low Conf. Int. Upp.
ATE -227.195618 27.372036 -8.300282 1.038645e-16 -280.843822 -173.547414
POM0 3403.250651 9.596571 354.631943 0.000000e+00 3384.441717 3422.059585
POM1 3176.055033 25.654642 123.800406 0.000000e+00 3125.772859 3226.337206
ipw_ra
coef std err z P>|z| Conf. Int. Low Conf. Int. Upp.
ATE -229.967078 26.629411 -8.635830 5.830196e-18 -282.159765 -177.774391
POM0 3403.335639 9.571288 355.577620 0.000000e+00 3384.576260 3422.095018
POM1 3173.368561 24.871955 127.588224 0.000000e+00 3124.620425 3222.116697
Results in Stata¶
The results in statsmodels are very close to the results in Stata because both packages use the same approach.
[13]:
for m, st in methods_st:
print("\n", m)
res = pd.DataFrame(st.table[:2, :6], index = ["ATE", "POM0"], columns=st.table_colnames[:6])
print(res)
ra
b se z pvalue ll ul
ATE -239.639211 23.824021 -10.058722 8.408247e-24 -286.333435 -192.944988
POM0 3403.242272 9.525207 357.288005 0.000000e+00 3384.573209 3421.911335
ipw
b se z pvalue ll ul
ATE -230.688638 25.815244 -8.936140 4.030006e-19 -281.285586 -180.091690
POM0 3403.462709 9.571369 355.587873 0.000000e+00 3384.703170 3422.222247
aipw
b se z pvalue ll ul
ATE -230.989201 26.210565 -8.812828 1.220276e-18 -282.360964 -179.617438
POM0 3403.355253 9.568472 355.684297 0.000000e+00 3384.601393 3422.109114
aipw_wls
b se z pvalue ll ul
ATE -227.195618 27.347936 -8.307597 9.765984e-17 -280.796587 -173.594649
POM0 3403.250651 9.596622 354.630065 0.000000e+00 3384.441618 3422.059684
ipw_ra
b se z pvalue ll ul
ATE -229.967078 26.626676 -8.636718 5.785117e-18 -282.154403 -177.779752
POM0 3403.335639 9.571260 355.578657 0.000000e+00 3384.576315 3422.094963
Treatment effects without inference¶
It is possible to compute POM and ATE without computing standard errors and inferential statistics. In this case the GMM model is not computed.
[14]:
for m in methods:
print("\n", m)
res = getattr(teff, m)(return_results=False)
print(res)
ra
(np.float64(-239.6392114643395), np.float64(3403.242271935487), np.float64(3163.6030604711477))
ipw
(np.float64(-230.6886377952617), np.float64(3403.4627086845567), np.float64(3172.7740708892948))
aipw
(np.float64(-230.98920111257803), np.float64(3403.3552531738355), np.float64(3172.3660520612575))
aipw_wls
(np.float64(-227.19561818674902), np.float64(3403.2506509757864), np.float64(3176.0550327890373))
ipw_ra
(np.float64(-229.96707793513224), np.float64(3403.3356393074205), np.float64(3173.3685613722882))
Treatment effect on the treated¶
Treatment effects on subgroups are not available for aipw
and aipw-wls
.
effect_group
choses the group for which treatement effect and potential outcomes are computed Options are “all” for sample average treatment effect, 1
for average treatment effect on the treated and 0
for average treatment effect on the untreated.
Note: The row labels in the pandas dataframe, POM and ATE, are the same even for treatment effect on subgroups.
[15]:
for m in methods:
if m.startswith("aipw"):
continue
res = getattr(teff, m)(effect_group=1)
print("\n", m)
print(res.summary_frame())
ra
coef std err z P>|z| Conf. Int. Low Conf. Int. Upp.
ATE -223.301651 22.742195 -9.818826 9.342545e-23 -267.875534 -178.727767
POM0 3360.961373 12.757489 263.450069 0.000000e+00 3335.957154 3385.965592
POM1 3137.659722 19.070923 164.525844 0.000000e+00 3100.281400 3175.038045
ipw
coef std err z P>|z| Conf. Int. Low Conf. Int. Upp.
ATE -225.179608 23.658119 -9.518069 1.764269e-21 -271.548669 -178.810546
POM0 3362.839334 14.197866 236.855264 0.000000e+00 3335.012028 3390.666640
POM1 3137.659726 19.070923 164.525845 0.000000e+00 3100.281404 3175.038049
ipw_ra
coef std err z P>|z| Conf. Int. Low Conf. Int. Upp.
ATE -223.545262 23.794008 -9.395023 5.720507e-21 -270.180660 -176.909864
POM0 3361.204984 14.465009 232.367989 0.000000e+00 3332.854088 3389.555880
POM1 3137.659722 19.070923 164.525844 0.000000e+00 3100.281400 3175.038045
Treatment effect on the untreated¶
Similar to ATT, we can compute average treatment effect on the untreated by using option effect_group=0
.
[16]:
for m in methods:
if m.startswith("aipw"):
# not available
continue
res = getattr(teff, m)(effect_group=0)
print("\n", m)
print(res.summary_frame())
ra
coef std err z P>|z| Conf. Int. Low Conf. Int. Upp.
ATE -243.375488 24.902030 -9.773319 1.465697e-22 -292.182569 -194.568406
POM0 3412.911593 9.283454 367.633804 0.000000e+00 3394.716358 3431.106829
POM1 3169.536106 23.128805 137.038471 0.000000e+00 3124.204480 3214.867731
ipw
coef std err z P>|z| Conf. Int. Low Conf. Int. Upp.
ATE -231.878176 27.699436 -8.371224 5.702294e-17 -286.168073 -177.588279
POM0 3412.911593 9.283454 367.633804 0.000000e+00 3394.716357 3431.106829
POM1 3181.033418 26.119760 121.786472 0.000000e+00 3129.839629 3232.227206
ipw_ra
coef std err z P>|z| Conf. Int. Low Conf. Int. Upp.
ATE -231.125972 28.813022 -8.021580 1.043933e-15 -287.598458 -174.653487
POM0 3412.911593 9.283454 367.633804 0.000000e+00 3394.716358 3431.106829
POM1 3181.785621 27.301318 116.543297 0.000000e+00 3128.276021 3235.295221
The docstring for the TreatmentEffect class and it’s methods can be obtained using help
help(teff)