差分の差分法

スライド

新しいタブで開く

セットアップ

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

pacman::p_load(tidyverse,     # Rの必須パッケージ
               summarytools,  # 記述統計
               modelsummary,  # 推定結果の要約
               estimatr,      # ロバストな回帰分析
               gsynth,        # 一般化SCM
               panelView,     # パネルデータのチェック
               CausalImpact)  # CausalImpactの実装

did_df <- read_csv("data/did_data5.csv")

did_df
# A tibble: 18,749 × 14
   county state  year shooting fatal_s…¹ non_f…² turnout demvote popul…³ non_w…⁴
    <dbl> <chr> <dbl>    <dbl>     <dbl>   <dbl>   <dbl>   <dbl>   <dbl>   <dbl>
 1   1001 01     1996        0         0       0    56.1    32.5   40207   0.190
 2   1001 01     2000        0         0       0    56.8    28.7   44021   0.189
 3   1001 01     2004        0         0       0    60.2    23.7   48366   0.196
 4   1001 01     2008        0         0       0    64.1    25.8   53277   0.205
 5   1001 01     2012        0         0       0    61.0    26.5   55027   0.212
 6   1001 01     2016        0         0       0    60.7    24.0   55416   0.228
 7   1003 01     1996        0         0       0    51.8    27.1  125412   0.123
 8   1003 01     2000        0         0       0    54.6    24.8  141342   0.122
 9   1003 01     2004        0         0       0    59.8    22.5  156266   0.122
10   1003 01     2008        0         0       0    62.4    23.8  175827   0.122
# … with 18,739 more rows, 4 more variables: change_unem_rate <dbl>,
#   county_f <dbl>, state_f <chr>, year_f <dbl>, and abbreviated variable names
#   ¹​fatal_shooting, ²​non_fatal_shooting, ³​population, ⁴​non_white
# ℹ Use `print(n = ...)` to see more rows, and `colnames()` to see all variable names

 データの詳細はスライドを参照すること。DID推定には時間(年)とカウンティー(郡)の固定効果を投入し、州レベルでクラスタリングした標準誤差を使う予定である。これらの変数を予めfactor化しておこう。factor化した変数は変数名の後ろに_fを付けて、新しい列として追加しておく。

did_df <- did_df |>
  mutate(county_f = factor(county),
         state_f  = factor(state),
         year_f   = factor(year))

did_df
# A tibble: 18,749 × 14
   county state  year shooting fatal_s…¹ non_f…² turnout demvote popul…³ non_w…⁴
    <dbl> <chr> <dbl>    <dbl>     <dbl>   <dbl>   <dbl>   <dbl>   <dbl>   <dbl>
 1   1001 01     1996        0         0       0    56.1    32.5   40207   0.190
 2   1001 01     2000        0         0       0    56.8    28.7   44021   0.189
 3   1001 01     2004        0         0       0    60.2    23.7   48366   0.196
 4   1001 01     2008        0         0       0    64.1    25.8   53277   0.205
 5   1001 01     2012        0         0       0    61.0    26.5   55027   0.212
 6   1001 01     2016        0         0       0    60.7    24.0   55416   0.228
 7   1003 01     1996        0         0       0    51.8    27.1  125412   0.123
 8   1003 01     2000        0         0       0    54.6    24.8  141342   0.122
 9   1003 01     2004        0         0       0    59.8    22.5  156266   0.122
10   1003 01     2008        0         0       0    62.4    23.8  175827   0.122
# … with 18,739 more rows, 4 more variables: change_unem_rate <dbl>,
#   county_f <fct>, state_f <fct>, year_f <fct>, and abbreviated variable names
#   ¹​fatal_shooting, ²​non_fatal_shooting, ³​population, ⁴​non_white
# ℹ Use `print(n = ...)` to see more rows, and `colnames()` to see all variable names

 連続変数(shootingからchange_unem_rateまで)の記述統計量を出力する。ここで一つ注意が必要だ。それはselect()関数の使い方である。具体的に言えば、使い方そのものは変わらない。問題は、現在、select()関数が2つ以上存在するという点だ。最初に読み込んだ{gsynth}パッケージは{MASS}パッケージにも依存しているため、自動的に{MASS}も読み込まれる。この{MASS}にもselect()関数が存在する。しかも、{gsynth}は{tidyverse}よりも先に読み込まれているので、select()のみ書くと、{MASS}のselect()関数として認識する。したがって、「どのパッケージのselect()関数か」を具体的に書く必要がある。書き方はパッケージ名::関数名()だ。{tidyverse}内には{dplyr}というパッケージがあり、列の抽出に使うselect()関数はこの{dplyr}が提供しているため、dplyr::select()と書く必要がある。

did_df |>
  dplyr::select(shooting:change_unem_rate) |>
  descr(stats = c("mean", "sd", "min", "max", "n.valid"),
        transpose = TRUE, order = "p")
Error : Can't find dplyr
Descriptive Statistics  

                               Mean     Std.Dev      Min           Max    N.Valid
------------------------ ---------- ----------- -------- ------------- ----------
                shooting       0.00        0.07     0.00          1.00   18749.00
          fatal_shooting       0.00        0.05     0.00          1.00   18749.00
      non_fatal_shooting       0.00        0.04     0.00          1.00   18749.00
                 turnout      57.94       10.32     1.08        300.00   18627.00
                 demvote      39.01       13.79     3.14         92.85   18627.00
              population   95121.69   306688.63    55.00   10137915.00   18749.00
               non_white       0.13        0.16     0.00          0.96   18749.00
        change_unem_rate      -0.39        2.34   -19.60         15.30   18749.00

Diff-in-Diff

 それでは差分の差分法の実装について紹介する。推定式は以下の通りである。

\[ \mbox{Outcome}_{i, t} = \beta_0 + \beta_1 \mbox{Shooting}_{i, t} + \sum_k \delta_{k, i, t} \mbox{Controls}_{k, i, t} + \lambda_{t} + \omega_{i} + \varepsilon_{i, t} \]

  • \(\mbox{Otucome}\): 応答変数
    • turnout: 投票率(大統領選挙)
    • demvote: 民主党候補者の得票率
  • \(\mbox{Shooting}\): 処置変数
    • shooting: 銃撃事件の発生有無
    • fatal_shooting: 死者を伴う銃撃事件の発生有無
    • non_fatal_shooting: 死者を伴わない銃撃事件の発生有無
  • \(\mbox{Controls}\): 統制変数
    • population: カウンティーの人口
    • non_white: 非白人の割合
    • change_unem_rate: 失業率の変化
    • 統制変数あり/なしのモデルを個別に推定
  • \(\lambda\): 年固定効果
  • \(\omega\): カウンティー固定効果

応答変数が2種類、処置変数が3種類、共変量の有無でモデルを分けるので、推定するモデルは計12個である。

モデル オブジェクト名 応答変数 処置変数 統制変数
モデル1 did_fit1 turnout shooting なし
モデル2 did_fit2 turnout shooting あり
モデル3 did_fit3 turnout fatal_shooting なし
モデル4 did_fit4 turnout fatal_shooting あり
モデル5 did_fit5 turnout non_fatal_shooting なし
モデル6 did_fit6 turnout non_fatal_shooting あり
モデル7 did_fit7 demvote shooting なし
モデル8 did_fit8 demvote shooting あり
モデル9 did_fit9 demvote fatal_shooting なし
モデル10 did_fit10 demvote fatal_shooting あり
モデル11 did_fit11 demvote non_fatal_shooting なし
モデル12 did_fit12 demvote non_fatal_shooting あり

 まずはモデル1を推定し、did_fit1という名のオブジェクトとして格納する。基本的には線形回帰分析であるため、lm()でも推定はできる。しかし、差分の差分法の場合、通常、クラスター化した頑健な標準誤差(cluster robust standard error)を使う。lm()単体ではこれが計算できないため、今回は{estimatr}パッケージが提供するlm_robust()関数を使用する。使い方はlm()同様、まず回帰式と使用するデータ名を指定する。続いて、固定効果をfixed_effects引数で指定する1。書き方は~固定効果変数1 + 固定効果変数2 + ...である。回帰式と違って、~の左側には変数がないことに注意すること。続いて、clusters引数でクラスタリングする変数を指定する。今回は州レベルでクラスタリングするので、state_fで良い。最後に標準誤差のタイプを指定するが、デフォルトは"CR2"となっている。今回のデータはそれなりの大きさのデータであり、"CR2"だと推定時間が非常に長くなる。ここでは推定時間が比較的早い"stata"とする。

did_fit1 <- lm_robust(turnout ~ shooting, 
                      data          = did_df, 
                      fixed_effects = ~year_f + county_f,
                      clusters      = state_f,
                      se_type       = "stata")

summary(did_fit1)

Call:
lm_robust(formula = turnout ~ shooting, data = did_df, clusters = state_f, 
    fixed_effects = ~year_f + county_f, se_type = "stata")

Standard error type:  stata 

Coefficients:
         Estimate Std. Error t value Pr(>|t|) CI Lower CI Upper DF
shooting  -0.5211     0.6457  -0.807   0.4235   -1.819   0.7764 49

Multiple R-squared:  0.8575 ,   Adjusted R-squared:  0.8288
Multiple R-squared (proj. model):  6.907e-05 ,  Adjusted R-squared (proj. model):  -0.2008 
F-statistic (proj. model): 0.6513 on 1 and 49 DF,  p-value: 0.4235

 処置効果の推定値は-0.521である。これは学校内銃撃事件が発生したカウンティーの場合、大統領選挙において投票率が約-0.521%p低下することを意味する。しかし、標準誤差がかなり大きく、統計的有意な結果ではない。つまり、「学校内銃撃事件が投票率を上げる(or 下げる)とは言えない」と解釈できる。決して「学校内銃撃事件が投票率を上げない(or 下げない)」と解釈しないこと。

 共変量を投入してみたらどうだろうか。たとえば、人口は自治体の都市化程度を表すこともあるので、都市化程度と投票率には関係があると考えられる。また、人口が多いと自然に事件が発生する確率もあがるので、交絡要因として考えられる。人種や失業率も同様であろう。ここではカウンティーの人口(population)、非白人の割合(non_white)、失業率の変化(change_unem_rate)を統制変数として投入し、did_fit2という名で格納する。

did_fit2 <- lm_robust(turnout ~ shooting + 
                        population + non_white + change_unem_rate, 
                      data          = did_df, 
                      fixed_effects = ~year_f + county_f,
                      clusters      = state_f,
                      se_type       = "stata")

summary(did_fit2)

