導入・入門から実践までのスライドと実践するためのコードを載せています。因果推論を実務でカジュアルに使えるまで上達させたいものです。
以下、
・統計的因果推論に関するスライド
・資料
・用語整理
・統計的因果推論に関するRの実行例
・書籍情報
からなります。
統計的因果推論に関するスライド
統計的因果推論の学習を進める上でのヒントが記されていると思われるスライドです。カジュアルに回帰分析を行うことに関する注意も書かれています。
因果推論の歴史や利用することのモチベーションが非常にわかりやすく書かれています。
多重共線性を引き合いに出していたりしているのも理解が捗ります。
以下の2つは星野先生の『調査観察データの統計科学』通称赤本の1~3章に関するスライドです。
Rによる実践に関してのスライドです。コードが書かれているので実行してみると良いでしょう。
傾向スコアの計算をRで実践しているスライドです。
資料
こちらの資料は定義が書かれていて助かります。
2016/01/23 (Sat) 第 3 回因果推論を学ぶ会
こちらは論文や教科書の紹介もしています。
プロペンシティスコア（Propensity score; PS）（１）－PSの正しい使い方
カーネギーメロン大学の機械学習における因果推論に関しての講義資料です。
Lecture Notes 17 Causal Inference
Googleのハル・バリアンがまとめたペーパーも良いです。
Causal Inference in Economics and Marketing
用語整理
- ATE（Average Treatment Effect：平均処置効果（因果効果））
- 例で述べるとするならば、母集団全てのユーザーにバナーを見せた場合のブランド名検索による訪問数の平均と、全てのユーザーにバナーを見せなかった場合のブランド名検索による訪問数の平均の差として表されます。バナーを見てしまったユーザーにとって、バナーを見なかったら、という反実仮想なデータは当然存在しないので、そのままでは計算できません。ただし、RCT（Randomized Control Trial）、無作為化比較対照実験、の状況ではバイアスなく推定できるとされています。RCTなケースは限られていると思いますが。
- ATT（Average Treatment Effect on the Treated：処置群における平均処置効果）
- バナーを見せたユーザーにおける、バナーを見せた場合と見せなかった場合の差の期待値。マーケティングにおける施策のROIを計算する際に使うことが望ましいとされています。ATEと同じくRCTにおいてバイアスなく推定できるとされています。
- ATU（Average Treatment Effect on the Untreated：対照群における平均処置効果）
- バナーを見せていないユーザーにおける、バナーを見せた場合と見せなかった場合の差の期待値。マーケティング施策を拡大させるか否かを判断する際に使うことができます。ATEと同じくRCTにおいてバイアスなく推定できるとされています。
- 強い意味での無視可能性
- 共変量に対し求める強い仮定のことで、「バナーを見たか見てないかのバイナリーな変数」や「ブランド名検索での訪問数」などに影響を与えるような共変量に対し、共変量自体で条件をつけて期待値をとると、「バナーを見たか見てないかのバイナリーな変数」と「潜在的なブランド名検索での訪問数」が独立するような特徴が求められています。「バナーを見たか見てないかのバイナリーな変数」が「過去のサイト訪問数（共変量）」や「特定ページへの接触（共変量）」で、配信対象を割り振られている場合は、そのバイナリーな変数は「潜在的なブランド名検索での訪問数」に影響を与えないとされています。
- マッチング
- バナーを見せられたユーザーの持つ、共変量（サイトへの訪問数や、見たページのカテゴリなど）の値と同じ（完全マッチング）、あるいは近い（距離を使ったマッチング）共変量を持っているが、バナーを見せられていない他のユーザーを「同じ人」と見なして、「バナーを見た・見てない」の与える「ブランド名検索での訪問数」への因果効果を推定します。
- 傾向スコア（Propensity score）
- 処置への割り当ての確率。つまり、上述の例でいうところの、バナーを見せられる確率。確率なので、当然0〜1の間の値をとります。推定には2項ロジットモデルが使われているようです。真の傾向スコアを推定できれば、ATE・ATT・ATUを計算することが可能になるそうです。この理屈はベイズの定理より導くことができるようです。詳しくは資料の”第 3 回因果推論を学ぶ会”を見てみてください。
統計的因果推論に関するRの実行例
“Rで学ぶ 傾向スコア解析入門 – 無作為割り当てが出来ない時の因果効果推定”で紹介されていたコードを以下に掲載します。
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 |
#Matchingパッケージを呼び出す library(Matching) #サンプルデータを呼び出す data(lalonde) #データの確認 head(lalonde) # age educ black hisp married nodegr re74 re75 re78 u74 u75 treat # 1 37 11 1 0 1 1 0 0 9930.05 1 1 1 # 2 22 9 0 1 0 1 0 0 3595.89 1 1 1 # 3 30 12 1 0 0 0 0 0 24909.50 1 1 1 # 4 27 11 1 0 0 1 0 0 7506.15 1 1 1 # 5 33 8 1 0 0 1 0 0 289.79 1 1 1 # 6 22 9 1 0 0 1 0 0 4056.49 1 1 1 #ロジスティック回帰モデルを用いて傾向スコアを推定する。 #re78に対する効果を調べたいので、説明変数から除外しておく。 #treatはCMなどを見せたかどうかのバイナリーデータ。 logi <- glm(treat~., data=lalonde[, -9], family = binomial) #因果効果を計算するために、マッチングを推定する。 #今回の効果を見るためのre78をYに指定して、CMなどを見せたかどうかのtreatをTrに指定する。 #共変量としてロジスティック回帰の予測値logi$fittedを用いている。 nsw1 <- Match(Y = lalonde$re78, Tr=lalonde$treat, X =logi$fitted) summary(nsw1) # Estimate... 2138.6 # AI SE...... 797.76 # T-stat..... 2.6807 # p.val...... 0.0073468 # # Original number of observations.............. 445 # Original number of treated obs............... 185 # Matched number of observations............... 185 # Matched number of observations (unweighted). 322 #マッチングのペアを確認する。 lalonde2 <- lalonde lalonde2$id <- 1:nrow(lalonde2) lalonde2$score <- logi$fitted pair.df <- cbind(lalonde2[nsw1$index.treated, c("id","score")], lalonde2[nsw1$index.control, c("id", "score")]) names(pair.df) <- c("t.id", "t.score", "c.id", "c.score") head(pair.df) # t.id t.score c.id c.score # 1 1 0.3927536 357 0.3935865 # 2 2 0.2271642 231 0.2242716 # 3 3 0.5313484 261 0.5313484 # 4 4 0.3285956 254 0.3285956 # 4.1 4 0.3285956 328 0.3286097 # 4.2 4 0.3285956 333 0.3286097 #確かに近そうだ。 #キャリパーマッチング（ペアが特定の距離以上になる時はマッチングしないマッチング） nsw2 <- Match(Y =lalonde2$re78, Tr=lalonde2$treat, X =logi$fitted, caliper = T) summary(nsw2) # Estimate... 2138.6 # AI SE...... 797.76 # T-stat..... 2.6807 # p.val...... 0.0073468 # # Original number of observations.............. 445 # Original number of treated obs............... 185 # Matched number of observations............... 185 # Matched number of observations (unweighted). 322 # # Caliper (SDs)........................................ TRUE # Number of obs dropped by 'exact' or 'caliper' 0 #IPW(Inverse Probability Weight)の計算 ivec1 <- lalonde$treat ivec2 <- rep(1, nrow(lalonde)) - ivec1 ivec <- cbind(ivec1, ivec2) head(ivec) # ivec1 ivec2 # [1,] 1 0 # [2,] 1 0 # [3,] 1 0 # [4,] 1 0 # [5,] 1 0 # [6,] 1 0 iestp1 <- (ivec1/logi$fitted)*(length(ivec1)/sum(ivec1)) iestp2 <- (ivec2/logi$fitted)*(length(ivec2)/sum(ivec2)) iestp <- iestp1 + iestp2 head(iestp) # 1 2 3 4 5 6 # 6.124465 10.588840 4.526984 7.320260 6.053291 6.745951 ipwe <- lm(re78~ivec - 1, weights=iestp, data = lalonde) summary(ipwe) # Call: # lm(formula = re78 ~ ivec - 1, data = lalonde, weights = iestp) # # Weighted Residuals: # Min 1Q Median 3Q Max # -17241 -9985 -4048 6201 129799 # # Coefficients: # Estimate Std. Error t value Pr(>|t|) # ivecivec1 6213.0 462.2 13.44 <2e-16 *** # ivecivec2 4589.4 436.4 10.52 <2e-16 *** # --- # Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 # # Residual standard error: 15090 on 443 degrees of freedom # Multiple R-squared: 0.3967, Adjusted R-squared: 0.3939 # F-statistic: 145.6 on 2 and 443 DF, p-value: < 2.2e-16 #因果効果 causal_effect <- ipwe$coefficients[1]- ipwe$coefficients[2] causal_effect # ivecivec1 # 1623.581 #標準誤差 std_causal_effect <- sqrt(summary(ipwe)$coefficients[3]^2 + summary(ipwe)$coefficients[4]^2) std_causal_effect # [1] 635.7069 |
書籍情報
データ分析の力 因果関係に迫る思考法 (光文社新書)
「原因と結果」の経済学―――データから真実を見抜く思考法
岩波データサイエンス Vol.3
調査観察データの統計科学―因果推論・選択バイアス・データ融合 (シリーズ確率と情報の科学)
