library(haven)
library(tidyverse)
library(estimatr)
library(stargazer)
クロスセクションデータ: | 異なる個体について、同時点の観測値を集めたデータ |
時系列データ: | 同じ個体について、異なる時点で観測した値を集めたデータ |
パネルデータ: | 複数の個体について、時系列の観測値を集めたデータ |
time-series data
)可視化してみる
差分の差分 (DID)法ではパネルデータを使う
Difference-in-Differences
)DID
, DD
, dif-in-dif
とも呼ばれる例):
・高速道路無料化する・・・処置群
・高速道路無料化しない・・・統制群
→ 高速道路無料化しない統制群の人々から「不公平だ」という不満が寄せられるため、RCT
による実験がなかなかできない
宮城県 | :ホテル・旅館の宿泊費半額サービスあり | |
山形県 | :ホテル・旅館の宿泊費半額サービスなし | |
都道府県名 | 宿泊者数半額サービス実施「後」(万人 | |
山形県 | 4000 | |
宮城県 | 6000 | |
差 | 2000 | |
・宿泊者数半額サービスが県外からの観光客数を 2000万人増やした
・宮城県と山形県は、処置「宿泊者数半額サービス」以外が同質ではない
→ そもそも宮城県の方が県外からの観光客数が多かっただけなのでは?
宿泊者数(万人) | 宿泊者数(万人) | |||
都道府県名 | サービス実施「前」 | サービス実施「後」 | 差(後−前) | |
宮城県 | 6500 | 6000 | −500 | |
・宿泊者数半額サービスが県外からの観光客数を 500万人減らした
宿泊者数(万人) | 宿泊者数(万人) | |||
都道府県名 | サービス実施「前」 | サービス実施「後」 | 差(後−前) | |
宮城県 | 6500 | 6000 | −500 | |
山形県 | 6000 | 4000 | -2000 | |
差 | 500 | 2000 | 1500 | |
Difference-in-differences: DID
)
を利用する「時点間の差」と「時点間の差」の個体間(=群間)の差を計算する
宮城と山形で after
の個体間を比較 =
「時点間の差」
宮城と山形で before
の個体間を比較 =
「時点間の差」
両者(=群間)の差を求める
\[δ_{DD} = (Y_{宮城, after} − Y_{山形,
after}) - (Y_{宮城, before} − Y_{山形, before})\\= (6000 - 6500) -
(4000 - 6000)\\= -500 + 2000 = 1500\]
「個体間(群間)の差」と「個体間(群間)の差」の時点間の差を計算する
宮城で処置の前後を比較 = 「個体間(群間)の差」
山形で処置の前後を比較 = 「個体間(群間)の差」
両者の時点間の差を求める
Counter factural
)を表しているparallel trend
) の仮定平行トレンドの仮定
・差分の差分法で因果推論を行うために必要な仮定
・「もし処置がなかったら、処置された個体の結果(あるいは処置群の結果の平均)と統制された個体の結果(あるいは統制群の結果の平均)は平行な時間的変化を示す」
・宮城県で東京オリンピック競技が開催された
・山形県では東京オリンピック競技が開催されなかった
・宮城県:緊急事態宣言の対象
・山形県:緊急事態宣言の対象
ceteris paribus
)」なくてもよい変数名 | 内容 | データの種類 |
\(Y_{kt}\) | :結果変数 | パネル |
\(k\) | :処置が適応される単位(例:国、自治体、企業、人など) | |
\(t\) | :時間 | |
\(D_k\) | :個体 \(k\) が処置群に属することを示すダミー変数 \(D = 0: 統制群、D = 1: 処置群\) | クロスセクション |
\(P_t\) | :時間 \(t\) が処置後
(post-treatment ) であることを示すダミー変数 \(P_t = 0: 処置前、P_t = 1: 処置後\) |
クロスセクション |
\[Y_{kt} = α + βD_k + γP_t + δ(D_k × P_t) + e_{kt}\]
\((D_k×P_t)\) の値が 1
になるのは \(D_k=1\) でかつ \(P_k=1\) の時
→ \(D_k=1: 処置群\) でありかつ \(P_t=1: 処置後\) の個体の時
→ 宿泊者数半額サービスを実施した 後の(\(P_t =
1\)) 宮城県(\(D_k=1\))
\(Y_{kt}\) ・・・時間 (\(t\)) と個体 (\(k\)) の値によって変化するパネルデータ
\(D_k\) ・・・ 個体 (\(k\)) に値によってのみ変化するクロスセクションデータ
\(P_t\) ・・・ 時間 (\(t\)) に値によってのみ変化するクロスセクションデータ
\(E[Y_{kt}|D_k = 1, P_t = 1] = α + β + γ + δ\) | (1) |
\(E[Y_{kt}|D_k = 1, P_t = 0] = α + β\) | (2) |
\(E[Y_{kt}|D_k = 0, P_t = 1] = α + γ\) | (3) |
\(E[Y_{kt}|D_k = 0, P_t = 0] = α\) | (4) |
\[E[Y_{kt}|D_k = 1, P_t = 1]- E[Y_{kt}|D_k = 1, P_t = 0] ・・・(5)\]
\[E[Y_{kt}|D_k = 0, P_t = 1] -
E[Y_{kt}|D_k = 0, P_t = 0]・・・(6) \]
・\((5)\) から \((6)\) を引く・・・前後比較の群間比較
\[δ = E[Y_{kt}|D_k = 1, P_t = 1]- E[Y_{kt}|D_k = 1, P_t = 0] \\ - E[Y_{kt}|D_k = 0, P_t = 1] - E[Y_{kt}|D_k = 0, P_t = 0]\]
\[Y_{kt} = α + βD_k + γP_t + δ(D_k × P_t) + e_{kt}\]
cluster-robust standard error
)
を使って推定の不確実性をチェックするlm()
ではなく lm_robust()
を使うDID回帰でクラスタ化した標準誤差を使う理由
・実際の DID 分析では2期間より長い期間を分析することが多い
・2期間より長い期間のパネルデータはある程度長くなければならない
①
ある程度長くなければ平行トレンドの仮定を確認できないから
→ ある程度長いデータがないと前後のトレンドを比較できない
② 処置のタイミングが異なることがあり得るから
・極端に早期に処置した個体と、極端に遅く処置した個体がある場合
→ 全体の期間を含むデータが必要なので、2期間より長いデータが必要
serial correlation
)があり相互に強い相関をもつため☆ 例)若者の死者数が多い州では、その前年も翌年も死者数も多い
・前の年のデータが分かれば次の年のデータも推測ができる
→ 一つ一つの観測値が独立に得られたとはみなせない
・時系列データは相関が強い
→ バラバラに集めたデータと比較すると、時系列データは情報量が少ない
→ 標本サイズを割り引いて考えないと、標準誤差の計算を間違える
→ 標準誤差が小さくなりすぎて、誤って統計的に有意だと判断してしまう
・ある固体内で、結果変数の値をたくさん観測しているとしても、たくさんある観測値のひとつひとつが独立だとは見なせない
・結果変数の値は固体の中でクラスタを形成している
→ それらを割り引いて標準誤差を計算する必要がある
・時系列データでは似ている個体がデータにたくさん含まれていることを想定する必要がある
・似ている個体データ(=クラスタ)が集まっている単位を指定し、標準誤差を計算する
lm()
ではなく lm_robust()
を使うここでは、Card and Krueger (1994)
が分析した、最低賃金の上昇が失業に与える影響の研究を例に考える
1992年4月にアメリカのニュージャージー州 (NJ)で最低賃金が引き上げられた
標準的な経済理論は、最低賃金が引き上げられれば、フルタイム労働者はレイオフされると想定
Counter factural
) として考えることが出来る処置前 | 処置後 | 差 | |
ニュージャージー | 29.7 | 32 | -2.3 |
ペンシルバニア | 31 | 27.2 | 3.8 |
差 | -1.3 | 4.8 | 6.1 |
Card and Krueger (1994)
が分析で使っているオリジナルデータを加工した card.csv
をダウンロードして RProject
フォルダの中にある
data
フォルダに入れる<- read_csv("data/card.csv", na = ".") card
names(card)
[1] "SHEET" "CHAIN" "CO_OWNED" "STATE" "SOUTHJ" "CENTRALJ"
[7] "NORTHJ" "PA1" "PA2" "SHORE" "NCALLS" "EMPFT"
[13] "EMPPT" "NMGRS" "WAGE_ST" "INCTIME" "FIRSTINC" "BONUS"
[19] "PCTAFF" "MEALS" "OPEN" "HRSOPEN" "PSODA" "PFRY"
[25] "PENTREE" "NREGS" "NREGS11" "TYPE2" "STATUS2" "DATE2"
[31] "NCALLS2" "EMPFT2" "EMPPT2" "NMGRS2" "WAGE_ST2" "INCTIME2"
[37] "FIRSTIN2" "SPECIAL2" "MEALS2" "OPEN2R" "HRSOPEN2" "PSODA2"
[43] "PFRY2" "PENTREE2" "NREGS2" "NREGS112"
numeric
)
に変換する$EMPFT = as.numeric(card$EMPFT)
card$EMPPT = as.numeric(card$EMPPT)
card$WAGE_ST = as.numeric(card$WAGE_ST)
card$EMPFT2 = as.numeric(card$EMPFT2)
card$EMPPT2 = as.numeric(card$EMPPT2)
card$WAGE_ST2 = as.numeric(card$WAGE_ST2) card
codebook
)
を読み、分析に使う変数を分かりやすい名前に変えて抜き出すmutate()
で、変数名の変換はrename()
でできるが、変換した変数だけをデータフレームに残したいので、transmute()
を使う<- card %>%
minwage transmute(
state = ifelse(STATE == 1, "NJ", "PA"), # NJ なら state = 1、PA なら state = 0
fulltime_before = EMPFT, # 最低時給上昇前のフルタイム労働者の数
parttime_before = EMPPT, # 最低時給上昇前のパートタイム労働者の数
wage_before = WAGE_ST, # 最低時給上昇前の賃金
fulltime_after = EMPFT2, # 最低時給上昇後のフルタイム労働者の数
parttime_after = EMPPT2, # 最低時給上昇後のパートタイム労働者の数
wage_after = WAGE_ST2, # 最低時給上昇後の賃金
full_prop_before = fulltime_before / (fulltime_before + parttime_before),
full_prop_after = fulltime_after / (fulltime_after + parttime_after)
%>%
) na.omit() # 説明のため、完全ケース分析にする
::datatable(minwage)%>%
DT::formatRound(columns=c('full_prop_before', 'full_prop_after'),
DTdigits=2)
Card and Krueger (1994)
のオリジナルデータ
・Card and Krueger (1994)
が分析で使っているオリジナルデータ
(Stata
用にフォーマットされている)は下のコマンドでダウンロードできる
・下のコマンドを実行する前にパソコンの RProject
フォルダ内にあらかじめ data
フォルダを作っておくこと
・データ (public.dat
) は data
フォルダ内にダウンロードされる
dir.create("tmp")
download.file(url = "http://davidcard.berkeley.edu/data_sets/njmin.zip",
destfile = "tmp/njmin.zip")
unzip("tmp/njmin.zip", exdir = "data")
<- read_table("data/public.dat",
card_org col_names = FALSE, na = ".")
注意:public.dat
は .dta
(Stata
形式) ではない
・このデータには変数名がついておらず、1行目から値が保存されているので、col_names = FALSE
と指定する
・欠測が 「.」
で表されているので、na = "."
を指定
・変数を確認する
names(card_org)
[1] "X1" "X2" "X3" "X4" "X5" "X6" "X7" "X8" "X9" "X10" "X11" "X12"
[13] "X13" "X14" "X15" "X16" "X17" "X18" "X19" "X20" "X21" "X22" "X23" "X24"
[25] "X25" "X26" "X27" "X28" "X29" "X30" "X31" "X32" "X33" "X34" "X35" "X36"
[37] "X37" "X38" "X39" "X40" "X41" "X42" "X43" "X44" "X45" "X46" "X47"
・実際に分析する前に、このオリジナルデータを分析者にとって使いやすくクリーニングする必要がある
・このセクション 4. DID の実例(最低賃金の影響)
では コードブック
(codebook
)
を読み、分析に使う変数を分かりやすい名前に変えたデータセット
(minwage
) を使っている
wage_before
)
と最低時給上昇「後」の賃金 (wage_after
) が
$5.05未満のファーストフード店の割合を計算する%>%
minwage group_by(state) %>%
summarize(before = mean(wage_before < 5.05),
after = mean(wage_after < 5.05),
.groups = "drop") %>%
::kable(digits = 4) knitr
state | before | after |
---|---|---|
NJ | 0.9107 | 0.0034 |
PA | 0.9403 | 0.9552 |
before
) には、どちらの州でも
90%以上の労働者の時給が $5.05未満 (91%: NJ, 94%: PA
)after
)にはニュージャージー (NJ) | :最低賃金の引き上げあり | |
ペンシルベニア (PA) | :最低賃金の引き上げなし | |
full_prop_after
)
を単純に比較してみる%>%
minwage group_by(state) %>%
summarize(fulltime = mean(full_prop_after),
.groups = "drop") %>%
::kable(digits = 4) knitr
state | fulltime |
---|---|
NJ | 0.3204 |
PA | 0.2723 |
・NJ と PA は、処置「最低賃金の上昇」以外が同質ではない
- このままでは因果推論はできない
→ そもそも NJ
の方がもともとフルタイム労働者の割合が高かったのでは?
full_prop_before
:
処置の「前」フルタイム労働者の割合full_prop_after
: 処置「後」のフルタイム労働者の割合%>%
minwage filter(state == "NJ") %>%
summarize(across(.cols = starts_with("full_"), mean)) %>%
::kable(digits = 4) knitr
full_prop_before | full_prop_after |
---|---|
0.2965 | 0.3204 |
full_prop_before
:
処置の「前」フルタイム労働者の割合full_prop_after
: 処置「後」のフルタイム労働者の割合<- minwage %>%
d_full group_by(state) %>%
summarize(across(.cols = starts_with("full_"), mean),
.groups = "drop")
::kable(d_full, digits = 3) knitr
state | full_prop_before | full_prop_after |
---|---|---|
NJ | 0.297 | 0.320 |
PA | 0.310 | 0.272 |
<- d_full %>%
p_did pivot_longer(cols = starts_with("full"),
names_to = "time",
names_prefix = "full_prop_",
values_to = "prop") %>%
mutate(time = factor(time, levels = c("before", "after"),
labels = c("処置前", "処置後"))) %>%
ggplot(aes(x = time, y = prop, group = state)) +
geom_line(aes(color = state)) +
geom_point(aes(shape = state)) +
ylim(0.25, 0.35) +
labs(x = "", y = "フルタイム労働者の平均割合") +
scale_color_discrete(name = "") +
scale_shape_discrete(name = "")
plot(p_did)
%>%
minwage group_by(state) %>%
summarize(across(.cols = starts_with("full_"), mean),
.groups = "drop") %>%
mutate(dif_ba = full_prop_after - full_prop_before) %>%
with(dif_ba[state == "NJ"] - dif_ba[state == "PA"])
[1] 0.06155831
処置前 | 処置後 | 差 | |
ニュージャージー | 29.7 | 32 | -2.4 |
ペンシルバニア | 31 | 27.2 | 3.8 |
差 | -1.3 | 4.8 | 6.1 |
変数名 | 内容 |
\(prop\) | :フルタイム労働者の割合(結果変数) |
\(D\) | :処置群を表すダミー変数 \(D = 0\): 統制群、\(D = 1\): 処置群 |
\(P\) | :処置後を表すダミー変数 \(P_t = 0\): 処置前、\(P_t = 1\): 処置後 |
\(DP\) | :\(D\) と \(P\) の交差項 |
minwage
) は横長
(wide
)で 処置前の結果変数 (full_prop_before
)
と処置後の結果変数 (full_prop_after
)
の値が異なる列にある形式prop
)
に統一する必要がある ::datatable(minwage)%>%
DT::formatRound(columns=c('full_prop_before', 'full_prop_after'),
DTdigits=2)
pivot_longer()
を使って縦長 (long
)
に変換し、必要な変数を作る<- minwage %>%
minwage_long ::select(state, starts_with("full_")) %>%
dplyrpivot_longer(cols = starts_with("full"),
names_to = "time",
names_prefix = "full_prop_",
values_to = "prop") %>%
mutate(D = ifelse(state == "NJ", 1, 0), # NJ なら D = 1, PA なら D = 0
P = ifelse(time == "after", 1, 0)) # 処置後なら P = 1, 処置前なら P = 0
::datatable(minwage_long) %>%
DT::formatRound(columns=c('prop'),
DTdigits=2)
<- lm(prop ~ D * P, data = minwage_long)
did_fit00 tidy(did_fit00) %>%
select(term, estimate) %>%
::kable(digits = 2) knitr
term | estimate |
---|---|
(Intercept) | 0.31 |
D | -0.01 |
P | -0.04 |
D:P | 0.06 |
Angrist and Pischke(2015)
アメリカの法定飲酒年齢 (MLDA)
州ごとによって異なる (50州+ワシントン D.C.)
MLDA: 18, 19, 20, 21
確認したい効果:
MLDA を 21歳 => 18歳 に引き下げると、18-20歳の死亡率が上がるのか?
上で紹介した宮城県と福島県における宿泊者数半額サービスや NJ と PA における最低賃金引き上げの事例とは次の 2 点において異なる
legal
) を作る0 ≦ legal ≦ 1
)MLDA
)が18歳ならlegal = 1
MLDA
)が21歳ならlegal = 0
MLDA
) の引き下げは、legal
の値を大きくするMLDA
の平均処置効果を推定するMLDA
の平均処置効果は、上で作成した legal
を使うことで、legal
の平均処置効果として代替する\[Y_{st} = α + δLEGAL_{st} +{\sum_{k=1}^{50}}{β_kSTATE_{ks}} + {\sum_{j=1971}^{1983}}{γ_jYEAR_{jt}} + e_{st}\]
\(s = {[1, 2, 3,..., 51]}\) | : 州(\(s = 51\)は D.C.) |
\(t = {[1970, 1972, ... , 1983]}\) | : 年(1970年は参照カテゴリ) |
\(Y_{st}\) | : \(t\) 年の \(s\) 州での 18-20歳の死者数(10万人あたり) |
\(STATE_{ks}\) | : \(k = s\) の時に 1 になる変数(D.C.を除く50州のダミー)例)CAなら \(STATE_{CA,s} = 1\)、その他の州は 0 |
\(YEAR_{jt}\) | : \(j = t\) の時に 1 になる変数(1970年を除く13個の年ダミー)例)1971年なら \(YEAR_{1971,s} = 1\)、その他の年は 0 |
\(LEGAL_{st}\) | :処置効果(18歳から20歳までの若者のうち、合法的に飲酒できる人の割合) |
download.file(url = "http://masteringmetrics.com/wp-content/uploads/2015/01/deaths.dta",
destfile = "data/deaths.dta")
Stata
形式のデータなので、haven::read_dta()
で読み込む<- read_dta("data/deaths.dta") MLDA
MLDA
が含む変数の確認names(MLDA)
[1] "year" "state" "legal1820" "dtype" "agegr"
[6] "count" "pop" "age" "legal" "beertaxa"
[11] "beerpercap" "winepercap" "spiritpercap" "totpercap" "mrate"
mrate
mrate
mrate
は次の 4 種類の死に関するデータを含んでいるdtype
によって分類されているすべての死 | All deaths | dtype = 1 |
自動車事故による死 | Motor vehicle accidents | dtype = 2 |
自殺 | Suicide | dtype = 3 |
内蔵疾患による死 | All internal causes | dtype = 6 |
dtype = 1
) のみを扱う注意:(7. Exercise では他の 3 種類の死に関しても取り扱う)
<- MLDA %>%
df5 filter(year <= 1983,
== 2, # 18-21 years old
agegr == 1) %>% # all deaths
dtype mutate(state = factor(state),
year_fct = factor(year))
state
) と年 (year
)
の組み合わせが1つひとつの行を構成する 州-年パネル
(state-year panel
)%>%
df5 summarize(across(.cols = c(state, year),
.fns = n_distinct))
# A tibble: 1 × 2
state year
<int> <int>
1 51 14
mrate
mrate : |
mortality rate: 死亡率(10万人あたりの死者数) |
mrate
の分布を確認する<- ggplot(df5, aes(x = mrate, y = after_stat(density))) +
hist_mrate geom_histogram(color = "black", binwidth = 10) +
labs(x = "死者数(10万人あたり)", y = "確率密度")
plot(hist_mrate)
<- ggplot(df5, aes(x = year, y = mrate, group = state)) +
ts_mrate geom_line() +
labs(x = "年", y = "死者数(10万人あたり)") +
theme(legend.position = "none")
plot(ts_mrate)
D : |
個体の差(処置群と統制群の区別)を表す | state : 州 |
P : |
時間(処置・施策が実施される前後の区別)を表す | yeat_fct : 年 |
D×P : |
処置を受けた後の観測値であることを示す | state:year_fct |
state
を使うyear_fct
を使うyear
をそのまま使うと数値として扱われてしまうので、factor
型に変換した year_fct
を使う)state:year_fct
)
では表現できないMLDA
の変更であるが、MLDA
は18歳, 19歳, 20歳,
21歳のいずれかであるためMLDA
が変更されるタイミングは、州によって異なるためMLDA
変更もあり得るためlegal
) を作るlegal
:
18歳から20歳までの若者のうち、合法的に飲酒できる人の割合
(0 ≦ legal ≦ 1
)MLDA
) の引き下げは、legal
の値を大きくするMLDA
が変更された場合には、その期間に応じて値が調整されるlegal
の値が大きいほど、若者(18-20歳)がアルコールにアクセスしやすいMLDA
を引き下げるという処置(施策)を実行すると、legal
の値が大きくなる\[Y_{st} = α + δLEGAL_{st} +{\sum_{k=1}^{50}}{β_kSTATE_{ks}} + {\sum_{j=1971}^{1983}}{γ_jYEAR_{jt}} + e_{st}\]
\(Y_{st}\): | 結果変数(死亡率:10万人あたりの死者数) | mrate |
\(LEGAL_{st}\): | 18歳から20歳までの若者のうち、合法的に飲酒できる人の割合 | legal |
\(STATE_{ks}\): | 個体の差(処置群と統制群の区別)を表す | state |
\(YEAR_{jt}\): | 時間(処置・施策が実施される前後の区別)を表す | year_fct |
<- df5 %>%
summ.stat select(mrate, legal, state, year)
stargazer(as.data.frame(summ.stat), type = "html")
Statistic | N | Mean | St. Dev. | Min | Max |
mrate | 714 | 136.442 | 39.118 | 52.474 | 331.887 |
legal | 714 | 0.557 | 0.439 | 0.000 | 1.000 |
year | 714 | 1,976.500 | 4.034 | 1,970 | 1,983 |
::datatable(summ.stat) %>%
DT::formatRound(columns=c('mrate', 'legal'),
DTdigits=2)
<- lm(mrate ~ 0
fit_ad0 + legal
+ state
+ year_fct, # year を factor 型に変換した変数
data = df5)
tidy(fit_ad0) %>%
select(term, estimate, std.error, p.value) %>%
filter(term == "legal") %>%
::kable(digits = 2) knitr
term | estimate | std.error | p.value |
---|---|---|---|
legal | 10.8 | 3.14 | 0 |
変数 legal
の効果が 10.80 → これが平均処置効果
18-20歳の人々がアルコールを摂取できない状況から摂取できる状況に変化すると(legal = 0 → legal = 1
)、10万人あたりの死者数が約11人増える
つまり、MLDA
を下げると、酒を飲めるようになった年齢の人たちのなかで死者が増える
Angrist and Pischke (2009) (p.196)
の表5.2 の
“All deaths”
の行の (1) の列にある数値と同じ
Angrist and Pischke (2009) (p.196)
の表5.2の標準誤差・・・
4.59serial correlation
) を考慮していないため州というクラスタを考慮に入れた標準誤差
(cluster-robust standard error
) を計算する
→ stimatr
パッケージの estimatr::lm_robust()
を使う
lm_robust()
でクラスタ標準誤差を得るためには、次の 2
つを指定
cluster
には、クラスタを表す変数を指定
(cluster = state
)
se_type
には、標準誤差を計算する方法を指定
通常は se_type = "CR2"
と指定するが、ここでは
Angrist and Pischke (2009)
と同じ結果を得るために
se_type = "stata"
と 指定
<- lm_robust(mrate ~ 0
fit_ad1 + legal
+ state
+ year_fct,
data = df5,
clusters = state,
se_type = "stata")
tidy(fit_ad1) %>%
select(term, estimate, std.error, p.value) %>%
filter(term == "legal") %>%
::kable(digits = 2) knitr
term | estimate | std.error | p.value |
---|---|---|---|
legal | 10.8 | 4.59 | 0.02 |
legal
の係数の推定値は先ほどとまったく同じ
(10.80)Angrist and Pischke (2009) (p.196)
の表5.2
の “All deaths” の行の (1) の列にある数値に一致DID
分析では、平行トレンドが仮定されているが、分析対象となる個体と期間の両者が多い場合には、仮定を緩めた分析をすることが可能
「傾きが同じ」という平行トレンドの仮定を「傾きが異なる」に緩める
ここでは次の DID 回帰式を推定する
\[Y_{st} = α + δLEGAL_{st} +{\sum_{k=1}^{50}}{β_kSTATE_{ks}} + {\sum_{j=1971}^{1983}}{γ_jYEAR_{jt}} + {\sum_{k=1}^{50}}{θ_k(STATE_{ks}×t)} + e_{st}\]
Angrist and Pischke (2009) (p.196)
の表5.2 の (2)
列では、州ごとの時間トレンドが考慮されているstate × year
(factor
型でなく、数値)で表すstate*year
と入力すると、sate
,
year
, state:year
の 3
つの項が全てモデルに含まれてしまうので、state:year
と入力する<- lm_robust(mrate ~ 0
fit_ad2 + legal
+ state
+ year_fct
+ state:year,
data = df5,
clusters = state,
se_type = "stata")
tidy(fit_ad2) %>%
select(term, estimate, std.error, p.value) %>%
filter(term == "legal") %>%
::kable(digits = 2) knitr
term | estimate | std.error | p.value |
---|---|---|---|
legal | 8.47 | 5.1 | 0.1 |
All deaths
” の行の (2)
の列と同じ結果が得られたState trends
の項目に Yes
と書かれていることが、平行トレンド (parallel trend
)
を考慮したという意味Texas
や California
などと、人口の少ない Wyoming
や Vermont
は異なるlm()
同様、lm_robust()
でも
weights
で重みが指定可能weighted least square: WLS
)
というpop
という変数に記録されているので、weights = pop
と指定する<- lm_robust(mrate ~ 0
fit_ad3 + legal
+ state
+ year_fct,
data = df5,
weights = pop,
clusters = state,
se_type = "stata")
tidy(fit_ad3) %>%
select(term, estimate, std.error, p.value) %>%
filter(term == "legal") %>%
::kable(digits = 2) knitr
term | estimate | std.error | p.value |
---|---|---|---|
legal | 12.41 | 4.6 | 0.01 |
All deaths
” の行の (3)
の列と同じ結果が得られた<- lm_robust(mrate ~ 0
fit_ad4 + legal
+ state
+ year_fct
+ state:year,
data = df5,
weights = pop,
clusters = state,
se_type = "stata")
tidy(fit_ad4) %>%
select(term, estimate, std.error, p.value) %>%
filter(term == "legal") %>%
::kable(digits = 2) knitr
term | estimate | std.error | p.value |
---|---|---|---|
legal | 9.65 | 4.64 | 0.04 |
All deaths
” の行の (4)
の列と同じ結果が得られた・4 つの結果のうち、州ごとのトレンドと人口の違いを考慮した最後のモデル (4) から、次のような結論が得られる
set.seed(2020-07-01)
<- tibble(year = 1990:2009)
sim_df <- rnorm(50, mean = 100, sd = 30)
init <- -3 + 1:50 * 0.2 + rnorm(50, mean = 0, sd = 1)
trend <- sapply(1:50,
outcome function(i) init[i] + trend[i] * 0:19)
colnames(outcome) <- str_c("state_", 1:50)
<- outcome %>%
sim_df as_tibble() %>%
mutate(year = 1990:2009) %>%
pivot_longer(cols = starts_with("state"),
names_to = "state",
names_prefix = "state_",
values_to = "y0")
::datatable(sim_df) %>%
DT::formatRound(columns=c('y0'),
DTdigits=2)
state : |
州を表す変数: (1,2, ..., 50) |
year : |
年を表す変数: (1990, 1991, ..., 2009) |
y0 : |
処置がない場合の結果変数 |
y0
がどのような変化示すか、図示してみるggplot(sim_df, aes(x = year, y = y0, group = state)) +
geom_line()
state
の型を調べるstr(sim_df$state)
chr [1:1000] "1" "2" "3" "4" "5" "6" "7" "8" "9" "10" "11" "12" "13" "14" ...
state
は文字型 (character
)<- sim_df %>%
sim_df mutate(D = ifelse(state %in% as.character(1:25), 1, 0))
<- sim_df %>%
sim_df mutate(P = ifelse(year > 2000, 1, 0))
delta
の処置効果を与えるdelta = 10
とする<- 10
delta <- rep(rnorm(50, mean = delta, sd = 2), each = 20)
effect <- sim_df %>%
sim_df arrange(state, year) %>%
mutate(outcome = case_when(
== 0 ~ y0,
D == 0 ~ y0,
P TRUE ~ y0 + effect
))
::datatable(sim_df) %>%
DT::formatRound(columns=c('y0', 'outcome'),
DTdigits=2)
outcome
) を可視化してみるggplot(sim_df, aes(x = year, y = outcome, group = state)) +
geom_line(aes(color = as.factor(D))) +
scale_color_discrete(name = "", labels = c("統制群", "処置群")) +
labs(x = "年", y = "結果変数") +
geom_vline(xintercept = 2001, linetype = "dotted") +
theme_set(theme_classic(base_size = 10,
base_family = "HiraginoSans-W3"))
DID
回帰を実行してみる<- lm_robust(outcome ~ D * P,
fit_s1 data = sim_df,
clusters = state,
se_type = "stata")
tidy(fit_s1) %>%
filter(term == "D:P") %>%
select(term, estimate, std.error, p.value) %>%
::kable(digits = 2) knitr
term | estimate | std.error | p.value |
---|---|---|---|
D:P | -42.94 | 4.99 | 0 |
state:year
をコントロールして、回帰分析を行ってみる<- lm_robust(outcome ~ D * P
fit_s2 + state:year,
data = sim_df,
clusters = state,
se_type = "stata")
tidy(fit_s2) %>%
filter(term == "D:P") %>%
select(term, estimate, std.error, p.value) %>%
::kable(digits = 2) knitr
term | estimate | std.error | p.value |
---|---|---|---|
D:P | 10.07 | 0.45 | 0 |
<- function(delta = 10,
did_trend_sim n_units = 50,
n_periods = 20,
n_treat = 25,
treat_timing = 10) {
if (n_treat > n_units) stop("n_treat must be smaller than n_units.")
if (treat_timing >= n_periods | treat_timing < 2) stop("treat_timing must be within time periods.")
<- tibble(year = 1:n_periods)
sim_df <- rnorm(n_units, mean = 100, sd = 20)
init <- -n_units/20 + 1:n_units * 0.2 + rnorm(50, mean = 0, sd = 1)
trend <- sapply(1:n_units, function(i) init[i] + trend[i] * ((1:n_periods) - 1))
outcome colnames(outcome) <- str_c("unit_", 1:n_units)
<- outcome %>%
sim_df as_tibble() %>%
mutate(time = 1:n_periods) %>%
pivot_longer(cols = starts_with("unit"),
names_to = "unit",
names_prefix = "unit_",
values_to = "y0") %>%
mutate(D = ifelse(unit %in% as.character(1:n_treat), 1, 0)) %>%
mutate(P = ifelse(time > treat_timing, 1, 0)) %>%
mutate(outcome = case_when(
== 0 ~ y0,
D == 0 ~ y0,
P TRUE ~ y0 + rnorm(n(), mean = delta, sd = 2)
))
## 標準誤差は記録しないので、lm() で推定する
<- lm(outcome ~ D * P, data = sim_df)
f1 <- lm(outcome ~ D * P + unit:time, data = sim_df)
f2 return(c(coef(f1)[4], coef(f2)[4]))
}
<- replicate(1000, did_trend_sim(delta = 10)) res
<- res %>%
res_df t() %>%
as_tibble(.name_repair = c("unique"))
names(res_df) <- c("no_trend", "with_trend")
<- ggplot(res_df, aes(x = no_trend, y = after_stat(density))) +
hist_res1 geom_histogram(color = "lightgrey") +
geom_vline(xintercept = 10, color = "tomato") +
labs(x = "推定値", y = "確率密度", title = "トレンドのコントロールなし")
hist_res1
<- ggplot(res_df, aes(x = with_trend, y = after_stat(density))) +
hist_res2 geom_histogram(color = "lightgrey") +
geom_vline(xintercept = 10, color = "tomato") +
labs(x = "推定値", y = "確率密度", title = "トレンドのコントロールあり")
hist_res2
シミュレーションの結論 ・平行トレンドの仮定が成り立たない場合、トレンドを考慮に入れないと正しい推定ができない
・トレンドをコントロールすることによって、平均処置効果の推定が可能になる
・ただし、これは作成したデータのトレンドが線形だった(そして私たちはそれを事前に知っていた)からうまくいっただけ
・非線形のトレンドがある場合、線形トレンドをコントロールしても、推定はうまくいかないことが予想される
・パネルデータの分析では、どのようなトレンドがあるのか見極めることが重要
Angrist and Pischke (2009)
の例について、以下を実行しなさい
Q1:「自動車事故による死」(Motor vehicle accidents: dtype = 2
)
を使って DID
回帰を実行し、上の表5.2と同じ結果が再現できることを確認しなさい
Q2:「自殺」(Suicide: dtype = 3
)
を使って DID
回帰を実行し、上の表5.2と同じ結果が再現できることを確認しなさい
Q3:「内蔵疾患による死」(All internal causes: dtype = 6
)
を使って DID
回帰を実行し、上の表5.2と同じ結果が再現できることを確認しなさい