Call:
lm_robust(formula = turnout ~ shooting + population + non_white + 
    change_unem_rate, data = did_df, clusters = state_f, fixed_effects = ~year_f + 
    county_f, se_type = "stata")

Standard error type:  stata 

Coefficients:
                   Estimate Std. Error t value Pr(>|t|)   CI Lower   CI Upper
shooting         -7.098e-01  6.237e-01  -1.138  0.26066 -1.963e+00  5.436e-01
population        8.029e-06  5.688e-06   1.411  0.16443 -3.402e-06  1.946e-05
non_white        -3.483e+01  1.558e+01  -2.236  0.02990 -6.613e+01 -3.534e+00
change_unem_rate  1.592e-01  6.584e-02   2.418  0.01938  2.688e-02  2.915e-01
                 DF
shooting         49
population       49
non_white        49
change_unem_rate 49

Multiple R-squared:  0.8598 ,   Adjusted R-squared:  0.8316
Multiple R-squared (proj. model):  0.01669 ,    Adjusted R-squared (proj. model):  -0.1811 
F-statistic (proj. model): 5.047 on 4 and 49 DF,  p-value: 0.001739

 処置効果の推定値は-0.710である。これは他の条件が同じ場合、学校内銃撃事件が発生したカウンティーは大統領選挙において投票率が約-0.710%p低下することを意味する。ちなみに、e-01\(\times 10^{-1}\)を、e-06\(\times 10^{-6}\)を、e+01\(\times 10^{1}\)意味する。今回も統計的に非有意な結果が得られている。

 これまでの処置変数は死者の有無と関係なく、学校内銃撃事件が発生したか否かだった。もしかしたら、死者を伴う銃撃事件が発生した場合、その効果が大きいかも知れない。したがって、これからは処置変数を死者を伴う学校内銃撃事件の発生有無(fatal_shooting)、死者を伴わない学校内銃撃事件の発生有無(non_fatal_shooting)に変えてもう一度推定してみよう。

did_fit3 <- lm_robust(turnout ~ fatal_shooting, 
                      data          = did_df, 
                      fixed_effects = ~year_f + county_f,
                      clusters      = state_f,
                      se_type       = "stata")

did_fit4 <- lm_robust(turnout ~ fatal_shooting + 
                        population + non_white + change_unem_rate, 
                      data          = did_df, 
                      fixed_effects = ~year_f + county_f,
                      clusters      = state_f,
                      se_type       = "stata")

did_fit5 <- lm_robust(turnout ~ non_fatal_shooting, 
                      data          = did_df, 
                      fixed_effects = ~year_f + county_f,
                      clusters      = state_f,
                      se_type       = "stata")

did_fit6 <- lm_robust(turnout ~ non_fatal_shooting + 
                        population + non_white + change_unem_rate, 
                      data          = did_df, 
                      fixed_effects = ~year_f + county_f,
                      clusters      = state_f,
                      se_type       = "stata")

 これまで推定してきた6つのモデルを比較してみよう。

modelsummary(list(did_fit1, did_fit2, did_fit3, 
                  did_fit4, did_fit5, did_fit6))
Model 1 Model 2 Model 3 Model 4 Model 5 Model 6
shooting −0.521 −0.710
(0.646) (0.624)
population 0.000 0.000 0.000
(0.000) (0.000) (0.000)
non_white −34.834 −34.844 −34.814
(15.575) (15.618) (15.614)
change_unem_rate 0.159 0.159 0.160
(0.066) (0.066) (0.066)
fatal_shooting −0.678 −0.918
(0.619) (0.658)
non_fatal_shooting −0.239 −0.327
(1.094) (1.145)
Num.Obs. 18627 18627 18627 18627 18627 18627
R2 0.857 0.860 0.857 0.860 0.857 0.860
R2 Adj. 0.829 0.832 0.829 0.832 0.829 0.832
AIC 103543.4 103237.2 103543.3 103237.1 103544.6 103239.4
BIC 103559.0 103276.4 103559.0 103276.3 103560.2 103278.6
RMSE 3.90 3.87 3.90 3.87 3.90 3.87
Std.Errors by: state_f by: state_f by: state_f by: state_f by: state_f by: state_f

 いずれのモデルも統計的に有意な処置効果は確認されていない。これらの結果を表として報告するには紙がもったいない気もする。これらの結果はOnline Appendixに回し、本文中には処置効果の点推定値と95%信頼区間を示せば良いだろう。

 {broom}のtidy()関数で推定結果のみを抽出し、それぞれオブジェクトとして格納しておこう。

tidy_fit1 <- tidy(did_fit1, conf.int = TRUE)
tidy_fit2 <- tidy(did_fit2, conf.int = TRUE)
tidy_fit3 <- tidy(did_fit3, conf.int = TRUE)
tidy_fit4 <- tidy(did_fit4, conf.int = TRUE)
tidy_fit5 <- tidy(did_fit5, conf.int = TRUE)
tidy_fit6 <- tidy(did_fit6, conf.int = TRUE)

 全て確認する必要はないので、tidy_fit1のみを確認してみる。

tidy_fit1
      term   estimate std.error  statistic   p.value conf.low conf.high df
1 shooting -0.5210841 0.6456718 -0.8070417 0.4235426 -1.81861  0.776442 49
  outcome
1 turnout

 以上の6つの表形式オブジェクトを一つの表としてまとめる。それぞれのオブジェクトには共変量の有無_処置変数の種類の名前を付けよう。共変量なしのモデルはM1、ありのモデルはM2とする。処置変数はshootingの場合はTr1fatal_shootingTr2non_fatal_shootingTr3とする。

did_est1 <- bind_rows(list("M1_Tr1" = tidy_fit1,
                           "M2_Tr1" = tidy_fit2,
                           "M1_Tr2" = tidy_fit3,
                           "M2_Tr2" = tidy_fit4,
                           "M1_Tr3" = tidy_fit5,
                           "M2_Tr3" = tidy_fit6),
                      .id = "Model")

did_est1
    Model               term      estimate    std.error  statistic    p.value
1  M1_Tr1           shooting -5.210841e-01 6.456718e-01 -0.8070417 0.42354260
2  M2_Tr1           shooting -7.097922e-01 6.237320e-01 -1.1379763 0.26066411
3  M2_Tr1         population  8.028568e-06 5.688141e-06  1.4114574 0.16442833
4  M2_Tr1          non_white -3.483443e+01 1.557543e+01 -2.2364985 0.02990378
5  M2_Tr1   change_unem_rate  1.592020e-01 6.584450e-02  2.4178480 0.01937671
6  M1_Tr2     fatal_shooting -6.779229e-01 6.191749e-01 -1.0948810 0.27892152
7  M2_Tr2     fatal_shooting -9.179995e-01 6.576559e-01 -1.3958659 0.16904785
8  M2_Tr2         population  8.026778e-06 5.747531e-06  1.3965611 0.16883977
9  M2_Tr2          non_white -3.484436e+01 1.561829e+01 -2.2309975 0.03029042
10 M2_Tr2   change_unem_rate  1.591126e-01 6.591001e-02  2.4140886 0.01955574
11 M1_Tr3 non_fatal_shooting -2.387205e-01 1.093832e+00 -0.2182423 0.82814672
12 M2_Tr3 non_fatal_shooting -3.269742e-01 1.145303e+00 -0.2854915 0.77647098
13 M2_Tr3         population  7.821482e-06 5.712400e-06  1.3692110 0.17717684
14 M2_Tr3          non_white -3.481384e+01 1.561438e+01 -2.2296020 0.03038921
15 M2_Tr3   change_unem_rate  1.595359e-01 6.590170e-02  2.4208169 0.01923637
        conf.low     conf.high df outcome
1  -1.818610e+00  7.764420e-01 49 turnout
2  -1.963229e+00  5.436442e-01 49 turnout
3  -3.402178e-06  1.945932e-05 49 turnout
4  -6.613442e+01 -3.534428e+00 49 turnout
5   2.688251e-02  2.915215e-01 49 turnout
6  -1.922202e+00  5.663557e-01 49 turnout
7  -2.239609e+00  4.036096e-01 49 turnout
8  -3.523318e-06  1.957687e-05 49 turnout
9  -6.623047e+01 -3.458237e+00 49 turnout
10  2.666148e-02  2.915637e-01 49 turnout
11 -2.436859e+00  1.959418e+00 49 turnout
12 -2.628546e+00  1.974598e+00 49 turnout
13 -3.658017e-06  1.930098e-05 49 turnout
14 -6.619211e+01 -3.435580e+00 49 turnout
15  2.710152e-02  2.919704e-01 49 turnout

 続いて、処置効果のみを抽出する。処置効果はterm列の値が"shooting""fatal_shooting""non_fatal_shooting"のいずれかと一致する行であるため、filter()関数を使用する。

did_est1 <- did_est1 |>
  filter(term %in% c("shooting", "fatal_shooting", "non_fatal_shooting"))

did_est1
   Model               term   estimate std.error  statistic   p.value  conf.low
1 M1_Tr1           shooting -0.5210841 0.6456718 -0.8070417 0.4235426 -1.818610
2 M2_Tr1           shooting -0.7097922 0.6237320 -1.1379763 0.2606641 -1.963229
3 M1_Tr2     fatal_shooting -0.6779229 0.6191749 -1.0948810 0.2789215 -1.922202
4 M2_Tr2     fatal_shooting -0.9179995 0.6576559 -1.3958659 0.1690478 -2.239609
5 M1_Tr3 non_fatal_shooting -0.2387205 1.0938325 -0.2182423 0.8281467 -2.436859
6 M2_Tr3 non_fatal_shooting -0.3269742 1.1453027 -0.2854915 0.7764710 -2.628546
  conf.high df outcome
1 0.7764420 49 turnout
2 0.5436442 49 turnout
3 0.5663557 49 turnout
4 0.4036096 49 turnout
5 1.9594182 49 turnout
6 1.9745977 49 turnout

 ちなみにgrepl()関数を使うと、"shooting"が含まれる行を抽出することもできる。以下のコードは上記のコードと同じ機能をする。

did_est1 <- did_est1 |>
  filter(grepl("shooting", term))

 つづいて、Model列をModelTreat列へ分割する。

did_est1 <- did_est1 |>
  separate(col  = Model,
           into = c("Model", "Treat"),
           sep  = "_")

did_est1
  Model Treat               term   estimate std.error  statistic   p.value
