UPDATE: 2024-12-07 23:07:21.769122

はじめに

このノートではRを使って因果推論に関する基礎的な理論から分析を実行するための方法をまとめている。ここでは差分の差分法の部分について、下記の矢内先生の資料(資料のライセンスに基づく)をもとに、差分の差分法をまとめている。私の記事は分析の方法メインなので、矢内先生の資料を読む方が数理も学べて良いと思う。

上記の資料のライセンスは下記の通り。

The text of this work is licensed under the Creative Commons Attribution 4.0 International License. The R Code in this work is licensed under the MIT License.

差分の差分法の数理

差分の差分法の数理については下記が詳しいので、ここでは扱わない。

準備

必要なライブラリを読み込んでおく。

library(tidyverse)
library(modelr)
library(modelsummary)
library(fastDummies)

ダミー変数と回帰モデル

01ダミー(2値ダミー)

差分の差分法を発展させると、ダミー変数などを取り込んだ回帰モデルがでてくるので、ダミー変数と回帰モデルのおさらいを行う。サンプルデータを用意する。何らかの時系列データである。ダミー変数の説明は下記を参考にしている。

d1 <- read_csv("~/Desktop/dummy1.csv")
d1
## # A tibble: 31 × 5
##     time     x     d     z     y
##    <dbl> <dbl> <dbl> <dbl> <dbl>
##  1  1981  2251     0  25.3 10204
##  2  1982  2230     0  25.8 10212
##  3  1983  2246     0  26.8 10308
##  4  1984  2290     0  27.3 11832
##  5  1985  2318     0  27.7 11613
##  6  1986  2280     0  26.3 11592
##  7  1987  2123     0  26.1 10571
##  8  1988  2087     0  26.2  9888
##  9  1989  2076     0  26.2 10297
## 10  1990  2055     0  27.5 10463
## # ℹ 21 more rows

変数の関係を可視化する。time=1993に何かがあって、値が少し外れている。

ggplot(d1, aes(x, y, col = as.factor(d))) + 
  geom_point() + 
  ylim(c(5000, 12000)) + 
  scale_color_brewer(palette = "Set1", name = "Dummy") +
  theme_classic()

さまざまな回帰モデルを作成してから、ダミー変数の役割を考える。fit1とfit3では、ダミー変数がモデルに含まれているかどうかが異なり、ダミー変数が1であれば切片である高さの差が−2750下がるという対応関係になっている。また、fit4ではダミー変数に該当するレコードを除外しているが、そのときに求める回帰係数は、ダミー変数を含むモデルfit3と同じ4.087である。

本来であれば、ダミー変数で処理をするのではなく、fit5のようにダミーが表している変数をモデルに含められる方がよい。ただ、このような変数が観察できないので、ダミー変数に影響を偽装させることになる。この例であれば、ダミー変数は1993年とその他の年を区別するものなので、z以外の影響もダミー変数は含んでいると考える。

fit1 <- lm(y ~ x, data = d1)
fit2 <- lm(y ~ d , data = d1)
fit3 <- lm(y ~ x + d, data = d1)
fit4 <- lm(y ~ x, data = d1 %>% filter(d != 1))
fit5 <- lm(y ~ x + z, data = d1)

regs <- list(fit1 = fit1, fit2 = fit2, fit3 = fit3, fit4 = fit4, fit5 = fit5)
msummary(regs, statistic = NULL, gof_omit = "AIC|BIC|Log.Lik.|F|RMSE")
fit1 fit2 fit3 fit4 fit5
(Intercept) 2323.662 9752.533 1868.249 1868.249 -10950.214
x 3.806 4.087 4.087 4.166
d -1941.533 -2750.056
z 471.359
Num.Obs. 31 31 31 30 31
R2 0.674 0.091 0.852 0.838 0.821
R2 Adj. 0.663 0.059 0.842 0.832 0.808

文字ではわかりにくいので可視化しておく。fit1では外れ値の影響を受けて、直線が下に引っ張られているが、fit4では外れ値を削除して計算しているので、ダミー変数を含むfit3と同じ傾きになっている。

sim <- d1 %>% 
  data_grid(x, d) %>%
  gather_predictions(fit1, fit2, fit3, fit4)

ggplot() +
  geom_point(data = d1, aes(x = x, y = y   , col = as.factor(d))) +
  geom_line(data = sim, aes(x = x, y = pred, col = as.factor(d))) +
  facet_grid( ~ model) + 
  ylim(c(5000, 12000)) + 
  scale_color_brewer(palette = "Set1", name = "Dummy") +
  theme_classic()

ここでは、2値ダミーを扱ったが、ダミー変数は多値(3種類以上のダミーはベースが必要)でもよいし、ダミー変数を組み合わせることも可能である。ただ、ダミー変数が行っていることは(係数ダミーを除く)、切片をコントロールして、ダミー変数の組み合わせ分の回帰直線を上下させているといえる。

四半期ダミー

時系列データの分析では四半期ダミーが使われることが多い。ここでは、四半期ダミーを含むデータを用意する。x,yは対数である。

d2 <- read_csv("~/Desktop/dummy2.csv")
d2
## # A tibble: 48 × 8
##     time quater     y     x    q1    q2    q3 trend
##    <dbl> <chr>  <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
##  1  2000 q1      9.81  11.7     1     0     0     1
##  2  2000 q2      9.56  11.7     0     1     0     2
##  3  2000 q3      9.72  11.7     0     0     1     3
##  4  2000 q4      9.66  11.7     0     0     0     4
##  5  2001 q1      9.82  11.7     1     0     0     5
##  6  2001 q2      9.61  11.7     0     1     0     6
##  7  2001 q3      9.72  11.7     0     0     1     7
##  8  2001 q4      9.57  11.7     0     0     0     8
##  9  2002 q1      9.74  11.7     1     0     0     9
## 10  2002 q2      9.52  11.7     0     1     0    10
## # ℹ 38 more rows

この四半期ダミーの役割は、季節による変動を処理するために利用される。四半期は4つのカテゴリーがあるので、ここではQ4を基準にダミー変数を作成する。また、季節ダミーのような短期的な変動ではなく長期的な変動を処理するためにトレンド変数が利用されることも多い。

xが1%増加すると、yが2.6%増加する。四半期ファミーの解釈は、Q4に比べると、Q1は9%増える。ややこしいがダミー変数は対数化してないので、ダミー変数が1単位増えると、yが9%(=100*β=0.09)増える。Q4に比べると、Q2は3.5%減る。

トレンド変数は１期間ごとの変化率に対応するので、毎四半期ごとに0.3%(=100*β=−0.003)減ることになる。

fit6 <- lm(y ~ x + q1 + q2 + q3 + trend, data = d2)
msummary(list(fit6 = fit6), statistic = NULL, gof_omit = "AIC|BIC|Log.Lik.|F|RMSE")
fit6
(Intercept) -20.749
x 2.601
q1 0.090
q2 -0.037
q3 0.126
trend -0.003
Num.Obs. 48
R2 0.900
R2 Adj. 0.888

