1 Data

Angrist and Pischke (2014) Mastering ’Metrics, Princeton University Press. が提供する 法廷飲酒年齢のRD分析のデータとコード を用いる. このデータは Carpenter and Dobkin (AEJ: Applied 2009) に基づく.

注:このページは Carpenter and Dobkin (2009) の reproduction/replication というわけではない.

上のサイトでは AEJfigs.dta というファイル名. Moodle には MMch4_RD.dta というファイル名でアップロードしている.

主な変数

  • outcome
    • all: mortality rate from all causes (per 100,000)
  • running variable
    • agecell: age cell = 19.07, 19.15, 19.23, 19.31, …, 22.93

つまり,個人レベルのデータではなく,年齢の cell を単位として集計されたデータである.

setwd("c://ws_stat")
library(tidyverse)
mm <- haven::read_dta("MMch4_RD.dta")
nrow(mm)
## [1] 50
mm <- mm %>%
  mutate(
    age = agecell - 21,  # centering 
    over21 = ifelse(agecell >= 21, 1, 0),
  )

1.1 Reproduction/Replication by others

2 Plot discontinuity

2.1 Naive fitting

閾値(21歳)の前後で切片が「ジャンプ」するように線形モデルを推定.

naive_linear <- lm(all ~ age + over21, data = mm)
naive_linear$coefficients[3]  # size of `jump` 
##   over21 
## 7.662709
mm$allfit <- predict(naive_linear, newdata = mm)  # 回帰直線を引くために fitted value を計算
ggplot(mm, aes(x = agecell, y = all)) +
  geom_point() +
  geom_line(aes(y = allfit)) 

縦軸(10万人当たり死亡率)が 8 ポイント近く「ジャンプ」している.

これがもっとも単純な RD 分析.

2.2 Flexible fitting

閾値(21歳)の前後で別々の滑らかな曲線を当てはめる.

  • LOESS (locally estimated scatterplot smoothing) と呼ばれる局所回帰の手法を使用してデータに当てはまる曲線を推定しているが(引数 method = loess),この演習では理解する必要はない
  • aescol として閾値の前と後を区別する factor を指定することで,閾値前後で別々の曲線が推定できる
ggplot(mm, aes(x = agecell, y = all, col = factor(over21))) + 
  geom_point() + 
  geom_smooth(method = loess, se = FALSE) +
  guides(color = FALSE)  # legend off 

縦軸(10万人当たり死亡率)が 10 ポイント近く「ジャンプ」している.

cf. Carpenter and Dobkin (2009) Figure 3

2.3 rdrobust::rdplot 関数を利用

RDD 分析は rdrobust パッケージ を用いて行うことができる.

  • 引数 c で cutoff を指定
  • デフォルトでは4次の多項式 (引数 p = 4) でフィッティング
  • デフォルトでは bin の数が小さすぎるので,nbins 引数で調整
# install.packages("rdrobust")
library(rdrobust)
rdplot(y = mm$all, x = mm$agecell, c = 21, nbins = 50)

3 Estimate the gap

3.1 Using linear/quadratic model

Baseline: Linear

\[ y = \beta_0 + \beta_1 1(\mbox{Age} \ge 21) + \beta_2 \mbox{Age} + \beta_3 1(\mbox{Age} \ge 21) \times \mbox{Age} +u \\ = \begin{cases} (\beta_0 + \beta_1) + (\beta_2 + \beta_3) \mbox{Age} +u \quad \mbox{if Age} \ge 21 \\ \beta_0 + \beta_2 \mbox{Age} +u \quad \mbox{if Age} < 21 \end{cases} \]

  • 不均一分散に頑健な標準誤差を計算するために fixest::feols を使う.
  • 興味があるのは over21 変数の係数.これが閾値(21歳)前後の「ジャンプ」すなわち因果効果を捕捉している.
  • over21:age は,age 変数の係数が閾値の前後で異なってよいことを表現するためのもの.21歳未満では age の係数,21歳以上では age の係数と over21:age の係数を足したもの.
library(fixest)
feols_linear <- feols(all ~ over21 + age + over21:age, data = mm, vcov = "hetero")
feols_quad <- feols(all ~ over21 + age + I(age^2) + over21 + over21:age + over21:I(age^2), 
                    data = mm, vcov = "hetero")
library(modelsummary)
modelsummary(list(feols_linear, feols_quad), 
             gof_omit = "AIC|BIC|RMSE", stars = TRUE)
 (1)   (2)
(Intercept) 93.618*** 93.073***
(0.628) (0.780)
over21 7.663*** 9.548***
(1.273) (1.830)
age 0.827 -0.831
(0.721) (2.850)
over21 × age -3.603** -6.017
(1.124) (4.528)
I(age^2) -0.840
(1.541)
over21 × I(age^2) 2.904
(2.257)
Num.Obs. 48 48
R2 0.668 0.682
R2 Adj. 0.645 0.644
Std.Errors Heteroskedasticity-robust Heteroskedasticity-robust
+ p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001

