下町データサイエンティストの日常

下町データサイエンティストの成果物

傾向スコアを用いた因果推論入門~実装編~

0. はじめに

理論編に引き続き、実装を行なっていきます。

pira-nino.hatenablog.com

実装はpythonとRの両方で行いました。

作業用のコードは以下のgithubで公開していますが、基本ブログ内のコードで完結することを心がけています。
また、基礎集計などの因果推論のメインでない部分はpythonのみの実装例を示します。

github.com

そもそも作業用コードなのでリファクタはちゃんとしていないことに加え、ところどころ怪しい箇所があるので指摘を貰えたら嬉しいです。普段私はpythonユーザでありRなんて滅多に触らないので、特にRに関してはあまり自信がないです。
加えて、pythonにはマッチングの便利パッケージがないので自作の部分が多いです。
一応、以下のパッケージは触ってみましたがあまりしっくりきませんでした*1。結局傾向スコアマッチングに関してはDoWhyパッケージに自作で手を加えた物を作りました汗

github.com

github.com

追記：ブログを書いている途中に思ったのですが、「CVを切る」「ROC曲線は」といった感じで機械学習の基礎的単語は断りなく使っているので察していただければ幸いです汗

1. 対象データ

みんな大好き岩波D3 vol.3のデータを使いました。

github.com

ストーリとしては、「CMを見た人は見ていない人に比べ本当にゲームアプリの利用時間が増えたか」といった問題に対して因果推論のアプローチをとるといった感じです。 詳しくは書籍を買って読んできださい。

www.iwanami.co.jp

個人の感想ですが、「このデータめっちゃ綺麗だな。。。」と思いました。
こんなに因果推論の諸々がうまくいく例はそうそうないという認識を持っていただければ幸いです。
実務で因果推論を行う際には、「強く無視できる割り当ては大丈夫か？」「そもそも割り当ての定義は本当に現状でいいか？」をよく吟味してから初めて因果推論を行なってください。

2. 先に各手法の結果を確認

この後諸々の手法の実装をしていきますが、先に結果を載せます。

手法 実装 因果効果 SD(TVwatch_day) 備考
傾向スコアマッチング
(ロジスティック回帰) 		python 272 0.087
傾向スコアマッチング
(ロジスティック回帰) 		R 181 0.0093
IPW
(ロジスティック回帰) 		python 470 -
IPW
(ロジスティック回帰) 		R 1,025 -
IPW
(書籍の結果) 		R 1,503 - 書籍と同じ変数を用いた結果の再現
ダブルロバスト
(ロジスティック回帰 /Ridge回帰) 		python 340 -
ダブルロバスト
(ロジスティック回帰 /Ridge回帰) 		R 927 -
Proximityマッチング
(RF) 		python 2,478 0.71 insampleのモデリング
Proximityマッチング
(RF) 		R 2,478 0.29

手法により結構値がブレました。。。
では、どの手法を使用すべきかという判断が難しいのも因果推論の一つの悩みです。例えば、SDが小さいマッチングを信用すべきなのは納得がいきますが、IPWとダブルロバストはそれぞれを構成する傾向スコアモデリングや回帰モデリングの精度に依存します。 個人的にはマッチングが無難に信用できる手法かなぁとも思っています。
また詳しい私の考えは最後に書きますが、上記の値を出す手順を以下で述べていきます。

3.準備

python

import lightgbm as lgb
import matplotlib.font_manager as fm
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import seaborn as sns
from sklearn.ensemble import RandomForestClassifier
from sklearn.linear_model import LinearRegression, LogisticRegression
from sklearn.metrics import log_loss, mean_squared_error, roc_auc_score
from sklearn.model_selection import KFold, train_test_split
from sklearn.metrics import roc_curve,auc
from sklearn.neighbors import NearestNeighbors


%matplotlib inline
# その他設定
pd.set_option('display.max_columns', 100)
plt.style.use('seaborn-darkgrid')
font = fm.FontProperties(fname='./TakaoPGothic.ttf')

dataset_url = 'https://raw.githubusercontent.com/iwanami-datascience/vol3/master/kato%26hoshino/q_data_x.csv'  
df = pd.read_csv(dataset_url)  
print(df.shape)

(10000, 35)

R

library(cvTools)
library(Metrics)
library(ROCR)
library(MLmetrics)
library(glmnet)
library(Matching)

df <- read.csv("https://raw.githubusercontent.com/iwanami-datascience/vol3/master/kato%26hoshino/q_data_x.csv")
print(dim(df))

[1] 10000    35

以下の変数達が分析で主に用いる変数ですので、頭に入れておいて下さい。

  • 目的変数(y): gamesecond
  • 割り当て(z): cm_dummy

4. 基礎集計~CMを見た人と見ていない人での目的変数の差~

fig = plt.figure(figsize=(7, 7))
ax = fig.add_subplot(1, 1, 1)
ax = sns.barplot(x='cm_dummy', y='gamesecond', data=df, ax=ax)
ax.set_xticklabels(['CM見てない(z=0)', 'CM見た(z=1)'], fontproperties=font)
ax.set_xlabel('')
ax.set_ylabel('平均ゲーム利用時間 [秒]', fontproperties=font)
fig.show()

f:id:pira_nino:20190714210821p:plain

上図より、単純に集計するとCMを見ていない人の方がゲームのプレイ時間が長いという解釈になってしまいます。
しかし、実際にはCMに接触している人ほど、女性が多い / 年齢層が高い...といったようにセレクションバイアス*2がかかっています。これ故に、CMを見た or 見ないで比較を行うことは妥当であるとは言えません。
そこで因果推論を用いてCMの効果を推定しようというのが導入となります。

5. 因果推論の流れのおさらい

理論編でも載せましたが以下の順序で分析を行います。 f:id:pira_nino:20190616163743p:plain

6. 傾向スコアモデリング

傾向スコアモデリングはzの2値分類モデルを作成し、確率P(z|x)を求めるタスクと同等です。

6.1 モデリング

まずpythonの実装をします。
今回は、「ロジスティック回帰」「ランダムフォレスト(RF) 」「LightGBM」の実装を行いました。
なお、そこまで効果が期待できないと思われるので、パラメータチューニングは行ないませんでした。 書籍では変数選択をしていますが、簡単のために全変数を用いています。 加えて、結果がそれっぽかったので、カテゴリカル変数(area_kanto等)を特にcategorical_featureとして扱うといった処理も面倒なので省略しています汗 また、念のためにCVを切って学習をさせます。

