回帰不連続デザイン

スライド

セットアップ

実習に必要なパッケージとデータセットを読み込む。

pacman::p_load(tidyverse, rdd, rdrobust, rddensity, 
               summarytools, BalanceR)

df <- read_csv("data/rdd_data.csv")

df

# A tibble: 1,266 × 7
    year prefname outcome      rv total_cand en_cand total_votes
   <dbl> <chr>      <dbl>   <dbl>      <dbl>   <dbl>       <dbl>
 1  2000 Hokkaido    34.6 -10.8            4    2.75      261382
 2  2003 Hokkaido    40.8 -10.4            3    2.29      259740
 3  2005 Hokkaido    36.8  -2.45           4    2.60      313909
 4  2009 Hokkaido    31.1  -8.72           4    2.29      337445
 5  2012 Hokkaido    39.6   1.09           5    4.17      276894
 6  2000 Hokkaido    35.3   2.52           6    4.19      244655
 7  2003 Hokkaido    44.5  -5.13           5    2.86      236430
 8  2005 Hokkaido    30.8  -0.407          3    2.42      285519
 9  2009 Hokkaido    35.0 -11.7            5    2.50      304810
10  2012 Hokkaido    38.8   5.87           5    4.17      239023
# ℹ 1,256 more rows

分析に入る前に記述統計を確認する。

df |>
  descr(stats = c("mean", "sd", "min", "max", "n.valid"),
        transpose = TRUE, order = "p")

Descriptive Statistics  
df  
N: 1266  

                         Mean    Std.Dev         Min         Max   N.Valid
----------------- ----------- ---------- ----------- ----------- ---------
             year     2006.18       4.20     2000.00     2012.00   1266.00
          outcome       47.02      11.65       13.26       84.50   1266.00
               rv        3.58      10.20      -23.76       36.31   1266.00
       total_cand        3.76       0.95        2.00        9.00   1266.00
          en_cand        2.55       0.59        1.37        5.35   1266.00
      total_votes   214025.95   42366.82   104398.00   339780.00   1266.00

処置効果の推定

　因果効果の推定は{rdd}パッケージのRDestimate()関数、あるいは{rdrobust}パッケージのrdrobust()を使う。機能面では{rdrobust}の方が優れているものの、パッケージの使いやすさとしては{rdd}の方が優れている。本講義では頑健な推定方法については紹介しなかったものの、近年は{rdrobust}がより使われているため、ここでも{rdrobust}を使用する。いずれのパッケージも自動的に最適バンド幅を設定し¹、交差項が含まれた局所線形回帰分析を行った結果を返してくれる²。また、デフォルトのカーネルは三角 (triangular)カーネル関数だ。

rdd_fit1 <- rdrobust(y = df$outcome, x = df$rv, c = 0)

summary(rdd_fit1)

Sharp RD estimates using local polynomial regression.

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

Number of Obs.                  479          787
Eff. Number of Obs.             318          381
Order est. (p)                    1            1
Order bias  (q)                   2            2
BW est. (h)                   8.064        8.064
BW bias (b)                  12.613       12.613
rho (h/b)                     0.639        0.639
Unique Obs.                     479          787

=============================================================================
        Method     Coef. Std. Err.         z     P>|z|      [ 95% C.I. ]       
=============================================================================
  Conventional     0.725     1.580     0.459     0.646    [-2.371 , 3.820]     
        Robust         -         -     0.528     0.598    [-2.692 , 4.677]     
=============================================================================

　最適バンド幅は8.064であり、処置効果は約0.725である。これは自民党候補者の投票率から非自民候補者の最高得票率を引いた値（rv）が-8.064から8.064までのデータを使うことを意味する。そして、これらのデータに対して交差項が含まれる線形回帰分析を行うことになる。また、閾値周辺に重みを付けるために三角カーネル関数による重み付けを行った。

　結果として現職は新人に比べ、約72.5%ポイント得票率が高いという結果が得られたが、標準誤差はかなり大きく、必ずしも現職が新人より得票するとは言えないだろう（$p$ = 0.646）。今回の推定結果から日本における現職効果について、統計的に有意な効果は確認できない。

