慶應義塾大学・株式会社Nospareの菅澤です．
今回は，回帰不連続デザインを用いてサブグループごとの因果効果の推定を考えたときに，階層ベイズ的アプローチを導入することで安定的な推定を実現 することができる最近の研究成果について紹介します．本記事の最後にも情報を記載しておりますが，本記事は最近公開したワーキングペーパーに基づいています．
回帰不連続デザイン（RDD）とは
回帰不連続デザイン（RDD: Regression Discontinuity Design）の導入と理論的な詳細については石原先生の記事（記事1，記事2）をご覧ください．本記事では，最低限必要な問題設定について紹介しておきます．
$Y_i$を目的変数，$W_i$を処置変数とします．RDDでは，ある変数$X_i$が存在して$X_i$がある閾値$c$を超えるか否かで処置を受けるか否かが決まることを想定します．すなわち，$W_i=I(X_i\geq c)$であることを仮定します．このような$X_i$はランニング変数と呼ばれます．（例えば，一定期間にある金額以上を購入したユーザーにキャンペーンを実施したり，ある成績以上の学生に奨学金を与えるような状況が考えられます．）
RDDにおける推定対象は次のような量です．
\tau\equiv \lim_{x \downarrow c} E[Y_i|X_i=x]- \lim_{x \uparrow c} E[Y_i|X_i=x]
これは，閾値上の個体に対する平均処置効果 になっており，閾値から離れた人の処置効果については何も分からないことに注意が必要です．
$\tau$の表現から，閾値の右側（左側）のデータだけを使って$\lim_{x \downarrow c} E[Y_i|X_i=x]$（$\lim_{x \uparrow c} E[Y_i|X_i=x]$）をノンパラメトリック推定することで$\tau$を推定することができます．以下では，局所多項式回帰を使ったノンパラメトリック推定のみに焦点を当てます．
サブグループごとの因果効果の推定
通常のRDDではランニング変数以外について平均化した因果効果を考えています．一方で，実際のデータには様々な個体の情報が付随しているケースが多々あります．例えば，目的変数とランニング変数以外に，性（男女）・年代（20歳未満, 20代, ..., 80代以上）・居住地域（8地方区分）の情報が得られている状況を考えてみます．この場合の自然な興味として「性・年代・居住地域による因果効果の違い」が出てくると考えられます．このようなケースでは，3つの変数の組み合わせで$2\times 8\times 8=128$個のサブグループが構成されるため，3つの変数による因果効果の違いを考えることは128個のサブグループごとに因果効果がどのように異なるかを探索することに相当します．
目的変数とランニング変数のみを用いて行う通常のRDDでは，このようなサブグループについて平均化をしたような因果効果が推定されますが，サブグループごとの因果効果（i.e. 属性別の因果効果）に興味があるケースは応用上とても多くあります．
以下では，RDDを用いてサブグループごとの因果効果を推定することを考えていきます．
サブグループごとのRDDと問題点
数式的な表現をわかりやすくするため，サブグループのインデックスを明示的に導入した形で設定を書き直します．
$G$をサブグループの数として，$g=1,\ldots,G$をサブグループのインデックスとします．さらに，$g$番目のサブグループにおけるサンプルサイズを$n_g$とし，サブグループ内の個体のインデックスを$i=1,\ldots,n_g$で表します．このとき，$\sum_{g=1}^G n_g$がトータルのサンプルサイズです．
また，$Y_{ig}$を目的変数，$X_{ig}$をランニング変数とし，$c$を閾値とします．（閾値がグループごとに異なるケースにも対応できますが，簡単のために共通なものとしておきます．）このとき，推定対象は
\tau_g\equiv \lim_{x \downarrow c} E[Y_{ig}|X_{ig}=x]- \lim_{x \uparrow c} E[Y_{ig}|X_{ig}=x], \ \ \ g=1,\ldots,G
となります．$\tau_g$を推定する最も単純な方法としては，各サブグループごとに別々にRDDを実行する方法が考えられますが，推定が非常に不安定になってしまう可能性があります．その理由としては
- RDDは閾値周辺のデータのみしか使わないので，因果効果推定に寄与する実質的なサンプルサイズは小さくなる傾向にある
- サブグループの数$G$が多くなると，トータルサンプルサイズが大きくてもサブグループごとのサンプルサイズ$n_g$は小さくなる
の2点が考えれます．
そこで，サブグループごとに別々に推定を行うのではなく，階層的な構造を導入してサブグループごとの因果効果を同時推定 することで，互いの情報を借りあって（いわゆる階層ベイズ法における"borrowing strength"の効果によって）推定精度を向上させることを考えます．
階層RDD
サブグループごとのRDDに対して階層的な構造を導入するためには，目的変数$Y_{ig}$に対する何らかの確率モデルが必要です．
キーとなる事実は，$\tau_g$が以下の損失関数最小化によって推定できる点です．
L(\tau_g, \beta_g; h_g)=\frac12\sum_{i=1}^{n_g}K\left(\frac{|X_{ig}-c|}{h_g}\right)\left(Y_{ig}-\tau_gW_{ig}-Z_{ig}^\top\beta_g\right)^2, \ \ \ \ g=1,\ldots,G.
ここで，$K(\cdot)$は重み付け関数，$\beta_g$は局所多項式の係数，$h_g$はバンド幅であり，$Z_{ig}$は
Z_{ig}=(1, (X_{ig}-c)_{-}, (X_{ig}-c)_{+},\dots,(X_{ig}-c)_{-}^{q}, (X_{ig}-c)_{+}^{q})
のように定義されるベクトルです．既に説明したように，RDDでは閾値の左右で局所多項式回帰を当てはめ，閾値上（$X_{ig}=c$上）で差分を取ることで因果効果を推定しますが，そのような手順によって計算できる$\tau_g$の推定量と，上記の損失関数を$\beta_g$とともに最小化して得られる$\tau_g$の推定量は一致します．（バンド幅$h_g$の設定方法については後ほど言及します．）
以下では，上記の損失関数を利用して階層モデルを組んでいきます．
一般に，階層モデルを与えるには$Y_{ig}$に対する確率モデルを定める必要がありますが，今の問題設定では$Y_{ig}$に依存する損失関数しか定められておりません．そこで，一般化ベイズという考えを用いて階層モデル を作っていきます．一般化ベイズ法についての解説はこちらの記事をご覧ください．
一般化ベイズ法の枠組みを用いて，$(\tau_g,\beta_g)$の一般化事後分布を以下のように定めます．
\pi(\tau_g,\beta_g|Y_{1g},\ldots,Y_{n_gg};h_g)\propto \pi(\tau_g,\beta_g)\exp(-\omega L(\tau_g, \beta_g; h_g))
ここで，$\pi(\tau_g,\beta_g)$は$(\tau_g,\beta_g)$に対する事前分布であり，$\omega$はlearning rateと呼ばれるパラメータです．この事後分布は$Y_{ig}$に対して以下のような擬似モデルを仮定した場合の事後分布に一致します．
p(Y_{ig} | \tau_g, \beta_g, \omega ) \propto \phi(Y_{ig}; \tau_gW_{ig}+Z_{ig}^\top\beta_g, \omega^{-1})^{k_{ig}}
ただし，$k_{ig}=K(h_g^{-1}|X_{ig}-c|)$です．
あとは，サブグループごとのパラメータ$\tau_g$および$\beta_g=(\beta_{g1},\ldots,\beta_{gp})$に対する事前分布を定める必要がありますが，異なるサブグループ間の情報を共有するため以下のような事前分布（モデル）を考えます．
\tau_g\sim N(m_\tau, \psi_\tau), \ \ \ \ \beta_{gk}\sim N(m_{\beta_k}, \psi_{\beta_k}), \ \ \ k=1,\ldots,p.
ここで，$m_\tau, \psi_\tau, m_{\beta_k}, \psi_{\beta_k}$は未知のパラメータです．これらのパラメータを全データから推定することで，サブグループごとのinformativeな事前分布が定まり最終的な$\tau_g$の推定を安定化させる ことができます．
以上が提案モデルの概要です．グループごとのRDDを階層モデルの枠組みに埋め込んだ方法ということで，階層RDD（HRDD: Hierarchical RDD） と名付けます．
MCMCについて
$m_\tau, \psi_\tau, m_{\beta_k}, \psi_{\beta_k}$に対してさらなる事前分布を導入し，MCMCによって事後分布の計算を行います．これらのパラメータは正規分布の平均・分散パラメータなので，標準的な正規・逆ガンマ事前分布を用います．また，$\omega$（learning rate）に対して逆ガンマ事前分布を導入します．
このような設定もとで，事後分布はギブスサンプラーで簡単に推定することができます．（詳細は論文をご覧ください．）
バンド幅の選択
HRDDを実行するにはサブグループごとのバンド幅$h_g$を設定する必要があります．そのために，まずはサブグループごとのRDDを考え，既存のバンド幅選択手法（例えばRのrdrobustパッケージなど）を用いて最適な$h_g$を計算することができます．
いくつかの拡張
これまでに説明したHRDDの手法は単なる特殊ケースに過ぎません．損失関数を変えたり事前分布の構造を変えることで様々な目的変数の形や用途に利用することが可能 です．論文中で扱っている具体例の概要を以下に列挙します．
-
外れ値に対するロバスト化 : 二乗損失は外れ値の影響を大きく受けてしまいますが，$\omega$を$\omega u_{ig}$に変えて外れ値の影響を吸収する潜在変数$u_{ig}$を導入することでロバストなHRDDを実行することができます．
-
目的変数が二値変数だった場合 : 重み付きロジスティック損失を損失関数として用いることで二乗損失のときと同様にサブグループごとの因果効果が推定できます．
-
縮小事前分布の利用 : 因果効果$\tau_g$が0か否かに興味がある場合，変数選択的な要素としてspike-and-slab型の事前分布$\tau_g\sim \pi N(0, \varepsilon \psi_\tau) + (1-\pi)N(m_\tau, \psi_\tau)$を用いることができます．
また，論文中では扱っていないですが，チェックロスを用いることでサブグループごとの分位点処置効果を推定することもできるなど，他にも様々な形のHRDDを行うことが可能です．
数値例
最後に簡単な数値実験の結果を紹介します．（データ生成の具体的なシナリオについては論文をご参照ください．）
以下の図はシミュレーションデータに対する結果例です．100個のサブグループに対する真値とHRDDによる推定値（図中のHRDD-N）およびサブグループごとにRDDを適用した推定値（図中のsRDD-r）が図示されています．また，赤い縦線はHRDDによる95%信用区間，青い縦線がサブグループごとにRDDを適用して得られる95%信頼区間を表しています．
この結果から次のようなことが確認できます．
- サブグループごとにRDDを適用させるよりもHRDDによって得られる点推定値の方が精度が高い．
- サブグループごとにRDDを適用して得られる95%信頼区間よりもHRDDによって得られる95%信用区間の方が幅が狭い（より効率的な区間推定ができている）．
どちらも階層構造を導入することで適切に"borrowing strength"が出来ていることに起因する恩恵 であると考えられます．
論文中では，Colombian scholarship dataという奨学金のデータをもとに106個のサブグループに対する処置効果を推定した例も紹介しておりますので，興味があればご覧ください．
まとめ
今回は階層回帰不連続デザインを用いたサブグループごとの因果効果推定手法（HRDD）について紹介しました．HRDDはサブグループごとの因果効果の異質性を適切に推定する手法ですが，一般的な連続な説明変数に対する異質性に対応することはできません．（もちろん，連続変数をある閾値で切ることで強制的に離散化してサブグループを定義することでHRDDを用いることは可能です．）
階層ベイズの考え方を使うと，連続的な変数によって変化する異質因果効果をRDDの枠組みで推定することも可能でして，具体的な方法論や実装については現在研究を進めています．
冒頭でも述べたように，本記事は以下の論文の内容に基づいておりますので，詳細については論文をご参照ください．
また，提案手法を実装したコードはGithubで公開しておりますので，興味があれば使ってみてください．
株式会社Nospareでは統計学の様々な分野を専門とする研究者が所属しております．統計アドバイザリーやビジネスデータの分析につきましては株式会社Nospareまでお問い合わせください．
