はじめに

岩波データサイエンス3巻に，加藤・星野「因果効果推定の応用」という解説がある．

因果推論の文献では，具体例が医学や社会科学から取られていることが多く， ビジネスで使ってみたいと思っている身としてはイマイチ実感が湧かなかったのだが， この解説はTVCMがスマホアプリ利用に与える因果効果を推定するという ビジネス上の話題を扱っており，興味深い． また，分析に使ったデータが公開されているので自分で分析してみることもできる．

ここでは，上記解説の分析を再現してみるとともに， 解説では使われていない因果推論の手法も同データで試してみることで， 因果推論への理解を深めたい．

問題設定

岩波DS3解説中で分析に使われているデータはgithubにある． データの項目定義はサポートサイトの補論に記載されているが，再掲しておく．

library(dplyr)
library(tidyr)
library(broom)
df <- read.csv("data/q_data_x.csv")
変数名 内容
cm_dummy CM接触
gamedummy アプリ利用
gamesecond アプリ利用秒数
gamecount アプリ利用回数
area_kanto 居住地が関東
area_keihan 居住地が京浜
area_tokai 居住地が東海
area_keihanshin 居住地が京阪神
age 年齢
sex 性別(男性=1)
marry_dummy 結婚の有無(既婚=1)
child_dummy 子供の有無(子供あり=1)
job_dummy1 職業:正社員・公務員
job_dummy2 職業:自営業・個人事業主
job_dummy3 職業:派遣社員・契約社員
job_dummy4 職業:その他職業
job_dummy5 職業:パート・アルバイト
job_dummy6 職業:専業主婦・専業主夫
job_dummy7 職業:学生
job_dummy8 職業:無職
inc 年間収入
pmoney 1月あたりのお小遣い
fam_str_dummy1 家族構成:単身
fam_str_dummy2 家族構成:夫婦のみ
fam_str_dummy3 家族構成:2世代同居
fam_str_dummy4 家族構成:3世代同居
fam_str_dummy5 家族構成:その他
TVwatch_day 1日当たりの平均TV視聴時間
T 13~19歳の男女
F1 20~34歳の女性
F2 35~49歳の女性
F3 50歳以上の女性
M1 20~34歳の男性
M2 35~49歳の男性
M3 50歳以上の男性

分析1：岩波DS3解説の再現

まずは，解説記事の分析結果を再現してみる． サンプルコードはgithubで提供されているが，なるべく丸写しにならないようにする．

分析の目的は以下の通りである．

平均値の比較

最初に，CM接触群とCM非接触群の平均的なアプリ利用状況を確認しておく．

df %>%
  group_by(cm_dummy) %>% 
  summarise(mean(gamedummy), mean(gamecount), mean(gamesecond)) %>% 
  knitr::kable(digits=3)
cm_dummy mean(gamedummy) mean(gamecount) mean(gamesecond)
0 0.073 10.048 3107.706
1 0.075 8.564 2478.066

単純な平均値の比較では，アプリ利用率はCM接触群の方がわずかに高いものの， 利用秒数と利用回数の平均値はCM接触群の方が小さいことがわかる．

セレクションバイアス

CMを見ている方がアプリを使わないなんておかしいなということで， セレクションバイアスを疑うことになる． そこで，各群における共変量の分布（平均値と標準偏差）を確認しておく．

