2014/12/18 実践編を更新 [計量経済学][時事ネタ?] 非連続回帰デザイン (RDD) 実践編 - ill-identified diary

概要

数式が多い
Sharp RD と Fuzzy RD の手法について前回より厳密かつ実用性のある説明をする.
R の rdd パッケージを使って実演する. → 理論の説明で長くなったので後半に分割

潜在効果モデル

ある施策によってでた結果の効果を表すには, その施策を受けた場合と受けなかった場合の結果を比較すればいい. ある人 (あるいは個体) に対して, $Y_i(1)$ , $Y_{i}(0)$ をそれぞれ施策を受けた/受けなかった人の結果の大きさを表すとすると, その差の期待値が平均処置効果 (ATE; Average Treatment Effect) と呼ばれる.

$\begin{equation} \mathit{ATE} =\mathrm{E}Y_{i}(1) -\mathrm{E}Y_{i}(0) \end{equation}$

また, もうひとつこの期待値を, 条件付き期待値とした

$\begin{equation} \mathit{ATT} = \mathrm{E}\left[ Y_{i}(1)-Y_{i}(0) | D_{i} = 1 \right] \end{equation}$

を ATT (Average Treatment effect on Treated ) と呼ぶ. $D_{i}$ は施策の有無を表すダミー変数で, 1 なら処置あり, ゼロならなしである. ATT は ATE と違い, 全ての人に処置を与えた場合, どれくらいの底上げ効果があるか, を意味すると解釈でき, ATE より ATT を「因果効果」として扱うことが多い. なお, 前回も言及した黒澤 (2005) では, ATE と ATT を別の物としてはっきり定義していたが, ATT も平均処置効果 ( average treatment effect) と呼ばれることが多い. これは, 次に説明するように, 一定の仮定のもとで, 両者を同一視できるからである.

さて, この ATT を測定する際に, 次のようなことが問題になる. 人ひとりに対しては, 施策を行った場合・行ってない場合のいずれか片方の状態しか観察できない, ということである. たとえば今日, 卵かけご飯を食べた一日, あるいは何も食べなかった一日を, あなたはどちらか一方しか選べないし, その結果どういう一日になったかは選んだほうしか経験できない. これが「潜在効果」である. 当たり前といえば当たり前だが, この事実が ATT の測定にどのような影響を与えるか. まず, $\mathrm{E}[Y_{i}(0)|D_{i} = 1 ]$ は観測できないため, これを避ける必要がある. ATT の式を変形すると,

$\begin{equation} \mathit{ATT} = \mathit{E}\left[ Y_{i}(1) - Y_{i}(0) | D_{i}=1\right] \\ = \mathrm{E}\left[Y_{i}(1)|D_{i}=1\right]-\mathrm{E}\left[Y_{i}(0)|D_{i}=0\right] \\ + \left\{ \mathrm{E}\left[Y_{i}(0)|D_{i}=0\right]-\mathrm{E}\left[Y_{i}(0)|D_{i}=1\right]\right\} \end{equation}$

(なぜかアンパサンドが文字化けするので数式の位置調整をしていません. なにが原因?)

となる. 第3項にある $\mathrm{E}[Y_{i}(0)|D_{i}=1]$ はやはり原理上観測できないが, もしこれが $\mathrm{E}[Y_{i}(0)|D_{i}=0$ ] と等しいなら, 第3項がゼロとなり, ATT を求めることができる.

第3項は, (現実に) 施策を受けなかった人の結果と, その人が (現実に反して) 施策を受けたとした場合の結果の差である. 例えば, この施策が本人の意志とは全く無関係にもたらされるのなら, 2つの期待値は等しいはずだ. つまり, 大数の法則の類推だ. 互いに独立なサンプルの平均値はある値に収束する. では逆に, 等しくならない場合とは何か. それがセレクションバイアス (自己選択バイアス) である.

セレクションバイアスの例として, 「任意参加の勉強会に参加した場合の学習効果」を考えてみよう. 勉強会は強制参加でもなく, 参加費もない. そしてその日あなたは特にスケジュールが埋まっていない, と仮定する. これなら, あなたは単純に, 怠けたいか, より自分を磨きたいかという好みの問題だけで勉強会の参加不参加を決めることができる. このような条件で, あなた以外にも勉強会に参加した数十人と, 勉強会に参加しなかった数十人を集め, 勉強会が終わった後でテストを課してそれぞれのグループの得点の平均点を比較してみよう. この平均点の差は, 勉強会による効果とイコールで繋げられるだろうか.

