効果検証入門 3章傾向スコアを用いた分析+αをPythonでやってみる

Last updated at 2024-05-24Posted at 2020-03-04

この記事は、「効果検証入門～正しい比較のための因果推論／計量経済学の基礎」の「3章傾向スコアを用いた分析」を読んで、学んだことや感じたことを整理したものです。

基礎編

傾向スコアとは

観測された共変量$X$が与えられた条件下で介入群に割り当てられる確率$e(X)=P(Z=1|X)$¹
- 全ての共変量が観測されることはない

回帰分析 v.s. 傾向スコア²

回帰分析の仮定:
- CIA(Conditional Independence Assumption): $\{Y^{(1)}, Y^{(0)}\}\perp Z|X$
- アウトカムと共変量の線形性(本書では省略？)
傾向スコアの仮定
- CIAをごにょごにょしたもの: $\{Y^{(1)}, Y^{(0)}\}\perp Z|e(X)$
  - 割り当てを共変量$X$で説明するのと共変量をバランシングする傾向スコア$e(X)$で説明するのは同じことだと考える³
- SUTVA(Stable Unit Treatment Value Assumption)⁴⁵^sutva-in-ra

CIA v.s. 強く無視できる割り当て

本書ではCIAが使われているが、厳密には強く無視できる割り当てが正しい⁶
強く無視できる割り当てはCIAに次の2つの仮定が加わる
- 共変量$X$は観測されたもののみ
- $0<P(Z=1|X)<1$
  - 要するにどのサンプルも、介入群/非介入群どちらにも割り当てられる可能性がないとNGということで、結構重要だと思う
  - 例えば男性が全員介入群に割り当てられるとすると、$P(Z=1|男性)=1$となるのでNG

傾向スコアの推定

傾向スコアの推定には多くの場合でロジスティック回帰が使われるが、機械学習でもOK
モデルには共変量だけでなく、アウトカムに影響を及ぼすその他の変数も含める⁷
機械学習を使う場合はprobability calibrationが必要⁸
- ランダムフォレストとSVMくらい？
一方、損失関数がloglossなら確率に収束するのでcalibration不要⁹

傾向スコアを使った介入効果の推定

ATT v.s. ATE

ATT(Average Treatment effect on Treated):
- $E[Y^{(1)}-Y^{(0)}|Z=1]$
- 介入を受けたサンプルにおける介入効果の期待値
- 本書ではマッチングによって求める
ATE(Average Treatment Effect):
- $E[Y^{(1)}-Y^{(0)}]$
- サンプル全体の介入効果の期待値
- 本書ではIPWによって求める

マッチング v.s. IPW¹⁰

マッチング:
- 介入群と非介入群から傾向スコアの値が近いサンプルのペアを取り出し、一方の欠測(反事実)をもう一方で補完する
IPW(Inverse Probability Weighting):
- 傾向スコアの逆数でサンプルを重み付けし、母集団の分布に近づける
- IPWはコモンサポート(後述)を考慮しないため、使うべきではない？⁷
  - コモンサポートを考慮すれば使ってもよいと思う

ATT in マッチング

介入群のサンプルに対し傾向スコアの値が近いサンプルを非介入群からマッチング、ペアの差の平均をATTとする

ATE in マッチング

非介入群のサンプルに対しても同様にマッチング、ペアの差の平均をATEとする

ATT in IPW

次式で求める

\frac{\sum_{i=1}^{N}Z_{i}Y_{i}}{\sum_{i=1}^{N}Z_{i}}
-\frac{\sum_{i=1}^{N}\frac{(1-Z_{i})\hat{e}(X_{i})}{1-\hat{e}(X_{i})}Y_{i}}{\sum_{i=1}^{N}\frac{(1-Z_{i})\hat{e}(X_{i})}{1-\hat{e}(X_{i})}}

ATE in IPW

次式で求める

\frac{\sum_{i=1}^{N}\frac{Z_{i}}{\hat{e}(X_{i})}Y_{i}}{\sum_{i=1}^{N}\frac{Z_{i}}{\hat{e}(X_{i})}}
-\frac{\sum_{i=1}^{N}\frac{1-Z_{i}}{1-\hat{e}(X_{i})}Y_{i}}{\sum_{i=1}^{N}\frac{1-Z_{i}}{1-\hat{e}(X_{i})}}

コモンサポートについて

