Regression discontinuities design (RDD)　

回帰不連続デザイン
回帰不連続デザイン（regression discontinuity design: RDD）は、介⼊の因果効果を推定するための疑似実験デザイン（quasi experimental design）の 1 つである.

治療群と対照群を分離するランダム・プロセスがなければ、治療群への割り当てが回帰不連続デザイン（RDD）に従う場合、治療効果が識別できる.
これは，ある閾値で，オブザベーションを治療群と対照群に分離する実行（ruuning variable）を必要とする．

Thistlethwaite，Campbell が Journal of Educational Psychology (1960) で初めて報告
特に 2010 年代に⼊ってからは、医療ビッグデータの広がりに併せて主要な医学系雑誌・疫学系雑誌での適⽤事例が増加
データ発生の機序に注目した準実験デザイン
非連続回帰デザインが適用できるデータは限定的
比較的弱い仮定のもとで因果効果を推定できる
Ruuning variable (連続変数)の値 (Z) が特定のカットオフ値よりも高いか低いかによって治療群 (T = 1) か対照群 (T = 0) になるかを決めて、介入効果を推定する
カットオフ値の前後の人々は本質的には同じ
偶然(as if as random)それぞれの群に割り振られた可能性が高い, 異なるのはどちらの群に属するかということだけ
この二つの群の間でアウトカムを比較

RDD type

RDDには2種類の亜種がある．

Sharp RDD
閾値が治療群と対照群を正確に分離する
Running varibaleが閾値を超えると必ず処置
例）20歳を超えると飲酒可能になる
Fazzy RDD
閾値は処理される確率に影響を与える． Running varibaleが閾値を超えると確率的に処置
これは実際には操作変数のアプローチ（LATEを推定する）である

閾値をわずかに下回る個人のoutomce (Y)の値は、欠落した反事実アウトカムに該当する.

推定方法

RDDを推定する方法には3つの方法がある

方法 1.
実行変数の値が閾値に近いサブサンプルを選択する
問題：サンプルが小さいほど標準誤差が大きい
方法2.
より大きなサンプルを選択し、パラメトリックに推定する
問題：これは関数形式と多項式に依存する
方法 3.
閾値に近いサブサンプルを選択して，パラメトリックに推定する拡張：カットオフの左右で異なる機能形態

RDDアプローチでは、いくつかの仮定をテストすることができる.
閾値に近い個体は、治療によって影響を受ける特性を除いて、ほぼ同一である.
治療の前に、治療と対照群の間でアウトカムが異なってはならない.
閾値を示す変数の分布は、このカットオフ値の周りにジャンプがないはずである.

Bandwidth

RDDではバンド（データの幅）の広さを決めて回帰する．例えば，飲酒と死亡の関係を年齢21歳以上，20歳以下を実行変数として使用する場合に，何歳ごとの幅で推定するかだ．

	サンプル数	精度	バイアス
広いBW	多く必要	高い	高い
狭いBW	少ない	低い	低い

Bandwidthの決め方はサンプル数，精度，バイアスのトレードオフになる．適切なBandwidthはサンプルサイズの影響を受けるが明確なルールはない．

事例1

前期の試験が悪かった生徒に補講授業を行うことで後期の試験があがるか（補講の因果効果）を推定したいとしよう．成績がわるかったものに対処しないのは非倫理的なため，RCTは行えない．そこでRDDを用いて考える．

ある大学では補講は，前期の点数が60/100点以下の赤点のもののみ受けるというルールがある．つまり
R: running variable　前期の試験が60点未満かどうか
T: treatment 補講　60点未満は強制的に補講（T=1），60点以上は補講なし（T=0)
Y: outcome 後期試験の点数

となる．

すべての学生は必ず，補講群，補講なし群に分類され，補講をさぼるものや，点数が高いのに補講に参加するものはいないと仮定する．

データの読み込み

まずはデータを読み込む
http://www.ner.takushoku-u.ac.jp/masano/class_material/waseda/keiryo/exam.score.csv

