傾向スコアと機械学習とprobability calibrationの話
はじめに
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を下記に示します。
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は以下のようになりました。
かなり理想的な線からずれています。そもそもランダムフォレストの予測値では0.1以下または0.7以上の値がありません。これを確率値として使ってしまうのは問題ありそうな気がします。 実際ATEは-280とマイナスの値が出ています。
これを補正するためにpythonだとCalibratedClassifierCVというパッケージがあるのでそれを使って補正を行います。 probability calibrationについては以下を参考にしています。
#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は以下のようになり、ほぼ理想的な直線と一致しました。
また補正後の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は以下のようになります。
ランダムフォレストと同様にマイルドな予測値が出るようになっており、0.4くらいの予測値でも実際は0.2だったり、0.6くらいの予測値の実際の値は0.8だったりしています。
これもCalibratedClassifierCVで補正するとこのような感じになります。
IPWで推定する場合、傾向スコアに0または1が入ると計算できないため、傾向スコアが0.01-0.99の間にあるサンプルに限定してATEを算出しました。ATEの結果は以下のようになりました。
補正前:610
補正後:2440
lightgbmだと補正後の値がかなり大きくなっています。 データの検証をしたところ、傾向スコアが0に近い(CMを視聴する可能性がかなり低い)人で、ゲームのプレイ時間が長い人が散見されました。
原因を探していて下記の記事を見つけたのですが、
IPW推定量だと傾向スコアが0.01と0.5の人がいた場合、0.01の人に50倍の重みをつけることになるため、その部分のサンプルに偏りがあると全体としてのATEが大きく変化してしまうようです。 傾向スコアが0.01とかだと実際はまず割り当てられない→反事実が存在しないサンプルに大きな重みを与えることになり、これはおかしな結果につながりそうです。
この影響を防ぐために傾向スコアが0.1以上の人に限定して再度ATEを算出したところ、1267と他の推定方法に近い値が出ました(これが正しい対処法かどうかはちょっと確証がないです)。層別マッチングでやっても傾向スコア0~0.1の部分のATEがかなり大きくなってしまうのでIPW推定量ほどではないにせよ影響は生じてしまいそうです。 (そもそもこのデータセットはCM接触者に多量のダミーデータが入っているのでそのせいかもしれませんが…)
ちなみに予測後の傾向スコアのヒストグラムは以下のような形になっており、 lightgbmは分類性能が高いために0と1付近と予測されているサンプルが多い傾向があるようです。
今回の結果をまとめると以下のようになりました。
ロジスティック回帰 | ランダムフォレスト | lightgbm | |
---|---|---|---|
AUC | 0.784 | 0.806 | 0.967 |
補正前ATE | 262 | -280 | 610 |
補正後ATE | - | 853 | 2440 1267(傾向スコア0.1以上に限定) |
まとめ
傾向スコアをロジスティック回帰、ランダムフォレスト、lightgbmを用いて推定しました。 ランダムフォレストとlightgbmについては予測確率を実際の確率に近づけるためにprobability calibrationを行い、 それぞれのCM接触効果のATEをIPW推定量を用いて算出しました。
calibrationを行うことでlightgbm等を用いても予測確率を実際の確率に近づけることができましたが、 傾向スコアの分類精度が高い影響で、傾向スコアからATEを算出する部分は少し工夫する必要がありました。
とはいえ傾向スコアの精度はロジスティック回帰よりもlightgbmの方が確実に高くなるので、 傾向スコアからATEを算出する部分をうまく行えれば、より実際に近い効果推定ができる可能性は高いのではないかと思います。
使用したコードは以下になります。