1    M1   Tr1           shooting -0.5210841 0.6456718 -0.8070417 0.4235426
2    M2   Tr1           shooting -0.7097922 0.6237320 -1.1379763 0.2606641
3    M1   Tr2     fatal_shooting -0.6779229 0.6191749 -1.0948810 0.2789215
4    M2   Tr2     fatal_shooting -0.9179995 0.6576559 -1.3958659 0.1690478
5    M1   Tr3 non_fatal_shooting -0.2387205 1.0938325 -0.2182423 0.8281467
6    M2   Tr3 non_fatal_shooting -0.3269742 1.1453027 -0.2854915 0.7764710
   conf.low conf.high df outcome
1 -1.818610 0.7764420 49 turnout
2 -1.963229 0.5436442 49 turnout
3 -1.922202 0.5663557 49 turnout
4 -2.239609 0.4036096 49 turnout
5 -2.436859 1.9594182 49 turnout
6 -2.628546 1.9745977 49 turnout

 可視化に入る前にModel列とTreat列の値を修正する。Model列の値が"M1"なら"County-Year FE"に、それ以外なら"County-Year FE + Covariates"とリコーディングする。戻り値が2種類だからif_else()を使う。Treat列の場合、戻り値が3つなので、recode()case_when()を使う。ここではrecode()を使ってリコーディングする。最後にModelTreatを表示順番でfactor化し(fct_inorder())、更に順番を逆転する(fct_rev())。

did_est1 <- did_est1 |>
  mutate(Model = if_else(Model == "M1",
                           "County-Year FE", 
                           "County-Year FE + Covariates"),
         Treat = recode(Treat,
                        "Tr1" = "Any Shooting (t-1)",
                        "Tr2" = "Fatal Shooting (t-1)",
                        "Tr3" = "Nonfatal Shooting (t-1)"),
         Model = fct_rev(fct_inorder(Model)),
         Treat = fct_rev(fct_inorder(Treat)))

did_est1
                        Model                   Treat               term
1              County-Year FE      Any Shooting (t-1)           shooting
2 County-Year FE + Covariates      Any Shooting (t-1)           shooting
3              County-Year FE    Fatal Shooting (t-1)     fatal_shooting
4 County-Year FE + Covariates    Fatal Shooting (t-1)     fatal_shooting
5              County-Year FE Nonfatal Shooting (t-1) non_fatal_shooting
6 County-Year FE + Covariates Nonfatal Shooting (t-1) non_fatal_shooting
    estimate std.error  statistic   p.value  conf.low conf.high df outcome
1 -0.5210841 0.6456718 -0.8070417 0.4235426 -1.818610 0.7764420 49 turnout
2 -0.7097922 0.6237320 -1.1379763 0.2606641 -1.963229 0.5436442 49 turnout
3 -0.6779229 0.6191749 -1.0948810 0.2789215 -1.922202 0.5663557 49 turnout
4 -0.9179995 0.6576559 -1.3958659 0.1690478 -2.239609 0.4036096 49 turnout
5 -0.2387205 1.0938325 -0.2182423 0.8281467 -2.436859 1.9594182 49 turnout
6 -0.3269742 1.1453027 -0.2854915 0.7764710 -2.628546 1.9745977 49 turnout

 それでは{ggplot2}を使ってpoint-rangeプロットを作成してみよう。

did_est1 |>
  ggplot() +
  # x = 0の箇所に垂直線を引く。垂直線は破線(linetype = 2)とする。
  geom_vline(xintercept = 0, linetype = 2) +
  geom_pointrange(aes(x = estimate, xmin = conf.low, xmax = conf.high,
                      y = Treat, color = Model),
                  position = position_dodge2(1/2)) +
  labs(x = "Change in Turnout (%p)", y = "", color = "") +
  # 色を指定する。
  # Modelの値が County-Year FE なら黒、
  # County-Year FE + Covariates ならグレー、
  scale_color_manual(values = c("County-Year FE" = "black", 
                                "County-Year FE + Covariates" = "gray50")) +
  # 横軸の下限と上限を-10〜10とする。
  coord_cartesian(xlim = c(-10, 10)) +
  theme_bw(base_size = 12) +
  theme(legend.position = "bottom")

 元の論文を見ると、点の上に点推定値が書かれているが、私たちもこれを真似してみよう。文字列をプロットするレイヤーはgeom_text()geom_label()annotate()があるが、ここではgeom_text()を使用する。文字列が表示される横軸上の位置(x)と縦軸上の位置(y)、そして出力する文字列(label)をマッピングする。点推定値は3桁まで出力したいので、sprintf()を使って、3桁に丸める。ただし、これだけだと点と文字が重なってしまう。vjust-0.75にすることで、出力する文字列を点の位置を上の方向へ若干ずらすことができる。

did_est1 |>
  ggplot() +
  geom_vline(xintercept = 0, linetype = 2) +
  geom_pointrange(aes(x = estimate, xmin = conf.low, xmax = conf.high,
                      y = Treat, color = Model),
                  position = position_dodge2(1/2)) +
  geom_text(aes(x = estimate, y = Treat, color = Model, 
                label = sprintf("%.3f", estimate)),
            position = position_dodge2(1/2),
            vjust = -0.75) +
  labs(x = "Change in Turnout (%p)", y = "", color = "") +
  scale_color_manual(values = c("County-Year FE" = "black", 
                                "County-Year FE + Covariates" = "gray50")) +
  coord_cartesian(xlim = c(-10, 10)) +
  theme_bw(base_size = 12) +
  theme(legend.position = "bottom")

 ちなみにこのコードを見ると、geom_pointrange()geom_text()xycolorを共有しているので、ggplot()内でマッピングすることもできる。

did_est1 |>
  ggplot(aes(x = estimate, y = Treat, color = Model)) +
  geom_vline(xintercept = 0, linetype = 2) +
  geom_pointrange(aes(xmin = conf.low, xmax = conf.high),
                  position = position_dodge2(1/2)) +
  geom_text(aes(label = sprintf("%.3f", estimate)),
            position = position_dodge2(1/2),
            vjust = -0.75) +
  labs(x = "Change in Turnout (%p)", y = "", color = "") +
  scale_color_manual(values = c("County-Year FE" = "black", 
                                "County-Year FE + Covariates" = "gray50")) +
  coord_cartesian(xlim = c(-10, 10)) +
  theme_bw(base_size = 12) +
  theme(legend.position = "bottom")

 続いて、民主党候補者の得票率(demvote)を応答変数として6つのモデルを推定し、同じ作業を繰り返す。

did_fit7 <- lm_robust(demvote ~ shooting, 
                      data          = did_df, 
                      fixed_effects = ~year_f + county_f,
                      clusters      = state_f,
                      se_type       = "stata")

did_fit8 <- lm_robust(demvote ~ shooting + 
                        population + non_white + change_unem_rate, 
                      data          = did_df, 
                      fixed_effects = ~year_f + county_f,
                      clusters      = state_f,
                      se_type       = "stata")

did_fit9 <- lm_robust(demvote ~ fatal_shooting, 
                      data          = did_df, 
                      fixed_effects = ~year_f + county_f,
                      clusters      = state_f,
                      se_type       = "stata")

did_fit10 <- lm_robust(demvote ~ fatal_shooting + 
                         population + non_white + change_unem_rate, 
                       data          = did_df, 
                       fixed_effects = ~year_f + county_f,
                       clusters      = state_f,
                       se_type       = "stata")

did_fit11 <- lm_robust(demvote ~ non_fatal_shooting, 
                       data          = did_df, 
                       fixed_effects = ~year_f + county_f,
                       clusters      = state_f,
                       se_type       = "stata")

did_fit12 <- lm_robust(demvote ~ non_fatal_shooting + 
                         population + non_white + change_unem_rate, 
                       data          = did_df, 
                       fixed_effects = ~year_f + county_f,
                       clusters      = state_f,
                       se_type       = "stata")

modelsummary(list("Model 7"  = did_fit7,  "Model 8"  = did_fit8, 
                  "Model 9"  = did_fit9,  "Model 10" = did_fit10, 
                  "Model 11" = did_fit11, "Model 12" = did_fit12))
Model 7 Model 8 Model 9 Model 10 Model 11 Model 12
shooting 4.513 2.364
(0.875) (0.737)
population 0.000 0.000 0.000
(0.000) (0.000) (0.000)
non_white 86.873 86.866 86.796
(17.863) (17.855) (17.855)
change_unem_rate −0.139 −0.139 −0.140
(0.127) (0.128) (0.127)
fatal_shooting 4.404 1.782
(1.092) (0.836)
non_fatal_shooting 4.217 2.930
(1.298) (1.037)
Num.Obs. 18627 18627 18627 18627 18627 18627
R2 0.882 0.894 0.882 0.894 0.882 0.894
R2 Adj. 0.859 0.873 0.859 0.873 0.859 0.873
AIC 110745.2 108755.4 110772.0 108768.2 110786.5 108762.1
BIC 110760.8 108794.6 110787.7 108807.3 110802.2 108801.3
RMSE 4.73 4.48 4.73 4.48 4.73 4.48
Std.Errors by: state_f by: state_f by: state_f by: state_f by: state_f by: state_f

 今回はいずれも統計的に有意な結果が得られている。例えば、モデル7(did_fit7)の場合、処置効果の推定値は4.513である。これは学校内銃撃事件が発生したカウンティーの場合、大統領選挙において民主党候補者の得票率が約4.513%p増加することを意味する。

 以上の結果を図としてまとめてみよう。

tidy_fit7  <- tidy(did_fit7)
tidy_fit8  <- tidy(did_fit8)
tidy_fit9  <- tidy(did_fit9)
tidy_fit10 <- tidy(did_fit10)
tidy_fit11 <- tidy(did_fit11)
tidy_fit12 <- tidy(did_fit12)

did_est2 <- bind_rows(list("M1_Tr1" = tidy_fit7,
                           "M2_Tr1" = tidy_fit8,
                           "M1_Tr2" = tidy_fit9,
                           "M2_Tr2" = tidy_fit10,
                           "M1_Tr3" = tidy_fit11,
                           "M2_Tr3" = tidy_fit12),
                      .id = "Model")

did_est2
    Model               term      estimate    std.error statistic      p.value
