条件独立下的处理效果¶
作者：Josef Perktold
此笔记本说明了 statsmodels 中新的处理效果功能的基本用法。
主要类是 statsmodels.treatment.treatment_effects.TreatmentEffect。
此类使用 5 种不同的方法估计处理效果和潜在结果，ipw、ra、aipw、aipw-wls、ipw-ra。最后三种方法需要处理或选择模型和结果模型。标准误差和推断基于选择或处理模型、结果模型和效应函数的联合 GMM 表示。推断方法遵循 Stata，但 Stata 支持更广泛的模型。估计和推断在条件独立或可忽略性下有效。
结果模型目前仅限于基于 OLS 的线性模型。处理目前仅限于二元处理，可以是 Logit 或 Probit。
该示例遵循 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 | 否 | 是 | 0 | 1 | 0 | 1 | 576.0 |
| 1 | 3260 | 未婚 | 0 | 0 | 1 | 0 | 0 | 20 | 10 | 0 | ... | 1 | 7 | 0 | 否 | 是 | 0 | 0 | 0 | 1 | 400.0 |
| 2 | 3572 | married | 0 | 0 | 1 | 0 | 0 | 22 | 9 | 30 | ... | 1 | 3 | 0 | 否 | 是 | 0 | 1 | 0 | 1 | 484.0 |
| 3 | 2948 | married | 0 | 0 | 0 | 0 | 0 | 26 | 12 | 30 | ... | 1 | 1 | 0 | 否 | 是 | 0 | 1 | 0 | 1 | 676.0 |
| 4 | 2410 | married | 0 | 0 | 0 | 0 | 0 | 20 | 12 | 21 | ... | 1 | 3 | 1 | 是 | 是 | 0 | 1 | 1 | 1 | 400.0 |
5 行 × 28 列
创建 TreatmentEffect 实例并计算 ipw¶
TreatmentEffect 类需要 - 结果模型的 OLS 模型实例，- 选择模型的结果实例，- 处理指示变量。
在以下示例中，我们使用 Probit 作为选择模型。也支持使用 Logit。
[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
创建 TreatmentEffect 实例后，我们可以调用 5 种方法中的任何一种来计算潜在结果、POM0、POM1 和平均处理效果 ATE。POM0 是未接受处理组的潜在结果，POM1 是接受处理组的潜在结果，处理效果是 POM1 - POM0。
例如 teff.ipw() 使用逆概率加权计算 POM 和 ATE。处理概率通常也称为倾向得分。估计的 summary 包括 POM 和 ATE 的标准误差和置信区间。
标准误差和其他推断统计量基于选择和结果模型的广义矩估计 (GMM) 表示以及结果统计量的矩条件。方法 ipw 使用选择模型但不使用结果模型。方法 ra 使用结果模型但不使用选择模型。双重稳健估计量 aipw、aipw-wls、ipw-ra 包括选择模型和结果模型，其中这两个模型中至少有一个必须被正确指定才能得到处理效果的一致估计。目标变量 POM0、POM1 和 ATE 的矩条件基于 POM0 和 ATE。剩余的 POM1 计算为 POM0 和 ATE 的线性组合。
内部 gmm 结果作为 results_gmm 附加到处理结果中。
默认情况下，处理效果方法计算平均处理效果，其中平均值取自样本观测值。选项 effect_group 可用于计算接受处理者的平均处理效果 ATT，使用 effect_group=1 或未接受处理者的平均处理效果，使用 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]:
| 系数 | 标准误差 | z | P>|z| | 置信区间下限 | 置信区间上限 | |
|---|---|---|---|---|---|---|
| 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: Thu, 03 Oct 2024
Time: 16:04:28
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
==============================================================================
接受处理者的平均处理效果
见下文
[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
==============================================================================
未接受处理者的平均处理效果
[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
==============================================================================
其他计算 ATE 的方法与 ipw 的方法相同或类似，例如回归调整 ra 和双重稳健 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]:
| 系数 | 标准误差 | z | P>|z| | 置信区间下限 | 置信区间上限 | |
|---|---|---|---|---|---|---|
| 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]:
| 系数 | 标准误差 | z | P>|z| | 置信区间下限 | 置信区间上限 | |
|---|---|---|---|---|---|---|
| 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 |
TreatmentEffect 中的所有方法¶
以下计算并打印所有方法的 ATE 和 POM。（我们包含对 TreatmentEffect 的调用作为提醒。）
[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
Stata 中的结果¶
statsmodels 中的结果非常接近 Stata 中的结果，因为这两个软件包使用相同的方法。
[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
没有推断的处理效果¶
可以不计算标准误差和推断统计量来计算 POM 和 ATE。在这种情况下，不会计算 GMM 模型。
[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))
处理对接受处理者的影响¶
子组的处理效果不适用于 aipw 和 aipw-wls。
effect_group 选择计算处理效果和潜在结果的组。选项为“all”表示样本平均处理效果，1 表示接受处理者的平均处理效果，0 表示未接受处理者的平均处理效果。
注意：即使是子组的处理效果，pandas 数据框中的行标签 POM 和 ATE 也是相同的。
[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
处理对未接受处理者的影响¶
与 ATT 类似，我们可以通过使用选项 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
可以使用 help 获取 TreatmentEffect 类及其方法的文档字符串
help(teff)
最后更新：2024 年 10 月 3 日