はじめに
ビジネスの現場では、施策の実施に伴ってその効果検証が求められる。A/Bテストのように施策の実施対象が介入群と対照群に分かれていれば、因果推論の枠組みで効果検証を行うことができる。
介入群と対照群がそれぞれ複数のユニットを含み、ユニットごとの時系列データが利用できるとき、代表的な分析手法として差の差法（Difference in Differences; DiD）が挙げられる。
DiDは手法のデザインが直感的で分かりやすい一方、いくつかの課題がある。
- 平行トレンド仮定を満たさない場合の対処
- 比較対象を選ぶルールが曖昧なため、cherry-pickingを誘発するおそれ
- ユニット・時間に特有の係数を入れられないこと
- 信頼区間の説明が実務的に難しいこと
これらの課題に対処するべく、今回はベイズ因果推論の枠組みを用いたマルチファクターモデルによる推定方法と、そのRパッケージであるbpCausalを紹介する。原論文はPang, Liu and Xu(2021)である。
ベイジアンファクターモデル
準備
ユニットを$i = 1,2, \cdots, N$、時間を$t = 1,2, \cdots, T$で表す。二値の処置変数を$w_{it}$で表し、一度介入群にアサインされたならば、そのあと対照群に移ることはないものとする（staggered adoption）。
各ユニットについて、介入のタイミングを確率変数と考え、$a_i$で表す（$a_i \in A=\lbrace 1,2,\cdots, T, c\rbrace$）。$a_i = c (>T)$のとき、ユニット$i$は観測されるデータの中で介入されない。
介入効果を以下で定義する。
$\delta_{it} = y_{it}(a_i)-y_{it}(c) \quad (a_i\leq t \leq T)$
つまり、介入効果は介入群のユニット$i$に対して、介入後のアウトカムと反実仮想のアウトカムの差とされる。
割り当てメカニズム
鍵となる仮定は潜在的な無視可能性（Latent ignorability）である。
処置前の時変・時変でない共変量を含むベクトルを$X_i$、ユニットレベルの異質性（ユニット固定効果など）とユニット特有の時間トレンドを$U_i$で表すと、以下が成立する。
$Pr(a_i|X_i, Y_i(0),U_i) = Pr(a_i|X_i, Y_i(0)^{mis}, Y_i(0)^{obs}, U_i)=Pr(a_i|X_i, U_i)$
つまり、$X_i$と$U_i$で条件付ければ、$Y_i(0)$なる時系列は割り当てメカニズムと独立である。
この仮定は固定効果モデルでしばしば仮定される強い外生性の仮定を拡張したものである（強い外生性の仮定とは、$U_i$で条件付ければ過去のアウトカムが現在・未来の処置に影響しないというもの）。
階層・動的ファクターモデル
ユニット$i$の時間$t$におけるポテンシャルアウトカムは以下で定義される。
$y_{it}(c)=X_{it}'\beta_{it} + \gamma_i'f_t + \epsilon_{it}$
$\beta_{it} = \beta + \alpha_i + \xi_t$
$\xi_t = \phi_{\xi}\xi_{t-1} + e_t$
$f_t = \phi_f f_{t-1} + \nu_t$
ただし、$X_{it}$は観測される共変量（時間不変、ユニット不変を許す）、$\gamma_{i}'f_t$は潜在的なマルチファクター項。
スパースモデリング
$\beta$の事前分布は以下のような階層構造にすることで、ベイズ縮小を可能にしている。
$\beta_k|\tau_k^2 \sim N(0, \tau_k^2) \quad \forall 1\leq k \leq p_1$
$\tau_k^2|\lambda_{\beta} \sim Exp(\frac{\lambda_{\beta}^2}{2})$
$\lambda_{\beta}^2 \sim \mathcal G(a_1, a_2)$
ただし$p_1$は共変量の個数。$\lambda_{\beta}$はLassoにおける正則化パラメータに相当。
他のパラメータ$\alpha_i$,$\xi_i$, $\gamma_i$についても縮小アプローチがとられている。具体的には、それぞれ$\alpha_i = w_{\alpha}\cdot \tilde\alpha_i$, $\xi_i = w_{\xi}\cdot \tilde\xi_i$, $\gamma_i = w_{\gamma}\cdot \tilde\gamma_i$とre-parameterizeし、それぞれの重み$w$が0に近似されるのであればモデルに含まれないようにする。
実装
ここでは、原論文で扱われているドイツ再統一のレプリケーションを試みる。1989年のドイツ再統一が西ドイツのGDPに与えた影響の大きさを推定する。データはHARVARD Dataverseより取得した。
前処理
国がドイツならば介入群、そうでなければ対照群とする変数$treat$と、1990年以降ならば介入後、そうでなければ介入前とする変数$D$を用意する。また、連続値の欠損保管を行った。
df <- foreign::read.dta("repgermany.dta") |>
dplyr::mutate(treat = dplyr::case_when(
index == 7 ~ 1,
T ~ 0
)) |>
dplyr::mutate(D = dplyr::case_when(
index == 7 & year > 1990 ~ 1,
T ~ 0
)) |>
dplyr::mutate(dplyr::across(c(invest60, invest70), ~ .x * 100))
df_mice <- mice::mice(df, m = 5, maxit = 50, method = "pmm", seed = 500)|>
mice::complete(1)
MCMC
マルコフ連鎖モンテカルロ（MCMC）を使い、ポテンシャルアウトカムのモデル$f(y_{it}(0)|X_{it}, \theta_{it})$を得る（$\theta_{it}$はパラメータ）。
set.seed(123456789)
out <- bpCausal::bpCausal(data = df_mice,
index = c("country", "year"),
Yname = "gdp",
Dname = "D",
Xname = c(),
Zname = c(),
Aname = c("infrate", "trade", "schooling", "invest60", "invest70", "invest80", "industry"),
re = "both",
ar1 = TRUE,
r = 8,
niter = 15000,
burn = 5000,
xlasso = 1,
zlasso = 1,
alasso = 1,
flasso = 1,
a1 = 0.001, a2 = 0.001,
b1 = 0.001, b2 = 0.001,
c1 = 0.001, c2 = 0.001,
p1 = 0.001, p2 = 0.001)
ATTの可視化
介入群に対する平均処置効果（Average Treatment Effect on Treated; ATT）を取り出し、可視化する。
eout <- bpCausal::effSummary(out,
usr.id = NULL,
cumu = FALSE,
rela.period = TRUE)
df_att <- eout |>
as.data.frame() |>
tibble::as.tibble() |>
dplyr::mutate(year = seq(from = as.Date("1960-01-01"),
to = as.Date("2003-01-01"),
by = "1 year"))
ggplot2::ggplot(df_att) +
ggplot2::geom_line(ggplot2::aes(x = year, y = est.eff.estimated_ATT)) +
ggplot2::geom_ribbon(ggplot2::aes(x = year, ymin = est.eff.estimated_ATT_ci_l, ymax = est.eff.estimated_ATT_ci_u), alpha = 0.2) +
ggplot2::geom_hline(yintercept = 0, linetype = 2) +
ggplot2::geom_vline(xintercept = as.Date("1989-01-01"), linetype = 2) +
ggplot2::labs(y = "ATT") +
ggplot2::theme_bw()
（ATTのサイズが若干大きくなってしまった…）
まとめ
ビジネスの現場における施策の効果検証では、ユニットごとの効果の異質性を示せることと、不確実性の解釈が容易であること（施策の効果サイズがXX%の確率でこの区間に収まっている）が重要だと考える。
今回紹介したベイズマルチファクターモデルおよびbpCausalは、（デザインこそDiDや合成統制法のように直感的ではないものの）DiDに代わって上記のようなニーズを満たせるのではないだろうか。
参考文献
Pang, X., Liu, L., & Xu, Y. (2021). A Bayesian alternative to synthetic control for comparative case studies. Political Analysis, 30(2), 269–288.