コモンサポートとは、介入群と非介入群の傾向スコアの分布が重なっている範囲
IPWの計算においてコモンサポート外のサンプルを使うと、傾向スコアの逆数(重み)が大きくなりすぎてうまく推定できない¹¹
- 「コモンサポート外であること」と「重みが大きくなりすぎること」は等価か？
コモンサポート外のサンプルは除外するのが一般的？
- 本書ではpropensity score trimmig/clippingを紹介
コモンサポート外を除外することで新たなバイアスを生むのでは？
- だから除外しなくてもいいように、比較可能な2群をちゃんと観測できるように実験計画をデザインすることが重要
- その上で生じたバイアスを傾向スコア等で調整すべき
コモンサポート外では反事実を補完できるサンプルがないため、そもそも「強く無視できる割り当て」を満たしていないのでは？
- $P(Z=1|X)=0$または$P(Z=1|X)=1$でないのはロジスティック回帰を使ってるから

傾向スコアの評価

傾向スコアによって介入群と非介入群の共変量の分布が2群間でバランシングされているかどうかが重要
- つまり、仮定を$\{Y^{(1)}, Y^{(0)}\}\perp Z|e(X)$から$\{Y^{(1)}, Y^{(0)}\}\perp Z|X$に戻すということだと思う
本書では標準化平均差(ASAM)で評価している
- ASAMが0.1以下ならOK(なぜ？)
旧来のAUCが0.8以上ならOK説はナンセンス？¹²
- AUCは0.8以上あればOKとされている¹³
- AUCが小さい(0.5に近い)とはどういうことか:
  - 割り当てが共変量に依存しない($Z\perp X$) → 傾向スコアで調整する必要がない → AUCはある程度大きくないと意味がない？
- AUCが大きい(1に近い)とはどういうことか:
  - コモンサポートが小さい → 介入効果を推定できない → AUCは大きければ良いというものでもない？

実践編

MineThatData E-Mail Analytics And Data Mining Challenge dataset¹⁴

後日追記予定

LaLonde(1986)

データの説明については本書参照
論文の説明についてはこちら参照

データ読み込み

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import warnings
warnings.filterwarnings('ignore')

df_cps1 = pd.read_stata('https://users.nber.org/~rdehejia/data/cps_controls.dta')
df_cps3 = pd.read_stata('https://users.nber.org/~rdehejia/data/cps_controls3.dta')
df_nsw = pd.read_stata('https://users.nber.org/~rdehejia/data/nsw_dw.dta')

データセット作成

NSW: RCTのデータセット
CPS1+NSW: NSWの非介入群をCPS1で置き換えたデータセット
CPS3+NSW: NSWの非介入群をCPS3で置き換えたデータセット

treat_label = {0: 'untreat', 1: 'treat'}
X_columns = ['age', 'education', 'black', 'hispanic', 'married', 'nodegree', 're74', 're75']
y_column = 're78'
z_column = 'treat'

CPS1+NSW

df_cps1_nsw = pd.concat([df_nsw[df_nsw[z_column] == 1], df_cps1], ignore_index=True)

CPS3+NSW

df_cps3_nsw = pd.concat([df_nsw[df_nsw[z_column] == 1], df_cps3], ignore_index=True)

介入群と非介入群の分布を確認

NSW

RCTを謳うだけあって均衡の取れたデータである

df_nsw[z_column].map(treat_label).value_counts().plot.bar(title='# of treatment in NSW')