exam <- read.csv("exam.score.csv")
head(round(exam, digits = 0))

散布図

前期試験点数を x 軸、後期試験点数を y 軸にして散布図を描く

exam %>% ggplot(aes(Z,Y))+
  geom_point()+
  xlab("Zi: Test Score (Spring)")+
  ylab("Yi: Test Score (Fall)")+
  ggtitle("Two Test Scores: Spring and Fall")+
  annotate("text",x=42,y=100,label="Treatment Group (T = 1)")+
  annotate("text",x=75,y=100,label="Control Group (T = 0)")+
  geom_vline(xintercept=60, colour="red")   # Zi = 60 に縦の点線を引く

前期試験の点数60点より左が Treatment Group (T = 1), 前期試験の点数60点より右が Control Group (T = 0), 前期試験の点数が高いほど、後期試験の点数が高い.

効果の識別

補講授業に参加したことが学力に与える因果効果が知りたいが，ある学生 i についての因果効果はYi(1)−Yi(0) は片方しか観測できない（片方は反事実）である．このため，集団の平均的な因果的効果を推定する．補講授業に参加すること（Ti）が、潜在的な点数 (Yi(1),Yi(0) ) と完全に独立であるならば，

$[Treatment群のYの平均値] − [Control群のYの平均値] = 平均介入効果 (ATE)$

と表せる．

介入群と対照群に属する学生において、前期試験の値が極めて 60に近い個人同士は比較可能である．つまり，補講授業への参加以外に、後期試験の点数Yi が不連続になる要因はないということである．

しかし、前期試験の点数ZiはTiとYiの両方に影響しているはずである．補講授業に参加することTiが、潜在的な点数Yiと完全に独立ではない.このためZiを条件付けすれば，独立させることができる．

データの要約と回帰式

まずはデータを要約

stargazer(exam,type="text")

## 
## ==============================================================
## Statistic  N   Mean  St. Dev.  Min   Pctl(25) Pctl(75)   Max  
## --------------------------------------------------------------
## Y         200 67.737  10.781  40.663  60.999   74.110  100.000
## Z         200 60.680  12.187  26.891  54.007   67.345  91.624 
## T         200 0.490   0.501     0       0        1        1   
## --------------------------------------------------------------

補講ありとなしの学生に関する回帰式を書く

fit1 <- lm(Y ~ Z, data = exam[exam$Z < 60, ])　#補講あり
fit2 <- lm(Y ~ Z, data = exam[exam$Z > 60, ])　#補講なし

補講を受けた人（前期試験が60点未満）のデータから回帰式を用いて前期試験が60点の人の点数を予測した値をy1
補講を受けない人（前期試験が60点以上）のデータから回帰式を用いて前期試験が60点の人の点数を予測した値をy2

z1.range <- c(60) #　補講あり
y1 <- predict(fit1, newdata = data.frame(Z = z1.range))
y1

##        1 
## 73.10664

z2.range <- c(60) # 60 to max　補講なし
y2 <- predict(fit2, newdata = data.frame(Z = z2.range))
y2

##        1 
## 60.06645

条件付きの平均因果効果は以下で計算できる
$補講を受けない学生の前期試験60点の場合の予測値−補講を受けた学生の前期試験60点での予測値$

y1-y2

##        1 
## 13.04019

つまり，前期試験の点数が Z = 60点という条件下で、補講を受けたことが後期試験の点数に与えた影響（＝条件付平均介入効果 τ ）の推定値がおよそ 13.04点ということである．

散布図からの理解

exam$T<-as.factor(exam$T)
exam %>% ggplot(aes(Z,Y,color=T))+
  geom_point()+
  xlab("Spring Exam Score")+
  ylab("Fall Exam Score")+
  ggtitle("Two Test Scores: Spring and Fall")+
  annotate("text",x=42,y=100,label="Treatment Group (T = 1)")+
  annotate("text",x=75,y=100,label="Control Group (T = 0)")+
  geom_vline(xintercept=60, colour="red")+   # Zi = 60 に縦の点線を引く
  geom_smooth(method="lm",se=TRUE,Color="black",size=1)

