【R】差の差法でイベントスタディやるときのコード
Rのfixestパッケージが先月末にアップデートされたようです。
このブログ投稿はこのバージョンの話ではないですが、3つ以上の固定効果を使うと、固定効果の推定がおかしくなるバグを直したらしいのでアプデ推奨のようです。
⚠ Hot Fix ⚠
— Laurent Bergé (@lrberge) April 1, 2022
Concerns #Rstats {fixest} users working in int. trade
When using 3+ fixed-effects, the value of the FE coefficients extracted with {fixef()} could be erroneous.
Pls update to 0.10.4 to fix this (ps: you need to build from source)https://t.co/mmEzLLlVEg
1/2
さて、本題です。
差の差法（Difference in Differences, DID）のイベントスタディ*1をやろうと思って推定しても、図書くのめんどいなということがあったのですが、fixestの関数使えば簡単にできるようになっていたのでメモ。*2
準備
fixestパッケージに入っているbase_didというデータセットで実演します。
# パッケージロード library(fixest) # データサマリー > skimr::skim(base_did) ── Data Summary ──────────────────────── Values Name base_did Number of rows 1080 Number of columns 6 _______________________ Column type frequency: numeric 6 ________________________ Group variables None ── Variable type: numeric ───────────────────────────────────────────────────────────────────────────────────────────────────────────── skim_variable n_missing complete_rate mean sd p0 p25 p50 p75 p100 hist 1 y 0 1 2.02 5.71 -15.1 -1.87 1.74 5.78 21.5 ▁▆▇▃▁ 2 x1 0 1 0.0353 2.98 -10.6 -2.04 0.0941 2.08 9.75 ▁▃▇▅▁ 3 id 0 1 54.5 31.2 1 27.8 54.5 81.2 108 ▇▇▇▇▇ 4 period 0 1 5.5 2.87 1 3 5.5 8 10 ▇▇▇▇▇ 5 post 0 1 0.5 0.500 0 0 0.5 1 1 ▇▁▁▁▇ 6 treat 0 1 0.509 0.500 0 0 1 1 1 ▇▁▁▁▇ # timing of policy: 政策の施行はピリオド6から table(base_did$period, base_did$post) # treated individuals: 1-55 政策が施行されるグループには1番から55番の個人が入る。全部で108番までいる。 table(base_did$id, base_did$treat)
2 x 2の差の差法
一番基本的な差の差法の推定である、施行前後ダミー・処置ダミー・政策ダミー（2つのダミーの相互作用）による2x2の推定は以下の通り
# 2 x 2 model1 <- feols(y ~ post*treat, data = base_did)
> model1 OLS estimation, Dep. Var.: y Observations: 1,080 Standard-errors: IID Estimate Std. Error t value Pr(>|t|) (Intercept) 0.322801 0.317345 1.017194 3.0929e-01 post 0.713468 0.448793 1.589747 1.1219e-01 treat 0.142340 0.444695 0.320084 7.4897e-01 post:treat 4.993390 0.628893 7.939967 5.0601e-15 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 RMSE: 5.15642 Adj. R2: 0.180566
固定効果モデル
2-way fixed effect modelは二元配置固定効果推定法というのかな？
時間と個人の固定効果をそれぞれ含めたモデルは以下の通り。
# 2 way fixed effects model2 <- feols(y ~ post:treat | id + period, data = base_did) # feols(y ~ post:treat, data = base_did, fixef = c("id","period")) でもOK
結果 > model2 OLS estimation, Dep. Var.: y Observations: 1,080 Fixed-effects: id: 108, period: 10 Standard-errors: Clustered (id) Estimate Std. Error t value Pr(>|t|) post:treat 4.99339 0.521434 9.57627 4.6502e-16 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 RMSE: 4.74941 Adj. R2: 0.222439 Within R2: 0.064601
イベントスタディ
事前トレンドの仮定を評価するために、イベントスタディモデルで推定したい場合は処置ダミーと時間ピリオドの相互作用項を使います。
i()関数で簡単にできる上に、ベースとなる時間ピリオドの指定もできます。
この場合は、時間ピリオドの6以降が政策施行の期間なので、直前の5をベースにしています。
# event study model3 <- feols(y ~ i(period,treat, ref = 5) | id+period, data = base_did)
> model3 OLS estimation, Dep. Var.: y Observations: 1,080 Fixed-effects: id: 108, period: 10 Standard-errors: Clustered (id) Estimate Std. Error t value Pr(>|t|) period::1:treat -2.015417 1.34235 -1.501414 1.3619e-01 period::2:treat -1.664474 1.38858 -1.198689 2.3330e-01 period::3:treat 0.504127 1.32338 0.380939 7.0400e-01 period::4:treat -0.884617 1.41588 -0.624784 5.3344e-01 period::6:treat 1.159133 1.22689 0.944775 3.4690e-01 period::7:treat 4.334751 1.30818 3.313584 1.2575e-03 ** period::8:treat 3.825722 1.51448 2.526102 1.2998e-02 * period::9:treat 4.640141 1.29328 3.587874 5.0439e-04 *** period::10:treat 6.946824 1.42559 4.872933 3.8321e-06 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 RMSE: 4.69194 Adj. R2: 0.234782 Within R2: 0.087104
いい感じの結果が出ています（作られたデータなんで当たり前ですが）。
この結果を図にするときに、わざわざ係数と標準誤差を抜き出して、データフレームにして加工してからggplotで描画、なんてしていましたが、
少なくともざっと結果を見る分にはその必要がありません！
係数をプロットする関数を用意してくれています。
係数プロット
coefplot(model3)
coefplot()関数を使うと、推定された係数を標準誤差による95%信頼区間のエラーバー付きで表示してくれますが、もし他に説明変数を入れていると、それも表示してしまいます。引数dropなどで編集はできますが、次のiplot()関数を使うと、イベントスタディとしてすごく見やすい図を作ってくれます。
iplot(model3)
ベースとなる時間ピリオドに縦線を勝手に引いてくれるので、政策前後がわかりやすくなっています。i()関数で推定した係数のみを表示するので、変数の選択などの引数も不要です。
fixestまじですごいな…以上メモでした。
*1:event study 9 Difference-in-Differences | Causal Inference
*2:この関数自体は昨年6月頃の0.9.0のアップデートで実装されてたようです。
