はじめに

RCTが使えない場合の因果推論の手法として傾向スコアを使う方法があります。傾向スコアの算出はロジスティック回帰を用いるのが一般的ですが、この部分は別にlightgbmとか機械学習的な手法でやってもいいのでは？と思っていましたが既にやっている記事がありました。

機械学習系の手法で算出した傾向スコアの値は、どうやら実際の確率とはずれが生じるようで、calibrationを行った方がよいということでした。この部分について、実際にそうなのかが気になったので試してみることにしました。

実装

データは上記記事と同様に岩波データサイエンスvol3のものを使用します。 CM接触に対するゲームのプレイ時間の因果効果を示すことが目的です。まずは一般的なロジスティック回帰を用いた傾向スコアの導出を行います。ここでの傾向スコアはCM接触の割り当てられやすさを示します。

#データの読み込み
df=pd.read_csv("q_data_x.csv")
X=df[['T','F1', 'F2', 'F3', 'M1', 'M2', 'M3', 'TVwatch_day',"child_dummy",'area_kanto', 'area_keihan', 'area_tokai',
       'area_keihanshin','marry_dummy', 'job_dummy1',
       'job_dummy2', 'job_dummy3', 'job_dummy4', 'job_dummy5', 'job_dummy6',
       'job_dummy7', 'job_dummy8']]
y=df["cm_dummy"].values

#ロジスティック回帰
from sklearn.linear_model import LogisticRegression
clf = LogisticRegression()
clf.fit(X, y)
#傾向スコア
pred=clf.predict_proba(X)[:,1]

ロジスティック回帰による予測のcalibration curveを下記に示します。

f:id:rmizutaa:20190330112058p:plain

calibration curveは、予測確率と実際の確率がどれだけ一致しているかを示します。青い点線が理想的な値（予測確率と実際の確率が一致）でオレンジ色の線がロジスティック回帰の予測値を示したものになります。例えば一番左下のオレンジの点は0-0.1と予測されたサンプルの予測値平均は0.083で、実際の値の平均値は0であることを示しています。

calibration curveについて、詳しくはこちらを参照ください。

https://scikit-learn.org/stable/auto_examples/calibration/plot_calibration_curve.html

ロジスティック回帰の場合はほぼ理想的なcurbration curveになっていることがわかります。またこのときのAUCは0.784で、IPW推定量で算出したATE（Average Treatment Score）は262でした。

次にランダムフォレストの場合です。 trainingに使ったデータを使ってそのまま出した予測値はオーバーフィット感があるのでout of foldで予測した値を使います。

#4foldでoof
pred=np.zeros((X.shape[0]))
sss = StratifiedKFold(n_splits=4, random_state=94, shuffle=True)
for i, (train_index, test_index) in enumerate(sss.split(X, y)):
    X_train, X_valid = X.iloc[train_index,:], X.iloc[test_index,:]
    y_train, y_valid = y[train_index], y[test_index]
    
    clf = RandomForestClassifier(n_estimators=100, max_depth=4,random_state=0)
    clf.fit(X_train, y_train)
    pred[test_index]=clf.predict_proba(X_valid)[:,1]

calibrationcurveは以下のようになりました。

f:id:rmizutaa:20190330112321p:plain

かなり理想的な線からずれています。そもそもランダムフォレストの予測値では0.1以下または0.7以上の値がありません。これを確率値として使ってしまうのは問題ありそうな気がします。実際ATEは-280とマイナスの値が出ています。

これを補正するためにpythonだとCalibratedClassifierCVというパッケージがあるのでそれを使って補正を行います。 probability calibrationについては以下を参考にしています。

1.16. Probability calibration — scikit-learn 0.20.3 documentation

