回帰非連続デザインについて

## はじめに
東北大学/株式会社Nospareの石原です．今回は，計量経済学などで研究されている**回帰非連続デザイン** (Regression Discontinuity Design; RDD) について紹介します．回帰非連続デザインとは，制度の非連続性などを利用して因果効果を識別できる研究デザインのことをいいます．本記事では，回帰非連続デザインの識別・推定方法について紹介します．

## 回帰非連続デザインとは
処置を受けたときに実現するであろう潜在的な目的変数を $Y_{i}(1)$，処置を受けなかったときに実現するであろう潜在的な目的変数を $Y_{i}(0)$ と表記し，$D_i$ を処置を受けたかどうかを表すダミー変数とします．このとき，観測できる結果変数 $Y_i$ は
```math
Y_i = D_i Y_i(1) + (1-D_i) Y_i(0)
```
と書くことができます．さらに，1次元の変数 $X_i$ が存在して，$X_i$ が閾値 $c$ を超えるかどうかで処置を受けるかどうかが決まると仮定します．つまり，
```math
D_i = 1 \{ X_i \geq c \} \hspace{1in} (1)
```
が成り立っていると仮定します．このような研究デザインを**回帰非連続デザイン** (Regression Discontinuity Design; RDD) といい，$X_i$ を**ランニング変数**と呼びます．

回帰非連続デザインの例としては，次のような状況が考えられます．ある試験の点数 ($X_i$) が合格点 ($c$) を超えた人だけ奨学金を受けることができるとすると，奨学金を受けるかどうかを表すダミー変数 $D_i$ は (1) 式のように決定されます．例えば，$Y_i$ を大学に進学したかどうかを表すダミー変数とすると，奨学金が大学進学率に与える因果効果を分析することができます．回帰非連続デザインでは，このように処置が非連続に決定される状況を用いて，処置が結果変数に与える因果効果を推定します．

## 処置効果の識別
回帰非連続デザインでは，**条件付き期待値 $E[Y_i(d)|X_i=x]$ が閾値 $c$ で連続であることが，処置効果の識別条件となります**．実際に，この仮定の下では
```math
\begin{align}
& \lim_{x \downarrow c} E[Y_i|X_i=x] - \lim_{x \uparrow c} E[Y_i|X_i=x] \\
= & \lim_{x \downarrow c} E[Y_i(1)|X_i=x] - \lim_{x \uparrow c} E[Y_i(0)|X_i=x] \\
= & E[Y_i(1)-Y_i(0)|X_i=c]
\end{align}
```
となり，ランニング変数が閾値の人に対する平均処置効果 $E[Y_i(1)-Y_i(0)|X_i=c]$ が識別できます．

下の図は回帰非連続デザインの識別方法のイメージです．ランニング変数が閾値を下回った場合は $Y_i(1)$ を観測することはできないので，オレンジ色の点線の部分の条件付き期待値は観測できません．同様に，青色の点線の部分も観測できません．したがって，観測できる条件付き期待値 $E[Y_i|X_i=x]$ は閾値より下では青色の実線，閾値より上ではオレンジ色の実線となり，閾値でジャンプすることになります．回帰非連続デザインでは，このジャンプ幅が因果効果となります．

![fig1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2618276/c3da3745-1854-e8af-848e-59c164534d91.png)

回帰非連続デザインでは，通常の平均処置効果 $ATE=E[Y_i(1)-Y_i(0)]$ や処置を受けた人に対する平均処置効果 $ATT=E[Y_i(1)-Y_i(0)|D_i=1]$ は識別できないことに注意する必要があります．回帰非連続デザインで識別されるのは $E[Y_i(1)-Y_i(0)|X_i=c]$ であり，閾値から離れた人の処置効果については何も分かりません．奨学金の例では，試験の点数が合格点を大きく上回った（または，下回った）人に対する奨学金の大学進学率への因果効果は分からないということになります．

