【パッケージ開発】fixesでイベントスタディを効率化

Staggered DiDにも対応しました。

R
公開

2025年5月25日

はじめに

本ページでは、two-way fixed effectsモデルのDiDにおけるイベントスタディを行うパッケージ、fixesについて説明します。

僕が今年の年始から地道に作成してきたパッケージで、一度tipsに書いたことがあるのですが、最近アップデートも重なってきたので再度書き直すことにしました。

ご意見待ってます

実戦経験が足りないパッケージなので、挙動がおかしかったり、こういう機能にも対応してほしいということがあれば、ページ下のコメントからどしどしお寄せください。こぢんまりしたページなので、いただいたものは余裕で全部対応できると思います。

これまでは処置個体において処置タイミングが同一のDiDにのみ対応していましたが、直近のアップデートで各個体で処置タイミングが異なる、いわゆるStaggered DiDにも対応しました。

果たしてあらゆるケースに対応できているか、個人開発のパッケージとしては不安ではありますが、いくつかのデータセットでテストした分には問題なかったので、そのテストコードと併せて解説していきます。

使用するパッケージ

使用するのはfixesです。インストールは以下のコマンドでできます。

pak::pak("fixes")

もしくは

install.packages("fixes")

です。開発版はGitHubからインストールできます。

pak::pak("yo5uke/fixes")

準備

パッケージを読み込んで、デモに使うデータを準備します。

library(fixes)
library(tidyverse)

データはfixestパッケージで用意されているものと、MixtapeのStaggered DiDの章で使用されているものを使います。

また、不足している変数を追加し、使わない変数もあるのであらかじめ除いておきます。

df_staggは未処置の個体の処置年、処置からの相対年をNAに変更し、df_castleに関してはMixtape内で行われている処理をしておきます1

df_base <- fixest::base_did
df_stagg <- fixest::base_stagg
castle <- haven::read_dta("https://github.com/scunning1975/mixtape/raw/master/castle.dta")


df_stagg <- df_stagg |> 
  mutate(
    year_treated = if_else(year_treated == 10000, NA_real_, year_treated), 
    time_to_treatment = if_else(time_to_treatment == -1000, NA_real_, time_to_treatment)
  )

# 以下Mixtape参考
dropped_vars <- c("r20004", "r20014",
                  "r20024", "r20034",
                  "r20044", "r20054",
                  "r20064", "r20074",
                  "r20084", "r20094",
                  "r20101", "r20102", "r20103",
                  "r20104", "trend_9", "trend_46",
                  "trend_49", "trend_50", "trend_51"
)

region <- castle |> 
  select(starts_with("r20")) |> 
  names() |> 
  setdiff(dropped_vars)

df_castle <- castle |> 
  select(state, sid, year, treatment_date, l_homicide, all_of(region), popwt) |> 
  mutate(is_treated = if_else(!is.na(treatment_date), 1, 0), .after = treatment_date)

各データフレームは以下のような感じになっています。

y x1 id period post treat
2.87530627 0.5365377 1 1 0 1
1.86065272 -3.0431894 1 2 0 1
0.09416524 5.5768439 1 3 0 1
3.78147485 -2.8300587 1 4 0 1
-2.55819959 -5.0443544 1 5 0 1
1.72873240 -0.6363849 1 6 1 1
id year year_treated time_to_treatment treated treatment_effect_true x1 y
90 1 2 -1 1 0 -1.0947021 0.01722971
89 1 3 -2 1 0 -3.7100676 -4.58084528
88 1 4 -3 1 0 2.5274402 2.73817174
87 1 5 -4 1 0 -0.7204263 -0.65103066
86 1 6 -5 1 0 -3.6711678 -5.33381664
85 1 7 -6 1 0 -0.3152137 0.49562631
state sid year treatment_date is_treated l_homicide r20001 r20002 r20003 r20011 r20012 r20013 r20021 r20022 r20023 r20031 r20032 r20033 r20041 r20042 r20043 r20051 r20052 r20053 r20061 r20062 r20063 r20071 r20072 r20073 r20081 r20082 r20083 r20091 r20092 r20093 popwt
Alabama 1 2000 2006 1 2.027356 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 4499293
Alabama 1 2001 2006 1 2.164867 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 4499293
Alabama 1 2002 2006 1 1.936334 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 4499293
Alabama 1 2003 2006 1 1.919567 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 4499293
Alabama 1 2004 2006 1 1.749841 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 4499293
Alabama 1 2005 2006 1 2.130440 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 4499293

パッケージの説明

基本的な使い方と一緒にパッケージについて説明していきます。関数はイベントスタディを実行するrun_es()関数とイベントスタディをggplotベースでプロットするplot_es()です。

run_es()

ひとまず、run_es()の基本的な引数を説明します。

run_es(
  data,
  outcome,
  treatment,
  time,
  staggered = FALSE,
  timing,
  lead_range = NULL,
  lag_range = NULL,
  covariates = NULL,
  fe,
  cluster = NULL,
  weights = NULL,
  baseline = -1,
  interval = 1,
  time_transform = FALSE,
  unit = NULL
)
  • data：使用するデータフレーム
  • outcome：アウトカム変数
  • treatment：処置されていれば（時間によらず）1をとるダミー変数2
  • time：時間を示す変数（yearなど）
  • staggered：処置タイミングが個体ごとに異なるモデルの場合はTRUE（デフォルトはFALSE）
  • timing：処置タイミングを表す変数（staggeredがFALSEの場合は数値、TRUEの場合は処置タイミングを示す変数）
  • lead_range：処置前の期間（処置年は含まず、デフォルトでは最大の期間をとる）
  • lag_range：処置後の期間（処置年は含まず、デフォルトでは最大の期間をとる）
  • covariates：共変量（~ cov1 + cov2のように書く）
  • fe：固定効果（~ id + yearのように書く）
  • cluster：クラスタリングの単位（~ idのように書く）
  • weights：ウェイト（~ weightのように書く）
  • baseline：基準となる相対年（デフォルトでは処置の1期前）
  • interval：時間のインターバル（デフォルトでは1（1年ごとのデータ）、5にすれば国勢調査に対応（5年ごと））
  • time_transform：TRUEで時間変数に通し番号を振って処理。時間変数が文字列だったりイレギュラーな間隔の場合に対応（デフォルトはFALSE）
  • unit：time_transformがTRUEの場合に個体の単位を指定（idなど）