1  M1_Tr1           shooting  4.513237e+00 8.752702e-01  5.156393 4.514820e-06
2  M2_Tr1           shooting  2.363839e+00 7.367346e-01  3.208536 2.353748e-03
3  M2_Tr1         population  3.174252e-05 6.778757e-06  4.682646 2.275534e-05
4  M2_Tr1          non_white  8.687344e+01 1.786287e+01  4.863352 1.234074e-05
5  M2_Tr1   change_unem_rate -1.389392e-01 1.274198e-01 -1.090405 2.808682e-01
6  M1_Tr2     fatal_shooting  4.404109e+00 1.091944e+00  4.033274 1.920110e-04
7  M2_Tr2     fatal_shooting  1.782393e+00 8.363792e-01  2.131083 3.812471e-02
8  M2_Tr2         population  3.206812e-05 6.801197e-06  4.715070 2.039958e-05
9  M2_Tr2          non_white  8.686629e+01 1.785535e+01  4.865000 1.227169e-05
10 M2_Tr2   change_unem_rate -1.392341e-01 1.275138e-01 -1.091914 2.802109e-01
11 M1_Tr3 non_fatal_shooting  4.216683e+00 1.298056e+00  3.248460 2.098382e-03
12 M2_Tr3 non_fatal_shooting  2.929866e+00 1.037217e+00  2.824737 6.825655e-03
13 M2_Tr3         population  3.229215e-05 6.717999e-06  4.806810 1.495567e-05
14 M2_Tr3          non_white  8.679618e+01 1.785514e+01  4.861131 1.243440e-05
15 M2_Tr3   change_unem_rate -1.400321e-01 1.274994e-01 -1.098296 2.774427e-01
        conf.low    conf.high df outcome
1   2.754316e+00 6.272159e+00 49 demvote
2   8.833157e-01 3.844363e+00 49 demvote
3   1.812010e-05 4.536494e-05 49 demvote
4   5.097665e+01 1.227702e+02 49 demvote
5  -3.949989e-01 1.171206e-01 49 demvote
6   2.209766e+00 6.598453e+00 49 demvote
7   1.016261e-01 3.463160e+00 49 demvote
8   1.840061e-05 4.573564e-05 49 demvote
9   5.098461e+01 1.227480e+02 49 demvote
10 -3.954828e-01 1.170146e-01 49 demvote
11  1.608142e+00 6.825225e+00 49 demvote
12  8.454996e-01 5.014232e+00 49 demvote
13  1.879182e-05 4.579247e-05 49 demvote
14  5.091493e+01 1.226774e+02 49 demvote
15 -3.962518e-01 1.161876e-01 49 demvote
did_est2 <- did_est2 |>
  filter(grepl("shooting", term))

did_est2
   Model               term estimate std.error statistic      p.value  conf.low
1 M1_Tr1           shooting 4.513237 0.8752702  5.156393 4.514820e-06 2.7543158
2 M2_Tr1           shooting 2.363839 0.7367346  3.208536 2.353748e-03 0.8833157
3 M1_Tr2     fatal_shooting 4.404109 1.0919440  4.033274 1.920110e-04 2.2097656
4 M2_Tr2     fatal_shooting 1.782393 0.8363792  2.131083 3.812471e-02 0.1016261
5 M1_Tr3 non_fatal_shooting 4.216683 1.2980560  3.248460 2.098382e-03 1.6081422
6 M2_Tr3 non_fatal_shooting 2.929866 1.0372173  2.824737 6.825655e-03 0.8454996
  conf.high df outcome
1  6.272159 49 demvote
2  3.844363 49 demvote
3  6.598453 49 demvote
4  3.463160 49 demvote
5  6.825225 49 demvote
6  5.014232 49 demvote
did_est2 <- did_est2 |>
  separate(col  = Model,
           into = c("Model", "Treat"),
           sep  = "_")

did_est2
  Model Treat               term estimate std.error statistic      p.value
1    M1   Tr1           shooting 4.513237 0.8752702  5.156393 4.514820e-06
2    M2   Tr1           shooting 2.363839 0.7367346  3.208536 2.353748e-03
3    M1   Tr2     fatal_shooting 4.404109 1.0919440  4.033274 1.920110e-04
4    M2   Tr2     fatal_shooting 1.782393 0.8363792  2.131083 3.812471e-02
5    M1   Tr3 non_fatal_shooting 4.216683 1.2980560  3.248460 2.098382e-03
6    M2   Tr3 non_fatal_shooting 2.929866 1.0372173  2.824737 6.825655e-03
   conf.low conf.high df outcome
1 2.7543158  6.272159 49 demvote
2 0.8833157  3.844363 49 demvote
3 2.2097656  6.598453 49 demvote
4 0.1016261  3.463160 49 demvote
5 1.6081422  6.825225 49 demvote
6 0.8454996  5.014232 49 demvote
did_est2 <- did_est2 |>
  mutate(Model = if_else(Model == "M1",
                           "County-Year FE", 
                           "County-Year FE + Covariates"),
         Treat = recode(Treat,
                        "Tr1" = "Any Shooting (t-1)",
                        "Tr2" = "Fatal Shooting (t-1)",
                        "Tr3" = "Nonfatal Shooting (t-1)"),
         Model = fct_rev(fct_inorder(Model)),
         Treat = fct_rev(fct_inorder(Treat)))

did_est2
                        Model                   Treat               term
1              County-Year FE      Any Shooting (t-1)           shooting
2 County-Year FE + Covariates      Any Shooting (t-1)           shooting
3              County-Year FE    Fatal Shooting (t-1)     fatal_shooting
4 County-Year FE + Covariates    Fatal Shooting (t-1)     fatal_shooting
5              County-Year FE Nonfatal Shooting (t-1) non_fatal_shooting
6 County-Year FE + Covariates Nonfatal Shooting (t-1) non_fatal_shooting
  estimate std.error statistic      p.value  conf.low conf.high df outcome
1 4.513237 0.8752702  5.156393 4.514820e-06 2.7543158  6.272159 49 demvote
2 2.363839 0.7367346  3.208536 2.353748e-03 0.8833157  3.844363 49 demvote
3 4.404109 1.0919440  4.033274 1.920110e-04 2.2097656  6.598453 49 demvote
4 1.782393 0.8363792  2.131083 3.812471e-02 0.1016261  3.463160 49 demvote
5 4.216683 1.2980560  3.248460 2.098382e-03 1.6081422  6.825225 49 demvote
6 2.929866 1.0372173  2.824737 6.825655e-03 0.8454996  5.014232 49 demvote
did_est2 |>
  ggplot() +
  geom_vline(xintercept = 0, linetype = 2) +
  geom_pointrange(aes(x = estimate, xmin = conf.low, xmax = conf.high,
                      y = Treat, color = Model),
                  position = position_dodge2(1/2)) +
  geom_text(aes(x = estimate, y = Treat, color = Model, 
                label = sprintf("%.3f", estimate)),
            position = position_dodge2(1/2),
            vjust = -0.75) +
  labs(x = "Change in Democratic Vote Share (%p)", y = "", color = "") +
  scale_color_manual(values = c("County-Year FE" = "black", 
                                "County-Year FE + Covariates" = "gray50")) +
  coord_cartesian(xlim = c(-10, 10)) +
  theme_bw(base_size = 12) +
  theme(legend.position = "bottom")

 最後に、これまで作成した2つの図を一つにまとめてみよう。bind_rows()関数を使い、それぞれの表に識別子(Outcome)を与える。

did_est <- bind_rows(list("Out1" = did_est1,
                          "Out2" = did_est2),
                     .id = "Outcome")

did_est
   Outcome                       Model                   Treat
1     Out1              County-Year FE      Any Shooting (t-1)
2     Out1 County-Year FE + Covariates      Any Shooting (t-1)
3     Out1              County-Year FE    Fatal Shooting (t-1)
4     Out1 County-Year FE + Covariates    Fatal Shooting (t-1)
5     Out1              County-Year FE Nonfatal Shooting (t-1)
6     Out1 County-Year FE + Covariates Nonfatal Shooting (t-1)
7     Out2              County-Year FE      Any Shooting (t-1)
8     Out2 County-Year FE + Covariates      Any Shooting (t-1)
9     Out2              County-Year FE    Fatal Shooting (t-1)
10    Out2 County-Year FE + Covariates    Fatal Shooting (t-1)
11    Out2              County-Year FE Nonfatal Shooting (t-1)
12    Out2 County-Year FE + Covariates Nonfatal Shooting (t-1)
                 term   estimate std.error  statistic      p.value   conf.low
1            shooting -0.5210841 0.6456718 -0.8070417 4.235426e-01 -1.8186103
2            shooting -0.7097922 0.6237320 -1.1379763 2.606641e-01 -1.9632286
3      fatal_shooting -0.6779229 0.6191749 -1.0948810 2.789215e-01 -1.9222015
4      fatal_shooting -0.9179995 0.6576559 -1.3958659 1.690478e-01 -2.2396086
5  non_fatal_shooting -0.2387205 1.0938325 -0.2182423 8.281467e-01 -2.4368592
6  non_fatal_shooting -0.3269742 1.1453027 -0.2854915 7.764710e-01 -2.6285461
7            shooting  4.5132372 0.8752702  5.1563930 4.514820e-06  2.7543158
8            shooting  2.3638394 0.7367346  3.2085357 2.353748e-03  0.8833157
9      fatal_shooting  4.4041092 1.0919440  4.0332738 1.920110e-04  2.2097656
10     fatal_shooting  1.7823931 0.8363792  2.1310825 3.812471e-02  0.1016261
11 non_fatal_shooting  4.2166835 1.2980560  3.2484602 2.098382e-03  1.6081422
12 non_fatal_shooting  2.9298657 1.0372173  2.8247368 6.825655e-03  0.8454996
   conf.high df outcome
1  0.7764420 49 turnout
2  0.5436442 49 turnout
3  0.5663557 49 turnout
4  0.4036096 49 turnout
5  1.9594182 49 turnout
6  1.9745977 49 turnout
7  6.2721585 49 demvote
8  3.8443631 49 demvote
9  6.5984529 49 demvote
10 3.4631600 49 demvote
11 6.8252248 49 demvote
12 5.0142319 49 demvote

 Outcome列のリコーディングし、factor化する。

did_est <- did_est |>
  mutate(Outcome = if_else(Outcome == "Out1",
                           "Change in Turnout (%p)",
                           "Change in Democratic Vote Share (%p)"),
         Outcome = fct_inorder(Outcome))

did_est
                                Outcome                       Model
1                Change in Turnout (%p)              County-Year FE
2                Change in Turnout (%p) County-Year FE + Covariates
3                Change in Turnout (%p)              County-Year FE
4                Change in Turnout (%p) County-Year FE + Covariates
5                Change in Turnout (%p)              County-Year FE
6                Change in Turnout (%p) County-Year FE + Covariates
7  Change in Democratic Vote Share (%p)              County-Year FE
8  Change in Democratic Vote Share (%p) County-Year FE + Covariates
9  Change in Democratic Vote Share (%p)              County-Year FE
10 Change in Democratic Vote Share (%p) County-Year FE + Covariates
11 Change in Democratic Vote Share (%p)              County-Year FE
12 Change in Democratic Vote Share (%p) County-Year FE + Covariates
                     Treat               term   estimate std.error  statistic