自由に過ごせる休みの日にあえて勉強会に参加した人間は, 間違いなく勉強熱心だ. 多分普段から勉強を欠かさない性格だろう. 一方, 参加しなかった人間は, 参加しなかった人間よりものぐさだ. 多分普段からあまり自主勉強をしないだろう. それならそもそも勉強会をする前から, 両者の潜在的な成績に差がある可能性が高い. これがセレクションバイアスの一例だ. このように, 参加者の意志が介在する実験では, 常にバイアスが発生する可能性がある. そこで, バイアスの解決方法の1つとして, 非連続回帰デザイン (RDD; Regression-Discontinuity Design) を用いる.

ちなみに, この例から「参加者・非参加者に対して, 勉強会の後だけでなく前にもテストをすれば勉強会の効果を測定できるのではないか」と思った人は鋭い. 実際, そのようなコンセプトに基づく手法があり, DID, Difference-In-Differences method と呼ばれる. これも訳語がまだ定着してないが, まれに直訳で「差の差」推定法などと呼ばれる. 今回は DID については書かないが, 暇があれば書く.

2019/10/14 追記: かなり長いが以下の記事で DID の考え方について言及した.

非連続回帰デザイン

これ以降は, 前回を先に読み大まかな方針を先に知っておいたほうがいいかもしれない.

ある施策がなされるかどうかは, 変数 $X_{i}$ の値が一定以上かどうかで決まる (一定以下でもいいが, ここでは単純化のため一定以上とする).

$\begin{equation} D_{i} = 1\{ X_{i} \geq c \}. \end{equation}$

補足しておくと, $1\{\cdot \}$ は指示関数で, カッコ内の条件を満たす時に 1, そうでないときにゼロとなる関数である. この場合, $X_{i}$ が $c$ 以上だと施策がなされる. この $c$ をカットオフポイントという. また, $X_{i}$ を running variable とか forcing variable とか呼ぶ. ここではスカラとして扱うが, 多変量でも良い.

このとき, 先の ATT の式の条件付き期待値には $D_{i}=1, 0$ と条件があったが, これを $X \lesseqgtr c$ に置き換えることができる. しかし, これだけではまだバイアスを消すことはできない. RDD は, 条件付き期待値の極限を取ることで, ATT を求める.

処置群 (施策がなされた) と対照群 (施策がなされてない) に分かれる瞬間であるカットオフポイントで, 両者の結果変数 $Y_{i}$ は大きく変わるはずである. 断層のように「非連続な」回帰曲線となることから, 「回帰非連続」デザインという. 参考として, $c=0$ の場合のプロット図とそれを作成した際の R のコードを用意した. ただしデータは乱数で適当に発生させただけである. 緑の直線は $y_{i} = \alpha + \beta x_{i} + \tau D_i +\varepsilon$ を線形回帰した直線である.

f:id:ill-identified:20141207155700p:plain

library(ggplot2)
set.seed(3)
N <- 300
df         <- data.frame(obs =1:N, x = rnorm(n = N, mean = 0), y=rnorm(n=N, mean=10))
df$y       <- ifelse(df$x >=0, df$y+rnorm(n=N, mean=2), df$y) + rnorm(n=N,mean=2)*(df$x)/10
df$treated <- ifelse(df$x >=0, T,F)
srd <- lm(y~x+treated, data=df)
# ggplot2 用に回帰直線のデータ
grid.x <- c(-2, 0, 2)
srd.df <- data.frame(x=rep(grid.x,2), variable=c(rep("y(0)",3),rep("y(1)",3)),
                     y=c(grid.x*srd$coefficients[2]+srd$coefficients[1], 
                         grid.x*srd$coefficients[2]+srd$coefficients[1]+srd$coefficients[3])
)
srd.df[srd.df$x==2  & srd.df$variable=="y(0)","y"] <- NA
srd.df[srd.df$x==-2 & srd.df$variable=="y(1)","y"] <- NA

