効果検証入門 - 第2章 介入効果を測るための回帰分析(前半) -
Ei Taguchi
Abstract
効果検証入門第２章の2.1~2.2についてまとめたものです2.1 回帰分析の導入
単回帰分析
目的変数Yと変数Xの関係性を分析する
- Xが1単位増減したときにYがどう変化するかがわかる
下の図は実際に単回帰分析を実行したときのイメージである
ここでX=0のときのYの値が切片(\(\beta_0\))となる
この直線の傾きが\(\beta_1\)となる
そして各点と直線との誤差は\(u_i\)であらわされる
よって, \(Y_i\)は\[Y_i = \beta_0 + \beta_1X_i + u_i\]となる
ここで\(\beta_0\), \(\beta_1\)はパラメータである(分析で推定するもの)
手持ちのデータから得られた値は\(\hat{\beta_0}\)と\(\hat{\beta_1}\)とあらわす
この回帰直線の誤差が小さくなれば各点の説明精度も高いことになるので, 最も小さくなるように線をひく→最小2乗法(OLS)
手元にあるデータから得られたパラメータ\(\hat{\beta}\)は, 母集団上で得られる\(\beta\)に対する推定値となる
ここで, \(Y = E[Y|X] + u\)(\(E[Y|X]\)はYの条件付き期待値)と分解できるので, Yに対する近似値を推定するだけでなく, Yの条件付き期待値の近似値を推定していることにもなる
#相関のある乱数を発生
set.seed(123)
rho <- 0.8
x <- rnorm(n=100,mean=0,sd=3)
e1 <- rnorm(n=100,mean=0,sd=3)
e2 <- rnorm(n=100,mean=0,sd=3)
y1 <- sqrt(rho)*x+sqrt(1-rho)*e1
y2 <- sqrt(rho)*x+sqrt(1-rho)*e2
cor(y1, y2)## [1] 0.7745571ggplot()+
aes(x = y1, y = y2) +
geom_point() +
geom_smooth(method = "lm",
se = FALSE) +
labs(x = "X", y = "Y")## `geom_smooth()` using formula 'y ~ x'
効果分析のための単回帰分析
ここでは3種類の変数が登場する
Y : 被説明変数
Z : 介入変数
X : 共変量
本著の例では, Yはメールマーケティングによる購買量, Zはメール配信をしたかどうか, Xは昨年の購買量や最後の購買日など分析者が選定するものである(基本Xは複数ある→重回帰分析)
重回帰分析においても単回帰分析と同様にパラメータの推定値と条件付き期待値が求まる
モデル式にはXの2乗など柔軟に関数の形を導入できるが, これは分析者が適切に設定する必要がある
回帰分析による効果の推定
分析で知りたいのは介入があったかどうかで条件付き期待値に差があるのかどうかである. それぞれ\[E[Y|X, Z = 1] = \beta_0 + \beta_1*X + \beta_2 * X^2 + \beta_3*1\] \[E[Y|X, Z = 0] = \beta_0 + \beta_1*X + \beta_2 * X^2 + \beta_3*0\]
となるので, この差をみると, \[E[Y|X, Z = 1] -E[Y|X, Z = 0] = \beta_3\]である
- つまりここでは介入変数である\(\beta_3\)に興味があるということになる
回帰分析における有意差検定
- 手に入れたデータからの推定値\(\hat{\beta_3}\)について有意差検定を行う. 結果の解釈については第1章と同様
Rによるメールメーケティングの分析(回帰編)
- 実際にバイアスのあるメールメーケティングデータに対して, 回帰分析を実行する
- 目的変数 : 購入額(spend)
- 介入変数 : メールの有無(treatment)
- 共変量 : 過去の購入額(history)
#まずは使用するパッケージの読み込み
library(tidyverse)## ─ Attaching packages ──────────────────── tidyverse 1.3.0 ─## ✓ tibble 3.0.6 ✓ dplyr 1.0.4
## ✓ tidyr 1.1.2 ✓ stringr 1.4.0
## ✓ readr 1.4.0 ✓ forcats 0.5.1
## ✓ purrr 0.3.4## ─ Conflicts ───────────────────── tidyverse_conflicts() ─
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()library(broom)#回帰分析を実行(biased_dataは第1章で作成済)
biased_reg <- lm(spend ~ treatment + history, data = biased_data)
#分析結果の表示
summary(biased_reg)##
## Call:
## lm(formula = spend ~ treatment + history, data = biased_data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.74 -1.46 -1.26 -0.48 497.74
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.3241996 0.1444390 2.245 0.02480 *
## treatment 0.9026109 0.1743057 5.178 2.25e-07 ***
## history 0.0010927 0.0003366 3.246 0.00117 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 15.36 on 31860 degrees of freedom
## Multiple R-squared: 0.001339, Adjusted R-squared: 0.001276
## F-statistic: 21.35 on 2 and 31860 DF, p-value: 5.406e-10Coefficientsに注目すると, それぞれの\(\beta\)の値がEstimateから読み取れる
また, p値については\(Pr(>|t|)\)で確認できる
ここでtreatmentに注目すると推定値は0.9026109となっており, p値も十分に小さいことがわかるので「メール配信の効果はない」という帰無仮説は棄却できる
- よってメール配信をすると平均0.9程度売上があがる, と解釈できる
#coefficientsのみに注目したい場合broomパッケージで結果を表示する
biased_reg_coef <- tidy(biased_reg)
#データフレームで出力してくれる
biased_reg_coef## # A tibble: 3 x 5
## term estimate std.error statistic p.value
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 (Intercept) 0.324 0.144 2.24 0.0248
## 2 treatment 0.903 0.174 5.18 0.000000225
## 3 history 0.00109 0.000337 3.25 0.00117効果検証のための回帰分析で行わないこと
基本的にはtreatmentの結果以外には注目しない(その他のパラメータの値を改善するような努力をここではしないから)
- 有意差検定の結果もtreatment以外は解釈をしない
回帰分析の結果は予測器としても使えるが, 効果検証という観点ではspendの予測が目的ではない
- 予測式としては以下のようになる\[spend = 0.324 + 0.903treatment + 0.001history\]
2.2 回帰分析におけるバイアス
共変量の追加による効果への作用
共変量を増やすことで起きる結果を確認していく
まず, RCTを行ったデータとバイアスのあるデータでメール配信の効果を確認
rct_reg <- lm(spend ~ treatment, data = male_df) rct_reg_coef <- summary(rct_reg) %>% tidy() rct_reg_coef## # A tibble: 2 x 5 ## term estimate std.error statistic p.value ## <chr> <dbl> <dbl> <dbl> <dbl> ## 1 (Intercept) 0.653 0.103 6.36 2.09e-10 ## 2 treatment 0.770 0.145 5.30 1.16e- 7nonrct_reg <- lm(spend ~ treatment, data = biased_data) nonrct_reg_coef <- summary(nonrct_reg) %>% tidy() nonrct_reg_coef## # A tibble: 2 x 5 ## term estimate std.error statistic p.value ## <chr> <dbl> <dbl> <dbl> <dbl> ## 1 (Intercept) 0.548 0.127 4.32 0.0000156 ## 2 treatment 0.979 0.173 5.67 0.0000000143それぞれのtreatmentの推定値が, 第1章で算出したメール配信有無での消費額平均の差と同様の結果である
やはりセレクションバイアスによって過剰推定されていることがわかる
セレクションバイアスは人為的に発生しているおり, ここでは以下を共変量として追加する
最後の購入からの経過月数(recency)
昨年の購入額(history)
昨年どのチャネルから購入したか(channel)
nonrct_mreg <- lm(spend ~ treatment + recency + history + channel, data = biased_data) nonrct_mreg_coef <- summary(nonrct_mreg) %>% tidy() nonrct_mreg_coef## # A tibble: 6 x 5 ## term estimate std.error statistic p.value ## <chr> <dbl> <dbl> <dbl> <dbl> ## 1 (Intercept) 0.502 0.379 1.32 0.185 ## 2 treatment 0.847 0.178 4.74 0.00000211 ## 3 recency -0.0403 0.0259 -1.55 0.121 ## 4 history 0.00103 0.000375 2.74 0.00608 ## 5 channelPhone -0.00178 0.304 -0.00585 0.995 ## 6 channelWeb 0.226 0.303 0.745 0.456
- treatmentの値は0.847と改善され, RCTにおける推定値0.770に近づいた
脱落変数バイアス(OVB)
追加すべき共変量 = 目的変数Yと介入変数Zに対して相関のある変数
セレクションバイアスのあるモデルAとセレクションバイアスの影響を取り除いたモデルBを以下とする\[Y_i = \alpha_0 + \alpha_1Z_i + u_i\] \[Y_i = \beta_0 + \beta_1Z_i + \beta_2X_{omit, i} + e_i\]
モデルAの誤差項は\(u_i = \beta_2X_{omit, i} + e_i\)となっている
ここで\(X_{omit, i}\)はバイアスを取り除くのに必要な共変量である→脱落変数
モデルAにおける介入効果を示すパラメータ\(\alpha_1\)は, \[\alpha_1 = \beta_1 + \gamma_1\beta_2\]
この式からわかるように\(\alpha_1\)はモデルBにおける介入効果Zの推定値\(\beta_1\)と, 何かしらの値を取る\(\gamma_1\beta_2\)の和である
- \(\gamma_1\beta_2\)は脱落変数バイアス(OVB)と呼ばれている→必要な共変量を加えることでOVBの影響を取り除くことができる
\(\beta_2\)はモデルBにおける\(X_{omit, i}\)と\(Y_i\)との相関にあたる
\(\gamma_1\)は\(X_{omit, i}\)に対して \(Z_i\) を回帰させてときの回帰係数である(\(Z_i\)と\(X_{omit, i}\)の相関)
よって, \(\gamma_1\beta_2\)は\(X_{omit, i}\)と\(Y_i\)との相関に, \(X_{omit, i}\)と \(Z_i\) をかけたものである→推定値が歪められる
ここで統計的に有意かどうかを考慮すると, 本来追加すべき共変量も取り除かれる可能性もあり, OVBを発生させてしまう
- よって, 介入変数以外のパラメータに関しては有意差検定の結果を気にする必要はない
RによるOVBの確認
昨年の購入額(history)が含まれているかどうかでモデル比較をする. historyは売上に対して強い相関をもつためOVBが発生する.
モデルA\[spend_i = \alpha_0+\alpha_1treatment_i+\alpha_2recency_i+\alpha_3channel_i+e_i\]
モデルB\[spend_i = \beta_0+\beta_1treatment_i+\beta_2recency_i+\beta_3channel_i+\beta_4history_i+e_i\]
また, OVBは以下のようになる\[\alpha_1-\beta_1 =\gamma_1\beta_4\]
\(\gamma_1\)は回帰式であるモデルCであらわされる\[history_i = \gamma_0 +\gamma_1treatmen_i+\gamma_2recency_i+\gamma_3channel_i+\epsilon_i\]
上記モデルA~Cの回帰分析を実行していく
#それぞれのモデルをベクトルに格納 formula_vec <- c(spend ~ treatment + recency + channel, spend ~ treatment + recency + channel + history, history ~ treatment + channel + recency) #名前をつける #LETTERSで大文字アルファベットを指定できる #namesでデータに名前がつく names(formula_vec) <- paste("reg", LETTERS[1:3], sep = "_") #enframeでデータフレーム化できる models <- formula_vec %>% enframe(name = "model_index", value = "formula") #まとめて回帰分析を実行する #mapでそれを実行している #.xで使うベクトルやリストを指定, .fで関数や式を指定する df_models <- models %>% #回帰分析の部分 mutate(model = map(.x = formula, .f = lm, data = biased_data)) %>% #必要な部分だけ取り出す mutate(lm_result = map(.x = model, .f = tidy)) #モデルの結果を整える df_results <- df_models %>% #formulaを文字列に直す mutate(formula = as.character(formula)) %>% #selectで取り出す列を指定(このままだとネストデータ) select(formula, model_index, lm_result) %>% #分析しやすいようにunnestする(colsで取り出したい部分を選ぶ) #lm_resultに格納されているデータが展開される unnest(cols = c(lm_result)) df_results## # A tibble: 16 x 7 ## formula model_index term estimate std.error statistic p.value ## <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl> ## 1 spend ~ treatment… reg_A (Inter… 1.10e+0 0.311 3.52 4.29e- 4 ## 2 spend ~ treatment… reg_A treatm… 8.75e-1 0.178 4.91 9.24e- 7 ## 3 spend ~ treatment… reg_A recency -5.52e-2 0.0254 -2.17 2.97e- 2 ## 4 spend ~ treatment… reg_A channe… -3.12e-1 0.282 -1.11 2.68e- 1 ## 5 spend ~ treatment… reg_A channe… -8.47e-2 0.282 -0.301 7.64e- 1 ## 6 spend ~ treatment… reg_B (Inter… 5.02e-1 0.379 1.32 1.85e- 1 ## 7 spend ~ treatment… reg_B treatm… 8.47e-1 0.178 4.74 2.11e- 6 ## 8 spend ~ treatment… reg_B recency -4.03e-2 0.0259 -1.55 1.21e- 1 ## 9 spend ~ treatment… reg_B channe… -1.78e-3 0.304 -0.00585 9.95e- 1 ## 10 spend ~ treatment… reg_B channe… 2.26e-1 0.303 0.745 4.56e- 1 ## 11 spend ~ treatment… reg_B history 1.03e-3 0.000375 2.74 6.08e- 3 ## 12 history ~ treatme… reg_C (Inter… 5.77e+2 4.65 124. 0. ## 13 history ~ treatme… reg_C treatm… 2.72e+1 2.66 10.2 1.39e- 24 ## 14 history ~ treatme… reg_C channe… -3.02e+2 4.21 -71.6 0. ## 15 history ~ treatme… reg_C channe… -3.02e+2 4.20 -71.8 0. ## 16 history ~ treatme… reg_C recency -1.45e+1 0.379 -38.2 1.97e-312ここで確認したいのは\(\alpha_1-\beta_1\)と\(\gamma_1\beta_4\)である.
まず, \(\alpha_1\), \(\beta_1\), \(\gamma_1\)を取り出す
#treatment行の中のestimateを取り出す treatment_coef <- df_results %>% filter(term == "treatment") %>% pull(estimate) treatment_coef## [1] 0.8746297 0.8465757 27.2395995- 1つ目が\(\alpha_1\), 2つ目が\(\beta_1\), 3つ目が\(\gamma_1\)をあらわす
次に, \(\beta_4\)を取り出す
history_coef <- df_results %>% filter(model_index == "reg_B", term == "history") %>% pull(estimate) history_coef## [1] 0.001029897
OVBは\(\alpha_1-\beta_1\), もしくは\(\gamma_1\beta_4\)なので,
#alpha-beta coef_gap = treatment_coef[1] - treatment_coef[2] #gamma*beta OVB <- treatment_coef[3] * history_coef coef_gap## [1] 0.02805398OVB## [1] 0.02805398よって共変量を追加したモデルと追加しなかったモデルとの推定される効果の差がOVBと一致することがわかる
OVBが与えてくれる情報
介入変数Zと目的変数Yの両方と関係性をもつような変数をモデルに含めるとバイアスを低減させることができるとわかった. このような変数を交絡因子と呼ぶ
Yを予測するという観点ではYとの相関がないものを選ぶべきだが, OVBの値は0になるのでその変数をモデルに加えても効果の推定値は変わらない
バイアスを減らすにはYとZどちらともある程度相関のある変数を加えることがポイント
バイアスを発生させるような変数(\(X_{omit, i}\))が手に入らなくても, その変数とY, Zの関係性が+, -どちらになるかを考えると推定が過剰か過少かが想定できる
Conditional Independence Assumption
効果検証において, モデルに含まれていない変数によるOVBが0になることが理想である
- これは介入Zがランダムに割り振られていることに等しい→CIA(Conditional Independence Assumption)
変数の選び方とモデルの評価
共変量の選択とCIA
分析者は以下のようなステップでモデルを作る
どのように介入が割り当てられているのかを考える
介入の割り当てに関係のありそうな共変量を選択する
選択した共変量とYとの関係性を考え関数を決める
このモデルがCIAを満たしている必要があるが, 2つの問題が考えられる
バイアスの評価ができないという問題
得られた効果の推定値がどの程度バイアスを持っているのかを評価できない
- OVBはモデル間のバイアスの変化を示しており, 残りのバイアスの大きさを示しているわけではない
必要な共変量がデータにはないという問題
手持ちのデータの中ではバイアスを最小にすることは可能でも, 他にバイアスを発生させる要因があると完全になくすことはできない
ドメイン知識でバイアスを推察することもできるがこれは主観的評価が強まる
これらの解決策としては操作変数法, 固定効果モデル, DIDなどが挙げられる
Sensitivity Analysis
手持ちのデータには含まれない変数がセレクションバイアスを引き起こしているかの評価方法→Sensitivity Analysis
分析者が重要だと認識している共変量以外の共変量をモデルから抜くことで, 効果の推定値が大きく変動するかどうかを確認する
その結果推定値が大きく変化しないのであれば, ほかの変数による影響を受けにくいことを示唆する
Post treatment bias
OVBが0でない変数であればなんでも入れていいわけではない
ここではサイト来訪履歴(visit)を扱う
cor_visit_treatment <- lm(treatment ~ visit + channel + recency + history, data = biased_data) %>% tidy() cor_visit_treatment## # A tibble: 6 x 5 ## term estimate std.error statistic p.value ## <chr> <dbl> <dbl> <dbl> <dbl> ## 1 (Intercept) 0.726 0.0112 65.0 0. ## 2 visit 0.144 0.00761 18.9 2.30e- 79 ## 3 channelPhone -0.0751 0.00948 -7.92 2.51e- 15 ## 4 channelWeb -0.0738 0.00947 -7.80 6.38e- 15 ## 5 recency -0.0292 0.000795 -36.7 3.36e-289 ## 6 history 0.000109 0.0000117 9.31 1.41e- 20結果から, 0.144程度相関があると考えられる
次に実際にモデルにvisitを追加してみる
bad_control_reg <- lm(spend ~ treatment + channel + recency + history + visit, data = biased_data) %>% tidy() bad_control_reg## # A tibble: 7 x 5 ## term estimate std.error statistic p.value ## <chr> <dbl> <dbl> <dbl> <dbl> ## 1 (Intercept) -0.438 0.376 -1.16 2.44e- 1 ## 2 treatment 0.294 0.177 1.66 9.68e- 2 ## 3 channelPhone 0.121 0.300 0.403 6.87e- 1 ## 4 channelWeb 0.117 0.299 0.392 6.95e- 1 ## 5 recency 0.00988 0.0257 0.385 7.00e- 1 ## 6 history 0.000525 0.000371 1.42 1.57e- 1 ## 7 visit 7.16 0.242 29.6 3.85e-190treatmentは0.294となり, 大きく減少した
メール配信されたグループは購買意欲の高いユーザも購買意欲の弱いユーザもメールがあるからサイトに来訪している
一方メール配信されなかったグループは元々購買意欲の高いユーザだけサイトに来訪していることになる
よって, サイト来訪有無で売上を比較するとメール配信されなかったグループの方が平均売上が高くなってしまう
介入によって影響を受けた変数を入れることによって起きるバイアスをPost Treatment Biasという
visitのように介入の後で決まるような値は分析から除外するべき
なんらかの施策を評価する際は利用者や購入者だけでなく, 対象となるユーザ全体のデータを用意する必要がある