頑健性の確認

　RDDで（局所）処置効果を推定する際、分析する側はバンド幅、カーネル関数、モデル（一次関数か、二次関数かなど）を決める必要がある。これらは恣意的なものであるため、これらを少し変更しても推定値が安定しているか、つまりどれほど頑健かを確認する必要がある。

バンド幅

　rdrobust()の場合、基本的には最適バンド幅を使うことになるが、h引数を使って任意のバンド幅を指定することもできる。たとえば、既に得られた最適バンド幅8.064を使って推定してみよう。

rdd_bw1 <- rdrobust(y = df$outcome, x = df$rv, c = 0, h = 8.064)

summary(rdd_bw1)

Sharp RD estimates using local polynomial regression.

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

Number of Obs.                  479          787
Eff. Number of Obs.             318          381
Order est. (p)                    1            1
Order bias  (q)                   2            2
BW est. (h)                   8.064        8.064
BW bias (b)                   8.064        8.064
rho (h/b)                     1.000        1.000
Unique Obs.                     479          787

=============================================================================
        Method     Coef. Std. Err.         z     P>|z|      [ 95% C.I. ]       
=============================================================================
  Conventional     0.725     1.580     0.459     0.646    [-2.371 , 3.820]     
        Robust         -         -     0.013     0.990    [-4.461 , 4.521]     
=============================================================================

　先ほどと同じ結果が得られている（Robust行はこの講義では無視する。Robust推定値についてはCalonico et al. (2015)を参照されたい³。）。頑健性を報告する際は最適バンド幅における処置効果に加え、最適バンド幅を半分にした場合、2倍にした場合の結果も報告するケースが多い。それではhを8.064の半分、2倍にしたモデルも推定してみよう。

rdd_bw2 <- rdrobust(y = df$outcome, x = df$rv, c = 0, h = 8.064 / 2)

summary(rdd_bw2)

Sharp RD estimates using local polynomial regression.

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

Number of Obs.                  479          787
Eff. Number of Obs.             188          185
Order est. (p)                    1            1
Order bias  (q)                   2            2
BW est. (h)                   4.032        4.032
BW bias (b)                   4.032        4.032
rho (h/b)                     1.000        1.000
Unique Obs.                     479          787

=============================================================================
        Method     Coef. Std. Err.         z     P>|z|      [ 95% C.I. ]       
=============================================================================
  Conventional     0.162     2.205     0.074     0.941    [-4.160 , 4.485]     
        Robust         -         -     0.316     0.752    [-5.326 , 7.375]     
=============================================================================

rdd_bw3 <- rdrobust(y = df$outcome, x = df$rv, c = 0, h = 8.064 * 2)

summary(rdd_bw3)

Sharp RD estimates using local polynomial regression.

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

Number of Obs.                  479          787
Eff. Number of Obs.             453          650
Order est. (p)                    1            1
Order bias  (q)                   2            2
BW est. (h)                  16.128       16.128
BW bias (b)                  16.128       16.128
rho (h/b)                     1.000        1.000
Unique Obs.                     479          787

=============================================================================
        Method     Coef. Std. Err.         z     P>|z|      [ 95% C.I. ]       
=============================================================================
  Conventional    -0.006     1.141    -0.005     0.996    [-2.242 , 2.229]     
        Robust         -         -     0.580     0.562    [-2.309 , 4.250]     
=============================================================================

　いずれも統計的に有意な処置効果は得られない。これらの結果をまとめると以下のよるになる。

bw_compare <- tibble(Bandwidth = c("Half", "Optimal", "Double"),
                     LATE      = c(0.162, 0.725, -0.006),
                     lower     = c(-4.160, -2.371, -2.242),
                     upper     = c(4.485, 3.820, 2.2249))

bw_compare

# A tibble: 3 × 4
  Bandwidth   LATE lower upper
  <chr>      <dbl> <dbl> <dbl>
1 Half       0.162 -4.16  4.49
2 Optimal    0.725 -2.37  3.82
3 Double    -0.006 -2.24  2.22

　これらをpoint-rangeプロットで可視化してみよう。

