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

著者: Josef Perktold

このノートブックは、statsmodelsの新しい処置効果機能の基本的な使用方法を説明します。

メインクラスはstatsmodels.treatment.treatment_effects.TreatmentEffectです。

このクラスは、ipw、ra、aipw、aipw-wls、ipw-raの5つの異なるメソッドを使用して、処置効果と潜在的な結果を推定します。最後の3つのメソッドでは、処置または選択モデルと結果モデルの両方が必要です。標準誤差と推論は、選択または処置モデル、結果モデル、および効果関数の結合GMM表現に基づいています。推論のアプローチはStataに従いますが、Stataはより幅広いモデルをサポートしています。推定と推論は、条件付き独立性または無視可能性の下で有効です。

結果モデルは現在、OLSに基づく線形モデルに限定されています。処置は現在、ロジットまたはプロビットのいずれかであるバイナリ処置に制限されています。

この例は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]:
出生時体重 母親が既婚かどうか 母親がヒスパニック系かどうか 父親がヒスパニック系かどうか 外国生まれかどうか アルコール摂取の有無 死亡した子供の数 母親の年齢 母親の教育水準 父親の年齢 ... 出生前ケア 誕生月 出生時体重（ポンド） 女児かどうか 出生前ケア（ダミー変数） 母親の喫煙状況 母親が既婚かどうか（ダミー変数） 女児かどうか(ダミー変数) 出生前ケア（ダミー変数） 母親の年齢の2乗
0 3459 既婚 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 既婚 0 0 1 0 0 22 9 30 ... 1 3 0 いいえ はい 0 1 0 1 484.0
3 2948 既婚 0 0 0 0 0 26 12 30 ... 1 1 0 いいえ はい 0 1 0 1 676.0
4 2410 既婚 0 0 0 0 0 20 12 21 ... 1 3 1 はい はい 0 1 1 1 400.0

5行×28列

TreatmentEffectインスタンスの作成とipwの計算

TreatmentEffectクラスには、- 結果モデルのOLSモデルインスタンス、- 選択モデルの結果インスタンス、- 処置指示変数が必要です。

次の例では、選択モデルとしてプロビットを使用します。ロジットの使用もサポートされています。

[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は、選択モデルと結果モデルの両方を含みます。これらのうち少なくとも1つが正しく指定されていれば、処置効果の一貫した推定値が得られます。ターゲット変数POM0、POM1、およびATEのモーメント条件は、POM0とATEに基づいています。残りのPOM1は、POM0とATEの線形結合として計算されます。

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

デフォルトでは、処置効果メソッドは、標本観測値の平均である平均処置効果を計算します。オプション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]:
係数 標準誤差 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を計算する他の方法（回帰調整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]:
係数 標準誤差 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

TreatmentEffectクラスとそのメソッドのdocstringは、helpを使用して取得できます。

help(teff)

最終更新日：2024年10月3日