ダミー変数のベースを変更したいときはファクタのレベルを調整する。

d2 %>%
  mutate(quater = fct_relevel(quater, "q1", "q2", "q3", "q4")) %>%
  lm(y ~ x + quater + trend, data = .) %>% 
  msummary(., statistic = NULL, gof_omit = "AIC|BIC|Log.Lik.|F|RMSE")
(1)
(Intercept) -20.658
x 2.601
quaterq2 -0.128
quaterq3 0.035
quaterq4 -0.090
trend -0.003
Num.Obs. 48
R2 0.900
R2 Adj. 0.888

時間ダミー

次は時間ダミーについてまとめておく。9県×3期間の合計27レコードがある。時系列としては3地点、クロスセクションとしては9県のデータ。

年次を見ると1999, 2004, 2009は3つのカテゴリーもっているので、これは時間ダミーとして利用できる。

d3 <- read_csv("~/Desktop/dummy3.csv")
d3
## # A tibble: 27 × 10
##        i pref      prefno  time     x     y   lny   lnx d2004 d2009
##    <dbl> <chr>      <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
##  1     1 tottori        1  1999   331  2166  7.68  5.80     0     0
##  2     2 shimane        2  1999   401  2575  7.85  5.99     0     0
##  3     3 okayama        3  1999   973  7410  8.91  6.88     0     0
##  4     4 hiroshima      4  1999  1481 11237  9.33  7.3      0     0
##  5     5 yamaguchi      5  1999   764  5662  8.64  6.64     0     0
##  6     6 tokushima      6  1999   405  2751  7.92  6.00     0     0
##  7     7 kagawa         7  1999   533  3714  8.22  6.28     0     0
##  8     8 ehime          8  1999   767  5114  8.54  6.64     0     0
##  9     9 kochi          9  1999   404  2474  7.81  6.00     0     0
## 10    10 tottori        1  2004   311  2319  7.75  5.74     1     0
## # ℹ 17 more rows

時間ダミーを加えると、年によって切片が異なる定数項ダミー、つまり時間ダミーを含む式になる。時間ダミーが有意なことからも、年が変わると、x以外の何らかの要因が変わることで、時間ダミーがyを増加させているとわかる。

時間ダミーは、地域の固有の要因だけではなく、サンプル全体に生じている変化の要因を含んでいる。時間ダミーは年次に共通している複数の説明要因を偽装している。そのため、回帰係数にはさまざまな効果が含まれている。このような理由から、年次によって、切片が異なる定数項ダミーの係数を時間効果(time effect)という。

fit7 <- lm(lny ~ lnx + d2004 + d2009, data = d3)
msummary(list(fit7 = fit7), stars = TRUE, statistic = NULL, gof_omit = "AIC|BIC|Log.Lik.|F|RMSE")
fit7
+ p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001
(Intercept) 1.278***
lnx 1.102***
d2004 0.100**
d2009 0.118***
Num.Obs. 27
R2 0.990
R2 Adj. 0.988

可視化しておく。

sim <- d3 %>% 
  data_grid(lnx, d2004, d2009) %>%
  gather_predictions(fit7)

ggplot() +
  geom_point(data = d3, aes(x = lnx, y = lny   , col = as.factor(time))) +
  geom_line(data = sim %>% filter(d2004 == 0 & d2009 == 0), aes(x = lnx, y = pred), col = "tomato") +
  geom_line(data = sim %>% filter(d2004 == 1 & d2009 == 0), aes(x = lnx, y = pred), col = "limegreen") +
  geom_line(data = sim %>% filter(d2004 == 0 & d2009 == 1), aes(x = lnx, y = pred), col = "royalblue") +
  labs(col = "時間ダミー") + 
  theme_classic()

クロスセクションとダミー変数

さきほどは時間ダミーに着目していたが、都道府県ももちろんダミー変数として扱える。さきほどのモデルに、都道府県ダミーを利用すると、地域固有の個別効果(individuals effects)を追加して推定できる。ここでは都道府県ダミーは\(E\)、時間ダミーを\(T\)とする。

\[ ln y_{i} = \alpha + \beta lnx_{i} + \sum_{r=2}^{9}\delta_{r}Er_{i} + \sum_{t=2}^{3}\gamma_{t}Tt_{i} \]

時間ダミーを含む場合と含まない場合のモデルを作成する。fit8ではlnxや個別効果が有意とされているが、fit9では時間ダミーを追加することで、個別効果が有意ではなくなっている。地域個別の要因はあまり関係なく、どちらかというと、時間ダミーが偽装している何らかの要因の変化がlnyと関係しているようである。

d33 <- d3 %>% 
  select(i, pref, time, lny, lnx) %>% 
  mutate(pref = fct_relevel(pref, 
                            "tottori", "shimane", "okayama", "hiroshima", "yamaguchi", "tokushima", "kagawa", "ehime", "kochi"),
         time = fct_relevel(as.character(time), "1999", "2004", "2009")
         )

fit8 <- lm(lny ~ lnx + pref, data = d33)
fit9 <- lm(lny ~ lnx + pref + time, data = d33)
msummary(list(fit8 = fit8, fit9 = fit9)
         , stars = TRUE, statistic = NULL, gof_omit = "AIC|BIC|Log.Lik.|F|RMSE")
fit8 fit9
+ p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001
(Intercept) 9.650*** 4.259*
lnx -0.337+ 0.593+
prefshimane 0.237*** 0.046
prefokayama 1.589*** 0.553
prefhiroshima 2.174*** 0.741
prefyamaguchi 1.258*** 0.453
preftokushima 0.349*** 0.144+
prefkagawa 0.702*** 0.243
prefehime 1.127*** 0.337
prefkochi 0.150** -0.036
time2004 0.077***
time2009 0.079*
Num.Obs. 27 27
R2 0.998 0.999
R2 Adj. 0.997 0.998

ここで、地域\(r\)を1から9として、時間\(t\)を1から3とすると、地域\(r\)の時間\(t\)での\(y\)は、\(y_{rt}\)と表現できる。\(x\)についも\(x_{rt}\)と表現できる。観測できない要因をダミー変数をお用いて、小屋いごとに特有で時間を通じて変化しない効果、例えば、地域ごとの地理的条件や歴史的背景などを偽装する。この場合、ダミー変数の係数は固定効果(fixed effects)とよばれる。さきほどと同じような発想で、モデル全体の切片をなくして、傾きは同じで、地域ごとに切片が異なる回帰式を推定する。

\[ ln y_{rt} = \beta lnx_{rt} + \alpha_{r} + u_{rt} \quad (r = 1...9, t = 1...3) \]

ようするに、展開すると下記のようになる。

\[ ln y_{rt} = \beta lnx_{rt} + \alpha_{1}E1_{r} + \alpha_{2}E2_{r} + ... + \alpha_{9}E9_{r} + u_{rt} \]