文字で見てもわかりにくいと思うので、読み込んだデータを使って関数を回してみましょう。

通常のDiDにおけるイベントスタディ

es <- run_es(
  data = df_base, 
  outcome = y, 
  treatment = treat, 
  time = period, 
  timing = 6, 
  fe = ~ id + period, 
  cluster = ~ id
)

基本的なケースであればこれくらいの指定で済みます。プロットは後のセクションでまとめて載せます。

Staggered DiD 1

es_stagg1 <- run_es(
  data = df_stagg, 
  outcome = y, 
  treatment = treated, 
  time = year, 
  staggered = TRUE, 
  timing = year_treated, 
  fe = ~ id + year, 
  cluster = ~ id
)

先ほどと違う点は、staggeredがTRUEになっており、timingが具体的な数値ではなく、処置年を示す変数を指定しているところです。

Staggered DiD 2

# 共変量の変数を式の形にしておく
covariates <- as.formula(paste("~", paste(region, collapse = "+")))

es_stagg2 <- run_es(
  data = df_castle, 
  outcome = l_homicide, 
  treatment = is_treated, 
  time = year, 
  staggered = TRUE, 
  timing = treatment_date, 
  covariates = covariates, 
  fe = ~ state + year, 
  cluster = ~ sid, 
  weights = ~ popwt
)

今回の例は共変量が多いので事前にas.formula()を使って処理してしまっていますが、関数内でcovariates = ~ cov1 + cov2のような指定の仕方で問題ありません。

plot_es()

こだわりがなければ、plot_es(es)の程度の書き方でプロットを出力することができます。引数の説明は後に回して、一旦先ほどのイベントスタディの結果をプロットしてみましょう。

通常のDiDにおけるイベントスタディ

plot_es(es)

後の2つもそうですが、ベースラインがデフォルトのままなので、処置前年（-1）が抜けて値が0をとっています。

Staggered DiD 1

plot_es(es_stagg1)

Staggered DiD 2

plot_es(es_stagg2)

Staggeredの方もいい感じにプロットできています。この図はMixtape内でlead変数とlag変数を手作業で作って実装しプロットされていますが、その図と見比べても同じ推移が示されていることがわかります。

使い方

基本的には上の使い方で十分ですが、ggplotベースで作ったこともあり、柔軟な対応も可能です。

plot_es(
  data,
  type = "ribbon",
  vline_val = 0,
  vline_color = "#000",
  hline_val = 0,
  hline_color = "#000",
  linewidth = 1,
  pointsize = 2,
  alpha = 0.2,
  barwidth = 0.2,
  color = "#B25D91FF",
  fill = "#B25D91FF"
)
  • data：run_es()で作成した結果のデータフレーム
  • type：ribbonとerrorbarが指定可能
  • vline_val：縦の破線の位置（デフォルトは0）
  • vline_color：縦の破線の色（デフォルトは#000（黒））
  • hline_val：横の破線の位置（デフォルトは0）
  • hline_color：縦の破線の色（デフォルトは#000（黒））
  • linewidth：折れ線の太さ（デフォルトは1）
  • pointsize：点の大きさ（デフォルトは2）
  • alpha：リボンの透明度（デフォルトは0.2）
  • color：線と点の色（デフォルトは#B25D91FF（ピンク））
  • fillリボンの色（デフォルトは#B25D91FF（ピンク））

plot_es()関数内でもある程度自由度をもって設定できます。
例えばプロットをエラーバーのタイプにしてみましょう。

plot_es(es_stagg2, type = "errorbar")

さらに色も変えてみます。

plot_es(es_stagg2, type = "errorbar", color = "#000")

リボンタイプとエラーバータイプがあり、さらにそこからある程度の自由度もある点が他のパッケージに対する優位性にもなっています。

ちなみに、ggplotベースということで、ggplotの関数をつなげることも可能です！

plot_es(es_stagg2, type = "errorbar", color = "#000") + 
  geom_hline(yintercept = .08, color = "red") + 
  annotate("text", x = 3, y = -.1, label = "DD Coefficient = 0.08\n(s.e. = 0.03)")

Mixtape風に係数を示す水平線と注釈を追加してみました。このようにどんどん要素を追加していくことができます。

おわりに

まだまだ開発途上ですが、staggeredまで対応できたので、これからしばらくはブラッシュアップを念頭にアップデートしていこうと思います。

冒頭にも書きましたが、もし使っていただいて「使いにくい」や「こんな機能が欲しい」等ありましたら、以下のコメントやGitHubのIssuesまでお願いします。

🔗 データ出典

データ出典：Causal Inference: The Mixtape（Scott Cunnningham, 2021）
データはGitHub上で公開されており、以下のリポジトリから入手可能です：

https://github.com/scunning1975/mixtape

また、データはMITライセンスのもとで提供されています。

注

  1. 参考：Mixtape↩︎

  2. TRUE/FALSEでも可↩︎