bw_compare |>
  mutate(Bandwidth = fct_inorder(Bandwidth)) |>
  ggplot() +
  geom_hline(yintercept = 0) +
  geom_pointrange(aes(x = Bandwidth, y = LATE, 
                      ymin = lower, ymax = upper)) +
  theme_bw(base_size = 12)

カーネル

　カーネル関数はkernel引数で指定することができる。指定しない場合、既定値としてkernel = "triangular"になるが、他にも"uniform"と"epanechnikov"がある。

rdd_kernel1 <- rdrobust(y = df$outcome, x = df$rv, c = 0, 
                        kernel = "triangular")
rdd_kernel2 <- rdrobust(y = df$outcome, x = df$rv, c = 0, 
                        kernel = "uniform")
rdd_kernel3 <- rdrobust(y = df$outcome, x = df$rv, c = 0, 
                        kernel = "epanechnikov")

summary(rdd_kernel1)

Sharp RD estimates using local polynomial regression.

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

Number of Obs.                  479          787
Eff. Number of Obs.             318          381
Order est. (p)                    1            1
Order bias  (q)                   2            2
BW est. (h)                   8.064        8.064
BW bias (b)                  12.613       12.613
rho (h/b)                     0.639        0.639
Unique Obs.                     479          787

=============================================================================
        Method     Coef. Std. Err.         z     P>|z|      [ 95% C.I. ]       
=============================================================================
  Conventional     0.725     1.580     0.459     0.646    [-2.371 , 3.820]     
        Robust         -         -     0.528     0.598    [-2.692 , 4.677]     
=============================================================================

summary(rdd_kernel2)

Sharp RD estimates using local polynomial regression.

Number of Obs.                 1266
BW type                       mserd
Kernel                      Uniform
VCE method                       NN

Number of Obs.                  479          787
Eff. Number of Obs.             245          264
Order est. (p)                    1            1
Order bias  (q)                   2            2
BW est. (h)                   5.651        5.651
BW bias (b)                  10.115       10.115
rho (h/b)                     0.559        0.559
Unique Obs.                     479          787

=============================================================================
        Method     Coef. Std. Err.         z     P>|z|      [ 95% C.I. ]       
=============================================================================
  Conventional     0.611     1.699     0.360     0.719    [-2.718 , 3.941]     
        Robust         -         -     0.516     0.606    [-2.878 , 4.933]     
=============================================================================

summary(rdd_kernel3)

Sharp RD estimates using local polynomial regression.

Number of Obs.                 1266
BW type                       mserd
Kernel                   Epanechnikov
VCE method                       NN

Number of Obs.                  479          787
Eff. Number of Obs.             298          349
Order est. (p)                    1            1
Order bias  (q)                   2            2
BW est. (h)                   7.420        7.420
BW bias (b)                  12.211       12.211
rho (h/b)                     0.608        0.608
Unique Obs.                     479          787

=============================================================================
        Method     Coef. Std. Err.         z     P>|z|      [ 95% C.I. ]       
=============================================================================
  Conventional     0.722     1.586     0.455     0.649    [-2.387 , 3.832]     
        Robust         -         -     0.560     0.575    [-2.633 , 4.740]     
=============================================================================

　以上の結果をまとめたものが以下である。

kernel_compare <- tibble(Kernel = c("Triangular", 
                                    "Unifrom", 
                                    "Epanechnikov"),
                         LATE   = c(0.725, 0.611, 0.722),
                         lower  = c(-2.371, -2.718, -2.387),
                         upper  = c(3.820, 3.941, 3.832))

kernel_compare

# A tibble: 3 × 4
  Kernel        LATE lower upper
  <chr>        <dbl> <dbl> <dbl>
1 Triangular   0.725 -2.37  3.82
2 Unifrom      0.611 -2.72  3.94
3 Epanechnikov 0.722 -2.39  3.83

kernel_compare |>
  mutate(Kernel = fct_inorder(Kernel)) |>
  ggplot() +
  geom_hline(yintercept = 0) +
  geom_pointrange(aes(x = Kernel, y = LATE, 
                      ymin = lower, ymax = upper)) +
  theme_bw(base_size = 12)

関数