graph.base <- ggplot(data=df) + geom_vline() + geom_point(aes(x,y, colour=treated))  + geom_line(data=srd.df, mapping = aes(x,y,group=variable), color="green", size=1.5 )

graph.base + theme_bw() + theme(title=element_text(size=20), axis.title= element_text(size=15), axis.title.y = element_text(angle = 0, hjust = .5)) + scale_color_discrete(name="", label=c("対照群","処置群"))

結果変数である $Y_{i}$ は, $X_{i}=c$ 上では, 左から近づけば $Y_{i}=Y_{i}(0)$ , 右から近づけば $Y_{i}=Y_{i}(1)$ となる. よって, 同じ $X_{i}=c$ でも, 左極限か右極限かで $\mathrm{E} Y_{i}|X_{i}$ の値は異なる.

$\begin{align} \lim_{x\uparrow c} \mathrm{E}\left[ Y_{i} | D_{i}=0 \right] = \lim_{x\uparrow c}\mathrm{E}\left[Y_{i}(0)|X_{i} = c\right], \\ \mathrm{E}\left[Y_{i}(1)|X_{i}=c\right] = \lim_{x \downarrow c}\mathrm{E}\left[Y_{i}(1)|X_{i}\right]. \end{align}$

よって, ATT は左右極限の差になる. この極限を, それぞれ $\mu_{+}$ , $\mu_{-}$ とおき,

$\begin{equation} \mathit{ATT} = \mu_{+} -\mu_{-} \end{equation}$

と表す. ここで, 簡単な線形回帰モデルを仮定し, 条件付き期待値関数を,

$\begin{equation} \mathrm{E}\left[ Y_{i} |X_{i}\right]= \alpha + \beta X_{i} +\tau D_{i} \end{equation}$

とする. $X_{i}=c$ において, $D_{i}$ は非連続だが, $X_{i}$ は連続である. よって, $\mu_{+}=\alpha + \beta X_{i} + \tau \cdot 1$ , $\mu_{-}=\alpha + \beta X_{i} + \tau \cdot 0$ となるので, 平均処置効果は $D_{i}$ の係数 $\tau$ と等しくなる.

もちろん回帰モデルはもう少し複雑な形状にすることができる. よく使わるのは $X_{i}$ の多項式である. また, 他の共変量を追加することも可能である.

順序が逆になったが, 数学的な厳密さのため書いておくが, この式が成り立つためには条件付き期待値 $\mathrm{E}Y_{i}(1)|X_{i}$ , $\mathrm{E}Y_{i}(0) |X_i$ がそれぞれ右・左連続であり, また $Y_{i}(1)|X_{i}$ , $Y_{i}(0)|X_{i}$ の条件付き分布もそれぞれ同様に連続である, という仮定をおいている.

Fuzzy RD

ここまでで説明した RDD に対して, Fuzzy RD は, 必ずしもカットオフポイント $X_{i}=c$ で変化が置きない場合である. よって, これまでの RDD (以降 Sharp RD と呼ぶ) と異なり条件付き期待値の, $D_{i}=0,1$ と $X_{i}$ の関係が成り立たない. 正確には,

$\begin{equation} \lim_{x\uparrow c}\mathrm{P}(D_{i}=1|X_{i}=x)\neq \lim_{x\downarrow c}\mathrm{P}(D_{i}=1|X_{i}=x), \end{equation}$

ということになる.

f:id:ill-identified:20141207153440p:plain

Sharp/Fuzzy RD の違いの図. 横軸が runnig variable $X_{i}$ , 縦軸が割り当て変数 $D_{i}=1$ である確率. Sharp ならカットオフポイント $c$ で即座に変化するが, Fuzzy では必ずしもそうならない.

Fuzzy RD では Sharp RD の式は成り立たず, 平均処置効果を求めるために新たな方法が必要となる. まず, これまでゼロか1の2値変数として扱ってきた $D_{i}$ を, 割り当て関数 $D_{i}(x)$ として扱う. この $D_{i}(x)$ は $x$ の単調関数と仮定する.

$Y_{i}=Y_{i}(0) + D_{i} \left( Y_{i}(1) - Y_{i}(0) \right)$ であるから,

