library(broom)     # For converting models into tables
library(haven)
library(huxtable)  # For side-by-side regression table
library(rdrobust)  # For robust nonparametric RDD  
library(rddensity) # For nonparametric RDD
library(stargazer)
library(tidyverse) # For ggplot
  • 作図用代用フォントの設定
theme_set(theme_gray(base_size = 10, base_family = "HiraginoSans-W3"))  # macOS用
theme_set(theme_gray(base_size = 10, base_family = "Meiryo"))        # Windows用

III. 現職優位効果の推定

  • 前のセッションでは補講を受けたか否かということが後期の試験点数に与える影響をもとめた
  • ここではある総選挙で当選したかどうかということが次の総選挙の結果に与える影響（現職優位）について検討してみる
  • 日本の衆院選の小選挙区で党派的現職優位効果がみられるかどうかを分析した Ariga et al. 2016 の研究を参考にする

リサーチ・クエスチョン

  • 選挙で当選した議員は次の選挙で有利なのか？　
  • 2017年総選挙で自民党候補者が僅差で勝った選挙区（処置群）
  • 2017年総選挙で自民党候補者が僅差で負けた選挙区（統制群）
  • 2021年総選挙においてこれら 2 種類の選挙区で、自民党候補者の得票率がどれだけ異なるのか？

仮説

もし党派的な現職優位効果が存在するなら、自民党候補者が僅差で勝った「処置群」の自民党得票率の方が高いはず 　　

現職優位を推定するリサーチ・デザイン

  • 「補講授業の処置効果」の推定とは異なり、政治家の現職優位効果の推定はやや複雑
  • 使うデータにおける閾値 \(c\) の意味に注意する必要がある　　
  • 分析で使うデータにの作成方法の詳細は下のコラム分析に必要なデータを作成を参照

8. 処置効果の推定手順（現職優位）

RDDの推定手順

Step 1: 処置群と統制群を分けるルールはあるか？ ルールの確認
Step 2: RDDデザインは sharp か fuzzy か？ 散布図で確認
Step 3: カットオフ付近での結果変数 \(Y\) の非連続性 散布図とRDプロットで確認
Step 4: カットオフ付近での割当変数 \(R\) の非連続性 理論的に検討
Step 5: 割当変数 \(R\) の値を自分では選べない ヒストグラム、rdplotdensity()関数、 rdensity()関数で検定
Step 6: 処置効果の推定（パラメトリック・ノンパラメトリック）
Step 7: 推定結果を比較する

分析に必要なデータの作成

  • 分析に必要な変数は次の 3 つ
\(Y_i\) 結果を表す 2021年総選挙での自民党候補者の得票率
\(D_i\) 原因を表す 2017総選挙で現職 (\(D_i=1\)) or 非現職 (\(D_i=0\))
\(R_i\) 割当変数 得票率の差（2017nen ）（自民党候補者ー次点者 or 当選者）

分析に必要なデータを作成

  • RProject フォルダ内に data という名称のフォルダを作成する

  • hr96-21.csv をクリックしてデータをパソコンにダウンロード　　

  • ダウンロードした hr96-21.csv を手動でRProject フォルダ内にある data フォルダに入れる

  • 次のいずれかの方法で hr96-21.csv を読み取る

読み取り方法 1
  • na = "."というコマンドは「欠損値をドットで置き換える」という意味
  • 欠損値を空欄のまま残すと、本来「数値 (numeric)」型のデータが「」文字型 (character)」として認識されるなど、エラーの原因になるため、読み取る時点で事前に対処する
df0 <- read_csv("data/hr96-21.csv",
               na = ".")
読み取り方法 2
  • 読み取った値の日本語が文字化けする場合
  • locale()関数を使って日本語エンコーディング (cp932) を指定する
df0 <- read.csv("data/hr96-21.csv",
               na = ".",
               locale = locale(encoding = "cp932"))
読み取り方法 3
df0 <- read.csv("data/hr96-21.csv",
               na = ".")
  • 読み取ったデータの変数を確認
names(df0)
 [1] "year"          "pref"          "ku"            "kun"          
 [5] "wl"            "rank"          "nocand"        "seito"        
 [9] "j_name"        "gender"        "name"          "previous"     
[13] "age"           "exp"           "status"        "vote"         
[17] "voteshare"     "eligible"      "turnout"       "seshu_dummy"  
[21] "jiban_seshu"   "nojiban_seshu"
  • ku と kun を統合して、district を作る
df3 <- df0 |> 
    mutate(
      district = str_c(ku, kun, sep = "_")
      ) |> 
  select(year, district, wl, rank, seito, j_name, status, voteshare)
names(df3)
[1] "year"      "district"  "wl"        "rank"      "seito"     "j_name"   
[7] "status"    "voteshare"

割当変数 \(R_i\)を作る

  • \(R_i\): 2017年総選挙での得票率の差（自民党候補者ー次点者 or 当選者）
自民党候補者が当選した場合
  • 2017年総選挙データだけを残す
  • 2017年総選挙における当選者 (rank = 1) と次点者 (rank = 2) だけを残す
  • voteshare の値を一つだけ上にずらした変数 vs_lag を作る
  • voteshare から vs_runnerup を引いた値 \(R_i\) を作る
df4 <- df3 |> 
  filter(year == 2017) |> 
  filter(rank < 3) |> 
  mutate(vs_lag = lead(voteshare)) |> 
  mutate(R = voteshare - vs_lag) 
DT::datatable(df4) |>
    DT::formatRound(columns=c('voteshare', 'R'),
                    digits=2)
  • 自民党候補者が小選挙区で当選した選挙区だけを残す
df4 <- df4 |> 
  filter(year == 2017) |> 
  filter(rank < 3) |> 
  mutate(vs_lag = lead(voteshare)) |> 
  mutate(R = voteshare - vs_lag) |> 
  filter(seito == "自民") |> 
  filter(rank == 1) |> 
  select(year, district, wl, rank, seito, j_name, status, voteshare, vs_lag, R)
非自民党候補者が当選した場合
  • 2017年総選挙データだけを残す
  • 2017年総選挙における当選者 (rank = 1) と次点者 (rank = 2) だけを残す
  • voteshare の値を一つだけ下にずらした変数 vs_runnerup を作る
  • voteshare から vs_lag を引いた値 \(R_i\) を作る
df4 <- df3 |> 
  filter(year == 2017) |> 
  filter(rank < 3) |> 
  mutate(vs_lag = lag(voteshare))|> 
  mutate(R = voteshare - vs_lag)
  • 自民党候補者が小選挙区で落選した選挙区だけを残す
df5 <- df4 |> 
  filter(year == 2017) |> 
  filter(rank < 3) |> 
  mutate(vs_lag = lag(voteshare))|> 
  mutate(R = voteshare - vs_lag) |> 
  filter(seito == "自民") |> 
  filter(rank == 2) |> 
  select(year, district, wl, rank, seito, j_name, status, vs_lag, R, voteshare)
データのマージ (df4 + df5)
  • 次の二つのデータフレームをマージする
    自民党候補者が当選した場合(df4)
    非自民党候補者が当選した場合(df5)
df6 <- bind_rows(df4, df5)

割当変数 \(D_i\)を作る

  • \(D_i\): 2017総選挙で現職 (\(D_i = 1\)) or 非現職(\(D_i = 0\))
  • 割当変数 \(R_i\)の値が 0 以上なら 現職 (\(D_i = 1\))、0以下なら非現職 (\(D_i = 0\)) を作る
df6 <- df6 |>
  mutate(D = ifelse(R > 0, "現職", "非現職"))
  • 必要な変数に絞る
df6 <- df6 |> 
  select(district, seito, j_name, R, D, voteshare)
DT::datatable(df6)

\(Y_i\)を作る

  • \(Y_i\): 2021年総選挙での自民党候補者の得票率

  • df_3 から 2021年総選挙データだけを抜きだす

df7 <- df3 |> 
  filter(year == 2021) |> 
  filter(seito == "自民") |> 
  rename(Y = voteshare) |> 
  select(district, seito, j_name, Y)

データのマージ (df6 + df7)

df8 <- full_join(df6, df7, by = "j_name") |> 
  select(district.x, seito.x, j_name, R, D, Y, voteshare)
names(df8)
[1] "district.x" "seito.x"    "j_name"     "R"          "D"         
[6] "Y"          "voteshare" 
df8 <- df8 |> 
  mutate(district = district.x,
         seto = seito.x,
         name = j_name,
         R = R,
         Y = Y) |> 
  select(district, name, Y, R, D, voteshare)
  • あちこちに欠損値がある
  • 欠損値を除外する