## `geom_smooth()` using formula 'y ~ x'

rdd package

ここまで，頑張って線形回帰をで書く，ということをおこなったが，簡単にもできる．詳細は後述する． Rのrddパッケージは、局所線形回帰によって限界平均治療効果(marginal average treatment effect)を推定する方法を含む、回帰不連続デザイン（RDD）の解析として使える.

bw <- with(exam, IKbandwidth(Z, Y, cutpoint = 60))
bw

## [1] 9.609246

rdd_simple <- RDestimate(Y ~ Z, data = exam, cutpoint = 60, bw = bw)
summary(rdd_simple)

## 
## Call:
## RDestimate(formula = Y ~ Z, data = exam, cutpoint = 60, bw = bw)
## 
## Type:
## sharp 
## 
## Estimates:
##            Bandwidth  Observations  Estimate  Std. Error  z value  Pr(>|z|) 
## LATE        9.609     125           -13.724   3.627       -3.783   1.547e-04
## Half-BW     4.805      71            -9.331   5.035       -1.853   6.385e-02
## Double-BW  19.218     173           -13.140   2.591       -5.071   3.967e-07
##               
## LATE       ***
## Half-BW    .  
## Double-BW  ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## F-statistics:
##            F       Num. DoF  Denom. DoF  p        
## LATE        9.954  3         121         1.299e-05
## Half-BW     9.000  3          67         8.683e-05
## Double-BW  11.082  3         169         2.232e-06

LATEとして-13.724と少し値が違う結果がでてきた．これの違いはBandwidthの設定によるものと思われる．

事例2

rdd package

##　シミュレーションデータ
次はシミュレーションデータでやってみよう．まずは変数を設定する．
R: running variable
T: treatment indicator
Y: outcome
X1: Rと相関する連続変数
X3: R，X1と独立する4つのレベルを持つカテゴリー共変量カテゴリカル変数

以下でシミレーションデータセットを作成する． Codeは以下から https://www.jepusto.com/rdd-interactions/#

set.seed(20210105)

simulate_RDD <- function(n = 2000, R = rnorm(n, mean = qnorm(.2))) {
  n <- length(R)
  T <- as.integer(R > 0)
  X1 <- 10 + 0.6 * (R - qnorm(.2)) + rnorm(n, sd = sqrt(1 - 0.6^2))
  X2 <- sample(LETTERS[1:4], n, replace = TRUE, prob = c(0.2, 0.3, 0.35, 0.15))
  Y0 <- 0.4 * R + 0.1 * (X1 - 10) + c(A = 0, B = 0.30, C = 0.40, D = 0.55)[X2] + rnorm(n, sd = 0.9)
  Y1 <- 0.35 + 0.3 * R + 0.18 * (X1 - 10) + c(A = -0.50, B = 0.30, C = 0.20, D = 0.60)[X2] + rnorm(n, sd = 0.9)
  Y <- (1 - T) * Y0 + T * Y1
  data.frame(R, T, X1, X2, Y0, Y1, Y)
}

RD_data <- simulate_RDD(n = 2000)

単純なRDD分析

sharp RDDの主な推定値は、限界平均治療効果（marginal ATE）、すなわち、適格性のしきい値のすぐ近くにある/近くにあるユニットに対する治療割り当ての平均効果である.
共変量X1, x2に依存する治療反応面をシミュレーションしたとしても、MATEを特定するために共変量をコントロールする必要はない.むしろ、実行変数、治療指標、およびそれらの相互作用に関する結果の局所線形回帰を使用することで十分である.

$Y_i = \beta_0 + \beta_1 R_i + \beta_2 T_i + \beta_3 R_i T_i + \epsilon_i$

通常，この回帰は，閾値のある帯域幅内のオブザベーションを用いて，いくつかのカーネルに基づいて定義された重みを用いて推定される．
rddパッケージのデフォルトは、ImbensとKalyanaramanによって提案された式を用いて選択された帯域幅で、三角エッジ・カーネルを使用する.