$\begin{align*} \mathrm{E}Y_{i}|X_{i} = \mathrm{E} \left[Y_{i}(0)|X_{i}\right]+\mathrm{E}\left[D_{i}|X_{i}\right]\mathrm{E}\left[Y_{i}(1)-Y_{i}(0)|X_{i} \right] \\ = \mathrm{E}\left[Y_{i}(0)|X_{i}\right]+\mathrm{E}\left[D_{i}|X_{i}\right]\tau_{\mathit{FRD}} \end{align*}$

ここで, また $X_{i}=c$ における左右の極限を取る.

$\begin{align*} \mu_{+} = \lim_{x\downarrow c}\mathrm{E}\left[Y_{i}|X_{i}=x\right] \\ = \lim_{x\downarrow c}\mathrm{E} \left[Y_{i}(0)|X_{i}=c\right]+\lim_{x\downarrow c}\mathrm{E} \left[D_{i}|X_{i}=c \right] \tau_{\mathit{FRD}}, \end{align*}$

よって,

$\begin{align*} \mu_{+} -\mu_{-} =\lim_{x\downarrow c} \mathrm{E}\left[D_{i}|X_{i}=c \right] \tau_{\mathit{FRD}} \\ -\lim_{x\uparrow c}\mathrm{E} \left[D_{i}|X_{i}=c \right] \tau_{\mathit{FRD}}. \end{align*}$

ここで, $\lim_{x\downarrow c} \mathrm{E} [Y_{i}|X_{i}=x ] \neq \lim_{x \uparrow c}\mathrm{E}[ Y_{i} |X_{i}=x$ ] であるから, 両辺をこの2つの期待値の極限の差で割ると,

$\begin{equation} \tau_{\mathit{FRD}} = \frac{\lim_{x\downarrow c}\mathrm{E}\left[Y_{i}(1)|X_{i}=x\right]-\lim_{x\uparrow c}\mathrm{E}\left[Y_{i}(0)|X_{i}=x\right]}{\lim_{x\downarrow c}\mathrm{E}\left[D_{i}|X_{i}=x\right]-\lim_{x\uparrow c}\mathrm{E}\left[D_{i}|X_{i}=x\right] } \\ = \frac{\mu_{+}-\mu_{-}}{D_{+}-D_{-}} \end{equation}$

ここで, $D_{-}$ , $D_{+}$ は $D_{i}$ の条件付き期待値の左右の極限を表す. Sharp RD の場合に当てはめると, $\lim_{x\downarrow c}\mathrm{E}\left[D_{i}|X_{i}=x\right]=1$ , $\lim_{x\uparrow c}\mathrm{E}\left[D_{i}|X_{i}=x\right]=0$ だから, 実はこの式は Sharp, Fuzzy いずれでも当てはまる一般形だと分かる.

これはいわゆるワルド推定量と呼ばれるもので, Hahn et al. (2001) "Identification and estimation of treatment effects with a regression‐discontinuity design" において, 2段階最小二乗法による推定量と等価であることが示されている. 単純な式の例を挙げると, 新たに外生的な割り当て変数 $Z_{i}$ を用意し, 1段階目の回帰式

$\begin{equation} D_{i} = \rho Z_{i} + \gamma X_{i} + \eta_{i} \end{equation}$

を推定し, そこから得た $D_{i}$ の予測値 $\hat{D_{i} }$ を用いて, 2段階目の

$\begin{equation} Y_{i} = \tau \hat{D_{i} } + \beta X_{i} + \varepsilon_{i} \end{equation}$

を最小二乗法を用いて $\tau$ を推定する. 係数の推定だけならば2回最小二乗法を行えばいいが, 2段階目の説明変数が予測値であるため, 検定を行う際は標準誤差に注意が必要である (2段階最小二乗法の詳しいやり方は他の適当な文献で). あるいは, 2段階最小二乗法というよりも傾向スコア (propensity score) 法の類推で考えるといいかもしれない.

推定方法

推定方法はこれまでに2種類のほうほうが提案されている. 1つはノンパラメトリックな手法から生まれた局所線形回帰 (local linear regression; LLR), もうひとつは多項式近似である. いずれの方法も Sharp RD, Fuzzy RD 両方に適用できる. ただし, それぞれ一長一短の特徴がある.

局所線形回帰 (LLR)

