差分の差分法

スライド

新しいタブで開く

セットアップ

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

pacman::p_load(tidyverse,     # Rの必須パッケージ
               summarytools,  # 記述統計
               modelsummary,  # 推定結果の要約
               estimatr)      # ロバストな回帰分析

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

did_df
# A tibble: 18,749 × 14
   county state  year shooting fatal_shooting non_fatal_shooting turnout demvote
    <dbl> <chr> <dbl>    <dbl>          <dbl>              <dbl>   <dbl>   <dbl>
 1   1001 01     1996        0              0                  0    56.1    32.5
 2   1001 01     2000        0              0                  0    56.8    28.7
 3   1001 01     2004        0              0                  0    60.2    23.7
 4   1001 01     2008        0              0                  0    64.1    25.8
 5   1001 01     2012        0              0                  0    61.0    26.5
 6   1001 01     2016        0              0                  0    60.7    24.0
 7   1003 01     1996        0              0                  0    51.8    27.1
 8   1003 01     2000        0              0                  0    54.6    24.8
 9   1003 01     2004        0              0                  0    59.8    22.5
10   1003 01     2008        0              0                  0    62.4    23.8
# ℹ 18,739 more rows
# ℹ 6 more variables: population <dbl>, non_white <dbl>,
#   change_unem_rate <dbl>, county_f <dbl>, state_f <chr>, year_f <dbl>

 データの詳細はスライドを参照すること。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_shooting non_fatal_shooting turnout demvote
    <dbl> <chr> <dbl>    <dbl>          <dbl>              <dbl>   <dbl>   <dbl>
 1   1001 01     1996        0              0                  0    56.1    32.5
 2   1001 01     2000        0              0                  0    56.8    28.7
 3   1001 01     2004        0              0                  0    60.2    23.7
 4   1001 01     2008        0              0                  0    64.1    25.8
 5   1001 01     2012        0              0                  0    61.0    26.5
 6   1001 01     2016        0              0                  0    60.7    24.0
 7   1003 01     1996        0              0                  0    51.8    27.1
 8   1003 01     2000        0              0                  0    54.6    24.8
 9   1003 01     2004        0              0                  0    59.8    22.5
10   1003 01     2008        0              0                  0    62.4    23.8
# ℹ 18,739 more rows
# ℹ 6 more variables: population <dbl>, non_white <dbl>,
#   change_unem_rate <dbl>, county_f <fct>, state_f <fct>, year_f <fct>

 連続変数(shootingからchange_unem_rateまで)の記述統計量を出力する。ここで一つ注意が必要だ。それはselect()関数の使い方である。具体的に言えば、使い方そのものは変わらない。

did_df |>
  select(shooting:change_unem_rate) |>
  descr(stats = c("mean", "sd", "min", "max", "n.valid"),
        transpose = TRUE, order = "p")
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))
 (1)   (2)   (3)   (4)   (5)   (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")

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