df9 <- df8 |> 
  na.omit()
  • 作成されたデータフレームをチェック
DT::datatable(df9) |>
    DT::formatRound(columns=c('Y', 'R', 'voteshare'),
                    digits=2)

df9 に含まれる変数一覧

\(district\) 小選挙区名 2021年総選挙で立候補している小選挙区名
\(name\) 立候補者名 2021年総選挙での自民党候補者の氏名
\(Y_i\) 結果を表す 2021年総選挙での自民党候補者の得票率
\(R_i\) 割当変数 2017年総選挙での得票率の差（自民党候補者ー次点者 or 当選者）
\(D_i\) 原因を表す 2017総選挙で現職か非現職か
\(voteshare\) 得票率 2017総選挙での自民党候補者の得票率
  • df9 の記述統計を確認する
stargazer::stargazer(as.data.frame(df9), 
                     type = "html")
Statistic N Mean St. Dev. Min Max
Y 280 49.899 12.645 23.480 84.070
R 280 10.346 17.810 -31.480 60.080
voteshare 280 48.906 10.623 28.110 85.720
  • サンプルサイズが 681 が 280 まで削減され、完全データができた

8.1 Step 1: 処置群と統制群を分けるルールはあるか？

  • 処置群と統制群が何らかのルールで決められているか？
    →　決められている:
    ・2017年総選挙で小選挙区当選（もしくは復活当選）したら「現職」、それ以外は「非現職」

8.2 Step 2: RDD の種類は sharp か fuzzy か？

  • RDD には sharp と fuzzy の 2 種類ある
  • ここで使う RDD が Sharp なのか Fuzzy なのかを判断するため、x 軸に「割当変数」\(R_i\)、y 軸に「処置変数」\(D_i\)（2017年総選挙で非現職か現職か）の散布図を描いて確かめる
## グラフの文字化け防止
theme_set(theme_bw(base_family = "HiraKakuProN-W3"))
ggplot(data = df9,
       mapping = aes(x = R, 
                     y = D, 
                     color = factor(D))) +
  geom_point(alpha = 0.5, 
             alpha = 0.5,
             position = position_jitter(width = 0, height = 0.15)) +
  geom_vline(xintercept = 0, 
             color = "purple", 
             linetype = 2) +
  scale_color_manual(values = c("blue", "red"))

  • カットオフ(0点) を境界線として「現職の自民党議員」と「非現職の自民党議員」が明確に分かれている
  • \(R_i\) が 0 点の議員は議員の体験がない
  • 数値で確認する
df9 |> 
  group_by(D, R < 0) |> 
  summarize(count = n())
# A tibble: 2 × 3
# Groups:   D [2]
  D      `R < 0` count
  <chr>  <lgl>   <int>
1 現職   FALSE     186
2 非現職 TRUE       94
  • 非現職の自民党議員が 94
  • 現職の自民党議員が 186
  • オーバーラップする自民党議員がいない → RDD デザインは sharp

8.3 Step 3: カットオフ付近での結果変数 (Y) の非連続性

  • カットオフ付近での結果変数 \(Y_i\) の非連続性を確認
  • 仮定 (1): カットオフによって処置の値が変わらない場合、結果変数の値が非連続的なジャンプをしないの確認
  • ここでは非連続性（= ジャンプ）があることを期待

散布図 (c=0)

ggplot(df9, aes(x = R, y = Y, color = factor(D)))+
  geom_vline(xintercept = 0, color = "purple", linetype = 2) +
  geom_point(size = 1, alpha = 0.5) + 
  geom_smooth(data = filter(df9, R <= 0),　method = "lm")　+ 
  geom_smooth(data = filter(df9, R > 0),　method = "lm") +
  scale_color_manual(values = c("blue", "red"))

  • 割当変数 \(R_i\)の カットオフ付近ではわずかながらジャンプが認められる
    →　補講の処置効果が期待できる
  • 「カットオフによって処置の値が変わらない」 = 「カットオフを作るルールはない」
  • 例えば、0 のところに「ジャンプ」があるのは、偶然なのではないかと考えてみる
  • 0 でなければ、0点でのジャンプはない」と仮定したい
    →　そのためには例えば「10 と 20 ではジャンプがない」ことを示せばよい
    →　そうすれば説得力が増す

RD プロット (c=0)

  • RDD では、上記のような通常の散布図に加えて RD プロットと呼ばれる特殊なグラフを使う

  • RD プロットとは、散布図の横軸をいくつかのビン (bin) に分割して、そのビンの中に入る \(Y\) の平均値をドットで示した図

  • rdrobustパッケージの rdplot()関数 を使って RD プロット を表示してみる

library(rdrobust)
rdplot(df9$Y, df9$R, c = 0)
[1] "Mass points detected in the running variable."

  • データをスムーズ化しているので、 \(Y\) にジャンプがあるかどうか確認しやすい
  • プラスのジャンプは認められず、むしろマイナスのジャンプが認められる

→ 仮定 (1): カットオフによって処置の値が変わらない場合、結果変数の値が非連続的なジャンプをしないは満たされたとは言えない

8.4 Step 4: カットオフ付近での割当変数 \(R_i\) の非連続性

  • 仮定(2): カットオフで『非連続的な変化』をするのは処置変数だけの確認

  • カットオフ付近で、結果に影響を与える他の「処置」がないことが必要

  • 例えば、この仮定が破られる例として・・・

  • 2021年総選挙（割当変数：\(R\)）で不正が行われたとする

  • この場合、カットオフである 0 を少し超える票を獲得する自民党当選者が多くなるはず

  • 不正がなければ、カットオフである 0 付近の候補者数はほぼ同じになるはず

  • 2017年総選挙で当選した自民党議員だけが2021年総選挙で不正を行ったとする
    →　この場合、2021年総選挙で多く得票したのは「現職」のせいなのか「選挙不正」のせいなのか区別できない

  • 2021年総選挙で選挙不正のような処置がない

→　仮定 (2): カットオフで「非連続的な変化」をするのは処置変数だけは満たされる

8.5 Step 5: 割当変数 \(R\) の値を自分では選べない

  • 仮定 (3): 割当変数の値を自分では選べないの確認

  • ヒストグラム、rdplotdensity()関数、rdensity()関数で検定

ヒストグラム

  • 自民党議員たちは当選するかどうかということを自分で選ぶことはできない
  • カットオフ付近で、割当変数の分布に不自然な変化がないかどうか確認する必要がある
  • 割当変数 \(R_i\)のカットオフ付近での非連続性を確認
  • ここでは非連続性（= ジャンプ）がないことを期待
  • 前期試験（割当変数：R）で不正が行われたとする
  • この場合、カットオフである 0 を少し超える票を獲得する自民党議員が多くなるはず（= 非連続になるはず）
  • 不正がなければ、カットオフである 0 付近の自民党議員数はほぼ同じ（= 連続）になるはず
  • もし、カットオフ前後で自民党議員数に大きな差がある（= 非連続）ようならそれは問題
  • カットオフ付近のヒストグラムを描く
ggplot(df9, aes(x = R, fill = D)) +
  geom_histogram(binwidth = 5, boundary = 0, color = "white") +
  geom_vline(xintercept = 0, color = "purple", linetype = 2) + scale_fill_manual(values = c("blue", "red"))

  • ギリギリで当選した自民党議員数（カットオフの右脇）の方が、ギリギリで落選した自民党議員数（カットオフの左脇）よりも圧倒的に少ない
  • 選挙不正はないと想定できる

rdplotdensity()関数で検定

  • 目視だけだと統計的有意性が分からないので、rdplotdensity()関数 を使って95%信頼区間を表示してみる
df9_modified <- df9 |> 
  filter(R > -30)　# マイナスの外れ値を除外  
rdplotdensity(rdd = rddensity(df9_modified$R, c = 0),
              X = df9_modified$R,
              type = "both") |> 
  summary()

        Length Class     Mode
Estl     4     lpdensity list
Estr     4     lpdensity list
Estplot 10     gg        list
  • カットオフ前後の 95%信頼区間がオーバーラップしている
    → カットオフである 0 付近の自民党候補者数に差がない
    → 連続している
  • カットオフである 0 を少し超える自民党候補者は増えていない
  • むしろ数が少ない
    →　選挙の不正操作はないと思われる

rddensity()関数で検定

  • 割当変数の操作の確認は rddensity パッケージの rddensity()関数を使っても確認できる
library(rddensity)
summary(rddensity(df9$R, c = 0))

Manipulation testing using local polynomial density estimation.

Number of obs =       280
Model =               unrestricted
Kernel =              triangular
BW method =           estimated
VCE method =          jackknife