## 処置効果の推定方法
識別の議論で確認したように，$E[Y_i(1)-Y_i(0)|X_i=c]$ を推定するには次の2つの値を推定すればよいことになります：
```math
\lim_{x \downarrow c} E[Y_i|X_i=x], \hspace{0.3in} \lim_{x \uparrow c} E[Y_i|X_i=x]
```
したがって，閾値の右側（左側）のデータだけを使って $\lim_{x \downarrow c} E[Y_i|X_i=x]$（$\lim_{x \uparrow c} E[Y_i|X_i=x]$）をノンパラメトリック推定すればよいことが分かります．以降では，右側の推定問題（$\lim_{x \downarrow c} E[Y_i|X_i=x]$ の推定問題）について考えます．左側の推定問題についても，右側と同じ議論を使うことができます．

### Local Constant 推定量
条件付き期待値関数のノンパラメトリック推定で最もよく知られている方法は **Nadaraya-Watson 推定量**です．$\mu(x) \equiv E[Y_i|X_i=x]$ とすると，$\mu(x)$ の Nadaraya-Watson 推定量は次のように書けます：
```math
\hat{\mu}(x) = \frac{\sum_{i=1}^n K(\frac{X_i-x}{h}) Y_i}{\sum_{i=1}^n K(\frac{X_i-x}{h})}
```
ここで，$K(\cdot)$ はカーネル関数と呼ばれる重み付け関数で $\int K(u) du =1$ を満たします．また，$h$ は**バンド幅**と呼ばれるパラメータです．例えば，$K(u)=\frac{1}{2} \times 1 \\{ -1 \leq u \leq 1 \\}$ とすると，Nadaraya-Watson 推定量は $|X_i-x| \leq h$ となるサンプルの標本平均と一致します．つまり，バンド幅 $h$ は $x$ にどれだけ近いデータを使うかを決定するパラメータとなります．

Nadaraya-Watson 推定量は次の最小化問題を解くことでも得ることができます：
```math
\hat{\mu}(x) = \text{arg} \min_{a} \sum_{i=1}^n K\left( \frac{X_i-x}{h} \right) (Y_i-a)^2 
```
カーネル関数で重みを付けて定数項に回帰していると見なすこともできるので，Nadaraya-Watson 推定量は **Local Constant 推定量**と呼ばれることもあります．

### Local Linear 推定量
通常の条件付き期待値の推定では Local Constant 推定量を使えば十分ですが，回帰非連続デザインの場合には注意しなければいけない問題があります．それは，推定したい対象がデータの端点にあるという問題です．例えば，$\lim_{x \downarrow c} E[Y_i|X_i=x]$ の推定では，$X_i \geq c$ のデータを用いて $X_i=c$ の人の条件付き期待値を推定することになるので，端点の条件付き期待値を推定するという問題が生じます．この問題点を解消するために，回帰非連続デザインでは，**Local Linear 推定量**を用いることが提案されています．

まず，$K(u)=\frac{1}{2} \times 1 \\{ -1 \leq u \leq 1 \\}$ として，Local Constant 推定量にどのような問題が生じるかを議論します．下の図は Local Constant 推定量で端点の条件付き期待値を推定するときのイメージです．青色の実線が真の条件付き期待値関数 $\mu(x)$ であり，真の $\mu(x)$ は右下がりの直線であるとします．このとき，$X_i < c$のデータは使わないので，$x=c$ での Local Constant 推定量は $c \leq X_i \leq c+h$ のサンプルの標本平均となります．したがって，Local Constant 推定量の期待値は $\mu(x)$ を $c \leq x \leq c+h$ の範囲で平均化したものとなります．つまり，Local Constant 推定量の期待値は図のオレンジ色の点線のように，$x=c$ での条件付き期待値より小さくなります．このように，Local Constant 推定量で端点の推定を行うとバイアスが大きくなってしまうという問題があります．

![fig2.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2618276/a62ccb61-8994-bdda-3181-3db3a4590085.png)