以下のコードは、共変量を制御せずにMATEを推定するためにrddを使用している.

library(rdd)
bw <- with(RD_data, IKbandwidth(R, Y, cutpoint = 0))
rdd_simple <- RDestimate(Y ~ R, data = RD_data, cutpoint = 0, bw = bw)
summary(rdd_simple)

## 
## Call:
## RDestimate(formula = Y ~ R, data = RD_data, cutpoint = 0, bw = bw)
## 
## Type:
## sharp 
## 
## Estimates:
##            Bandwidth  Observations  Estimate  Std. Error  z value  Pr(>|z|)    
## LATE       1.048      1143          0.3417    0.1303      2.622    0.008742  **
## Half-BW    0.524       615          0.2483    0.1782      1.393    0.163518    
## Double-BW  2.096      1791          0.3251    0.1018      3.193    0.001407  **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## F-statistics:
##            F      Num. DoF  Denom. DoF  p        
## LATE       33.09  3         1139        0.000e+00
## Half-BW    12.50  3          611        1.225e-07
## Double-BW  89.85  3         1787        0.000e+00

1.09の帯域幅を使用して、推定された限界平均治療効果は0.303です.下の図は、不連続性を示しています.　　

library(tidyverse)
RD_data$T<-as.factor(RD_data$T)
RD_data %>% ggplot(aes(R,Y,color=T)) + 
  geom_point()+
  scale_y_continuous(limits=c(-3,3))+
  scale_x_continuous(limits=c(-1.75,1.75))+
  geom_smooth(method="lm",se=TRUE,Color="black",size=1)

## Warning: Ignoring unknown parameters: Color

## `geom_smooth()` using formula 'y ~ x'

## Warning: Removed 350 rows containing non-finite values (stat_smooth).

## Warning: Removed 350 rows containing missing values (geom_point).

共変量を持つRDD

実際には、分析者がモデルに追加の共変量を含めることは非常に一般的である.
これは治療効果の特定には必要ではないが、精度を向上させる目的で有用である.
例えば、Cortes, Goodman, Nomi(2015)は、RDDを用いて、成績の低い9年生を代数学の補講に割り当てた場合の効果を推定している.彼らの主なモデルには、生徒の性別、人種/民族、無料/割引価格のランチの状況などのコントロールが含まれている.

ここには単純な無作為化実験との直接的な類似性がある：基本的な平均値の差は、サンプルの平均的な治療効果の無作為化に偏りのない推定値を提供するが、実際には、共変量を追加したモデルからの推定値を使用することは非常に有用である.

シミュレーションした例に戻るが、次の表は、共変量のどちらか一方、または両方をコントロールした場合のRDestimateによって生成された推定値である.

RD_est <- function(mod, covariates) {
  RD_fit <- RDestimate(as.formula(paste(mod, covariates)), 
                       data = RD_data, cutpoint = 0)
  with(RD_fit, c(est = est[[1]], se = se[1], p = p[1]))
}

covariates <- list("No covariates" = "",
                "X1 only" = "| X1",
                "X2 only" = "| X2",
                "X1 + X2" = "| X1 + X2")

## ------------------------------------------------------------------------------

## You have loaded plyr after dplyr - this is likely to cause problems.
## If you need functions from both plyr and dplyr, please load plyr first, then dplyr:
## library(plyr); library(dplyr)

## ------------------------------------------------------------------------------

## 
## Attaching package: 'plyr'

## The following objects are masked from 'package:dplyr':
## 
##     arrange, count, desc, failwith, id, mutate, rename, summarise,
##     summarize

## The following object is masked from 'package:purrr':
## 
##     compact

アルコールと死亡の関係実際のデータ

とりあえず、RStudioでデータを読み込んでいます.データセットには、回答者の年齢に応じた集計値が含まれています.

データは以下からもらいます http://masteringmetrics.com/wp-content/uploads/2015/01/AEJfigs.dta