inner_join(
  df %>% 
    group_by(cm_dummy) %>% 
    summarise_each(funs(mean), -starts_with("game")) %>% 
    gather(var, mean, -cm_dummy) %>%
    spread(cm_dummy, mean, sep="_mean_"),
  df %>% 
    group_by(cm_dummy) %>% 
    summarise_each(funs(sd), -starts_with("game")) %>% 
    gather(var, sd, -cm_dummy) %>%
    spread(cm_dummy, sd, sep="_sd_"),
  by="var"
) %>% select(var, ends_with("0"), ends_with("1")) %>% 
knitr::kable(digits=3)
var cm_dummy_mean_0 cm_dummy_sd_0 cm_dummy_mean_1 cm_dummy_sd_1
age 40.187 10.636 41.767 10.148
area_kanto 0.063 0.243 0.131 0.337
area_keihan 0.509 0.500 0.701 0.458
area_keihanshin 0.303 0.460 0.075 0.263
area_tokai 0.124 0.330 0.093 0.291
child_dummy 0.422 0.494 0.424 0.494
F1 0.138 0.345 0.113 0.317
F2 0.136 0.343 0.226 0.418
F3 0.050 0.218 0.055 0.229
fam_str_dummy1 0.150 0.358 0.145 0.352
fam_str_dummy2 0.128 0.334 0.168 0.374
fam_str_dummy3 0.619 0.486 0.622 0.485
fam_str_dummy4 0.084 0.277 0.051 0.219
fam_str_dummy5 0.019 0.136 0.014 0.117
inc 369.246 264.487 341.697 270.695
job_dummy1 0.601 0.490 0.518 0.500
job_dummy2 0.057 0.232 0.050 0.218
job_dummy3 0.069 0.253 0.086 0.280
job_dummy4 0.011 0.106 0.014 0.115
job_dummy5 0.091 0.287 0.156 0.363
job_dummy6 0.092 0.289 0.111 0.314
job_dummy7 0.046 0.210 0.031 0.172
job_dummy8 0.033 0.179 0.035 0.184
M1 0.167 0.373 0.103 0.304
M2 0.342 0.474 0.311 0.463
M3 0.154 0.361 0.180 0.384
marry_dummy 0.630 0.483 0.670 0.470
pmoney 3.558 3.360 3.544 3.403
sex 0.671 0.470 0.597 0.491
T 0.013 0.115 0.013 0.112
TVwatch_day 5714.982 5690.371 11461.881 8851.091

居住地，性別，年齢，職業，TV視聴時間といった項目で， CM接触群と非接触群の間にかなり違いがあることが見てとれる．

平均処置効果の推定

セレクションバイアスを補正して因果効果を推定する． 解説に従い，傾向スコアをロジスティック回帰で推定してIPW推定量を計算する．

参考のため，平均処置効果（ATE, Average Treatment Effect）と 処置群での平均処置効果（ATT, Average Treatment effect on the Treated）の IPW推定量の式を書いておく．

\[ ATE = \frac{\sum_{i=1}^{N}\frac{z_i y_i}{e(x_i)}}{\sum_{i=1}^{N}\frac{z_i}{e(x_i)}} - \frac{\sum_{i=1}^{N}\frac{(1 - z_i) y_i}{1 - e(x_i)}}{\sum_{i=1}^{N}\frac{(1 - z_i)}{1 - e(x_i)}} \\ ATT = \frac{\sum_{i=1}^{N} z_i y_i}{\sum_{i=1}^{N} z_i} - \frac{\sum_{i=1}^{N}\frac{(1 - z_i) e_i y_i}{1 - e(x_i)}}{\sum_{i=1}^{N}\frac{(1 - z_i)e_i}{1 - e(x_i)}} \]

# 傾向スコア推定
ps_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,
  data=df, family=binomial()
)
df$ps <- ps_model$fitted.values

# 傾向スコアの逆数の重み
## ATE
df$ipw_ate  <- with(
  df, ifelse(cm_dummy == 1, 1 / ps,  1 / (1 - ps))
)
## ATT
df$ipw_att  <- with(
  df, ifelse(cm_dummy == 1, 1, ps / (1 - ps))
)

