LoginSignup
12
14

More than 3 years have passed since last update.

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

Last updated at Posted at 2020-03-04

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

基礎編

傾向スコアとは

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

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

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

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

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

傾向スコアの推定

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

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

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. IPW11

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

image.png

ATT in マッチング

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

image.png

ATE in マッチング

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

image.png

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の計算においてコモンサポート外のサンプルを使うと、傾向スコアの逆数(重み)が大きくなりすぎてうまく推定できない12
    • 「コモンサポート外であること」と「重みが大きくなりすぎること」は等価か?
  • コモンサポート外のサンプルは除外するのが一般的?
    • 本書では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説はナンセンス?13
    • AUCは0.8以上あればOKとされている14
    • AUCが小さい(0.5に近い)とはどういうことか:
      • 割り当てが共変量に依存しない($Z\perp X$) → 傾向スコアで調整する必要がない → AUCはある程度大きくないと意味がない?
    • AUCが大きい(1に近い)とはどういうことか:
      • コモンサポートが小さい → 介入効果を推定できない → AUCは大きければ良いというものでもない?

実践編

MineThatData E-Mail Analytics And Data Mining Challenge dataset15

  • 後日追記予定

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')

image.png

_, 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]))

image.png

CPS1+NSW

  • なかなか不均衡なデータである
df_cps1_nsw[z_column].map(treat_label).value_counts().plot.bar(title='# of treatment in CPS1+NSW')

image.png

_, 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]))

image.png

CPS3+NSW

  • NSWほどではないが、CPS1+NSWよりは均衡の取れたデータである
df_cps3_nsw[z_column].map(treat_label).value_counts().plot.bar(title='# of treatment in CPS3+NSW')

image.png

_, 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]))

image.png

回帰分析

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]

image.png

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]

image.png

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]

image.png

結果

データセット 介入効果[$]
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]

image.png

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))

image.png

  • 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]

image.png

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))

image.png

  • 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)

image.png

  • 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]

image.png

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]

image.png

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')

image.png

  • マッチングにおける共変量の分布を確認
    • 平均は真実を語らないので分布も見る
    • 元の分布に比べてちゃんとバランシングできている
_, 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]))

image.png

  • 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]))

image.png

  • 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]))

image.png

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)

image.png

  • 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]

image.png

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]

image.png

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')

image.png

  • マッチングにおける共変量の分布を確認
_, 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]))

image.png

  • 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]))

image.png

  • 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]))

image.png

結果

  • 共変量はバランシングできているが、介入効果の推定値は$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)

image.png

  • 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]

image.png

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]

image.png

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')

image.png

  • マッチングにおける共変量の分布を確認
_, 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]))

image.png

  • 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]))

image.png

  • 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]))

image.png

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)

image.png

  • 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]

image.png

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]

image.png

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')

image.png

  • マッチングにおける共変量の分布を確認
_, 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]))

image.png

  • 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]))

image.png

  • 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]))

image.png

結果

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

所感

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

  1. 本書では傾向スコアは$P(X)$と表記されているが、便宜上$e(X)$とする 

  2. 少し変な比較だが、本書からはこのように読み取れる 

  3. 星野崇宏(2009),調査観察データの統計科学,岩波書店. 

  4. https://www.krsk-phs.com/entry/counterfactual_assumptions 

  5. https://pira-nino.hatenablog.com/entry/casual_inference 

  6. 傾向スコアの文脈でしか見たことがないが、回帰分析でもSUTVAを仮定する? 

  7. 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. 

  8. https://healthpolicyhealthecon.com/2015/05/07/propensity-score-2/ 

  9. https://rmizutaa.hatenablog.com/entry/2019/03/30/224854 

  10. 門脇大輔,阪田隆司,保坂桂佑,平松雄司(2019),Kaggleで勝つデータ分析の技術,技術評論社. 

  11. 傾向スコアを使った介入効果の推定手法はマッチングとIPW以外にもあるが、本書ではこの2つを取り扱っている 

  12. https://healthpolicyhealthecon.com/2015/05/04/propensity-score-1/ 

  13. https://speakerdeck.com/tomoshige_n/causal-inference-and-data-analysis 

  14. 岩波データサイエンス刊行委員会(2016),岩波データサイエンス Vol.3,岩波書店. 

  15. https://blog.minethatdata.com/2008/03/minethatdata-e-mail-analytics-and-data.html 

12
14
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
12
14