この問題に対処するために，回帰非連続デザインでは Local Linear 推定量が使われています．$\mu(x)$ の Local Linear 推定量は，次の最小化問題を解くことで得られます：
```math
\left( \hat{a}(x), \hat{b}(x) \right) = \text{arg} \min_{a, b} \sum_{i=1}^n K\left( \frac{X_i-x}{h} \right) \left\{Y_i-a-b(X_i-x) \right\}^2
```
ここで，$\hat{a}(x)$ が $\mu(x)$ の Local Linear 推定量となります．Local Constant 推定量との違いは，最小化問題に線形項が追加されていることです．$K(u)=\frac{1}{2} \times 1 \\{ -1 \leq u \leq 1 \\}$ として $x=c$ での条件付き期待値の推定を考えると，Local Linear 推定量は次の最小化問題を解くことになります：
```math
\left( \hat{a}(c), \hat{b}(c) \right) = \text{arg} \min_{a, b} \sum_{i: \, c \leq X_i \leq c+h} \left\{Y_i-a-b(X_i-x) \right\}^2
```
したがって，Local Linear 推定量を求めることは，$c \leq X_i \leq c+h$ を満たすデータを用いて線形回帰を行うことに対応します．Local Constant 推定量の場合と異なり，$\mu(x)$ を線形近似することになるので $\mu(x)$ が直線の場合にはバイアスが出ないことが分かります．このように，Local Linear 推定量を使うことで，端点の推定に対してもバイアスを小さくすることができると考えられます．

### バンド幅の選択
Local Linear 推定量では，バンド幅をどのように選択するかが重要になります．バンド幅を大きくすると，使用するデータが増えるので推定量の分散は小さくなりますが，閾値から離れたデータを使うことになるので推定量のバイアスは大きくなります．一方で，バンド幅を小さくすると，閾値に近いデータだけを使うのでバイアスは小さくなりますが，使用するデータが減るので分散は大きくなります．このように，バンド幅の選択にはバイアスと分散のトレードオフがあります．そのため，適切なバンド幅を選択することが重要になります．詳しくは紹介しませんが，回帰非連続デザインでのバンド幅の選択については，Imbens and Kalyanaraman (2012), Calonico, Cattaneo, and Titiunik (2014), Arai and Ichimura (2018) 等の論文で提案されています．

## 最後に
今回は，回帰非連続デザインの識別・推定方法について紹介しました．回帰非連続デザインは経済学だけでなく政治学や教育学などの様々な分野の実証研究で使われています．"[A Practical Introduction to Regression Discontinuity Designs](https://www.amazon.co.jp/Introduction-Regression-Discontinuity-Quantitative-Computational/dp/1108710204)" というテキストなどもあるので，より詳しく知りたい方はテキストやサーベイ論文などを読んでみてください．また，R には "rdrobust" というパッケージがあるので，簡単に実装することができます．

株式会社Nospareでは統計学の様々な分野を専門とする研究者が所属しております．統計アドバイザリーやビジネスデータの分析につきましては[株式会社Nospare](https://bit.ly/38WtGZH)までお問い合わせください．
<a href="https://bit.ly/38WtGZH">
<img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1023506/d5c625c0-d6f4-680b-d702-3c40f11d6c39.jpeg" alt="株式会社Nospare" width=400>
</a>

## 参考文献
1. [Arai, Y., & Ichimura, H. (2018). Simultaneous selection of optimal bandwidths for the sharp regression discontinuity estimator. Quantitative Economics, 9(1), 441-482.](https://onlinelibrary.wiley.com/doi/abs/10.3982/QE590)
2. [Calonico, S., Cattaneo, M. D., & Titiunik, R. (2014). Robust nonparametric confidence intervals for regression‐discontinuity designs. *Econometrica*, 82(6), 2295-2326.](https://onlinelibrary.wiley.com/doi/abs/10.3982/ECTA11757)
3. [Imbens, G., & Kalyanaraman, K. (2012). Optimal bandwidth choice for the regression discontinuity estimator. *The Review of economic studies*, 79(3), 933-959.](https://academic.oup.com/restud/article-abstract/79/3/933/1533189?login=false)