# IPW推定値
ipwe_models <- list(
  ## ATE
  ATE_gamedummy  = lm(gamedummy  ~ cm_dummy + I(1-cm_dummy) + 0, weights=ipw_ate, data=df),
  ATE_gamecount  = lm(gamecount  ~ cm_dummy + I(1-cm_dummy) + 0, weights=ipw_ate, data=df),
  ATE_gamesecond = lm(gamesecond ~ cm_dummy + I(1-cm_dummy) + 0, weights=ipw_ate, data=df),
  ## ATT
  ATT_gamedummy  = lm(gamedummy  ~ cm_dummy + I(1-cm_dummy) + 0, weights=ipw_att, data=df),
  ATT_gamecount  = lm(gamecount  ~ cm_dummy + I(1-cm_dummy) + 0, weights=ipw_att, data=df),
  ATT_gamesecond = lm(gamesecond ~ cm_dummy + I(1-cm_dummy) + 0, weights=ipw_att, data=df)
)
lapply(ipwe_models, tidy) %>%
  bind_rows(.id="outcome")%>%
  select(outcome, term, estimate) %>% 
  spread(term, estimate) %>%
  mutate(causal_effect = cm_dummy - `I(1 - cm_dummy)`) %>% 
  knitr::kable(digits=3)
outcome cm_dummy I(1 - cm_dummy) causal_effect
ATE_gamecount 13.636 8.317 5.320
ATE_gamedummy 0.094 0.062 0.032
ATE_gamesecond 4143.330 2639.412 1503.918
ATT_gamecount 8.564 6.249 2.315
ATT_gamedummy 0.075 0.050 0.026
ATT_gamesecond 2478.066 2080.266 397.800

IPW推定の結果，アプリの利用有無・利用回数・利用秒数のすべてにおいて， ATE，ATTのどちらも正の因果効果があることがわかる． また，ATTの方がATEよりも小さくなっている．

調整効果の推定

個人の属性等による調整効果を推定する． IPWのウェイトを付けて回帰するだけなので簡単だ．

# アプリ利用秒数への調整効果
ipwe_mod_models <- list(
  gamesecond_treated = lm(gamesecond ~ child_dummy + area_kanto +area_tokai +
                             area_keihanshin + T + F1 + F2 + F3 + M1 + M2,
                           weights=ipw_ate, data=df[df$cm_dummy==1,]),
  gamesecond_untreated = lm(gamesecond ~ child_dummy + area_kanto +area_tokai +
                               area_keihanshin + T + F1 + F2 + F3 + M1 + M2,
                             weights=ipw_ate, data=df[df$cm_dummy==0,])
)
lapply(ipwe_mod_models, tidy) %>%
  bind_rows(.id="outcome")%>%
  knitr::kable(digits=1)
outcome term estimate std.error statistic p.value
gamesecond_treated (Intercept) 498.6 819.6 0.6 0.5
gamesecond_treated child_dummy -1938.9 666.3 -2.9 0.0
gamesecond_treated area_kanto -2729.2 1083.7 -2.5 0.0
gamesecond_treated area_tokai -3030.9 996.9 -3.0 0.0
gamesecond_treated area_keihanshin 7833.9 807.3 9.7 0.0
gamesecond_treated T 6413.3 2904.3 2.2 0.0
gamesecond_treated F1 -1042.6 1221.1 -0.9 0.4
gamesecond_treated F2 1142.6 1047.0 1.1 0.3
gamesecond_treated F3 -265.1 1539.0 -0.2 0.9
gamesecond_treated M1 1247.3 1122.2 1.1 0.3
gamesecond_treated M2 9624.8 929.9 10.4 0.0
gamesecond_untreated (Intercept) 631.9 596.7 1.1 0.3
gamesecond_untreated child_dummy 1827.3 511.6 3.6 0.0
gamesecond_untreated area_kanto -28.3 846.7 0.0 1.0
gamesecond_untreated area_tokai 41.5 646.1 0.1 0.9
gamesecond_untreated area_keihanshin -1091.4 607.3 -1.8 0.1
gamesecond_untreated T 7141.7 2363.3 3.0 0.0
gamesecond_untreated F1 382.8 858.5 0.4 0.7
gamesecond_untreated F2 944.7 874.0 1.1 0.3
gamesecond_untreated F3 -282.5 1035.3 -0.3 0.8
gamesecond_untreated M1 7164.8 847.5 8.5 0.0
gamesecond_untreated M2 981.1 699.3 1.4 0.2