1       Any Shooting (t-1)           shooting -0.5210841 0.6456718 -0.8070417
2       Any Shooting (t-1)           shooting -0.7097922 0.6237320 -1.1379763
3     Fatal Shooting (t-1)     fatal_shooting -0.6779229 0.6191749 -1.0948810
4     Fatal Shooting (t-1)     fatal_shooting -0.9179995 0.6576559 -1.3958659
5  Nonfatal Shooting (t-1) non_fatal_shooting -0.2387205 1.0938325 -0.2182423
6  Nonfatal Shooting (t-1) non_fatal_shooting -0.3269742 1.1453027 -0.2854915
7       Any Shooting (t-1)           shooting  4.5132372 0.8752702  5.1563930
8       Any Shooting (t-1)           shooting  2.3638394 0.7367346  3.2085357
9     Fatal Shooting (t-1)     fatal_shooting  4.4041092 1.0919440  4.0332738
10    Fatal Shooting (t-1)     fatal_shooting  1.7823931 0.8363792  2.1310825
11 Nonfatal Shooting (t-1) non_fatal_shooting  4.2166835 1.2980560  3.2484602
12 Nonfatal Shooting (t-1) non_fatal_shooting  2.9298657 1.0372173  2.8247368
        p.value   conf.low conf.high df outcome
1  4.235426e-01 -1.8186103 0.7764420 49 turnout
2  2.606641e-01 -1.9632286 0.5436442 49 turnout
3  2.789215e-01 -1.9222015 0.5663557 49 turnout
4  1.690478e-01 -2.2396086 0.4036096 49 turnout
5  8.281467e-01 -2.4368592 1.9594182 49 turnout
6  7.764710e-01 -2.6285461 1.9745977 49 turnout
7  4.514820e-06  2.7543158 6.2721585 49 demvote
8  2.353748e-03  0.8833157 3.8443631 49 demvote
9  1.920110e-04  2.2097656 6.5984529 49 demvote
10 3.812471e-02  0.1016261 3.4631600 49 demvote
11 2.098382e-03  1.6081422 6.8252248 49 demvote
12 6.825655e-03  0.8454996 5.0142319 49 demvote

 図の作り方はこれまでと変わらないが、ファセット分割を行うため、facet_wrap()レイヤーを追加する。

did_est |>
  ggplot() +
  geom_vline(xintercept = 0, linetype = 2) +
  geom_pointrange(aes(x = estimate, xmin = conf.low, xmax = conf.high,
                      y = Treat, color = Model),
                  position = position_dodge2(1/2)) +
  geom_text(aes(x = estimate, y = Treat, color = Model, 
                label = sprintf("%.3f", estimate)),
            position = position_dodge2(1/2),
            vjust = -0.75) +
  labs(x = "Treatment Effects", y = "", color = "") +
  scale_color_manual(values = c("County-Year FE" = "black", 
                                "County-Year FE + Covariates" = "gray50")) +
  coord_cartesian(xlim = c(-10, 10)) +
  facet_wrap(~Outcome, ncol = 2) +
  theme_bw(base_size = 12) +
  theme(legend.position = "bottom")

 以上の結果から「学校内銃撃事件の発生は投票参加を促すとは言えないものの、民主党候補者の得票率を上げる」ということが言えよう。

SCM

 ここでは、差分の差分法の応用としてSynthetic Control Method(SCM)を{gsynth}パッケージの使い方に重点を起きながら説明する。Synthetic Control MethodをRに実装したパッケージは{Synth}、{gsynth}、{bpCausal}などがある。SCMの代表的な論文の一つであるAbadie et al.(2015)は{Synth}パッケージを使っているが、使い方はかなり複雑である。したがって、ここでは使い方が最も簡単な{gsynth}パッケージについて説明する2

 ここで使用するのは歴代参院選(第7回以降)における都道府県の投票率である3。処置変数は第24回参院選から導入された「合区」だ。具体的には鳥取選挙区と島根選挙区が「鳥取・島根選挙区」に、徳島選挙区と高知選挙区が「徳島・高知選挙区」として合区されたか否かである。合区によって自分の1票の価値がほぼ半分に下がり、投票参加を妨げたのではないかという議論もあるが、本当だろうか。

scm_df <- read_csv("data/did_data6.csv")

scm_df
# A tibble: 938 × 9
   Election PrefCode PrefName_E PrefNam…¹ Eligi…² Voters EffVote Magni…³ Candi…⁴
      <dbl>    <dbl> <chr>      <chr>       <dbl>  <dbl>   <dbl>   <dbl>   <dbl>
 1        7        1 Hokkaido   北海道    2903802 1.83e6 1683764       4       7
 2        7        2 Aomori     青森県     823336 5.09e5  478380       1       5
 3        7        3 Iwate      岩手県     849109 6.20e5  591977       1       3
 4        7        4 Miyagi     宮城県    1049109 6.95e5  666205       1       4
 5        7        5 Akita      秋田県     771983 5.91e5  558603       1       3
 6        7        6 Yamagata   山形県     784438 6.37e5  615414       1       3
 7        7        7 Fukushima  福島県    1123975 8.93e5  855609       2       5
 8        7        8 Ibaraki    茨城県    1246254 7.85e5  750181       2       4
 9        7        9 Tochigi    栃木県     909245 6.38e5  614249       2       6
10        7       10 Gunma      群馬県     990083 7.20e5  690448       2       6
# … with 928 more rows, and abbreviated variable names ¹​PrefName_J, ²​Eligible,
#   ³​Magnitude, ⁴​Candidates
# ℹ Use `print(n = ...)` to see more rows
scm_df |>
  dplyr::select(Eligible:Candidates) |>
  descr(stats = c("mean", "sd", "min", "max", "n.valid"),
        transpose = TRUE, order = "p")
Error : Can't find dplyr
Descriptive Statistics  

                         Mean      Std.Dev         Min           Max   N.Valid
---------------- ------------ ------------ ----------- ------------- ---------
        Eligible   1933435.87   1867850.63   364184.00   11454822.00    938.00
          Voters   1139531.46   1029039.25   226580.00    6477709.00    938.00
         EffVote   1096934.10    997876.11   217637.00    6298466.00    938.00
       Magnitude         1.60         0.87        1.00          6.00    938.00
      Candidates         5.49         5.42        2.00         72.00    938.00
変数名 説明 備考
Election 選挙 XX回参議院議員通常選挙
PrefCode 都道府県 JIS規格コード
PrefName_E 都道府県 英語
PrefName_J 都道府県 日本語
Eligible 有権者数
Voters 投票者数
EffVote 有効投票数
Magnitude 定数
Candidates 候補者数

 合区された選挙区における定数および候補者数の扱いについてだが、たとえば島根選挙区の場合、合区前の定数は1である。合区後の鳥取・島根選挙区も定数1だ。この場合、合区後の島根「県」の定数を1にするか0.5にするか、そして候補者数も合区における候補者数にすべきか、2に割るかが問題となるが、ここでは定数を1とし、候補者数も合区における候補者数とした。このようなコーディングが実証上、正しいかどうかは分からないが、実習の段階では問題ないだろう。

 それでは、データを以下のように加工する。

  1. Treatという変数を作成する。PrefNameが「鳥取県」、「島根県」、「徳島県」、「高知県」であり、Electionが24以上なら1、それ以外のケースは0とする。
  2. Turnoutという変数を作成する。投票者数(Voters)を有権者数(Eligible)で割り、100を掛ける。
  3. Spoiltという変数を作成する。有効投票数(EffVote)を有権者数(Eligible)で割った値を1から引き、100を掛ける。
  4. 都道府県名(PrefName_E)をfactor変数に変換し、PrefName_Fという列として追加します。水準(levels)の順番はデータの出現順とする。
  5. 都道府県名(PrefName_J)をfactor化し、要素の順番はscm_df内における表示順番とする。
scm_df <- scm_df |>
  mutate(Treat      = if_else(PrefName_E %in% c("Tottori", "Shimane",
                                                "Tokushima","Kochi") &
                                Election >= 24, 1, 0),
         Turnout    = Voters / Eligible * 100,
         Spoilt     = (1 - (EffVote / Voters)) * 100,
         PrefName_J = fct_inorder(PrefName_J))

 本格的な分析に入る前に、各ケースの処置有無を可視化してみよう。パネルデータの確認には{panelView}パッケージのpanelview()関数を使う。第一引数は応答変数 ~ 処置変数であり、data引数にデータフレームのオブジェクト名を指定する。また、indexにはユニットと時間を表す変数名を長さ2のcharacter型ベクトルとして指定します。ここでは"PrefName_J""Election"である。最後に、pre.post = TRUEを指定すると、処置前後を色分けしてくれるので、見やすくなる。

panelview(Turnout ~ Treat, data = scm_df, 
          index = c("PrefName_J", "Election"), pre.post = TRUE,
          xlab = "参院選", ylab = "都道府県")

 これを見ると処置を受けるケースが4県であり、どれも第24回参院選から処置を受けることが分かる。また、沖縄県の場合、第7・8回参院選のデータが欠損していることが分かる。

 次は投票率の変化を時系列的に示してみよう。使い方は先ほどとほぼ同じだが、panelView()内にtype = "outcome"引数を追加する。