カーネル法を用いて, ノンパラメトリックな推定をすることが可能である. しかし, カーネル法では推定結果にバイアスをもたらす「境界問題」が存在する. そのため, Hahn et al. (2001) では Fan (1991) "Design-adaptive Nonparametric Regression" を引き合いに出し, カーネル法のなかでも特にこの問題を回避できる Local Linear Regression (LLR) を使用することを提案している. LLR は, 矩形カーネルを用いて $c$ 周辺を推定する. 平たく言えば, カットオフポイント $c$ の両側の $h/2$ の範囲だけサンプルとして, $\begin{equation} Y_{i} = \alpha + \beta (X_{i} -c) + \tau D_{i} +\varepsilon_{i} \end{equation}$ を推定する. Fan (1991) では, より普及しているカーネル密度推定量の Nadaraya-Watson 推定量より, LLR のほうが境界値のバイアスが少ないことを指摘している. 結局, シンプルな線形回帰モデルに立ち返った方がうまくいくということである.

カーネル密度推定の話は日本語文献も充実しているので詳しい話はしない. PRML の訳本とか, RjpWiki の『ヒストグラムと密度の推定』とか, あるいは日本経済学会のノンパラメトリック, セミパラメトリック計量経済分析のセッション資料スライドも簡潔で概要をつかむには良い. 英語でもいいなら, 計量経済学のテキストで, Cameron and Trivedi (2005) "Microeconometrics: Methods and Applications" (9章など) や Li and Racine (2006) "Nonparametric Econmetrics: Theory and Practice" などが内容が豊富だ (Wooldridge (2010) にはあまりノンパラメトリックの話は書かれていない.).

Microeconometrics: Methods and Applications

作者:Cameron, A. Colin,Trivedi, Pravin K.
発売日: 2005/05/09
メディア: ハードカバー

Nonparametric Econometrics: Theory and Practice

作者:Li, Q.,Racine, Jeffrey Scott
発売日: 2006/11/27
メディア: ハードカバー

バンド幅の選び方

ただし, モデルの見かけ上単純な線形モデルだが, バンド幅 $h$ の選び方よって誤差が変わる. カーネル法では平均時乗誤差 (MSE) の最小化からバンド幅を決定する. しかし, この場合に限って言えば, この方法には問題があることを Imbens and Kalyanaraman (2012) "Optimal Bandwidth Choice for the Regression Discontinuity Estimator" で指摘されている. 後ほど紹介する rdd パッケージでは, Imbens and Kalyanaraman (2012) で提案された方法でバンド幅を決定する関数が用意されている. これは Fuzzy でも使用可能である.

LLRの弱点

LLR の問題として, $\tau$ の有意性の検定が難しいことがある. McCrary (2008) "Manipulation of the running variable in the regression discontinuity design: A density test" において, 検定方法が提案されているが, この方法は効果全ての人に対し同じであるという前提に依存している. より厳密には, 人によって効果に上下がある場合, 集計されることによってその効果が均されてしまい (いわゆる集計バイアス), 非連続な変化の大きさに対する正確な検定ができなくなる, ということである. 個人レベルのデータで効果の差異がないという状況はまれなので, かなり強い仮定であり, 注意が必要である.

多項式近似による推定

カーネル法を用いない場合, カットオフポイントの左右それぞれで, 個別に多項式近似で推定するという方法もある. Lee and Lemiuex (2010) "Regression discontinuity designs in economics" で紹介された方法では, $X$ の多項式 $P(X)=\beta_{0} +\beta_{1}X_{i} + \beta_{2}X_{i}^2 + \cdots + \beta_{p}X_{i}^p$ を用いて,

$\begin{equation} Y_{i} = P(X_{i}) + \tau D_{i} +\varepsilon_{i} \end{equation}$

とする. こちらはノンパラメトリックでないため, 普通に有意検定ができる *1.計量経済学の教科書にも載っているような初歩的な手法の応用なので分かりやすいが, LLR と違い, カットオフポイントから離れた場所のサンプルからも影響を受けるため, 当てはまりの問題が生じることがある (逆に多項式近似でバンド幅を設けた場合, オーバーフィットになる可能性が高い).

どうやって非連続点を見つけるか

RDD を応用できそうな非連続な変化を見つける最もポピュラーな方法は, running variable と結果変数をプロットして目視で判断することである. rdd パッケージは, plot関数に RDD の推定結果を与えるとプロットしてくれる機能を提供してくれる.