子供の有無や居住地，性別，年齢の効果がけっこうリアリティが感じられる結果に なっており，おもしろい．

分析2：本に載っていないことを試す

ここから先は，[星野2009]（赤本）を参考にしつつ，岩波DS3には載っていない分析方法を試してみたい．

傾向スコアによる層別解析

傾向スコアによる層別で平均因果効果を推定してみる． 層の数をいくつに設定するのがよいのか分からないので，とりあえず5, 10, 20層の 3パターンを試してみることにする．

# 層別する（層ごとのデータ数が等しくなるようにする）
df$strata5  <- cut(df$ps, c(0, quantile(df$ps, seq(0,1,0.2))[2:5], 1))
df$strata10 <- cut(df$ps, c(0, quantile(df$ps, seq(0, 1, 0.1))[2:10], 1))
df$strata20 <- cut(df$ps, c(0, quantile(df$ps, seq(0, 1, 0.05))[2:20], 1))

# 平均因果効果
te_strata <- list()
for(strata in c("strata5", "strata10", "strata20")) {
  te_strata[[strata]] <- df %>% 
    group_by_(strata) %>% 
    summarise(te_gamedummy  = sum(cm_dummy*gamedummy)  / sum(cm_dummy) -
                sum((1-cm_dummy)*gamedummy) / sum(1-cm_dummy),
              te_gamecount  = sum(cm_dummy*gamecount)  / sum(cm_dummy) - 
                sum((1-cm_dummy)*gamecount) / sum(1-cm_dummy),
              te_gamesecond = sum(cm_dummy*gamesecond) / sum(cm_dummy) -
                sum((1-cm_dummy)*gamesecond) / sum(1-cm_dummy),
              n=n()) %>% 
    summarise(ate_gamedummy  = sum(te_gamedummy*n) / sum(n),
              ate_gamecount  = sum(te_gamecount*n) / sum(n),
              ate_gamesecond = sum(te_gamesecond*n) / sum(n))
}
knitr::kable(te_strata$strata5,  digits=3)
ate_gamedummy ate_gamecount ate_gamesecond
0.022 1.499 68.263 
knitr::kable(te_strata$strata10, digits=3)
ate_gamedummy ate_gamecount ate_gamesecond
0.03 4.564 1317.082 
knitr::kable(te_strata$strata20, digits=3)
ate_gamedummy ate_gamecount ate_gamesecond
0.056 16.048 5482.291

層の数によってかなり異なる結果になってしまった． 10層の場合がIPW推定値に最も近いが，だからといってこれが正しいのかどうか分からない．

おそらく実際に使う場合には共変量がバランスしているかどうかなど， きちんとチェックして適切な層数を決める必要がありそうだが， とりあえず色々な手法を試してみるのが目的なので，次に進む．

傾向スコアによるマッチング

マッチングはMatchingパッケージで簡単に実行できる．

library(Matching)
match_result <- with(df, list(
  gamedummy  = Match(Y=gamedummy,  Tr=cm_dummy, X=ps),
  gamecount  = Match(Y=gamecount,  Tr=cm_dummy, X=ps),
  gamesecond = Match(Y=gamesecond, Tr=cm_dummy, X=ps)
))