panelview(Turnout ~ Treat, data = scm_df, 
          index = c("PrefName_J", "Election"), 
          type  = "outcome",
          main  = "投票率の推移",
          xlab  = "参院選",
          ylab  = "投票率 (%)")

 これを見ると、都道府県内における投票率の変化にはバラツキがあるが、傾向としてはかなり似通っていることが分かる。また、処置後の変化(濃い青の線)を見ると、他の都道府県よりかなり落ちているようにも見える。これは合区によって何かの変化が生じた可能性があることを示唆している。

 それではSCMをやってみよう。基本的な使い方はlm()関数に近いが、それでも必要な引数がそれなりにある。ここではその一部を紹介する。

  • 第1引数は応答変数 ~ 処置変数 + 統制変数1 + 統制変数2 + ...であり、統制変数は必須ではない。ここでは統制変数として定数(Magnitude)と候補者数(Candidates)を使用する。
  • dataにはデータフレームのオブジェクト名を指定する。
  • indexpanelview()と同じだが、一つ注意が必要です。それは、ユニットを表す変数をこれまではfactor化された"PrefName_J"を使ってきたが、factor化されていない"PrefName_E"にすることだ。まだ開発途上のパッケージということもあり、ユニット変数がfactor型だと、正しく推定できない。PrefName_Eはfactor型でなくcharacter型ですので、問題ない(あえてPrefName_Eをfactor化しなかったのはこのためである)。また、numeric型であるPrefCodeを指定すれば良い。。
  • forceは固定効果をユニットレベルにするか、時間レベルにするか、両方にするか、あるいはなしにするかを意味する。既定値はユニットレベル("unit")だが、ここでは両方("two-way")にする。
  • seは標準誤差を計算するか否かであり、既定値はFALSEである。ここでは標準誤差も計算するためにTRUEに指定する。標準誤差を計算しない場合、計算が早く終わる。
  • inferenceは推定方法を意味し、既定値はノンパラメトリック推定("nonparametric")です。ただし、処置ケースが少ない場合、パラメトリック推定が推奨されているため、ここでは"parametric"を指定する。
  • nbootsは標準誤差を計算する際に使用されるブートストラッピングの回数である。既定値は200ですが、ここでは500回にする。
  • parallelは並列計算の有無を意味する。既定値はTRUEで、このままでも通常は問題ない。ただし、R Markdown上で{gsynth}を使用する場合は並列計算ができないため、FALSEにしておく(処理時間が倍以上になる)。
scm_fit <- gsynth(Turnout ~ Treat + Magnitude + Candidates, 
                  data = scm_df, index = c("PrefName_E", "Election"), 
                  force = "two-way", se = TRUE, inference = "parametric",
                  nboots = 500, parallel = TRUE)
print(scm_fit)
Call:
gsynth.formula(formula = Turnout ~ Treat + Magnitude + Candidates, 
    data = scm_df, index = c("PrefName_E", "Election"), force = "two-way", 
    se = TRUE, nboots = 500, inference = "parametric", parallel = FALSE)

Average Treatment Effect on the Treated:
        Estimate  S.E. CI.lower CI.upper   p.value
ATT.avg   -6.563 1.463    -9.43   -3.696 7.242e-06

   ~ by Period (including Pre-treatment Periods):
         ATT   S.E. CI.lower CI.upper   p.value n.Treated
-16  0.09054 1.2301  -2.3205  2.50155 9.413e-01         0
-15 -0.56609 1.1603  -2.8403  1.70812 6.256e-01         0
-14  3.32790 1.0443   1.2811  5.37466 1.439e-03         0
-13  1.44804 0.9972  -0.5064  3.40244 1.465e-01         0
-12 -1.53058 1.0356  -3.5602  0.49908 1.394e-01         0
-11 -1.21442 0.8801  -2.9393  0.51048 1.676e-01         0
-10 -2.51664 1.2571  -4.9806 -0.05272 4.530e-02         0
-9  -1.69683 1.0860  -3.8253  0.43163 1.182e-01         0
-8   1.12683 0.9586  -0.7520  3.00562 2.398e-01         0
-7  -1.03078 1.0960  -3.1789  1.11730 3.470e-01         0
-6   3.60723 1.5520   0.5653  6.64915 2.011e-02         0
-5  -0.85728 1.2777  -3.3615  1.64695 5.022e-01         0
-4   0.34762 0.9497  -1.5137  2.20893 7.143e-01         0
-3  -0.38462 0.9695  -2.2848  1.51553 6.916e-01         0
-2   0.26204 0.9226  -1.5463  2.07033 7.764e-01         0
-1   1.56199 0.8999  -0.2017  3.32571 8.260e-02         0
0   -1.97495 0.9112  -3.7608 -0.18913 3.019e-02         0
1   -6.57256 1.4040  -9.3243 -3.82083 2.849e-06         4
2   -6.87683 2.0956 -10.9841 -2.76959 1.032e-03         4
3   -6.24002 1.7173  -9.6059 -2.87413 2.795e-04         4

Coefficients for the Covariates:
               beta    S.E. CI.lower CI.upper   p.value
Magnitude   2.64853 0.79438   1.0916   4.2055 0.0008558
Candidates -0.05892 0.03945  -0.1362   0.0184 0.1353069

 平均値な処置効果(ATT)は-6.563であり、統計的にも有意な結果が得られた。これは合区が行われた選挙区は、もし合区しなかった場合の投票率に比べ、約-6.563%p低いことを意味します。また、~ by Period以下の欄では各選挙ごとのATTが表示されます。1以降が合区以降の時期であり、それぞれ約-6.573%p、-6.877%p、-6.240%pである。上記の-6.563はこの3つの数値の平均値である。ちなみに処置群の県が4つにも関わらず、ATTが一つだけ表示されるのは、4つの平均値を出しているからだ。それぞれの県ごとに効果を確認する方法は後ほど解説する。とりあえず、この結果から合区は投票率を下げたという解釈ができよう。

 以上の結果を可視化するためにはplot()関数を使う。タイトル、横軸、縦軸のラベルはmainxlabylabで指定できる。

plot(scm_fit, 
     main = "Estimated ATT", 
     xlab = "Election (0 = 2013 Election)", 
     ylab = "ATT (%p)")

 処置を受けるは架空の4県とその実際の4県の間に投票率の差はあまりなかったが、処置を受けてから差が広がることが分かる。差分でなく、合成された架空の4県と実際の4県のトレンドを見るためにはtype = "counterfactual"、もしくはtype = "ct"を指定します。

plot(scm_fit, 
     type = "counterfactual", 
     main = "",
     xlab = "Election", 
     ylab = "Turnout (%)")

 処置を受けたのは4県であるが、線が1本になっている。ここからも処置前は架空の4県と実際の4県はほぼ同じトレンドを示していますが、処置後は傾向が変わったことが確認できる。

 もし、4つの都道府県を個別に示したい場合はどうすれば良いだろうか。この場合、plot()内にid引数を使えば良い。たとえば、鳥取県の処置効果(ATT)を確認するためにはid = "Tottori"と指定する。

plot(scm_fit, id = "Tottori")

plot(scm_fit, id = "Tottori", type = "ct")

 4つの都道府県のデータを個別のファセットで全て出すためにはゼロベースで作図した方が効率的かも知れない。また、この場合は図のカスタマイズの幅も広がる。{gsynth}から得られたオブジェクトには架空の鳥取、島根、徳島、高知のデータが個別に格納されているため、それを抽出すれば個別の図を作成することもできる。たとえば、処置群(4県)の観察済みのデータはオブジェクト名$Y.trで抽出できる。

scm_fit$Y.tr
      Kochi  Shimane Tokushima  Tottori
7  74.63966 83.49050  69.39499 81.24437
8  75.69312 84.78405  69.88441 83.69962
9  71.37616 78.14244  63.66989 77.03502
10 75.69410 86.92480  78.45754 84.82852
11 70.90031 85.15588  63.40844 83.37384
12 72.83410 87.58638  77.05070 84.86225
13 61.27247 75.36695  49.81165 74.79201
14 72.01185 86.89197  70.26960 85.81677
15 70.52034 82.31831  65.59127 78.71451
16 54.89730 73.79272  48.42299 67.28538
17 50.63857 67.09330  47.14380 67.57092
18 56.21296 73.26601  56.91326 70.03591
19 58.38653 68.61257  57.24322 66.68463
20 57.29701 68.87462  54.59808 64.16773
21 58.39998 71.80893  58.46591 67.67363
22 58.49266 71.70006  58.23675 65.76520
23 49.88987 60.89257  49.29437 58.87696
24 45.51522 62.19714  46.97798 56.28266
25 46.34053 54.04000  38.59305 49.97934
26 47.36311 56.37101  45.72428 48.92585

 このデータは行列構造となっているため、as_tibble()で表形式へ変更し、treat_dfという名で格納する。

treat_df <- as_tibble(scm_fit$Y.tr)

treat_df
# A tibble: 20 × 4
   Kochi Shimane Tokushima Tottori
   <dbl>   <dbl>     <dbl>   <dbl>
 1  74.6    83.5      69.4    81.2
 2  75.7    84.8      69.9    83.7
 3  71.4    78.1      63.7    77.0
 4  75.7    86.9      78.5    84.8
 5  70.9    85.2      63.4    83.4
 6  72.8    87.6      77.1    84.9
 7  61.3    75.4      49.8    74.8
 8  72.0    86.9      70.3    85.8
 9  70.5    82.3      65.6    78.7
10  54.9    73.8      48.4    67.3
11  50.6    67.1      47.1    67.6
12  56.2    73.3      56.9    70.0
13  58.4    68.6      57.2    66.7
14  57.3    68.9      54.6    64.2
15  58.4    71.8      58.5    67.7
16  58.5    71.7      58.2    65.8
17  49.9    60.9      49.3    58.9
18  45.5    62.2      47.0    56.3
19  46.3    54.0      38.6    50.0
20  47.4    56.4      45.7    48.9

 このtreat_dfだけは各投票率がいつの投票率かが分からない。Election列をKochi列前に追加し(mutate()内の最後に.before = Kochiを指定すると、mutate()で生成・修正された列がKochi列の前へ移動する)、7から26までを入れる。

treat_df <- treat_df |>
  mutate(Election = 7:26, .before = Kochi)

treat_df
# A tibble: 20 × 5
   Election Kochi Shimane Tokushima Tottori
      <int> <dbl>   <dbl>     <dbl>   <dbl>
 1        7  74.6    83.5      69.4    81.2
 2        8  75.7    84.8      69.9    83.7
 3        9  71.4    78.1      63.7    77.0
 4       10  75.7    86.9      78.5    84.8
 5       11  70.9    85.2      63.4    83.4
 6       12  72.8    87.6      77.1    84.9
 7       13  61.3    75.4      49.8    74.8
 8       14  72.0    86.9      70.3    85.8
 9       15  70.5    82.3      65.6    78.7
10       16  54.9    73.8      48.4    67.3
11       17  50.6    67.1      47.1    67.6
12       18  56.2    73.3      56.9    70.0
13       19  58.4    68.6      57.2    66.7
14       20  57.3    68.9      54.6    64.2
15       21  58.4    71.8      58.5    67.7
16       22  58.5    71.7      58.2    65.8
17       23  49.9    60.9      49.3    58.9
18       24  45.5    62.2      47.0    56.3
19       25  46.3    54.0      38.6    50.0
20       26  47.4    56.4      45.7    48.9

 続いて、{tidyr}のpivot_longer()関数を使ってtreat_dfをlong型データへ変換する。pivot_*()関数の使い方については『私たちのR』第15章「整然データ構造」を参照されたい。

treat_df <- treat_df |>
  pivot_longer(cols      = Kochi:Tottori,
               names_to  = "Pref",
               values_to = "Turnout")

