条件付き独立性の下での処置効果

著者: Josef Perktold

このノートブックは、statsmodelsにおける新しい処置効果機能の基本的な使用法を示しています。

主要なクラスは statsmodels.treatment.treatment_effects.TreatmentEffect です。

このクラスは、ipw、ra、aipw、aipw-wls、ipw-raの5つの異なる方法を使用して処置効果と潜在的な結果を推定します。最後の3つの方法は、処置または選択モデルと結果モデルの両方を必要とします。 標準誤差と推論は、選択または処置モデル、結果モデル、および効果関数の共同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 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

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 インスタンスを作成した後、五つのメソッドのいずれかを呼び出して、潜在的な結果である POM0、POM1 、平均処置効果（ATE）を計算できます。POM0 は処置なしのグループの潜在的な結果、POM1 は処置グループの潜在的な結果、処置効果はPOM1 - POM0です。

例えば、 teff.ipw() は逆確率加重法を使用してPOMとATEを計算します。処置の確率は、一般に傾向スコアとも呼ばれます。推定の summary には、POMとATEの標準誤差と信頼区間が含まれます。

標準誤差およびその他の推論統計量は、選択モデルと結果モデルの一般化方法（GMM）表現、および結果統計量のモーメント条件に基づいています。 メソッド ipw は選択モデルを使用しますが、結果モデルは使用しません。 メソッド ra は結果モデルを使用しますが、選択モデルは使用しません。 二重ロバスト推定量である aipw 、 aipw-wls 、 ipw-ra は選択モデルと結果モデルの両方を含み、そのうち少なくとも1つが正しく指定されていないと処置効果の一貫した推定は得られません。 ターゲット変数であるPOM0、POM1、およびATEのモーメント条件は、POM0とATEに基づいています。残りのPOM1は、POM0とATEの線形結合として計算されます。

内部のGMM結果は、 results_gmm として処置結果に添付されます。

デフォルトでは、処置効果メソッドは平均処置効果（ATE）を計算し、サンプル観測値の平均を取ります。 オプション effect_group を使用すると、 effect_group=1 で処置を受けたグループに対する平均処置効果（ATT）を計算したり、 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.048541e-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, 23 Mar 2024
Time:                        10:53:17
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
==============================================================================

処置を受けた対象に対する平均処置効果（ATT）

下記をご覧ください。

[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を計算する他の方法は、例えば回帰調整（ra）や二重ロバスト（ipw_ra）のように、 ipw と同様の方法で機能します。

[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.048541e-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.038644e-16     -280.843822      -173.547414
POM0  3403.250651   9.596571  354.631943  0.000000e+00     3384.441717      3422.059585
POM1  3176.055033  25.654641  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
(-239.6392114643363, 3403.2422719354845, 3163.603060471148)

 ipw
(-230.6886377952615, 3403.462708684556, 3172.774070889295)

 aipw
(-230.98920111257803, 3403.355253173835, 3172.366052061257)

 aipw_wls
(-227.19561818674856, 3403.2506509757864, 3176.055032789038)

 ipw_ra
(-229.96707793513042, 3403.335639307419, 3173.3685613722887)

処置群への処置効果

aipw および aipw-wls では、サブグループへの処置効果は利用できません。

effect_group は、処置効果と潜在的な結果が計算されるグループを選択します。
オプションは次の通りです：
「all」 はサンプル平均処置効果（SATE）を、
1 は処置群への平均処置効果（ATT）を、
0 は非処置群への平均処置効果（ATU）を指定します。

注：pandas のDataFrame、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.702293e-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.786473  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.633805  0.000000e+00     3394.716358      3431.106829
POM1  3181.785621  27.301318  116.543297  0.000000e+00     3128.276021      3235.295221

TreatmentEffect クラスとそのメソッドのドクストリングは、次のコマンドを使用して取得できます。

help(teff)

最終更新日: 2025年01月28日