データの成型

library(haven)
library(stargazer)
hw10data <- read_dta(file = "AEJfigs.dta")
summary(hw10data)

##     agecell           all           allfitted         internal    
##  Min.   :19.07   Min.   : 88.43   Min.   : 91.71   Min.   :15.98  
##  1st Qu.:20.08   1st Qu.: 92.79   1st Qu.: 93.04   1st Qu.:18.60  
##  Median :21.00   Median : 95.69   Median : 95.18   Median :20.29  
##  Mean   :21.00   Mean   : 95.67   Mean   : 95.80   Mean   :20.29  
##  3rd Qu.:21.92   3rd Qu.: 98.03   3rd Qu.: 97.79   3rd Qu.:21.98  
##  Max.   :22.93   Max.   :105.27   Max.   :102.89   Max.   :24.37  
##                  NA's   :2                         NA's   :2      
##  internalfitted     external     externalfitted     alcohol      
##  Min.   :16.74   Min.   :71.34   Min.   :73.16   Min.   :0.6391  
##  1st Qu.:18.67   1st Qu.:73.04   1st Qu.:74.06   1st Qu.:0.9962  
##  Median :20.54   Median :74.81   Median :74.74   Median :1.2119  
##  Mean   :20.28   Mean   :75.39   Mean   :75.52   Mean   :1.2573  
##  3rd Qu.:21.66   3rd Qu.:77.24   3rd Qu.:76.06   3rd Qu.:1.4701  
##  Max.   :24.04   Max.   :83.33   Max.   :81.78   Max.   :2.5193  
##                  NA's   :2                       NA's   :2       
##  alcoholfitted       homicide     homicidefitted     suicide     
##  Min.   :0.7943   Min.   :14.95   Min.   :16.26   Min.   :10.89  
##  1st Qu.:1.0724   1st Qu.:16.61   1st Qu.:16.54   1st Qu.:11.61  
##  Median :1.2471   Median :16.99   Median :16.99   Median :12.20  
##  Mean   :1.2674   Mean   :16.91   Mean   :16.95   Mean   :12.35  
##  3rd Qu.:1.4455   3rd Qu.:17.29   3rd Qu.:17.25   3rd Qu.:12.82  
##  Max.   :1.8174   Max.   :18.41   Max.   :17.76   Max.   :14.83  
##                   NA's   :2                       NA's   :2      
##  suicidefitted        mva          mvafitted         drugs      
##  Min.   :11.59   Min.   :26.86   Min.   :27.87   Min.   :3.202  
##  1st Qu.:11.61   1st Qu.:30.12   1st Qu.:30.17   1st Qu.:3.755  
##  Median :12.25   Median :31.64   Median :31.73   Median :4.314  
##  Mean   :12.36   Mean   :31.62   Mean   :31.68   Mean   :4.250  
##  3rd Qu.:13.04   3rd Qu.:33.10   3rd Qu.:33.40   3rd Qu.:4.756  
##  Max.   :13.55   Max.   :36.39   Max.   :34.82   Max.   :5.565  
##                  NA's   :2                       NA's   :2      
##   drugsfitted    externalother    externalotherfitted
##  Min.   :3.449   Min.   : 7.973   Min.   : 8.388     
##  1st Qu.:3.769   1st Qu.: 9.149   1st Qu.: 9.347     
##  Median :4.323   Median : 9.561   Median : 9.690     
##  Mean   :4.255   Mean   : 9.599   Mean   : 9.610     
##  3rd Qu.:4.679   3rd Qu.:10.122   3rd Qu.: 9.939     
##  Max.   :5.130   Max.   :11.483   Max.   :10.353     
##                  NA's   :2

hw10data <- subset(hw10data, select = c("agecell", "all", "allfitted",
                                        "internal", "internalfitted",
                                        "external","externalfitted"))

head(hw10data)

hw10data$over21 =  as.numeric(hw10data$agecell>=21)
#length(hw10data$agecell)
#length(hw10data$all[!is.na(hw10data$all)])
hw10data <- na.omit(hw10data)
summary(hw10data)