d333 <- d3 %>% 
  dummy_cols(select_columns = c("pref"), remove_first_dummy = FALSE) %>% 
  select(time, lny, lnx, starts_with("pref_")) %>% 
  mutate(time = fct_relevel(as.character(time), "2009", "2004", "1999"))

fit10 <- lm(lny ~ 
              0 + 
              lnx + 
              pref_tottori + 
              pref_shimane + 
              pref_okayama + 
              pref_hiroshima + 
              pref_yamaguchi + 
              pref_tokushima + 
              pref_kagawa + 
              pref_ehime + 
              pref_kochi + 
              time, 
            data = d333)
msummary(list(fit10 = fit10), stars = TRUE, statistic = NULL, gof_omit = "AIC|BIC|Log.Lik.|F|RMSE")
fit10
+ p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001
lnx 0.593+
pref_tottori 4.259*
pref_shimane 4.305*
pref_okayama 4.812*
pref_hiroshima 4.999+
pref_yamaguchi 4.712*
pref_tokushima 4.403*
pref_kagawa 4.502*
pref_ehime 4.595+
pref_kochi 4.222*
time2009 0.079*
time2004 0.077***
Num.Obs. 27
R2 1.000
R2 Adj. 1.000

差分の差分法

差分の差分法の基礎

Card and Krueger (1994)の論文で示されている最低賃金の上昇が失業に与える影響の研究データを利用する。この研究は、最低賃金の引き上げと雇用の関係を調べたもの。1992年にアメリカのニュージャージー州(NJ)で最低時給が4.25ドルから5.05ドルに上昇したのに対し、隣のペンシルバニア州(PA)では最低時給の上昇がなかったという事実を対照群として利用することで分析している。

# fulltime_before:最低時給上昇前のフルタイム労働者の数
# parttime_before:最低時給上昇前のパートタイム労働者の数
# wage_before    :最低時給上昇前の賃金
# fulltime_after :最低時給上昇後のフルタイム労働者の数
# parttime_after :最低時給上昇後のパートタイム労働者の数
# wage_after     :最低時給上昇後の賃金
# full_prop_before = fulltime_before / (fulltime_before + parttime_before),
# full_prop_after  = fulltime_after  / (fulltime_after + parttime_after)
df <- read_csv("~/Desktop/public_mini.csv")
df
## # A tibble: 358 × 9
##    state fulltime_before parttime_before wage_before fulltime_after
##    <chr>           <dbl>           <dbl>       <dbl>          <dbl>
##  1 PA               20                20        5                 0
##  2 PA                6                26        5.5              28
##  3 PA               50                35        5                15
##  4 PA               10                17        5                26
##  5 PA                2                 8        5.25              3
##  6 PA                2                10        5                 2
##  7 PA                2.5              20        5                 1
##  8 PA               40                30        5                 9
##  9 PA                8                27        5                 7
## 10 PA               10.5              30        5.5              18
## # ℹ 348 more rows
## # ℹ 4 more variables: parttime_after <dbl>, wage_after <dbl>,
## #   full_prop_before <dbl>, full_prop_after <dbl>

ニュージャージー州(NJ)とペンシルバニア州(PA)の地理的関係は下記の通り。

引用: 在ニューヨーク日本国領事館
引用: 在ニューヨーク日本国領事館

基礎集計として、最低時給上昇前の賃金wage_beforeと最低時給上昇後の賃金wage_afterが5.05ドル未満のお店の割合を計算する。ペンシルバニア州は最低賃金の引き上げがなされてないので、前後であまり差がないが、ニュージャージー州のほとんどの店は最低賃金を5.05ドル以上にあげていることがわかる。

df %>% 
  group_by(state) %>% 
  summarize(before = mean(wage_before < 5.05),
            after = mean(wage_after < 5.05),
            .groups = "drop")
## # A tibble: 2 × 3
##   state before   after
##   <chr>  <dbl>   <dbl>
## 1 NJ     0.911 0.00344
## 2 PA     0.940 0.955

単純に法律施行後のニュージャージー州とペンシルバニア州のフルタイム労働者の割合の差を計算すると、ニュージャージー州の方がフルタイム労働者の割合が4ポイント程度高いが、これはニュージャージー州はそもそも前後で上昇トレンドがあったからかもしれない。つまり、法律は関係なく、勝手にフルタイム労働者の割合が上がっていた可能性もある。

df %>% 
  group_by(state) %>% 
  summarize(fulltime = mean(full_prop_after),
            .groups = "drop") %>% 
  pivot_wider(names_from = state, values_from = fulltime) %>% 
  set_names(c("NJ_after","PA_after")) %>% 
  mutate(diff_after = NJ_after - PA_after)
## # A tibble: 1 × 3
##   NJ_after PA_after diff_after
##      <dbl>    <dbl>      <dbl>
## 1    0.320    0.272     0.0481

次に法律施行前後のニュージャージー州のフルタイム労働者の割合の差を計算する。ニュージャージー州のフルタイム労働者の前後比較では、2.4ポイント増加している。これも、法律は関係なく、勝手にフルタイム労働者の割合が上がっていた可能性もある。

df %>% 
  filter(state == "NJ") %>% 
  # summarize(across(.cols = starts_with("full_"), mean))でもOK
  summarise(full_prop_before = mean(full_prop_before),
            full_prop_after = mean(full_prop_after),
            .groups = "drop") %>% 
  mutate(diff_before_after = full_prop_after - full_prop_before)
## # A tibble: 1 × 3
##   full_prop_before full_prop_after diff_before_after
##              <dbl>           <dbl>             <dbl>
## 1            0.297           0.320            0.0239

つまり、法律は関係なく、勝手にフルタイム労働者が増加している問題を対処しないといけない。このためにペンシルバニア州のデータを利用して差分の差分法を行う。法律は関係なく、勝手に上がっているならば、ペンシルバニア州も勝手に上がってると考えられるので、ペンシルバニア州の上がった分を考慮した上で、ニュージャージー州のフルタイム労働者の増加を検討する。

df2 <- df %>% 
  group_by(state) %>%
  summarize(across(.cols = starts_with("full_"), mean),
            .groups = "drop") %>% 
  mutate(diff_before_after = full_prop_after - full_prop_before)
df2
## # A tibble: 2 × 4
##   state full_prop_before full_prop_after diff_before_after
##   <chr>            <dbl>           <dbl>             <dbl>
## 1 NJ               0.297           0.320            0.0239
## 2 PA               0.310           0.272           -0.0377

2点間の可視化を行うとわかるが、個体間比較や前後比較の推定値より大きくなったは、ペンシルバニア州の前後で、後ろの方が割合が下がっているため、差分の差分法による推定値が、個体間比較や前後比較の推定値より大きくなる。

df2 %>% 
  select(- diff_before_after) %>% 
  pivot_longer(cols = starts_with("full"),
               names_to = "before_after",
               names_prefix = "full_prop_",
               values_to = "prop") %>% 
  ggplot(aes(x = before_after, y = prop, group = state)) +
  geom_line(aes(color = state)) +
  geom_point(aes(shape = state)) +
  geom_text(aes(x = before_after, y = prop+0.005, label = round(prop,3))) + 
  ylim(0.25, 0.35) +
  theme_bw()