　パラメトリック推定、またはセミパラメトリック推定の場合、応答変数と強制変数間の関係をある関数として仮定する必要がある。セミパラメトリック推定は関数設定の影響力が比較的小さいが、バンド幅内のケース数が少なくなると、関数の形に影響を受けやすい。ここでは1次関数から4次関数までモデルを変えながら、推定値が安定しているかを確認してみよう。rdrobust()で関数の次数を指定するためにはp引数を使用する。既定値は1であるが、p次関数を使う場合はpを指定する必要がある。

rdd_p1 <- rdrobust(y = df$outcome, x = df$rv, c = 0, p = 1)
rdd_p2 <- rdrobust(y = df$outcome, x = df$rv, c = 0, p = 2)
rdd_p3 <- rdrobust(y = df$outcome, x = df$rv, c = 0, p = 3)
rdd_p4 <- rdrobust(y = df$outcome, x = df$rv, c = 0, p = 4)

summary(rdd_p1)

Sharp RD estimates using local polynomial regression.

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

Number of Obs.                  479          787
Eff. Number of Obs.             318          381
Order est. (p)                    1            1
Order bias  (q)                   2            2
BW est. (h)                   8.064        8.064
BW bias (b)                  12.613       12.613
rho (h/b)                     0.639        0.639
Unique Obs.                     479          787

=============================================================================
        Method     Coef. Std. Err.         z     P>|z|      [ 95% C.I. ]       
=============================================================================
  Conventional     0.725     1.580     0.459     0.646    [-2.371 , 3.820]     
        Robust         -         -     0.528     0.598    [-2.692 , 4.677]     
=============================================================================

summary(rdd_p2)

Sharp RD estimates using local polynomial regression.

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

Number of Obs.                  479          787
Eff. Number of Obs.             358          433
Order est. (p)                    2            2
Order bias  (q)                   3            3
BW est. (h)                   9.205        9.205
BW bias (b)                  12.622       12.622
rho (h/b)                     0.729        0.729
Unique Obs.                     479          787

=============================================================================
        Method     Coef. Std. Err.         z     P>|z|      [ 95% C.I. ]       
=============================================================================
  Conventional     0.260     2.139     0.122     0.903    [-3.931 , 4.452]     
        Robust         -         -    -0.030     0.976    [-4.804 , 4.661]     
=============================================================================

summary(rdd_p3)

Sharp RD estimates using local polynomial regression.

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

Number of Obs.                  479          787
Eff. Number of Obs.             379          476
Order est. (p)                    3            3
Order bias  (q)                   4            4
BW est. (h)                  10.311       10.311
BW bias (b)                  13.458       13.458
rho (h/b)                     0.766        0.766
Unique Obs.                     479          787

=============================================================================
        Method     Coef. Std. Err.         z     P>|z|      [ 95% C.I. ]       
=============================================================================
  Conventional    -0.389     2.676    -0.145     0.884    [-5.633 , 4.855]     
        Robust         -         -    -0.240     0.811    [-6.462 , 5.055]     
=============================================================================

summary(rdd_p4)

Sharp RD estimates using local polynomial regression.

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

Number of Obs.                  479          787
Eff. Number of Obs.             420          567
Order est. (p)                    4            4
Order bias  (q)                   5            5
BW est. (h)                  12.653       12.653
BW bias (b)                  15.671       15.671
rho (h/b)                     0.807        0.807
Unique Obs.                     479          787

=============================================================================
        Method     Coef. Std. Err.         z     P>|z|      [ 95% C.I. ]       
=============================================================================
  Conventional    -0.472     3.002    -0.157     0.875    [-6.356 , 5.411]     
        Robust         -         -    -0.124     0.902    [-6.776 , 5.972]     
=============================================================================

　以上の結果をまとめたものが以下である。

order_compare <- tibble(Order = 1:4,
                        LATE  = c(0.725, 0.260, -0.389, -0.472),
                        lower = c(-2.371, -3.931, -5.633, -6.356),
                        upper = c(3.820, 4.452, 4.855, 5.411))

order_compare

# A tibble: 4 × 4
  Order   LATE lower upper
  <int>  <dbl> <dbl> <dbl>