##     agecell           all           allfitted         internal    
##  Min.   :19.07   Min.   : 88.43   Min.   : 91.71   Min.   :15.98  
##  1st Qu.:20.03   1st Qu.: 92.79   1st Qu.: 93.01   1st Qu.:18.60  
##  Median :21.00   Median : 95.69   Median : 95.18   Median :20.29  
##  Mean   :21.00   Mean   : 95.67   Mean   : 95.71   Mean   :20.29  
##  3rd Qu.:21.97   3rd Qu.: 98.03   3rd Qu.: 97.69   3rd Qu.:21.98  
##  Max.   :22.93   Max.   :105.27   Max.   :102.59   Max.   :24.37  
##  internalfitted     external     externalfitted      over21   
##  Min.   :16.74   Min.   :71.34   Min.   :73.23   Min.   :0.0  
##  1st Qu.:18.61   1st Qu.:73.04   1st Qu.:74.07   1st Qu.:0.0  
##  Median :20.51   Median :74.81   Median :74.74   Median :0.5  
##  Mean   :20.27   Mean   :75.39   Mean   :75.44   Mean   :0.5  
##  3rd Qu.:21.72   3rd Qu.:77.24   3rd Qu.:75.89   3rd Qu.:1.0  
##  Max.   :24.04   Max.   :83.33   Max.   :81.48   Max.   :1.0

carpenter_dobkin_2009 <- as.data.frame(hw10data)

では、21歳で発生する閾値（＝最低飲酒年齢）を見てみましょう.
21歳以降の死亡率には顕著な跳躍がありそうだ．

carpenter_dobkin_2009 %>% 
  ggplot(aes(x = agecell, y = all)) + 
  geom_point() +
  geom_vline(xintercept = 21, color = "red", size = 1, linetype = "dashed") + 
  annotate("text", x = 20.4, y = 105, label = "Minimum Drinking Age",
           size=4) +
  labs(y = "Mortality rate (per 100.000)",
       x = "Age (binned)")

最初にある個人がカットオフ値より下か上かを示すダミー変数（閾値）を計算する．

閾値のダミー変数は、21年のカットオフ以下のオブザベーションではゼロに等しく、21年のカットオフより上のオブザベーションでは1に等しいように作成．

100,000人あたりのすべての死亡（すべて）を閾値ダミーと年齢の閾値（21歳）を中心とする回答者の年齢に回帰させるために、関数lm()で線形モデルを指定する．

なお%$%はデータフレームから直接変数を渡す演算子でmagrittrから使える．

これは、各年齢ビンからカットオフを差し引くことで、関数I()で行う．

lm_same_slope <- carpenter_dobkin_2009 %>%
  mutate(threshold = ifelse(agecell >= 21, 1, 0)) %$%
  lm(all ~ threshold + I(agecell - 21))
summary(lm_same_slope)

## 
## Call:
## lm(formula = all ~ threshold + I(agecell - 21))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -5.0559 -1.8483  0.1149  1.4909  5.8043 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      91.8414     0.8050 114.083  < 2e-16 ***
## threshold         7.6627     1.4403   5.320 3.15e-06 ***
## I(agecell - 21)  -0.9747     0.6325  -1.541     0.13    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.493 on 45 degrees of freedom
## Multiple R-squared:  0.5946, Adjusted R-squared:  0.5765 
## F-statistic: 32.99 on 2 and 45 DF,  p-value: 1.508e-09

ダミーの前変数閾値の係数は平均的な治療効果である.
平均では、最低飲酒年齢に達した人の100.000人当たりの死亡率は7.66ポイント高い.

RDD の適用に関連した様々な関数を含む R パッケージ rddtoolsがあったがVersionの関係で4以降は対応していない．https://github.com/bquast/rddtools rddを使用する．

散布図を使って、閾値の両側で同じ傾きを示す回帰の適合線を描きます.