ニュージャージー州とペンシルバニア州の前後差diff_before_afterを使って更に、差を計算することで、差分の差分法による推定値が計算できる。最低賃金の引き上げは、フルタイム労働者の割合を6.2ポイント増加させる結果となった。

df2 %>% 
  select(state, diff_before_after) %>% 
  pivot_wider(names_from = state, values_from = diff_before_after) %>% 
  set_names(c("NJ_after","PA_after")) %>% 
  mutate(did = NJ_after - PA_after)
## # A tibble: 1 × 3
##   NJ_after PA_after    did
##      <dbl>    <dbl>  <dbl>
## 1   0.0239  -0.0377 0.0616

これは回帰係数の交互作用の回帰係数と一致する。

df %>% 
  select(state, starts_with("full_")) %>% 
  pivot_longer(cols = starts_with("full"),
               names_to = "time",
               names_prefix = "full_prop_",
               values_to = "prop") %>% 
  mutate(D = if_else(state == "NJ", 1, 0),
         P = if_else(time == "after", 1, 0)) %>% 
  lm(prop ~ D * P, data = .) %>% 
  tidy() %>% 
  select(term, estimate)
## # A tibble: 4 × 2
##   term        estimate
##   <chr>          <dbl>
## 1 (Intercept)   0.310 
## 2 D            -0.0134
## 3 P            -0.0377
## 4 D:P           0.0616

簡単に図解すると下記のようなイメージ。黄色が介入効果、赤色が時間効果、緑色が介入効果と時間効果の効果となるので、交互作用項の回帰係数がDID推定値と解釈できる。

図: DIDと回帰係数
図: DIDと回帰係数

以上より、1992年4月にアメリカのニュージャージー州で最低賃金が引き上げられたため、フルタイム労働者の割合は6.1%ポイント増えたと考えられる。

差分の差分法の発展

Mastering ’MetricsのAngrist and Pischke(2009)の5.2節(pp.191-203)にある法定飲酒年齢 (minimum legal drinking age; MLDA)の変更が18歳から20歳までの若者の死亡率との関連を調べた研究のデータを利用する。このデータは州と年のパネルデータ。

# すべての死 All deaths  dtype = 1
# 自動車事故による死 Motor vehicle accidents dtype = 2
# 自殺    Suicide dtype = 3
# 内蔵疾患による死  All internal causes dtype = 6

df <- read_csv("~/Desktop/deaths.csv")
df2 <- df %>% 
  filter(year <= 1983,
         agegr == 2,      # 18-21 years old
         dtype == 1) %>%  # all deaths      
  mutate(state = factor(state),
         year_fct = factor(year))
df2 %>% 
  select(year, state, age, mrate, legal) %>% 
  filter(state == 1)
## # A tibble: 14 × 5
##     year state   age mrate legal
##    <dbl> <fct> <dbl> <dbl> <dbl>
##  1  1970 1      19.0  154. 0    
##  2  1971 1      19.0  162. 0    
##  3  1972 1      19.0  159. 0    
##  4  1973 1      19.0  141. 0    
##  5  1974 1      19.0  143. 0    
##  6  1975 1      19.0  148. 0.294
##  7  1976 1      19.0  133. 0.665
##  8  1977 1      19.0  146. 0.668
##  9  1978 1      19.0  139. 0.667
## 10  1979 1      19.0  135. 0.668
## 11  1980 1      19.0  120. 0.672
## 12  1981 1      19.0  106. 0.679
## 13  1982 1      19.0  118. 0.676
## 14  1983 1      19.0  128. 0.680

死亡率の変化を可視化しておく。可視化した結果を見る限り死亡率は下がっているようにも見える。

ggplot() +
  geom_line(data = df2, aes(x = year, y = mrate, group = state)) +
  geom_line(data = df2 %>% group_by(year) %>% summarise(mean_mrate = mean(mrate)), 
            aes(x = year, y = mean_mrate), col = "tomato", size = 2) +
  theme_classic()

このデータは単純に差分の差分法を利用できない。なぜならば、処置のパターンが1つではなく、処置のタイミングも州によって異なるためである。処置のパタンについては、下記の通り。詳細は参照先のこちらを確認のこと。

  • 21歳から18歳(3歳引き下げ)
  • 21歳から20歳(1歳引き下げ)
  • 18歳から21歳(3歳引き上げ)

回帰分析を利用してDID回帰を行うが、さきほどの回帰モデルには下記が必要だった。

  • 処置群と統制群の個体差の区別をD: state
  • 時間前後の区別を表すP: year_fct
  • 処置後の観測値を表すD×P: legal

この分析では、処置後の観測値であるかどうかは、これら2つの交差項では表現できない。そのため、研究データでは、この問題を解決するために、18歳から20歳の人が飲酒できる割合を表す変数legalを作ることで対処している。

この変数は、t年のs州で、18歳から20歳までの人口のうち何割が合法的に飲酒できるかを表している。MLDAが21歳だとすると、20歳以下で合法的に飲酒できる人はいないので、この変数の値は0になる。MLDAが18歳なら、18歳以上の全員が合法的に飲酒できるので、この変数の値は1になる。この値が大きいほど、18歳から20歳の人が飲酒しやすいことを意味する。MLDAを引き下げるという処置を実行すると、legalの値は大きくなる。

なんとなくわかるようでわからない感じのわたしの理解力助けるためにメモを残しておく。正確かどうかはわからない。下記の州(state=1)では、1975年から飲酒年齢を引き下げている。おそらく19歳から飲酒可能にしていると思われる。1970年から1974年までは、「18歳から20歳までの人口」のうち何割が合法的に飲酒できるかを表すlegalは0である、つまり、この期間に「19歳から20歳」も人はそもそも飲酒が法律的に認められない。

1975年は何かあったのか、年の途中から飲酒年齢が引き下げられたのかわからないが、数値が他の年と比べて小さいため、年の途中でMLDAが変更されたために調整されていると思われる。その他の年では、だいたい67%なので、「19歳から20歳」の人口を「18歳から20歳」の人口で割っていると思われる。超絶簡単なイメージで説明すると、各年代100人いて、毎年100人が18歳(in)になり、100人が21歳(out)になるのであれば、18歳から20歳という範囲termには、常に一定inoutが発生することで、legalは200(19-20)/300(18-20)=0.66となる。浅野先生の資料が詳しいので、こちらから画像をお借りする。

#     year state mrate legal
#    <dbl> <fct> <dbl> <dbl>
#  1  1970 1      154. 0    
#  2  1971 1      162. 0    
#  3  1972 1      159. 0    
#  4  1973 1      141. 0    
#  5  1974 1      143. 0    
#  6  1975 1      148. 0.29
#  7  1976 1      133. 0.66
#  8  1977 1      146. 0.66
#  9  1978 1      139. 0.66
# 10  1979 1      135. 0.66
# 11  1980 1      120. 0.67
# 12  1981 1      106. 0.67
# 13  1982 1      118. 0.67
# 14  1983 1      128. 0.68