c = 0                 Left of c           Right of c          
Number of obs         94                  186                 
Eff. Number of obs    66                  55                  
Order est. (p)        2                   2                   
Order bias (q)        3                   3                   
BW est. (h)           13.088              11.35               

Method                T                   P > |T|             
Robust                -1.3672             0.1715              


P-values of binomial tests (H0: p=0.5).

Window Length              <c     >=c    P>|T|
1.310     + 1.310          16       5    0.0266
2.619     + 2.426          19       9    0.0872
3.927     + 3.541          27      10    0.0076
5.236     + 4.657          35      14    0.0038
6.545     + 5.772          43      19    0.0032
7.853     + 6.888          46      22    0.0049
9.162     + 8.003          52      26    0.0043
10.470    + 9.119          58      35    0.0220
11.779    + 10.234         66      45    0.0572
13.088    + 11.350         66      55    0.3634
  • ここでの帰無仮説：「カットオフ前後の平均値は連続している」
  • p-value = 0.1715 →　帰無仮説は棄却されず
    → 連続している　→　割り当て変数 \(R\) に操作は行われていない

・仮定 (3): 割当変数の値を自分では選べないは満たされる

9. パラメトリック RD（現職優位）

  • Step 6: 処置効果を推定するは「9. パラメトリック RD（現職優位）」と「10. ノンパラメトリック RD（現職優位）」で実行する
  • ここでは、現職優位効果の回帰関数として以下の 4 つの関数形を考えて、まずパラメトリック回帰をやってみる

9.1 推定に使う関数形（現職優位）

\(1. Y_i = α + ρD_i + βR_i + ε_i\)
\(2. Y_i = α + ρD_i + βR_i + γ(D_i×R_i) + ε_i\)
\(3. Y_i = α + ρD_i + β_1R_i + β_2R_i^2 + ε_i\)
\(4. Y_i = α + ρD_i + β_1R_i + β_2R_i^2+ γ_1(D_i×R_i)+ γ_2(D_i×R_i^2) + ε_i\)

  • 通常、パラメトリック RD では、割当変数 (R) を中心化した変数 RC を作る
    →　しかし、ここでは割当変数 (R) のカットオフ値は 0
    →　中心化する必要がない
    →　回帰分析では割当変数 (R) を使う
  • R の値を 2 乗した R_sq を作る　　
  • 処置変数 D の値（「非現職」「現職」）をそれぞれ「0」「1」に変更したダミー変数 D1を作る
  • ここで作成した 2 つの変数 R_sq と D1をデータフレーム df9 に追加する
df9 <- df9 |> 
  mutate(R_sq = R^2) |>                  # R を2乗した R_sqの作成
  mutate(D1 = ifelse(D == "現職", 1, 0)) # 現職ダミー D1 の作成
  • データの中身を確認する
stargazer(as.data.frame(df9), type = "html")
Statistic N Mean St. Dev. Min Max
Y 280 49.899 12.645 23.480 84.070
R 280 10.346 17.810 -31.480 60.080
voteshare 280 48.906 10.623 28.110 85.720
R_sq 280 423.107 606.427 0.001 3,609.606
D1 280 0.664 0.473 0 1
  • ①から④の回帰式を lm推定する formula としてまとめて MODELS_9 と名前を付ける
MODELS_9 <- alist(
  model1 = Y ~ D1 + R,
  model2 = Y ~ D1 * R,
  model3 = Y ~ D1 + R + R_sq,
  model4 = Y ~ D1 * R + D1 * R_sq) |> 
  enframe(name = "modelname", value = "formula")
  • MODELS_9 の中身を確認
MODELS_9
modelnameformula
model1Y ~ D1 + R
model2Y ~ D1 * R
model3Y ~ D1 + R + R_sq
model4Y ~ D1 * R + D1 * R_sq
  • 4 つの回帰式を、一挙に推定
  • 推定結果を rdd_9 として保存
rdd_9 <- MODELS_9 |> 
  mutate(model = map(.x = formula, .f = lm, data = df9)) |> 
  mutate(pred = map(.x = model, .f = predict),   # 予測値を計算する
         result = map(.x = model, .f = tidy))
names(rdd_9)
[1] "modelname" "formula"   "model"     "pred"      "result"
  • rdd_9 に保存されている結果を 1 つずつ確認する

9.1.1 model1 の推定結果

  • model1 による推定結果を図示してみる
p1_rdd2 <- df9 |> 
  mutate(pred = rdd_9$pred[[1]]) |> 
  ggplot(aes(x = R, color = as.factor(D1))) +
  geom_point(aes(y = Y)) +
  geom_line(aes(y = pred)) +
  geom_vline(xintercept = 0, color = "gray")+
  labs(x = "得票率の差（2017年総選挙）", y = "得票率（2021年総選挙）") +
  theme(legend.position = "right") +
  scale_color_manual(values = c("blue", "red")) 
plot(p1_rdd2)

  • 割当変数 \(R\) のカットオフ 0 の時に負のジャンプがある
  • この分析は、処置変数と割当変数のみを説明変数として利用する、最も単純な パラメトリック RD
  • そのため、予測値は「直線」であり、割当変数が 0 未満と 0 以上について直線の傾きが同じ
  • このときの処置効果の推定値と標準誤差を計算してみる
rdd_9 |> 
  unnest(cols = result) |> 
  filter(modelname == "model1", term == "D1") |> 
  select(term, estimate, std.error, p.value) |> 
  knitr::kable(digits = 2)
term estimate std.error p.value
D1 -1.67 2.02 0.41
  • 処置効果は - 1.67 という結果を得た

9.1.2 model2 の推定結果

  • 処置変数と割当変数の交差項を含む model2 の結果を図示してみる
p2_rdd2 <- df9 |> 
  mutate(pred = rdd_9$pred[[2]]) |> 
  ggplot(aes(x = R, color = as.factor(D1))) +
  geom_point(aes(y = Y)) +
  geom_line(aes(y = pred)) +
  geom_vline(xintercept = 0, color = "gray") +
  labs(x = "得票率の差（2017年総選挙）", y = "得票率（2021年総選挙）") +
  theme(legend.position = "right") +
  scale_color_manual(values = c("blue", "red")) 
plot(p1_rdd2)

  • 処置効果の推定値と標準誤差を計算してみる
rdd_9 |> 
  unnest(cols = result) |> 
  filter(modelname == "model2", term == "D1") |> 
  select(term, estimate, std.error, p.value) |> 
  knitr::kable(digits = 2)
term estimate std.error p.value
D1 0.92 2.15 0.67
  • 交差項を含んだ model2 の処置効果は 0.92 という結果を得た

9.1.3 model3 の推定結果

  • model1 に割当変数の 2 乗項を加えた model3 の推定結果を図示してみる
p3_rdd2 <- df9 |> 
  mutate(pred = rdd_9$pred[[3]]) |> 
  ggplot(aes(x = R, color = as.factor(D1))) +
  geom_point(aes(y = Y)) +
  geom_line(aes(y = pred)) +
  geom_vline(xintercept = 0, color = "gray") +
  labs(x = "得票率の差（2017年総選挙）", y = "得票率（2021年総選挙）") +
  theme(legend.position = "right") +
  scale_color_manual(values = c("blue", "red")) 
plot(p3_rdd2)

  • 2 乗項を使った結果、予測値が直線ではなく「曲線」として推定された
  • 交差項がないので、曲線は現職と非現職が同じ傾きと仮定
  • この関数形ではジャンプが確認できる
  • 処置効果の推定値と標準誤差を計算してみる
rdd_9 |> 
  unnest(cols = result) |> 
  filter(modelname == "model3", term == "D1") |> 
  select(term, estimate, std.error, p.value) |> 
  knitr::kable(digits = 2)
term estimate std.error p.value
D1 2.73 2.41 0.26
  • 2 乗項を含む model3 の処置効果は 0.92 という結果を得た

9.1.4 model4 の推定結果

  • 最後に、交差項と 2 乗項を含む model4 の推定結果を可視化してみる

  • 比較のため、model1 で推定した直線も点線として表示する　　

p4_rdd2 <- df9 |> 
  mutate(pred0 = rdd_9$pred[[1]],
         pred = rdd_9$pred[[4]]) |> 
  ggplot(aes(x = R, color = as.factor(D1))) +
  geom_point(aes(y = Y)) +
  geom_line(aes(y = pred0), linetype = "dashed") +
  geom_line(aes(y = pred)) +
  geom_vline(xintercept = 0, color = "gray") +
  labs(x = "得票率の差（2017年総選挙）", y = "得票率（2021年総選挙）") +
  theme(legend.position = "right") +
  scale_color_manual(values = c("blue", "red")) 

