「効果検証入門：DIDの章」を読んで気になったのでDIDのテクをいろいろ調べて試してみる
はじめに
前回はこちら
DIDの理論の話は以前にも shinomiya-note.hatenablog.com
回帰分析を利用したDIDと自己相関、標準誤差
DIDを利用する条件として、平行トレンド仮定があります。
DIDは目的変数の期間における時間トレンドが、仮に介入が発生しなかった場合（連続値の場合影響値が0の場合）変わらないという前提です。
特に複数期間を扱うデータだと、自己相関の問題も考えてやるべきでしょう。
また、サンプルによっても時間トレンドが異なると考えられます。これも考慮すべきでしょう。
DIDを複数期間、介入を連続値と二値、自己相関の検討、サンプルの効果（混合効果）といろいろ条件を変えてやってみます。またCausalImpactも試してみます。
データについて
高知工科大学の授業ページが本書のことも踏まえて解説しています。これはいい。参考にさせていただきます。
法定飲酒年齢（MLDA）の変更が18歳から20歳までの若者の死に与える影響の分析
このデータは、州 (state) と年 (year) の組み合わせが1つひとつの行を構成する 州-年パネルデータ (state-year panel)
この分析で考える処置（介入）はMLDAの変更であるが、MLDA は18歳, 19歳, 20歳, 21歳のいずれかである。 また、MLDA が変更されるタイミングは、州によって異なる。さらに、（少なとも理論的には）複数回の MLDA 変更も可能である。
そこで、このデータセットには、処置変数として legal が用意されている。 この変数は、特定の年の特定の州で、18歳から20歳までの人口のうち何割が合法的に飲酒できるかを表す。例えば、 t 年のs 州でのMLDAが21歳だとすると、20歳以下で合法的に飲酒できる者はいないので、この変数の値は0になる。MLDAが18歳なら、18歳以上の全員が合法的に飲酒できるので、この変数の値は1になる
つまり
- 介入変数は連続値で0～1を取る
- 各サンプルは州番号
- 期間は複数
library("tidyverse") library("haven") # library("broom") library("estimatr") #Stata 形式のデータなので、haven::read_dta() で読み込む。 url <- ("http://masteringmetrics.com/wp-content/uploads/2015/01/deaths.dta") MLDA <- read_dta(url) # 分析に必要なのはデータの一部なので、その部分を抜き出す myd <- MLDA %>% dplyr::filter(year <= 1983, agegr == 2, # 18-21 years old dtype == 1) %>% # all deaths mutate(state = factor(state), year_fct = factor(year)) myd %>% summarize(across(.cols = c(state, year), .fns = n_distinct)) # 死亡率の時系列変化を可視化 ts_mrate <- ggplot(myd, aes(x = year, y = mrate, group = state)) + geom_line() + labs(x = "年", y = "死者数（10万人あたり）") + theme(legend.position = "none") plot(ts_mrate
二値データに要約して試してみる
本来は正しくないやり方ですが、参考にむりやり期間とlegalを二値データに要約して試してみます
# 二値データに変換 # 1970年に0,1983年に0.5以上のstateに絞り込む pre <- myd %>% filter(year == 1970 & legal == 0) %>% select(state) post0 <- myd %>% filter(year == 1983 & legal == 0) %>% select(state) %>% mutate(z=0) post1 <- myd %>% filter(year == 1983 & legal > 0.5) %>% select(state) %>% mutate(z=1) target_state <- pre %>% inner_join(bind_rows(post0,post1), by='state') %>% select(state,z) myd01 <- myd %>% inner_join(target_state, by='state')%>% filter(year == 1970 | year == 1983) %>% mutate(time =if_else(year==1970,0,1),year_fct=factor(year),z_f=factor(z)) # plot myd01_means <- myd01 %>% group_by(year,year_fct,time,z_f,z) %>% summarise(mrate_mean = mean(mrate)) ts_mrate01 <- ggplot(myd01_means, aes(x=year_fct,y=mrate_mean,group=z_f)) + geom_line(aes(color=z_f)) + geom_point(aes(shape=z_f)) + labs(x = "年", y = "死者数（10万人あたり）") plot(ts_mrate01)
# 二値データで推定 fit_ap01l <- lm(mrate ~z + time + z:time, data = myd01) tidy(fit_ap01l) %>% select(term, estimate, std.error) %>% filter(term == "z:time") %>% knitr::kable(digits = 2) |term | estimate| std.error| |:------|--------:|---------:| |z:time | 11.55| 16.97|
18-20歳の人々がアルコールを摂取できない状況から摂取できる状況に変化すると、10万人あたりの死者数が11.55人増えそうです。
全体の効果
要約していないデータでやってみます。介入は連続値、期間も複数です。
# 回帰分析を利用したDID fit_ap <- lm(mrate ~ legal + year_fct, data = myd) tidy(fit_ap) %>% select(term, estimate, std.error) %>% filter(term == "legal") %>% knitr::kable(digits = 2) |term | estimate| std.error| |:-----|--------:|---------:| |legal | -3.53| 3.38|
二値の要約データでの結果と異なり、マイナスになってしまいました。直感的には反する結果ですね。
地域の特徴を入れる
各州の特徴をダミー扱いの共変量（固定効果）として導入して推定します。
# 回帰分析を利用したDID # stateを共変量に（固定効果） fit_ap0 <- lm(mrate ~ legal + state + year_fct, data = myd) tidy(fit_ap0) %>% select(term, estimate, std.error) %>% filter(term == "legal") %>% knitr::kable(digits = 2) |term | estimate| std.error| |:-----|--------:|---------:| |legal | 10.8| 3.14|
各州の状況を考慮すると10.8人になります。
自己相関を考慮する（クラスター標準誤差）
期間毎にデータを取るパネルデータは同一サンプルから複数データを取るので、自己相関が想定されます。
- 自己相関を持つと誤差項の分散が小さくなる
- 回帰分析のパラメータの標準誤差は誤差項の分散を利用すため、パラメータの標準誤差が過小になる
そのため、クラスター（かつロバストな）標準誤差を使います。本書ではmiceaddsパッケージの:lm.cluster関数を使っていましたが、estimatrパッケージの方が使い勝手がよさそう。
# stateを共変量に（固定効果、クラスターロバスト標準誤差） fit_ap1 <- lm_robust(mrate ~ 0+ legal + state + year_fct, data = myd, clusters = state, se_type = "stata") tidy(fit_ap1) %>% select(term, estimate, std.error) %>% filter(term == "legal") %>% knitr::kable(digits = 2) |term | estimate| std.error| |:-----|--------:|---------:| |legal | 10.8| 4.59|
推定値は変わりませんが、過小評価されていたパラメータの標準誤差が大きくなりました（3.14→4.59）。確かに過小評価されていたようです*1。
誤差項の分散とパラメータの標準誤差の関係
回帰分析においては以下のような関係があるからです \begin{align} & d_{i} = y_{i} − (\hat{a} +\hat{b}x_{i})・・・残差 \\ \\ & s^{2} = \frac{\sum d^{2}_{i} }{n-2} ・・・誤差項の分散 \\ & \\ & s.e.(\hat{b}) = \sqrt{ \frac{s^{2}}{\sum ( x_{i}-\bar{x})^{2}} } ・・・パラメータの標準誤差 \\ \end{align}
ロバスト標準誤差
回帰分析における条件として均一分散の仮定があります。がからなずしも均一ではありません。それを考慮してより頑健な標準誤差を算出しています。
クラスター標準誤差
クラスター間（今回はstate間）の誤差項の相関を考慮した標準誤差。
サンプル間の期間トレンドの違いを考慮する（トレンド仮定を緩める）
DID 分析では、平行トレンドが仮定されているが、分析対象となる個体と期間の両者が多い場合には、仮定を緩めた分析をすることができる。個体ごとの時間トレンドの違いを説明する変数を回帰分析に加えることにより、トレンドが完全に平行でなくても、平均置効果を推定することができるようになるのである。
stateによってyearの効果が異なる、とする交互作用項をさらに加えます。
fit_ap2 <- lm_robust(mrate ~ 0+ legal + state + year_fct + state:year, data = myd, clusters = state, se_type = "stata") tidy(fit_ap2) %>% select(term, estimate, std.error) %>% filter(term == "legal") %>% knitr::kable(digits = 2) |term | estimate| std.error| |:-----|--------:|---------:| |legal | 8.47| 5.1|
8.47人。stateの固定効果のみより少し少ないですね。
階層モデルでサンプル間の期間トレンドの違いを考慮する
”個体ごとの時間トレンドの違いを説明する変数を回帰分析に加える”という意味であれば、stateによって異なる回帰直線を仮定し、その代表値を固定効果とする階層モデルでもよいんじゃないでしょうか？（誰か教えてください）
# stateを共変量に（HLM階層モデル） hlm_ap0 <- lmer(mrate ~ legal + (legal|state)+ year_fct , data = myd) # 固定効果 fixef(hlm_ap0)[2] legal 8.832869
8.83人。ほぼ同じですね。
さらに人口の違いを考慮
さらにstaeごとのトレンドと人口の違いを考慮し、重みづけしてやります。
fit_ap4 <- lm_robust(mrate ~ 0 + legal + state + year_fct + state:year, data = myd, weights = pop, clusters = state, se_type = "stata") tidy(fit_ap4) %>% select(term, estimate, std.error) %>% filter(term == "legal") %>% knitr::kable(digits = 2) |term | estimate| std.error| |:-----|--------:|---------:| |legal | 9.65| 4.64|
9.65人。
結果
様々なモデル試してみた結果、18-20歳の死者は10万人あたり約8～10人増えると言えそうですね。こうした様々なモデルを試して結果に大きな違いが出ないことが重要です。
*1:因果推論においてパラメータの標準誤差自体は検定結果が有意差かどうかに影響するのみ。またサンプル数が大きくなれば検定自体は有意になりやすくなるので、サンプル数が大きければそこまで気にしなくても問題ないか