飲酒引き上げをいつから行ったかを調べると、legalが14年に渡りすべて0の州が12州あるので、おそらくこれらの州は飲酒年齢を引き下げていない。

df2 %>% 
  group_by(state) %>% 
  select(year, state, mrate, legal) %>% 
  mutate(tmp = ifelse(legal > 0, 1, 0),
         cond_start = if_else(tmp == 1, year, NA_real_),
         is_raw = if_else(all(is.na(cond_start)), FALSE, TRUE),
         start_year = if_else(is_raw == FALSE, min(year), min(cond_start, na.rm = TRUE))
  ) %>%
  ungroup() %>% 
  distinct(start_year, state, is_raw) %>% 
  group_by(start_year, is_raw) %>% 
  summarise(cnt = n(),
            state_group = paste0(state, collapse = ",")) %>% 
  arrange(start_year)
## # A tibble: 8 × 4
## # Groups:   start_year [7]
##   start_year is_raw   cnt state_group                                  
##        <dbl> <lgl>  <int> <chr>                                        
## 1       1970 FALSE     12 5,6,18,21,29,32,35,38,41,42,49,53            
## 2       1970 TRUE      16 2,8,11,16,20,22,23,28,31,36,37,39,45,46,54,55
## 3       1971 TRUE       2 47,50                                        
## 4       1972 TRUE       8 4,9,10,13,15,19,26,44                        
## 5       1973 TRUE       9 12,17,25,27,30,33,34,48,56                   
## 6       1974 TRUE       2 24,51                                        
## 7       1975 TRUE       1 1                                            
## 8       1976 TRUE       1 40

データへの理解が深まったところで、DID回帰を行う。モデルは下記の通り。

\[ y_{st} = \alpha + \delta LEGAL_{st} + \sum_{k=1}^{50}\beta_{k}STATE_{ks} + \sum_{j=1971}^{1983}\gamma_{j}YEAR_{jt} + e_{st} \]

記号 内容
\(s\)=[1,2,3,…,51] 州(\(s\)=51は D.C.)
\(t\)=[1970,1972,…,1983] 年(1970年は参照カテゴリ)
\(Y_{st}\) \(t\)年の\(s\)州での18から20歳の死者数(10万人あたり)
\(LEGAL_{st}\) 処置効果(18歳から20歳の飲酒可能な人の割合)
\(STATE_{ks}\) \(k=s\)の時に1(D.C.を除く50州のダミー)
\(YEAR_{jt}\) \(j=t\)の時に1(1970年を除く13個の年ダミー)

ダミー変数も明示して、どのようなモデルなのかわかりやすくしておく。

lm(mrate ~ 0 + legal + state + year_fct, data = df2) %>% 
  tidy() %>% 
  select(term, estimate, std.error) %>% 
  as.data.frame()
##            term   estimate std.error
## 1         legal  10.804141  3.137576
## 2        state1 153.669154  5.200633
## 3        state2 261.555750  5.310643
## 4        state4 191.472585  5.248148
## 5        state5 159.695144  5.279205
## 6        state6 152.931870  5.279205
## 7        state8 141.285486  5.629773
## 8        state9 108.455150  5.390877
## 9       state10 121.624050  5.192291
## 10      state11 157.957094  5.629773
## 11      state12 158.734019  5.319487
## 12      state13 153.589183  5.371284
## 13      state15 118.168566  5.460636
## 14      state16 168.354686  5.271968
## 15      state17 138.629450  5.191376
## 16      state18 136.155379  5.279205
## 17      state19 133.625685  5.332511
## 18      state20 135.131848  5.629773
## 19      state21 149.327154  5.279205
## 20      state22 163.227278  5.629773
## 21      state23 139.547863  5.265282
## 22      state24 130.551487  5.304303
## 23      state25 111.193610  5.248715
## 24      state26 137.091383  5.226174
## 25      state27 131.932044  5.271285
## 26      state28 153.314569  5.629773
## 27      state29 152.345279  5.279205
## 28      state30 215.685715  5.292326
## 29      state31 138.077315  5.238494
## 30      state32 244.794469  5.279205
## 31      state33 130.600105  5.242543
## 32      state34 113.800240  5.329778
## 33      state35 220.562702  5.279205
## 34      state36 125.282941  5.600285
## 35      state37 141.257107  5.623048
## 36      state38 156.074033  5.279205
## 37      state39 119.946170  5.592418
## 38      state40 154.479735  5.222120
## 39      state41 155.908843  5.279205
## 40      state42 123.763698  5.279205
## 41      state44  97.864771  5.342273
## 42      state45 136.233762  5.629773
## 43      state46 146.162173  5.555958
## 44      state47 161.672377  5.417643
## 45      state48 161.490118  5.330628
## 46      state49 134.920186  5.279205
## 47      state50 139.243233  5.483938
## 48      state51 124.343565  5.314165
## 49      state53 145.959715  5.279205
## 50      state54 149.895274  5.626146
## 51      state55 132.134166  5.629773
## 52      state56 245.722502  5.234076
## 53 year_fct1971  -5.501817  3.434680
## 54 year_fct1972  -6.561624  3.467123
## 55 year_fct1973  -5.057369  3.573380
## 56 year_fct1974 -14.686167  3.659049
## 57 year_fct1975 -19.981426  3.686623
## 58 year_fct1976 -23.397556  3.697394
## 59 year_fct1977 -21.096330  3.710735
## 60 year_fct1978 -23.066713  3.695079
## 61 year_fct1979 -21.617661  3.636559
## 62 year_fct1980 -20.248734  3.600693
## 63 year_fct1981 -33.263357  3.580025
## 64 year_fct1982 -41.553328  3.567483
## 65 year_fct1983 -45.018810  3.542118

平均処置効果legalの効果は10.8と推定された。18歳から20歳の人が飲酒できるようになると、10万人あたりの死者数が約11人増えたことを意味する。つまり、MLDAを下げると、飲酒可能な人の中で死亡率があがることにあんる。

このような時系列データで回帰モデルを計算する際は、同一州内の観測値が類似する(系列相関)という問題が発生し、標準誤差を過小に推定してしまう。このような州の類似性をクラスタとして調整して計算するクラスタ標準誤差(cluster-robust standard error)を利用するのが望ましい。estimatrパッケージのlm_robust関数で計算できる。state56などが出てきて変な感じもするが、これはstateが連番ではなく欠番があるため。

library(estimatr)
# mrate ~ 0 + legal + state + year_fct のように書いてもよい
lm_robust(mrate ~ legal,
          fixed_effects = state + year_fct,
          data = df2, 
          #clusters = state, 
          # se_type = "stata"
          ) 
##       Estimate Std. Error  t value     Pr(>|t|) CI Lower CI Upper  DF
## legal 10.80414   2.724831 3.965068 8.153258e-05 5.453592 16.15469 649

