効果検証入門 3.傾向スコアを用いた分析

前回：効果検証入門 2.介入効果を測るための回帰分析

回帰分析では、共変量で条件づけられた介入の決定が完全にランダム(CIA)であるようなモデルがセレクションバイアスを軽減した。今回はCIAを構成することが(介入の割り当て確率が入手できないといった理由から)困難である場合について、傾向スコアを使った分析を学ぶ。

介入変数Zを目標変数としてロジスティック回帰することで、共変量Xが与えられた時の介入の確率(傾向スコア。交絡因子ポイントといったところか。)を推定する。介入確率をP(X_i)を表記することにする。

\hat{P}(X_i)=\sigma(\hat{\beta}X_i)

これにより、複数のデータの共変量が与えられた時、その傾向スコアが一致するならばそのグループ内では介入はランダムだと仮定することができる。このように、CIAよりも緩和された条件を用いて分析ができる。

※傾向スコアが類似するデータをクラスタリングすれば良いので、決定木のような機械学習手法を用いても良い。(大体の論文はロジスティック回帰で算出していた。もしかすると傾向スコアの算出方法についての研究はあまり進んでいないのかも。わざわざニューラルネットワークで傾向スコア算出したりしても牛刀割鶏だからか？？時間を作って検証したい。パッと思いつく手法↓)

• 主成分回帰：結局やりたいことは「多くの交絡スコアを一次元に集約し調整する」ということなので次元削減してからロジスティック回帰するのは有効なのでは。そもそも今回のサンプルサイズに対して説明変数8個でロジスティック回帰するのは微妙。参考論文
• ランダムフォレスト：誤って介入決定に全く影響しない変数をモデルに入れてしまった場合、線形性を表現できないロジスティック回帰では厳しい。決定木は特徴量の軸と直行する決定境界を作るので、そのような場合には対応しやすい。
• ニューラルネット：私の不勉強で弱点をよく知らないので載せておく。データ数が多ければニューラルネットでなんとかなるイメージなのですが有識者の方教えてください。

※傾向スコアの評価方法としては、Hosmer-Lemeshow検定やAUCがあるらしい。分類器の指標なので普通に混合行列などでも良いのでは？

3.1 傾向スコアマッチング

• ATT(Average Treatment effect on Treated)：介入を受けたサンプルにおける介入効果の期待値。一般に、介入を受けたサンプルでは介入がなかったと仮定した時のYの値は得られない。そこで傾向スコアが類似している介入を受けていないサンプルを持ってきて、両者の差をとることでATTを推定する。

3.1.1 LaLondeデータセット

LaLondeデータセットを分析する。LoLondeデータセットはdtaファイルなので読み込みが若干面倒。こちらのサイトを参考にした。まずはターミナルでStatsModelsをインストールする。(数分かかった)

$ pip install git+https://github.com/statsmodels/statsmodels

StatsModelsを使ってdtaファイルを読み込む。

import statsmodels.iolib.foreign as smio
import pandas as pd
import matplotlib.pyplot as plt
import random
from sklearn.linear_model import LogisticRegression

CPS_data = pd.DataFrame.from_records(smio.genfromdta('cps_controls.dta'))
NSW_data = pd.DataFrame.from_records(smio.genfromdta('nsw_dw.dta'))

このデータセットは失業者への就職支援の効果を調べるもの。'treat'が介入変数であり、これが1なら就職支援を行った、ということになる。Yは're78'で、就職支援によって年収が変化するかを調査する。他にも'education'(学歴)や'married'(結婚)などの変数が含まれている。

NSWはRCT(就職支援の対象をランダムに選択)である。それに対し、CPSは実験の外で得られた調査データであり、これを非介入群とすることでセレクションバイアスを擬似的に発生させることができる。

1. NSWのうち介入群だけを取り出す。
2. 介入群と非介入群(CPS)について、're78'(年収)のヒストグラムを作成

#---
#　NSWからの介入群の抽出
#---

random.seed(1234)
Treated = NSW_data[NSW_data['treat'] == 1].copy()
index = sorted( random.sample(range(len(CPS_data)), k = len(NSW_treated)) ) #介入と非介入のサンプル数を揃える
Not_Treated = CPS_data.iloc[index,:].copy()

#---
#　プロット
#---

plt.hist([Treated['re78'], Not_Treated['re78']], bins = 20,
         histtype='barstacked', label = ['Treated', 'Not Treated'])
plt.xlabel('re78')
plt.legend()
plt.show()

(なぜ.copy()してるのかについてはこちら参照。copy()しないと後でたくさんWarningが出る。)

ヒストグラムから分かるように、非介入群の方が年収が高い、という結果になっている。これは「就職支援を行うと年収が下がる」という不自然な結論を導きかねない。しかし、ここでは「元々就職支援を必要としない人」が非介入群に含まれており、単純に就職支援の効果を分析することができない。実際、NSWでは