20190805追記
LightGBMのCVの切り方 / 予測方法にご指摘を頂きましたので修正します。

u++さんありがとうございます！！！

upura.hatenablog.com

学習はこんな感じでCVを切って行いました。

f:id:pira_nino:20190805205250p:plain

尚、テストデータの予測結果のAUCは変わらず1.00であったので以下のコード以降は変更していないです。
また、ダブルロバストの回帰については、ちゃんと考えると実装に手間がかかり、かつ最終結果にRidge回帰を用いているので修正を行なっていません。

準備

z_col = 'cm_dummy'
x_cols = [
    col for col in df.columns
    if col not in ['gamecount', 'gamedummy', 'gamesecond', 'cm_dummy']
]
df_x = df[x_cols]
df_z = df[z_col]

モデリングのメイン部分

test_ratio = 0.2
oof_preds_lgb = np.zeros(df.shape[0])
oof_preds_lr = np.zeros(df.shape[0])
oof_preds_rf = np.zeros(df.shape[0])
insample_logloss_list_lgb = []
insample_logloss_list_lr = []
insample_logloss_list_rf = []
feature_importance_df = pd.DataFrame()

kf = KFold(n_splits=5, shuffle=True, random_state=0)
for n_fold, (train_index, valid_index) in enumerate(kf.split(df_x)):
    print('CV_{}'.format(n_fold))
    X_train, z_train = df_x.values[train_index], df_z.values[train_index]
    X_test, z_test = df_x.values[valid_index], df_z.values[valid_index]

    # train -> TRAIN / VALに分割
    num_test = int(len(X_train) * test_ratio)
    index_all = np.random.choice(len(X_train),  len(X_train), replace=False)
    index_VAL = index_all[0:num_test]
    index_TRAIN = index_all[num_test:len(X_train)]
    X_TRAIN, z_TRAIN = X_train[index_TRAIN], z_train[index_TRAIN]
    X_VAL, z_VAL = X_train[index_VAL], z_train[index_VAL]
    lgb_train = lgb.Dataset(X_TRAIN, label=z_TRAIN, free_raw_data=False)
    lgb_val = lgb.Dataset(X_VAL, label=z_VAL, free_raw_data=False)

    # LightGBM
    print('LightGBM_fitting...')
    params = {
        'task': 'train',
        'objective': 'binary',
        'metric': 'binary',
    }
    model_lgb = lgb.train(
        params,
        lgb_train,
        valid_sets=[lgb_train, lgb_val],
        valid_names=['train', 'val'],
        num_boost_round=100000,
        early_stopping_rounds=100,
        verbose_eval=1000)
    # LogisticRegression
    print('LogisticRegression_fitting...')
    model_lr = LogisticRegression()
    model_lr.fit(X_train, z_train)

    #RF
    print('RandomForestClassifier_fitting...')
    model_rf =  RandomForestClassifier(n_estimators=500)
    model_rf.fit(X_train, z_train)
    
    # pred_for_insample
    # NOTE: 正確にはLightGBMではTrain ->(TRAIN., VAL)に分割している
    pred_insample_lgb = model_lgb.predict(X_train, num_iteration=model_lgb.best_iteration)
    pred_insample_lr =   model_lr.predict_proba(X_train)[:, 1]
    pred_insample_rf =   model_rf.predict_proba(X_train)[:, 1]
    insample_logloss_lgb = (log_loss(z_train, pred_insample_lgb))
    insample_logloss_lr = (log_loss(z_train, pred_insample_lr))
    insample_logloss_rf = (log_loss(z_train, pred_insample_rf))
    print('insample_logloss_of_LightGBM:{}'.format(insample_logloss_lgb))
    print('insample_logloss_of_LogisticRegression:{}'.format(insample_logloss_lr))
    print('insample_logloss_of_RandomForestClassifier:{}'.format(insample_logloss_rf))
    insample_logloss_list_lgb.append(insample_logloss_lgb)
    insample_logloss_list_lr.append(insample_logloss_lr)
    insample_logloss_list_rf.append(insample_logloss_rf)

    # pred_for_oof
    oof_pred_lgb = model_lgb.predict(X_test, num_iteration=model_lgb.best_iteration)
    oof_pred_lr = model_lr.predict_proba(X_test)[:, 1]
    oof_pred_rf = model_rf.predict_proba(X_test)[:, 1]
    
    oof_preds_lgb[valid_index] = oof_pred_lgb
    oof_preds_lr[valid_index] = oof_pred_lr
    oof_preds_rf[valid_index] = oof_pred_rf

    oof_logloss_lgb = (log_loss(z_test, oof_pred_lgb))
    oof_logloss_lr = (log_loss(z_test, oof_pred_lr))
    oof_logloss_rf = (log_loss(z_test,oof_pred_rf))
    print('oof_logloss_of_LightGBM:{}'.format(oof_logloss_lgb))
    print('oof_logloss_of_LogisticRegression:{}'.format(oof_logloss_lr))
    print('oof_logloss_of_RandomForestClassifier:{}'.format(oof_logloss_rf))


    # 重要度の算出
    fold_importance_df = pd.DataFrame()
    fold_importance_df["feature"] = x_cols
    fold_importance_df["importance"] = np.log1p(
        model_lgb.feature_importance(
        importance_type='gain', iteration=model_lgb.best_iteration))
    fold_importance_df["fold"] = n_fold + 1
    feature_importance_df = pd.concat(
        [feature_importance_df, fold_importance_df], axis=0)

print('all_CV_done')
print('CV_insample_logloss_of_LightGBM:{}'.format(np.mean(insample_logloss_list_lgb)))
print('CV_insample_logloss_of_LogisticRegression:{}'.format(np.mean(insample_logloss_list_lr)))
print('CV_insample_logloss_of_RandomForestClassifier:{}'.format(np.mean(insample_logloss_list_rf)))
print('CV_oof_logloss_of_LightGBM:{}'.format(log_loss(df_z, oof_preds_lgb)))
print('CV_oof_logloss_of_LogisticRegression:{}'.format(log_loss(df_z, oof_preds_lr)))
print('CV_oof_logloss_of_RandomForestClassifier:{}'.format(log_loss(df_z, oof_preds_rf)))

# #予測値をdfに追加
df['pred_lgb'] = oof_preds_lgb
df['pred_lr'] = oof_preds_lr
df['pred_rf'] = oof_preds_rf

出力

........(略)

