岩波データサイエンスvol.3
最近、業務の中で統計的因果推論の手法を適用するにあたって勉強をしており、その一環として各手法や実装方法の理解を深めるために巷に出回っているデータを用いて実験的にコードを書いて試してみたりしています。
本稿はまさにタイトルの通りなのですが、そういった勉強を進める中で『岩波データサイエンスvol.3』で取り上げられているCM接触の因果効果について、本文中で紹介されている方法とは別のやり方でいくつか因果効果の推定を行ってみた、という記事になります*1。具体的には、本文中では主にIPW推定量を用いた回帰分析によって因果効果の推定を行っていますが、本ブログでは「傾向スコアマッチング」と「Causal Forest」の両手法を同じデータに適用してみるとどうなるかを試してみたい思います。まずは紙幅の関係上、本稿では「傾向スコアマッチング」に関して取り上げられればと思います。
データの準備
本稿で扱うデータは上述の『岩波データサイエンスvol.3』中の加藤・星野の章で取り上げられた、CM接触のアプリゲームプレイへの因果効果を分析する際のデータを用います。データの詳細に関してはここでは割愛しますが、詳細を把握されたい方は上記の章を確認してください。
また、傾向スコアの推定に用いる共変量や Causal Forest の説明変数に用いるデータは上記加藤・星野の章で取り上げられたものと同じものを使用することとします。
手法の説明
「傾向スコアマッチング」と「Causal Forest」に関しては、以下先達のまとめが非常にわかりやすくて参考になります。
より本格的に参考にしたい場合は、傾向スコアであれば G. Imbens & D. Rubin (2015)、Causal Forest であれば S. Wager & S. Athey (2017) あたりが理論も含めてしっかり押さえることができるかと思います。
因果効果の推定
以下、傾向スコアマッチングを用いて因果効果の推定を行っていこうと思います。分析に際して、今回は R言語を用いることにします。
まずは、上述のサンプルデータをGithubからダウンロードしてくることにします。
> library(devtools) > filepath <- "https://raw.githubusercontent.com/iwanami-datascience/vol3/master/kato%26hoshino/q_data_x.csv" > df <- read.csv(url(filepath))
傾向スコアマッチング
傾向スコアマッチングを行うに際して、今回は Matching と cobalt ライブラリを主に使用していきたいと思います。前者が傾向スコアマッチング関連のライブラリで*2、後者はマッチング後の共変量のバランスの可視化や確認に便利な機能を揃えています*3。
最初にロジスティック回帰を行って傾向スコアの推定を行いますが、ここは特段オリジナルの加藤・星野のものと変わるところはありません。
# ロジスティック回帰による傾向スコア推定 > ps_mod <- glm( formula = 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, data = df, family = binomial(link = "logit") ) > ps <- ps_mod$fitted.values
次に、Matching ライブラリ中の Match 関数を用いて実際に傾向スコアマッチングを行っていこうと思います。推定量は ATT(Average Treatment Effect on the Treated = 処置群に対する平均処置効果)とします。
> ps_match <- Match( Y=df$gamesecond, Tr=df$cm_dummy, X=ps, estimand="ATT", M=1, caliper=TRUE, ties=FALSE, replace=TRUE )
また、マッチング後の共変量のバランスを見る際には cobalt ライブラリ中の bal.tab 関数を用い、それを love.plot 関数で可視化すると便利です。
> bal_m <- bal.tab( ps_match, 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, data=df ) > love.plot(bal_m, starts="std", binary="std", grid=F, limits=c(0,1), threshold=.1)
するとこのような形で共変量の調整前後のバランスを簡単に可視化することができます。
今回はダミー変数が多く、調整後のほうが調整前よりもバランスが悪くなっている変数もありますが、調整によってどの共変量も概ね両群でバランスがとれていると判断し、ATTの確認を行うこととします。
> summary(ps_match) Estimate... 147.54 SE......... 320.4 T-stat..... 0.46047 p.val...... 0.64518 Original number of observations.............. 10000 Original number of treated obs............... 4144 Matched number of observations............... 4144 Matched number of observations (unweighted). 4144 Caliper (SDs)........................................ TRUE Number of obs dropped by 'exact' or 'caliper' 0
Estimate の値が ATT の推定量になります。今回だと、CM接触の因果効果は 147.54 (秒) と推定されています。結果としては『岩波データサイエンス』の中で述べられている推定値よりも低い値が出てしまっているようです。
ちなみに上記を何度か試行錯誤する中で分かったのですが、上記 Match 関数で ties = FALSE とした場合、マッチング対象として選ばれる組み合わせが変化するらしく、毎回同じような結果が出るわけではないようです。これはある種マッチングの実行時に選択バイアスが発生しているとも言える訳で、それをそのまま放置することはよくないため、ATTの推定量を点ではなく「ばらつき」で考えることにします。具体的には、単純ですが上記のマッチングを反復処理し、得られた推定量の要約統計量を見ることで得られた推定量がどれくらいばらつきうるのかを確認することとします。
コード部分は単純なので省略して結果だけ記載します。以下、mean_boot はマッチング処理を200回繰り返した際にえられた推定量を格納した配列です。
> mean(mean_boot) [1] 192.7477 > quantile(mean_boot, probs=c(.05, .25, .5, .75, .95)) 5% 25% 50% 75% 95% -38.21432 108.72116 186.61245 289.68641 410.14094
こうして結果を見てみると興味深いのが、オリジナルの加藤・星野の章で記載されていた IPW推定量によるATTの推定値との差分がそこそこ大きいのと、マッチングによって得られたATTのばらつきを見た際に95%信頼区間が0をまたがって正負の双方に分布しているといったあたりでしょうか。この結果はATTの推定量を点ではなくばらつきとして捉えたからこそ分かったことでもあります。そもそも手法が違うので結果を単純比較することはできませんが、ATT推定量のばらつきとしては一応オリジナルの加藤・星野の結果である 397.8 (秒) も 95%信頼区間の中に入っていると言えば入っているようです。
ただ、後者に関してはオリジナルの結果とは異なり、CM接触による因果効果は「必ずしもあるとは言えない」とも言えるのかもしれません。このあたり、自分の解釈ややり方が悪かったりするかもしれないのですが、観察データから因果効果を推定するにあたって難しいところだなと思った次第でした。