#4foldでoof
pred=np.zeros((X.shape[0]))
sss = StratifiedKFold(n_splits=4, random_state=94, shuffle=True)
for i, (train_index, test_index) in enumerate(sss.split(X, y)):
    X_train, X_valid = X.iloc[train_index,:], X.iloc[test_index,:]
    y_train, y_valid = y[train_index], y[test_index]
    
    clf = RandomForestClassifier(n_estimators=100, max_depth=4,random_state=0)
    calibrated_clf = CalibratedClassifierCV(clf, method='isotonic', cv=5)
    calibrated_clf.fit(X_train, y_train)
    calibratedpred[test_index] = calibrated_clf.predict_proba(X_valid)[:,1]

補正後のcalibration curveは以下のようになり、ほぼ理想的な直線と一致しました。

f:id:rmizutaa:20190330112420p:plain

また補正後のATEは853になり、補正を入れることでCM接触効果がプラスになりました。

最後にlightgbmを使った場合です。

lgb_pred=np.zeros((X.shape[0]))
lgb_calibratedpred=np.zeros((X.shape[0]))
sss = StratifiedKFold(n_splits=4, random_state=94, shuffle=True)
for i, (train_index, test_index) in enumerate(sss.split(X, y)):
    X_train, X_valid = X.iloc[train_index,:], X.iloc[test_index,:]
    y_train, y_valid = y[train_index], y[test_index]
    
    model = lgb.LGBMClassifier(boosting_type= 'gbdt',objective = 'binary')
    model.fit(X_train, y_train)
    lgb_pred[test_index]=model.predict_proba(X_valid)[:,1]

    lgb_calibrated_model = CalibratedClassifierCV(model, method='isotonic', cv=5)
    lgb_calibrated_model.fit(X_train, y_train)
    lgb_calibratedpred[test_index] = lgb_calibrated_model.predict_proba(X_valid)[:,1]

補正前のcalibraion_curveは以下のようになります。

f:id:rmizutaa:20190330112632p:plain

ランダムフォレストと同様にマイルドな予測値が出るようになっており、0.4くらいの予測値でも実際は0.2だったり、0.6くらいの予測値の実際の値は0.8だったりしています。

これもCalibratedClassifierCVで補正するとこのような感じになります。

f:id:rmizutaa:20190330112703p:plain

IPWで推定する場合、傾向スコアに0または1が入ると計算できないため、傾向スコアが0.01-0.99の間にあるサンプルに限定してATEを算出しました。ATEの結果は以下のようになりました。

補正前：610
補正後：2440

lightgbmだと補正後の値がかなり大きくなっています。データの検証をしたところ、傾向スコアが0に近い（CMを視聴する可能性がかなり低い）人で、ゲームのプレイ時間が長い人が散見されました。

原因を探していて下記の記事を見つけたのですが、

プロペンシティスコア（Propensity score; PS）（２）－PSに関する5つの誤解 – 医療政策学×医療経済学

IPW推定量だと傾向スコアが0.01と0.5の人がいた場合、0.01の人に50倍の重みをつけることになるため、その部分のサンプルに偏りがあると全体としてのATEが大きく変化してしまうようです。傾向スコアが0.01とかだと実際はまず割り当てられない→反事実が存在しないサンプルに大きな重みを与えることになり、これはおかしな結果につながりそうです。

この影響を防ぐために傾向スコアが0.1以上の人に限定して再度ATEを算出したところ、1267と他の推定方法に近い値が出ました(これが正しい対処法かどうかはちょっと確証がないです)。層別マッチングでやっても傾向スコア0~0.1の部分のATEがかなり大きくなってしまうのでIPW推定量ほどではないにせよ影響は生じてしまいそうです。 (そもそもこのデータセットはCM接触者に多量のダミーデータが入っているのでそのせいかもしれませんが…)

ちなみに予測後の傾向スコアのヒストグラムは以下のような形になっており、 lightgbmは分類性能が高いために0と1付近と予測されているサンプルが多い傾向があるようです。

f:id:rmizutaa:20190330113009p:plain