Angrist and Pischke (2009)の論文の数値とestimateが一致している。std.errorは一致していないので、修正する。

同一州内の観測値が類似するという系列相関が発生し、クラスター内で相関してしまうので、これを修正するクラスター標準誤差を計算することで、標準誤差も修正できる。クラスタ化を考慮しない分析では、標準誤差が過小評価されてしまう。

lm_robust(mrate ~ legal,
          fixed_effects = state + year_fct,
          data = df2, 
          clusters = state, 
          se_type = "stata"
          ) 
##       Estimate Std. Error  t value   Pr(>|t|) CI Lower CI Upper DF
## legal 10.80414   4.592205 2.352713 0.02261311 1.580427 20.02786 50

次は重み付きDID回帰の結果を確認する。州の人口を考慮した重み付き回帰を実行することで、人口の多い州、人口の少ない州の重みを考慮できる。人口の多い州の影響を反映し、人口が小さい州の過度な影響に左右しないようにできる。

lm_robust(mrate ~ legal,
          fixed_effects = state + year_fct,
          data = df2, 
          clusters = state, 
          se_type = "stata",
          weights = pop
          ) 
##       Estimate Std. Error  t value    Pr(>|t|) CI Lower CI Upper DF
## legal 12.41251   4.599346 2.698755 0.009467674 3.174449 21.65056 50

こちらも数字が一致していることが確認できる。

差分の差分法では、平行トレンドが仮定されるが、今回のように個体と期間の数が多い場合、個体ごとの時間トレンドの違いを説明する変数をモデルに加えることで、並行トレンド仮定を緩めて平均置効果を推定する。トレンドが完全に平行でなくても、平均置効果を推定することができるようになるらしい。

\[ y_{st} = \alpha + \delta LEGAL_{st} + \sum_{k=1}^{50}\beta_{k}STATE_{ks} + \sum_{j=1971}^{1983}\gamma_{j}YEAR_{jt} + \sum_{k=1}^{50}\theta_{k}(STATE_{ks} * time ) + e_{st} \]

下記の分析では、各州の時間トレンドは線形、つまり、処置がないときの死亡率は、変化しない、単調増加、単調減少のどれかであることを仮定する。州ごとの時間トレンドはstate × yearで表している。year_fctではないので注意。stateは個体効果、year_fctは時間効果、state:yearは州ごとのトレンド効果という感じ。

# 時間トレンドを含める場合、fixed_effects = state + year_fct + state:yearのようには書けない。実行はできるが結果が意図通りではない。
lm_robust(mrate ~ 0 + legal + state + year_fct + state:year,
          data = df2, 
          clusters = state, 
          se_type = "stata"
          ) 
