條件獨立下的處理效應¶
作者：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 | notmarried | 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]:
| 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: 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]:
| 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 |
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 類別及其方法的 docstring
help(teff)
上次更新：2024 年 10 月 03 日