treat_df
# A tibble: 80 × 3
   Election Pref      Turnout
      <int> <chr>       <dbl>
 1        7 Kochi        74.6
 2        7 Shimane      83.5
 3        7 Tokushima    69.4
 4        7 Tottori      81.2
 5        8 Kochi        75.7
 6        8 Shimane      84.8
 7        8 Tokushima    69.9
 8        8 Tottori      83.7
 9        9 Kochi        71.4
10        9 Shimane      78.1
# … with 70 more rows
# ℹ Use `print(n = ...)` to see more rows

 同じ作業を架空の処置群についても行う。架空の処置群はオブジェクト名$Y.ctで抽出できる。今回は全ての作業をパイプ演算子で繋ぎ、コードを効率化する。

counter_df <- scm_fit$Y.ct |>
  as_tibble() |>
  mutate(Election = 7:26, .before = Kochi) |>
  pivot_longer(cols      = Kochi:Tottori,
               names_to  = "Pref",
               values_to = "Turnout")

counter_df
# A tibble: 80 × 3
   Election Pref      Turnout
      <int> <chr>       <dbl>
 1        7 Kochi        73.5
 2        7 Shimane      84.0
 3        7 Tokushima    67.6
 4        7 Tottori      83.3
 5        8 Kochi        75.3
 6        8 Shimane      85.7
 7        8 Tokushima    70.7
 8        8 Tottori      84.7
 9        9 Kochi        65.6
10        9 Shimane      76.1
# … with 70 more rows
# ℹ Use `print(n = ...)` to see more rows

 最後にtreat_dfcounter_dfを結合する。treat_dfの行なら"観測値"の値が、counter_dfの行なら"反実仮想"の値が付けられたTypeという列を追加し、このType列をfactor化する。要素の順番は"反実仮想""観測値"の順とする。

tr_ct_df <- bind_rows(list("観測値"   = treat_df, 
                           "反実仮想" = counter_df),
                      .id = "Type") |>
  mutate(Type = factor(Type, levels = c("反実仮想", "観測値")))

 データが揃ったので{ggplot2}を使って折れ線グラフを作ってみよう。

tr_ct_df |>
  # Pref列をリコーディング & factor化
  mutate(Pref = recode(Pref,
                       "Tottori"   = "鳥取",
                       "Shimane"   = "島根",
                       "Tokushima" = "徳島",
                       "Kochi"     = "高知"),
         Pref = factor(Pref, levels = c("鳥取", "島根", "徳島", "高知"))) |>
  ggplot() +
  # x = 23の箇所に垂直線の破線を引く
  geom_vline(xintercept = 23, linetype = 2) +
  geom_line(aes(x = Election, y = Turnout, linetype = Type, color = Type),
            size = 1) +
  # 線の色
  scale_color_manual(values = c("観測値"   = "black", 
                                "反実仮想" = "orange")) +
  # 線のタイプ(1 = 実線; 2 = 破線)
  scale_linetype_manual(values = c("観測値"   = 1, 
                                   "反実仮想" = 2)) +
  # colorの凡例は残し、linetypeの凡例は無くす
  guides(linetype = "none") +
  labs(x = "選挙", y = "投票率 (%)", color = "") +
  # 都道府県ごとにファセット分割
  facet_wrap(~Pref, ncol = 2) +
  theme_bw(base_size = 12) +
  theme(legend.position = "bottom")

 ただし、この図だけだと処置効果(ATT)が分かりにくい。処置効果は観測値(=処置群)と反実仮想(=架空の統制群)の差分である。tr_ct_dfからもATTの計算はできるが、ATTの不確実性(標準誤差や信頼区間など)までは分からない。処置効果(ATT)を抽出はオブジェクト名$est.indで出来るが、3次元配列(array構造)になっているため、オブジェクト名$est.ind[, , "処置群名"]で抽出する必要がある。処置群名は今回の場合だとPrefName_E上の名前である。たとえば、鳥取における処置効果は以下のように抽出する。これらもまた表形式でなく、行列構造であるので、as_tibble()を使って表形式にしよう。

as_tibble(scm_fit$est.ind[, , "Tottori"])
# A tibble: 20 × 5
        Eff  S.E. CI.lower CI.upper  p.value
      <dbl> <dbl>    <dbl>    <dbl>    <dbl>
 1  -2.03    2.33  -6.59       2.53 0.382   
 2  -0.960   2.13  -5.14       3.22 0.653   
 3   1.69    2.07  -2.36       5.74 0.415   
 4   1.07    2.01  -2.87       5.00 0.596   
 5   0.721   2.00  -3.21       4.65 0.719   
 6  -0.704   1.74  -4.11       2.71 0.686   
 7  -0.926   2.33  -5.49       3.63 0.691   
 8   0.816   2.08  -3.27       4.90 0.695   
 9   0.247   1.90  -3.48       3.98 0.897   
10  -3.24    2.23  -7.61       1.13 0.147   
11   5.88    3.03  -0.0512    11.8  0.0520  
12   0.143   2.28  -4.33       4.61 0.950   
13   0.0189  1.89  -3.69       3.73 0.992   
14  -1.48    1.77  -4.94       1.98 0.401   
15  -0.0566  1.83  -3.65       3.53 0.975   
16   0.456   1.79  -3.06       3.97 0.799   
17  -1.63    1.76  -5.07       1.81 0.353   
18  -6.90    2.86 -12.5       -1.30 0.0157  
19  -9.00    3.94 -16.7       -1.28 0.0223  
20 -10.9     3.23 -17.2       -4.56 0.000750

 ここでEff列が架空の鳥取と実際の鳥取の差、つまり処置効果である。そして、CI.lowerCI.upperがそれぞれ95%信頼区間の下限と上限だ。これら4つの県のデータを結合し、Election変数を追加、都道府県をfactor化したものをatt_dfという名で格納する。

att_df <- bind_rows(list("鳥取" = as_tibble(scm_fit$est.ind[, , "Tottori"]),
                         "島根" = as_tibble(scm_fit$est.ind[, , "Shimane"]),
                         "徳島" = as_tibble(scm_fit$est.ind[, , "Tokushima"]),
                         "高知" = as_tibble(scm_fit$est.ind[, , "Kochi"])),
                    .id = "Pref") |>
  mutate(Election = rep(7:26, 4), # 7~26を4回繰り返し
         .before = Pref) |>       # Election列をPref列の前に
  mutate(Pref = fct_inorder(Pref))

att_df
# A tibble: 80 × 7
   Election Pref     Eff  S.E. CI.lower CI.upper p.value
      <int> <fct>  <dbl> <dbl>    <dbl>    <dbl>   <dbl>
 1        7 鳥取  -2.03   2.33    -6.59     2.53   0.382
 2        8 鳥取  -0.960  2.13    -5.14     3.22   0.653
 3        9 鳥取   1.69   2.07    -2.36     5.74   0.415
 4       10 鳥取   1.07   2.01    -2.87     5.00   0.596
 5       11 鳥取   0.721  2.00    -3.21     4.65   0.719
 6       12 鳥取  -0.704  1.74    -4.11     2.71   0.686
 7       13 鳥取  -0.926  2.33    -5.49     3.63   0.691
 8       14 鳥取   0.816  2.08    -3.27     4.90   0.695
 9       15 鳥取   0.247  1.90    -3.48     3.98   0.897
10       16 鳥取  -3.24   2.23    -7.61     1.13   0.147
# … with 70 more rows
# ℹ Use `print(n = ...)` to see more rows

 それでは折れ線グラフを作ってみよう。ただし、今回は点推定値(折れ線グラフ)だけでなく、信頼区間も示す必要があるのでgeom_ribbon()レイヤーを追加する。

att_df |>
  ggplot() +
  geom_hline(yintercept = 0, linetype = 2) +
  geom_vline(xintercept = 23, linetype = 2) +
  # 信頼区間
  geom_ribbon(aes(x = Election, ymin = CI.lower, ymax = CI.upper),
              alpha = 0.25) +
  geom_line(aes(x = Election, y = Eff)) +
  # 以下のgeom_point()なあってもなくても良い
  geom_point(aes(x = Election, y = Eff), 
             size = 3, shape = 21, fill = "black", color = "white") +
  labs(x = "Election", y = "Gap (%p)") +
  facet_wrap(~Pref, ncol = 2) +
  theme_bw(base_size = 12) +
  theme(legend.position = "bottom")

 個別に見ると、合区直後、合区によって投票率を下がったのは鳥取、徳島、高知である。しかし、合区から3回目の選挙となる2022年(第26回)では鳥取のみとなる。

 SCMは統制群から合成された架空の鳥取、島根、徳島、高知を作成し、こちらを実際の統制群として用いる手法である。その際に43都道府県を重み付け合成されるが、それぞれの重みはいくつだろうか。gsynth()から得られたオブジェクトからwgt.impliedを抽出すればその重みが分かる。通常のSCMの場合、重みは正の値であるが、一般化SCMの場合は負も正もあり得る。

scm_fit$wgt.implied
Kochi Shimane Tokushima Tottori
Aichi 0.022 −2.169 −2.755 0.448
Akita −0.151 −0.161 −0.798 0.170
Aomori −0.069 −1.188 −1.613 0.225
Chiba 0.309 −2.096 −1.885 0.327
Ehime 0.238 0.471 0.851 −0.009
Fukui −0.630 4.081 4.711 −1.122
Fukuoka −0.541 0.634 0.342 −0.368
Fukushima 0.947 −0.235 0.669 0.402
Gifu 0.354 0.149 0.566 0.097
Gunma 0.318 2.085 2.968 −0.291
Hiroshima −0.254 1.622 2.146 −0.570
Hokkaido −0.553 −3.254 −5.161 0.643
Hyogo 0.117 −3.045 −3.918 0.723
Ibaraki 0.378 −0.483 0.902 −0.263
Ishikawa −0.480 1.476 1.294 −0.432
Iwate 0.006 −2.615 −4.020 0.833
Kagawa −0.391 3.227 3.960 −0.902
Kagoshima 0.098 0.202 −0.102 0.200
Kanagawa 0.506 −3.426 −3.741 0.828
Kumamoto −0.791 1.767 1.083 −0.490
Kyoto 0.397 −4.141 −5.373 1.200
Mie 0.016 −1.162 −1.447 0.231
Miyagi 0.712 −0.122 1.072 0.070
Miyazaki −1.071 5.630 6.465 −1.699
Nagano 0.498 −1.419 −1.244 0.446
Nagasaki −0.921 1.910 1.572 −0.763
Nara 0.378 −2.473 −2.658 0.589
Nigata 0.856 −2.080 −1.318 0.535
Oita −0.593 2.253 1.749 −0.444
Okayama −0.264 0.966 1.560 −0.557
Okinawa 0.958 0.416 0.697 0.637
Osaka 0.387 −3.404 −4.601 1.118
Saga −0.913 4.668 5.637 −1.552
Saitama 0.217 −0.975 −0.328 −0.031
Shiga −0.012 −0.073 −0.613 0.235
Shizuoka 0.599 0.961 2.654 −0.330
Tochigi 0.409 0.338 1.366 −0.143
Tokyo 0.619 −4.424 −5.420 1.301
Toyama −0.819 4.384 5.389 −1.467
Wakayama −0.050 −1.103 −1.730 0.323
Yamagata 0.652 −1.077 −0.815 0.508
Yamaguchi −1.048 1.821 0.402 −0.382
Yamanashi −0.444 2.064 1.486 −0.274