plot(p4_rdd2)

  • 交差項と 2 乗項の両方を使ったので、予測値が曲線になりカットオフ前後で異なる曲線が描かれている
  • ここでは正のジャンプがある
  • 処置効果の推定値と標準誤差を計算してみる
rdd_9 |> 
  unnest(cols = result) |> 
  filter(modelname == "model4", term == "D1") |> 
  select(term, estimate, std.error, p.value) |> 
  knitr::kable(digits = 2)
term estimate std.error p.value
D1 3.37 3 0.26
  • model1 〜 model4 の推定結果のまとめ 推定した 4 つのモデルの推定結果を比べてみる
rdd_9 |> 
  unnest(cols = result) |> 
  filter(term == "D1") |> 
  select(modelname, term, estimate, std.error, p.value) |> 
  knitr::kable(digits = 2)
modelname term estimate std.error p.value
model1 D1 -1.67 2.02 0.41
model2 D1 0.92 2.15 0.67
model3 D1 2.73 2.41 0.26
model4 D1 3.37 3.00 0.26

10. ノンパラメトリック RD（現職優位）

  • データの中身を確認する
stargazer(as.data.frame(df9), type = "html")
Statistic N Mean St. Dev. Min Max
Y 280 49.899 12.645 23.480 84.070
R 280 10.346 17.810 -31.480 60.080
voteshare 280 48.906 10.623 28.110 85.720
R_sq 280 423.107 606.427 0.001 3,609.606
D1 280 0.664 0.473 0 1
  • ここで使う変数は観察可能な 2 つの変数 \(R_i, Y_i\)　
\(Y\) 結果を表す 2021年総選挙での自民党候補者の得票率
\(R\) 割当変数 2017年総選挙での得票率の差（自民党候補者ー次点者 or 当選者）

10.1 散布図と RD プロットを描く

  • 最初に、割り付け変数 \(R\) を \(x\) 軸、結果変数 \(Y\) を \(y\) 軸とした散布図と RDプロットを描く
df9 |> 
  ggplot(aes(R, Y)) +
  geom_point() +
  stat_smooth()  # 線形ではなく LOWESS 曲線を描く

  • データをスムーズ化し、\(Y\) にジャンプがあるかどうか確認しやすい RDプロットを描いてみる
library(rdrobust)
rdplot(df9$Y, df9$R, c = 0)
[1] "Mass points detected in the running variable."

  • x = 0 のカットオフでマイナスのジャンプが確認できる

10.2 rdrobust パッケージを使った分析

  • RDD を実行する前に必要な仮定の確認は終わっているので、ここでは Step 6-2 を実行する
Step 1: 処置群と統制群を分けるルールはあるか？ 8.1 で確認済み
Step 2: RDDデザインは sharp か fuzzy か？ 8.2 で確認済み
Step 3: カットオフ付近での結果変数 \(Y\) の非連続性 8.3 で確認済み
Step 4: カットオフ付近での割当変数 \(R\) の非連続性 8.4 で確認済み
Step 5: 割当変数 \(R\) の値を自分では選べない 8.5 で確認済み
Step 6-1: 処置効果の推定（パラメトリック） 9 で実行済み
Step 6-2: 処置効果の推定（ノンパラメトリック）
Step 7: 推定結果を比較する
  • ここでは次の複数のバンド幅をモデリングして推定する
    ・IK バンド幅
    ・CER 最適バンド幅
    ・MSE 最適化バンド幅
  • まず、Imbens and Kalyanarman (2012) が推奨するバンド幅を確認してみる
IK バンド幅
IKband2 <- rdbwselect_2014(df9$Y, df9$R, c = 0, bwselect = "IK")
IKband2$bws
            h        b
[1,] 15.04033 10.33074
  • Imbens and Kalyanarman (2012) は 15.04033 が最適なバンド幅だと推奨している
  • では、この推奨されたバンド幅を使って、補講授業の処置効果を推定してみる
model_IK2 <- rdrobust(df9$Y, df9$R, c = 0, h = IKband2$bws[1,1], all = TRUE)
[1] "Mass points detected in the running variable."
summary(model_IK2)
Call: rdrobust

Number of Obs.                  280
BW type                      Manual
Kernel                   Triangular
VCE method                       NN

Number of Obs.                   94          186
Eff. Number of Obs.              72           79
Order est. (p)                    1            1
Order bias  (q)                   2            2
BW est. (h)                  15.040       15.040
BW bias (b)                  15.040       15.040
rho (h/b)                     1.000        1.000
Unique Obs.                      49          184

=============================================================================
        Method     Coef. Std. Err.         z     P>|z|      [ 95% C.I. ]       
=============================================================================
  Conventional     1.636     2.432     0.673     0.501    [-3.130 , 6.403]     
Bias-Corrected    -2.087     2.432    -0.858     0.391    [-6.854 , 2.680]     
        Robust    -2.087     2.787    -0.749     0.454    [-7.549 , 3.375]     
=============================================================================
  • IK によって推奨されたバンド幅 (15.04033) の時の推定値は 1.636
CER 最適バンド幅
model_cer2 <- rdrobust(df9$Y, df9$R, c = 0, bwselect = "cerrd", all = TRUE)
[1] "Mass points detected in the running variable."
summary(model_cer2)
Call: rdrobust

Number of Obs.                  280
BW type                       cerrd
Kernel                   Triangular
VCE method                       NN

Number of Obs.                   94          186
Eff. Number of Obs.              19            9
Order est. (p)                    1            1
Order bias  (q)                   2            2
BW est. (h)                   3.006        3.006
BW bias (b)                   6.593        6.593
rho (h/b)                     0.456        0.456
Unique Obs.                      49          184

=============================================================================
        Method     Coef. Std. Err.         z     P>|z|      [ 95% C.I. ]       
=============================================================================
  Conventional    -8.049     4.185    -1.923     0.054   [-16.251 , 0.154]     
Bias-Corrected    -7.965     4.185    -1.903     0.057   [-16.167 , 0.238]     
        Robust    -7.965     4.681    -1.701     0.089   [-17.140 , 1.210]     
=============================================================================
MSE 最適化バンド幅
model_mse2 <- rdrobust(df9$Y, df9$R, c = 0, bwselect = "mserd", all = TRUE)
[1] "Mass points detected in the running variable."
summary(model_mse2)
Call: rdrobust

Number of Obs.                  280
BW type                       mserd
Kernel                   Triangular
VCE method                       NN

Number of Obs.                   94          186
Eff. Number of Obs.              29           11
Order est. (p)                    1            1
Order bias  (q)                   2            2
BW est. (h)                   3.985        3.985
BW bias (b)                   6.593        6.593
rho (h/b)                     0.604        0.604
Unique Obs.                      49          184

=============================================================================
        Method     Coef. Std. Err.         z     P>|z|      [ 95% C.I. ]       
=============================================================================
  Conventional    -5.183     3.731    -1.389     0.165   [-12.495 , 2.129]     
Bias-Corrected    -5.758     3.731    -1.543     0.123   [-13.070 , 1.554]     
        Robust    -5.758     4.583    -1.256     0.209   [-14.739 , 3.224]     
=============================================================================

10.3 分析結果のまとめ

  • 「補講の処置効果」に関して、ここまで実施した分析結果をまとめてみる

ノンパラメトリックRDD の推定結果

手法 推定値 標準誤差 p-value
パラメトリックmodel1 -1.67 2.02 0.41
パラメトリックmodel2 0.92 2.15 0.67
パラメトリックmodel3 2.73 2.41 0.26
パラメトリックmodel4 3.37 3.00 0.26
ノンパラIK バンド幅 1.636 2.432 0.501
ノンパラCER -8.049 4.185 0.054
ノンパラMSE -5.183 3.731 0.165
  • パラメトリック RD もノンパラメトリック RDも得られた推定値も符号もバラバラで統計的にも有意でない
現職優位の処置効果の推定結果 2017年と2021年総選挙における自民党候補者に関しては、現職優位があるとはいえない

IV. 法定飲酒年齢効果の推定

  • ここでは法定飲酒年齢 (minimum legal drinking age; MLDA) が21歳の誕生日付近の死亡者数に与える影響を分析する

Research Question:

法定飲酒年齢の引き下げは、18歳から20歳の若者の死亡率を上げたか？
  • Angrist and Pischke (2009)による研究
  • アメリカの法定飲酒年齢 (MLDA: minimum legal drinking age)
  • 州ごとによって異なる (50州＋ワシントン D.C.)
  • MLDA: 18, 19, 20, 21歳
  • ここでは MLDA が21歳の場合のみを考える
  • 確認したい効果：
  • 法律で飲酒が認められている州では、18-20歳の死亡率が上がるのか?