1     1  0.725 -2.37  3.82
2     2  0.26  -3.93  4.45
3     3 -0.389 -5.63  4.86
4     4 -0.472 -6.36  5.41

order_compare |>
  ggplot() +
  geom_hline(yintercept = 0) +
  geom_pointrange(aes(x = Order, y = LATE, 
                      ymin = lower, ymax = upper)) +
  labs(x = "Order of Local Polynomial Regression") +
  theme_bw(base_size = 12)

　結果は大きく変わらず、安定していることが分かる。

可視化

　{rdrobust}パッケージには可視化に便利なrdplot()関数が用意されている。使い方はrdrobust()とほぼ同じで、yには応答変数を、xには強制変数、cには閾値を指定すれば良い。

rdplot(y = df$outcome, x = df$rv, c = 0)

　観測値が少ないように見えるが、これは強制変数を区間に分け、区間内の平均値を示したものである。通常のRDDはサンプルサイズが大きいため、散布図+回帰直線（曲線）だと線が見えなかったり、傾向が見にくくなる傾向がある。サンプルサイズが数百程度なら全観測値を見せても良いだろうが、今回は1200以上であり、このような見せ方が効果的である。

　1つ注意すべき点は、表示される回帰直線（曲線）の場合、2次関数が使用される。1次関数にフィットさせるためには、rdrobust()同様、p = 1を指定すれば良い。

rdplot(y = df$outcome, x = df$rv, c = 0, p = 1)

　また、この回帰直線の場合、カーネル関数は矩形関数である。三角形関数にするためには、更にkernel = "triangular"を指定する。また、x.lab、y.lab、title引数でラベルを修正することもできる。

rdplot(y = df$outcome, x = df$rv, c = 0, p = 1,
       kernel = "triangular",
       x.label = "Vote Margin in Election t",
       y.label = "Vote Share in Election t+1",
       title = "")

　以上の例は、バンド幅を設定せず、全観測値を利用したものである。rdrobust()のようにバンド幅を指定することはできないため、subset引数を使って使用するデータを制限することができる。たとえば、rvが-15より大きく、15より小さいケースのみを使う場合、subset = (df$rv > -15 & df$rv < 15)と指定すれば良い。

rdplot(y = df$outcome, x = df$rv, c = 0, p = 1,
       kernel = "triangular", 
       subset = (df$rv > -15 & df$rv < 15),
       x.label = "Vote Margin in Election t",
       y.label = "Vote Share in Election t+1",
       title = "")

仮定の確認

交絡要因の連続性

　RDDの重要な仮定の一つとして、交絡要因の連続性がある。交絡要因として考えられる要因が、処置群に割り当てられることでジャンプした場合、観察される処置効果がが処置によるものか、交絡要因のジャンブによるものかが識別できないからだ。今回の例では処置効果が見られていないが、それでもこの仮定は確認する価値がある。処置による効果（$X \rightarrow Y$）と交絡要因による効果（$Z \rightarrow Y$）が両方存在するケースを考えてみよう。もしこの2つの効果の符号が逆である場合、処置効果（$X \rightarrow Y$）が交絡要因による効果（$Z \rightarrow Y$）に相殺される可能性もあるからだ。

　確認する方法は簡単だ。もう一度RDDをするだけだ。ただし、応答変数が得票率（outcome）でなく、交絡要因に代わるだけだ。今回は候補者数（total_cand）、有効候補者数（en_cand）、得票数（total_votes）に対してRDDを行ってみよう。

assumption_fit1 <- rdrobust(y = df$total_cand, x = df$rv)
assumption_fit2 <- rdrobust(y = df$en_cand, x = df$rv)
assumption_fit3 <- rdrobust(y = df$total_votes, x = df$rv)

summary(assumption_fit1)

Sharp RD estimates using local polynomial regression.

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

Number of Obs.                  479          787
Eff. Number of Obs.             243          261
Order est. (p)                    1            1
Order bias  (q)                   2            2
BW est. (h)                   5.547        5.547
BW bias (b)                   8.621        8.621
rho (h/b)                     0.643        0.643
Unique Obs.                     479          787