carpenter_dobkin_2009 %>%
  select(agecell, all) %>%
  mutate(threshold = as.factor(ifelse(agecell >= 21, 1, 0))) %>%
  ggplot(aes(x = agecell, y = all)) +
  geom_point(aes(color = threshold)) +
  geom_smooth(method = "lm", se = FALSE) +
  scale_color_brewer(palette = "Accent") +
  guides(color = FALSE) +
  geom_vline(xintercept = 21, color = "red",
    size = 1, linetype = "dashed") +
  labs(y = "Mortality rate (per 100.000)",
       x = "Age (binned)")

## `geom_smooth()` using formula 'y ~ x'

Estimation strategy: different slopes

関数 lm()では，カットオフ値を中心とした閾値ダミーと年齢との間の相互作用を指定することで，これを達成することができます．

lm_different_slope <- carpenter_dobkin_2009 %>%
  mutate(threshold = ifelse(agecell >= 21, 1, 0)) %$%
  lm(all ~ threshold + I(agecell - 21) + threshold:I(agecell - 21))

summary(lm_different_slope)

## 
## Call:
## lm(formula = all ~ threshold + I(agecell - 21) + threshold:I(agecell - 
##     21))
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -4.368 -1.787  0.117  1.108  5.341 
## 
## Coefficients:
##                           Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                93.6184     0.9325 100.399  < 2e-16 ***
## threshold                   7.6627     1.3187   5.811  6.4e-07 ***
## I(agecell - 21)             0.8270     0.8189   1.010  0.31809    
## threshold:I(agecell - 21)  -3.6034     1.1581  -3.111  0.00327 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.283 on 44 degrees of freedom
## Multiple R-squared:  0.6677, Adjusted R-squared:  0.645 
## F-statistic: 29.47 on 3 and 44 DF,  p-value: 1.325e-10

このアプローチは治療効果の解釈を変えるものではない! 平均して、最低飲酒年齢に達した人の10万人当たりの死亡率は7.66ポイント高い.

一応rddtoolsが復活したときようのコード

rdd_data(y = carpenter_dobkin_2009$all, 
         x = carpenter_dobkin_2009$agecell, 
         cutpoint = 21) %>% 
  rdd_reg_lm(slope = "separate") %>% 
  summary()

散布図で傾きの違いを見てみましょう.カットオフの右側の傾きが負であるのに対し、左側の傾きは正です.

carpenter_dobkin_2009 %>%
  select(agecell, all) %>%
  mutate(threshold = as.factor(ifelse(agecell >= 21, 1, 0))) %>%
  ggplot(aes(x = agecell, y = all, color = threshold)) +
  geom_point() +
  geom_smooth(method = "lm", se = FALSE) +
  scale_color_brewer(palette = "Accent") +
  guides(color = FALSE) +
  geom_vline(xintercept = 21, color = "red",
    size = 1, linetype = "dashed") +
  labs(y = "Mortality rate (per 100.000)",
       x = "Age (binned)")

## `geom_smooth()` using formula 'y ~ x'

RDD を適用する際には、機能形態の仕様に特に注意が必要です.

以下、私は年齢と100.000人あたりの死亡率（全）との間の二次関係をモデル化しています.二次項は、関数I()を介して式に入ります.前節と同様に、カットオフ付近で異なる傾きを使用しています.

lm_quadratic <- carpenter_dobkin_2009 %>% 
mutate(threshold = ifelse(agecell >= 21, 1, 0)) %$% 
  lm(all ~ threshold + I(agecell - 21) + I((agecell -21)^2) + threshold:I(agecell - 21) +
       threshold:I((agecell - 21)^2))

summary(lm_quadratic)

