回帰不連続デザイン
UPDATE: 2022-10-28 05:31:59
はじめに
このノートではRを使って因果推論に関する基礎的な理論から分析を実行するための方法をまとめている。ここでは回帰不連続デザインについて、下記の矢内先生の資料(資料のライセンスに基づく)を参考にしながら、回帰不連続デザインをまとめている。私の記事は分析の方法メインなので、矢内先生の資料を読む方が数理も学べて良いと思う。
上記の資料のライセンスは下記の通り。
The text of this work is licensed under the Creative Commons Attribution 4.0 International License. The R Code in this work is licensed under the MIT License.
回帰不連続デザインの数理
回帰不連続デザインの数理については下記が詳しいので、ここでは扱わない。
回帰不連続デザイン
回帰不連続デザインとカーネル密度推定
回帰不連続デザインとバンド幅
準備
必要なライブラリを読み込んでおく。
library(tidyverse)
library(broom)
library(rdd)
library(raincloudplots)
library(rdrobust)ここで使用するサンプルデータは、参考資料と同じく法定飲酒年齢の引き下げと21歳の誕生日付近の死亡者数に関する研究のデータ。
このデータにはMLDAが21歳の場合のみが含まれており、このデータの1レコードは「誕生日によって区切られた年齢グループ」とのことで、１レコードが19歳の誕生日からおよそ30日単位のブロックになっているデータ。
# data: http://masteringmetrics.com/wp-content/uploads/2015/01/AEJfigs.dta
# agecell: 19歳の誕生日からおよそ30日単位のブロック(0-30/31-60/61-90/...)のようなイメージ
# all: 死者数
df <- read_csv('~/Desktop/AEJfigs.csv', show_col_types = FALSE)
df## # A tibble: 48 × 2
## agecell all
## <dbl> <dbl>
## 1 19.1 92.8
## 2 19.2 95.1
## 3 19.2 92.1
## 4 19.3 88.4
## 5 19.4 88.7
## 6 19.5 90.2
## 7 19.6 96.2
## 8 19.6 89.6
## 9 19.7 93.4
## 10 19.8 90.9
## # … with 38 more rows
## # ℹ Use `print(n = ...)` to see more rows回帰不連続デザインの解説で、(個人的に)よく見かける変数変換を行う。変数を中心化することで、基準値からのプラス・マイナスのズレで表現できる。
df <- df %>%
mutate(age = agecell - 21, # 21歳が基準
over = as.integer(agecell > 21), # 21歳以上を表す処置変数
age_sq = age^2)
df## # A tibble: 48 × 5
## agecell all age over age_sq
## <dbl> <dbl> <dbl> <int> <dbl>
## 1 19.1 92.8 -1.93 0 3.73
## 2 19.2 95.1 -1.85 0 3.42
## 3 19.2 92.1 -1.77 0 3.12
## 4 19.3 88.4 -1.68 0 2.84
## 5 19.4 88.7 -1.60 0 2.57
## 6 19.5 90.2 -1.52 0 2.31
## 7 19.6 96.2 -1.44 0 2.07
## 8 19.6 89.6 -1.36 0 1.84
## 9 19.7 93.4 -1.27 0 1.62
## 10 19.8 90.9 -1.19 0 1.42
## # … with 38 more rows
## # ℹ Use `print(n = ...)` to see more rows21という値が0に位置するので、イメージとしてはこのようになる。
df_1x1 <- data_1x1(
array_1 = df %>% pull(agecell),
array_2 = df %>% pull(age),
jit_distance = .09,
jit_seed = 1)
raincloud_1x1_repmes(
data = df_1x1,
colors = (c('dodgerblue', 'darkorange')),
fills = (c('dodgerblue', 'darkorange')),
line_color = 'gray',
# line_alpha = .3,
# size = 1,
# alpha = .5,
align_clouds = FALSE) + theme_bw()
パラメトリックRD
パラメトリックRDでは、回帰式の関数形を仮定し、割当変数(running variable)の最小から最大までの範囲について回帰式の予測値を計算することで局所平均処置効果を計算する。 つまり、回帰不連続デザインでは回帰モデルの関数形を正しく特定する必要がある。誤った関数形で効果があったように見えても、正しい関数形では効果はないかもしれないため。ここでは参考資料に従い、4つのモデルを作成する。
models <- alist(
model1 = all ~ over + age,
model2 = all ~ over * age,
model3 = all ~ over + age + age_sq,
model4 = all ~ over * age + over * age_sq) %>%
enframe(name = "model_name", value = "formula")
models## # A tibble: 4 × 2
## model_name formula
## <chr> <list>
## 1 model1 <language>
## 2 model2 <language>
## 3 model3 <language>
## 4 model4 <language>用意したモデルをmap()を利用して、一気に回帰係数を推定し、予測、結果の整形を行う。map()が理解し難い場合はfor-loopで4回繰り返したと思えばOK。
rdd <- models %>%
mutate(
model = map(.x = formula, .f = function(x){lm(formula = x, data = df)}),
pred = map(.x = model, .f = function(x){predict(x)}),
result = map(.x = model, .f = function(x){tidy(x)})
)
rdd## # A tibble: 4 × 5
## model_name formula model pred result
## <chr> <list> <list> <list> <list>
## 1 model1 <language> <lm> <dbl [48]> <tibble [3 × 5]>
## 2 model2 <language> <lm> <dbl [48]> <tibble [4 × 5]>
## 3 model3 <language> <lm> <dbl [48]> <tibble [4 × 5]>
## 4 model4 <language> <lm> <dbl [48]> <tibble [6 × 5]>1つ目のモデルの結果を可視化すると、基準点を境に死亡数が高くなっていることがわかる。
df %>%
mutate(pred = rdd$pred[[1]]) %>%
ggplot(aes(x = agecell, color = as.factor(over))) +
geom_point(aes(y = all)) +
geom_line(aes(y = pred)) +
geom_vline(xintercept = 21, color = "gray") +
theme_bw() +
theme(legend.position = "none") 
この場合の局所平均処置効果はoverの係数に該当するので、回帰係数を取り出すと、7.66ほど10万人あたりの死亡数が高くなっていることになる。
rdd %>%
unnest(cols = result) %>%
filter(model_name == "model1", term == "over") %>%
select(term, estimate, std.error)## # A tibble: 1 × 3
## term estimate std.error
## <chr> <dbl> <dbl>
## 1 over 7.66 1.44他のモデルも可視化して、基準点でジャンプがあるか確認してみると、すべての関数形でジャンプが確認できるため、21歳を基準に10万人あたりの死亡数が高くなっていると思われる。
df %>%
select(agecell, over, all) %>%
mutate(
model1 = rdd$pred[[1]],
model2 = rdd$pred[[2]],
model3 = rdd$pred[[3]],
model4 = rdd$pred[[4]],
) %>%
pivot_longer(
cols = model1:model4,
names_to = "model",
values_to = "pred"
) %>%
ggplot(aes(x = agecell, color = as.factor(over))) +
geom_point(aes(y = all)) +
geom_line(aes(y = pred)) +
geom_vline(xintercept = 21, color = "gray") +
facet_wrap( ~ model) +
theme_bw() +
theme(legend.position = "none") 
推定した4つのモデルの結果を比べると、4つ目のモデルの局所平均処置効果は9.54となっており、関数形の違いで結果が異なってしまうことが回帰不連続デザインで起こりうる。つまり、誤った関数形を仮定すると、誤った推定結果を得ることにつながる。
rdd %>%
unnest(cols = result) %>%
filter(term == "over") %>%
select(term, estimate, std.error)## # A tibble: 4 × 3
## term estimate std.error
## <chr> <dbl> <dbl>
## 1 over 7.66 1.44
## 2 over 7.66 1.32
## 3 over 7.66 1.34
## 4 over 9.55 1.99ノンパラメトリックRD
ノンパラメトリックRDでは回帰式の関数形を考えず、推定対象となる範囲(バンド幅)を基準値の周辺に限定し、処置群と統制群の平均値を比較する。つまり、バンド幅の選び方によって推定値が変わる。ノンパラメトリックRDではバンド幅の選び方が重要。ただ、基準値付近の値を利用するので、サンプルサイズが大きい場合、バンド幅が狭いほどバイアスは小さくなる。
暫定的にバンド幅を下記の通り0.25刻みで用意し、
df2 <- tibble(bandwidth = seq(from = 0.5, to = 1.5, by = 0.25))
df2## # A tibble: 5 × 1
## bandwidth
## <dbl>
## 1 0.5
## 2 0.75
## 3 1
## 4 1.25
## 5 1.5このバンド幅でデータをフィルタリングする。
df2 %>%
mutate(
banding_data = map(.x = bandwidth, .f = function(x){df %>% filter(age > -1*x, age < x)})
) %>% unnest(banding_data) %>%
group_by(bandwidth) %>%
summarise(min = min(agecell),
max = max(agecell),
diff = max(agecell) - min(agecell),
sanplesize = n()
)## # A tibble: 5 × 5
## bandwidth min max diff sanplesize
## <dbl> <dbl> <dbl> <dbl> <int>
## 1 0.5 20.5 21.5 0.904 12
## 2 0.75 20.3 21.7 1.40 18
## 3 1 20.1 21.9 1.89 24
## 4 1.25 19.8 22.2 2.38 30
## 5 1.5 19.6 22.4 2.88 36フィルタリングしたデータをもとに回帰係数を計算すると、バンド幅が小さいほど標準誤差が大きく、バンド幅が大きくなると標準誤差は小さくなるし、バンド幅によって回帰係数の推定値も変化する。
df3 <- df2 %>%
mutate(
# formula = "all ~ over + age"として、map2で渡すときにas.formulaでも良い
formula = alist(all ~ over + age),
banding_data = map(.x = bandwidth, .f = function(x){df %>% filter(age > -1*x, age < x)}),
model = map2(.x = banding_data, .y = formula,
.f = function(x, y){lm(formula = y, data = x)}),
result = map(.x = model, .f = function(x){tidy(x)}),
pred = map2(.x = banding_data, .y = model,
.f = function(x, y){predict(y, newdata = x)}),
)
df3 %>%
unnest(result) %>%
filter(term == "over") %>%
select(bandwidth, term, estimate, std.error)## # A tibble: 5 × 4
## bandwidth term estimate std.error
## <dbl> <chr> <dbl> <dbl>
## 1 0.5 over 8.88 2.45
## 2 0.75 over 9.85 2.20
## 3 1 over 9.75 1.93
## 4 1.25 over 9.44 1.70
## 5 1.5 over 8.70 1.55各バンド幅でフィルタリングされたデータや予測値を可視化すると各モデルの設定がわかりやすい。
df3 %>%
select(bandwidth, banding_data, pred) %>%
unnest(c(banding_data, pred)) %>%
ggplot(aes(x = agecell, color = as.factor(over))) +
geom_point(aes(y = all)) +
geom_line(aes(y = pred)) +
geom_vline(xintercept = 21, color = "gray") +
facet_wrap( ~ bandwidth) +
theme_bw() +
theme(legend.position = "none") 
rdrobustパッケージ
回帰不連続デザインを行えるRのrdrobustパッケージを利用して同じ分析を行う。バンド幅を0.5で行うと、先程と同じく、基準値前後で6サンプル(Eff.
Number of
Obs.)が計算対象になっており、Coef.=8.878局所局所平均処置効果を表している。
rdd <- rdrobust(df$all, df$age, c = 0, h = 0.5 , kernel = "uniform")
summary(rdd)## Sharp RD estimates using local polynomial regression.
##
## Number of Obs. 48
## BW type Manual
## Kernel Uniform
## 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.500 0.500
## BW bias (b) 0.500 0.500
## rho (h/b) 1.000 1.000
## Unique Obs. 24 24
##
## =============================================================================
## Method Coef. Std. Err. z P>|z| [ 95% C.I. ]
## =============================================================================
## Conventional 8.878 3.029 2.931 0.003 [2.942 , 14.815]
## Robust - - 2.055 0.040 [0.504 , 21.338]
## =============================================================================カーネル関数とバンド幅
rdrobustパッケージでは、カーネル関数を選択できる。回帰不連続デザインでのカーネル関数の役割はバンド幅のデータに対する重み付けである。下記のようなカーネル関数を利用して重み付けを行う。uniformは重みが1なので、手計算で行ってきた数値と一致する。デフォルトでは三角形の関数が利用され、線形で重みが付けられる。
layout(matrix(1:4, 2, 2, byrow = TRUE))
plot(density(0, bw = 1, kernel = "gaussian"), main = "gaussian")
plot(density(0, bw = 1, kernel = "rectangular"), main = "uniform")
plot(density(0, bw = 1, kernel = "triangular"), main = "triangular")
plot(density(0, bw = 1, kernel = "epanechnikov"), main = "epanechnikov")
バンド幅の推奨される設定方法として、Imbens and kalyanaramanが考案したIKバンドもしくは、Coverage Error rate Regression Discontinuity(CER)がある。詳細は下記が詳しい。
IKバンドはrdbwselect_2014()のhで計算でき、
list(
rdbwselect_2014(df$all, df$age, c = 0, bwselect = "IK")
)## [[1]]
## Call:
## rdbwselect_2014(y = df$all, x = df$age, c = 0, bwselect = "IK")
##
## BW Selector IK
## Number of Obs 48
## NN Matches 3
## Kernel Type Triangular
##
## Left Right
## Number of Obs 24 24
## Order Loc Poly (p) 1 1
## Order Bias (q) 2 2
##
## h b
## [1,] 1.645225 0.9166231CERバンド幅を使用するときは、bwselect = "cerrd"と指定する。
rdd2 <- rdrobust(df$all, df$age, c = 0, bwselect = "cerrd")
summary(rdd2)## Sharp RD estimates using local polynomial regression.
##
## Number of Obs. 48
## BW type cerrd
## Kernel Triangular
## VCE method NN
##
## Number of Obs. 24 24
## Eff. Number of Obs. 5 5
## Order est. (p) 1 1
## Order bias (q) 2 2
## BW est. (h) 0.406 0.406
## BW bias (b) 0.780 0.780
## rho (h/b) 0.521 0.521
## Unique Obs. 24 24
##
## =============================================================================
## Method Coef. Std. Err. z P>|z| [ 95% C.I. ]
## =============================================================================
## Conventional 10.179 4.106 2.479 0.013 [2.131 , 18.226]
## Robust - - 2.171 0.030 [0.994 , 19.495]
## =============================================================================RDプロット
回帰不連続デザインではRDプロットと呼ばれる可視化方法がある。これはx軸をビニングしてサンプルを分割し、そのビンでの平均値を計算し、多項式を当てはめることで、分布形を調べることができる。
rdplt <- rdplot(df$all, df$age, c = 0, p = 4)
RDプロットの設定詳細はsumamry()で確認でき、引数を調整することで、設定は変更可能である。
summary(rdplt)## Call: rdplot
##
## Number of Obs. 48
## Kernel Uniform
##
## Number of Obs. 24 24
## Eff. Number of Obs. 24 24
## Order poly. fit (p) 4 4
## BW poly. fit (h) 1.932 1.932
## Number of bins scale 1 1
##
## Bins Selected 4 5
## Average Bin Length 0.483 0.386
## Median Bin Length 0.483 0.386
##
## IMSE-optimal bins 8 5
## Mimicking Variance bins 4 5
##
## Relative to IMSE-optimal:
## Implied scale 0.500 1.000
## WIMSE variance weight 0.889 0.500
## WIMSE bias weight 0.111 0.500連続性の仮定
回帰不連続デザインを実行するにあたり条件付き回帰関数の連続性を満たす必要がある。グラフ化して確認する方法と検定で行う方法がある。
グラフで確認する場合は、基準値の前後がなめらかになっているかどうかを確認する。基準値の後に大きな変化があるようであれば、回帰不連続デザインの仮定が満たされない可能性がある。
set.seed(1989)
x <- rnorm(500, 0, 5)
hist(x, breaks = 30)
abline(v = 0)
検定を行う際は、Cattaneo, Jansson and Maが提案する方法で検定できる。この検定では、帰無仮説が右連続と左連続が一致する、つまり連続であるが設定され、対立仮説は連続ではない、である。つまり、連続であるためには、帰無仮説が棄却されると困るので、p値が大きければ帰無仮説を棄却できにくくなるので、そのほうが望ましい。
さきほどのxで検定すると、p=0.2925なので、帰無仮説が棄却できず、帰無仮説を採択することになる。つまり、連続性の仮定を満たしていることになる。
library(rddensity)
res <- rddensity(x, c = 0)
summary(res)##
## Manipulation testing using local polynomial density estimation.
##
## Number of obs = 500
## Model = unrestricted
## Kernel = triangular
## BW method = estimated
## VCE method = jackknife
##
## c = 0 Left of c Right of c
## Number of obs 262 238
## Eff. Number of obs 142 88
## Order est. (p) 2 2
## Order bias (q) 3 3
## BW est. (h) 3.734 2.629
##
## Method T P > |T|
## Robust -1.0526 0.2925
##
##
## P-values of binomial tests (H0: p=0.5).
##
## Window Length / 2 <c >=c P>|T|
## 0.226 12 8 0.5034
## 0.452 24 19 0.5424
## 0.678 35 22 0.1112
## 0.903 42 28 0.1196
## 1.129 49 33 0.0970
## 1.355 51 39 0.2461
## 1.581 57 45 0.2760
## 1.807 65 55 0.4114
## 2.033 76 61 0.2315
## 2.258 82 72 0.4684検定結果を可視化する際は、rdplotdensity()を利用する。これをみると連続ではないようにも見える。P-values of binomial testsの結果は、どれも0.05を基準とすると下回っていないため、連続性の仮定が満たされていると思われる。
plt <- rdplotdensity(res, X = x, type = "both")