2017-07-30

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 ！そちらのサポートページで公開されているデータを使いました。

github.com

岩波データサイエンス Vol.3

作者: 岩波データサイエンス刊行委員会
出版社/メーカー: 岩波書店
発売日: 2016/06/10
メディア: 単行本（ソフトカバー）
この商品を含むブログ (4件) を見る

　
10000行で35変数ありますが，「CM視聴の有無でゲーム利用回数はどの程度異なるか」をアウトカムにして，23変数を使うことにしました。イメージ的にはこんな感じ（あんまりいい図ではないけど）。

f:id:hikaru1122:20170730181643p:plain 　
上の図は 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列で並べて，それぞれマッチング前とマッチング後のラベルをつけ，共通の凡例を表示させます。

　
f:id:hikaru1122:20170730162121p:plain

　
確認しても，これで共変量が調整されているかなんとも言えないです･･･。なお，この並べ方は次のページを参考にしました。

http://www.sthda.com/english/wiki/ggplot2-easy-way-to-mix-multiple-graphs-on-the-same-page#.V5F9h53_d2E.twitter

因果効果の推定結果を確認

　
参考資料には「あとは回帰分析すればいいぜ」と書いてあったので，やってみます。なんか怪しい。まぁ，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」とのことです。

2016-10-29

統計的因果推論勉強会第5回を開催しました！

統計的因果推論統計分析ノンプログラマーのためのR入門

経営学系統計分析エンドユーザーのため統計的因果推論勉強会第5回を行いました。最近週休3日制を打ち出した某有名大企業の新しいオフィス内での開催で，とても楽しかったです。参加者は7名（クローズド）。

speakerdeck.com

　
今回は宮川本の第5章（の特にバックドア基準）と星野本第4章「共変量選択と無視できない欠測」が勉強範囲でした。目的は2つ。1つは強く無視できる割り当て条件とバックドア基準が似たようなものであることを理解すること。もう1つはRでの実習。Rnotebookを使うと便利でした。配布する資料もカンタンに作れるので。

試用したものはこちらにあります。

github.com

　
Rでの実習は傾向スコアマッチングを星野本の巻末付録に載っているもので写経（使用したデータはいろいろねつ造したフィクションのものです），IPW法を岩波データサイエンス第3巻のサポートページに載っているもので写経しました。

　
c統計量（このcは大文字のほうがいいのか，小文字なのかわかりません･･･）については，私たちはお手軽に rms パッケージの lrm() でやったほうがラクです。「C」と記載されているところがc統計量です。

　
f:id:hikaru1122:20161029205730p:plain

　
星野本第4章までで，論文を読む・書くツールはだいたい出そろっているように思います。しかし，いろいろ謎は残ります。最大の謎は推定した因果効果の大きさ・小ささをどう判断したらよいか，です。効果量（effect size）のようにある程度の基準があれば便利なのですが･･･。

　
さて，最新の『行動計量学』には一般化傾向スコアを用いた論文が出ています。今度はそれもみんなで読んでみたいですね。

2016-08-28

統計的因果推論の勉強会第4回を開催しました！

統計的因果推論

毎月下旬恒例，統計的因果推論勉強会（第4回）を開催しました。今回の範囲は宮川本第4章と星野本第3章「IPW推定量」のところ。スライドはこちらです。

　
speakerdeck.com

　
口頭での説明を加えて進めています（スライドだけではいろいろ不足しているかと）。今日は宮川本と星野本以外に，岩波データサイエンス第3巻も使っています。ちゃんと買って参加してくれる人も嬉しい限り。また，"Causal Inference in Statistice: A Primer" も使いました。

　
宮川本も星野本もどんどん難しくなっている･･･。1ヶ月に1度なので，傾向スコアがなんだとか，反事実がどうだとか，みんなで思い出しながらやりました。やっと傾向スコアのイメージが全体で共有できたように思います。

　
それにしても，IPW推定量はまだまだピンとこない。腹落ちするのに，もう少し時間がかかりそう。そして，数学の勉強も必要･･･。先は長い。

　
今回，メインに使用した本の一覧はこちら。

Causal Inference in Statistics: A Primer

作者: Judea Pearl,Madelyn Glymour,Nicholas P. Jewell
出版社/メーカー: Wiley
発売日: 2016/03/07
メディア: ペーパーバック
この商品を含むブログを見る

岩波データサイエンス Vol.3

作者: 岩波データサイエンス刊行委員会
出版社/メーカー: 岩波書店
発売日: 2016/06/10
メディア: 単行本（ソフトカバー）
この商品を含むブログ (3件) を見る