_, ax = plt.subplots(3, 3, figsize=(16, 12))
for i, c in enumerate(X_columns + [y_column]):
    (df_nsw[c].groupby(df_nsw[z_column].map(treat_label))
     .plot.hist(density=True, alpha=0.5, title=f'{c} densities', legend=True, ax=ax[i // 3, i % 3]))

CPS1+NSW

なかなか不均衡なデータである

df_cps1_nsw[z_column].map(treat_label).value_counts().plot.bar(title='# of treatment in CPS1+NSW')

_, ax = plt.subplots(3, 3, figsize=(16, 12))
for i, c in enumerate(X_columns + [y_column]):
    (df_cps1_nsw[c].groupby(df_cps1_nsw[z_column].map(treat_label))
     .plot.hist(density=True, alpha=0.5, title=f'{c} densities', legend=True, ax=ax[i // 3, i % 3]))

CPS3+NSW

NSWほどではないが、CPS1+NSWよりは均衡の取れたデータである

df_cps3_nsw[z_column].map(treat_label).value_counts().plot.bar(title='# of treatment in CPS3+NSW')

_, ax = plt.subplots(3, 3, figsize=(16, 12))
for i, c in enumerate(X_columns + [y_column]):
    (df_cps3_nsw[c].groupby(df_cps3_nsw[z_column].map(treat_label))
     .plot.hist(density=True, alpha=0.5, title=f'{c} densities', legend=True, ax=ax[i // 3, i % 3]))

回帰分析

import statsmodels.api as sm

NSW

results = sm.OLS(df_nsw[y_column], sm.add_constant(df_nsw[[z_column] + X_columns])).fit()
results.summary().tables[1]

CPS1+NSW

results = sm.OLS(df_cps1_nsw[y_column], sm.add_constant(df_cps1_nsw[[z_column] + X_columns])).fit()
results.summary().tables[1]

CPS3+NSW

results = sm.OLS(df_cps3_nsw[y_column], sm.add_constant(df_cps3_nsw[[z_column] + X_columns])).fit()
results.summary().tables[1]

結果

データセット	介入効果[$]
NSW	1676.3
CPS1+NSW	699.1
CPS3+NSW	1548.2

ゲームスタート

NSWの介入効果は約$1600
バイアスを含んだ状態から傾向スコアを用いてどれだけ$1600に近づけるかチャレンジ

傾向スコアの推定

本書のI(re74^2)とI(re75^2)はどこから出てきた？

CPS1+NSW

傾向スコアをロジスティック回帰で推定

ps_model = sm.Logit(df_cps1_nsw[z_column], sm.add_constant(df_cps1_nsw[X_columns])).fit()
ps_model.summary().tables[1]

df_cps1_nsw['ps'] = 1 / (1 + np.exp(-ps_model.fittedvalues))

傾向スコアの分布を確認

(df_cps1_nsw.groupby(df_cps1_nsw[z_column].map(treat_label))['ps']
 .plot.hist(density=True, alpha=0.5, title=f'propensity score densities', legend=True))

AUCを計算

from sklearn.metrics import roc_auc_score
roc_auc_score(df_cps1_nsw[z_column], df_cps1_nsw['ps'])

0.9708918648513447

CPS3+NSW

傾向スコアをロジスティック回帰で推定

ps_model = sm.Logit(df_cps3_nsw[z_column], sm.add_constant(df_cps3_nsw[X_columns])).fit()
ps_model.summary().tables[1]

df_cps3_nsw['ps'] = 1 / (1 + np.exp(-ps_model.fittedvalues))

傾向スコアの分布を確認

(df_cps3_nsw.groupby(df_cps3_nsw[z_column].map(treat_label))['ps']
 .plot.hist(density=True, alpha=0.5, title=f'propensity score densities', legend=True))

AUCを計算

from sklearn.metrics import roc_auc_score
roc_auc_score(df_cps3_nsw[z_column], df_cps3_nsw['ps'])

0.8742266742266742

ATTの推定

CPS1+NSW

マッチング

ペアを作る関数(再帰ジェネレータ)を定義
- こちらを参考にさせていただいた

from sklearn.neighbors import NearestNeighbors
def pairs_generator(s0: pd.Series, s1: pd.Series, threshold: np.float) -> pd.DataFrame:
    """
    K近傍法を使ってペアを作る再帰ジェネレータ。
    
    Parameters
    ----------
    s0, s1 : pd.Series
    threshold : np.float
        ペアを作る距離の閾値。最近傍点でも閾値を超えた場合はペアとならない。
    
    Returns
    -------
    pairs : pd.DataFrame
        s0のサンプルと、最も近いs1のサンプルで作られたペアのインデックス。
        列lがs0のインデックス、列rがs1のインデックスとなる。
    
    Notes
    -----
    複数のペアに同じサンプルが使われることはない。
    
    Examples
    --------
    s0 = df[df[z_column] == 0]['ps']
    s1 = df[df[z_column] == 1]['ps']
    threshold = df['ps'].std() * 0.1
    pairs = pd.concat(pairs_generator(s1, s0, threshold), ignore_index=True)
    """
    neigh_dist, neigh_ind = NearestNeighbors(1).fit(s0.values.reshape(-1, 1)).kneighbors(s1.values.reshape(-1, 1))
    pairs = pd.DataFrame(
        {'d': neigh_dist[:, 0], 'l': s0.iloc[neigh_ind[:, 0]].index},
        index=s1.index,
    ).query(f'd < {threshold}').groupby('l')['d'].idxmin().rename('r').reset_index()
    if len(pairs) > 0:
        yield pairs
        yield from pairs_generator(s0.drop(pairs['l']), s1.drop(pairs['r']), threshold)

ペアを作る

s1 = df_cps1_nsw[df_cps1_nsw[z_column] == 1]['ps']
s0 = df_cps1_nsw[df_cps1_nsw[z_column] == 0]['ps']
threshold = df_cps1_nsw['ps'].std() * 0.1
pairs1_cps1_nsw = pd.concat(pairs_generator(s1, s0, threshold), ignore_index=True)
pairs0_cps1_nsw = pd.concat(pairs_generator(s0, s1, threshold), ignore_index=True) # ATEの計算で使う

介入群のペアを確認
- オレンジは捨てられたサンプル
- 本来は図のような交差が発生しないようにペアを作るのが望ましい(と思う)

ax = df_cps1_nsw.loc[pairs1_cps1_nsw.unstack()].plot.scatter(x='ps', y=z_column, label='matched', alpha=0.3, figsize=(16, 4), title='pairs on treated')
df_cps1_nsw.loc[~df_cps1_nsw.index.isin(pairs1_cps1_nsw.unstack())].plot.scatter(x='ps', y=z_column, label='unmatched', alpha=0.3, c='tab:orange', ax=ax)
for x in zip(df_cps1_nsw['ps'].iloc[pairs1_cps1_nsw['l']].values, df_cps1_nsw['ps'].iloc[pairs1_cps1_nsw['r']].values):
    plt.plot(x, [1, 0], c='tab:blue', alpha=0.3)

ATTを計算
- 標準偏差がかなり大きい

(pairs1_cps1_nsw['l'].map(df_cps1_nsw[y_column]) - pairs1_cps1_nsw['r'].map(df_cps1_nsw[y_column])).agg(['mean', 'std'])

mean    1929.083008
std     9618.346680
dtype: float64

IPW

重みを計算

df_cps1_nsw['w'] = df_cps1_nsw[z_column] / df_cps1_nsw['ps'] + (1 - df_cps1_nsw[z_column]) / (1 - df_cps1_nsw['ps'])

重み付き線形回帰でATTを推定

att_model = sm.WLS(df_cps1_nsw[y_column], df_cps1_nsw[[z_column]].assign(untreat=1-df_cps1_nsw[z_column]), weights=df_cps1_nsw['w'] * df_cps1_nsw['ps']).fit()
att_model.summary().tables[1]

att_model.params['treat'] - att_model.params['untreat']

1180.4078220960955

IPW(コモンサポート)

傾向スコアの分布が重なっている範囲のサンプルを抽出

df_cps1_nsw_cs = df_cps1_nsw[df_cps1_nsw['ps'].between(*df_cps1_nsw.groupby(z_column)['ps'].agg(['min', 'max']).agg({'min': 'max', 'max': 'min'}))].copy()

重み付き線形回帰でATTを推定

att_model = sm.WLS(df_cps1_nsw_cs[y_column], df_cps1_nsw_cs[[z_column]].assign(untreat=1-df_cps1_nsw_cs[z_column]), weights=df_cps1_nsw_cs['w'] * df_cps1_nsw_cs['ps']).fit()
att_model.summary().tables[1]

att_model.params['treat'] - att_model.params['untreat']

1243.8640545001808

共変量の2群間の標準化平均差を確認
- 調整後は概ね赤い点線の下に収まっている
- 本書のグラフと一致しない？(特に調整前のblack)

pd.DataFrame([{
    'unadjusted': np.subtract(*df_cps1_nsw.groupby(z_column)[c].mean()) / df_cps1_nsw[c].std(),
    'adjusted(matching)': np.subtract(*df_cps1_nsw.loc[pairs1_cps1_nsw.unstack()].groupby(z_column)[c].mean()) / df_cps1_nsw[c].std(),
    'adjusted(weighting)': np.subtract(*df_cps1_nsw.groupby(z_column).apply(lambda x: np.average(x[c], weights=x['w'] * x['ps']))) / df_cps1_nsw[c].std(),
    'adjusted(weighting in common support)': np.subtract(*df_cps1_nsw_cs.groupby(z_column).apply(lambda x: np.average(x[c], weights=x['w'] * x['ps']))) / df_cps1_nsw_cs[c].std(),
} for c in X_columns], index=X_columns).abs().plot(marker='o', alpha=0.5, grid=True, ylim=(0, 100), logy=True, title='covariates balance')
plt.axhline(y=0.1, linestyle='--', c='r')

マッチングにおける共変量の分布を確認
- 平均は真実を語らないので分布も見る
- 元の分布に比べてちゃんとバランシングできている

_, ax = plt.subplots(3, 3, figsize=(16, 12))
for i, c in enumerate(X_columns + [y_column]):
    (df_cps1_nsw.loc[pairs1_cps1_nsw.unstack()][c].groupby(df_cps1_nsw[z_column].map(treat_label))
     .plot.hist(density=True, alpha=0.5, title=f'{c} densities', legend=True, ax=ax[i // 3, i % 3]))

IPWにおける共変量の分布を確認
- 元の分布に比べてちゃんとバランシングできている

_, ax = plt.subplots(3, 3, figsize=(16, 12))
df = df_cps1_nsw.sample(1000000, replace=True, weights=df_cps1_nsw['w'] * df_cps1_nsw['ps'], random_state=2020)
for i, c in enumerate(X_columns + [y_column]):
    (df[c].groupby(df_cps1_nsw[z_column].map(treat_label))
     .plot.hist(density=True, alpha=0.5, title=f'{c} densities', legend=True, ax=ax[i // 3, i % 3]))

IPW(コモンサポート)における共変量の分布を確認
- 元の分布に比べてちゃんとバランシングできている

_, ax = plt.subplots(3, 3, figsize=(16, 12))
df = df_cps1_nsw_cs.sample(1000000, replace=True, weights=df_cps1_nsw_cs['w'] * df_cps1_nsw_cs['ps'], random_state=2020)
for i, c in enumerate(X_columns + [y_column]):
    (df[c].groupby(df_cps1_nsw_cs[z_column].map(treat_label))
     .plot.hist(density=True, alpha=0.5, title=f'{c} densities', legend=True, ax=ax[i // 3, i % 3]))

CPS3+NSW

マッチング

ペアを作る

s1 = df_cps3_nsw[df_cps3_nsw[z_column] == 1]['ps']
s0 = df_cps3_nsw[df_cps3_nsw[z_column] == 0]['ps']
threshold = df_cps3_nsw['ps'].std() * 0.1
pairs1_cps3_nsw = pd.concat(pairs_generator(s1, s0, threshold), ignore_index=True)
pairs0_cps3_nsw = pd.concat(pairs_generator(s0, s1, threshold), ignore_index=True) # ATEの計算で使う

介入群のペアを確認

ax = df_cps3_nsw.loc[pairs1_cps3_nsw.unstack()].plot.scatter(x='ps', y=z_column, label='matched', alpha=0.3, figsize=(16, 4), title='pairs on treated')
df_cps3_nsw.loc[~df_cps3_nsw.index.isin(pairs1_cps3_nsw.unstack())].plot.scatter(x='ps', y=z_column, label='unmatched', alpha=0.3, c='tab:orange', ax=ax)
for x in zip(df_cps3_nsw['ps'].iloc[pairs1_cps3_nsw['l']].values, df_cps3_nsw['ps'].iloc[pairs1_cps3_nsw['r']].values):
    plt.plot(x, [1, 0], c='tab:blue', alpha=0.3)

ATTを計算

(pairs1_cps3_nsw['l'].map(df_cps3_nsw[y_column]) - pairs1_cps3_nsw['r'].map(df_cps3_nsw[y_column])).agg(['mean', 'std'])

mean    1463.115845
std     8779.598633
dtype: float64

IPW

重みを計算

df_cps3_nsw['w'] = df_cps3_nsw[z_column] / df_cps3_nsw['ps'] + (1 - df_cps3_nsw[z_column]) / (1 - df_cps3_nsw['ps'])

重み付き線形回帰でATTを推定

att_model = sm.WLS(df_cps3_nsw[y_column], df_cps3_nsw[[z_column]].assign(untreat=1-df_cps3_nsw[z_column]), weights=df_cps3_nsw['w'] * df_cps3_nsw['ps']).fit()
att_model.summary().tables[1]

att_model.params['treat'] - att_model.params['untreat']

1214.0711733143507

IPW(コモンサポート)

傾向スコアの分布が重なっている範囲のサンプルを抽出

df_cps3_nsw_cs = df_cps3_nsw[df_cps3_nsw['ps'].between(*df_cps3_nsw.groupby(z_column)['ps'].agg(['min', 'max']).agg({'min': 'max', 'max': 'min'}))].copy()

重み付き線形回帰でATTを推定

att_model = sm.WLS(df_cps3_nsw_cs[y_column], df_cps3_nsw_cs[[z_column]].assign(untreat=1-df_cps3_nsw_cs[z_column]), weights=df_cps3_nsw_cs['w'] * df_cps3_nsw_cs['ps']).fit()
att_model.summary().tables[1]

att_model.params['treat'] - att_model.params['untreat']

988.8657151185471

共変量の2群間の標準化平均差を確認

pd.DataFrame([{
    'unadjusted': np.subtract(*df_cps3_nsw.groupby(z_column)[c].mean()) / df_cps3_nsw[c].std(),
    'adjusted(matching)': np.subtract(*df_cps3_nsw.loc[pairs1_cps3_nsw.unstack()].groupby(z_column)[c].mean()) / df_cps3_nsw[c].std(),
    'adjusted(weighting)': np.subtract(*df_cps3_nsw.groupby(z_column).apply(lambda x: np.average(x[c], weights=x['w'] * x['ps']))) / df_cps3_nsw[c].std(),
    'adjusted(weighting in common support)': np.subtract(*df_cps3_nsw_cs.groupby(z_column).apply(lambda x: np.average(x[c], weights=x['w'] * x['ps']))) / df_cps3_nsw_cs[c].std(),
} for c in X_columns], index=X_columns).abs().plot(marker='o', alpha=0.5, grid=True, ylim=(0, 100), logy=True, title='covariates balance')
plt.axhline(y=0.1, linestyle='--', c='r')

マッチングにおける共変量の分布を確認

_, ax = plt.subplots(3, 3, figsize=(16, 12))
for i, c in enumerate(X_columns + [y_column]):
    (df_cps3_nsw.loc[pairs1_cps3_nsw.unstack()][c].groupby(df_cps3_nsw[z_column].map(treat_label))
     .plot.hist(density=True, alpha=0.5, title=f'{c} densities', legend=True, ax=ax[i // 3, i % 3]))

IPWにおける共変量の分布を確認

_, ax = plt.subplots(3, 3, figsize=(16, 12))
df = df_cps3_nsw.sample(1000000, replace=True, weights=df_cps3_nsw['w'] * df_cps3_nsw['ps'], random_state=2020)
for i, c in enumerate(X_columns + [y_column]):
    (df[c].groupby(df_cps3_nsw[z_column].map(treat_label))
     .plot.hist(density=True, alpha=0.5, title=f'{c} densities', legend=True, ax=ax[i // 3, i % 3]))

IPW(コモンサポート)における共変量の分布を確認

_, ax = plt.subplots(3, 3, figsize=(16, 12))
df = df_cps3_nsw_cs.sample(1000000, replace=True, weights=df_cps3_nsw_cs['w'] * df_cps3_nsw_cs['ps'], random_state=2020)
for i, c in enumerate(X_columns + [y_column]):
    (df[c].groupby(df_cps3_nsw_cs[z_column].map(treat_label))
     .plot.hist(density=True, alpha=0.5, title=f'{c} densities', legend=True, ax=ax[i // 3, i % 3]))

結果

共変量はバランシングできているが、介入効果の推定値は$1600に対して結構外れている気がする

データセット	マッチング	IPW	IPW(コモンサポート)
CPS1+NSW	1929.1	1180.4	1243.9
CPS3+NSW	1463.1	1214.1	988.9

ATEの推定

CPS1+NSW

マッチング

ペアを確認

pairs_cps1_nsw = pd.concat([pairs1_cps1_nsw, pairs0_cps1_nsw.rename({'l': 'r', 'r': 'l'}, axis=1)], ignore_index=True)

ax = df_cps1_nsw.loc[pairs_cps1_nsw.unstack()].plot.scatter(x='ps', y=z_column, label='matched', alpha=0.3, figsize=(16, 4), title='pairs')
df_cps1_nsw.loc[~df_cps1_nsw.index.isin(pairs_cps1_nsw.unstack())].plot.scatter(x='ps', y=z_column, label='unmatched', alpha=0.3, c='tab:orange', ax=ax)
for x in zip(df_cps1_nsw['ps'].iloc[pairs_cps1_nsw['l']].values, df_cps1_nsw['ps'].iloc[pairs_cps1_nsw['r']].values):
    plt.plot(x, [1, 0], c='tab:blue', alpha=0.3)

ATEを計算

(pairs_cps1_nsw['l'].map(df_cps1_nsw[y_column]) - pairs_cps1_nsw['r'].map(df_cps1_nsw[y_column])).agg(['mean', 'std'])

mean    1836.516968
std     9634.636719
dtype: float64

IPW

重み付き線形回帰でATEを推定

ate_model = sm.WLS(df_cps1_nsw[y_column], df_cps1_nsw[[z_column]].assign(untreat=1-df_cps1_nsw[z_column]), weights=df_cps1_nsw['w']).fit()
ate_model.summary().tables[1]

ate_model.params['treat'] - ate_model.params['untreat']

-6456.300804768925

IPW(コモンサポート)

重み付き線形回帰でATEを推定

ate_model = sm.WLS(df_cps1_nsw_cs[y_column], df_cps1_nsw_cs[[z_column]].assign(untreat=1-df_cps1_nsw_cs[z_column]), weights=df_cps1_nsw_cs['w']).fit()
ate_model.summary().tables[1]

ate_model.params['treat'] - ate_model.params['untreat']

545.3149727568589

共変量の2群間の標準化平均差を確認

pd.DataFrame([{
    'unadjusted': np.subtract(*df_cps1_nsw.groupby(z_column)[c].mean()) / df_cps1_nsw[c].std(),
    'adjusted(matching)': np.subtract(*df_cps1_nsw.loc[pairs_cps1_nsw.unstack()].groupby(z_column)[c].mean()) / df_cps1_nsw[c].std(),
    'adjusted(weighting)': np.subtract(*df_cps1_nsw.groupby(z_column).apply(lambda x: np.average(x[c], weights=x['w']))) / df_cps1_nsw[c].std(),
    'adjusted(weighting in common support)': np.subtract(*df_cps1_nsw_cs.groupby(z_column).apply(lambda x: np.average(x[c], weights=x['w']))) / df_cps1_nsw_cs[c].std(),
} for c in X_columns], index=X_columns).abs().plot(marker='o', alpha=0.5, grid=True, ylim=(0, 100), logy=True, title='covariates balance')
plt.axhline(y=0.1, linestyle='--', c='r')

マッチングにおける共変量の分布を確認

_, ax = plt.subplots(3, 3, figsize=(16, 12))
for i, c in enumerate(X_columns + [y_column]):
    (df_cps1_nsw.loc[pairs_cps1_nsw.unstack()][c].groupby(df_cps1_nsw[z_column].map(treat_label))
     .plot.hist(density=True, alpha=0.5, title=f'{c} densities', legend=True, ax=ax[i // 3, i % 3]))

IPWにおける共変量の分布を確認

_, ax = plt.subplots(3, 3, figsize=(16, 12))
df = df_cps1_nsw.sample(1000000, replace=True, weights=df_cps1_nsw['w'], random_state=2020)
for i, c in enumerate(X_columns + [y_column]):
    (df[c].groupby(df_cps1_nsw[z_column].map(treat_label))
     .plot.hist(density=True, alpha=0.5, title=f'{c} densities', legend=True, ax=ax[i // 3, i % 3]))

IPW(コモンサポート)における共変量の分布を確認

_, ax = plt.subplots(3, 3, figsize=(16, 12))
df = df_cps1_nsw_cs.sample(1000000, replace=True, weights=df_cps1_nsw_cs['w'], random_state=2020)
for i, c in enumerate(X_columns + [y_column]):
    (df[c].groupby(df_cps1_nsw_cs[z_column].map(treat_label))
     .plot.hist(density=True, alpha=0.5, title=f'{c} densities', legend=True, ax=ax[i // 3, i % 3]))

CPS3+NSW

マッチング

ペアを確認

pairs_cps3_nsw = pd.concat([pairs1_cps3_nsw, pairs0_cps3_nsw.rename({'l': 'r', 'r': 'l'}, axis=1)], ignore_index=True)

ax = df_cps3_nsw.loc[pairs_cps3_nsw.unstack()].plot.scatter(x='ps', y=z_column, label='matched', alpha=0.3, figsize=(16, 4), title='pairs')
df_cps3_nsw.loc[~df_cps3_nsw.index.isin(pairs_cps3_nsw.unstack())].plot.scatter(x='ps', y=z_column, label='unmatched', alpha=0.3, c='tab:orange', ax=ax)
for x in zip(df_cps3_nsw['ps'].iloc[pairs_cps3_nsw['l']].values, df_cps3_nsw['ps'].iloc[pairs_cps3_nsw['r']].values):
    plt.plot(x, [1, 0], c='tab:blue', alpha=0.3)

ATEを計算

(pairs_cps3_nsw['l'].map(df_cps3_nsw[y_column]) - pairs_cps3_nsw['r'].map(df_cps3_nsw[y_column])).agg(['mean', 'std'])

mean    1486.608276
std     8629.639648
dtype: float64

IPW

重み付き線形回帰でATEを推定

ate_model = sm.WLS(df_cps3_nsw[y_column], df_cps3_nsw[[z_column]].assign(untreat=1-df_cps3_nsw[z_column]), weights=df_cps3_nsw['w']).fit()
ate_model.summary().tables[1]

ate_model.params['treat'] - ate_model.params['untreat']

224.6762634665365

IPW(コモンサポート)

重み付き線形回帰でATEを推定

ate_model = sm.WLS(df_cps3_nsw_cs[y_column], df_cps3_nsw_cs[[z_column]].assign(untreat=1-df_cps3_nsw_cs[z_column]), weights=df_cps3_nsw_cs['w']).fit()
ate_model.summary().tables[1]

ate_model.params['treat'] - ate_model.params['untreat']

801.0736196669177

共変量の2群間の標準化平均差を確認

pd.DataFrame([{
    'unadjusted': np.subtract(*df_cps3_nsw.groupby(z_column)[c].mean()) / df_cps3_nsw[c].std(),
    'adjusted(matching)': np.subtract(*df_cps3_nsw.loc[pairs_cps3_nsw.unstack()].groupby(z_column)[c].mean()) / df_cps3_nsw[c].std(),
    'adjusted(weighting)': np.subtract(*df_cps3_nsw.groupby(z_column).apply(lambda x: np.average(x[c], weights=x['w']))) / df_cps3_nsw[c].std(),
    'adjusted(weighting in common support)': np.subtract(*df_cps3_nsw_cs.groupby(z_column).apply(lambda x: np.average(x[c], weights=x['w']))) / df_cps3_nsw_cs[c].std(),
} for c in X_columns], index=X_columns).abs().plot(marker='o', alpha=0.5, grid=True, ylim=(0, 100), logy=True, title='covariates balance')
plt.axhline(y=0.1, linestyle='--', c='r')

マッチングにおける共変量の分布を確認

_, ax = plt.subplots(3, 3, figsize=(16, 12))
for i, c in enumerate(X_columns + [y_column]):
    (df_cps3_nsw.loc[pairs_cps3_nsw.unstack()][c].groupby(df_cps3_nsw[z_column].map(treat_label))
     .plot.hist(density=True, alpha=0.5, title=f'{c} densities', legend=True, ax=ax[i // 3, i % 3]))

IPWにおける共変量の分布を確認

_, ax = plt.subplots(3, 3, figsize=(16, 12))
df = df_cps3_nsw.sample(1000000, replace=True, weights=df_cps3_nsw['w'], random_state=2020)
for i, c in enumerate(X_columns + [y_column]):
    (df[c].groupby(df_cps3_nsw[z_column].map(treat_label))
     .plot.hist(density=True, alpha=0.5, title=f'{c} densities', legend=True, ax=ax[i // 3, i % 3]))

IPW(コモンサポート)における共変量の分布を確認

_, ax = plt.subplots(3, 3, figsize=(16, 12))
df = df_cps3_nsw_cs.sample(1000000, replace=True, weights=df_cps3_nsw_cs['w'], random_state=2020)
for i, c in enumerate(X_columns + [y_column]):
    (df[c].groupby(df_cps3_nsw_cs[z_column].map(treat_label))
     .plot.hist(density=True, alpha=0.5, title=f'{c} densities', legend=True, ax=ax[i // 3, i % 3]))

結果

データセット	マッチング	IPW	IPW(コモンサポート)
CPS1+NSW	1836.5	-6456.3	545.3
CPS3+NSW	1486.6	224.7	801.0

所感

全然「バイアス」を取り除けてない…

本書では傾向スコアは$P(X)$と表記されているが、便宜上$e(X)$とする ↩
少し変な比較だが、本書からはこのように読み取れる ↩
星野崇宏(2009)，調査観察データの統計科学，岩波書店． ↩
https://www.krsk-phs.com/entry/counterfactual_assumptions ↩
https://pira-nino.hatenablog.com/entry/casual_inference ↩
Paul R. Rosenbaum and Donald B. Rubin. "The Central Role of the Propensity Score in Observational Studies for Causal Effects." Biometrika, Vol. 70, pp. 41-55. ↩
https://healthpolicyhealthecon.com/2015/05/07/propensity-score-2/ ↩ ↩²
https://rmizutaa.hatenablog.com/entry/2019/03/30/224854 ↩
門脇大輔，阪田隆司，保坂桂佑，平松雄司(2019)，Kaggleで勝つデータ分析の技術，技術評論社． ↩
傾向スコアを使った介入効果の推定手法はマッチングとIPW以外にもあるが、本書ではこの2つを取り扱っている ↩
https://healthpolicyhealthecon.com/2015/05/04/propensity-score-1/ ↩
https://speakerdeck.com/tomoshige_n/causal-inference-and-data-analysis ↩
岩波データサイエンス刊行委員会(2016)，岩波データサイエンス Vol．3，岩波書店． ↩
https://blog.minethatdata.com/2008/03/minethatdata-e-mail-analytics-and-data.html ↩

You get articles that match your needs
You can efficiently read back useful information
You can use dark theme

What you can do with signing up

効果検証入門 3章 傾向スコアを用いた分析+αをPythonでやってみる

基礎編

傾向スコアとは

回帰分析 v.s. 傾向スコア2

CIA v.s. 強く無視できる割り当て

傾向スコアの推定

傾向スコアを使った介入効果の推定

ATT v.s. ATE

マッチング v.s. IPW10

ATT in マッチング

ATE in マッチング

ATT in IPW

ATE in IPW

コモンサポートについて

傾向スコアの評価

実践編

MineThatData E-Mail Analytics And Data Mining Challenge dataset14

LaLonde(1986)

データ読み込み

データセット作成

CPS1+NSW

CPS3+NSW

介入群と非介入群の分布を確認

NSW

CPS1+NSW

CPS3+NSW

回帰分析

NSW

CPS1+NSW

CPS3+NSW

結果

ゲームスタート

傾向スコアの推定

CPS1+NSW

CPS3+NSW

ATTの推定

CPS1+NSW

マッチング

IPW

IPW(コモンサポート)

CPS3+NSW

マッチング

IPW

IPW(コモンサポート)

結果

ATEの推定

CPS1+NSW

マッチング

IPW

IPW(コモンサポート)

CPS3+NSW

マッチング

IPW

IPW(コモンサポート)

結果

所感

効果検証入門 3章傾向スコアを用いた分析+αをPythonでやってみる

回帰分析 v.s. 傾向スコア²

マッチング v.s. IPW¹⁰

MineThatData E-Mail Analytics And Data Mining Challenge dataset¹⁴