all_CV_done
CV_insample_logloss_of_LightGBM:0.0016683426433416388
CV_insample_logloss_of_LogisticRegression:0.5396573440875929
CV_insample_logloss_of_RandomForestClassifier:0.0017502169196193698
CV_oof_logloss_of_LightGBM:0.0032708288644122943
CV_oof_logloss_of_LogisticRegression:0.5421404243283353
CV_oof_logloss_of_RandomForestClassifier:0.004715561791194972

見やすく表にまとめます。

手法 insample_logloss oof_logloss
ロジスティック回帰 0.540 0.542
RF 0.001 0.004
LightGBM 0.002 0.003

余談ですが、こんな感じでモデリングのテンプレ関数を持っておくと日々の業務で役立ちます。
話は戻り、出力を見る感じ各手法ともinsampleとoofの結果に大きな差がないあたりからCVの切り方は大丈夫そうです。
また、各手法の比較を行うと「LightGBM >= RF > ロジスティック回帰」といった結果が伺えます。

ちなみに、LightGBMの重要度はこんな感じです。

def display_importances(feature_importance_df_, png_path=None):
    cols = feature_importance_df_[[
        "feature", "importance"
    ]].groupby("feature").mean().sort_values(
        by="importance", ascending=False)[:100].index
    best_features = feature_importance_df_.loc[
        feature_importance_df_.feature.isin(cols)]

    plt.figure(figsize=(12, 15))
    sns.barplot(
        x="importance",
        y="feature",
        data=best_features.sort_values(by="importance", ascending=False))
    plt.title('LightGBM Features (avg over folds)')
    plt.tight_layout()
    if png_path is not None:
        plt.savefig(png_path)

display_importances(feature_importance_df)

f:id:pira_nino:20190714214619p:plain

次にROC曲線 / AUC(c統計量)を求めます。

fig = plt.figure(figsize=[6, 6])
ax = fig.add_subplot(111,aspect = 'equal')
for model in ['lgb', 'lr', 'rf']:
    prob = df['pred_{}'.format(model)]
    fpr, tpr, t = roc_curve(df['cm_dummy'], prob)
    roc_auc = auc(fpr, tpr)
    print('model:{}_auc:{:.6f}'.format(model, roc_auc))
    ax.plot(fpr, tpr, lw=2, alpha=0.3, label='{}_ROC(AUC = {:.6f})'.format(model, roc_auc))
ax.legend(prop={'size': 10})
fig.show()

f:id:pira_nino:20190714215213p:plain

これを見ると、LightGBMとRFはカンスト状態で、ロジスティック回帰はそこそこという感じであることがわかります。
また、書籍のロジスティック回帰の結果がAUC=0.795なので私の実装の結果と大きな差異がなく、概ね実装は間違いなさそうなことも伺えます。

6.2 テスト / コントロールの傾向スコアの可視化

モデルごとにテスト / コントロールの傾向スコアの分布を可視化します。

def visualize_ps(model, png_path=None):
    fig = plt.figure(figsize=[6, 6])
    ax = fig.add_subplot(111)
    bins = np.linspace(0, 1, 100)
    ax.hist(
        df[df['cm_dummy'] == 0]['pred_{}'.format(model)],
        bins=bins,
        color='red',
        width=0.01,
        alpha=0.2,
        label='ctr')
    ax.hist(
        df[df['cm_dummy'] == 1]['pred_{}'.format(model)],
        bins=bins,
        color='blue',
        width=0.01,
        alpha=0.2,
        label='test')
    ax.set_xlim(0, 1.0)
    ax.legend(prop={'size': 10})
    ax.set_title('{}_PS'.format(model), fontdict={'fontsize': 10})
    fig.show()
    if png_path is not None:
        fig.savefig(png_path)


for model in ['lgb', 'lr', 'rf','lgb_calibrated']:
    visualize_ps(model)

f:id:pira_nino:20190714215603p:plainf:id:pira_nino:20190714215607p:plainf:id:pira_nino:20190714215611p:plain

上からロジスティック回帰 / RF / LightGBMの結果です。
傾向スコアの分布は、テスト / コントロールで左右対称になることが理想でした。
RFとLightGBMの分布も左右対称といえば左右対称ですが極端に左右に偏っていることが伺えます。 これだと、追々の傾向スコアマッチングの際に「傾向スコアの近いテスコンのマッチングが望めない」ことは明らかです。
加えて、IPWの際にも「傾向スコアの逆数を用いるのでinfに吹き飛ぶ」ことが懸念されます。 理論上、データ数が無限にあれば例えば、確率0.5っぽいデータが存在し左右に偏ることなく理想的な形になるはずですが、データ数が少なくモデルが賢すぎると極端に左右に偏ることが頻繁に起こります。
ではなぜこのようにテスト / コントロールで極端に左右に偏ったのか？を考えてみましょう。上記結果は、LightGBM / RFといった賢いモデルを用いると、共変量をもとにテスト / コントロールの判別が完璧にできるということを意味しています。これは「テスト / コントロールで似たような背景(共変量)を持つデータは存在しない」ということを示唆します。そもそもマッチングを行うモチベーションは「似たような背景を持つテスト / コントロールの比較」ですので、現状況では似たような背景を持つテスト / コントロールのデータは存在せず、マッチングができないということに繋がります。
このように、賢すぎるモデルを用いるとテスト / コントロールの判別が完璧にできてしまい、マッチングどころでないということが因果推論あるあるです。。。そもそも賢いモデルを用いても、傾向スコアが近しいデータがテスト / コントロールで存在することが理想ですが、こうなってしまったら共変量を考え直すなど問題設定のフェーズに戻る必要があります。
本記事では色々妥協して、ほどほどに当たっていて、傾向スコアの分布もいい感じにマッチングに使えそうな雰囲気を出しているロジスティック回帰の結果を用いて、マッチングを行います。
本当は、RFやLightGBMによる得られる傾向スコアの分布がいい感じになるまで問題設定を見直すのが理想だと思います。

また、モデリングで明らかになった重要度の高い変数に対して、テスト / コントロールの確率の逆数をかけた補正を行い、注目変数の分布のBefore / Afterを確認します。
今回は最も重要度の高かったTVwatch_dayについてロジスティック回帰で得られた傾向スコアの逆数を掛けた補正を行います。

col = 'TVwatch_day'
model = 'lr'
z_col = 'cm_dummy'
ps_col = 'pred_{}'.format(model)

df_test = df.loc[df[z_col] == 1][[col, ps_col]]
df_ctr = df.loc[df[z_col] == 0][[col, ps_col]]
x_test_inv = df_test[col] / df_test[ps_col]
x_ctr_inv = df_ctr[col] / (1 - df_ctr[ps_col])