=============================================================================
        Method     Coef. Std. Err.         z     P>|z|      [ 95% C.I. ]       
=============================================================================
  Conventional    -0.066     0.209    -0.315     0.753    [-0.475 , 0.344]     
        Robust         -         -    -0.415     0.678    [-0.593 , 0.385]     
=============================================================================

summary(assumption_fit2)

Sharp RD estimates using local polynomial regression.

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

Number of Obs.                  479          787
Eff. Number of Obs.             280          314
Order est. (p)                    1            1
Order bias  (q)                   2            2
BW est. (h)                   6.629        6.629
BW bias (b)                  10.114       10.114
rho (h/b)                     0.655        0.655
Unique Obs.                     479          787

=============================================================================
        Method     Coef. Std. Err.         z     P>|z|      [ 95% C.I. ]       
=============================================================================
  Conventional     0.032     0.111     0.289     0.772    [-0.186 , 0.251]     
        Robust         -         -     0.149     0.881    [-0.242 , 0.282]     
=============================================================================

summary(assumption_fit3)

Sharp RD estimates using local polynomial regression.

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

Number of Obs.                  479          787
Eff. Number of Obs.             291          341
Order est. (p)                    1            1
Order bias  (q)                   2            2
BW est. (h)                   7.131        7.131
BW bias (b)                  11.445       11.445
rho (h/b)                     0.623        0.623
Unique Obs.                     479          787

=============================================================================
        Method     Coef. Std. Err.         z     P>|z|      [ 95% C.I. ]       
=============================================================================
  Conventional  1222.634  7020.544     0.174     0.862[-12537.380 , 14982.648] 
        Robust         -         -     0.042     0.967[-16145.541 , 16846.362] 
=============================================================================

　いずれも統計的に有意なジャンプは見られない。以上の検定結果から「仮定は満たされている」ことは主張できないものの、「仮定が満たされていないとは言えない」までは主張できるはずだ。

バランスチェック

　ノンパラメトリックRDDの場合、バンド幅内であれば、処置群と統制群の性質はほぼ同じであると仮定する。つまり、処置変数を除く共変量が処置群と統制群の間において均質であることを意味する。それでは{BalanceR}を使って、候補者数（total_cand）、有効候補者数（en_cand）、得票数（total_votes）が処置群と統制群の間に差があるかを確認してみよう。

df |>
  # 処置の有無を示す treat 変数を作成
  mutate(treat = if_else(rv > 0, "yes", "no")) |>
  BalanceR(group = treat,
           cov = total_cand:total_votes) |>
  plot(abs = TRUE) +
  scale_y_discrete(label = c("total_cand" = "Total number of candidates",
                             "en_cand" = "Effective number of candidates",
                             "total_votes" = "Total votes"))

　得票数（total_votes）の場合、標準化差分が非常に大きいことが分かる。それではバンド幅内のサンプルに限定すればどうだろうか。filter()を使ってrvが-8.064より大きく、8.064より小さいサンプルに絞っってバランスチェックをしてみよう。

df |>
  mutate(treat = if_else(rv > 0, "yes", "no")) |>
  filter(rv > -8.064 & rv < 8.064) |>
  BalanceR(group = treat,
           cov = total_cand:total_votes) |>
  plot(abs = TRUE) +
  scale_y_discrete(label = c("total_cand" = "Total number of candidates",
                             "en_cand" = "Effective number of candidates",
                             "total_votes" = "Total votes"))

　有効候補者数（en_cand）のバランスがむしろ悪くなったものの、他の2つの変数のバランスは改善されていることが分かる。

強制変数の操作可能性

　RDDのもう一つの重要な仮定として、閾値周辺において強制変数の操作が行われてはいけない。得票率の差を操作することは極めて困難なので、今回は問題はないと考えられるが、たとえばフランス地方議会選挙のように人口によって制度が変わる場合、特定の選挙制度を採用するために人口を操作することは不可能ではないだろう。

　この仮定を確認、検定する手法がMcCrayの密度検定 (density test)だ (McCray 2006)⁴。簡単に説明すると、強制変数の密度関数が閾値周辺においてジャンプしているか否かを確認する方法である。もし、操作が行われているとしたら、密度関数が断絶するだろう。