となっている。そこで、傾向スコアマッチングによる分析を行う。

3.1.2 傾向スコアマッチングの実装

1. TreatedとNot_Treatedを連結したjoin_dataを作成する。
2. join_dataに対してtreatを被説明変数としたロジスティック回帰を実施する。
3. TreatedとNot_Treatedに傾向スコアの列を追加し、傾向スコアに基づきソートする。
4. Treatedの傾向スコアと一致するNot_Treatedのデータと順探索でマッチングさせる。
5. 共変量のバランス(二群における各共変量の平均の差)を確認する。
6. マッチング後のデータに対し、ATTを算出する。

#---
#　TreatedとNot_Treatedを連結
#---

join_data = pd.concat([Treated, Not_Treated])
X_list = ['age', 'education', 'black', 'hispanic', 'married', 'nodegree', 're74', 're75']
join_X = join_data[X_list].copy()
join_Z = join_data[['treat']].copy()

#---
#　ロジスティック回帰の実施
#---

model = LogisticRegression(max_iter = 1000, multi_class = 'ovr')
model.fit(join_X, join_Z.values.ravel())

#---
#　傾向スコアの算出、傾向スコアによるソート
#---

Treated['PropensityScore'] = model.predict(Treated[X_list])
Not_Treated['PropensityScore'] = model.predict(Not_Treated[X_list])
sorted_Treated = Treated.sort_values('PropensityScore', ascending=False)
sorted_Not_Treated = Not_Treated.sort_values('PropensityScore', ascending=False)

#---
#　マッチング
#---

matched_list_Treated, matched_list_Not_Treated = [],[]
top = 0
for i in range(0, len(Treated)):
    if sorted_Treated.iloc[i,:]['PropensityScore'] == 0:
        break
    for j in range(top, len(Not_Treated)):
        if sorted_Treated.iloc[i,:]['PropensityScore'] == sorted_Not_Treated.iloc[j,:]['PropensityScore']:
            matched_list_Treated.append(i)
            matched_list_Not_Treated.append(j)
            top = j+1
            break
        else:
            break
matched_Treated = sorted_Treated.iloc[matched_list_Treated,:]
matched_Not_Treated = sorted_Not_Treated.iloc[matched_list_Not_Treated,:]
            
#---
#　共変量のバランスをプロット
#---

mean_Treated = Treated[X_list].describe().loc['mean'].values
mean_Not_Treated = Not_Treated[X_list].describe().loc['mean'].values
var_Treated = Treated[X_list].describe().loc['std'].values ** 2
var_Not_Treated = Not_Treated[X_list].describe().loc['std'].values ** 2
Unadjusted_Balance = abs(mean_Treated - mean_Not_Treated) / np.sqrt((var_Treated + var_Not_Treated) / 2)

mean_matched_Treated = matched_Treated[X_list].describe().loc['mean'].values
mean_matched_Not_Treated = matched_Not_Treated[X_list].describe().loc['mean'].values
var_matched_Treated = matched_Treated[X_list].describe().loc['std'].values ** 2
var_matched_Not_Treated = matched_Not_Treated[X_list].describe().loc['std'].values ** 2
Adjusted_Balance = abs(mean_matched_Treated - mean_matched_Not_Treated) / np.sqrt((var_matched_Treated + var_matched_Not_Treated) / 2)


plt.plot(X_list, Unadjusted_Balance, marker = 'o', linestyle = '-', label = 'Unadjusted')
plt.plot(X_list, Adjusted_Balance, marker = 'o', linestyle = '-', label = 'Adjusted')
plt.legend()
plt.ylabel('Covariate Balance')
plt.show()

#---
#　ATTの算出
#---

ATE_NSW = NSW_data[NSW_data['treat'] == 1].copy()['re78'].mean() - NSW_data[NSW_data['treat'] == 0].copy()['re78'].mean()
ATT = matched_Treated['re78'].mean() - matched_Not_Treated['re78'].mean()
print(f'ATE of NSW data = {ATE_NSW:.2f}, ATT of biased Data = {ATT:.2f}')

###出力###
ATE of NSW data = 1794.35, ATT of biased Data = 1605.57

このように、セレクションバイアスのあるデータに対してRCTデータに近い精度で介入効果を測定することができた。しかし、気になる点としては傾向スコアが0か1の値しか取っていないこと。Xを標準化したところ明らかに変な結果になってしまったのでやめたが、普通対称群の傾向スコアが重なる部分のサンプルをマッチングするため今回は傾向スコア=1のサンプルのみマッチングすることになってしまった。

3.3 逆確率重み付き推定

介入群の結果の平均を取る際に、介入確率(傾向スコア)の逆数で重み付することによって、セレクションバイアスによって真の平均を推定する。しかし、今回のように傾向スコアが大変低いデータが多数存在する場合は、重みが発散してしまい有効ではない。(良いデータが見つかれば実装したい)