fig = plt.figure(figsize=[12, 6])
ax = fig.add_subplot(211)
bins = 100
ax.hist(df_ctr[col], bins=bins, color='red', alpha=0.2, label='ctr', density=True)
ax.hist(df_test[col], bins=bins, color='blue', alpha=0.2, label='test', density=True)
ax.legend(prop={'size': 10})
ax.set_title('before_inv', fontdict={'fontsize': 10})
ax.set_xlim(0, 100000)
ax.set_ylim(0, 0.00025)

ax = fig.add_subplot(212)
ax.hist(x_ctr_inv, bins=bins, color='red', alpha=0.2, label='ctr', density=True)
ax.hist(x_test_inv, bins=bins, color='blue', alpha=0.2, label='test', density=True)
ax.legend(prop={'size': 10})
ax.set_title('after_inv', fontdict={'fontsize': 10})
ax.set_xlim(0, 100000)
ax.set_ylim(0, 0.00025)

fig.show()

f:id:pira_nino:20190727155625p:plain

補正後の分布が真っ平らになることが理想です。
上図を見るとbeforeよりはafterの方が平に近づいていることが分かります。

これらの結果を踏まえて後々の分析にとって都合の良いロジスティック回帰を中心に話を進めていきます。*3

6.3Rでの傾向スコアモデリングの実装

さくっと以下のような感じで実装します。

# 読み込み / 設定
df <- read.csv("https://raw.githubusercontent.com/iwanami-datascience/vol3/master/kato%26hoshino/q_data_x.csv")
print(dim(df))
not_use_cols <- c("gamecount", "gamedummy", "gamesecond")
df_model <- df[, !(colnames(df) %in% not_use_cols)]
print(dim(df_model))
# 5CVで学習
k <- 5
folds <- cvFolds(NROW(df), K = k)
df$oof_pred <- rep(0, nrow(df))
for (i in 1:k) {
  print(paste0("CV", i))
  train <- df_model[folds$subsets[folds$which != i], ]
  validation <- df_model[folds$subsets[folds$which == i], ]
  # NOTE: 本来はglmnetなどを用いて、正則化した方が良い
  # NOTE: 書籍では1行後のように、共変量を絞っている
  model <- glm(formula = as.factor(cm_dummy) ~ ., data = train, family=binomial)
  # model <-glm(cm_dummy ~ TVwatch_day + age + sex + marry_dummy + child_dummy + inc + pmoney + area_kanto +area_tokai + area_keihanshin + job_dummy1 + job_dummy2 + job_dummy3 + job_dummy4 + job_dummy5 + job_dummy6 + job_dummy7  + fam_str_dummy1 + fam_str_dummy2 + fam_str_dummy3 + fam_str_dummy4, family=binomial(link="logit") , data = df_model)
  valid_pred <- predict(model, validation, type="response")
  df[folds$subsets[folds$which == i], ]$oof_pred <- valid_pred
}

# loglossとaucの確認
print(LogLoss(df$oof, df$cm_dummy))
oof_train_auc <- prediction(df$oof, df$cm_dummy)
tmp.auc <- performance(oof_train_auc, "auc")
print(as.numeric(tmp.auc@y.values))

本題に戻り上記を回すと、oofのloglossが0.54(pythonでも0.54)、aucが0.79(pythonでも0.79)とpythonでの結果と一致していることが伺えます。
実装に関してですが、ランク落ちした状態で当てはめを行っているため、誤った予測をしている可能性がありますと出ますが今回はスルーしています。勿論ほんとはダメです。
対処法として、変数を絞ってモデリングを行うか、glmnetを使い正則化項を用いることで解消されます。

www.rdocumentation.org

7. 傾向スコアマッチング

7.1 結構面倒なpython

後述のRではサクッとできるライブラリがあるのですが、pythonにはいい感じのライブラリがありません。
ということで、自分はここ最近マイクロソフトが公開したDoWhyライブラリに手を加えました。DoWhyライブラリはモデリング → マッチングを一気通貫してできるのですが、Dagを作る必要がある / モデリング本体を取り出して、上記のような分析ができないという点から今回は本格採用は避けました。
DoWhyの使い方に関しては以下のブログに日本語で分かりやすくまとまっています。

www.krsk-phs.com

DoWhyライブラリを漁ると、class PropensityScoreMatchingEstimator(CausalEstimator)という今回やりたい感じの実装が見つかりましたで、今回の分析に都合がいいように以下のように手を加えます。

def ps_matching(df, y_col, z_col, ps_col, caliper=None):
    treated = df.loc[df[z_col] == 1]
    control = df.loc[df[z_col] == 0]
    control_neighbors = (NearestNeighbors(
        n_neighbors=1,
        algorithm='ball_tree').fit(control[ps_col].values.reshape(-1, 1)))
    distances, indices = control_neighbors.kneighbors(
        treated[ps_col].values.reshape(-1, 1))
    print('num_of_test_{}'.format(len(treated)))
    if caliper is not None:
        less_caliper_index = np.where(distances.reshape(-1) <= caliper)
        distances = distances[less_caliper_index]
        indices = indices[less_caliper_index]

    # estimate ATE on treated by summing over difference between matched neighbors
    ate = 0
    numtreatedunits = len(indices)
    print('num_of_pair_{}'.format(numtreatedunits))

    df_pair_test =  pd.DataFrame(columns=df.columns)
    df_pair_ctr = pd.DataFrame(columns=df.columns)
    for i in range(numtreatedunits):
        treated_outcome = treated.iloc[i][y_col].item()
        control_outcome = control.iloc[indices[i]][y_col].item()
        ate += treated_outcome - control_outcome
        df_pair_ctr = df_pair_ctr.append(control.iloc[indices[i], :])
        df_pair_test = df_pair_test.append(treated.iloc[i, :])
    ate /= numtreatedunits

    # df_pair: マッチしたテスコンを1行に持つDataFrame
    df_pair_ctr.reset_index(inplace=True)
    df_pair_test.reset_index(inplace=True)
    df_pair_ctr.columns = ["ctr_" + colname for colname in df_pair_ctr.columns]
    df_pair_test.columns = ["test_" + colname for colname in df_pair_test.columns]

    df_pair = pd.concat([df_pair_test, df_pair_ctr], axis=1)
    print(ate)
    return ate, df_pair

ate, df_pair = ps_matching(df, y_col='gamesecond',  z_col='cm_dummy', ps_col='pred_lr', caliper=0.2)