summary(match_result$gamedummy)
## 
## Estimate...  0.020597 
## AI SE......  0.0093214 
## T-stat.....  2.2097 
## p.val......  0.027129 
## 
## Original number of observations..............  10000 
## Original number of treated obs...............  4144 
## Matched number of observations...............  4144 
## Matched number of observations  (unweighted).  53955
summary(match_result$gamecount)
## 
## Estimate...  1.6133 
## AI SE......  1.7147 
## T-stat.....  0.94086 
## p.val......  0.34678 
## 
## Original number of observations..............  10000 
## Original number of treated obs...............  4144 
## Matched number of observations...............  4144 
## Matched number of observations  (unweighted).  53955
summary(match_result$gamesecond)
## 
## Estimate...  179.44 
## AI SE......  551.51 
## T-stat.....  0.32535 
## p.val......  0.74491 
## 
## Original number of observations..............  10000 
## Original number of treated obs...............  4144 
## Matched number of observations...............  4144 
## Matched number of observations  (unweighted).  53955

結果は，gamedummyとgamesecondについてはIPW推定値と桁は合っている． gamesecondは一桁ずれている． だがgamecountとgamesecondは誤差が大きく，そもそも因果効果が有意でないようだ．

この例に限らず，マーケティング施策というのは，処置群でも非処置群でも 結局ほとんどの人は望ましいアクションを取ってくれないというデータになることが 多いので，効果推定のばらつきが大きいのは納得できるように思う． （ただ，それだと因果効果がんばって推定してもなあという気もしてくるが……．）

# CMを見てもほとんどの人はゲームをやらない
data.frame(
  treated   = quantile(df$gamecount[df$cm_dummy==0], seq(0,1,0.1)),
  untreated = quantile(df$gamecount[df$cm_dummy==1], seq(0,1,0.1))
) %>% t() %>% knitr::kable()
0% 10% 20% 30% 40% 50% 60% 70% 80% 90% 100%
treated 0 0 0 0 0 0 0 0 0 0 765
untreated 0 0 0 0 0 0 0 0 0 0 688

傾向スコアによるカーネルマッチング

# バンド幅（赤本に載っているおまじないの値）
bw1 <- 1.06 * nrow(df[df$cm_dummy==1,])^(-0.2) * sd(df$ps[df$cm_dummy==1])
bw0 <- 1.06 * nrow(df[df$cm_dummy==0,])^(-0.2) * sd(df$ps[df$cm_dummy==0])

km_result <- list()
for(col in c("gamedummy", "gamecount", "gamesecond")) {
  # カーネル回帰
  ## CM非接触群の、仮にCMに接触していた場合の結果変数
  y1_est <- ksmooth(x=df$ps[df$cm_dummy==1], y=df[df$cm_dummy==1, col],
                    x.points=df$ps[df$cm_dummy==0], bandwidth=bw1, kernel="normal")
  ## CM接触群の、仮にCMに接触していなかった場合の結果変数
  y0_est <- ksmooth(x=df$ps[df$cm_dummy==0], y=df[df$cm_dummy==0, col],
                    x.points=df$ps[df$cm_dummy==1], bandwidth=bw0, kernel="normal")

  # 平均因果効果
  ate <- structure(
    mean(c(df[df$cm_dummy==1, col], y1_est$y) - c(y0_est$y, df[df$cm_dummy==0, col])),
    names=col
  )
  km_result[col] <- ate
}

knitr::kable(as.data.frame(km_result), digits=3)
gamedummy gamecount gamesecond
0.038 9.251 2838.06

これもまあ，桁はIPW推定と合っている．

傾向スコアをロジスティック回帰以外の方法で推定する

傾向スコアを精度よく推定するのが大事というなら， ロジスティック回帰よりもっと複雑な機械学習モデルを使う方が良いのでは， という疑問が湧いてくるのが人情だと思うので，実際に試してみる．

モデルは何でもよいのだが，ここでは適当に学習させてもそこそこ良い精度が出そうな ランダムフォレストとGBMを使う．

library("randomForest")
library("gbm")