11. 処置効果の推定手順（法定飲酒年齢）

RDDの推定手順

Step 1: 処置群と統制群を分けるルールはあるか？ ルールの確認
Step 2: RDDデザインは sharp か fuzzy か？ 散布図で確認
Step 3: カットオフ付近での結果変数 \(Y\) の非連続性 散布図とRDプロットで確認
Step 4: カットオフ付近での割当変数 \(R\) の非連続性 理論的に検討
Step 5: 割当変数 \(R\) の値を自分では選べない ヒストグラム、rdplotdensity()関数、 rdensity()関数で検定
Step 6: 処置効果の推定（パラメトリック・ノンパラメトリック）
Step 7: 推定結果を比較する

分析に必要なデータの作成

分析に必要なデータを作成 

download.file(url = "http://masteringmetrics.com/wp-content/uploads/2015/01/AEJfigs.dta",
              dest = "data/AEJfigs.dta")
  • Stata 形式のデータなので、haven::read_dta() で読み込む
MLDA <- read_dta("data/AEJfigs.dta")
  • データフレーム MLDA が含む変数の確認
names(MLDA)
 [1] "agecell"             "all"                 "allfitted"          
 [4] "internal"            "internalfitted"      "external"           
 [7] "externalfitted"      "alcohol"             "alcoholfitted"      
[10] "homicide"            "homicidefitted"      "suicide"            
[13] "suicidefitted"       "mva"                 "mvafitted"          
[16] "drugs"               "drugsfitted"         "externalother"      
[19] "externalotherfitted"
  • 分析に必要な変数だけに限定する
df10 <- MLDA |> 
  select(agecell, all)
agecell： 年齢を記録した変数
all： 死者数（全ての死因）

  • 「28. 差分の差分法:DID」で使っているデータとは異なり、このデータには MLDA が21歳の場合のみが含まれている

  • このデータの観測単位は「誕生日によって区切られた年齢グループ」

  • 例えば、「19歳の誕生日になってから1ヶ月未満の人たち」「19歳の誕生日になってから1ヶ月以上2ヶ月未満の人たち」といいうように、年齢（約30日単位）が1つの観測点を構成している

  • 回帰分析を行うために、age、over_21, age_sq とい 3 つの変数を作成する

df10 <- MLDA |> 
  mutate(age = agecell - 21,    # 21歳を基準に age を中心化した変数
    over_21 = as.integer(agecell > 21), # 21歳以上を示すダミー変数
    age_sq = age^2) |> 
  select(agecell, age, age_sq, over_21, all)
summary(df10)
    agecell           age                 age_sq          over_21    
 Min.   :19.07   Min.   :-1.9315071   Min.   :0.0000   Min.   :0.00  
 1st Qu.:20.08   1st Qu.:-0.9246578   1st Qu.:0.2044   1st Qu.:0.00  
 Median :21.00   Median :-0.0000048   Median :0.8934   Median :0.00  
 Mean   :21.00   Mean   :-0.0000002   Mean   :1.2446   Mean   :0.48  
 3rd Qu.:21.92   3rd Qu.: 0.9246578   3rd Qu.:2.0689   3rd Qu.:1.00  
 Max.   :22.93   Max.   : 1.9315071   Max.   :3.7307   Max.   :1.00  
                                                                     
      all        
 Min.   : 88.43  
 1st Qu.: 92.79  
 Median : 95.69  
 Mean   : 95.67  
 3rd Qu.: 98.03  
 Max.   :105.27  
 NA's   :2
  • 死者数 (all) に欠測 (NA) がある 2 つの観測個体を除外する
df10 <- df10 |> 
  filter(complete.cases(all))
agecell： 年齢を記録した変数
age： 21歳を中心に年齢を中心化（21歳の誕生日からのズレを測る変数）
age_sq： age を 2 乗した変数
over_21： 処置変数：21歳以上を示すダミー変数 
DT::datatable(df10) |>
    DT::formatRound(columns=c('agecell', 'age', 'age_sq', 'all'),
                    digits=2)
stargazer(as.data.frame(df10), 
          type = "html")
Statistic N Mean St. Dev. Min Max
agecell 48 21.000 1.151 19.068 22.932
age 48 0.000 1.151 -1.932 1.932
age_sq 48 1.296 1.171 0.002 3.731
over_21 48 0.500 0.505 0 1
all 48 95.673 3.831 88.428 105.268
all： 死者数（全ての死因）
over_21： 処置変数：21歳以上を示すダミー変数
age： 21歳を中心に年齢を中心化（21歳の誕生日からのズレを測る変数）
agecell： 年齢を記録した変数
age_sq： age を 2 乗した変数

11.1 Step 1: 処置群と統制群を分けるルール確認

  • 処置群と統制群が何らかのルールで決められているかどうか
    →　決められている:
    ・21歳以上なら「法律で飲酒が認められている \(D=1\)」
    ・21歳以下なら「法律で飲酒が認められていない\(D=0\)」

11.2 Step 2: RDD の種類は sharp か fuzzy か？

  • RDD には sharp と fuzzy の 2 種類ある
  • ここで使う RDD が Sharp なのか Fuzzy なのかを判断するため、x 軸に「割当変数」\(R_i\)、y 軸に「処置変数」\(D_i\)（法律で飲酒が認められているか）の散布図を描いて確かめる
## グラフの文字化け防止
theme_set(theme_bw(base_family = "HiraKakuProN-W3"))
ggplot(data = df10,
       mapping = aes(x = agecell, 
                     y = factor(over_21), 
                     color = factor(over_21))) +
  geom_point(alpha = 0.5, 
             alpha = 0.5,
             position = position_jitter(width = 0, height = 0.15)) +
  geom_vline(xintercept = 21, 
             color = "purple", 
             linetype = 2) +
  scale_color_manual(values = c("blue", "red"))

  • カットオフ(21歳) を境界線として「法律で飲酒が認められている 人」（水色）と「法律で飲酒が認められていない人」（トマト色）が明確に分かれている
  • 数値で確認する
df10 |> 
  group_by(over_21, age < 0) |> 
  summarize(count = n())
# A tibble: 2 × 3
# Groups:   over_21 [2]
  over_21 `age < 0` count
    <int> <lgl>     <int>
1       0 TRUE         24
2       1 FALSE        24
  • 21歳未満がが 24人
  • 21歳以上がが 24人
  • オーバーラップする人がいない → RDD デザインは sharp

11.3 Step 3: カットオフ付近での結果変数 (Y) の非連続性

  • カットオフ付近での結果変数 \(Y_i\) の非連続性を確認
  • 仮定 (1): カットオフによって処置の値が変わらない場合、結果変数の値が非連続的なジャンプをしないの確認
  • ここでは非連続性（= ジャンプ）があることを期待

散布図 (c=21)

ggplot(df10, aes(x = agecell, y = all, color = factor(over_21)))+
  geom_vline(xintercept = 21, color = "purple", linetype = 2) +
  geom_point(size = 1, alpha = 0.5) + 
  geom_smooth(data = filter(df10, over_21 <= 0),　method = "lm")　+ 
  geom_smooth(data = filter(df10, over_21 > 0),　method = "lm") +
  scale_color_manual(values = c("blue", "red"))

  • 割当変数 \(R_i\)の カットオフ付近ではジャンプが認められる
    →　法定飲酒年齢の処置効果が期待できる
  • 「カットオフによって処置の値が変わらない」 = 「カットオフを作るルールはない」
  • 例えば、21歳のところに「ジャンプ」があるのは、偶然なのではないかと考えてみる
  • 21歳でなければ、21歳でのジャンプはない」と仮定したい
    →　そのためには例えば「20歳と22歳ではジャンプがない」ことを示せばよい
    →　そうすれば説得力が増す

RD プロット (c=21)

  • RDD では、上記のような通常の散布図に加えて RD プロットと呼ばれる特殊なグラフを使う

  • RD プロットとは、散布図の横軸をいくつかのビン (bin) に分割して、そのビンの中に入る \(Y\) の平均値をドットで示した図

  • rdrobustパッケージの rdplot()関数 を使って RD プロット を表示してみる

library(rdrobust)
rdplot(df10$all, df10$agecell, c = 21)

  • データをスムーズ化しているので、 \(Y\) にジャンプがあるかどうか確認しやすい
  • プラスのジャンプが認められる
  • 「カットオフによって処置の値が変わらない場合」の一例として 20 歳でのジャンプの有無を確認してみる

