Proximity Matching による因果効果の推定？
2017年7月30日のKazutan.Rで「Proximity Matching による因果効果の推定」というタイトルで LT をしてきました。スライドは↓ですが，R のスクリプトは少なめなので，その補足説明をここに残しておこうと思います。
https://hikarugoto.github.io/kazutanR20170730/
LT の動機
今回の発表のモチベーションは，たまたま次の論考（何かの予稿集かな？）を見かけ，傾向スコアに代わる何かがあるのだったら，調べてみようかなと思った次第です。ランダムフォレストと傾向スコアについては，＠TJO さんや＠nakamichi さんなどがすでにブログで書いていらっしゃいますが，この論考では Proximity というのを使ってマッチングを行うものなので，ちょっと方法が違うようです。
Zhao, P., Su, X., Ge, T., & Fan, J. (2016). Propensity score and proximity matching using random forest. Contemporary Clinical Trials, 47, 85-92.
ランダムフォレストから Proximity を出して，マッチングを行う例は，すでに “Causal inference with observational data in R: A walk-through from some recent research” という YouTube 動画と資料があったので今回の LT ではそれをなぞってみました。
使用したデータ
統計的因果推論といえば，岩波データサイエンス vol.3 ！ そちらのサポートページで公開されているデータを使いました。
- 作者: 岩波データサイエンス刊行委員会
- 出版社/メーカー: 岩波書店
- 発売日: 2016/06/10
- メディア: 単行本（ソフトカバー）
- この商品を含むブログ (4件) を見る
10000行で35変数ありますが，「CM視聴の有無でゲーム利用回数はどの程度異なるか」をアウトカムにして，23変数を使うことにしました。イメージ的にはこんな感じ（あんまりいい図ではないけど）。
上の図は DiagrammeR で書きました。
grViz(" digraph boxes_and_circles { # a 'graph' statement graph [overlap = true, fontsize = 10] # several 'node' statements node [shape = box, fontname = Helvetica] 原因変数（CM視聴の有無）; ゲーム利用回数 node [shape = circle, width = 1] // sets as circles 共変量 # several 'edge' statements 共変量->{原因変数（CM視聴の有無） ゲーム利用回数} 原因変数（CM視聴の有無）->{ゲーム利用回数} } ")
使用したパッケージ
使ったパッケージは次のとおりです。
library(ggplot2) library(randomForest) library(dplyr) library(magrittr) library(broom) library(knitr) library(DiagrammeR) library(ggpubr)
分析の手順
まず，使用する変数を元のデータから取り出します。select で -（マイナス）を使い，不要な列を削除しました。
CI_data %<>% select(-gamedummy, -area_keihan, -job_dummy8, -fam_str_dummy5, -T:-M3, -gamesecond) CI_data %>% glimpse() Observations: 10,000 Variables: 23 $ cm_dummy <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ... $ area_kanto <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, ... $ area_tokai <int> 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ... $ area_keihanshin <int> 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 1, 1, 0, 0, 0, 1, 0, 0, 0, 0, ... $ age <dbl> 44.5, 34.5, 24.5, 44.5, 34.5, 24.5, 19.0, 54.5, 44.5, 34.5, 34.5, 34.5, ... ＜以下，省略＞
%<>% はあまり見かけませんが，左からデータフレームを右に流していって，その後，同じ名前のオプジェクトに置き換える破壊的な操作をするパイプです。今回は軽い分析なので使ってみました。 2行目はチラ見関数です。str() よりもモダンな感じ？
ランダムフォレストの実行
proximity を得るために，ランダムフォレストを実行します。
cm_dummy = CI_data$cm_dummy==1 # 論理ベクトル rf = randomForest(cm_dummy ~ ., data = CI_data, proximity = TRUE) # 1分くらいかかる
1行目は CM視聴した有無を論理値ベクトルで入手し，それを randamForest で使います。なぜ，論理値ベクトルにする必要があるのか，私には不明です。proximity = T のオプションを忘れないのが重要です。
Proximity でマッチングさせたデータを作る
参考資料では，以下のようにして，マッチングを行っています。
prox = rf$proximity - diag(nrow(rf$proximity)) # 何やってるの？ prox.true = prox[cm_dummy, !cm_dummy] # 何やってるの？ untre.samples = apply(prox.true, 1, which.max) # 行ごとに proximity が最大の列を抜き出し？ tre = CI_data[cm_dummy,] # 処置群のみ取り出し untre = CI_data[!cm_dummy,] # 対照群のデータのみ取り出し untre.samples2 = unique(untre.samples) # 重複している行を削除。 untre2 = untre[untre.samples2,] # 処置群と似ている対照群を抜き出している？ whole.sample.u = rbind(tre, untre2) # それぞれのデータを結合。
もっとも重要なところなのでしょうが，勉強不足でよくわからん･･･。目をつぶって先に進みます。
共変量のバランスは取れたのかチェックしてみる。
せっかくなので，CM見た群と見なかった群の共変量の分布が近づいたかどうか，TVwatch_day と pmoney について見てみます。マッチング前後で見比べるため，ggplot で作図したものを並べます。ggplot を並べる方法はいくつかありますが，今回は ggpubr の ggarrange() を使います。
a1 = CI_data %>% ggplot(aes(x = TVwatch_day)) + geom_density(aes(colour = factor(cm_dummy) ,fill = factor(cm_dummy)), alpha = 0.5) a2 = whole.sample.u %>% ggplot(aes(x = TVwatch_day)) + geom_density(aes(colour = factor(cm_dummy) ,fill = factor(cm_dummy)), alpha = 0.5) b1 = CI_data %>% ggplot(aes(x = pmoney)) + geom_density(aes(colour = factor(cm_dummy) ,fill = factor(cm_dummy)), alpha = 0.5) b2 = whole.sample.u %>% ggplot(aes(x = pmoney)) + geom_density(aes(colour = factor(cm_dummy) ,fill = factor(cm_dummy)), alpha = 0.5)
TVwatch_day のマッチング前の密度分布（a1），マッチング後（a2）として，pmoney も同様にマッチング前後（b1 と b2）を ggplot() で作ります。並べるには次のようにggpub パッケージの ggarrange()を使います。
ggarrange(a1, a2, b1, b2, ncol = 2, nrow = 2, labels = c("Before Matching","After Matching","Before Matching","After Matching"), common.legend = TRUE)
図を2行2列で並べて，それぞれマッチング前とマッチング後のラベルをつけ，共通の凡例を表示させます。
確認しても，これで共変量が調整されているかなんとも言えないです･･･。なお，この並べ方は次のページを参考にしました。
因果効果の推定結果を確認
参考資料には「あとは回帰分析すればいいぜ」と書いてあったので，やってみます。なんか怪しい。まぁ，RCT に近づいたからやっていいという考え方だとは思うのですが･･･。
lm(gamecount ~ ., data = whole.sample.u) %>% tidy() %>% head(3) %>% kable()
lm() で得た結果を2行目でデータフレームにして，3行目ではじめの数行だけ表示させ，4行目で下の表を出すためのマークダウンを得ました。
| term | estimate | std.error | statistic | p.value |
|---|---|---|---|---|
| (Intercept) | 46.467720 | 7.462751 | 6.226621 | 0.0000000 |
| cm_dummy | 5.080259 | 1.499469 | 3.388039 | 0.0007088 |
| area_kanto | -3.428637 | 2.041777 | -1.679242 | 0.0931606 |
上の表から，CMを見た場合の平均処置効果（ATE）は 5.08 であることがわかりました。岩波DS vol.3 でIPWによってATEを出したときの結果が 5.32 ですから（96ページ），近いと言えば近い*1。
| ATE | |
|---|---|
| 岩波DS vol.3 | 5.32 |
| RF で IPW @nakamichi*2 | 5.32 |
| Proximity でマッチング | 5.08 |
分析してみた感想
やっていることをよく理解できず，今回の結果を信じるには至っていません。ちなみに傾向スコアマッチングは使わないほうがいいという論考が出ています*3。今回の方法は傾向スコアの弱点をクリアしているものの*4，IPW法など他のやり方があるので，わざわざマッチングをしなくてもいいように思いました。
最後に
やはり，Kazutan.R に参加してよかったです。主催者 kazutan 先生やハンズオンの MrUnadadon氏，LTされた方々に感謝。私にはデータ分析をする人が身近にいないので，情報の入手先は twitter やこういう勉強会がすべてなので（書籍の情報も twitter で知ることが多い），助かります。
*1:ちなみに，スライドで2.32と書いたのは ATT でまちがいでした。スミマセン
*2:https://rpubs.com/nakamichi_/study-iwanami-ds3
*3:King and Nielsen （2016）"Why Propensity Scores Should Not Be Used for Matching"
*4:Zhao らによれば，傾向スコアの限界は「model misspecifications, categorical variables with more than two levels, difficulties in handling missing data, and nonlinear relationships」とのことです。