3.2 rdrobust::rdrobust 関数を利用

何らかのカーネル関数を用いて推定する. デフォルトは三角カーネル(kernel = "triangular").

summary(rdrobust(y = mm$all, x = mm$agecell, c = 21))
## Sharp RD estimates using local polynomial regression.
## 
## Number of Obs.                   48
## BW type                       mserd
## Kernel                   Triangular
## VCE method                       NN
## 
## Number of Obs.                   24           24
## Eff. Number of Obs.               6            6
## Order est. (p)                    1            1
## Order bias  (q)                   2            2
## BW est. (h)                   0.493        0.493
## BW bias (b)                   0.780        0.780
## rho (h/b)                     0.632        0.632
## Unique Obs.                      24           24
## 
## =============================================================================
##         Method     Coef. Std. Err.         z     P>|z|      [ 95% C.I. ]       
## =============================================================================
##   Conventional     9.595     3.591     2.672     0.008     [2.557 , 16.633]    
##         Robust         -         -     2.205     0.027     [1.077 , 18.300]    
## =============================================================================

9.6 ポイントの「ジャンプ」が認められる.

3.3 Honest confidence interval (おまけ)

定式化によって推定値は大きく影響を受ける可能性があり,定式化の誤りを考慮して推定値の信頼区間を計算する手法が開発されている.

Kolesár and Rothe (AER 2018) に基づく “honest” な信頼区間を計算する.

この計算は RDHonest::RDHonest で計算できるが,この パッケージ RDHonest は CRAN には含まれないため,パッケージ作成者個人の GitHub からダウンロードする. GitHub からパッケージをダウンロードするためには remotes パッケージが必要.

# install.packages("remotes")
# remotes::install_github("kolesarm/RDHonest")
library(RDHonest)
RDHonest(all ~ agecell, data = mm, cutoff = 21)
## 
## Call:
## RDHonest(formula = all ~ agecell, data = mm, cutoff = 21)
## 
## Inference for Sharp RD parameter (using Holder class), confidence level 95%:
##              Estimate Std. Error Maximum Bias  Confidence Interval
## I(agecell>0) 11.43565   4.903254     1.981793 (1.090396,  21.7809)
## 
## Onesided CIs:  (-Inf, 21.48258), ( 1.38872, Inf)
## Number of effective observations: 6.757664
## Maximal leverage for sharp RD parameter: 0.4602486
## Smoothness constant M: 156.7686
## P-value: 0.03002816 
## 
## Based on local regression with bandwidth: 0.3189381, kernel: triangular
## Regression coefficients:
##         I(agecell>0)  I(agecell>0):agecell           (Intercept)  
##               11.436               -43.082                94.592  
##              agecell  
##                9.042  
## 2 observations with missing values dropped

4 McCrary test (おまけ)

RDD の識別の重要な仮定は running variable の閾値前後における counterfactual outcome の条件付き期待値の連続性である. しかしながら,たとえば individual 自身が running variable を「操作」できる場合 (sorting),この仮定は満たされない場合がある.

そこで,McCrary (JE 2008) が提案した密度検定を行いこの仮定を検証する. 閾値の前後で running variable が不連続に多くなっていれば何らかのメカニズムによって sorting が生じていることになり,連続性の仮定を満たさない.

検定には rdd パッケージの DCdensity 関数を用いる.

注:Carpenter and Dobkin (2009) の running variable は cell であるためこの検定方法は不適であるし,そもそも年齢を直接操作することはできない.

以下にパッケージの help に書かれているサンプルコードを示す. Running variable running_variable が処置を受ける閾値を 0 とし,その前後で density が異なるようにデータを生成する.

# install.packages("rdd")
library(rdd)
set.seed(2)
running_variable <- runif(1000, -1, 1)  # running variable 
running_variable <- running_variable + 2 * (runif(1000, -1, 1) > 0 & running_variable < 0)
hist(running_variable, breaks = 100, freq = FALSE)
abline(v = 0, col = 2, lwd = 2, lty = 2)  # red dashed vertical line 

rdd::DCdensity(runvar = running variable, cutpoint = cutoff) で検定を行う.

出力されるのは「Sorting は生じていない」という帰無仮説のもとでの p 値. p 値が有意水準(5%など)を下回る場合,帰無仮説は棄却され,「Sorting が生じているだろう」と結論付けられる.

rdd::DCdensity(runvar = running_variable, cutpoint = 0)

## [1] 0.004311685

縦軸を相対頻度(密度)に設定したヒストグラムと比較すると,この検定が何を行っているかをイメージできるだろう.

.