統計的因果推論―回帰分析の新しい枠組み (シリーズ・予測と発見の科学)

作者: 宮川雅巳
出版社/メーカー: 朝倉書店
発売日: 2004/04
メディア: 単行本
購入: 6人クリック: 67回
この商品を含むブログ (23件) を見る

調査観察データの統計科学―因果推論・選択バイアス・データ融合 (シリーズ確率と情報の科学)

作者: 星野崇宏
出版社/メーカー: 岩波書店
発売日: 2009/07/29
メディア: 単行本
購入: 29人クリック: 285回
この商品を含むブログ (26件) を見る

　
岩崎版『統計的因果推論』もIPWについて，1章分扱っています。やっぱり難しい･･･。

統計的因果推論 (統計解析スタンダード)

作者: 岩崎学
出版社/メーカー: 朝倉書店
発売日: 2015/11/10
メディア: 単行本（ソフトカバー）
この商品を含むブログ (2件) を見る

　
森田『実証分析入門』の16章と17章は統計的因果推論の復習にオススメ。

実証分析入門データから「因果関係」を読み解く作法

作者: 森田果
出版社/メーカー: 日本評論社
発売日: 2014/06/09
メディア: 単行本
この商品を含むブログ (8件) を見る

2016-07-23

統計的因果推論の勉強会第3回を開催しました！

統計的因果推論統計分析参加報告

最近，ブログは月1更新。仕事と学生はなんとか両立できていると信じたい･･･。

　
さて，先日，「経営学系統計学エンドユーザー*1のための統計的因果推論　勉強会」第3回を実施しました。これで3回目。クローズド*2。

　
今回の内容は宮川本と星野本の第3章です。宮川本はパス解析，星野本は傾向スコアについて。星野本は核心に入ってため，難しく，一気に第3章全部はつらいため，60〜69頁をやりました。数式が増えて，本当につらい。

　
宮川本は構造方程式モデリング（SEM）の話を中心に組んでみました。みんなが SEM に慣れているわけではないので3つ実例を見ながら進行。また，みかけの相関と選択バイアスについても復習をしました。

　
ちなみに因果ダイアグラムを用いたみかけの相関と選択バイアスについては『岩波DS Vol.3』の中の「相関と因果と丸と矢印のはなし　はじめてのバックドア基準」が感動的にわかりやすいです。これは次回の勉強会でみんなで読んでみようと思っています。

岩波データサイエンス Vol.3

作者: 岩波データサイエンス刊行委員会
出版社/メーカー: 岩波書店
発売日: 2016/06/10
メディア: 単行本（ソフトカバー）
この商品を含むブログ (3件) を見る

　
スライドはこちらです。

　
speakerdeck.com

　
なお，星野本第3章の詳解は↓に載っています。私にはかなり難しくて，理解が厳しい･･･。 d.hatena.ne.jp

　
また，こちらのブログ「調査観察データにおける因果推論」のシリーズも詳しく解説されています。 smrmkt.hatenablog.jp

　
星野本でやれることはやったと思い，『岩波DS Vol.3』を読み直してみるのがいい気がしてきました。次回は岩波DSを題材にもう少し傾向スコアについて勉強を進めていこうと思います。

*1:主に文系で数学的な内容を学んでいくよりは，実際に社会科学の論文を読み書きするできるようになるのが目的

*2:会場が関係者しか入れないところなので

2016-06-25

統計的因果推論の勉強会第2回を開催しました！

統計的因果推論統計分析

宮川本と星野本の第2章が今回の内容です。毎月それぞれ1章ずつ読み進めています。参加者は私を含めて7人くらい。減ると思って心配していましたが，前回より増えてよかったです。

　
使用したスライドはこちらです。

　
speakerdeck.com

　
第2章はどちらも導入部の終わり，といった感じです。キーワードは

潜在反応モデル（反実仮想モデル）
共変量
傾向スコア

　
でしょうか。例題として，簡単な傾向スコアによる分析もやってみました。R 仲間が増えないかな。我々は文系統計学エンドユーザーなので，数式のところは理解が不足してしまいます。宮川本も星野本も，読んでいて「？」なところは多く出てきます。それでも，統計的因果推論は今後の役に立つかも，ということでコツコツ学んでいます。

　
今回は「強く無視できる割り当て条件（強い意味での無視可能性）」への理解がイマイチでした。どうして傾向スコアで条件付ければ交絡のことを考えなくていいのか，というのがボンヤリしています。

　
きっとこれは星野本3章を読んでいけば解決できるはず，という期待を持ちながら来月に向けて勉強を続けていきます。

　
P.S.
スライドのアップ先を Slideshare から　Speaker Deck に移りました。Slideshare は Mac で作ったPDFで日本がうまく表示されなかったので（解決策はあるみたいですが，面倒ですね）。