また, グラフによる確認が必要なのは $Y$ の非連続点だけではない. running variable $X$ の分布もプロットし, カットオフポイント周辺の分布の形状が不自然でないかも確認する必要がある. もしカットオフポイントの前後どちらかで不自然に度数が多い場合, 個人が意図してどちらのグループに入るか選ぶことができる状況にある可能性があるため, バイアスが発生する. 前回から繰り返し言っているように, ある人に施策がなされるかなされないかは, 本人の意志とは関係なく無作為に決定されなければ正しく推定できない.そのような場合は Sharp RDD ではなく Fuzzy を検討したほうがいいだろう.

まとめ

RDD を応用できそうなポイントは, まずデータをプロットして目視確認する. また, running variable の分布からセレクションバイアスがないかも確認する.
Sharp RDD はカーネル法の一種である local linear regression をカットオフポイント周辺で適用することで推定するのがよい. または, 多項式近似を用いる. 前者は有意性の検定をするのに強い仮定が必要で, 後者は当てはまりの問題があるため, 一長一短である.
Fuzzy RDD は2段階最小二乗法を用いる.
誤差を最小化するバンド幅は Imbens and Kalyanaraman (2012) で提案されている.

長くなったのでRで実演する話は分割します.

[計量経済学][時事ネタ?] 非連続回帰デザイン (RDD) 実践編 - ill-identified diary

参考文献

黒澤昌子. (2005). 積極的労働政策の評価―レビュー. フィナンシャル・レビュー, (77).
西口慶彦 (2014) 『ノンパラメトリック，セミパラメトリック計量経済分析：概観』日本経済学会・春季大会数量的経済・政策分析分科会「ノンパラメトリック, セミパラメトリック計量経済分析」セッション資料
Angrist, J. D., & Pischke, J.-S. (2008). Mostly Harmless Econometrics: An Empricist’s Companion. Princeton University Press. 邦訳大森義明・田中隆一・野口晴子・小原美紀訳『ほとんど無害な計量経済学』エヌティティ出版
Cameron, A. C., & Trivedi, P. K. (2005). Microeconometrics: Methods and Applications . Cambridge University Press.
---, Gelbach, J. B., & Miller, D. L. (2006). Robust inference with multi-way clustering. NBER Working Paper, 0327.
Imbens, G. W., & Lemieux, T. (2008). Regression discontinuity designs: A guide to practice. Journal of Econometrics, 142(2), 615–635. doi:10.1016/j.jeconom.2007.05.001
Imbens, G. W., & Kalyanaraman, K. (2012). Optimal Bandwidth Choice for the Regression Discontinuity Estimator. Review of Economic Studies, 79(3), 933–959. doi:10.1093/restud/rdr043
Li, Q., & Racine, J. S. (2006). Nonparametric Econometrics: Theory and Practice. Princeton University Press.
Fan, J. (1992). Design-adaptive Nonparametric Regression. Journal of American Statistical Association, 87(420), 998–1004.
Hahn, J., Todd, P., & Van der Klaauw, W. (1999). Evaluating the effect of an antidiscrimination law using a regression-discontinuity design. NBER Working Paper, (7131).
Hahn, J., Todd, P., & Van der Klaauw, W. (2001). Identification and estimation of treatment effects with a regression‐discontinuity design. Econometrica, 69(1), 201–209.
McCrary, J. (2008). Manipulation of the running variable in the regression discontinuity design: A density test. Journal of Econometrics, (December). およびそのSTATAコード資料
Wooldridge, J. M. (2010). Econometric Analysis of Cross Section and Panel Data (2nd ed.). The MIT Press.
ill-identified (2014) 『[時事ネタ] 少人数教育は実際どれくらい意味があるのか』
--- (2014) 『非連続回帰デザイン(Regression-Discontinuity Design)』
Drew Dimmery (2013). rdd: Regression Discontinuity Estimation. R package version 0.56. http://CRAN.R-project.org/package=rdd
RjpWiki 『ヒストグラムと密度の推定』

*1:ただし, 標準誤差は不均一分散に対してロバストなものを用いることが推奨される

ill-identified diary

所属組織の見解などとは一切関係なく小難しい話しかしません

[計量経済学] 非連続回帰デザイン (RDD) 理論編

概要