　密度検定は{rdd}のDCdensity()で簡単に行うことができる。第一引数は強制変数を、cutpointには閾値を指定する（既定値は0であるため、今回は省略可能）。

DCdensity(df$rv, cutpoint = 0)

[1] 0.3002424

　図と長さ1のnumeric型ベクトルが出力されるが⁵、図は密度分布を可視化したものであり、数値は「密度関数は連続している」という帰無仮説に対する$p$値である。これが$\alpha$（通常、$\alpha = 0.05$）を下回る場合、帰無仮説は棄却され、密度関数が断絶していると判断できる。つまり、RDDの仮定を満たしていないことを意味する。

　{rdd}のDCdensity()以外にも、密度検定専用のパッケージ{rddenstiy}のrddensity()を使うことも可能だ。検定方法は基本的に同じだが、検定の際に使用するパラメーターや標準誤差計算のアルゴリズムが異なるため、結果はやや異なる。使い方はXに強制変数を、cに閾値を指定すれば良い。他にも十数種類のパラメーターが指定できるが詳細はコンソール上で?rddensityを入力し、ヘルプを参照すること。

Density_Test <- rddensity(X = df$rv, c = 0)
summary(Density_Test)


Manipulation testing using local polynomial density estimation.

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

c = 0                 Left of c           Right of c          
Number of obs         479                 787                 
Eff. Number of obs    259                 289                 
Order est. (p)        2                   2                   
Order bias (q)        3                   3                   
BW est. (h)           6.112               6.078               

Method                T                   P > |T|             
Robust                -0.818              0.4134              


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

Window Length / 2          <c     >=c    P>|T|
0.448                      20      24    0.6516
0.896                      41      41    1.0000
1.343                      65      59    0.6536
1.791                      89      79    0.4876
2.239                     108      98    0.5307
2.687                     129     116    0.4434
3.134                     151     135    0.3751
3.582                     168     163    0.8260
4.030                     188     185    0.9175
4.478                     200     204    0.8814

　密度検定の結果（$p$値）は中間辺りにある# Robust行の0.4134だ。ここでも帰無仮説は棄却されず、強制変数の操作が行われているとは言えない。これらの結果を可視化の際はrdplotdensity()関数を使う。第一引数はrddensity()から得られたオブジェクト名を指定し、Xには強制変数を指定する。その他の引数についてはヘルプ（コンソール上で?rdplotdensity）を参照すること。

Density_Plot <- rdplotdensity(Density_Test, X = df$rv, 
                              type = "both", lwd = 1, pwd = 3, pty = 19)

　また、要約結果の下段にあるBinomial testsは密度分布に代わるもう一つの検定手法だ。1行目はケースが20個入る範囲と、その中での処置群と統制群の大きさ、そしてその差の検定である。ここでは0.215だが、これはrvが-0.215から0.215の間に20個のケースがあるということを意味する。統制群は7ケース、処置群は13ケースである。もし、強制変数の操作が行われなかったのであれば、処置群の割合は0.5になるはずである。右のP>|T|列は、$p = 0.5$を帰無仮説とした二項検定における$p$値である。もし、この値が$\alpha$を下回ると、閾値周辺において何らかの操作が行われた可能性があることを示唆する。

　2行目は1行目の幅を2倍に、3行目は1行目の幅を3倍に、…したものである。いずれも$p$値は0.05以上であり、強制変数の操作が行われたとは言えない。

脚注

もし、自分でバンド幅を指定したい場合、bw = ...の引数を加える。↩︎
閾値のデフォルトは0だ。もし、閾値が0ではない場合、cutpoint = ...の引数を設定する。↩︎
Calonico, S., M. D. Cattaneo, and R. Titiunik. 2015b. “rdrobust: An R Package for Robust Nonparametric Inference in Regression-Discontinuity Designs,” R Journal, 7(1): 38-51.↩︎
McCray, Justin. 2008. “Manipulation of the running variable in the regression discontinuity design: A density test,” Journal of Econometrics, 142(2): 698-714.↩︎
図が不要ならplot = FALSEを指定する。↩︎