キャリパーの設定はTODOになっていたので自分で追加しました。
また、後述のSDの計算の際にどのデータ同士がマッチングしたかのデータが必要なのでdf_pairも出力するようにしました。

出力

num_of_test_4144
num_of_pair_4144
271.882722007722

これにより、CMを見せた時の因果効果は272秒といえます。

7.2 R (Matchingライブラリ)

Rはマッチングライブラリを用いることで簡単に実装できます。

sekhon.berkeley.edu

# 傾向スコアを用いたマッチング
caliper <- 0.2
M <- 1
ties <- FALSE
match <- Match(Y=df$gamesecond, 
               Tr=df$cm_dummy,
               X=df$oof_pred,
               estimand="ATT",
               M=M,
               caliper=caliper,
               ties=ties)
print(summary(match))

出力

Estimate...  181.49 
SE.........  321.06 
T-stat.....  0.56529 
p.val......  0.57188 

Original number of observations..............  10000 
Original number of treated obs...............  4144 
Matched number of observations...............  4144 
Matched number of observations  (unweighted).  4144 

Caliper (SDs)........................................   0.2 
Number of obs dropped by 'exact' or 'caliper'  0

この結果、CMを見せたときの因果効果は 181秒といえます。pythonの結果が271秒であったので多少ブレましたね汗

7.3マッチングの評価

「テストとコントロールのマッチングがイケてるマッチングか？」を定量的に評価する手法として　standardized difference(SD)が知られています。

SD = \frac{\bar{x}_{test} - \bar{x}_{ctr} }{\sqrt{ \frac{s^2_{test} + s^2_{ctr} }{2} }}

SDは、汎用的(universal)な閾値はないけど0.1以下ならば共変量が揃ったマッチングを行えていると解釈できます。
An Introduction to Propensity Score Methods for Reducing the Effects of Confounding in Observational Studies

まずpythonの実装を行います。

def cal_sd(df_pair, colname):
    x_test = df_pair['test_{}'.format(colname)]
    x_ctr = df_pair['ctr_{}'.format(colname)]
    sd = (np.mean(x_test) - np.mean(x_ctr)) / np.sqrt((np.var(x_test) + np.var(x_ctr)) / 2)
    return abs(sd)

ここでは、LightGBMで算出された重要度が高い3つの変数についてSDの計算を行いました。

for colname in ['TVwatch_day', 'inc', 'pmoney']:
    sd = cal_sd(df_pair, colname)
    print('SD_of_{}_is:{}'.format(colname, sd))

出力

SD_of_TVwatch_day_is:0.0866092467165027
SD_of_inc_is:0.038831545007124726
SD_of_pmoney_is:0.04211676096041157

この結果、各変数のSDが0.1を切っているので、傾向スコアマッチングはそれなりに間違っていなさそうであることがわかりました。

次に同様にRでの実装を行います。

# マッチングしたデータのペアのデータフレームの作成
df_tmp<- df
df_tmp$id <- 1:nrow(df_tmp)
df_pair <- cbind(df_tmp[match$index.treated, c('id', colnames(df))],
                 df_tmp[match$index.control, c('id',colnames(df))])
test_cols <- list()
ctr_cols <- list()
for (i in 1 :length(colnames(df))){
  new_name <- paste0("test_", colnames(df)[i])
  test_cols[i] <- new_name
  new_name <- paste0("ctr_", colnames(df)[i])
  ctr_cols[i] <- new_name
}
colnames(df_pair) = c('test_id', test_cols, 'ctr_id', ctr_cols)

# SDの計算
for (col in c("TVwatch_day", "inc", "pmoney")){
  test_colname <- paste0("test_", col)
  ctr_colname <- paste0("ctr_", col)
  mean_test <- mean(df_pair[, test_colname])
  mean_ctr <- mean(df_pair[, ctr_colname])
  var_test <- var(df_pair[, test_colname])
  var_ctr <- var(df_pair[, ctr_colname])
  SD <- abs((mean_test - mean_ctr) /  sqrt((var_test + var_ctr) / 2))
  print(paste0("SD_of_", col, "_is:", SD))
}

出力

[1] "SD_of_TVwatch_day_is:0.0093254842027378"
[1] "SD_of_inc_is:0.0703479485090115"
[1] "SD_of_pmoney_is:0.0548384777267239"

SDの値もpythonと若干ブレた値が算出されましたが、どの変数も0.1を切っているので悪くないマッチングを行えたと解釈できます。

8.IPW

IPWは以下のように実装しました。
それっぽいRライブラリを見つけましたが、これくらい簡単なので実装しちゃいます。

\begin{aligned} W_0 &= \sum_{i=1}^N \frac{1-z_i}{1-e_i} \\ W_1 &= \sum_{i=1}^N \frac{z_i}{e_i} \end{aligned}
\begin{aligned} \hat{E}(y_0) &= \frac{1}{W_0}\sum_{i=1}^N \frac{y_i(1-z_i)}{1-e_i} \\ \hat{E}(y_1) &=\frac{1}{W_1} \sum_{i=1}^N \frac{z_i y_i}{e_i} \end{aligned}

python

Y = df['gamesecond']
Z = df['cm_dummy']

PS = df['pred_lr']

IPW1 = sum(Z * Y / PS) / sum(Z / PS)
IPW0 = sum(((1 - Z) * Y) / (1 - PS)) / sum((1 - Z) / (1 -PS))
effect = IPW1 - IPW0
print('IPW_effect:{}'.format(effect))

R

Y <- df$gamesecond
Z <- df$cm_dummy

PS <- df$oof_pred

IPW1 <-sum(Z * Y / PS) / sum(Z / PS)
IPW0 <- sum(((1 - Z) * Y) / (1 - PS)) / sum((1 - Z) / (1 -PS))
effect <- IPW1 - IPW0
print(paste0("effect", effect))

結果