##                   Estimate   Std. Error      t value     Pr(>|t|)      CI Lower
## legal         8.466624e+00 5.097812e+00   1.66083508 1.030056e-01 -1.772632e+00
## state1        8.301347e+03 3.044177e+02  27.26959027 1.115951e-31  7.689906e+03
## state2        4.827761e+03 8.955706e+02   5.39070993 1.903719e-06  3.028955e+03
## state4        1.302659e+04 5.537238e+02  23.52543262 1.063173e-28  1.191441e+04
## state5        1.048942e+04 1.001069e+03  10.47822266 3.260953e-14  8.478716e+03
## state6        6.196865e+03 1.001069e+03   6.19024870 1.107054e-07  4.186159e+03
## state8        5.716430e+03 1.006167e+03   5.68139477 6.815124e-07  3.695485e+03
## state9        6.674982e+02 4.572720e+02   1.45974008 1.506164e-01 -2.509596e+02
## state10       3.631577e+03 7.786071e+02   4.66419638 2.343847e-05  2.067698e+03
## state11       1.227551e+04 1.006167e+03  12.20027891 1.323410e-16  1.025457e+04
## state12       5.720604e+03 4.436259e+02  12.89510890 1.597872e-17  4.829555e+03
## state13       1.256791e+04 6.256728e+02  20.08703292 1.363659e-25  1.131121e+04
## state15       4.658527e+02 4.416915e+02   1.05470160 2.966318e-01 -4.213107e+02
## state16       9.041919e+03 7.574949e+02  11.93660615 3.000397e-16  7.520445e+03
## state17       6.054395e+03 1.027004e+03   5.89520363 3.183282e-07  3.991598e+03
## state18       5.723788e+03 1.001069e+03   5.71767694 5.991048e-07  3.713082e+03
## state19       7.885272e+03 6.622383e+02  11.90700073 3.291062e-16  6.555127e+03
## state20       4.221631e+03 1.006167e+03   4.19575781 1.112077e-04  2.200686e+03
## state21       4.612981e+03 1.001069e+03   4.60805593 2.833135e-05  2.602275e+03
## state22       2.587834e+03 1.006167e+03   2.57197409 1.313064e-02  5.668893e+02
## state23       4.365182e+03 1.298806e+03   3.36091918 1.494978e-03  1.756454e+03
## state24       7.441929e+03 1.666805e+02  44.64787424 5.897876e-42  7.107142e+03
## state25       5.836604e+03 8.822989e+02   6.61522257 2.400711e-08  4.064454e+03
## state26       8.909473e+03 1.470647e+03   6.05819991 1.777314e-07  5.955592e+03
## state27       7.495299e+03 5.554100e+02  13.49507448 2.705652e-18  6.379725e+03
## state28       6.802820e+03 1.006167e+03   6.76112719 1.419155e-08  4.781875e+03
## state29       6.205279e+03 1.001069e+03   6.19865439 1.074156e-07  4.194574e+03
## state30       4.353277e+03 5.200036e+02   8.37162864 4.451934e-11  3.308819e+03
## state31       7.278361e+03 1.046177e+03   6.95710270 7.003839e-09  5.177053e+03
## state32       1.542771e+04 1.001069e+03  15.41123735 1.253420e-20  1.341700e+04
## state33       4.332938e+03 8.371362e+02   5.17590600 4.037874e-06  2.651501e+03
## state34       6.315848e+03 5.903023e+02  10.69934418 1.574307e-14  5.130191e+03
## state35       8.165233e+03 1.001069e+03   8.15651531 9.535557e-11  6.154527e+03
## state36       7.600656e+03 1.104340e+03   6.88253025 9.162855e-09  5.382523e+03
## state37       8.278937e+03 1.028541e+03   8.04920701 1.396128e-10  6.213052e+03
## state38       7.808874e+03 1.001069e+03   7.80053666 3.387820e-10  5.798168e+03
## state39       4.541523e+03 1.127148e+03   4.02921533 1.905153e-04  2.277579e+03
## state40                 NA           NA           NA           NA            NA
## state41       6.217085e+03 1.001069e+03   6.21044760 1.029636e-07  4.206379e+03
## state42       4.168779e+03 1.001069e+03   4.16432834 1.231816e-04  2.158073e+03
## state44       5.072893e+03 9.060537e+02   5.59888829 9.131299e-07  3.253031e+03
## state45       7.886137e+03 1.006167e+03   7.83780465 2.965612e-10  5.865192e+03
## state46       1.110184e+04 7.813317e+02  14.20887059 3.468779e-19  9.532490e+03
## state47       6.600759e+03 9.424862e+02   7.00355980 5.924474e-09  4.707719e+03
## state48       1.663799e+03 3.722708e+02   4.46932571 4.511288e-05  9.160715e+02
## state49       5.634743e+03 1.001069e+03   5.62872671 8.215202e-07  3.624037e+03
## state50       5.042159e+03 5.147369e+02   9.79560352 3.206043e-13  4.008279e+03
## state51       7.508273e+03 1.081705e+02  69.41149735 2.203044e-51  7.291007e+03
## state53       4.129961e+03 1.001069e+03   4.12555127 1.396872e-04  2.119255e+03
## state54       5.137581e+03 1.018217e+03   5.04566542 6.347492e-06  3.092432e+03
## state55       4.983706e+03 1.006167e+03   4.95316142 8.736521e-06  2.962761e+03
## state56       1.322648e+04 4.629416e+02  28.57052165 1.246070e-32  1.229664e+04
## year_fct1971 -2.198310e+00 3.289304e+00  -0.66832057 5.070027e-01 -8.805072e+00
## year_fct1972  2.824255e-01 3.287292e+00   0.08591432 9.318776e-01 -6.320295e+00
## year_fct1973  5.409201e+00 4.081651e+00   1.32524841 1.911129e-01 -2.789036e+00
## year_fct1974 -7.715550e-01 5.133235e+00  -0.15030581 8.811281e-01 -1.108196e+01
## year_fct1975 -2.765714e+00 4.883387e+00  -0.56635153 5.736883e-01 -1.257429e+01
## year_fct1976 -2.916808e+00 5.815619e+00  -0.50154734 6.181873e-01 -1.459782e+01
## year_fct1977  2.654118e+00 6.043832e+00   0.43914484 6.624486e-01 -9.485277e+00
## year_fct1978  3.895767e+00 6.968034e+00   0.55909132 5.785951e-01 -1.009994e+01
## year_fct1979  8.463277e+00 7.433437e+00   1.13854162 2.603203e-01 -6.467220e+00
## year_fct1980  1.299063e+01 6.703434e+00   1.93790698 5.829206e-02 -4.736118e-01
## year_fct1981  3.166470e+00 6.011374e+00   0.52674642 6.006989e-01 -8.907730e+00
## year_fct1982 -1.914022e+00 6.934406e+00  -0.27601810 7.836718e-01 -1.584219e+01
## year_fct1983 -2.209433e+00 6.524369e+00  -0.33864309 7.362966e-01 -1.531401e+01
## state1:year  -4.132824e+00 1.550839e-01 -26.64895054 3.280975e-31 -4.444319e+00
## state2:year  -2.320492e+00 4.535085e-01  -5.11675532 4.960534e-06 -3.231390e+00
## state4:year  -6.504238e+00 2.808405e-01 -23.15990174 2.180789e-28 -7.068323e+00
## state5:year  -5.237299e+00 5.085935e-01 -10.29761283 5.937506e-14 -6.258839e+00
## state6:year  -3.068924e+00 5.085935e-01  -6.03413898 1.937215e-07 -4.090464e+00
## state8:year  -2.830560e+00 5.085935e-01  -5.56546630 1.027804e-06 -3.852100e+00
## state9:year  -2.929629e-01 2.314801e-01  -1.26560717 2.115196e-01 -7.579044e-01
## state10:year -1.786546e+00 3.953350e-01  -4.51906973 3.820321e-05 -2.580600e+00
## state11:year -6.140660e+00 5.085935e-01 -12.07380628 1.957577e-16 -7.162200e+00
## state12:year -2.824228e+00 2.248165e-01 -12.56237037 4.363797e-17 -3.275785e+00
## state13:year -6.291108e+00 3.167424e-01 -19.86190900 2.252743e-25 -6.927304e+00
## state15:year -1.859363e-01 2.233996e-01  -0.83230357 4.091954e-01 -6.346475e-01
## state16:year -4.499853e+00 3.838151e-01 -11.72401503 5.842944e-16 -5.270769e+00
## state17:year -3.003727e+00 5.209485e-01  -5.76588061 5.047252e-07 -4.050083e+00
## state18:year -2.838061e+00 5.085935e-01  -5.58021498 9.755391e-07 -3.859601e+00
## state19:year -3.932112e+00 3.353737e-01 -11.72456686 5.832800e-16 -4.605730e+00
## state20:year -2.077388e+00 5.085935e-01  -4.08457429 1.594564e-04 -3.098928e+00
## state21:year -2.269390e+00 5.085935e-01  -4.46208976 4.621468e-05 -3.290930e+00
## state22:year -1.236562e+00 5.085935e-01  -2.43133680 1.866969e-02 -2.258102e+00
## state23:year -2.148272e+00 6.577204e-01  -3.26623935 1.972150e-03 -3.469342e+00
## state24:year -3.709409e+00 8.475952e-02 -43.76392330 1.564213e-41 -3.879654e+00
## state25:year -2.907116e+00 4.470776e-01  -6.50248531 3.602992e-08 -3.805098e+00
## state26:year -4.448780e+00 7.448898e-01  -5.97239970 2.416031e-07 -5.944935e+00
## state27:year -3.735778e+00 2.815748e-01 -13.26744370 5.279196e-18 -4.301338e+00
## state28:year -3.374128e+00 5.085935e-01  -6.63423280 2.241807e-08 -4.395668e+00
## state29:year -3.073478e+00 5.085935e-01  -6.04309337 1.876100e-07 -4.095018e+00
## state30:year -2.103671e+00 2.635673e-01  -7.98153227 1.776348e-10 -2.633061e+00
## state31:year -3.622991e+00 5.300505e-01  -6.83518137 1.086756e-08 -4.687629e+00
## state32:year -7.692745e+00 5.085935e-01 -15.12552587 2.716740e-20 -8.714285e+00
## state33:year -2.136541e+00 4.242632e-01  -5.03588699 6.565986e-06 -2.988699e+00
## state34:year -3.148105e+00 2.989881e-01 -10.52919656 2.755537e-14 -3.748640e+00
## state35:year -4.030592e+00 5.085935e-01  -7.92497753 2.172924e-10 -5.052133e+00
## state36:year -3.792000e+00 5.583276e-01  -6.79171234 1.271051e-08 -4.913434e+00
## state37:year -4.127069e+00 5.199279e-01  -7.93777064 2.076060e-10 -5.171374e+00
## state38:year -3.882922e+00 5.085935e-01  -7.63462710 6.132997e-10 -4.904462e+00
## state39:year -2.246956e+00 5.698843e-01  -3.94282728 2.509888e-04 -3.391602e+00
## state40:year  6.770603e-02 1.265590e-03  53.49759282 8.382872e-46  6.516402e-02
## state41:year -3.077648e+00 5.085935e-01  -6.05129270 1.821823e-07 -4.099189e+00
## state42:year -2.057582e+00 5.085935e-01  -4.04563158 1.807408e-04 -3.079122e+00
## state44:year -2.527280e+00 4.586960e-01  -5.50970713 1.251701e-06 -3.448598e+00
## state45:year -3.930868e+00 5.085935e-01  -7.72889980 4.376565e-10 -4.952409e+00
## state46:year -5.552889e+00 3.950030e-01 -14.05784107 5.329073e-19 -6.346276e+00
## state47:year -3.267904e+00 4.768916e-01  -6.85250930 1.020971e-08 -4.225769e+00
## state48:year -7.702950e-01 1.886740e-01  -4.08267766 1.604342e-04 -1.149258e+00
## state49:year -2.793634e+00 5.085935e-01  -5.49286218 1.328390e-06 -3.815174e+00
## state50:year -2.490604e+00 2.602957e-01  -9.56836529 6.945134e-13 -3.013424e+00
## state51:year -3.746099e+00 5.512076e-02 -67.96166578 6.261489e-51 -3.856812e+00
## state53:year -2.026712e+00 5.085935e-01  -3.98493480 2.195001e-04 -3.048252e+00
## state54:year -2.533342e+00 5.146979e-01  -4.92199770 9.725675e-06 -3.567143e+00
## state55:year -2.464472e+00 5.085935e-01  -4.84566180 1.263773e-05 -3.486012e+00
## state56:year -6.577963e+00 2.349943e-01 -27.99200943 3.266704e-32 -7.049963e+00
##                   CI Upper DF
## legal         1.870588e+01 50
## state1        8.912788e+03 50
## state2        6.626568e+03 50
## state4        1.413878e+04 50
## state5        1.250013e+04 50
## state6        8.207571e+03 50
## state8        7.737375e+03 50
## state9        1.585956e+03 50
## state10       5.195455e+03 50
## state11       1.429646e+04 50
## state12       6.611653e+03 50
## state13       1.382461e+04 50
## state15       1.353016e+03 50
## state16       1.056339e+04 50
## state17       8.117193e+03 50
## state18       7.734494e+03 50
## state19       9.215417e+03 50
## state20       6.242577e+03 50
## state21       6.623687e+03 50
## state22       4.608780e+03 50
## state23       6.973911e+03 50
## state24       7.776717e+03 50
## state25       7.608753e+03 50
## state26       1.186335e+04 50
## state27       8.610873e+03 50
## state28       8.823766e+03 50
## state29       8.215985e+03 50
## state30       5.397735e+03 50
## state31       9.379670e+03 50
## state32       1.743841e+04 50
## state33       6.014376e+03 50
## state34       7.501505e+03 50
## state35       1.017594e+04 50
## state36       9.818789e+03 50
## state37       1.034482e+04 50
## state38       9.819580e+03 50
## state39       6.805466e+03 50
## state40                 NA NA
## state41       8.227791e+03 50
## state42       6.179485e+03 50
## state44       6.892756e+03 50
## state45       9.907082e+03 50
## state46       1.267119e+04 50
## state47       8.493798e+03 50
## state48       2.411527e+03 50
## state49       7.645449e+03 50
## state50       6.076038e+03 50
## state51       7.725540e+03 50
## state53       6.140666e+03 50
## state54       7.182729e+03 50
## state55       7.004651e+03 50
## state56       1.415633e+04 50
## year_fct1971  4.408453e+00 50
## year_fct1972  6.885146e+00 50
## year_fct1973  1.360744e+01 50
## year_fct1974  9.538851e+00 50
## year_fct1975  7.042858e+00 50
## year_fct1976  8.764206e+00 50
## year_fct1977  1.479351e+01 50
## year_fct1978  1.789147e+01 50
## year_fct1979  2.339378e+01 50
## year_fct1980  2.645487e+01 50
## year_fct1981  1.524067e+01 50
## year_fct1982  1.201414e+01 50
## year_fct1983  1.089515e+01 50
## state1:year  -3.821329e+00 50
## state2:year  -1.409593e+00 50
## state4:year  -5.940154e+00 50
## state5:year  -4.215759e+00 50
## state6:year  -2.047384e+00 50
## state8:year  -1.809020e+00 50
## state9:year   1.719786e-01 50
## state10:year -9.924926e-01 50
## state11:year -5.119120e+00 50
## state12:year -2.372671e+00 50
## state13:year -5.654912e+00 50
## state15:year  2.627750e-01 50
## state16:year -3.728938e+00 50
## state17:year -1.957371e+00 50
## state18:year -1.816521e+00 50
## state19:year -3.258494e+00 50
## state20:year -1.055848e+00 50
## state21:year -1.247850e+00 50
## state22:year -2.150220e-01 50
## state23:year -8.272019e-01 50
## state24:year -3.539165e+00 50
## state25:year -2.009134e+00 50
## state26:year -2.952625e+00 50
## state27:year -3.170219e+00 50
## state28:year -2.352588e+00 50
## state29:year -2.051938e+00 50
## state30:year -1.574280e+00 50
## state31:year -2.558353e+00 50
## state32:year -6.671204e+00 50
## state33:year -1.284384e+00 50
## state34:year -2.547569e+00 50
## state35:year -3.009052e+00 50
## state36:year -2.670566e+00 50
## state37:year -3.082763e+00 50
## state38:year -2.861382e+00 50
## state39:year -1.102309e+00 50
## state40:year  7.024805e-02 50
## state41:year -2.056108e+00 50
## state42:year -1.036042e+00 50
## state44:year -1.605962e+00 50
## state45:year -2.909328e+00 50
## state46:year -4.759502e+00 50
## state47:year -2.310039e+00 50
## state48:year -3.913322e-01 50
## state49:year -1.772094e+00 50
## state50:year -1.967785e+00 50
## state51:year -3.635385e+00 50
## state53:year -1.005172e+00 50
## state54:year -1.499541e+00 50
## state55:year -1.442932e+00 50
## state56:year -6.105963e+00 50
#  %>% tidy() %>% 
#  select(term, estimate, std.error) %>% 
#  as.data.frame()

論文の結果とも一致している。

最後に時間トレンドと重みを考慮したDID回帰モデルの数字を確認する。

# 時間トレンドを含める場合、fixed_effects = state + year_fct + state:yearのようには書けない。実行はできるが結果が意図通りではない。
lm_robust(mrate ~ 0 + legal + state + year_fct + state:year,
          data = df2, 
          clusters = state, 
          se_type = "stata",
          weights = pop
          ) %>% 
  tidy()%>% 
  select(term, estimate, std.error, p.value) %>% 
  filter(term == "legal") 
##    term estimate std.error   p.value
## 1 legal 9.649383  4.642592 0.0428203

他にもDID回帰は下記のように発展版がある。画像は立教大学の安藤先生の下記のスライドより引用している。

DID回帰の発展版
DID回帰の発展版

上記に関しては下記が参考になるかも。