View slides in full screen
セットアップ
実習ではMeyersson（2014）の（一部の）再現を試みる。本研究は「イスラム系政党の市長が誕生することは、女性の教育参加（高校修了率など）に対して抑圧的に働くのか、あるいはエンパワーメントとして働くのか」という問いをRDDを使って分析した研究である。分析の結果、イスラム政党の市長が当選した自治体では、女性における高校修了率が大幅に上昇したことが明らかになった面白い研究である。分析は1994年の地方選挙を使用において、僅差で当選・落選したイスラム系首長に注目し、2000年の高校卒業者の割合を比較したものである。それでは今回使用するパッケージとデータセットを読み込んでおこう。
:: p_load (tidyverse,<- read_csv ("_data/Myersson_2014.csv" )
# A tibble: 3,204 × 18
province district is_district_center is_province_center is_sub_metro is_metro iwm ivs
<chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 Adana Aladag 1 0 0 0 -0.356 0.0154
2 Adana Aladag 0 0 0 0 -0.543 0.0251
3 Adana Buyuksehir 1 0 0 1 -0.204 0.111
4 Adana Ceyhan 0 0 0 0 -0.201 0.102
5 Adana Ceyhan 0 0 0 0 -0.430 0.0206
6 Adana Ceyhan 0 0 0 0 -0.450 0.0559
7 Adana Ceyhan 0 0 0 0 -0.519 0.00249
8 Adana Ceyhan 1 0 0 0 -0.207 0.0806
9 Adana Ceyhan 0 0 0 0 -0.514 0.00656
10 Adana Ceyhan 0 0 0 0 -0.578 0.00706
# ℹ 3,194 more rows
# ℹ 10 more variables: n_parties <dbl>, log_pop <dbl>, female_hs_2000 <dbl>, male_hs_2000 <dbl>,
# female_hs_1990 <dbl>, male_hs_1990 <dbl>, hh_size <dbl>, pop_young <dbl>, pop_old <dbl>,
# gender_ratio <dbl>
今回使用するデータはMeyersson（2014）が公開した再現用データの一部を抜粋したものである。本来は236変数で構成されており、本研究の本文・付録に掲載されているすべての分析を再現しようとすれば、これらのほとんどが必要だが、ここではMeyersson（2014）の主要な分析結果のみを再現するので、必要な18変数のみ利用する。
province州/県（İl）
district地区/群（ilçe）
is_province_center
is_district_center
is_metro大都市ダミー
is_submetro副大都市ダミー
iwmイスラム系候補者の惜敗率
1994年地方選挙
ivsイスラム系候補者の得票率
1994年地方選挙
n_parties得票した政党の数
1994年地方選挙
log_pop人口（対数）
2000年
female_hs_200015-20歳女性の高卒者の割合
2000年
male_hs_200015-20歳男性の高卒者の割合
2000年
female_hs_199015-20歳女性の高卒者の割合
1990年
male_hs_199015-20歳男性の高卒者の割合
1990年
hh_size平均世帯人数
2000年
pop_young19歳以下人口の割合
2000年
pop_old60歳超え人口の割合
2000年
gender_ratio男女比
2000年
分析に入る前に簡単なデータハンドリングをしておく。まず、処置変数treatを作成する。実際の推定の際は不要な変数だが図表を作成する際に使うこともあるので、割当変数（running variable; RV）が閾値を超えたら1、それ以外は0をとる変数を作る。ここでは1994年地方選挙におけるイスラム系候補者の惜敗率（iwm）が割当変数だ。これが0.1なら、イスラム系政党候補者が2位の候補者に比べ10%ポイント得票率が高かった、つまり、当選したことを意味する（2位より得票率が高いことは1位ということになる）。一方、-0.05なら1位の候補者に比べ、得票率が5%ポイント低く、落選したことになる。したがって、iwmが0より大きい場合は1、それ以外は0の値をとるtreat変数をif_else()関数で作成する。
<- df |> mutate (treat = if_else (iwm > 0 , 1 , 0 ))
続いて、メイン分析の推定に使用する変数が1つでも欠損している行を除外する。たとえば、df |> drop_na(x)は変数xの値が欠損（NA）している場合、その行 を丸ごと落とすことを意味する。もし、drop_na(x)の変わりにdrop_na(x, y)にすれば変数xか y、1つでも欠損していればその行を落とすことになる。もったいないと思うかも知れないが、どうせ分析の時に欠損値を含む行は自動的に落とされる。むしろ、共変量あり/なしバージョンを推定する場合は、事前に欠損値の整理をしておかないと大変なことになるので、事前に落としておく必要がある。
たとえば、同じデータセット（ n = 100 ）を使って、y ~ dとy ~ d + xのモデルを推定するとする、xに10個の欠損値が含まれているとしよう。もし欠損値を含む行を整理せずに推定を続けると、モデル1のサンプルサイズは100、モデル2のサンプルサイズは90になる。ここで2つの推定結果（たとえば、dの係数）が異なるとすれば、これは共変量調整が原因と言えるだろうか。もしかしたら、脱落した10行が原因だったかも知れない。2つのモデルは使用しているデータセットがことなるため、モデル間比較ができない。xが欠損している10行分を予め削除しておくと、2つのモデルはいずれも n = 90だから、推定結果の違いが見られたらそれは共変量調整によるものだと判断することができる。
<- df |> drop_na (female_hs_2000, male_hs_2000, treat, ivs, n_parties,
# A tibble: 2,629 × 19
province district is_district_center is_province_center is_sub_metro is_metro iwm ivs
<chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 Adana Aladag 1 0 0 0 -0.356 0.0154
2 Adana Aladag 0 0 0 0 -0.543 0.0251
3 Adana Buyuksehir 1 0 0 1 -0.204 0.111
4 Adana Ceyhan 0 0 0 0 -0.201 0.102
5 Adana Ceyhan 0 0 0 0 -0.430 0.0206
6 Adana Ceyhan 0 0 0 0 -0.450 0.0559
7 Adana Ceyhan 0 0 0 0 -0.519 0.00249
8 Adana Ceyhan 1 0 0 0 -0.207 0.0806
9 Adana Ceyhan 0 0 0 0 -0.514 0.00656
10 Adana Ceyhan 0 0 0 0 -0.578 0.00706
# ℹ 2,619 more rows
# ℹ 11 more variables: n_parties <dbl>, log_pop <dbl>, female_hs_2000 <dbl>, male_hs_2000 <dbl>,
# female_hs_1990 <dbl>, male_hs_1990 <dbl>, hh_size <dbl>, pop_young <dbl>, pop_old <dbl>,
# gender_ratio <dbl>, treat <dbl>
約600行分が削除され、残り2629行となっている。
記述統計量
まず、記述統計量を確認してみよう。Meyersson（2014）のTable 1では全体の平均値と標準偏差、グループ（treatの値が0か1か）ごとに平均値と標準誤差を掲載している。select()関数でdfからTable 1に登場した共変量のみを抽出し、{modelsummary}のdatasummary_skim()で全体（pooled）の記述統計量を確認する。小数点3桁まで出力するためにfmt = fmt_decimal(digits = 3)を指定知る。
|> select (female_hs_2000, male_hs_2000, treat, ivs, n_parties, |> datasummary_skim (fmt = fmt_decimal (digits = 3 ))
Unique
Missing Pct.
Mean
SD
Min
Median
Max
Histogram
female_hs_2000
2291
0
0.163
0.096
0.000
0.155
0.680
male_hs_2000
2271
0
0.192
0.077
0.000
0.187
0.683
treat
2
0
0.120
0.325
0.000
0.000
1.000
ivs
2492
0
0.139
0.154
0.000
0.070
0.995
n_parties
14
0
5.541
2.192
1.000
5.000
14.000
log_pop
1991
0
7.840
1.188
5.493
7.479
15.338
pop_young
2628
0
0.405
0.083
0.065
0.397
0.688
pop_old
2625
0
0.092
0.040
0.017
0.085
0.272
gender_ratio
2611
0
1.073
0.253
0.750
1.032
10.336
hh_size
2627
0
5.835
2.360
2.823
5.274
33.634
is_district_center
2
0
0.345
0.475
0.000
0.000
1.000
is_province_center
2
0
0.023
0.149
0.000
0.000
1.000
is_sub_metro
2
0
0.022
0.146
0.000
0.000
1.000
次は処置群（treat = 1）と統制群（treat = 0）の記述統計量だ。ここでは{modelsummary}のdatasummary_balance()を使う。datasumamry_skim()とは違って、第1引数がformula（~ 処置変数名）になっている。データフレーム名はdata引数に当てる。パイプ演算子を使う場合はプレイスホルダー（_）を使う。
|> select (female_hs_2000, male_hs_2000, treat, ivs, n_parties, |> datasummary_balance (~ treat, data = _, fmt = fmt_decimal (digits = 3 ))
0
1
Mean
Std. Dev.
Mean
Std. Dev.
Diff. in Means
Std. Error
female_hs_2000
0.166
0.096
0.140
0.090
-0.026
0.005
male_hs_2000
0.192
0.078
0.196
0.076
0.004
0.005
ivs
0.101
0.116
0.415
0.114
0.313
0.007
n_parties
5.494
2.050
5.889
3.019
0.395
0.175
log_pop
7.775
1.070
8.315
1.767
0.540
0.102
pop_young
0.400
0.082
0.445
0.075
0.046
0.005
pop_old
0.095
0.040
0.073
0.031
-0.022
0.002
gender_ratio
1.073
0.266
1.076
0.117
0.003
0.009
hh_size
5.752
2.376
6.445
2.147
0.693
0.131
is_district_center
0.338
0.473
0.394
0.489
0.056
0.029
is_province_center
0.017
0.129
0.067
0.250
0.050
0.014
is_sub_metro
0.014
0.119
0.076
0.266
0.062
0.015
記述統計を見るとイスラム系首長のいる自治体の場合、高校を卒業した15-20歳女性の割合が約14.0%、非イスラム系首長のいる自治体のそれは約16.6である。この差分は約-2.6%ポイントだ。男性のそれが0.4%ポイントであることを考えると、やはり通念通り、イスラム政党の政権掌握は女性の権利を後退させるように見える。
しかし、これは単純比較から得られた結果だ。観察データには常に内生性が残っている。単に、イスラム信者の多い地域においてイスラム系首長が誕生しやすく、女性に教育を受けさせない家庭が多い可能性もある。もしそうであれば、これは家庭の問題であって首長が原因とは言えないだろう。やはり、より厳密なリサーチデザインが必要だろう。ここで登場することが回帰不連続デザイン（regression discontinuity design; RDD）だ。イスラム系首長がギリギリ当選した自治体とギリギリ落選した自治体であれば、自治体のしての特徴はそこまで変わらないだろう。Meyersson（2014）はこの点に着目し、イスラム政党による政治的統制が女性のエンパワーメントに影響を与えるかどうかを検証した研究成果である。
推定
それではLATEを推定してみよう。Meyersson（2014）のFigure 4とTable 2が主要推定結果である。RDDを使った研究は図と表の組み合わせが定番である。図で大体の・・・ 雰囲気を見せ、表で具体的な処置効果を示す。Figure 4は4つの図で、Table 2は計18のモデルで構成されているが、ここでは紙幅（…?）の関係上、すべての結果を再現せず、メインの結果のみを再現する。
LATEの推定
本稿のメインの推定結果はTable 2のPanel Aの左から4列目とPanel Bの左から4列目である。いずれも2000年における高卒の割合を結果変数だが、Panel Aが女性、Panel Bが男性である。
まず、女性の高卒割合を結果変数とし、共変量を入れていないモデル（Panel Aの3列目）を推定してみよう。使用する関数は{rdrobust}のrdrobust()だ。必須引数はy、x、c、dataだ。yとxは結果変数と割当変数名、dataはxとyで指定した変数が格納されているデータフレーム名だ。最後にcは処置を受ける閾値（c utpoint）だ。今回はiwmが0より大きくなると処置群に割り当てられるからc = 0になる。ただし、cは既定値が0であり、ここの実習では指定しなくても問題ない（一応、最初は付けてみよう）。
rdrobust (y = female_hs_2000, x = iwm, c = 0 , data = df)
Call: rdrobust
Sharp RD estimates using local polynomial regression.
Number of Obs. 2629
BW type mserd
Kernel Triangular
VCE method NN
Left Right
Number of Obs. 2314 315
Eff. Number of Obs. 529 266
Order est. (p) 1 1
Order bias (q) 2 2
BW est. (h) 0.172 0.172
BW bias (b) 0.286 0.286
rho (h/b) 0.603 0.603
Unique Obs. 2312 315
何か表示されるが、推定結果はどこにあるだろうか。推定結果を見るためにはrdrobust()で推定したオブジェクトをsummary()関数に通す必要がある。
rdrobust (y = female_hs_2000, x = iwm, data = df) |> summary (all = TRUE )
Call: rdrobust
Sharp RD estimates using local polynomial regression.
Number of Obs. 2629
BW type mserd
Kernel Triangular
VCE method NN
Left Right
Number of Obs. 2314 315
Eff. Number of Obs. 529 266
Order est. (p) 1 1
Order bias (q) 2 2
BW est. (h) 0.172 0.172
BW bias (b) 0.286 0.286
rho (h/b) 0.603 0.603
Unique Obs. 2312 315
=======================================================================================
Method Coef. Std. Err. z P>|z| [ 95% C.I. ]
=======================================================================================
Conventional 0.030 0.014 2.116 0.034 [0.002 , 0.058]
Bias-Corrected 0.030 0.014 2.090 0.037 [0.002 , 0.058]
Robust 0.030 0.017 1.776 0.076 [-0.003 , 0.063]
=======================================================================================
{broom}のtidy()関数を使えばLATEの推定結果をデータフレーム型で抽出できる。ただし、この表をそのまま使うのは危険である。点推定値はConventional行を、それ以外はRobust行を使うことが推奨されている（後述）。
rdrobust (y = female_hs_2000, x = iwm, data = df) |> tidy ()
term estimate std.error statistic p.value conf.low conf.high
1 Conventional 0.03019525 0.01427059 2.115907 0.03435268 0.002225406 0.05816509
2 Bias-Corrected 0.02983235 0.01427059 2.090477 0.03657494 0.001862504 0.05780219
3 Robust 0.02983235 0.01679899 1.775842 0.07575901 -0.003093068 0.06275776
計3つの推定結果（Conventional、Bias-Corrected、Robust）が表示され、結果も微妙に異なる。どの結果が、Meyersson（2014）の結果と一致するだろうか。論文を見直してみても良いが、結論からいうと一致するものは「ない」。Meyersson（2014）の推定結果は点推定値が0.032、標準誤差が0.010になっている。この相違を説明するのはかなり面倒だ。簡単にいえば、Meyersson（2014）のRDDと{rdrobust}がデフォルトで提供しているRDDは同じRDDとはいえ、中身の設定が相当異なると言える。まず、使用しているバンド幅が異なる。Meyersson（2014）はImbens and Kalyanaraman（2012）、通常IK、またはIK12バンド幅を使用しているが、{rdrobust}はCalonico et al. （2014）のバンド幅を使用している。数年前までは引数でIK12が指定できたが、今はできなくなっている。
それではバンド幅を変えれば同じ結果が得られるだろうか。IK12バンド幅の値は0.2397612だ。これを使ってみよう。rdrobust()でバンド幅を手動で指定したい時にはh引数を使用する。
rdrobust (y = female_hs_2000, x = iwm, data = df,h = 0.2397612 ) |> summary (all = TRUE )
Call: rdrobust
Sharp RD estimates using local polynomial regression.
Number of Obs. 2629
BW type Manual
Kernel Triangular
VCE method NN
Left Right
Number of Obs. 2314 315
Eff. Number of Obs. 727 293
Order est. (p) 1 1
Order bias (q) 2 2
BW est. (h) 0.240 0.240
BW bias (b) 0.240 0.240
rho (h/b) 1.000 1.000
Unique Obs. 2312 315
=======================================================================================
Method Coef. Std. Err. z P>|z| [ 95% C.I. ]
=======================================================================================
Conventional 0.030 0.013 2.362 0.018 [0.005 , 0.054]
Bias-Corrected 0.027 0.013 2.187 0.029 [0.003 , 0.052]
Robust 0.027 0.018 1.541 0.123 [-0.007 , 0.062]
=======================================================================================
それでも結果が全然合わない。元の論文と{rdrobust}の違いはバンド幅の選び方だけではない。カーネル関数も異なる。講義ではRDDにおけるカーネル関数の事実上の標準は三角形カーネル関数（traiangular）と述べた（はず）が、{rdrobust}でもこれを採用している。しかし、Meyersson（2014）はなぜか一様カーネル関数（uniform / rectangular）を使用している。それではkernel引数でカーネル関数を一様カーネルに変えてみよう。
rdrobust (y = female_hs_2000, x = iwm, data = df,h = 0.2397612 , kernel = "uniform" ) |> summary (all = TRUE )
Call: rdrobust
Sharp RD estimates using local polynomial regression.
Number of Obs. 2629
BW type Manual
Kernel Uniform
VCE method NN
Left Right
Number of Obs. 2314 315
Eff. Number of Obs. 727 293
Order est. (p) 1 1
Order bias (q) 2 2
BW est. (h) 0.240 0.240
BW bias (b) 0.240 0.240
rho (h/b) 1.000 1.000
Unique Obs. 2312 315
=======================================================================================
Method Coef. Std. Err. z P>|z| [ 95% C.I. ]
=======================================================================================
Conventional 0.032 0.012 2.784 0.005 [0.009 , 0.055]
Bias-Corrected 0.026 0.012 2.242 0.025 [0.003 , 0.048]
Robust 0.026 0.017 1.557 0.120 [-0.007 , 0.058]
=======================================================================================
Conventional行で点推定値（0.032）が一致した。点推定値のズレの謎はこれで解決したが、標準誤差がややことなる。論文では0.010、{rdrobust}の推定結果は0.012だ。結論からいうと、これはどうしようもない。{rdrobust}ではAbadie and Imbens（2006）の不均一分散に頑健な最近傍分散推定量（heteroskedasticity-robust nearest neighbor variance estimator）が採用されており、Meyersson（2014）は分散不均一性を考慮した標準誤差を採用していない。問題は{rdrobust}では頑健ではない通常の標準誤差は使えないということだ。つまり、{rdrobust}で本論文の結果を100%再現することはできない。誤解しないで頂きたいが、これはMeyersson（2014）がわざと頑健ではない手法を採用したとか、方法論的に間違った手法を使ったとか、あるいは時代遅れだったとか、そういうものではない。この論文が執筆された時期を考えるとある意味、これが標準的な作法だったともいえる。実際に論文を読んで見れば分かるが、「ここまでする?」と思わせるレベルで方法論的に真摯に向かってる良い論文だ。
話を戻し、Conventional、Bias-corrected、Robustのどの推定結果を使うべきかについて考えてみよう。なるべく論文の正確に再現したいのであればConventionalだが、自分の論文であればより頑健なRobustを使うのが良いだろう。{rdrobust}パッケージでの推定方法等に関してはCalonico et al. （2014）を参照すること。
それでは共変量を投入した推定方法を紹介する。ここからは{rdrobust}のデフォルト設定（CCTバンド幅、頑健標準誤差、三角形カーネル関数）を使用し、推定結果はRobust列を使う。rdrobust()で共変量はcovs引数で指定する。具体的にはcovs ~ 共変量1 + 共変量2 + ...のように書く。使用する共変量はTable 2の下段に書かれている。また、論文内では推定後、県（province）単位でクラスタリングした標準誤差を使用していることを明記しているので、cluster = provinceを追加し、province変数でクラスタリングした頑健標準誤差を使用する。クラスター標準誤差の計算方法としてCR1、CR2、CR3が提供されており、cluster引数を指定すると自動的にCR1が適用される。CR2やCR3に変更したい場合はvce = "cr2"かvce = "cr3"を引数として追加すれば良い。CR3はやや大きめの標準誤差を出す保守的な計算方法だが、一般的に採用されている方法はCR2である（Imbens and Kolesár 2016）。ここでもCR2を使うとする。推定したオブジェクトはそのままsummary()関数に渡し、LATE推定結果を確認してみよう。
rdrobust (y = female_hs_2000, x = iwm, data = df, covs = ~ ivs + n_parties + pop_young + pop_old + gender_ratio + log_pop + + is_district_center + is_sub_metro + is_metro + province, cluster = province, vce = "CR2" ) |> summary (all = TRUE )
Call: rdrobust
Covariate-adjusted Sharp RD estimates using local polynomial regression. Std. errors are clustered (120 clusters).
Number of Obs. 2629
BW type mserd
Kernel Triangular
VCE method CR2
Left Right
Number of Obs. 2314 315
Eff. Number of Obs. 282 187
Order est. (p) 1 1
Order bias (q) 2 2
BW est. (h) 0.096 0.096
BW bias (b) 0.208 0.208
rho (h/b) 0.461 0.461
Clusters (g) 67 53
Unique Obs. 2312 315
=======================================================================================
Method Coef. Std. Err. z P>|z| [ 95% C.I. ]
=======================================================================================
Conventional 0.025 0.008 3.096 0.002 [0.009 , 0.041]
Bias-Corrected 0.026 0.008 3.231 0.001 [0.010 , 0.042]
Robust 0.026 0.009 2.938 0.003 [0.009 , 0.043]
=======================================================================================
これから共変量を投入したモデルを推定していくわけだが、毎回共変量のを書くのは面倒だろう。formula()関数を使えば、これらの共変量（~ ivs + n_parties + ...）を一つのオブジェクトとして扱うことができる。共変量リストをcov_formulaと名付けて格納してみよう。
<- formula (~ ivs + n_parties + pop_young + pop_old + gender_ratio + log_pop + + is_province_center + is_sub_metro + + province)
これからはcovs = cov_formulaと書くだけで、上記の11個の共変量リスト作成の手間が省ける。
rdrobust (y = female_hs_2000, x = iwm, data = df, covs = cov_formula, cluster = province, vce = "CR2" ) |> summary (all = TRUE )
Call: rdrobust
Covariate-adjusted Sharp RD estimates using local polynomial regression. Std. errors are clustered (120 clusters).
Number of Obs. 2629
BW type mserd
Kernel Triangular
VCE method CR2
Left Right
Number of Obs. 2314 315
Eff. Number of Obs. 282 186
Order est. (p) 1 1
Order bias (q) 2 2
BW est. (h) 0.096 0.096
BW bias (b) 0.208 0.208
rho (h/b) 0.461 0.461
Clusters (g) 67 53
Unique Obs. 2312 315
=======================================================================================
Method Coef. Std. Err. z P>|z| [ 95% C.I. ]
=======================================================================================
Conventional 0.025 0.008 3.095 0.002 [0.009 , 0.041]
Bias-Corrected 0.026 0.008 3.229 0.001 [0.010 , 0.042]
Robust 0.026 0.009 2.937 0.003 [0.009 , 0.043]
=======================================================================================
それでは本稿のメインのモデルを推定してみよう。結果変数の2種類であり、（1）2000年におけるの15〜20歳女性の高校卒業者割合（female_hs_2000）と（2）2000年におけるの15〜20歳男性の高校卒業者割合（male_hs_2000）だ。推定結果はそれぞれrdd_female_2000、rdd_male_2000と名付ける。
<- rdrobust (y = female_hs_2000, x = iwm, data = df, covs = cov_formula, cluster = province, vce = "CR2" )<- rdrobust (y = male_hs_2000, x = iwm, data = df, covs = cov_formula, cluster = province, vce = "CR2" )
これらの結果を{broom}のtidy()関数でLATEの推定値のみ表形式で抽出できるが、一つ注意点が必要だ。たとえば、rdd_female_2000の推定結果をtidy()で抽出してみよう
term estimate std.error statistic p.value conf.low conf.high
1 Conventional 0.02482684 0.008021680 3.094967 0.001968345 0.009104635 0.04054904
2 Bias-Corrected 0.02590521 0.008021680 3.229400 0.001240502 0.010183010 0.04162742
3 Robust 0.02590521 0.008820529 2.936923 0.003314864 0.008617295 0.04319313
summary()関数では点推定値が約0.025で、p 値が約0.003になっているが、これと一致する結果がない。これには理由がある。LATEの点推定値はConventional行を、統計的仮説検定に関わる推定量（標準誤差、\(z\) 統計量、p 値、信頼区間）はRobust行を使うからだ（Cattaneo et al. 2020; Valentim et al. 2021）。tidy()から得られた結果をよく見ると、たしかに2つの結果が混在している。
今のtidy()関数だけでは私たちに欲しい表が出来上がらないので、Robust行だけ残し、estimateの列だけはConventional行のものを取ってくるrdd_tidy()という新しい関数を追加してみよう。関数の作り方については『私たちのR』の「Rプログラミングの基礎 」と「関数の自作 」を参照するとし、ここでは関数の完成したコードだけ示す。
# 関数名は rdd_tidy で、必要な引数は x のみ <- function (x) {# tidy(x) の結果を temp と名付ける <- tidy (x)# temp の3行目だけを抽出し result と名付ける # また、1列目も不要だから[3, -1]とする <- temp[3 , - 1 ]# reuslt の estimate 列に temp の estimate 列の1番目の値を入れる $ estimate <- temp$ estimate[1 ]# result を返す
それではやってみよう。
rdd_tidy (rdd_female_2000)
estimate std.error statistic p.value conf.low conf.high
3 0.02482684 0.008820529 2.936923 0.003314864 0.008617295 0.04319313
Robust行をベースにConventionalの点推定値で上書きされた推定結果が出力される。これをbind_rows()関数で一つのデータフレームにまとめてみよう。bind_rows()関数の使い方は「Rの復習 」を参照すること。
list ("Female" = rdd_tidy (rdd_female_2000),"Male" = rdd_tidy (rdd_male_2000)) |> bind_rows (.id = "gender" )
gender estimate std.error statistic p.value conf.low conf.high
1 Female 0.02482684 0.008820529 2.936923 0.003314864 0.008617295 0.04319313
2 Male 0.01171863 0.009998157 1.177157 0.239132872 -0.007826628 0.03136543
以下はこの表からRobustの行のみ残し、pointrangeプロットで可視化するコードだ。
list ("Female" = rdd_tidy (rdd_female_2000),"Male" = rdd_tidy (rdd_male_2000)) |> bind_rows (.id = "gender" ) |> ggplot () + geom_hline (yintercept = 0 , linetype = "dashed" ) + geom_pointrange (aes (x = gender, y = estimate, ymin = conf.low, ymax = conf.high)) + labs (x = "Gender" , y = "LATE with 95% CIs" ) + theme_classic ()
イスラム系首長の誕生は女性の世俗的な高校教育を統計的有意に増加させたことが分かった。一方、男性のそれは効果量が小さく、統計的有意な処置効果は確認できなかった。
可視化
続いて、Figure 4を再現してみよう。Figure 4は計4つの図で構成されされているが、1つでも作れるようになったら、あとは使用する変数と図のスケールを少しいじるだけが。まずはメインの分析結果である2000年におけるの15〜20歳女性の高校卒業者割合（female_hs_2000）を応答変数とした図（Figure 4の左上）を作ってみよう。
まずは割当変数（iwm）と結果変数（female_hs_2000）の散布図を作ってみよう。論文のFigure 4ではiwmの値が-0.5以上、0.5以下のケースのみを使用しているため、dfをggplot()に渡す前にfilter()関数でiwmの値が-0.5以上、0.5以下のケースのみに絞る。
|> filter (iwm >= - 0.5 & iwm <= 0.5 ) |> ggplot (aes (x = iwm, y = female_hs_2000)) + geom_point ()
続いて、曲線をオーバラップしてみよう。この曲線は平滑化曲線 （smoothing line）と呼ばれるものであり、2つの変数間の関係を示す滑らかな曲線（または直線）である。よく使う回帰直線もできるが、RDD界隈では一般的にLOESS（locally estimated scatterplot smoothing）と呼ばれる局所回帰で推定した線が用いられる。RDDのLATEを推定する際の局所回帰に似たようなものであるが、簡単にいえば、「データ全体に1本のカーブを当てはめるのではなく、各点の周りだけを見て、その近所のデータだけで簡単な回帰線を引く」を、全部の点で繰り返して、それらをつなげたものである。そうすることで直線ではなく、滑らかな曲線が得られる。LOESSの詳細はCleveland（1979）をネット上の解説記事等を参照されたい。
LOESS線はgeom_smooth()レイヤーで足せる。method = "loess"をすればLOESSが、"lm"だと線形回帰、"gam"だと一般化加法モデルが使われる（自分で作った関数を使用することもできる）。また、マッピングではxとyだけでなく、groupも必要だ。groupを指定しないと一本の線が表示されるが、group = treatにしておくと、統制群と処置群、それぞれのLOESS曲線が計算され、出力される。
|> filter (iwm >= - 0.5 & iwm <= 0.5 ) |> ggplot () + geom_point (aes (x = iwm, y = female_hs_2000)) + geom_smooth (aes (x = iwm, y = female_hs_2000, group = treat), method = "loess" )
コードを見ると、geom_point()とgeom_smooth()はxとyのマッピングを共有しているため、この2つはggplot()内でマッピングしておくとより簡潔なコードになる。
|> filter (iwm >= - 0.5 & iwm <= 0.5 ) |> ggplot (aes (x = iwm, y = female_hs_2000)) + geom_point () + geom_smooth (aes (group = treat), method = "loess" )
線がジャンプしていることが分かる。実はこれだけで目的は達成したが、Figure 4を見ると散布図が非常にシンプルになっている。RDDにおける可視化は情報量よりも傾向 が重要だ。ここまで点が多いと、線はジャンプしていても、点の群れがジャンプしているかどうかは分かりにくい。ここでは以下の図のようにデータを4%ポイントの区間に分け、その区間内の結果変数（female_hs_2000）の平均値をプロットしてみよう。
任意の区間（bin）内で変数を要約（summary）し、プロットする関数はstat_summary_bin()だ。ここで要約方法は平均値（fun = "mean"）、出力方法は散布図（geom = "point"）とする。区間の幅は0.04（bindwidth = 0.04）にする。ここでさらにboundary = 0を指定しておくと、区間が0にかからなくなり、Figure 4のような図が作れるようになる。後はgeom_point()同様、size、shape、color、fill等で点のカスタマイズをする。ここでは点の大きさは3、形は枠線付きの円（21番）、枠線の色は黒、中身の色塗りはグレーとする。
|> filter (iwm >= - 0.5 & iwm <= 0.5 ) |> ggplot (aes (x = iwm, y = female_hs_2000, group = treat)) + stat_summary_bin (fun = "mean" , geom = "point" , binwidth = 0.04 , boundary = 0 ,size = 3 , shape = 21 , color = "black" , fill = "gray" )
ここでさらにLOESS曲線と割当変数の閾値に垂直線（geom_vline()）を乗せる。線の色は黒（color = "black"）、太さはやや細めの0.5（linewidth = 0.5）とする。
|> filter (iwm >= - 0.5 & iwm <= 0.5 ) |> ggplot (aes (x = iwm, y = female_hs_2000)) + geom_vline (xintercept = 0 , linetype = "dashed" , color = "gray" ) + stat_summary_bin (fun = "mean" , geom = "point" , binwidth = 0.04 , boundary = 0 ,size = 3 , shape = 21 , color = "black" , fill = "gray" ) + geom_smooth (aes (group = treat), method = "loess" , color = "black" , linewidth = 0.5 )
geom_smooth()内のspan引数で曲線の滑らかさは調整できる。既定値は0.75だが、これを大きくすればするほど滑らかな曲線となり、小さいするとより点（データポイント）にフィットするギザギザな線になる。
Error in `theme_jay()`:
! could not find function "theme_jay"
図 1: span引数による曲線の滑らかさの調節
あとはlabs()で軸タイトル、図のタイトルを修正し、coord_cartesian()で座標軸の下限と上限を調節する。今回は縦軸の下限と上限をそれぞれ-0.05、0.25にしたいからylim = c(-0.05, 0.25)とする。テーマはMeyersson（2014）内の図と似た雰囲気のクラシック（theme_classic()）を使う。
|> filter (iwm >= - 0.5 & iwm <= 0.5 ) |> ggplot (aes (x = iwm, y = female_hs_2000)) + geom_vline (xintercept = 0 , linetype = "dashed" , color = "gray" ) + stat_summary_bin (fun = "mean" , geom = "point" , binwidth = 0.04 , boundary = 0 ,size = 3 , shape = 21 , color = "black" , fill = "gray" ) + geom_smooth (aes (group = treat), method = "loess" , color = "black" , linewidth = 0.5 ) + labs (x = "Islamic win margin in 1994" , y = "Share of 15-20-year-olds \n with high school degree" ,title = "(a) Women in 2000" ) + coord_cartesian (ylim = c (- 0.05 , 0.25 )) + theme_classic ()
同じ要領で縦軸（結果変数）を2000年におけるの15〜20歳男性の高校卒業者割合（male_hs_2000）、1990年におけるの15〜20歳女性の高校卒業者割合（female_hs_1990）、1990年におけるの15〜20歳男性の高校卒業者割合（male_hs_1990）にした図を作成してみよう（以下の「▶ コード」をクリックすればコードが表示される）。
コード
|> filter (iwm >= - 0.5 & iwm <= 0.5 ) |> ggplot (aes (x = iwm, y = female_hs_2000)) + geom_vline (xintercept = 0 , linetype = "dashed" , color = "gray" ) + stat_summary_bin (fun = "mean" , geom = "point" , binwidth = 0.04 , boundary = 0 ,size = 3 , shape = 21 , color = "black" , fill = "gray" ) + geom_smooth (aes (group = treat), method = "loess" , color = "black" , linewidth = 0.5 ) + labs (x = "Islamic win margin in 1994" , y = "Share of 15-20-year-olds \n with high school degree" , title = "(a) Women in 2000" ) + coord_cartesian (ylim = c (- 0.05 , 0.25 )) + theme_classic ()|> filter (iwm >= - 0.5 & iwm <= 0.5 ) |> ggplot (aes (x = iwm, y = male_hs_2000)) + geom_vline (xintercept = 0 , linetype = "dashed" , color = "gray" ) + stat_summary_bin (fun = "mean" , geom = "point" , binwidth = 0.04 , boundary = 0 ,size = 3 , shape = 21 , color = "black" , fill = "gray" ) + geom_smooth (aes (group = treat), method = "loess" , color = "black" , linewidth = 0.5 ) + labs (x = "Islamic win margin in 1994" , y = "Share of 15-20-year-olds \n with high school degree" , title = "(b) Men in 2000" ) + coord_cartesian (ylim = c (- 0.05 , 0.25 )) + theme_classic ()|> filter (iwm >= - 0.5 & iwm <= 0.5 ) |> ggplot (aes (x = iwm, y = female_hs_1990)) + geom_vline (xintercept = 0 , linetype = "dashed" , color = "gray" ) + stat_summary_bin (fun = "mean" , geom = "point" , binwidth = 0.04 , boundary = 0 ,size = 3 , shape = 21 , color = "black" , fill = "gray" ) + geom_smooth (aes (group = treat), method = "loess" , color = "black" , linewidth = 0.5 ) + labs (x = "Islamic win margin in 1994" , y = "Share of 15-20-year-olds \n with high school degree" , title = "(c) Women in 1990" ) + coord_cartesian (ylim = c (- 0.05 , 0.25 )) + theme_classic ()|> filter (iwm >= - 0.5 & iwm <= 0.5 ) |> ggplot (aes (x = iwm, y = male_hs_1990)) + geom_vline (xintercept = 0 , linetype = "dashed" , color = "gray" ) + stat_summary_bin (fun = "mean" , geom = "point" , binwidth = 0.04 , boundary = 0 ,size = 3 , shape = 21 , color = "black" , fill = "gray" ) + geom_smooth (aes (group = treat), method = "loess" , color = "black" , linewidth = 0.5 ) + labs (x = "Islamic win margin in 1994" , y = "Share of 15-20-year-olds \n with high school degree" , title = "(d) Men in 1990" ) + coord_cartesian (ylim = c (- 0.05 , 0.25 )) + theme_classic ()
図 2: Meyersson（2014）のFigure 4の再現
イスラム系首長の誕生は女性の世俗的な高校教育を増加させるという傾向が確認できた。また、1990年における高校卒業者の割合を縦軸にした分析では線の大きなジャンプが確認出来ていない。1994年の選挙結果前にジャンプが見られないことは、逆の因果や第3の要因（交絡要因）がないことを強く示唆している。
しかし、こういった図からでは統計的有意性検定はできず、あくまでも「傾向」を確認することが主な目的である。曲線のジャンプの幅や信頼区間はLOESSの設定によっていくらでも変わるものであり、この図から効果の有無を語ることはできない。処置効果の有無やその効果量の大きさを語る際には図ではなく、推定結果そのものの図表を使おう。
{rdrobust}を使ったrdplot()を使えばより簡単に上記の図が作れる。rdplot()のy、x、c、data引数はrdrobust()のそれと同じだ。ここにnbins引数を使うことで、閾値を中心とし左右に何個の点を出力するかが調節できる。
<- df |> filter (iwm >= - 0.5 & iwm < 0.5 ) |> rdplot (y = female_hs_2000, x = iwm, nbins = 15 , data = _)
出力された図は{ggplot2}で作られたものだ。図オブジェクト名$rdplotでオブジェクトの中から図だけを抽出できるため、ここから+演算子を使うとさらにレイヤーを追加できる。
$ rdplot + labs (x = "Islamic win margin in 1994" , y = "Share of 15-20-year-olds \n with high school degree" , title = "(a) Women in 2000" )
RDDの仮定
RDDの識別条件として重要なのは以下の3点です。
閾値周辺において潜在的結果が連続であること
閾値周辺において応答変数に影響を与える共変量（処置変数を除く）が連続すること
閾値周辺において割当変数が連続であること
これらの仮定はそれぞれプラセボ検定、バランスチェック、密度検定を使うことで仮定の違反有無が確認できる。ただ、仮定の違反 有無を確認することが限界であり、「仮定が満たされていない」ことは確認できても、「仮定が満たされている」は確認できない（統計的仮説検定の仕組みを考えてみよう）。せいぜい「「仮定が満たされていない」とは言えない」くらいまでしか確認できない。それでも仮定違反が見つかるだけでも十分に意義はあろう。
プラセボ検定
潜在的結果の連続性は「もし、閾値が今とは異なる別のポイントなら、今の閾値周辺において直線（または曲線）は連続しているだろう」ということである。飲酒年齢の例で考えると、交通事故率は20歳を境界にジャンプしていると考えられる。もし、法改正により飲酒可能な年齢が18歳に引き下げられた場合、今の20歳周辺で見られるジャンプは見られなくなる、これが潜在的結果が連続している状況である。しかし、我々は潜在的結果のうち、反事実は観察できない。したがって、この連続性を直接立証することはできない。
ここで発想を転換してみよう。「処置が発生しない場所でも、同じようなジャンプが生じていないか」を確認することはどうだろうか。もし、本当の閾値以外の場所でも、ジャンプが頻繁に見つかるなら、「このデータでは、処置がなくても割当変数に沿って不連続な変化が発生しているのでは?」という疑念が生じることになる。これがプラセボテスト（placebo test）の考え方だ。偽の閾値を与えた場合、統計的に有意なLATEが得られるかどうかを確認することである。プラセボテストの方法はいくつか考えられるが、ここではImbens and Lemieux（2008）が提案した方法を紹介する。
まず、filter()関数を使用し、データセットを処置群と統制群で分割する。統制群のデータセットはdf_c、処置群のそれはdf_tとする。
<- df |> filter (treat == 0 )<- df |> filter (treat == 1 )
つづいて、任意の閾値を決める。Imbens and Lemieux（2008）は閾値を中心に左側に個体の50%、右側に50%が入るよう閾値を提案している。簡単にいえば、中央値（= 第2四分位数）だ。median()関数を使えば簡単に中央値が計算できる。むろん、必ず中央値を使う必要はない。たとえば閾値を中心に左側に個体の40%、右側に60%が入るような閾値を決めても良い。その場合はqunatile()関数を使用すれば良い（quantile(データフレーム名$iwm, probs = 0.4)）。
あとはこの2つのデータセットと閾値を使用し、同じRDDを行えば良い。今回はrdd_tidy()関数でLATEの推定値だけ確認しよう。
rdrobust (y = female_hs_2000, x = iwm, data = df_c, c = median (df_c$ iwm)) |> rdd_tidy ()
estimate std.error statistic p.value conf.low conf.high
3 -0.008077541 0.01506885 -0.6359359 0.5248182 -0.03911724 0.01995159
rdrobust (y = female_hs_2000, x = iwm, data = df_t, c = median (df_t$ iwm)) |> rdd_tidy ()
estimate std.error statistic p.value conf.low conf.high
3 -0.0303105 0.04551953 -0.7126352 0.4760715 -0.1216555 0.05677781
いずれも統計的有意なジャンプが見られなかった。これが潜在的結果の連続性を保証するわけではないが、少なくとも潜在的結果の非連続性の証拠はまだ 見つかってないことを意味している。
以下のように通常とRDD同様、プラセボテストの図を添えても良いだろう。本文に載せるには字数制限がもったいないので、付録で良いかも知れない。
プラセボテストにはいくつかのやり方が考えられる。たとえばデータを分割せず、全体のデータセットに対し、割当変数の閾値を少しずつ動かしながらLATEを推定することも考えられる。rdrobust()のc引数だけを修正すれば良いし、簡単にできよう。
密度検定
割当変数の連続する密度検定の結果はFigure 2として掲載されている。割当変数の連続性が確認されない場合、割当変数に人為的な操作が加わっている可能性が示唆される。変数の人為的操作は日常的に観察できる。
以下は日本で販売されている乗用車（除く軽自動車）の車幅の分布を示したものである（AIで適当に集めたデータなので、そこまで信頼できるものではない）。日本で販売されている軽自動車以外の乗用車の車幅の平均値は約1,800mm、標準偏差は約85mmであり、これは赤い線で示している（もう一度強調するが、AIで適当に集めたデータである）。傾向としては分布とおりだが、異常な箇所がいくつかある。それは車幅が1,700mmと1,800mmにギリギリいっていない車が異常に多いことだ。車幅はメーカー側で決めるものであり、意図的にこうなっている。維持費が比較的安いと言われている通称「５ナンバー」車両の条件の一つが車幅1,700mm以下であるため、多くの車の車幅が1,695mmになっている。また、1,800mmにも理由がある。都心部のマンションにある機械式駐車場の多くが1,800mmまで入庫できるからだ。今は1,850mmの機械式駐車場が標準で、タワマンのように富裕層の多いマンションだと2,000mmまで入庫できる機械式駐車場も多いが、まだ多くが1,800mmが限界になっている。「 売れそうな」車幅の車種をたくさん作っていることが、このような車幅の分布密度の断絶の理由でもある（むろん、一番多いのは車幅は1,475mmであり、これは日本の軽自動車規格である1,480mmギリギリの数値である）。
もし割当変数にこのような人為的操作が出来るのであれば、その割当変数は慎重に扱わねばならない。むろん、人為的操作が可能でも操作が行われている可能性もあり、この場合は割当変数とて使っても良いだろう。その人為的操作可能性を調べるのが密度検定だ。その意味で今回の割当変数（iwm）は問題なさそうだ。iwmが0を超えると、それはイスラム系候補者が当選したことになるが、これが操作可能であれば、当選することが至上目標である政治家はみんなが操作してしまうだろう。むろん、自分の得票率を上げるために何かの不正をする可能性はあり、それによって得票率が上がるかも知れない。しかし、ギリギリな線で勝つかどうかは別である。これはどう考えても不可能だろう。つまり、イスラム系候補者が「0.1%ポイント差で勝つようにしよう!」というモチベーションで何かをすることは考えにくい（本当に僅差で勝利したらそっちの方が色々と大変だし）。
理論的には問題なさそうだが、政府や与党、その他の政治エリートによる選挙結果を操作している可能性も排除できないので、やはりデータで確認する必要がある。Figure 2を上段をみるとヒストグラムを使用しているが、同じものを作ってみよう。
|> ggplot () + geom_histogram (aes (x = iwm, fill = iwm > 0 ),color = "white" , boundary = 0 , binwidth = 0.02 ) + geom_vline (xintercept = 0 , linetype = "dashed" , color = "firebrick" ) + labs (x = "Isalmic party vote margin in 1994 (running variable)" ,y = "Number of municipalities" ,fill = "Mayor party" ) + scale_fill_manual (labels = c ("FALSE" = "Non-Islamic" , "TRUE" = "Islamic" ),values = c ("FALSE" = "#5F7387" , "TRUE" = "#E74C3C" )) + theme_classic () + theme (legend.position = "bottom" )
ギザギザしているが、不自然なものはないように見える。ただし、そう「見える」だけであって、統計的有意な不自然さがあるかも知れない。そこで必要なのが密度検定だ。この分布を適切に示す密度曲線を作成し、その密度曲線が閾値でジャンプしているかどうかを決めるものである。密度検定の詳細についてはMcCrary（2008）を参照して頂くとし、ここでは{rddapp}パッケージを利用した密度検定の実装方法を紹介する。
{rddapp}はdc_test()という密度検定用の関数を提供している。必須引数は2つだ。1つ目は割当変数のベクトル、2つ目は閾値（cutpoint）だ。割当変数は$演算子を使用し、dfからiwm列を抽出すれば良い。
dc_test (df$ iwm, cutpoint = 0 )
Binwidth:
0.008626
Bandwidth:
0.1649
Estimate for log difference in heights:
Estimate Std. Error lower.CL upper.CL z value Pr(>|z|)
-0.09298 0.14856 -0.38415 0.19819 -0.62587 0.53140
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Confidence interval used: 0.95
密度検定の結果で注目するのは密度曲線の「Pr(>|z|) 」列だ。ここが予め設定した有意水準（\(\alpha\) ）を下回れば、割当変数が閾値周辺にて非 連続だと判定する。今回の p 値は約0.531であり、よく使う基準である0.05より大きいことから割当変数が非連続である証拠は見つからなかったことになる。少なくとも仮定を違反しているとは言えないということだろう。
また、dc_test()を使うと自動的に図も返してくれるが、この図はカスタマイズが非常に難しい。{ggplot2}で作成した図ならばレイヤーを上書きする形で調整できるが、これはBase Rで作成されたものなので、すでに表示されているものを削除したり、修正することはできず、この図の上に何かを追加することだけである。もし、自分で作図をしたい場合は密度検定の歳、ext.out = TRUEを付けておくと、散布図に使用したデータフレームが抽出できる。さらに図が不要であればplot = FALSEも付けておこう。
<- dc_test (df$ iwm, cutpoint = 0 , plot = FALSE , ext.out = TRUE )
密度検定のオブジェクト$dataで散布図用のデータセットが抽出できる。
head (dens_test$ data) # 長いので最初の6行だけ確認
cellmp cellval
1 -0.9963213 0.08819053
2 -0.9876952 0.00000000
3 -0.9790690 0.00000000
4 -0.9704429 0.04409527
5 -0.9618167 0.00000000
6 -0.9531905 0.04409527
後はこのデータを利用し、散布図（geom_point()） + 平滑化曲線（geom_smooth()）を作るだけだ。
$ data |> mutate (treat = if_else (cellmp > 0 , 1 , 0 )) |> ggplot (aes (x = cellmp, y = cellval, group = treat)) + geom_point () + geom_smooth (method = "loess" , color = "black" , linewidth = 0.75 , span = 0.5 ) + labs (x = "Islamic win margin" , y = "# of districts" ) + coord_cartesian (xlim = c (- 0.5 , 0.5 )) + theme_classic ()
Figure 2の下段のものとはやや異なるが、傾向自体は同じである。ここでも図は傾向を示すことが主目的であり、非連続性の検討は密度検定の具体的な数値に基づいて行うべきであろう。
バランスチェック
バランスチェックはMeyersson（2014）のFigure 3として掲載されている。RDDでは割当変数の値が閾値に達するか、閾値を超えると線のジャンプが生じ、そのジャンプの幅がLATEとなる。しかし、同じタイミングで別の変数におてもジャンプが発生したら、LATEが識別できなくなる。処置によるものか、あるいはそのジャンプした別の変数から影響されたのかが分からなくなるからだ。しかし、バランスチェックはこれまでのランダム化比較試験（RCT）や傾向スコア推定のような標準化平均差（standardized mean difference; SMD）を使うことはできない。RDDにおけるバランスチックは共変量のジャンプが重要であって、その平均値は重要ではない。たとえば、 図 4 (a) の場合、SMDが約3.5である。緩い基準でも0.25を使っても非常に大きい数値となっている。しかし、図を見れば分かるが、線のジャンプはない。一方、 図 4 (b) の場合、SMDは0.05であり、数値だけ見ればバランスが取れているかのように見えるが、図を見れば明らかにジャンプしていることが分かる。このようにRDDにおけるバランスチェックはSMDでなく、可視化（+推定）を用いることが推奨される。
可視化はLATEの可視化と全く同じで、縦軸（結果変数）の変数だけ替えれば良い。LATEの推定でも述べた通り、RDDの図は傾向を確認するものであって、統計的有意なジャンプがあったかどうかは図では判定できない。たとえば、今回は適合曲線としてLOESSを採用したが、これはパラメーター（span）の値に応じて見た目が大きく変わる。したがって、論文には載っていなかったが、以下のようにLATEの点推定値および p 値も一緒に掲載しておくと、より説得力のあるバランスチェックになるだろう。
以下はFigure 3の最初の4つだけ再現した例だ（以下の「▶ コード」をクリックすればコードが表示される）。
コード
|> filter (iwm >= - 0.32 & iwm <= 0.32 ) |> ggplot (aes (x = iwm, y = ivs, group = treat)) + geom_vline (xintercept = 0 , linetype = "dashed" , color = "gray" ) + stat_summary_bin (fun = "mean" , geom = "point" , binwidth = 0.08 , boundary = 0 ,size = 3 , shape = 21 , color = "black" , fill = "gray" ) + geom_smooth (method = "loess" , color = "black" , linewidth = 0.75 , span = 1 ) + labs (x = "Islamic win margin" , y = "Islamic vote share" ,subtitle = "LATE = 0.006 (p = 0.713)" ) + coord_cartesian (ylim = c (0 , 0.65 )) + theme_classic ()|> filter (iwm >= - 0.32 & iwm <= 0.32 ) |> ggplot (aes (x = iwm, y = pop_old, group = treat)) + geom_vline (xintercept = 0 , linetype = "dashed" , color = "gray" ) + stat_summary_bin (fun = "mean" , geom = "point" , binwidth = 0.08 , boundary = 0 ,size = 3 , shape = 21 , color = "black" , fill = "gray" ) + geom_smooth (color = "black" , linewidth = 0.75 , method = "loess" , span = 1 ) + labs (x = "Islamic win margin" , y = "Age 60+" ,subtitle = "LATE = -0.004 (p = 0.445)" ) + coord_cartesian (ylim = c (0.05 , 0.11 )) + theme_classic ()|> filter (iwm >= - 0.32 & iwm <= 0.32 ) |> ggplot (aes (x = iwm, y = pop_young, group = treat)) + geom_vline (xintercept = 0 , linetype = "dashed" , color = "gray" ) + stat_summary_bin (fun = "mean" , geom = "point" , binwidth = 0.08 , boundary = 0 ,size = 3 , shape = 21 , color = "black" , fill = "gray" ) + geom_smooth (color = "black" , linewidth = 0.75 , method = "loess" , span = 1 ) + labs (x = "Islamic win margin" , y = "Age 19-" ,subtitle = "LATE = 0.008 (p = 0.567)" ) + coord_cartesian (ylim = c (0.35 , 0.55 )) + theme_classic ()|> filter (iwm >= - 0.32 & iwm <= 0.32 ) |> ggplot (aes (x = iwm, y = log_pop, group = treat)) + geom_vline (xintercept = 0 , linetype = "dashed" , color = "gray" ) + stat_summary_bin (fun = "mean" , geom = "point" , binwidth = 0.08 , boundary = 0 ,size = 3 , shape = 21 , color = "black" , fill = "gray" ) + geom_smooth (color = "black" , linewidth = 0.75 , method = "loess" , span = 1 ) + labs (x = "Islamic win margin" , y = "Log polulation" ,subtitle = "LATE = 0.002 (p = 0.996)" ) + coord_cartesian (ylim = c (6 , 10 )) + theme_classic ()
いずれの共変量においても目に見えるジャンプは見えず、推定値ベースでみても統計的有意なジャンプは確認できなかった。
感度分析
感度分析とはデータ分析やモデリングにおいて、パラメーターやモデルを変化させたときに、出力結果がどの程度変化するかを調べる手法だ。たとえば、今回のメインの結果はバンド幅が約0.096、三角形カーネル関数、対数人口、若年層と高齢層の人口割合などの共変量を入れたモデルから得られた結果だ。もし、バンド幅が2倍になったらどうなるだろうか。バンド幅が半分になったら?カーネル関数をエパネチニコフ（Epanechnikov）関数にしたら?まったく同じ結果にはならないだろうが、これらの設定を少しいじるだけで結果が大きく変わったら（統計的有意性が変わったり、符号が変わったり）、今回の結果は「この設定でたまたま得られたもの」となり、信頼性の低い推定値になる。
したがって、RDDの推定結果が分析方法やモデルを少し変えてもどの程度安定しているかを確認することは非常に重要である（実はあらゆるデータ分析において重要である）。ここでは論文等でよく見られる感度分析の手法をMeyersson（2014）のデータを使ってやってみる。具体的にはバンド幅、局所回帰における多項式の次数、カーネル関数、共変量、そして近年注目されているドーナツRDDについて解説する。
バンド幅
RDDの感度分析の定番はバンド幅である。まず、メインの分析結果（rdd_female_2000）のバンド幅を確認してみよう。summary(rdrobustオブジェクト名)でもバンド幅は出力されるが、小数点3桁まで表示されない。rdrobustオブジェクト名$bwsと入力すれば、より正確なバンド幅が分かる。
left right
h 0.09593258 0.09593258
b 0.20821751 0.20821751
左側と右側のバンド幅が2列で表示されるが、通常のRDDだとバンド幅は左右対称だからどこを見ても良い。そして、hとb行があるが、通常、我々が言っているバンド幅はhだ。bはバイアス補正のために使用するもう一つのバンド幅であり、論文等で報告するのは基本的にhの方だ。今回のバンド幅は約0.096である。
それではこのバンド幅をh関数の実引数として割り当ててみよう。rdd_female_2000と同じ結果が得られるだろうか。
rdrobust (y = female_hs_2000, x = iwm, data = df, covs = cov_formula,cluster = province, vce = "CR2" , h = 0.09593258 ) |> summary ()
Call: rdrobust
Covariate-adjusted Sharp RD estimates using local polynomial regression. Std. errors are clustered (120 clusters).
Number of Obs. 2629
BW type Manual
Kernel Triangular
VCE method CR2
Left Right
Number of Obs. 2314 315
Eff. Number of Obs. 282 186
Order est. (p) 1 1
Order bias (q) 2 2
BW est. (h) 0.096 0.096
BW bias (b) 0.096 0.096
rho (h/b) 1.000 1.000
Clusters (g) 67 53
Unique Obs. 2312 315
=====================================================================
Point Robust Inference
Estimate z P>|z| [ 95% C.I. ]
---------------------------------------------------------------------
RD Effect 0.025 2.615 0.009 [0.009 , 0.063]
=====================================================================
rdd_female_2000の点推定値は0.025で、p 値は0.003だった。一致しない。なぜだろうか。それはバイアス補正に使われるbが原因である。
rdd_female_2000のbは0.208であるが、hを指定した場合、bはhと同じ値をとる。だから、hだけを動かしながら感度分析することは推奨しない。
バイアス補正のためにバンド幅が2つ用意される。推定用の\(h\) とバイアス補正用の\(b\) 。もしhのみ調整するとbが機能したい。だから、計算された最適バンド幅を入れても結果が一致しない。なので\(\rho = \frac{h}{b}\) を計算し、この\(\rho\) が一定になるように\(b\) を計算し、b引数で指定する（\(b = \frac{h}{\rho}\) ）。今回は\(\rho\) が0.461だから、新しい\(h\) に0.461を割った値をbに指定する。つまり、\(\rho\) を0.461に維持しながら\(h\) を動かすと、\(b\) も自然に動くことになる。
# h と rho の正確な値は rdd_female_2000$bws をベースに計算した rdrobust (y = female_hs_2000, x = iwm, covs = cov_formula,data = df, cluster = province, vce = "CR2" , h = 0.09593258 , b = 0.09593258 / 0.4607325 ) |> summary ()
Call: rdrobust
Covariate-adjusted Sharp RD estimates using local polynomial regression. Std. errors are clustered (120 clusters).
Number of Obs. 2629
BW type Manual
Kernel Triangular
VCE method CR2
Left Right
Number of Obs. 2314 315
Eff. Number of Obs. 282 186
Order est. (p) 1 1
Order bias (q) 2 2
BW est. (h) 0.096 0.096
BW bias (b) 0.208 0.208
rho (h/b) 0.461 0.461
Clusters (g) 67 53
Unique Obs. 2312 315
=====================================================================
Point Robust Inference
Estimate z P>|z| [ 95% C.I. ]
---------------------------------------------------------------------
RD Effect 0.025 2.937 0.003 [0.009 , 0.043]
=====================================================================
それではバンド幅を最適ハンド幅の半分にして推定するためにはどうすれば良いだろうか。まず、hは最適バンド幅である0.09593258を半分にすれば良い。つづいて、bはその新しいhを0.4607325で割った値を指定する。
rdrobust (y = female_hs_2000, x = iwm, covs = cov_formula,data = df, cluster = province, vce = "CR2" , h = 0.09593258 / 2 , b = (0.09593258 / 2 ) / 0.4607325 ) |> rdd_tidy ()
estimate std.error statistic p.value conf.low conf.high
3 0.03695309 0.01330636 2.928097 0.003410436 0.01288232 0.06504229
同じやり方でバンド幅を2倍にして推定してみよう。
rdrobust (y = female_hs_2000, x = iwm, covs = cov_formula,data = df, cluster = province, vce = "CR2" , h = 0.09593258 * 2 , b = (0.09593258 * 2 ) / 0.4607325 ) |> rdd_tidy ()
estimate std.error statistic p.value conf.low conf.high
3 0.0255538 0.007854997 3.253294 0.001140753 0.0101591 0.04095012
バンド幅を半分にしても、2倍にしても、LATEの点推定値は変動するものの、LATEの符号および統計的有意性は変わらなかった。RDDにおいて推定結果に最も大きく影響する要因はバンド幅であるため、バンド幅の感度分析は必須である（少なくとも最適バンド幅の半分と2倍は検証すべき）。
以下の図はバンド幅を0.05から0.2まで0.01ずつ大きくしながらLATEとその95%信頼区間を示した図である。また、最適バンド幅の箇所にはpointrangeを加えた。以下の図を効率よく作るためには{purrr}のmap()関数の使い方を理解する必要がある（知らなくても作れるが、かなりコードが長くなる）。興味のある人は『私たちのR』の「反復処理 」を読んでみよう。
コード
seq (0.05 , 0.2 , by = 0.01 ) |> enframe () |> mutate (rdd = map (value, \(x) rdrobust (y = female_hs_2000, x = iwm, data = df, covs = cov_formula, h = x, b = x / 0.4607325 ,cluster = province, vce = "CR2" )),rdd = map (rdd, \(x) rdd_tidy (x))) |> unnest (rdd) |> ggplot (aes (x = value, y = estimate, ymin = conf.low, ymax = conf.high)) + geom_hline (yintercept = 0 ) + geom_ribbon (alpha = 0 , linetype = "dashed" , color = "black" ) + geom_line () + annotate ("pointrange" , color = "tomato" ,x = 0.096 , y = 0.02482684 , ymin = 0.008617295 , ymax = 0.04319313 ) + labs (x = "Bandwidth" , y = "Estimated LATE with 95% CI" ) + coord_cartesian (ylim = c (- 0.025 , 0.1 )) + theme_classic ()
高次多項式
通常、RDDのLATEはバンド幅内のケースを対象とし、交差項（処置変数 \(\times\) 割当変数）付きの回帰分析を行う。{rdrobust}のrdrobust()の場合、以下のようなモデルを推定する（\(Y\) 、\(D\) 、\(R\) はそれぞれ結果変数、処置変数、割当変数）。
\[
\textsf{for}\,\, -h \leq X \leq h:\,\,\,\,\widehat{Y}_i = \alpha + \tau D_i + \beta R_i + \gamma D_i \cdot R_i \quad
\]
バンド幅内では割当変数と結果変数間の関係がシンプルな関数で表現できるケースが多いため、ほとんどの場合、一次関数で十分だろう。もし、一次関数だけではバンド幅内の関係性がうまく説明できなさそうであれば、二次関数や三次関数のような多項式モデルを推定することになる。もう一つのケースとしては感度分析としての多項式だ。バンド幅内を設定すると、必然的にサンプルサイズが小さくなってしまう。モデルのバイアス—分散トレードオフ関係を考えると、バイアスは小さくなると考えられるが、その一方、モデルの分散が大きくなりがちだ。若干のパラメーターの変動が結果を歪ませる可能性は否定できない。だから、感度分析を通じて、多項式でも推定値が安定することを示す必要がある。
割当変数の二乗項を投入した局所回帰のモデルは以下の通りだ。
\[
\textsf{for}\,\, -h \leq X \leq h:\,\,\,\,\widehat{Y}_i = \alpha + \tau D_i + \beta_1 R_i + \beta_2 R_i^2 + \gamma_1 D_i \cdot R_i + \gamma_2 D_i \cdot R_i \quad
\]
2つを比較すると以下のようなものになる。
ある程度の変動は見られるが、推定値が大きく変わるレベルの変動ではないことがわかる。{rdrobust}のrdrobust()で局所回帰分析における多項式の次数を変えるためにはp引数を使用する。既定値は1だが、これを2にすると二次（quadratic）関数、3にすると三次（cubic）関数、4にすると四次（quartic）関数が回帰モデルとして使われる。以下はバンド幅内で二次関数 + 交差項でLATEを推定する例だ。
rdrobust (y = female_hs_2000, x = iwm, data = df, covs = cov_formula, p = 2 ,cluster = province, vce = "CR2" ) |> summary ()
Call: rdrobust
Covariate-adjusted Sharp RD estimates using local polynomial regression. Std. errors are clustered (120 clusters).
Number of Obs. 2629
BW type mserd
Kernel Triangular
VCE method CR2
Left Right
Number of Obs. 2314 315
Eff. Number of Obs. 365 217
Order est. (p) 2 2
Order bias (q) 3 3
BW est. (h) 0.117 0.117
BW bias (b) 0.186 0.186
rho (h/b) 0.632 0.632
Clusters (g) 67 53
Unique Obs. 2312 315
=====================================================================
Point Robust Inference
Estimate z P>|z| [ 95% C.I. ]
---------------------------------------------------------------------
RD Effect 0.031 2.326 0.020 [0.005 , 0.057]
=====================================================================
LATEが若干大きくなったが、符号および統計的有意性は変わらなかった。このような手順で次数を3、または4まで1ずつ増やしながら推定値の傾向が安定していることを示そう（どこまで試すかに関する議論はないが、3〜4程度で十分だ）。
3〜4パターンの次数しか試さないので表や脚注の文章でも十分だが、以下のように可視化することもできる。
コード
enframe (1 : 4 ) |> mutate (rdd = map (value, \(x) rdrobust (y = female_hs_2000, x = iwm, data = df, covs = cov_formula, p = x,cluster = province, vce = "CR2" )),rdd = map (rdd, \(x) rdd_tidy (x))) |> unnest (rdd) |> ggplot (aes (x = value, y = estimate, ymin = conf.low, ymax = conf.high)) + geom_hline (yintercept = 0 ) + geom_pointrange () + labs (x = "Order of Local Polynomial Regression" , y = "Estimated LATE with 95% CI" ) + coord_cartesian (ylim = c (- 0.025 , 0.1 )) + theme_classic ()
カーネル関数
局所回帰分析に使用するカネール関数はrdrobust()関数のkernel引数で変更できる。既定値は三角形（"triangular"）カーネル関数だが、他にもエパネチニコフ（"epanechnikov"）カーネル関数、長方形・一様（"uniform"）カーネル関数も使える。一様カーネル関数を使う例は一回登場したので、ここではエパネチニコフカーネル関数を使ってみよう。
rdrobust (y = female_hs_2000, x = iwm, data = df, covs = cov_formula, cluster = province, vce = "CR2" , kernel = "epanechnikov" ) |> summary ()
Call: rdrobust
Covariate-adjusted Sharp RD estimates using local polynomial regression. Std. errors are clustered (120 clusters).
Number of Obs. 2629
BW type mserd
Kernel Epanechnikov
VCE method CR2
Left Right
Number of Obs. 2314 315
Eff. Number of Obs. 262 178
Order est. (p) 1 1
Order bias (q) 2 2
BW est. (h) 0.090 0.090
BW bias (b) 0.178 0.178
rho (h/b) 0.506 0.506
Clusters (g) 67 53
Unique Obs. 2312 315
=====================================================================
Point Robust Inference
Estimate z P>|z| [ 95% C.I. ]
---------------------------------------------------------------------
RD Effect 0.023 2.498 0.012 [0.005 , 0.042]
=====================================================================
結果はほぼ変わらない。一様カーネル関数を除くカーネル関数は閾値を頂点とした左右対称の山形である。以下はrdrobust()が提供する三種類のカーネル関数の形を可視化したものである。三角形とエパネチニコフの形はかなり似ていることが分かる。これにより、カーネル関数が選択がLATEの推定結果に影響しにくいことが分かるだろう（一様カーネル関数を除く）。
コード
tibble (x = seq (- 2 , 2 , by = 0.1 ),t = if_else (abs (x) <= 1 , 1 - abs (x), 0 ),e = if_else (abs (x) <= 1 , (3 / 4 ) * (1 - x^ 2 ), 0 ),u = if_else (abs (x) <= 1 , 0.5 , NA )) |> ggplot () + geom_line (aes (x = x, y = t, color = "Triangular" ), linewidth = 1 ) + geom_line (aes (x = x, y = e, color = "Epanechnikov" ), linewidth = 1 ) + geom_line (aes (x = x, y = u, color = "Uniform/Rectangular" ), linewidth = 1 ) + annotate ("segment" , x = - 2 , xend = - 1 , y = 0 , yend = 0 , linewidth = 1 ) + annotate ("segment" , x = 1 , xend = 2 , y = 0 , yend = 0 , linewidth = 1 ) + annotate ("segment" , x = - 1 , xend = - 1 , y = 0 , yend = 0.5 , linewidth = 1 ) + annotate ("segment" , x = 1 , xend = 1 , y = 0 , yend = 0.5 , linewidth = 1 ) + labs (x = "割当変数" , y = "重み" , color = NULL ) + scale_x_continuous (breaks = c (- 1 , 0 , 1 ), labels = c ("-h" , "閾値" , "h" )) + scale_y_continuous (breaks = NULL ) + scale_color_manual (values = c ("Triangular" = "tomato" ,"Epanechnikov" = "royalblue" ,"Uniform/Rectangular" = "black" )) + theme_classic () + theme (legend.position = "bottom" )
一様カーネル関数はRDDではあまり使われていいたいので、カーネル関数の感度分析は省略しても良いかも知れない。もし、付録等で掲載するのであれば、以下のような図で良いだろう。
コード
c ("triangular" , "epanechnikov" , "uniform" ) |> enframe () |> mutate (value = fct_relevel (value, "triangular" , "epanechnikov" , "uniform" )) |> mutate (rdd = map (value, \(x) rdrobust (y = female_hs_2000, x = iwm, data = df, covs = cov_formula, kernel = x,cluster = province, vce = "CR2" )),rdd = map (rdd, \(x) rdd_tidy (x))) |> unnest (rdd) |> ggplot (aes (x = value, y = estimate, ymin = conf.low, ymax = conf.high)) + geom_hline (yintercept = 0 ) + geom_pointrange () + labs (x = "Kernel functions" , y = "Estimated LATE with 95% CI" ) + coord_cartesian (ylim = c (- 0.01 , 0.06 )) + theme_classic ()
共変量の有無
共変量の有無が推定値に影響を与えることもある。本研究のメインの分析結果は共変量を入れており、具体的な推定値は以下の通りだ。
rdd_tidy (rdd_female_2000)
estimate std.error statistic p.value conf.low conf.high
3 0.02482684 0.008820529 2.936923 0.003314864 0.008617295 0.04319313
共変量を入れなかった場合のLATEはどうだろうか。
rdrobust (y = female_hs_2000, x = iwm, data = df,cluster = province, vce = "CR2" ) |> rdd_tidy ()
estimate std.error statistic p.value conf.low conf.high
3 0.02950802 0.01973107 1.479906 0.1388982 -0.009472052 0.06787233
このように結果が変わる。しかも今回は共変量がなかったら、LATEは統計的に有意ではなくなる。これは本研究の知見が共変量の有無に敏感であることを意味する。
また、共変量の有無以外にも、どの共変量を入れるかに関する感度分析もある。「\(X_1\) と\(X_2\) で統制すれば統計的有意性が変わらないけど、さらに\(X_3\) も入れると統計的有意性が変わる」みたいな議論もできる。また、データ分析全般においてよく採用される感度分析では、架空の共変量を想定するケースもあるが、これもまた共変量の感度分析の一種だとも言える。
回帰分析で代表されえる多変量解析を用いた因果推論の際、検証不可能でありながら、かなり致命的な懸念事項が欠落変数バイアス（omitted variable bias; OVB）だ。処置変数と結果変数両方に影響する未特定・未観察の変数がある場合、OVBが発生する。しかし、頑張って交絡要因を特定し、観察、測定しても私たちが想像もできなかった重要な交絡要因があるかも知れない。ここで感度分析が登場する。大雑把にいえば処置変数と結果変数両方と相関を持つ架空の共変量を仮定し、この相関の強さがどれくらいになれば、今回得られた結果が否定されるかを調べる方法だ。そこで「私が見落とした変数があるとしても、相関の相当強いものじゃない限り、今回の結果が覆ることはないだろう」という結果がベストだ。むろん、「私が見落とした変数があり、もしそれが処置・結果両方に少しでも関係するものであれば、今回同様の結果は得られないだろう」と結論づけ、結果の解釈を慎重に行わざるを得ないケースもあろう。
ドーナツRDD
RDDは閾値周辺の個体の重みが高く付くカーネル関数を使っているため、閾値周辺の異常な値（outlier）がある場合、その異常値の影響を強く受けている可能性をある。これは、今回得られたLATEが閾値周辺の異常値の影響を受けており、この異常値がなければ異なる結果が得られる可能性を示唆する。この閾値周辺の異常値の影響を検証することがDount RDDである。
Donut RDDはBarreca et al. （2011）で提案され、Noack and Rothe（2023）によって方法論的にまとめられるようになったが、簡単にいえば、以下のように、閾値周辺の個体を除外した状態でRDD推定する手法だ。むろん、通常のRDD同様、バンド幅外の個体も推定には使用しない。
問題はどれくらいの個体を落とすかだ。このドーナツ・ホール（donut hole）の大きさに関するガイドラインはなく、通常、最適バンド幅の\(\bigcirc\) %という基準を採用することが多い。今回は最適バンド幅の5%のデータを落としてみよう。iwmの絶対値（abs(iwm)）が最適バンド幅の5%（0.09593258 \(\times\) 0.05）より大きい個体だけを残し、df_donutと名付ける。
<- df |> filter (abs (iwm) >= 0.09593258 * 0.05 )dim (df_donut)
25個の個体が除外された（dfの行数は2629行）。 図 6 は-0.1〜0.1の範囲内のiwmをヒストグラムで示したものである。棒（bin）の幅は0.005である。iwm = 0周辺に若干の個体が残っているが、今回のドーナツ・ホールの大きさが0.0047966であるため、-0.005〜0.005内に個体が若干残っている。
それではこのdf_donutを使用し、LATEを推定してみよう。推定の際はバンド幅を最適バンド幅に固定することには注意されたい。dfとdf_donutは別のデータセットだから、hとb引数でバンド幅を固定しないと、バンド幅が変わってしまう。もしDonut RDDの結果、従来の結果に反するものが出たらこれが異常値のせいなのか、新しいバンド幅のせいなのかが識別できない。ただし、閾値周辺の個体をドロップしてからもう一度バンド幅を計算するアプローチもある。大事なのは再現可能性のために、どの手順・方法でDonut RDDを実行したかをしっかりと示すことだろう。
# rdd_female_2000$bws で正確なバンド幅の値を抽出した rdrobust (y = female_hs_2000, x = iwm, data = df_donut, covs = cov_formula, h = 0.09593258 , b = 0.20821751 ,cluster = province, vce = "CR2" ) |> summary ()
Call: rdrobust
Covariate-adjusted Sharp RD estimates using local polynomial regression. Std. errors are clustered (120 clusters).
Number of Obs. 2604
BW type Manual
Kernel Triangular
VCE method CR2
Left Right
Number of Obs. 2302 302
Eff. Number of Obs. 270 173
Order est. (p) 1 1
Order bias (q) 2 2
BW est. (h) 0.096 0.096
BW bias (b) 0.208 0.208
rho (h/b) 0.461 0.461
Clusters (g) 67 53
Unique Obs. 2300 302
=====================================================================
Point Robust Inference
Estimate z P>|z| [ 95% C.I. ]
---------------------------------------------------------------------
RD Effect 0.025 2.762 0.006 [0.008 , 0.046]
=====================================================================
従来のLATEの点推定値が0.025、95%信頼区間は[0.009, 0.043]だつたが、これに比べるとDonut RDDの結果は大きく変わらず推定結果がかなり安定していることが分かる。
Donut RDDを行う際、最適バンド幅の\(\bigcirc\) %の穴を開けるが、このドーナツ・ホールの大きさも恣意的に決めるものである。できれば最適バンド幅の1%、2%、3%、…でDonut RDDを実行し、結果を報告するのは良いだろう。
コード
c (0.096 * seq (0.01 , 0.1 , by = 0.01 )) |> enframe () |> mutate (data = map (value, \(x) filter (df, abs (iwm) >= x)),rdd = map (data, \(x) rdrobust (y = female_hs_2000, x = iwm, data = x, covs = cov_formula,h = 0.09593258 , b = 0.20821751 ,cluster = province, vce = "CR2" )),rdd = map (rdd, \(x) rdd_tidy (x))) |> unnest (rdd) |> ggplot (aes (x = value, y = estimate, ymin = conf.low, ymax = conf.high,color = p.value <= 0.05 )) + geom_hline (yintercept = 0 ) + geom_pointrange () + labs (x = "Donut hole size (% of bandwidth)" , y = "Estimated LATE with 95% CI" ) + scale_x_continuous (breaks = c (0.096 * seq (0.01 , 0.1 , by = 0.01 )), labels = 1 : 10 ) + scale_color_brewer (palette = "Set1" , direction = - 1 ) + theme_classic () + theme (legend.position = "none" )
一部において統計的非有意なLATEが得られているが、全体的な傾向（トレンド）を考慮するとかなり頑健な結果だと考えられる。