# ランダムフォレスト（共変量はロジスティック回帰と同じ，パラメータは適当）
ps_rf <- randomForest(as.factor(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)
ps_rf$fitted.values <- ps_rf$votes[,"1"]

# GBM（共変量はロジスティック回帰と同じ，パラメータは適当）
ps_gbm <- gbm(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, distribution = "bernoulli", n.trees=1000, interaction.depth = 5)
ps_gbm$fitted.values <- 1 / (1 + exp(-ps_gbm$fit))

# 傾向スコア推定のAUC（c統計量）確認
c(pROC::auc(df$cm_dummy, ps_model$fitted.values),
  pROC::auc(df$cm_dummy, ps_rf$fitted.values),
  pROC::auc(df$cm_dummy, ps_gbm$fitted.values))
## [1] 0.7917029 0.9936573 0.7996351

ものの本によるとAUCが0.8というのが精度の目安らしい．

ロジスティック回帰とGBMの精度は大体0.8で，GBMの方が少し高い． ランダムフォレストの精度は0.99と非常に高い． オーバーフィッティングだと思うのだが，傾向スコア推定は オーバーフィッティングしていてもよいのか悪いのか，分からない．

とりあえず，IPW推定量を計算してみる．

# IPW推定量
df$ps_rf <- ps_rf$fitted.values
df$ipw_rf  <- with(df, ifelse(cm_dummy == 1, 1 / ps_rf,  1 / (1 - ps_rf)))
df$ps_gbm <- ps_gbm$fitted.values
df$ipw_gbm  <- with(df, ifelse(cm_dummy == 1, 1 / ps_gbm,  1 / (1 - ps_gbm)))

ipwe_models2 <- list(
  lr_gamedummy  = lm(gamedummy  ~ cm_dummy + I(1-cm_dummy) + 0, weights=ipw_ate, data=df),
  lr_gamecount  = lm(gamecount  ~ cm_dummy + I(1-cm_dummy) + 0, weights=ipw_ate, data=df),
  lr_gamesecond = lm(gamesecond ~ cm_dummy + I(1-cm_dummy) + 0, weights=ipw_ate, data=df),
  rf_gamedummy   = lm(gamedummy  ~ cm_dummy + I(1-cm_dummy) + 0, weights=ipw_rf,  data=df),
  rf_gamecount   = lm(gamecount  ~ cm_dummy + I(1-cm_dummy) + 0, weights=ipw_rf,  data=df),
  rf_gamesecond  = lm(gamesecond ~ cm_dummy + I(1-cm_dummy) + 0, weights=ipw_rf,  data=df),
  gbm_gamedummy  = lm(gamedummy  ~ cm_dummy + I(1-cm_dummy) + 0, weights=ipw_gbm, data=df),
  gbm_gamecount  = lm(gamecount  ~ cm_dummy + I(1-cm_dummy) + 0, weights=ipw_gbm, data=df),
  gbm_gamesecond = lm(gamesecond ~ cm_dummy + I(1-cm_dummy) + 0, weights=ipw_gbm, data=df)
)
lapply(ipwe_models2, tidy) %>%
  bind_rows(.id="outcome") %>%
  dplyr::select(outcome, term, estimate) %>% 
  spread(term, estimate) %>%
  mutate(causal_effect = cm_dummy - `I(1 - cm_dummy)`) %>% 
  knitr::kable(digits=3)
outcome cm_dummy I(1 - cm_dummy) causal_effect
gbm_gamecount 9.898 9.409 0.489
gbm_gamedummy 0.082 0.071 0.011
gbm_gamesecond 2756.185 2968.089 -211.904
lr_gamecount 13.636 8.317 5.320
lr_gamedummy 0.094 0.062 0.032
lr_gamesecond 4143.330 2639.412 1503.918
rf_gamecount 12.875 10.092 2.783
rf_gamedummy 0.088 0.072 0.016
rf_gamesecond 3462.625 3164.428 298.197

モデルによって結構違う結果になった． （マッチングで出したSEを見返すと，まあこれくらいはぶれるのかなという気もするが……．） さて，どのモデルで推定した傾向スコアが良いのだろうか．AUCが高いモデルが良いのだろうか．

ここでIPW推定量の定義を思い出すと，推定された傾向スコアで重み付けをするのだから， 推定値が真の傾向スコアの値に近くないと，過少あるいは過大な重みを付けることになり， 推定精度が悪くなりそうである．

ところがAUCというのは推定された傾向スコアの大小関係だけで決まる精度指標であって， 傾向スコアの値は見ていない． なのでAUCだけで決めてはダメなのではないかという疑問が出てくる．

そこで，推定された傾向スコアの値が真の傾向スコアに近いか確認するため， calibration plot を描いてみる．

# 傾向スコア推定値のプロット
caret::calibration(as.factor(cm_dummy) ~ ps+ps_rf+ps_gbm, class="1", data=df) %>%
  lattice::xyplot(auto.key = list(columns = 2))

ロジスティック回帰が最もよくフィットしているように見える． ランダムフォレストは傾向スコアが小さい領域では過大な推定で， 傾向スコアが大きい領域では過少な推定になっている． GBMは傾向スコアが80%以上に推定される人がいない （が，それ以外の領域ではランダムフォレストより良いように見える）．

ということは結局，ロジスティック回帰を使うのが良いのだろうか．

ただ，ランダムフォレストも傾向スコア推定値のフィットが悪いわりには， 因果効果の推定値がロジスティック回帰に結構近いような気もしてきて，分からない……．

決定木による因果効果の異質性の推定

最近，決定木で因果効果の異質性を推定する方法が提案された[Athey2015]． ありがたいことにgithubでRパッケージcausalTreeが公開されているので，これも試してみる．

アルゴリズムの詳細は論文を見てほしいが（自分はまだ細部まできちんと理解できていない）， 基本的には，普通の決定木が二乗誤差が小さくなるように枝を分割していくのに対して， causalTreeでは因果効果の誤差が小さくなるように枝を分割していくというイメージだ．

使い方は，線形モデルで調整効果を推定した例と同じで， + あらかじめ傾向スコアを推定しておき， + 傾向スコアの逆数で重み付けしてcausalTreeを実行 とすればよい（と思う．間違っていたら，すみません……）．

#devtools::install_github("susanathey/causalTree")
library(causalTree)
library(rpart.plot)

ct_gamesecond <- causalTree(
  gamesecond ~ child_dummy + area_kanto +area_tokai + area_keihanshin + area_keihan +
    T + F1 + F2 + F3 + M1 + M2 + M3,
  data = df[,1:35], treatment = df$cm_dummy, weights = df$ipw_ate,
  split.Rule = "CT", cv.option = "CT", 
  split.Honest = TRUE, cv.Honest = TRUE, split.Bucket = FALSE, 
  xval = 5, cp = 0, minsize = 100
)

cp_opt <- ct_gamesecond$cptable[,1][which.min(ct_gamesecond$cptable[,4])]
ct_gamesecond_opt <- prune(ct_gamesecond, cp_opt)
rpart.plot(ct_gamesecond_opt, extra=101, nn=TRUE)

これは割とおもしろい結果ではないかと思う．

というのも，岩波DS3本編の調整効果の分析では，子どもありの集団に対しては TVCMはネガティブな因果効果を及ぼしており，子供あり集団を除いて再分析しても 因果効果はあまり変わらないという結論だったのだが， 木を見ると，子どもありの集団に対する因果効果も一様ではないことが 一目で分かるからだ（と言っても，別に岩波DS3の分析と矛盾するわけではないけど）．

たとえば，

  • ノード27 :「M2層（男性35～49歳）・京浜居住・子どもあり」ではATE = 6268
  • ノード330:「F1層（女性20～34歳）・京阪居住・子どもあり」ではATE = -2625

といったノードを見ると，同じ居住地で子どもありでも性別・年齢で効果が正反対だったりする．

参考文献