散布図 (c=20)

df10 |>
  mutate(threshold = as.factor(ifelse(agecell > 20, 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 = 20, color = "black", linetype = 2) +
  scale_color_manual(values = c("blue", "red"))

RD プロット (c=20)

rdplot(df10$all, df10$agecell, c = 20)

- 20歳 では割当変数 \(R_i\)の カットオフ付近よりも小さいジャンプが認められる

  • もう一つ「カットオフによって処置の値が変わらない場合」の例として 22歳 でのジャンプの有無を確認してみる

散布図 (c=22)

df10 |>
  mutate(threshold = as.factor(ifelse(agecell > 22, 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 = 22, color = "black", linetype = 2) +
  scale_color_manual(values = c("blue", "red"))

  • 22歳では割当変数 \(R_i\)の カットオフ付近とは異なり、むしろ負のジャンプが認められる

RD プロット (c=22)

rdplot(df10$all, df10$agecell, c = 22)

  • プラスのジャンプが認められる
    →　カットオフによって処置の値が変わらない場合（20歳と22歳の場合）
    →　結果変数の値が非連続的なジャンプをせず、21歳 の時に最も大きな「ジャンプ」がある

→ 仮定 (1):カットオフによって処置の値が変わらない場合、結果変数の値が非連続的なジャンプをしないは満たされる

11.4 Step 4: カットオフ付近での割当変数 \(Ri\) の非連続性

カットオフで「非連続的な変化」をするのは処置変数だけ

  • カットオフ付近で、結果に影響を与える他の「処置」がないことが必要

  • 例えば、この仮定が破られる例として・・・

  • 酒を飲んでいいのは 21 歳以上と法律で決められているのに、多くの若者が法律を守らずに酒を飲むという事実が認められた場合
    →　この場合、若者の死者数に対する法定飲酒年齢の処置効果が正しく推定できない

  • このようなことがなければ・・・
    →　カットオフで「非連続的な変化」をするのは処置変数だけという仮定は満たされる

11.5 Step 5: 割当変数 \(R\) の値を自分では選べない

  • 仮定 (3): 割当変数の値を自分では選べないの確認

  • ヒストグラム、rdplotdensity()関数、rdensity()関数で検定

ヒストグラム

  • カットオフ付近のヒストグラムを描く
  • 変数の型を確認する
str(df10)
tibble [48 × 5] (S3: tbl_df/tbl/data.frame)
 $ agecell: num [1:48] 19.1 19.2 19.2 19.3 19.4 ...
  ..- attr(*, "label")= chr "Age Cell"
  ..- attr(*, "format.stata")= chr "%9.0g"
 $ age    : num [1:48] -1.93 -1.85 -1.77 -1.68 -1.6 ...
  ..- attr(*, "label")= chr "Age Cell"
  ..- attr(*, "format.stata")= chr "%9.0g"
 $ age_sq : num [1:48] 3.73 3.42 3.12 2.84 2.57 ...
  ..- attr(*, "label")= chr "Age Cell"
  ..- attr(*, "format.stata")= chr "%9.0g"
 $ over_21: int [1:48] 0 0 0 0 0 0 0 0 0 0 ...
 $ all    : num [1:48] 92.8 95.1 92.1 88.4 88.7 ...
  ..- attr(*, "label")= chr "All"
  ..- attr(*, "format.stata")= chr "%9.0g"
  • 変数 over_21 の型が numeric なので数値型 (factor) に変換する
df10$over_21 <- as.factor(df10$over_21)
ggplot(df10, aes(x = agecell, fill = over_21)) +
  geom_histogram(binwidth = 0.3, boundary = 0, color = "white") +
  geom_vline(xintercept = 21, color = "purple", linetype = 2) + scale_fill_manual(values = c("blue", "red"))

  • 21歳前後の若者の数が極端に多い（あるいは少ない）ということはないとわかる

rdplotdensity()関数で検定

rdplotdensity(rdd = rddensity(df10$agecell, c = 21),
              X = df10$agecell,
              type = "both") |> 
  summary()

        Length Class     Mode
Estl     4     lpdensity list
Estr     4     lpdensity list
Estplot 10     gg        list
  • カットオフ前後の 95%信頼区間がオーバーラップしている
    → カットオフである 0 付近の自民党候補者数に差がない
    → 連続している
  • カットオフである 0 を少し超える自民党候補者は増えていない
  • むしろ数が少ない
    →　選挙の不正操作はないと思われる

rddensity()関数で検定

  • 割当変数の操作の確認は rddensity パッケージの rddensity()関数を使っても確認できる
library(rddensity)
summary(rddensity(df10$agecell, c = 21))

Manipulation testing using local polynomial density estimation.

Number of obs =       48
Model =               unrestricted
Kernel =              triangular
BW method =           estimated
VCE method =          jackknife

c = 21                Left of c           Right of c          
Number of obs         24                  24                  
Eff. Number of obs    24                  24                  
Order est. (p)        2                   2                   
Order bias (q)        3                   3                   
BW est. (h)           1.932               1.932               

Method                T                   P > |T|             
Robust                0                   1                   


P-values of binomial tests (H0: p=0.5).

Window Length / 2          <c     >=c    P>|T|
0.781                      10      10    1.0000
0.909                      11      11    1.0000
1.037                      13      13    1.0000
1.164                      14      14    1.0000
1.292                      16      16    1.0000
1.420                      17      17    1.0000
1.548                      19      19    1.0000
1.676                      20      20    1.0000
1.804                      22      22    1.0000
1.932                      24      24    1.0000
  • ここでの帰無仮説：「カットオフ前後の平均値は連続している」
  • p-value = 1 →　帰無仮説は棄却されず
    → 連続している　→　割り当て変数 \(R\) に操作は行われていない

・仮定 (3): 割当変数の値を自分では選べないは満たされる

12. パラメトリック RD（法定飲酒年齢）

  • ここでは、現職優位効果の回帰関数として以下の 4 つの関数形を考えて、まずパラメトリック回帰をやってみる

推定に使う関数形（法定飲酒年齢）

\(1. Y_i = α + ρD_i + βR_i + ε_i\)
\(2. Y_i = α + ρD_i + βR_i + γ(D_i×R_i) + ε_i\)
\(3. Y_i = α + ρD_i + β_1R_i + β_2R_i^2 + ε_i\)
\(4. Y_i = α + ρD_i + β_1R_i + β_2R_i^2+ γ_1(D_i×R_i)+ γ_2(D_i×R_i^2) + ε_i\)

  • 割当変数 \(R\) を中心化すれば、どの関数形においても \(ρ\) （ロー）の推定値が処置効果の推定値となる
  • ここでは割当変数 \(R = agecell\) を中心化した変数として \(age\) を使っている
DT::datatable(df10) |>
    DT::formatRound(columns=c('agecell', 'age', 'age_sq', 'all'),
                    digits=2)

分析に必要な変数

\(Y_i\) 結果を表す 21歳の誕生日付近の死亡者数（10万人あたり） all
\(D_i\) 原因を表す 法律で飲酒が認められている (\(D_i=1\)) over_21
\(R_i\) 割当変数 年齢：21歳を基準に中心化 age
  • 変数の型を確認する
str(df10)
tibble [48 × 5] (S3: tbl_df/tbl/data.frame)
 $ agecell: num [1:48] 19.1 19.2 19.2 19.3 19.4 ...
  ..- attr(*, "label")= chr "Age Cell"
  ..- attr(*, "format.stata")= chr "%9.0g"
 $ age    : num [1:48] -1.93 -1.85 -1.77 -1.68 -1.6 ...
  ..- attr(*, "label")= chr "Age Cell"
  ..- attr(*, "format.stata")= chr "%9.0g"
 $ age_sq : num [1:48] 3.73 3.42 3.12 2.84 2.57 ...
  ..- attr(*, "label")= chr "Age Cell"
  ..- attr(*, "format.stata")= chr "%9.0g"
 $ over_21: Factor w/ 2 levels "0","1": 1 1 1 1 1 1 1 1 1 1 ...
 $ all    : num [1:48] 92.8 95.1 92.1 88.4 88.7 ...
  ..- attr(*, "label")= chr "All"
  ..- attr(*, "format.stata")= chr "%9.0g"
  • 変数 over_21 の型が factor なので数値型 (numeric) に変換する
df10$over_21 <- as.numeric(df10$over_21)
  • ①から④の回帰式を lm推定する formula としてまとめて MODELS_12 と名前を付ける
MODELS_12 <- alist(
  model1 = all ~ over_21 + age,
  model2 = all ~ over_21 * age,
  model3 = all ~ over_21 + age + age_sq,
  model4 = all ~ over_21 * age + over_21 * age_sq) |> 
  enframe(name = "modelname", value = "formula")
  • MODELS_12 の中身を確認
MODELS_12
modelnameformula
model1all ~ over_21 + age
model2all ~ over_21 * age
model3all ~ over_21 + age + age_sq
model4all ~ over_21 * age + over_21 * age_sq
  • 4 つの回帰式を、一挙に推定
  • 推定結果を rdd_12 として保存
rdd_12 <- MODELS_12 |> 
  mutate(model = map(.x = formula, .f = lm, 
                     data = df10)) |> 
  mutate(pred = map(.x = model, .f = predict),   # 予測値を計算する
         result = map(.x = model, .f = tidy))
  • 結果は、rdd_12 に保存されているので、結果を1つずつ確認する

12.1 model1 の推定結果

  • 推定結果を図示してみる
p1_rdd_12 <- df10 |> 
  mutate(pred = rdd_12$pred[[1]]) |> 
  ggplot(aes(x = agecell, color = as.factor(over_21))) +
  geom_point(aes(y = all)) +
  geom_line(aes(y = pred)) +
  geom_vline(xintercept = 21, color = "gray") +
  ylim(80, 115) +
  labs(x = "年齢", y = "10万人あたりの死者数（すべての死因）") +
  theme(legend.position = "none")  +
  scale_color_manual(values = c("blue", "red")) 

plot(p1_rdd_12)

  • Angrist and Pischke (2009)の図4.2 (p.150) と同じ図が確認できる
  • 21歳の誕生日に、死者数のジャンプがある様子がわかる
  • この分析は、処置変数と割当変数のみを説明変数として利用する、最も単純な パラメトリック RD
  • そのため、予測値は「直線」
  • 21歳未満と21歳以上の両者について直線の傾きが同じ
  • このとき、処置効果の推定値と標準誤差は次のように計算できる
rdd_12 |> 
  unnest(cols = result) |> 
  filter(modelname == "model1", term == "over_21") |> 
  select(term, estimate, std.error, p.value) |> 
  knitr::kable(digits = 2)
term estimate std.error p.value
over_21 7.66 1.44 0

12.2 model2 の推定結果

  • 処置変数と割当変数の交差項をもつ model2 の結果を図示してみる
p2_rdd_12 <- df10 |> 
  mutate(pred = rdd_12$pred[[2]]) |> 
  ggplot(aes(x = agecell, color = as.factor(over_21))) +
  geom_point(aes(y = all)) +
  geom_line(aes(y = pred)) +
  geom_vline(xintercept = 21, color = "gray") +
  ylim(80, 115) +
  labs(x = "年齢", y = "10万人あたりの死者数（すべての死因）") +
  theme(legend.position = "none")  +
  scale_color_manual(values = c("blue", "red")) 

plot(p2_rdd_12)

  • 交差項を使った結果、21歳未満と21歳以上の直線の傾きが異なることがわかる
  • この関数形でもやはり、21歳の誕生日に死者数のジャンプがある様子が確認できる　　
  • 処置効果の推定値と標準誤差を計算してみる
rdd_12 |> 
  unnest(cols = result) |> 
  filter(modelname == "model2", term == "over_21") |> 
  select(term, estimate, std.error, p.value) |> 
  knitr::kable(digits = 2)
term estimate std.error p.value
over_21 7.66 1.32 0

12.3 model3 の推定結果

  • model1 に割当変数の 2 乗項を加えた model3 の推定結果を図示する
p3_rdd_12 <- df10 |> 
  mutate(pred = rdd_12$pred[[3]]) |> 
  ggplot(aes(x = agecell, color = as.factor(over_21))) +
  geom_point(aes(y = all)) +
  geom_line(aes(y = pred)) +
  geom_vline(xintercept = 21, color = "gray") +
  ylim(80, 115) +
  labs(x = "年齢", y = "10万人あたりの死者数（すべての死因）") +
  theme(legend.position = "none") +
  scale_color_manual(values = c("blue", "red")) 

plot(p3_rdd_12)

  • 2 乗項を使った結果、予測値が直線ではなく「曲線」として推定されている
  • 交差項がないので、曲線は21歳未満と21歳以上で同じ →　ジャンプをなくせば滑らかに繋がる
  • この関数形でもやはり、21歳の誕生日に死者数のジャンプがある様子が確認できる
  • 処置効果の推定値と標準誤差を計算してみる
rdd_12 |> 
  unnest(cols = result) |> 
  filter(modelname == "model3", term == "over_21") |> 
  select(term, estimate, std.error, p.value) |> 
  knitr::kable(digits = 2)
term estimate std.error p.value
over_21 7.66 1.34 0

12.4 model4 の推定結果

  • 交差項と 2 乗項を含む model4 の推定結果を可視化する
  • 比較のため、model1 で推定した直線も点線として表示する
p4_rdd_12 <- df10 |> 
  mutate(pred0 = rdd_12$pred[[1]],
         pred = rdd_12$pred[[4]]) |> 
  ggplot(aes(x = agecell, color = as.factor(over_21))) +
  geom_point(aes(y = all)) +
  geom_line(aes(y = pred0), linetype = "dashed") +
  geom_line(aes(y = pred)) +
  geom_vline(xintercept = 21, color = "gray") +
  ylim(80, 115) +
  labs(x = "年齢", y = "10万人あたりの死者数（すべての死因）") +
  theme(legend.position = "none") +
  scale_color_manual(values = c("blue", "red")) 

plot(p4_rdd_12)

  • Angrist and Pischke (2009)の図4.4 (p.158) と同じ図ができた
  • 交差項と 2 乗項の両方を使ったので、予測値が曲線になり、21歳未満と21歳以上では異なる曲線が描かれた
  • 21歳の誕生日に死者数のジャンプがある様子が確認できる
  • 処置効果の推定値と標準誤差を計算してみる
rdd_12 |> 
  unnest(cols = result) |> 
  filter(modelname == "model4", term == "over_21") |> 
  select(term, estimate, std.error, p.value) |> 
  knitr::kable(digits = 2)
term estimate std.error p.value
over_21 9.55 1.99 0
  • 推定した 4 つのモデルの推定結果を比べてみる
rdd_12 |> 
  unnest(cols = result) |> 
  filter(term == "over_21") |> 
  select(modelname, term, estimate, std.error, p.value) |> 
  knitr::kable(digits = 2)
modelname term estimate std.error p.value
model1 over_21 7.66 1.44 0
model2 over_21 7.66 1.32 0
model3 over_21 7.66 1.34 0
model4 over_21 9.55 1.99 0
  • model1 から model3 までの結果はほぼ同じ
  • model4 の推定値だけ大きい
    →　model1 から model3 か、model4 のうち、少なくとも一方の推定にバイアスが生じている
  • 私たちはどの関数形が正しいかを知らない
  • バイアスがどちらにあるか、あるいはどちらにもあるのかはわからない
  • パラメトリックRDの推定結果は関数形に依存する
  • 誤った関数形を使うと、誤った推定結果を得る
  • 我々は正しい式を知らないので、関数形の誤りの可能性は常に存在する
    →　経済学、経営学、心理学、政治学などの理論に基づいて妥当な関数形（回帰式）を選ぶ必要がある
  • 理論的に妥当性が高い複数の関数形によって推定した結果がどれも実質的に同じ
    →　推定結果に対する信頼性が高まる　　
  • 例えば、上の model4 の推定結果も model1 から model3 と同じ 7.66 であれば、処置効果が約7.7という推定結果の信頼性は高い
  • パラメトリックRDを行う場合には、推定結果が関数形に（あまり）依存しないということを示すことが重要
  • 関数形によって推定結果が大きく異なるようでは、（正しい関数形を知らない限り）推定結果を信用することはできない

13. ノンパラメトリック RD（法定飲酒年齢）

  • データの中身を確認する
stargazer(as.data.frame(df10), type = "html")
Statistic N Mean St. Dev. Min Max
agecell 48 21.000 1.151 19.068 22.932
age 48 0.000 1.151 -1.932 1.932
age_sq 48 1.296 1.171 0.002 3.731
over_21 48 1.500 0.505 1 2
all 48 95.673 3.831 88.428 105.268
  • ここで使う変数は観察可能な 2 つの変数 \(R_i, Y_i\)　
\(Y\) 結果を表す 21歳の誕生日付近の死亡者数（10万人あたり）
\(R\) 割当変数 年齢：21歳を基準に中心化

13.1 散布図と RD プロットを描く

  • 最初に、割り付け変数 \(R\) を \(x\) 軸、結果変数 \(Y\) を \(y\) 軸とした散布図と RDプロットを描く
df10 |> 
  ggplot(aes(agecell, all)) +
  geom_point() +
  stat_smooth()  # 線形ではなく LOWESS 曲線を描く

  • データをスムーズ化し、\(Y\) にジャンプがあるかどうか確認しやすい RDプロットを描いてみる
library(rdrobust)
rdplot(df10$all, df10$agecell, c = 21)

  • x = 21 のカットオフでのジャンプが確認できる

13.2 rdrobust パッケージを使った分析

  • RDD を実行する前に必要な仮定の確認は終わっているので、ここでは Step 6-2 を実行する
Step 1: 処置群と統制群を分けるルールはあるか？ 11.1 で確認済み
Step 2: RDDデザインは sharp か fuzzy か？ 11.2 で確認済み
Step 3: カットオフ付近での結果変数 \(Y\) の非連続性 11.3 で確認済み
Step 4: カットオフ付近での割当変数 \(R\) の非連続性 11.4 で確認済み
Step 5: 割当変数 \(R\) の値を自分では選べない 11.5 で確認済み
Step 6-1: 処置効果の推定（パラメトリック） 12 で実行済み
Step 6-2: 処置効果の推定（ノンパラメトリック）
Step 7: 推定結果を比較する
  • ここでは次の複数のバンド幅をモデリングして推定する
    ・IK バンド幅
    ・CER 最適バンド幅
    ・MSE 最適化バンド幅
  • まず、Imbens and Kalyanarman (2012) が推奨するバンド幅を確認してみる
IK バンド幅
IKband3 <- rdbwselect_2014(df10$all, df10$agecell, 
                           c = 21, 
                           bwselect = "IK")

IKband3$bws
            h         b
[1,] 1.645225 0.9166231
  • Imbens and Kalyanarman (2012) は 1.645225 が最適なバンド幅だと推奨している
  • では、この推奨されたバンド幅を使って、補講授業の処置効果を推定してみる
model_IK3 <- rdrobust(df10$all, df10$agecell, 
                      c = 21, 
                      h = IKband3$bws[1,1], 
                      all = TRUE)
summary(model_IK3)
Call: rdrobust

Number of Obs.                   48
BW type                      Manual
Kernel                   Triangular
VCE method                       NN

Number of Obs.                   24           24
Eff. Number of Obs.              20           20
Order est. (p)                    1            1
Order bias  (q)                   2            2
BW est. (h)                   1.645        1.645
BW bias (b)                   1.645        1.645
rho (h/b)                     1.000        1.000
Unique Obs.                      24           24

=============================================================================
        Method     Coef. Std. Err.         z     P>|z|      [ 95% C.I. ]       
=============================================================================
  Conventional     9.022     1.743     5.177     0.000     [5.606 , 12.438]    
Bias-Corrected    10.255     1.743     5.885     0.000     [6.840 , 13.671]    
        Robust    10.255     2.757     3.719     0.000     [4.851 , 15.659]    
=============================================================================
  • IK によって推奨されたバンド幅 (1.645225) の時の推定値は 9.022
CER 最適バンド幅
model_cer3 <- rdrobust(df10$all, df10$agecell, 
                       c = 21, 
                       bwselect = "cerrd", 
                       all = TRUE)

summary(model_cer3)
Call: rdrobust

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.181     4.107     2.479     0.013     [2.130 , 18.231]    
Bias-Corrected    10.247     4.107     2.495     0.013     [2.197 , 18.297]    
        Robust    10.247     4.721     2.170     0.030     [0.994 , 19.500]    
=============================================================================
MSE 最適化バンド幅
model_mse3 <- rdrobust(df10$all, df10$agecell, 
                       c = 21, 
                       bwselect = "mserd", 
                       all = TRUE)

summary(model_mse3)
Call: rdrobust

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]    
Bias-Corrected     9.689     3.591     2.698     0.007     [2.650 , 16.727]    
        Robust     9.689     4.394     2.205     0.027     [1.077 , 18.300]    
=============================================================================

13.3 分析結果のまとめ

  • 「法定飲酒年齢」に関して、ここまで実施した分析結果をまとめてみる

RDD の推定結果

手法 推定値 標準誤差 p-value 95%信頼区間
パラメトリックmodel1 7.66 1.44 0
パラメトリックmodel2 7.66 1.32 0
パラメトリックmodel3 7.66 1.34 0
パラメトリックmodel4 9.55 1.99 0
ノンパラIK バンド幅 9.022 1.743 0.000 5.606,12.438
ノンパラCER 10.181 4.107 0.013 2.130, 18.231
ノンパラCERバイアス修正 10.247 4.107 0.013 2.197, 18.297
ノンパラMSE 9.595 3.591 0.008 2.557, 16.633
ノンパラMSEバイアス修正 9.689 3.591 0.007 2.650, 16.727
ノンパラMSEロバスト 9.595 4.394 0.027 1.077, 18.300
  • パラメトリック RD もノンパラメトリック RDも得られた推定値は 10 近くで統計的にも有意
法定飲酒年齢の処置効果の推定結果 法定飲酒年齢を21歳に引き下げると、21歳の誕生日付近の死亡者数が（10万人あたり）約10人増える

V. Exercise

  • Angrist and Pischke (2009) によれば「飲酒可能な年齢になると死者数が急に増えるのは、飲酒運転による事故が増えるため」であると思われる
  • これを踏まえて、次の問題にこたえなさい
  • 下図右側の TABLE 4.1 は Angrist and Pischke (2009)からの引用である

  • 上図左側の表はこのセッションで得られたパラメトリック RDD 回帰の結果である

  • 上図左側の表にある model1 の推定値 7.66は右側の Age19-22 (1) の結果を示している

  • 上図左側の表にある model4 の推定値 9.55は右側の Age19-22 (2) の結果を示している

  • Angrist and Pischke (2009) の分析で使っているオリジナルデータ AEJfigs.dta に含まれている「自動車事故による死」(mva)と「自殺」(suicide) という変数を使って次の問題にこたえなさい

Q1:「自動車事故による死」(mva) を使ってパラメトリック RDD 回帰を実行し、TABLE 4.1 と同じ推定値 (4.53 & 4.66) が再現できることを確認しなさい　　

Q2:「自動車事故による死」(mva) を使ってノンパラメトリック RDD 回帰を実行し、Q1 で得られた結果と比較して、結論を述べなさい

Q3: 「自殺」(suicide) を使ってパラメトリック RDD 回帰を実行し、TABLE 4.1 と同じ推定値 (1.79 & 1.81) が再現できることを確認しなさい　

Q4: 「自殺」(suicide)を使ってノンパラメトリック RDD 回帰を実行し、Q3 で得られた結果と比較して、結論を述べなさい

Q5: 以上のことから、法定飲酒年齢と死者数の関係の結論としてどのようなことが導き出せるか？

  • 分析に必要なデータ（Stata形式）は次の URL からダウンロードし、haven::read_dta() を使って R で読み込むこと
download.file(url = "http://masteringmetrics.com/wp-content/uploads/2015/01/AEJfigs.dta",
              dest = "data/AEJfigs.dta")
Reference
  • Angrist, Joshua D., and Jörn-Steffen Pischke. 2009. Mostly Harmless Econometrics: An Empiricist’s Companion. Princeton: Princeton University Press.
  • Guido Imbens and Karthik Kalyanaraman (2012), Optimal Bandwidth Choice for the Regression Discontinuity Estimator Get access Arrow, The Review of Economic Studies, Volume 79.
  • Andrew Gelman & Guido Imbens (2019), Why High-Order Polynomials Should Not Be Used in Regression Discontinuity Designs, Journal of Business & Economic Statistics, Volume 37
  • PMAP 8521 • (10) Regression discontinuity I: (2) Drawing lines and measuring gaps
  • Philipp Leppert, R Tutorial: Regression Discontinuity Design
  • 矢内勇生「KUT計量経済学応用」
  • 高橋将宣『統計的因果推論の理論と実装』、共立出版、2022年
  • 清水昌平『統計的因果探索』、講談社、2017年
  • 宋財泫 (Jaehyun Song)・矢内勇生 (Yuki Yanai)「私たちのR: ベストプラクティスの探究」
  • 浅野正彦・ 矢内勇生『Rによる計量政治学』オーム社、2018年
  • 松林哲也 『政治学と因果推論』、2021年
  • Winston Chang, R Graphics Cookbook, O’Reilly Media, 2012.
  • Kosuke Imai, Quantitative Social Science: An Introduction, Princeton University Press, 2017