TWFE
それでは二元配置固定効果（Two-way fixed effect; TWFE ）モデルで差分の差分法を実装する方法について紹介する。推定式は以下の通りである。
\[
\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}
\]
\(i\) ：カウンティ\(t\) ：年\(\mbox{Otucome}\) ：応答変数
turnout：投票率（大統領選挙）demvote：民主党候補者の得票率 \(\mbox{Shooting}\) ：処置変数
shooting：銃撃事件の発生有無fatal_shooting：死者を伴う銃撃事件の発生有無non_fatal_shooting：死者を伴わない銃撃事件の発生有無 \(\mbox{Controls}_k\) ：統制変数
\(\mbox{Controls}_1\) （population）：カウンティーの人口\(\mbox{Controls}_2\) （non_white）：非白人の割合\(\mbox{Controls}_3\) （change_unem_rate）：失業率の変化統制変数あり/なしのモデルを個別に推定
\(\lambda\) ：年固定効果\(\omega\) ：カウンティー固定効果
今回は応答変数が2種類、処置変数が3種類、共変量の有無でモデルを分けるので、推定するモデルは計12個である。
モデル1
did_fit1turnoutshooting-
モデル2
did_fit2turnoutshooting\(\bigcirc\)
モデル3
did_fit3turnoutfatal_shooting-
モデル4
did_fit4turnoutfatal_shooting\(\bigcirc\)
モデル5
did_fit5turnoutnon_fatal_shooting-
モデル6
did_fit6turnoutnon_fatal_shooting\(\bigcirc\)
モデル7
did_fit7demvoteshooting-
モデル8
did_fit8demvoteshooting\(\bigcirc\)
モデル9
did_fit9demvotefatal_shooting-
モデル10
did_fit10demvotefatal_shooting\(\bigcirc\)
モデル11
did_fit11demvotenon_fatal_shooting-
モデル12
did_fit12demvotenon_fatal_shooting\(\bigcirc\)
まずはモデル1を推定し、did_fit1という名のオブジェクトとして格納する。基本的には線形回帰分析であるため、lm()でも推定はできる。しかし、差分の差分法の場合、通常、クラスター化した頑健な標準誤差（cluster robust standard error）を使う。lm()単体ではこれが計算できないため、今回は{fixest}パッケージが提供するfeols()関数を使用する。使い方はlm()に似ている。まず、回帰式を書く。たとえば、モデル1の場合、trunout ~ shootingである。ただし、feols()で固定効果を投入するには、回帰式の右側に|を追加し、その右側に固定効果として投入する変数を指定する必要がある。たとえば、county_f変数とyear_f変数が固定効果として投入される場合、trunout ~ shooting | county_f + year_fが回帰式である。最後に標準誤差をクラスタリングする単位はcluster引数で指定する。たとえば、州（state_f）単位で標準誤差をクラスタリングしたい時はcluster = ~ state_fと書く。
本研究のデータは1行がCounty \(\times\) Yearを意味するパネルデータだが、クラスター標準誤差を推定する際、カウンティでなく、州でクラスタリングしている。これについては明確な理由は述べられていない。もしかしたら、同一州内のカウンティ間の政治的トレンドの相関が理論的に懸念されたかも知れない。また、より保守的にアプローチし、より粗い方でクラスタイングをした可能性もある（Abadie et al. 2023）。
そもそも、本研究は政治学トップジャーナルであるAmerican Political Science Review（APSR ）に掲載された論文であるが、APSRへの掲載が本論文の方法論的正しさを保証するわけではない（むろん、トップジャーナルであれば相対的にクオリティーの良い論文が見つかりやすいだろう）。実際、本研究に対して方法論的に辛辣な批判をした研究もあり、彼らの研究によると、García-Montoya et al. （2022）の結果は過大評価されており、実際の因果効果はほぼゼロだという（Hassell and Holbein 2025）。また、Hassell and Holbein（2025）では州単位でクラスタリングした標準誤差が使用されており、州単位でのクラスタリングが適切だった可能性もある。。
<- feols (turnout ~ shooting | county_f + year_f, data = did_df, cluster = ~ state_f)summary (did_fit1)
OLS estimation, Dep. Var.: turnout
Observations: 18,620
Fixed-effects: county_f: 3,104, year_f: 6
Standard-errors: Clustered (state_f)
Estimate Std. Error t value Pr(>|t|)
shooting -0.521084 0.589288 -0.88426 0.38087
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
RMSE: 3.89827 Adj. R2: 0.828867
Within R2: 6.907e-5
処置効果の推定値は-0.521である。これは学校内銃撃事件が発生したカウンティーの場合、大統領選挙において投票率が約-0.521パーセンテージポイント（pp）低下することを意味する。しかし、標準誤差がかなり大きく、統計的有意な結果ではない。つまり、「学校内銃撃事件が投票率を上げる（or 下げる）とは言えない」と解釈できる。決して「学校内銃撃事件が投票率を上げない（or 下げない）」と解釈しないこと。
共変量を投入してみたらどうだろうか。たとえば、人口は自治体の都市化程度を表すこともあるので、都市化程度と投票率には関係があると考えられる。また、人口が多いと自然に事件が発生する確率もあがるので、交絡要因として考えられる。人種や失業率も同様であろう。ここではカウンティーの人口（population）、非白人の割合（non_white）、失業率の変化（change_unem_rate）を統制変数として投入し、did_fit2という名で格納する。統制変数は|の左側に+演算して加える。
<- feols (turnout ~ shooting + + non_white + change_unem_rate | county_f + year_f, data = did_df, cluster = ~ state_f)summary (did_fit2)
OLS estimation, Dep. Var.: turnout
Observations: 18,620
Fixed-effects: county_f: 3,104, year_f: 6
Standard-errors: Clustered (state_f)
Estimate Std. Error t value Pr(>|t|)
shooting -0.70979224 0.56925534 -1.24688 0.218369
population 0.00000803 0.00000519 1.54653 0.128411
non_white -34.83442620 14.21507415 -2.45053 0.017882 *
change_unem_rate 0.15920199 0.06009365 2.64923 0.010828 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
RMSE: 3.86574 Adj. R2: 0.831678
Within R2: 0.016687
処置効果の推定値は-0.710である。これは他の条件が同じ場合、学校内銃撃事件が発生したカウンティーは大統領選挙において投票率が約-0.710pp低下することを意味する。今回も統計的に有意な因果効果は確認されなかった。
これまでの処置変数は死者の有無と関係なく、学校内銃撃事件が発生したか否かだった。もしかしたら、死者を伴う銃撃事件が発生した場合、その効果が大きいかも知れない。したがって、これからは処置変数を死者を伴う学校内銃撃事件の発生有無（fatal_shooting）、死者を伴わない学校内銃撃事件の発生有無（non_fatal_shooting）に変えてもう一度推定してみよう。
<- feols (turnout ~ fatal_shooting | county_f + year_f, data = did_df, cluster = ~ state_f)<- feols (turnout ~ fatal_shooting + + non_white + change_unem_rate | county_f + year_f, data = did_df, cluster = ~ state_f)<- feols (turnout ~ non_fatal_shooting | county_f + year_f, data = did_df, cluster = ~ state_f)<- feols (turnout ~ non_fatal_shooting + + non_white + change_unem_rate | county_f + year_f, data = did_df, cluster = ~ state_f)
これまで推定してきた6つのモデルを比較してみよう。
modelsummary (list (did_fit1, did_fit2, did_fit3,
(1)
(2)
(3)
(4)
(5)
(6)
shooting
-0.521
-0.710
(0.589)
(0.569)
population
0.000
0.000
0.000
(0.000)
(0.000)
(0.000)
non_white
-34.834
-34.844
-34.814
(14.215)
(14.254)
(14.251)
change_unem_rate
0.159
0.159
0.160
(0.060)
(0.060)
(0.060)
fatal_shooting
-0.678
-0.918
(0.565)
(0.600)
non_fatal_shooting
-0.239
-0.327
(0.998)
(1.045)
Num.Obs.
18620
18620
18620
18620
18620
18620
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
R2 Within
0.000
0.017
0.000
0.017
0.000
0.017
R2 Within Adj.
0.000
0.016
0.000
0.016
-0.000
0.016
AIC
109727.5
109421.4
109727.4
109421.3
109728.6
109423.6
BIC
134085.0
133802.4
134084.9
133802.3
134086.1
133804.6
RMSE
3.90
3.87
3.90
3.87
3.90
3.87
FE: county_f
X
X
X
X
X
X
FE: year_f
X
X
X
X
X
X
いずれのモデルも統計的に有意な処置効果は確認されていない。これらの結果を表として報告するには紙がもったいない気もする。これらの結果はOnline Appendixに回し、本文中には処置効果の点推定値と95%信頼区間を示せば良いだろう。
{broom}のtidy()関数で推定結果のみを抽出し、それぞれオブジェクトとして格納しておこう。
<- tidy (did_fit1, conf.int = TRUE )<- tidy (did_fit2, conf.int = TRUE )<- tidy (did_fit3, conf.int = TRUE )<- tidy (did_fit4, conf.int = TRUE )<- tidy (did_fit5, conf.int = TRUE )<- tidy (did_fit6, conf.int = TRUE )
全て確認する必要はないので、tidy_fit1のみを確認してみる。
# A tibble: 1 × 7
term estimate std.error statistic p.value conf.low conf.high
<chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 shooting -0.521 0.589 -0.884 0.381 -1.71 0.663
以上の6つの表形式オブジェクトを一つの表としてまとめる。それぞれのオブジェクトには共変量の有無_処置変数の種類の名前を付けよう。共変量なしのモデルはM1、ありのモデルはM2とする。処置変数はshootingの場合はTr1、fatal_shootingはTr2、non_fatal_shootingはTr3とする。
<- 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" )
# A tibble: 15 × 8
Model term estimate std.error statistic p.value conf.low conf.high
<chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 M1_Tr1 shooting -0.521 0.589 -0.884 0.381 -1.71 0.663
2 M2_Tr1 shooting -0.710 0.569 -1.25 0.218 -1.85 0.434
3 M2_Tr1 population 0.00000803 0.00000519 1.55 0.128 -0.00000240 0.0000185
4 M2_Tr1 non_white -34.8 14.2 -2.45 0.0179 -63.4 -6.27
5 M2_Tr1 change_unem_rate 0.159 0.0601 2.65 0.0108 0.0384 0.280
6 M1_Tr2 fatal_shooting -0.678 0.565 -1.20 0.236 -1.81 0.458
7 M2_Tr2 fatal_shooting -0.918 0.600 -1.53 0.133 -2.12 0.288
8 M2_Tr2 population 0.00000803 0.00000525 1.53 0.132 -0.00000251 0.0000186
9 M2_Tr2 non_white -34.8 14.3 -2.44 0.0181 -63.5 -6.20
10 M2_Tr2 change_unem_rate 0.159 0.0602 2.65 0.0109 0.0382 0.280
11 M1_Tr3 non_fatal_shooting -0.239 0.998 -0.239 0.812 -2.24 1.77
12 M2_Tr3 non_fatal_shooting -0.327 1.05 -0.313 0.756 -2.43 1.77
13 M2_Tr3 population 0.00000782 0.00000521 1.50 0.140 -0.00000266 0.0000183
14 M2_Tr3 non_white -34.8 14.3 -2.44 0.0182 -63.5 -6.18
15 M2_Tr3 change_unem_rate 0.160 0.0601 2.65 0.0107 0.0387 0.280
続いて、処置効果のみを抽出する。処置効果はterm列の値が"shooting"、"fatal_shooting"、"non_fatal_shooting"のいずれかと一致する行であるため、filter()関数を使用する。
<- did_est1 |> filter (term %in% c ("shooting" , "fatal_shooting" , "non_fatal_shooting" ))
# A tibble: 6 × 8
Model term estimate std.error statistic p.value conf.low conf.high
<chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 M1_Tr1 shooting -0.521 0.589 -0.884 0.381 -1.71 0.663
2 M2_Tr1 shooting -0.710 0.569 -1.25 0.218 -1.85 0.434
3 M1_Tr2 fatal_shooting -0.678 0.565 -1.20 0.236 -1.81 0.458
4 M2_Tr2 fatal_shooting -0.918 0.600 -1.53 0.133 -2.12 0.288
5 M1_Tr3 non_fatal_shooting -0.239 0.998 -0.239 0.812 -2.24 1.77
6 M2_Tr3 non_fatal_shooting -0.327 1.05 -0.313 0.756 -2.43 1.77
ちなみにgrepl()関数を使うと、"shooting"が含まれる行を抽出することもできる。以下のコードは上記のコードと同じ機能をする。
<- did_est1 |> filter (grepl ("shooting" , term))
つづいて、Model列をModelとTreat列へ分割する。
separate()関数は列を任意の文字列を基準に複数の列へ分割する関数だ。{tidyr}に含まれている関数だが、{tidyr}は{tidyverse}を構成する一部だから、{tidyverse}を読み込んでおいたのであれば、別途読み込む必要はない。
必要な必須引数は3つ、col、into、sepだ。colは分割する列名、into分割後に列名（c()でベクトルとして指定する）、sepが分割の基準となる文字列だ。以下の例ではModel列を「-」（ハイフン）を基準にModel、Treatの2列に分割するコードだ。
|> separate (col = Model,into = c ("Model" , "Treat" ),sep = "-" )
Error:
! object 'データ' not found
これを利用すれば、2026-08-19のような値が格納されている列を年、月、日に分割できる。
<- did_est1 |> separate (col = Model,into = c ("Model" , "Treat" ),sep = "_" )
# A tibble: 6 × 9
Model Treat term estimate std.error statistic p.value conf.low conf.high
<chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 M1 Tr1 shooting -0.521 0.589 -0.884 0.381 -1.71 0.663
2 M2 Tr1 shooting -0.710 0.569 -1.25 0.218 -1.85 0.434
3 M1 Tr2 fatal_shooting -0.678 0.565 -1.20 0.236 -1.81 0.458
4 M2 Tr2 fatal_shooting -0.918 0.600 -1.53 0.133 -2.12 0.288
5 M1 Tr3 non_fatal_shooting -0.239 0.998 -0.239 0.812 -2.24 1.77
6 M2 Tr3 non_fatal_shooting -0.327 1.05 -0.313 0.756 -2.43 1.77
可視化に入る前にModel列とTreat列の値を修正する。Model列の値が"M1"なら"County-Year FE"に、それ以外なら"County-Year FE + Covariates"とリコーディングする。戻り値が2種類だからif_else()を使う。Treat列の場合、戻り値が3つなので、recode()かcase_when()を使う。ここではrecode()を使ってリコーディングする。最後にModelとTreatを表示順番でfactor化し（fct_inorder()）、更に順番を逆転する（fct_rev()）。
fct_rev()関数はfactor型の要素の順番を逆転する関数だ。{forcats}に含まれている関数だが、{forcats}は{tidyverse}を構成する一部だから、{tidyverse}を読み込んでおいたのであれば、別途読み込む必要はない。
<- 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)))
# A tibble: 6 × 9
Model Treat term estimate std.error statistic p.value conf.low conf.high
<fct> <fct> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 County-Year FE Any Sho… shoo… -0.521 0.589 -0.884 0.381 -1.71 0.663
2 County-Year FE + Covariates Any Sho… shoo… -0.710 0.569 -1.25 0.218 -1.85 0.434
3 County-Year FE Fatal S… fata… -0.678 0.565 -1.20 0.236 -1.81 0.458
4 County-Year FE + Covariates Fatal S… fata… -0.918 0.600 -1.53 0.133 -2.12 0.288
5 County-Year FE Nonfata… non_… -0.239 0.998 -0.239 0.812 -2.24 1.77
6 County-Year FE + Covariates Nonfata… non_… -0.327 1.05 -0.313 0.756 -2.43 1.77
それでは{ggplot2}を使ってpointrangeプロットを作成してみよう。
|> ggplot () + # x = 0の箇所に垂直線を引く。垂直線は破線（dashed）とする。 geom_vline (xintercept = 0 , linetype = "dashed" ) + 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 (pp)" , 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_classic (base_size = 12 ) + theme (legend.position = "bottom" )
元の論文を見ると、点の上に点推定値が書かれているが、私たちもこれを真似してみよう。文字列をプロットするレイヤーはgeom_text()とgeom_label()、annotate()があるが、ここではgeom_text()を使用する。文字列が表示される横軸上の位置（x）と縦軸上の位置（y）、そして出力する文字列（label）をマッピングする。点推定値は3桁まで出力したいので、sprintf()を使って、3桁に丸める。ただし、これだけだと点と文字が重なってしまう。vjustを-0.75にすることで、出力する文字列を点の位置を上の方向へ若干ずらすことができる。
|> 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 (pp)" , y = "" , color = "" ) + scale_color_manual (values = c ("County-Year FE" = "black" , "County-Year FE + Covariates" = "gray50" )) + coord_cartesian (xlim = c (- 10 , 10 )) + theme_classic (base_size = 12 ) + theme (legend.position = "bottom" )
ちなみにこのコードを見ると、geom_pointrange()とgeom_text()はx、y、colorを共有しているので、ggplot()内でマッピングすることもできる。
|> 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 (pp)" , y = "" , color = "" ) + scale_color_manual (values = c ("County-Year FE" = "black" , "County-Year FE + Covariates" = "gray50" )) + coord_cartesian (xlim = c (- 10 , 10 )) + theme_classic (base_size = 12 ) + theme (legend.position = "bottom" )
続いて、民主党候補者の得票率（demvote）を応答変数として6つのモデルを推定し、同じ作業を繰り返す。
<- feols (demvote ~ shooting | county_f + year_f, data = did_df, cluster = ~ state_f)<- feols (demvote ~ shooting + + non_white + change_unem_rate | county_f + year_f, data = did_df, cluster = ~ state_f)<- feols (demvote ~ fatal_shooting | county_f + year_f, data = did_df, cluster = ~ state_f)<- feols (demvote ~ fatal_shooting + + non_white + change_unem_rate | county_f + year_f, data = did_df, cluster = ~ state_f)<- feols (demvote ~ non_fatal_shooting | county_f + year_f, data = did_df, cluster = ~ state_f)<- feols (demvote ~ non_fatal_shooting + + non_white + change_unem_rate | county_f + year_f, data = did_df, cluster = ~ state_f)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.799)
(0.672)
population
0.000
0.000
0.000
(0.000)
(0.000)
(0.000)
non_white
86.873
86.866
86.796
(16.303)
(16.296)
(16.296)
change_unem_rate
-0.139
-0.139
-0.140
(0.116)
(0.116)
(0.116)
fatal_shooting
4.404
1.782
(0.997)
(0.763)
non_fatal_shooting
4.217
2.930
(1.185)
(0.947)
Num.Obs.
18620
18620
18620
18620
18620
18620
R2
0.882
0.894
0.882
0.894
0.882
0.894
R2 Adj.
0.859
0.873
0.859
0.873
0.858
0.873
R2 Within
0.004
0.105
0.002
0.104
0.001
0.104
R2 Within Adj.
0.003
0.105
0.002
0.104
0.001
0.104
AIC
116926.6
114937.5
116953.4
114950.3
116967.9
114944.2
BIC
141284.1
139318.5
141310.9
139331.3
141325.3
139325.2
RMSE
4.73
4.48
4.73
4.48
4.73
4.48
FE: county_f
X
X
X
X
X
X
FE: year_f
X
X
X
X
X
X
今回はいずれも統計的に有意な結果が得られている。例えば、モデル7（did_fit7）の場合、処置効果の推定値は4.513である。これは学校内銃撃事件が発生したカウンティーの場合、大統領選挙において民主党候補者の得票率が約4.513pp増加することを意味する。
以上の結果を図としてまとめてみよう。これまでのコードと同じ手順で作図するため、説明は割愛する。
<- tidy (did_fit7, conf.int = TRUE )<- tidy (did_fit8, conf.int = TRUE )<- tidy (did_fit9, conf.int = TRUE )<- tidy (did_fit10, conf.int = TRUE )<- tidy (did_fit11, conf.int = TRUE )<- tidy (did_fit12, conf.int = TRUE )<- 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" )
# A tibble: 15 × 8
Model term estimate std.error statistic p.value conf.low conf.high
<chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 M1_Tr1 shooting 4.51 0.799 5.65 0.000000807 2.91 6.12
2 M2_Tr1 shooting 2.36 0.672 3.52 0.000956 1.01 3.72
3 M2_Tr1 population 0.0000317 0.00000619 5.13 0.00000493 0.0000193 0.0000442
4 M2_Tr1 non_white 86.9 16.3 5.33 0.00000248 54.1 120.
5 M2_Tr1 change_unem_rate -0.139 0.116 -1.19 0.238 -0.373 0.0948
6 M1_Tr2 fatal_shooting 4.40 0.997 4.42 0.0000548 2.40 6.41
7 M2_Tr2 fatal_shooting 1.78 0.763 2.34 0.0237 0.248 3.32
8 M2_Tr2 population 0.0000321 0.00000621 5.17 0.00000436 0.0000196 0.0000445
9 M2_Tr2 non_white 86.9 16.3 5.33 0.00000247 54.1 120.
10 M2_Tr2 change_unem_rate -0.139 0.116 -1.20 0.237 -0.373 0.0946
11 M1_Tr3 non_fatal_shooting 4.22 1.18 3.56 0.000838 1.84 6.60
12 M2_Tr3 non_fatal_shooting 2.93 0.947 3.10 0.00325 1.03 4.83
13 M2_Tr3 population 0.0000323 0.00000613 5.27 0.00000308 0.0000200 0.0000446
14 M2_Tr3 non_white 86.8 16.3 5.33 0.00000250 54.0 120.
15 M2_Tr3 change_unem_rate -0.140 0.116 -1.20 0.235 -0.374 0.0938
<- did_est2 |> filter (grepl ("shooting" , term))
# A tibble: 6 × 8
Model term estimate std.error statistic p.value conf.low conf.high
<chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 M1_Tr1 shooting 4.51 0.799 5.65 0.000000807 2.91 6.12
2 M2_Tr1 shooting 2.36 0.672 3.52 0.000956 1.01 3.72
3 M1_Tr2 fatal_shooting 4.40 0.997 4.42 0.0000548 2.40 6.41
4 M2_Tr2 fatal_shooting 1.78 0.763 2.34 0.0237 0.248 3.32
5 M1_Tr3 non_fatal_shooting 4.22 1.18 3.56 0.000838 1.84 6.60
6 M2_Tr3 non_fatal_shooting 2.93 0.947 3.10 0.00325 1.03 4.83
<- did_est2 |> separate (col = Model,into = c ("Model" , "Treat" ),sep = "_" )
# A tibble: 6 × 9
Model Treat term estimate std.error statistic p.value conf.low conf.high
<chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 M1 Tr1 shooting 4.51 0.799 5.65 0.000000807 2.91 6.12
2 M2 Tr1 shooting 2.36 0.672 3.52 0.000956 1.01 3.72
3 M1 Tr2 fatal_shooting 4.40 0.997 4.42 0.0000548 2.40 6.41
4 M2 Tr2 fatal_shooting 1.78 0.763 2.34 0.0237 0.248 3.32
5 M1 Tr3 non_fatal_shooting 4.22 1.18 3.56 0.000838 1.84 6.60
6 M2 Tr3 non_fatal_shooting 2.93 0.947 3.10 0.00325 1.03 4.83
<- 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)))
# A tibble: 6 × 9
Model Treat term estimate std.error statistic p.value conf.low conf.high
<fct> <fct> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 County-Year FE Any Sho… shoo… 4.51 0.799 5.65 8.07e-7 2.91 6.12
2 County-Year FE + Covariates Any Sho… shoo… 2.36 0.672 3.52 9.56e-4 1.01 3.72
3 County-Year FE Fatal S… fata… 4.40 0.997 4.42 5.48e-5 2.40 6.41
4 County-Year FE + Covariates Fatal S… fata… 1.78 0.763 2.34 2.37e-2 0.248 3.32
5 County-Year FE Nonfata… non_… 4.22 1.18 3.56 8.38e-4 1.84 6.60
6 County-Year FE + Covariates Nonfata… non_… 2.93 0.947 3.10 3.25e-3 1.03 4.83
|> 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 (pp)" , y = "" , color = "" ) + scale_color_manual (values = c ("County-Year FE" = "black" , "County-Year FE + Covariates" = "gray50" )) + coord_cartesian (xlim = c (- 10 , 10 )) + theme_classic (base_size = 12 ) + theme (legend.position = "bottom" )
最後に、これまで作成した2つの図を一つにまとめてみよう。bind_rows()関数を使い、それぞれの表に識別子（Outcome）を与える。
<- bind_rows (list ("Out1" = did_est1,"Out2" = did_est2),.id = "Outcome" )
# A tibble: 12 × 10
Outcome Model Treat term estimate std.error statistic p.value conf.low conf.high
<chr> <fct> <fct> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 Out1 County-Year FE Any … shoo… -0.521 0.589 -0.884 3.81e-1 -1.71 0.663
2 Out1 County-Year FE + Cov… Any … shoo… -0.710 0.569 -1.25 2.18e-1 -1.85 0.434
3 Out1 County-Year FE Fata… fata… -0.678 0.565 -1.20 2.36e-1 -1.81 0.458
4 Out1 County-Year FE + Cov… Fata… fata… -0.918 0.600 -1.53 1.33e-1 -2.12 0.288
5 Out1 County-Year FE Nonf… non_… -0.239 0.998 -0.239 8.12e-1 -2.24 1.77
6 Out1 County-Year FE + Cov… Nonf… non_… -0.327 1.05 -0.313 7.56e-1 -2.43 1.77
7 Out2 County-Year FE Any … shoo… 4.51 0.799 5.65 8.07e-7 2.91 6.12
8 Out2 County-Year FE + Cov… Any … shoo… 2.36 0.672 3.52 9.56e-4 1.01 3.72
9 Out2 County-Year FE Fata… fata… 4.40 0.997 4.42 5.48e-5 2.40 6.41
10 Out2 County-Year FE + Cov… Fata… fata… 1.78 0.763 2.34 2.37e-2 0.248 3.32
11 Out2 County-Year FE Nonf… non_… 4.22 1.18 3.56 8.38e-4 1.84 6.60
12 Out2 County-Year FE + Cov… Nonf… non_… 2.93 0.947 3.10 3.25e-3 1.03 4.83
Outcome列のリコーディングし、factor化する。
<- did_est |> mutate (Outcome = if_else (Outcome == "Out1" ,"Change in Turnout (pp)" ,"Change in Democratic Vote Share (pp)" ),Outcome = fct_inorder (Outcome))
# A tibble: 12 × 10
Outcome Model Treat term estimate std.error statistic p.value conf.low conf.high
<fct> <fct> <fct> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 Change in Turnout (pp) Coun… Any … shoo… -0.521 0.589 -0.884 3.81e-1 -1.71 0.663
2 Change in Turnout (pp) Coun… Any … shoo… -0.710 0.569 -1.25 2.18e-1 -1.85 0.434
3 Change in Turnout (pp) Coun… Fata… fata… -0.678 0.565 -1.20 2.36e-1 -1.81 0.458
4 Change in Turnout (pp) Coun… Fata… fata… -0.918 0.600 -1.53 1.33e-1 -2.12 0.288
5 Change in Turnout (pp) Coun… Nonf… non_… -0.239 0.998 -0.239 8.12e-1 -2.24 1.77
6 Change in Turnout (pp) Coun… Nonf… non_… -0.327 1.05 -0.313 7.56e-1 -2.43 1.77
7 Change in Democratic V… Coun… Any … shoo… 4.51 0.799 5.65 8.07e-7 2.91 6.12
8 Change in Democratic V… Coun… Any … shoo… 2.36 0.672 3.52 9.56e-4 1.01 3.72
9 Change in Democratic V… Coun… Fata… fata… 4.40 0.997 4.42 5.48e-5 2.40 6.41
10 Change in Democratic V… Coun… Fata… fata… 1.78 0.763 2.34 2.37e-2 0.248 3.32
11 Change in Democratic V… Coun… Nonf… non_… 4.22 1.18 3.56 8.38e-4 1.84 6.60
12 Change in Democratic V… Coun… Nonf… non_… 2.93 0.947 3.10 3.25e-3 1.03 4.83
図の作り方はこれまでと変わらないが、ファセット分割を行うため、facet_wrap()レイヤーを追加する。
|> 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" )
以上の結果から「学校内銃撃事件の発生は投票参加を促すとは言えないものの、民主党候補者の得票率を上げる」ということが言えよう。
Staggered DID
以上がGarcía-Montoya, Arjona, and Lacombe（2022） の主要結果の再現であるが、政治学トップジャーナルであるAPSRに掲載されたことが、本研究が完全無欠ということを意味するわけではない。なぜなら、TWFEによる推定値が信頼できる推定値であるためには、並行トレンドの仮定のような差分の差分法に関する仮定以外にもいくつかの条件がある。その一つが「処置が段階的に導入されない」ことだ。処置が段階的に導入されていないことは 図 1 のように処置群が同時期に処置を受け、一度処置を受けるとその状態が続くことを意味する。
一方、 図 2 は処置が段階的に導入されている例である。このようなデータを使ってTWFE推定を行うと、バイアスが生じることが知られている。学校内の銃撃事件は全米において同時多発的に起こるものではないので、処置群のカウンティが同時期に処置を受けることは考えられない。したがって、以上のTWFE推定量にはバイアスが含まれている可能性が高い。
しかし、銃撃事件の場合は単なる段階的導入とは考えにくい。 図 2 では処置を受けるタイミングはユニットによって異なるが、一度処置を受けるとその状態は継続する。しかし、銃撃事件の場合、一度銃撃事件が起きたからといって、継続的に銃撃事件が発生するわけではない。つまり、García-Montoya, Arjona, and Lacombe（2022） のデータは @図 3 のような状態であると考えられる。
このように処置がOn-Offでスイッチできるようになっていることを、可逆的（reversible）な処置と呼ぶ。全カウンティを俯瞰することは難しいから、ここではコロラド州のみのデータを眺めてみよう。{panelView}パッケージのpanelview()関数を使えば良い。第一引数として応答変数 ~ 処置変数を指定し、data引数で使用するデータフレームのオブジェクト名、indexではユニットと時間を意味する変数を指定する。
|> filter (state == "06" ) |> panelview (turnout ~ shooting, data = _, index = c ("county" , "year" ),xlab = "Year" , ylab = "County ID in Colorado" )
予想通り、銃撃発生事件の時期も同じでなく、処置状態が可逆的であることが確認できる。したがって、ここでは{fect}のfect()関数を使用し、より厳密に処置効果を推定してみよう。方法の詳細はLiu et al. （2024） を参照されたい。fect()が提供する反事実予測アプローチは処置が段階的に導入されているデータ（ 図 2 ）でも、可逆的な場合（ 図 3 ）でも使える。
ここではdid_fit8と同じモデルをfect()を使って推定し、did_fit8_fectという名のオブジェクトとして保存してみよう。第1引数は回帰式であるが、固定効果で使用する変数（今回はcounty_fとyear_f）は回帰式でなく、index = c("ユニット変数名", "時間変数名")引数で別途指定する必要がある。順番は必ずユニット、時間の順にし、"で囲む必要があることに注意すること。あとはdataでデータフレームのオブジェクト名を指定するだけだが、このままだと推定値の不確実性（標準誤差や信頼区間など）が推定できないため、se = TRUEで標準誤差も計算するように別途指定しておこう（少し時間がかかる）。標準誤差はブートストラップ法で計算されるが、そのためのサンプリング回数が既定値だと200になっている。ブートストラップとしては少し少なめではある。もし、1,000回にしたい場合はnboot = 1000引数を追加してみよう。ただし、その分推定時間が長くなるので注意。
<- fect (demvote ~ shooting + + non_white + change_unem_rate, data = did_df,index = c ("county_f" , "year_f" ),se = TRUE ,nboot = 1000 )
Call:
fect.formula(formula = demvote ~ shooting + population + non_white +
change_unem_rate, data = did_df, index = c("county_f", "year_f"),
se = TRUE, nboots = 1000)
Estimator: fe
Fixed effects: county_f (unit) + year_f (time)
ATT:
ATT S.E. CI.lower CI.upper p.value
Tr obs equally weighted 2.155 0.7429 0.6990 3.611 0.0037217
Tr units equally weighted 2.279 0.6707 0.9642 3.593 0.0006799
Covariates:
Coef S.E. CI.lower CI.upper p.value
population 3.526e-05 5.531e-06 2.442e-05 0.0000461 1.819e-10
non_white 8.493e+01 5.678e+00 7.380e+01 96.0632949 0.000e+00
change_unem_rate -1.375e-01 2.797e-02 -1.923e-01 -0.0826865 8.830e-07
処置効果が2つ出力されるが、基本的には1行目のTr obs equally weighted行を解釈すれば良い。今回の場合、処置効果の点推定値は2.155である。つまり、学校内銃撃事件が発生する場合、民主党候補者の得票率が約2.155pp増加することを意味する。 p 値は0.004であり、 p \(\leq\) 0.05のことから統計的に有意な因果効果であると解釈できる。
TWFEから得られた推定結果（did_fit8）との違いも確認してみよう。
fixest::feols()2.364
1.046
3.682
fect::fect()2.155
0.699
3.611
{fect}パッケージを利用した反事実予測アプローチから得られた点推定量がやや0に近く、標準誤差も大きいことが分かる。分析する側として不確実性の推定値を使いたい気持ちは分かるが、このように処置が段階的に導入されたり、可逆的な場合、通常のTWFEではバイアスが発生することが知られているため、別のアプローチを試すべきだろう。{fect}には様々な機能が提供されており、詳細は公式ホームページ を参照されたい。
また、fect()から得られた処置効果はATTfeols()を使用した場合でも、内部にsunab()で処置の段階的導入を補正した場合や、i()関数でイベントスタディを実行した場合は、ATTが推定対象である。