2016-05-28

統計的因果推論の勉強会の1回目を開催しました！

統計的因果推論統計分析参加報告

2016年5月28日に統計的因果推論のクローズドな勉強会を開催しました。参加者は5名でした。私を含め全員が非専門家ですので「学び合い」がキーワードです。

　
対象は下のスライドにあるとおり，経営学を学ぶ統計学エンドユーザーです。進行は，私が宮川本と星野本を1章ずつまとめたものを発表し，それについて不明点の質疑応答を行うというようにしました。印象としては，潜在反応モデル／反実仮想モデルが腑に落ちるかどうかが重要のようです。

統計的因果推論勉強会　第1回 from Hikaru GOTO

www.slideshare.net

　
やはり，こういうのは発表する人がいちばん勉強になるような気がします。いろいろと質問を受けて，刺激にもなりました。2冊とも第1章はガイダンスみたいなものなので，本番は翌月の第2章からです。来月もがんばります。

2016-05-12

統計的因果推論の勉強会の前準備

社会人博士書籍紹介読んだ論文統計的因果推論

きっかけは，先月の月1ゼミでした。3時間のゼミのうち，はじめの1時間は輪読をしています。その中で私が「統計的因果推論というものがあるらしい」と情報共有をして，その後「日本社会心理学会　春の方法論セミナー」のページを紹介したところ，先生が興味を示されました。そして，5月からゼミ前の90分間を使って，自由参加で統計的因果推論の勉強会（1回につき1章ずつ）をスタートすることにしました。私が音頭を取って…*1。

　
私を含め，参加者となるのは経営学（主にマネジメント）を勉強・研究しに来ている社会人学生なので，基本，文系が多いです。統計分析の基本知識を一緒に復習しながら，勉強会を進めていく予定です。

　
まず，統計的因果推論勉強会の前準備をするために，資料にあたりました。そのまとめメモとして，書き残しておきます。書籍，ブログ・スライド，Cinii で検索した日本語論文の3タイプに分けます。

書籍

代表的なのは次の2冊。必ず紹介されています。もう紹介も不要なレベル。

統計的因果推論―回帰分析の新しい枠組み (シリーズ・予測と発見の科学)

作者: 宮川雅巳
出版社/メーカー: 朝倉書店
発売日: 2004/04
メディア: 単行本
購入: 6人クリック: 67回
この商品を含むブログ (23件) を見る

調査観察データの統計科学―因果推論・選択バイアス・データ融合 (シリーズ確率と情報の科学)

作者: 星野崇宏
出版社/メーカー: 岩波書店
発売日: 2009/07/29
メディア: 単行本
購入: 29人クリック: 285回
この商品を含むブログ (26件) を見る

　
ただし，数学を学部のときに学んでいない人にはきついです。宮川本は6章より先はサッパリ（目は通した）。現時点では，宮川本・星野本も合わせて3割くらい理解できたかどうか，というところ。私のような文系には次の森田本でイメージをつかむのがよさそうです。16章に説明があります。ただし，アニメが好きな人に限ります。

実証分析入門データから「因果関係」を読み解く作法

作者: 森田果
出版社/メーカー: 日本評論社
発売日: 2014/06/09
メディア: 単行本
この商品を含むブログ (7件) を見る

　
もう1つ日本語ではタイトルど直球の本があります。上記2冊よりはやさしい印象ですが，数式はけっこう出てきます。この本は4割くらい理解できたかもしれません。

統計的因果推論 (統計解析スタンダード)

作者: 岩崎学
出版社/メーカー: 朝倉書店
発売日: 2015/11/11
メディア: 単行本
この商品を含むブログ (2件) を見る

　
もっと文系にやさしく統計的因果推論を説明している本はないかと，洋書もチェックしました。次の2つが読みやすそうでした。Kindle のサンプルをチェック後，まず私は ”Primer” のほうを購入して読み進めています。今年出たばかりだし，著者の1人が Judea Pearl なので，大きなまちがいはないだろう，そして薄い（印刷版だと160ページくらい）というのが選定理由です。とりあえず，1章まではついていけてます。