R による実証分析: 第12章

学部：鶴岡ゼミ

星野匡郎, 田中久稔, & 北川梨津. (2023). R による実証分析: 回帰分析から因果分析へ. 株式会社 オーム社.

12章を担当した際の資料

Transcript

1. 第 12 章 : パネルデータ分析 星野匡郎 , 田中久稔 , & 北川梨津 . (2023). R による実証 分析 : 回帰分析から因果分析へ . 株式会社 オーム社 . 花澤楓 最終更新 : 2024-01-18 21:22:50.958964
2. 12.1 パネルデータ 2 / 48
3. パネルデータとクロスセクションデータ パネルデータ: 同一の個体を複数の時点で観察し収集したデータ クロスセクションデータ: 各個体についてある同一の1 時点での情報を収集したデータ Table: パネルデータ ID 2018 2019 2020 2021 A 10 15 18 22 B 20 25 30 35 Table: クロスセクションデータ ID 2020 A 18 B 30 3 / 48
4. パネルデータの特徴 変数X について2 つのインデックスが存在する。 個人を表す （ ） 観測時点を表す （ ） 変数 について3 つのタイプが考えられる。 1. 個人 と時点 の両方に依存: 年収、家賃、労働時間など 2. 個人間で異なるが、時点では固定: 性別、母国語、個人の先天的資質 3. 時点によって異なり、個人間では共通: 消費税率、パンデミック化の行動制限（母集団全体に一律 に影響する変数） が観察できない変数であっても、上記のいずれかに分類される。 i i = 1, … , n t t = 1, … , T Xi,t i t Xi,t 4 / 48
5. パネルデータにおける線形回帰モデル 回帰モデルは以下のようにかける。 ここで、誤差項 についても、 として、各タイプに対応する誤差項に分解できる。 もし、 が各説明変数と相関しないのであれば、OLS でバイアスなしの係数を推定可能 （Pooled OLS という）。 Y i,t = β 0 + X (1) i,t β 1 + ⋯ + X (k) i,t β k  タイプ1 + Z (1) i α 1 + ⋯ + Z (l) i α l  タイプ2 + W 1 t γ 1 + ⋯ + W (m) t γ m  タイプ3 + ε i,t for i = 1, … , n, t = 1, … , T . ε i,t ε i,t = ϵ i,t + e i + u t ϵ i,t + e i + u t 5 / 48
6. パネルデータの種類 パネルデータの場合、標本全体の大きさは となる。大数の法則や中心極限定理を適用するた めに、標本の大きさを無限に大きくするには、 1. のみを無限に大きくする (i.e., ): 公的調査統計 2. のみを無限に大きくする (i.e., ): 金融データ Big Data 3. と のどちらも無限に大きくする (i.e., ) どの枠組みを採用するかによって、選択可能な分析方法も変化する。 Unbalanced panel data: 個人ごとに観察できる時期や期間が異なる場合 Balanced panel data: 全ての個人について同一の時期・期間で観察できる場合 本書ではタイプ1 & Balanced panel data について考える。 n × T n n → ∞, T < ∞ T n < ∞, T → ∞ ⊆ n T n → ∞, T → ∞ 6 / 48
7. 12.2 パネルデータ回帰分析 7 / 48
8. モデル 説明変数の数を限定し、モデルを簡単化する。 例： : 地域 の年 における新規求人数 : 人口増加数 真の を正しく推定するには と が無相関である必要。ここで は、 年 に依存せず新規求人数に影響を与える、観察できない地域 固有の要因 「地域の潜在的能力」「交通の利便性」など パネルデータが利用可能であれば、IV を使用せずとも を正しく推定することができる。 Yi,t = β0 + Xi,t β1 + Zi α + Wt γ + ϵi,t + ei + ut for i = 1, … , n, t = 1, … , T . Yi,t i t Xi,t β1 Xi,t ϵi,t + ei + ut ei t i β1 8 / 48
9. 固定効果モデル 時点 について不変の 固有の効果と、定数項 をまとめる。 時点 にのみ依存し、個人間で共通の効果を表す項をまとめる。 よって、回帰モデルは ここで、指示関数を とすると、回帰モデルは t i β0 ϕi = β0 + Zi α + ei t ψt = Wt γ + ut Yi,t = Xi,t β1 + ϕi + ψt + ϵi,t for i = 1, … , n, t = 1, … , T . Ii,j = 1{i = j}, δt,k = 1{t = k} Y i,t = X i,t β 1 + n ∑ j=1 I i,j ϕ j + T ∑ k=1 δ t,k ψ k + ϵ i,t for i = 1, … , n, t = 1, … , T . 9 / 48
10. 固定効果モデル においては、 被説明変数: 説明変数: として重回帰モデルと解釈できる。 式(12.2) は、固定効果モデル（fixed effects model ）と呼ばれる。 特に、個人効果と時点効果の両方を含むモデルを2 方向固定効果モデル（two-way fixed effects model ）という。 Yi,t = Xi,t β1 + n ∑ j=1 Ii,j ϕj + T ∑ k=1 δt,k ψk + ϵi,t for i = 1, … , n, t = 1, … , T . (12.2) Yi,t Xi,t , Ii,1 , … , Ii,n , δt,1 , … , δt,T 10 / 48
11. 固定効果モデル : 注意点 1 係数 と について、「相対的大きさ」のみ推定可能であること。 かつ であることに注意すると、任意の定数 について、 が成立する。これは、 と が識別できないことを意味する。よって、値を基準化することが必 要。例えば、 を説明変数から外す など。 (ϕ1 , … , ϕn ) (ψ, … , ψT ) ∑ n j=1 Ii,j = 1 ∑ T k=1 δt,k = 1 c n ∑ j=1 Ii,j ϕj + T ∑ k=1 δt,k ψk = n ∑ j=1 Ii,j (ϕj + c) + T ∑ k=1 δt,k (ψk − c) ϕi ψt ϕ1 = 0 δt,1 11 / 48
12. 固定効果モデル : 注意点 2,3 説明変数として含むことができるのは、タイプ1 （個人 と時点 の両方に依存）に該当する変 数のみ。 もし興味のある変数が時点間で不変、もしくは個人間で一定の場合、固定効果に含まれるため。 個人効果 について大標本理論が使えない 「 が増え、 は一定」というデータを考えているため、個人効果パラメータ はn と同じだけ増えていく 各 について、 となるタイミングは が成立する 回のみ 各個人効果 を推定するのに、それぞれ 個のデータしか使えない は一定という設定のもとでは、 の推定値に対して、大数の法則や中心極限定 理を適用できる十分なデータがない 各時点効果パラメータ については、 個のデータより推定されるので、 の係数 と同様に 検定などを実行できる。 i t ϕj n T (ϕ1 , … , ϕn ) j Ii,j = 1 {i = j} T ϕ j T T (ϕ1 , … , ϕn ) ψk n Xi,t β1 t 12 / 48
13. 固定効果モデルの推定 簡単化のため、時点固有の効果をモデルから落として考える。 各 における全 期間の各変数の平均値を、 とする。ここで、 は に依存しないので、 期間平均は となる。 Yi,t = Xi,t β1 + n ∑ j=1 Ii,j ϕj + ϵi,t i = 1, … , n, t = 1, … , T . (12.3) i T ¯ ¯ ¯ ¯ ¯ Yi = T ∑ t=1 Yi,t , ¯ ¯ ¯ ¯ ¯ ¯ Xi = T ∑ t=1 Xi,t , ¯ ¯ ¯ ¯ ϵi = T ∑ t=1 ϵi,t . 1 T 1 T 1 T ∑ n j=1 Ii,j ϕj t T ∑ n j=1 Ii,j ϕj (= ϕi ) 13 / 48
14. 固定効果モデルの推定 式（12.3 ）の両辺について に関する平均をとる。 この式を式（12.3 ）から辺々引くことで となる。式（12.4 ）において、個人効果がモデルから除去されていることがわかる。 この式からモデルを推定することを、within 推定と呼ぶ。 t ¯ ¯ ¯ ¯ ¯ Yi = ¯ ¯ ¯ ¯ ¯ ¯ Xi β1 + n ∑ j=1 Ii,j ϕj + ¯ ¯ ¯ ¯ ϵi i = 1, … , n. Y i − ¯ ¯ ¯ ¯ ¯ Y i      被説明変数 = {X i,t − ¯ ¯ ¯ ¯ ¯ ¯ X i      説明変数 }β 1 + ϵ i,t − ¯ ¯ ¯ ¯ ϵ i    誤差項 i = 1, … , n. (12.4) 14 / 48
15. 固定効果モデルの推定 重回帰モデル: within 推定: 上記2 つのモデルから得られる推定値 は等しくなることを確かめる。重回帰モデルにおいて、係 数パラメータ のもとでの平均2 乗誤差MSE は、 Yi,t = Xi,t β1 + n ∑ j=1 Ii,j ϕj + ϵi,t i = 1, … , n, t = 1, … , T . (12.3) Y i − ¯ ¯ ¯ ¯ ¯ Y i = {X i,t − ¯ ¯ ¯ ¯ ¯ ¯ X i }β 1 + ϵ i,t − ¯ ¯ ¯ ¯ ϵ i i = 1, … , n. (12.4) ^ β 1 (β 1 , ϕ 1 , … , ϕ n ) MSE = n ∑ i=1 T ∑ t=1 (Yi,t − Xi,t β1 − n ∑ j=1 Ii,j ϕj) 2 1 nT 15 / 48
16. に関する一階条件より、 ϕj n ∑ i=1 T ∑ t=1 2 (Yi,t − Xi,t β1 − n ∑ j=1 Ii,j ϕj) = 0 n ∑ i=1 T ∑ t=1 n ∑ j=1 Ii,j ϕj = n ∑ i=1 T ∑ t=1 (Yi,t − Xi,t β1 ) n ∑ i=1 ϕi = n ∑ i=1 ( ¯ ¯ ¯ ¯ ¯ Yj − ¯ ¯ ¯ ¯ ¯ ¯ Xj β1) (∵ n ∑ j=1 Ii,j ϕj = ϕi) ∴ ^ ϕ j = ¯ ¯ ¯ ¯ ¯ Y j − ¯ ¯ ¯ ¯ ¯ ¯ X j ^ β 1 ∀j 1 nT 1 nT 1 nT 1 n 1 n 16 / 48
17. また、 を使って、 に関する一階条件より (∑(Yi − ¯ Y )(Xi − ¯ X) = ∑(Yi − ¯ Y )Xi) β1 0 = n ∑ i=1 T ∑ t=1 Xi,t (Yi,t − Xi,t ^ β 1 − n ∑ j=1 Ii,j ^ ϕ j ) 0 = n ∑ i=1 T ∑ t=1 Xi,t (Yi,t − Xi,t ^ β 1 − ^ ϕ i ) 0 = n ∑ i=1 T ∑ t=1 Xi,t (Yi,t − Xi,t ^ β 1 − ¯ ¯ ¯ ¯ ¯ Yj + ¯ ¯ ¯ ¯ ¯ ¯ Xj ^ β 1 ) 0 = n ∑ i=1 T ∑ t=1 X i,t ({Y i,t − ¯ ¯ ¯ ¯ ¯ Y j } − {X i,t − ¯ ¯ ¯ ¯ ¯ ¯ X j } ^ β 1 ) 0 = n ∑ i=1 T ∑ t=1 {Xi,t − ¯ ¯ ¯ ¯ ¯ ¯ Xj } ({Yi,t − ¯ ¯ ¯ ¯ ¯ Yj } − {Xi,t − ¯ ¯ ¯ ¯ ¯ ¯ Xj } ^ β 1 )  within 推定量の 一階条件と等しい ⇔ ^ β 1 はwithin 推定量と一致 1 nT 1 nT 1 nT 1 nT 1 nT 17 / 48
18. R による固定効果モデルの推定 prefdata <- read_csv("prefecture.csv") datasummary_skim(prefdata, output = "kableExtra") Unique (#) Missing (%) Mean SD Min Median Max year 23 0 2008.0 6.6 1997.0 2008.0 2019.0 suicide 214 0 22.5 4.9 12.4 22.2 44.6 unemp 66 0 3.8 1.2 1.1 3.7 8.4 pref: 都道府県名。個体 year: 年度。時点 suicide: 自殺死亡率（単位: 自殺者/10 万人）。被説明変数 unemp: 完全失業率（単位: % ）。説明変数 i t Yit Xit 18 / 48
19. まず OLS lm(suicide ~ -1 + unemp + pref, data = prefdata) |> # -1 で定数項を含めない modelsummary(stars = TRUE, coef_omit = c(-1) , gof_omit = "IC|Adj|F|RMSE|Log", output = "kableExtra") (1) unemp 3.067*** (0.086) Num.Obs. 1081 R2 0.988 + p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001 19 / 48
20. within 推定 : 推定値が OLS と同じことを確認 prefdata |> group_by(pref) |> mutate( suicidebar = mean(suicide), unempbar = mean(unemp), suicide2 = suicide - suicidebar, unemp2 = unemp - unempbar ) |> (\(x) lm(suicide2 ~ -1 + unemp2, data = x))() |> modelsummary(stars = TRUE, coef_omit = c(-1) , gof_omit = "IC|Adj|F|RMSE|Log", output = "kableExtr (1) unemp2 3.067*** (0.084) Num.Obs. 1081 R2 0.551 + p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001 20 / 48
21. 2 方向固定効果モデルの推定 都道府県効果、時点効果をコントロールして分析してみる。 feols(suicide ~ unemp | pref + year, data = prefdata) |> # `|` の後に固定効果の項を入れる modelsummary(stars = TRUE, coef_omit = c(-1) , gof_omit = "IC|Adj|F|RMSE|Log", output = "kableExtr (1) unemp 0.771* (0.299) Num.Obs. 1081 R2 0.871 R2 Within 0.023 Std.Errors by: pref + p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001 21 / 48
22. 各モデルの比較 OLS within 推定 2FE unemp 3.067 3.067 0.771* (0.086) (0.084) (0.299) Num.Obs. 1081 1081 1081 R2 0.988 0.551 0.871 R2 Within 0.023 Std.Errors by: pref + p < 0.1, p < 0.05, p < 0.01, p < 0.001 unemp の係数の推定量は、都道府県固定効果のモデルの結果と比べて大きく小さくなっているこ とがわかる。 時点効果をモデルに含めないことで上方バイアスがかかっていた 22 / 48
23. 固定効果で比較 以下の4 つのモデルを比較する。 (1): 都道府県効果も時点効果もコントロールしない (2): 時点効果のみをコントロール (3): 都道府県効果のみをコントロール (4): 2 方向固定効果モデル models <- list( "POLS" = feols(suicide ~ unemp, data = prefdata), # (1) "prefFE" = feols(suicide ~ unemp | pref, data = prefdata), # (2) "yearFE" = feols(suicide ~ unemp | year, data = prefdata), # (3) "TFE" = feols(suicide ~ unemp | pref + year, data = prefdata) # (4) ) 23 / 48
24. 結果 各モデルの比較 POLS prefFE yearFE TFE unemp 2.017*** 3.067*** 0.513** 0.771* (0.109) (0.165) (0.146) (0.299) Num.Obs. 1081 1081 1081 1081 R2 0.243 0.727 0.484 0.871 R2 Within 0.551 0.015 0.023 Std.Errors IID by: pref by: year by: pref + p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001 24 / 48
25. 12.3 差の差分析 (difference-in-differences; DID) 25 / 48
26. モデル 2 期間のパネルデータ: 関心のある処置変数: 処置群に属する場合、時点 と時点 の間に処置を受ける。 「 のときの潜在的な結果変数」: 実際に観察される結果変数は、 関心のある尺度: 「 における処置群に対する平均処置効果 」 t ∈ {1, 2} Di ∈ {0, 1} t = 1 t = 2 Di = d Yi,t (d) Yi,t = Di Yi,t (1) + (1 − Di )Yi,t (0) t = 2 ATET2 ATET2 = E[Yi,2 (1) − Yi,2 (0) ∣ Di = 1] = E[Yi,2 (1) − Yi,1 (1) ∣ Di = 1] − E[Yi,2 (0) − Yi,1 (1) ∣ Di = 1] = E[Yi,2 − Yi,1 ∣ Di = 1]  処置郡におけるY i,2 とY i,1 の 差の 平均値 − E[Yi,2 (0) − Yi,1 (1) ∣ Di = 1] (12.5) 26 / 48
27. 仮定 1. である全ての について、 : no anticipatory effect 2. : parallel trends 仮定1: 「処置を実際に受ける以前 に処置の影響はない」 このとき、 となる 例: 政府「来週から飲食店の営業時間を短縮します」 一部の人「会食の予定を前倒しだ」 仮定2: 「処置群が仮に処置を受けなかった場合の結果変数の変動 = 対照群の結果変数の変動」 と仮定2 を書き換えることができる。 Di = 1 i Yi,1 (1) − Yi,1 (0) = 0 E[Yi,2 (0) − Yi,1 (0) ∣ Di = 1] = E[Yi,2 (0) − Yi,1 (0) ∣ Di = 0] t = 1 Yi,1 = Yi,1 (1) → E[Yi,2 (0) ∣ Di = 1]                    データから直接推定できない = E[Yi,1 ∣ Di = 1] + E[Yi,2 − Yi,1 ∣ Di = 1] 27 / 48
28. ATET 2 = E[Y i,2 − Y i,1 ∣ D i = 1] − E[Y i,2 (0) − Y i,1 (1) ∣ D i = 1] = E[Yi,2 − Yi,1 ∣ Di = 1] − E[Yi,2 (0) − Yi,1 (0) ∣ Di = 1] (∵ 仮定1) = E[Y i,2 − Y i,1 ∣ D i = 1] − E[Y i,2 (0) − Y i,1 (0) ∣ D i = 0] (∵ 仮定2) = E[Yi,2 − Yi,1 ∣ Di = 1]  処置郡におけるY i,2 とY i,2 の 差の 期待値 − E[Yi,2 − Yi,1 ∣ Di = 0]  対照郡におけるY i,2 とY i,2 の 差の 期待値 28 / 48
29. の推定 グループごとの期待値をそれぞれの標本平均に置き換えて推定する。 データ ここで、 AT ET 2 {(Yi,1 , Yi,2 , Di ) : i = 1, … , n} ˆ ATET2 = ( ^ Y 2,D=1 − ^ Y 1,D=1 ) − ( ^ Y 2,D=0 − ^ Y 1,D=0 ) ^ Y t,D=1 = , ^ Y t,D=0 = ∑ n i=1 Di Yi,t ∑ n i=1 Di ∑ n i=1 (1 − Di )Yi,t ∑ n i=1 1 − Di 29 / 48
30. 回帰分析による推定 この重回帰モデルの最小二乗法から係数の推定値 を得る。結果としては、 ここで、 . 従って、 Yi,t = (δt,2 Di )β1 + n ∑ i=1 Ii,j ϕj + δt,2 ψ + ϵi,t , i = 1, … , n, t = 1, 2 ( ^ β 1 , ^ ϕ 1 , … , ^ ϕ n , ^ ψ) ^ Y 1,D=1 = ^ ϕ D=1 (D = 1, t = 1) ^ Y 2,D=1 = ^ β 1 + ^ ϕ D=1 + ^ ψ (D = 1, t = 2) ^ Y 1,D=0 = ^ ϕ D=0 (D = 0, t = 1) ^ Y 2,D=0 = ^ ϕ D=0 + ^ ψ (D = 0, t = 2) ^ ϕ D=1 = ∑ n i Di ^ ϕ i / ∑ n i Di , ^ ϕ D=0 = ∑ n i (1 − Di ) ^ ϕ i / ∑ n i (1 − Di ) ^ Y 2,D=1 − ^ Y 1,D=1 = ^ β 1 + ^ ψ, ^ Y 2,D=0 − ^ Y 1,D=0 = ^ ψ 30 / 48
31. ここで、残差を とする。一階条件より、 式(1) より、以下の式を得る。 MSE(β 1 , ϕ 1 , … , ϕ n , ψ) = n ∑ i=1 2 ∑ t=1 (Y i,t − (δ t,2 D i )β 1 − n ∑ i=1 I i,j ϕ j − δ t,2 ψ) 2 1 2n ^ e i,t = Y i,t − (δ t,2 D i ) ^ β 1 − ^ ϕ i − δ t,2 ^ ψ 0 = n ∑ i=1 2 ∑ t=1 (δt,2 Di )^ ei,t = n ∑ i : D i =1 ^ ei,2 0 = n ∑ i=1 2 ∑ t=1 I i,j ^ e i,t = 2 ∑ t=1 ^ e j,t 0 = n ∑ i=1 2 ∑ t=1 δt,2 ^ ei,t = n ∑ i=1 ^ ei,2 (1) (2) (3) 1 n 1 n 1 n 1 n 1 n 1 n ^ Y 2,D=1 = ^ β 1 + ^ ϕ D=1 + ^ ψ 31 / 48
32. 処置群の割合を、 、対照群の割合を、 とする。 式(3) より、 式(2) について、 のとき、 pT = ∑ n i=1 Di /n p C = ∑ n i=1 (1 − Di )/n 0 = p T ( ^ Y 2,D=1 − ^ β 1 − ^ ϕ D=1 − ^ ψ) + p C ( ^ Y 2,D=0 − ^ ϕ D=0 − ^ ψ) 0 = pC ( ^ Y 2,D=0 − ^ ϕ D=0 − ^ ψ) ∴ ^ Y 2,D=0 = ^ ϕ D=0 + ^ ψ 0 = 2 ∑ t=1 (Yj,t − (δt,2 Dj ) ^ β 1 − ^ ϕ j − δt,2 ^ ψ) 1 n Dj = 0 0 = Yj,1 + Yj,2 − 2 ^ ϕ j − ^ ψ 0 = ^ Y 1,D=0 + ^ Y 2,D=0 − 2 ^ ϕ j − ^ ψ ( 平均をとった) 32 / 48
33. より、 のとき、 となる。 を代入して、 以上の結果より、DID 分析による の推定量は に一致することがわかる。 ^ Y 2,D=0 = ^ ϕ D=0 + ^ ψ ^ Y 1,D=0 = 2 ^ ϕ D=0 + ^ ψ − ^ ϕ D=0 − ^ ψ = ^ ϕ D=0 Dj = 1 0 = Y j,1 + Y j,2 − ^ β 1 − 2 ^ ϕ j − ^ ψ 0 = ^ Y 1,D=1 + ^ Y 2,D=1 − ^ β 1 − 2 ^ ϕ j − ^ ψ ( 平均をとった) ^ Y 2,D=1 = ^ β 1 + ^ ϕ D=1 + ^ ψ ^ Y 1,D=1 = ^ β 1 + 2 ^ ϕ D=1 + ^ ψ − ^ β 1 − ^ ϕ D=1 − ^ ψ = ^ ϕ D=1 ATET2 ^ β 1 ^ Y 1,D=1 = ^ ϕ D=1 , ^ Y 2,D=1 = ^ β 1 + ^ ϕ D=1 + ^ ψ ^ Y 1,D=0 = ^ ϕ D=0 , ^ Y 2,D=0 = ^ ϕ D=0 + ^ ψ ∴ ( ^ Y 2,D=1 − ^ Y 1,D=1 ) − ( ^ Y 2,D=0 − ^ Y 1,D=0 ) = ^ β 1 33 / 48
34. DID 分析 説明変数を追加したモデルが使用される。 仮定1,2 が成立するならば、説明変数を追加しなくとも係数の推定量に影響はないが、追加すること で係数 の標準誤差を小さくできる。 「処置群ダミー と個人属性の交差項」を含むことで、平行トレンドの成否が何らかの属性に 依存する場合に対応できる Yi,t = (δt,2 Di )β1 + X (2) i,t β2 + ⋯ + X (k) i,t βk + n ∑ j=1 Ii,j ϕj + δt,2 ψ + ϵi,t for i = 1, … + n, t = 1, 2. β 1 D i 34 / 48
35. 平行トレンド仮定の検証 平行トレンドの仮定をデータから直接検証することはできない（潜在変数に依存しているため） 次善策: 処置以前のデータを複数期間取得し、事前の機関に平行トレンドが成立しているかを検証 3 期間データの場合: 処置以前の ~ 間のデータから をテストし、両グループで事前トレンドに差がないことを確認する 処置後の 期においても、同様の傾向が継続しているのではないかと考える t − 1 t E[Yi,t − Yi,t−1 ∣ Di = 1] = E[Yi,t − Yi,t−1 ∣ Di = 0] t + 1 35 / 48
36. 平行トレンド仮定の検証 : モデル . 処置は で実施される 処置効果と時点の固定効果は0 とする ここで、 : 実際のデータで処置郡に属しているかどうか : 潜在的な処置 このとき、 について以下が成り立つ。 t ∈ {0, 1, 2} t = 2 Yi,t (d) = ϕi + tDi + ϵi,t , i = 1, … , n, t = 0, 1, 2, d = 0, 1. Di d t = 0, 1 E[Yi,t+1 (0) − Yi,t (0) ∣ Di = 1] = 1 E[Yi,t+1 (0) − Yi,t (0) ∣ Di = 0] = 0 36 / 48
37. 処置群の結果変数は、処置を受ける以前から上昇傾向にあるため、平行トレンドの仮定が満たさ れない 一種のセレクションバイアス（人口が増加傾向にある地域に駅を新設するなど） データから事前トレンドの有無を検定するために が十分0 に近い 事前トレンドなし の推定に進む と は以下の回帰モデルの のOLS 推定量と一致する モデルの設定より、平行トレンドは満たされていない。 は有意に0 と異なるはず 事前トレンドテストを誤って通過することも（type-I error ） △ ˆ Trend = ( ^ Y 0,D=1 − ^ Y 1,D=1) − ( ^ Y 0,D=0 − ^ Y 1,D=0) △ ˆ Trend ⇔ ⇒ ATET2 △ ˆ Trend ˆ ATET2 β1 , β2 Yi,t = (δt,0 Di )β1 + (δt,2 Di )β2 + n ∑ j=1 Ii,j ϕj + ∑ k≠1 δt,k ψk + ϵi,t i = 1, … , n, t = 0, 1, 2. β1 ⇒ 37 / 48
38. 数値シミュレーションで確かめてみる n <- 80 # 標本サイズ ID <- rep(1:n, each = 3) # 個人ID t <- rep(c(0,1,2), times = n) # 時点 D <- rep(c(1,0), each = 3*(n/2)) # トリートメント phi <- rep(runif(n), each = 3) # 個人効果: [0,1] 上の一様分布 のp 値が10% 以上の時に のp 値が10% 以上か以下かを調べたい。 sim <- function(){ Y <- phi + D * t + rnorm(3 * n) t0D <- D * (t == 0) t2D <- D * (t == 2) df <- tibble(ID, Y, D, t, t0D, t2D) reg <- feols(Y ~ t0D + t2D | ID + t, data = df) pval <- reg$coeftable[, 4] return(c(pval[1] < 0.1, pval[2] < 0.1)) } △ ˆ Trend ˆ ATET 2 38 / 48
39. nrep <- 500 res <- matrix(0, nrep, 2) for(i in 1:nrep){ res[i, ] <- sim() } # 事前トレンド非有意、ATET 非有意 sum(res[, 1] == FALSE & res[, 2] == FALSE) ## [1] 0 # 事前トレンド非有意、ATET 有意 sum(res[, 1] == FALSE & res[, 2] == TRUE) ## [1] 40 事前トレンドのp 値が10% 以上である場合、全てのケースで が10% 水準で有意な結果 に。 真の は0 なので、これは良くない 事前トレンドテストを通過しても、推定値が誤っている可能性がある ATET2 ATET2 39 / 48
40. 時点 における処置群と対照群の結果変数の差は、 より となる。これの標本対応は . 従って、サンプルサイズが十分あれば が期待される。また、 より、 と表せる。どちらも の引き算が含まれるため、両者の値は を通じて相関している。 t Yi,t (d) = ϕi + tDi + ϵi,t , i = 1, … , n, t = 0, 1, 2, d = 0, 1. E[Yi,t ∣ Di = 1] − E[Yi,t ∣ Di = 0] = t △ ^ Y t = ^ Y t,D=1 − ^ Y t,D=0 △ ^ Y t ≃ t △ ˆ Trend = ( ^ Y 0,D=1 − ^ Y 1,D=1) − ( ^ Y 0,D=0 − ^ Y 1,D=0) ˆ ATET 2 = ( ^ Y 2,D=1 − ^ Y 1,D=1 ) − ( ^ Y 2,D=0 − ^ Y 1,D=0 ) △ ˆ Trend = △ ^ Y 0 − △ ^ Y 1 , ˆ ATET2 = △ ^ Y 2 − △ ^ Y 1 △ ^ Y 1 △ ^ Y 1 40 / 48
41. 平行トレンドの検証 : まとめ 平均的に、 の値が0 に近いほど は0 に近くなり、 は2 に近くなる。 「 が偶然小さな値に下振れしたデータ」を引き当てるまでモデルの設定やサンプリングを 変えることで、真の処置効果が0 であった場合でも、 を有意（事前トレンドの差が有意 ではないのに）にできる。 p-hacking となる。 △ ^ Y 1 △ ˆ Trend ˆ ATET 2 △ ^ Y 1 ATET2 ⇒ 41 / 48
42. Card and Krueger (1994) を R でやってみる Card and Krueger (1994) では、1992 年にNew Jersey で最低賃金が引き上げられた一方で、ペンシル ベニアではそうならなかったことを使って、最低賃金が雇用量に与える影響をDID で調べた store_id: 店舗のID post: 最低賃金引き上げ後ダミー emp: 従業員数 RJ: NJ 州ダミー price: ミールセットの価格 42 / 48
43. データ Unique (#) Missing (%) Mean SD Min Median Max store_id 357 0 249.7 153.9 1.0 238.0 999.0 post 2 0 0.5 0.5 0.0 0.5 1.0 emp 112 0 21.0 9.5 0.0 20.0 85.0 price 187 6 3.3 0.7 2.1 3.1 5.9 nj 2 0 0.8 0.4 0.0 1.0 1.0 43 / 48
44. 2 方向固定効果モデルによる DID mwdata %>% feols(emp ~ post:nj | store_id + post) |> modelsummary(stars = TRUE, coef_omit = c(-1) , gof_omit = "IC|Adj|F|RMSE|Log", output = "kableExtra") (1) post × nj 2.326 (1.452) Num.Obs. 714 R2 0.787 R2 Within 0.011 Std.Errors by: store_id + p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001 44 / 48
45. 手動 DID 回帰分析と推定値が等しいことを確認する。 mwdata |> group_by(nj) |> summarise(mean_post = mean(emp[post == 1]), mean_pre = mean(emp[post == 0])) |> mutate(diff = mean_post - mean_pre) |> summarise(did = diff[nj == 1] - diff[nj == 0]) ## # A tibble: 1 × 1 ## did ## <dbl> ## 1 2.33 45 / 48
46. ミールセットの価格の対数をアウトカムに mwdata |> filter(is.na(price)) |> print(n = 4) ## # A tibble: 42 × 5 ## store_id post emp price nj ## <dbl> <dbl> <dbl> <dbl> <dbl> ## 1 234 0 8 NA 1 ## 2 324 0 48 NA 1 ## 3 326 0 45.5 NA 1 ## 4 375 0 34 NA 1 ## # ℹ 38 more rows price に欠損値があるため、Unbalanced panel data になっている。欠損地があるものは除外してバ ランスさせる。 46 / 48
47. mwdata |> tidyr::drop_na(price) |> group_by(store_id) |> filter(n() == 2) ## # A tibble: 634 × 5 ## # Groups: store_id [317] ## store_id post emp price nj ## <dbl> <dbl> <dbl> <dbl> <dbl> ## 1 6 0 5 4.72 1 ## 2 14 0 16 4.4 1 ## 3 26 0 41.5 2.95 1 ## 4 27 0 13 4.25 1 ## 5 30 0 10 4.65 1 ## 6 31 0 15 4.65 1 ## 7 33 0 26 3.04 1 ## 8 36 0 18 4.5 1 ## 9 77 0 35 3.42 1 ## 10 78 0 21 3.16 1 ## # ℹ 624 more rows このデータにfeols() を適用して分析する。 47 / 48
48. (1) post × nj 0.033* (0.014) Num.Obs. 634 R2 0.925 R2 Within 0.016 Std.Errors by: store_id + p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001 最低賃金の引き上げは、統計的に有意にミールセットの価格を約3.25% 上昇させた Card and Krueger (1994) では、雇用主が最低賃金引き上げによるコスト増加の分を雇用量によっ て調整するのではなく、消費者に価格転嫁した可能性があると指摘 48 / 48