という感じになりました。Rとpythonでそれなりに値がブレましたね汗
書籍の完全再現はできたので、やってることは間違いでなく変数選択やハイパラの問題だと思います。((詳細の検証はしていませんが、sklearnのLogisticRegressionはデフォルトでL2正則化項がついている等の違いだと思います。 https://scikit-learn.org/stable/modules/generated/sklearn.linear_model.LogisticRegression.html))

9. ダブルロバスト

9.1 回帰モデル

まず、ダブルロバストに用いる回帰モデルを作成します。
z=0のデータのみで学習を行い、割り当てがない場合の目的変数yを予測するg_0(x)を作成します。
同様にz=1のデータのみで割り当てがある場合の目的変数を予測するg_1(x)も作成します。 この2つのモデルを用いて、全データに対して仮に z=0, 1だったときのyを予測を行います。
ここでは、「リッジ回帰」「ランダムフォレスト(RF) 」「LightGBM」の3種類のモデルを作成します。
ここでもCVを切った学習を行い、CVごとに全データへの予測を行い得られたCVごとの予測の平均値を最終予測値とします。

y_col = 'gamesecond'
z_col = 'cm_dummy'
x_cols = [
    col for col in df.columns if col not in [
        'gamecount', 'gamedummy', 'gamesecond', 'cm_dummy', 
        'pred_lgb', 'pred_lr', 'pred_rf', 
        'z0_reg_pred_lgb', 'z1_reg_pred_lgb',
        'z0_reg_pred_rf', 'z1_reg_pred_rf',
        'z0_reg_pred_rg', 'z1_reg_pred_rg'
    ]
]
df_x_all = df[x_cols]
df_y_all = df[y_col]
df_x_all.shape, df_y_all.shape
reg_preds_lgb = np.zeros(df.shape[0])
reg_preds_rf = np.zeros(df.shape[0])
reg_preds_rg = np.zeros(df.shape[0])

for z in [0, 1]:
    print('z={}'.format(z))
    df_tmp = df.loc[df[z_col] == z]
    df_x = df_tmp[x_cols]
    df_y = df_tmp[y_col]
    print(df_x.shape, df_y.shape)

    oof_preds_lgb = np.zeros(df_tmp.shape[0])
    oof_preds_rf = np.zeros(df_tmp.shape[0])
    oof_preds_rg = np.zeros(df_tmp.shape[0])
    insample_mse_list_lgb = []
    insample_mse_list_rf = []
    insample_mse_list_rg = []
    feature_importance_df = pd.DataFrame()
    n_splits = 5
    kf = KFold(n_splits=n_splits, shuffle=True, random_state=0)
    for n_fold, (train_index, valid_index) in enumerate(kf.split(df_x)):
        print('CV_{}'.format(n_fold))
        X_train, y_train = df_x.values[train_index], df_y.values[train_index]
        X_valid, y_valid = df_x.values[valid_index], df_y.values[valid_index]
        lgb_train = lgb.Dataset(X_train, label=y_train, free_raw_data=False)
        lgb_test = lgb.Dataset(X_valid, label=y_valid, free_raw_data=False)

        # LightGBM
        print('LightGBM_fitting...')
        params = {
            'task': 'train',
            'objective': 'regression',
            'metric': 'mse',
        }
        model_lgb = lgb.train(
            params,
            lgb_train,
            valid_sets=[lgb_train, lgb_test],
            valid_names=['train', 'test'],
            num_boost_round=100000,
            early_stopping_rounds=1000,
            verbose_eval=1000)

        # RF
        print('RandomForestRegressor_fitting...')
        model_rf = RandomForestRegressor(n_estimators=500)
        model_rf.fit(X_train, y_train)

        # Ridge
        print('Ridge_fitting...')
        model_rg = Ridge()
        model_rg.fit(X_train, y_train)

        # pred_for_insample
        pred_insample_lgb = model_lgb.predict( X_train, num_iteration=model_lgb.best_iteration)
        pred_insample_rf = model_rf.predict(X_train)
        pred_insample_rg = model_rg.predict(X_train)

        insample_mse_lgb = (mean_squared_error(y_train, pred_insample_lgb))
        insample_mse_rf = (mean_squared_error(y_train, pred_insample_rf))
        insample_mse_rg = (mean_squared_error(y_train, pred_insample_rg))
        print('insample_mse_of_LightGBM:{}'.format(insample_mse_lgb))
        print('insample_mse_of_RandomForestRegressor:{}'.format(insample_mse_rf))
        print('insample_mse_of_Ridge:{}'.format(insample_mse_rg))
        insample_mse_list_lgb.append(insample_mse_lgb)
        insample_mse_list_rf.append(insample_mse_rf)
        insample_mse_list_rg.append(insample_mse_rg)

        # pred_for_oof
        oof_pred_lgb = model_lgb.predict(X_valid, num_iteration=model_lgb.best_iteration)
        oof_pred_rf = model_rf.predict(X_valid)
        oof_pred_rg = model_rg.predict(X_valid)
    
        oof_preds_lgb[valid_index] = oof_pred_lgb
        oof_preds_rf[valid_index] = oof_pred_rf
        oof_preds_rg[valid_index] = oof_pred_rg
        
        oof_mse_lgb = (mean_squared_error(y_valid, oof_pred_lgb))
        oof_mse_rf = (mean_squared_error(y_valid, oof_pred_rf))
        oof_mse_rg = (mean_squared_error(y_valid, oof_pred_rg))
        print('oof_mse_of_LightGBM:{}'.format(oof_mse_lgb))
        print('oof_mse_of_RandomForestRegressor:{}'.format(oof_mse_rf))
        print('oof_mse_of_Ridge:{}'.format(oof_mse_rg))

        #全データへの予測値
        reg_preds_lgb += model_lgb.predict(df_x_all.values, num_iteration=model_lgb.best_iteration)
        reg_preds_rf += model_rf.predict(df_x_all.values)
        reg_preds_rg += model_rg.predict(df_x_all.values)

        # 重要度の算出
        fold_importance_df = pd.DataFrame()
        fold_importance_df["feature"] = x_cols
        fold_importance_df["importance"] = np.log1p(
            model_lgb.feature_importance(
                importance_type='gain', iteration=model_lgb.best_iteration))
        fold_importance_df["fold"] = n_fold + 1
        feature_importance_df = pd.concat(
            [feature_importance_df, fold_importance_df], axis=0)

    print('z={}_all_CV_done'.format(z))
    print('CV_insample_mse_of_LightGBM:{}'.format(np.mean(insample_mse_list_lgb)))
    print('CV_insample_mse_of_RandomForestRegressor:{}'.format(np.mean(insample_mse_list_rf)))
    print('CV_insample_mse_of_Ridge:{}'.format(np.mean(insample_mse_list_rg)))
    
    print('CV_oof_mse_of_LightGBM:{}'.format(mean_squared_error(df_y, oof_preds_lgb)))
    print('CV_oof_mse_of_RandomForestRegressor:{}'.format(mean_squared_error(df_y, oof_preds_rf)))
    print('CV_oof_mse_of_Ridge:{}'.format(mean_squared_error(df_y, oof_preds_rg)))
    #予測値をdfに追加
    df['z{}_reg_pred_lgb'.format(z)] = reg_preds_lgb / n_splits
    df['z{}_reg_pred_rf'.format(z)] = reg_preds_rf / n_splits
    df['z{}_reg_pred_rg'.format(z)] = reg_preds_rg / n_splits

出力

........(略)
z=0_all_CV_done
CV_insample_mse_of_LightGBM:1308.2350804468886
CV_insample_mse_of_RandomForestRegressor:3992318.2751389993
CV_insample_mse_of_Ridge:367626558.77199274
CV_oof_mse_of_LightGBM:797820.3214861164
CV_oof_mse_of_RandomForestRegressor:21480526.954195246
CV_oof_mse_of_Ridge:371428163.63210666

........(略)
z=1_all_CV_done
CV_insample_mse_of_LightGBM:13.106690796586198
CV_insample_mse_of_RandomForestRegressor:28172.902385427773
CV_insample_mse_of_Ridge:227173843.09575287
CV_oof_mse_of_LightGBM:5904.294258803237
CV_oof_mse_of_RandomForestRegressor:96039.45127396613
CV_oof_mse_of_Ridge:232637068.45961764

見やすく表にまとめます。

z 手法 insample_mse oof_mse
0 リッジ回帰 367,626,558 371,428,163
1 リッジ回帰 227,173,843 232,637,068
0 RF 3,992,318 21,480,526
1 RF 28,172 96,039
0 LightGBM 1,308 797,820
1 LightGBM 13 5,904

また、LightGBMによる重要度はこんな感じです。 f:id:pira_nino:20190717205230p:plain

insampleとoofの結果を比較すると、Ridge回帰以外過学習していることが見て取れます。
CVの切り方がイケてないことなど色々思い当たる節はありますが、簡単のためにこの先はリッジ回帰の話に絞ります。

同様にRでも実装を行います。 手間を省くのを目的としてRでは、cv.glmnet()を用いてリッジ回帰の学習を行います。

# 回帰モデル
x_col <- colnames(df)
x_col <- x_col[-which(x_col %in% c("gamesecond", "gamecount", "cm_dummy","gamedummy", "oof_pred"))]
x_col
# z=0, 1でデータセットを分割
df_0 <- df[df$cm_dummy==0, ]
df_1 <- df[df$cm_dummy==1, ]

# 学習
model_0 <- cv.glmnet(y=df_0$gamesecond, x=as.matrix(df_0[x_col]), family ="gaussian", nfolds=5,  alpha=0)
model_1 <- cv.glmnet(y=df_1$gamesecond, x=as.matrix(df_1[x_col]), family ="gaussian", nfolds=5,  alpha=0)

# 0, 1への予測とmseの算出
pred_00 <- predict(model_0, s="lambda.min", newx=as.matrix(df_0[x_col]))
print(mean((df_0$gamesecond - pred_00)^2))
pred_11 <-  predict(model_1, s="lambda.min", newx=as.matrix(df_1[x_col]))
print(mean((df_1$gamesecond - pred_11)^2))

この結果、z=0のデータに対して 368,126,906、z=1のデータに対して227,766,399とpythonの結果と概ね同等の値が得られましたので、問題なさそうです。

9.2ダブルロバスト

数式と照らし合わせて、実装します。

\begin{aligned} \hat{E}^{DR}(y_1) &= \frac{1}{N} \sum_{i=1}^{N}[ \frac{z_{i}y_{i1}}{e(x_i)} + (1-\frac{z_{i}}{e(x_i)}) g_1(x_i) ] \\ &=\frac{1}{N} \sum_{i=1}^{N}[y_{i1} +\frac{z_{i} - e(x)}{e(x)} (y_{i1} - g_1(x_i)) ] \\ \hat{E}^{DR}(y_0) &= \frac{1}{N} \sum_{i=1}^{N} [ \frac{(1-z_i)y_{i0}}{1-e(x_i)} +(1 - \frac{1-z_{i}}{1-e(x_i)} )g_0(x_i) ] \end{aligned} 
Y = df['gamesecond']
Z = df['cm_dummy']
PS = df['pred_lr']
reg_0 = df['z0_reg_pred_rg']
reg_1 = df['z1_reg_pred_rg']

# E(y_1)
dr_1 = np.mean(Z * Y / PS + (1 - Z / PS) * reg_1)

# E(y_0)
dr_0 = np.mean((1 - Z) * Y / (1 - PS) + (1 - (1 - Z) / (1 - PS)) * reg_0)

print('E(y_1):{}_E(y_0):{}_effect:{}'.format(dr_1, dr_0, dr_1 - dr_0))

出力

E(y_1):3102.2829932265536_E(y_0):2762.315547310952_effect:339.96744591560173

これより、ダブルロバストで算出された因果効果は340秒となりました。

同様にRの実装を行います。

# ダブルロバスト
Y <- df$gamesecond
Z <- df$cm_dummy
PS <- df$oof_pred
reg_0 <- predict(model_0, s="lambda.min", newx=as.matrix(df[x_col]))
reg_1 <- predict(model_1, s="lambda.min", newx=as.matrix(df[x_col]))

# E(y_1)
dr_1 <- mean(Z * Y / PS + (1 - Z / PS) * reg_1)
# E(y_0)
dr_0 <- mean((1 - Z) * Y / (1 - PS) + (1 - (1 - Z) / (1 - PS)) * reg_0)

print(paste0("E(y_1):", dr_1, "_E(y_0):", dr_0, "_effect:", dr_1-dr_0))

出力

[1] "E(y_1):3670.7507373371_E(y_0):2744.02213942598_effect:926.728597911123"

因果効果は927秒と算出されました。またしてもpythonの結果とブレましたね。。。。

10.Proximityマッチング

sklearnのRandomForestClassifierにはproximityマトリックスを計算するメソッドが実装されていないので、自作で実装する必要があります。
ここでは、以下の実装を参考にしました。

qiita.com

流れとしては、applyメソッドを用いて各データがどの葉っぱに落ちるかを取得した後に、データどうし同じ葉っぱに落ちた木の数をカウントしていくイメージです。 愚直に実装しているため、すごい量のfor文を回すので計算時間が結構時間がかかります。

# 学習
model_rf = RandomForestClassifier(n_estimators=500)
model_rf.fit(df_x.values, df_z.values)

# proximityマトリックスの算出
def proximity(data):
    prox_matrix = np.zeros((len(data), len(data)))
    n_estimators = len(data[0])

    for e, est in enumerate(np.transpose(np.array(data))):
        for n, n_node in enumerate(est):
            for k, k_node in enumerate(est):
                if n_node == k_node:
                    prox_matrix[n][k] += 1

    prox_matrix = 1.0 * np.array(prox_matrix) / n_estimators

    return prox_matrix

app = model_rf.apply(df_x.values)
prox_matrix = proximity(app)

念のため、loglossを計算すると0.0004であったのでモデリングは正しく行えていそうです。
上のコードを実行するとのデータ数×データ数のサイズのproximityマトリックスが算出されました。
これを用いて、テストデータの行に対して、最も類似度の高いコントロール群のデータのマッチングを行います。

# 対角成分の削除
prox_matrix_remove_diag = prox_matrix - np.diagflat(np.diag(prox_matrix))
# テスト群の行のみ抽出
prox_matrix_remove_diag_test = prox_matrix_remove_diag[df.loc[df[z_col] ==1].index, :]
# コントロール群の列のみ抽出
prox_matrix_remove_diag_test = prox_matrix_remove_diag_test[:, df.loc[df[z_col] == 0].index]
prox_matrix_remove_diag_test.shape
#  各テスト群のデータ'(行)に対して最も類似度の高いコントロール群(列)の抽出
max_index = np.argmax(prox_matrix_remove_diag_test, axis=1)

# 　マッチング
df_pair_test = df.loc[df[z_col] == 1]
df_pair_ctr = df.iloc[max_index, :]
df_pair_ctr.reset_index(inplace=True)
df_pair_test.reset_index(inplace=True)
df_pair_ctr.columns = ["ctr_" + colname for colname in df_pair_ctr.columns]
df_pair_test.columns = ["test_" + colname for colname in df_pair_test.columns]
df_pair = pd.concat([df_pair_test, df_pair_ctr], axis=1)
# 因果効果の算出
print(np.mean(df_pair['test_gamesecond'] - df_pair['ctr_gamesecond']) )

この結果、因果効果は2,478秒と推定されました。
次に、傾向スコアマッチングと同様にSDの算出を行います。

SD_of_TVwatch_day_is:0.7061676939448971
SD_of_inc_is:0.12851003077841508
SD_of_pmoney_is:0.5529597735167537

傾向スコアマッチングに比べ大幅に悪化し、かつSDの値が0.1を大きく超えており、望ましいマッチングができたとは言えません。
(なんでここまで悪い値になったかの原因がわかりません、、、)

さて、Rの実装も行います。

# Proximityマッチング
# モデリング -> Proximityマトリックスの算出
model <- ranger(formula = as.factor(cm_dummy) ~ ., data = df, num.trees = 500)
proximity_matrix <- extract_proximity(model, newdata=df)

# 対角成分を削除して、testのみの行 / ctrのみの列にする
prox_matrix_remove_diag <- proximity_matrix - diag(nrow(proximity_matrix))
test_index = df$cm_dummy == 1
prox_matrix_remove_diag_test <- prox_matrix_remove_diag[test_index, !test_index]
print(dim(prox_matrix_remove_diag_test))

#  各テスト群のデータ(行)に対して最も類似度の高いコントロール群(列)の抽出
max_index <- apply(prox_matrix_remove_diag_test, 1, which.max)
df_pair_test <- df[test_index, ]
df_pair_ctr <- df[max_index, ]
# 因果効果
estimate <- mean(df_pair_test$gamesecond) - mean(df_pair_ctr$gamesecond)
print(estimate)

この結果2,478秒とpythonと一致する結果が得られました。
次にSDを算出します。

# SDの計算
for (col in c("TVwatch_day", "inc", "pmoney")){
  mean_test <- mean(df_pair_test[, col])
  mean_ctr <- mean(df_pair_ctr[, col])
  var_test <- var(df_pair_test[, col])
  var_ctr <- var(df_pair_ctr[, col])
  SD <- abs((mean_test - mean_ctr) /  sqrt((var_test + var_ctr) / 2))
  print(paste0("SD_of_", col, "_is:", SD))
}

出力

[1] "SD_of_TVwatch_day_is:0.292491351380813"
[1] "SD_of_inc_is:0.0604999241643371"
[1] "SD_of_pmoney_is:0.151723869457804"

pythonと同様にTVwatch_dayに関しては0.1を超えたSDが算出されました。
これらの結果より、本データに対してはproximityマッチングはイケてない可能性が高いことが示唆されました。

11.改めて、各種手法による推定因果効果のまとめ

ここまで行なってきた各種手法により算出された因果効果をまとめてみます。

手法 実装 因果効果 SD(TVwatch_day) 備考
傾向スコアマッチング
(ロジスティック回帰) 		python 272 0.087
傾向スコアマッチング
(ロジスティック回帰) 		R 181 0.0093
IPW
(ロジスティック回帰) 		python 470 -
IPW
(ロジスティック回帰) 		R 1,025 -
IPW
(書籍の結果) 		R 1,503 - 書籍と同じ変数を用いた結果の再現
ダブルロバスト
(ロジスティック回帰 /Ridge回帰) 		python 340 -
ダブルロバスト
(ロジスティック回帰 /Ridge回帰) 		R 927 -
Proximityマッチング
(RF) 		python 2,478 0.71 insampleのモデリング
Proximityマッチング
(RF) 		R 2,478 0.29

因果推論あるあるですが、手法により算出される因果効果の値がブレます。
SDが計算できるマッチング系の手法は、結果を信頼していいかを定量的に判断できるものの、IPWやダブルロバストはどこまで信頼していいかを定量的に求める手法がないので評価は難しいです。
また、今回は傾向スコア算出に関しLightGBMやRFではテスト / コントロールを完璧に判別できるという事実がわかりながらも妥協してロジスティック回帰を用いてるので、そもそもの分析設計がイケていない可能性も十分にありえます。
因果推論難しいですね、、、、

12.まとめ

想定以上の長くなってしまいましたが、これでdoneです。
理論編をざっと理解し、上記のコードをコピペすればRubinの因果推論の王道は概ね試せるといった感じになっていると思います。
ただ、「強く無視できる割り付け」という強い制約 / 傾向スコアの分布がいい感じ という大きな壁を超えて実際にに初めて使える手法なので、くれぐれも「安易に因果推論で効果測定という発想をしないでください」という一言を最後にこのブログは終了です(汗

*1:モデル間比較ができなかったのが主な要因です。もしかしたらこれらのパッケージを使う方が最短でできるかもしれません。あくまでも個人の好みの感想です

*2:詳細は書籍を確認してください

*3:LightGBM / RFでテスコンの区別が完璧にできているので本当は妥協しちゃダメだと思います汗