## 
## Call:
## lm(formula = all ~ threshold + I(agecell - 21) + I((agecell - 
##     21)^2) + threshold:I(agecell - 21) + threshold:I((agecell - 
##     21)^2))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.3343 -1.3946  0.1849  1.2848  5.0817 
## 
## Coefficients:
##                               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                    93.0729     1.4038  66.301  < 2e-16 ***
## threshold                       9.5478     1.9853   4.809 1.97e-05 ***
## I(agecell - 21)                -0.8306     3.2901  -0.252    0.802    
## I((agecell - 21)^2)            -0.8403     1.6153  -0.520    0.606    
## threshold:I(agecell - 21)      -6.0170     4.6529  -1.293    0.203    
## threshold:I((agecell - 21)^2)   2.9042     2.2843   1.271    0.211    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.285 on 42 degrees of freedom
## Multiple R-squared:  0.6821, Adjusted R-squared:  0.6442 
## F-statistic: 18.02 on 5 and 42 DF,  p-value: 1.624e-09

carpenter_dobkin_2009 %>%
  select(agecell, all) %>%
  mutate(threshold = as.factor(ifelse(agecell >= 21, 1, 0))) %>%
  ggplot(aes(x = agecell, y = all, color = threshold)) +
  geom_point() +
  geom_smooth(method = "lm",
              formula = y ~ x + I(x ^ 2),
              se = FALSE) +
  scale_color_brewer(palette = "Accent") +
  guides(color = FALSE) +
  geom_vline(xintercept = 21, color = "red",
    size = 1, linetype = "dashed") +
  labs(y = "Mortality rate (per 100.000)",
       x = "Age (binned)")

平均して、最低飲酒年齢に達した人の10万人当たりの死亡率は9.55ポイント高くなっている.

library(rdd)
bw <- with(carpenter_dobkin_2009, IKbandwidth(agecell, all, cutpoint = 21))
rdd_simple <- RDestimate(all~ agecell, data = carpenter_dobkin_2009, cutpoint = 21, bw = 2)
summary(rdd_simple)

## 
## Call:
## RDestimate(formula = all ~ agecell, data = carpenter_dobkin_2009, 
##     cutpoint = 21, bw = 2)
## 
## Type:
## sharp 
## 
## Estimates:
##            Bandwidth  Observations  Estimate  Std. Error  z value  Pr(>|z|) 
## LATE       2          48            8.381     1.345       6.232    4.603e-10
## Half-BW    1          24            9.700     1.932       5.022    5.112e-07
## Double-BW  4          48            7.889     1.274       6.192    5.940e-10
##               
## LATE       ***
## Half-BW    ***
## Double-BW  ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## F-statistics:
##            F      Num. DoF  Denom. DoF  p        
## LATE       36.57  3         44          1.024e-11
## Half-BW    26.74  3         20          6.711e-07
## Double-BW  31.85  3         44          8.449e-11

References

Carpenter, Christopher, and Carlos Dobkin. 2009. “The Effect of Alcohol Consumption on Mortality: Regression Discontinuity Evidence from the Minimum Drinking Age.” American Economic Journal: Applied Economics 1 (1): 164–82. https://doi.org/10.1257/app.1.1.164.

https://rpubs.com/phle/r_tutorial_regression_discontinuity_design

https://www.jepusto.com/rdd-interactions/#

http://www.ner.takushoku-u.ac.jp/masano/class_material/waseda/keiryo/16_rdd.html

RDD practice

Ryohei Yamamoto

2021/01/05

Regression discontinuities design (RDD)

RDD type

推定方法

Bandwidth

事例1

データの読み込み

散布図

効果の識別

データの要約と回帰式

散布図からの理解

rdd package

事例2

rdd package

単純なRDD分析

共変量を持つRDD

アルコールと死亡の関係実際のデータ

データの成型

Estimation strategy: different slopes

References

RDD practice

Ryohei Yamamoto

2021/01/05

Regression discontinuities design (RDD)

RDD type

推定方法

Bandwidth

事例1

データの読み込み

散布図

効果の識別

データの要約と回帰式

散布図からの理解

rdd package

事例2

rdd package

単純なRDD分析

共変量を持つRDD

アルコールと死亡の関係 実際のデータ

データの成型

Estimation strategy: different slopes

References

Regression discontinuities design (RDD)　

アルコールと死亡の関係実際のデータ