CausalImpact

 最後にCausalImpactをRで実装した{CausalImpact}パッケージを使ってSCMの例を再現してみよう。この{CausalImpact}パッケージの場合、いくつかデータの制約がある。

  • 処置群の値は1列目 & 2列目以降は統制群の値
  • 処置群は1つのみ
  • 欠損値があってはならない。

 以上の2点である。scm_dfの場合、沖縄県に2つの欠損値(第7・8回参院選)が存在するため、第9回以降のデータのみを使用する。また、共変量などは使わないので選挙(Election)、都道府県名(PrefName_E)、投票率(Turnout)のみを残し、ci_dfという名で格納する。

ci_df <- scm_df |>
  # 第7・8回参院選の場合、沖縄県が欠損しているため、除外
  filter(Election >= 9) |>
  dplyr::select(Election, PrefName_E, Turnout)

ci_df
# A tibble: 846 × 3
   Election PrefName_E Turnout
      <dbl> <chr>        <dbl>
 1        9 Hokkaido      59.7
 2        9 Aomori        53.0
 3        9 Iwate         61.4
 4        9 Miyagi        60.9
 5        9 Akita         67.7
 6        9 Yamagata      70.0
 7        9 Fukushima     74.0
 8        9 Ibaraki       56.5
 9        9 Tochigi       62.0
10        9 Gunma         69.7
# … with 836 more rows
# ℹ Use `print(n = ...)` to see more rows

 続いて、{tidyr}のpivot_wider()関数を使ってci_dfをwide型データへ変換する。pivot_*()関数の使い方については『私たちのR』第15章「整然データ構造」を参照されたい。

ci_df <- ci_df |> 
  pivot_wider(names_from  = PrefName_E,
              values_from = Turnout)

ci_df
# A tibble: 18 × 48
   Election Hokkaido Aomori Iwate Miyagi Akita Yamagata Fukush…¹ Ibaraki Tochigi
      <dbl>    <dbl>  <dbl> <dbl>  <dbl> <dbl>    <dbl>    <dbl>   <dbl>   <dbl>
 1        9     59.7   53.0  61.4   60.9  67.7     70.0     74.0    56.5    62.0
 2       10     75.9   64.4  74.2   72.8  78.0     83.5     84.0    72.2    76.0
 3       11     73.7   66.1  71.9   71.4  77.4     79.9     75.2    62.5    70.2
 4       12     76.3   73.9  75.4   75.2  80.5     82.7     80.8    73.0    73.9
 5       13     57.1   54.0  57.8   53.9  63.7     65.2     68.7    50.7    55.3
 6       14     74.0   74.3  76.4   71.7  81.2     82.5     78.4    71.7    70.7
 7       15     70.9   61.4  70.9   61.2  74.5     75.2     71.7    61.5    62.1
 8       16     59.0   43.9  59.3   48.5  61.6     61.5     59.1    36.6    53.8
 9       17     46.9   46.1  56.4   41.1  57.0     59.1     51.7    36.9    35.9
10       18     59.9   63.1  65.1   54.3  64.1     64.4     65.2    51.0    56.8
11       19     58.5   51.0  66.0   55.6  60.7     63.1     60.9    50.2    53.8
12       20     61.7   53.9  63.3   53.9  65.3     61.7     60.3    50.1    51.0
13       21     62.4   53.9  63.4   55.8  67.7     67.3     61.6    54.0    56.7
14       22     61.9   54.6  60.4   53.3  65.1     64.0     61.6    55.1    56.6
15       23     54.4   46.2  57.5   50.7  56.2     60.8     54.5    49.7    49.7
16       24     56.8   55.3  57.8   52.4  60.9     62.2     57.1    50.8    51.4
17       25     53.8   42.9  56.5   51.2  56.3     60.7     52.4    45.0    44.1
18       26     54.0   49.5  55.4   48.8  55.6     61.9     53.4    47.2    47.0
# … with 38 more variables: Gunma <dbl>, Saitama <dbl>, Chiba <dbl>,
#   Tokyo <dbl>, Kanagawa <dbl>, Nigata <dbl>, Toyama <dbl>, Ishikawa <dbl>,
#   Fukui <dbl>, Yamanashi <dbl>, Nagano <dbl>, Gifu <dbl>, Shizuoka <dbl>,
#   Aichi <dbl>, Mie <dbl>, Shiga <dbl>, Kyoto <dbl>, Osaka <dbl>, Hyogo <dbl>,
#   Nara <dbl>, Wakayama <dbl>, Tottori <dbl>, Shimane <dbl>, Okayama <dbl>,
#   Hiroshima <dbl>, Yamaguchi <dbl>, Tokushima <dbl>, Kagawa <dbl>,
#   Ehime <dbl>, Kochi <dbl>, Fukuoka <dbl>, Saga <dbl>, Nagasaki <dbl>, …
# ℹ Use `colnames()` to see all variable names

 {CausalImpact}の場合、処置群は1つのみとなる。ci_dfには4県の処置群があるが、ここでは鳥取(PrefName_E == "Tottori")におけるATTを推定してみよう。relocate()関数を使って、Tottori列をElection列より前へ移動させる。続いて、不要になったElection関数とその他の統制群列を除外する。

ci_df <- ci_df |> 
  # Tottori列を最先端へ(= Election列より前)
  relocate(Tottori, .before = Election) |>
  # その他の処置群とElection列を除外
  dplyr::select(-c(Election, Shimane, Tokushima, Kochi))

 それでは実装してみよう。使う関数はCausalImpact()関数だ。第1引数はデータである。他に必要な引数としてはpre.periodpost.periodがあり、それぞれ統制の時期と初期が行われた時期を意味する。処置前の時期というのはci_dfでいうと、第9回から第23回参院選までだ。ci_dfの1行目が第9回参院選、15行目は第23回参院選である。16行目から18行目は処置が行われてからの時期だ。それぞれpre.periodpost.periodで指定する。

ci_fit <- ci_df |>
  CausalImpact(pre.period = c(1, 15), post.period = c(16, 18))

ci_fit
Posterior inference {CausalImpact}

                         Average         Cumulative   
Actual                   52              155          
Prediction (s.d.)        61 (3.6)        183 (10.9)   
95% CI                   [54, 68]        [162, 204]   
                                                      
Absolute effect (s.d.)   -9.4 (3.6)      -28.1 (10.9) 
95% CI                   [-16, -2.2]     [-49, -6.5]  
                                                      
Relative effect (s.d.)   -15% (5.9%)     -15% (5.9%)  
95% CI                   [-27%, -3.6%]   [-27%, -3.6%]

Posterior tail-area probability p:   0.01308
Posterior prob. of a causal effect:  98.692%

For more details, type: summary(impact, "report")

 推定されたATTはAbsolute effect行のAverageから確認できる。これは鳥取県の場合、合区によって投票率が約-9.364%p低下したことを意味する。ちなみにSCMで指定した鳥取県のATTはいくらだろうか。

tottori_att <- scm_fit |>
  cumuEff(cumu   = FALSE,     # TRUEの場合、累積ATT
          id     = "Tottori", # PrefName_E上の名前
          period = c(1, 3))   # 処置が行われた時期を1する。

tottori_att
$catt
[1]  -6.903364  -8.997478 -10.885151

$est.catt
        CATT     S.E.  CI.lower  CI.upper p.value
1  -6.903364 2.856475 -12.35897 -1.634203   0.020
2  -8.997478 3.937604 -16.08604 -1.068505   0.032
3 -10.885151 3.229445 -16.71687 -4.236510   0.000

 SCMから得られた結果、第24回におけるATTは-6.903、第25回は-8.997、第26回は-10.885であり、平均値は-8.929である。0.435%pの差はあるが、それなりに近い結果が得られている。

 分析の手法はことなるものの、SCMもCausalImpactも、架空の統制群を作るという点では共通している。したがって、SCMで示したような図も作成できよう。{CausalImpact}もplot()関数で作図ができる。

plot(ci_fit)

 上の図が架空の鳥取県と実際の鳥取県を示したもの、真ん中の図はその2つの差分であり、ATTを意味する。下の図はATTの累積である。

 それでは既に作成したSCMの図と比較してみよう。

 CausalImpactの方が信頼区間が広いものの、傾向としてはほとんど同じである。今回は{CausalImpact}の基本的な使い方のみを紹介したが、{bsts}パッケージの使用経験がある場合、ある程度のモデリングもできるため、より精度の良い推定ができるかも知れない。宋はCausalImpactの便利さ(=統制群すら必要としない点、{CausalImpact}パッケージの使いやすさ)は認めるものの、まだまだ改善の余地があると考えているので、ここまでに紹介させて頂きたい。詳細な紹介についてはill-identifiedさんの記事コードが詳しい。

  1. fixed_effects引数でなく、回帰式に説明変数として指定しても結果は同じである。しかし、回帰式に書く場合、固定効果の推定値も全て出力され、推定結果が非常に長くなる。しかし、固定効果の推定値は論文内で報告することもない。fixed_effectsで指定すると、それらの結果は省略される。↩︎

  2. {gsynth}は通常のSCMでなく、一般化SCM(Generalized Synthetic Control Method)を実装したパッケージである。通常のSCMは処置群は1つのみに使えるが、一般化SCMは処置群が複数あっても使える。一般化SCMについてはYiqing Xu. 2017. “Generalized Synthetic Control Method: Causal Inference with Interactive Fixed Effects Models,” Political Analysis, 25 (1): 57-76.を参照されたい。↩︎

  3. 第7〜23回参院選のデータは参議院議員通常選挙データベースからダウンロードしたものを架空したものであり、第24〜26回データは宋が手入力したものである。↩︎