はじめに

こちらの記事でSurrogate indexについて初めて知りました。

developers.cyberagent.co.jp

Netflixでも活用されていて、一定の成果を上げているようです。

arxiv.org

業務で扱う課題を解決してくれる可能性があったため、理解しておきたいと思いました。

論文等を読みつつRで挙動などを確認していきます。

はじめに
概要
- 課題感
- 方針
- 仮定
  - 仮定1 強く無視できる割り当て条件
  - 仮定2 代理性
  - 仮定3 互換性
定義
- Surrogate index：代理指数
- 代理スコア：Surrogate Score
- Sampling Score
- 潜在的条件付き期待値
関係性
因果効果の推定
論文での適用事例
Rで確認
- データの生成
- 推定
- 所感

はじめに

こちらでStaggered DIDを初めて知りました。

blog.brainpad.co.jp

介入タイミングが異なる場合のDIDということで、実務でも利用シーンがありそうな問題設計だと思いました。

簡単にまとめたいと思います。

概論

元の論文は以下になります。

arxiv.org

また、こちらの資料も参考になった。

DIDの計量経済手法の近年の展開～問題と対策～, 小西祥文慶應義塾大学

わかりやすい説明はこちらに譲ります。

blog.brainpad.co.jp

通常の差分の差分法(DID)では、介入タイミングが単一であり介入直後のLATEを推定することはできる。

介入タイミングが複数あり、動的な効果を知りたいが知りたい場合の拡張が議論されている。

Two Way Fixed Effect Event Study (TWFE Event Study)

$Y_{i,t}$ のモデルを以下のように定義する。

$\displaystyle{ Y_{i,t} = \alpha_i +\lambda_t + \sum_l \mu_l \mathbb{I} \{t − E_i = l \} +v_{i,t} }$

ここで $\alpha_i$ がグループiの固定効果、 $\lambda_t$ が時間固定効果、 $\mu_l$ は介入タイミングから時刻lだけずれた時点における介入効果。

$\mathbb{I} \{ t − E_i = l \}$ は介入タイミングからのずれを表現したインジケータである。 $v_{i,t}$ は誤差項である。

TWFEは介入タイミングによって介入効果が異なる場合/異質性がある場合は、大きなバイアスが発生してしまうことがある。

Staggered DiD

介入タイミングのグループ項を入れることで、介入タイミングが異なるコーホート同士の比較を適切に行えるようになる。

$\displaystyle{ Y_{i,t} = \alpha_i +\lambda_t + \sum_{c \in C} \sum_l \delta_{c,l} \mathbb{I} \{C_i = c \} \mathbb{I} \{t − E_i = l \} +\epsilon_{i,t} }$

ここで $\delta_{c,l}$ がコホートcと介入タイミングlの介入効果、 $\mathbb{I} \{C_i = c \}$ はiがコホート $c$ に所属しているかどうかのインジケータである。

$\epsilon_{i,t}$ は誤差項である。

Dynamic treatment effects $ATT_l$ は次のように推定される。

$\displaystyle{ ATT_l = \sum_c w_{l,c} \hat{\delta}_{l,c} }$

平均処置効果は次のように推定される。

$\displaystyle{ \theta = \sum_c w_{c} \hat{\delta}_{l,c} }$

ここで $w_{l,c}$ はコホートcと介入タイミングlに該当するサンプルの割合(重み)である。

$\displaystyle{ w_{l,c} = \frac{E [ \tilde{D}_{it} \mathbb{I} (t,k) ] } {E [ \tilde{D}^2_{it} ] } }$

$\tilde{D}_{it}$ は残差であり次にように分解できる。

$\displaystyle{ \tilde{D}_{it} = D_{it} - \hat{D}_{it} = D_{it} - \bar{D}_{i} - \bar{D}_{t} + \bar{\bar{D}}_{i} }$

期間中の処置確率が高いグループほど $D_i$ の値が高くなり，処置確率の高い期間ほど $D_T$ の値が高くなる。

これにより $\tilde{D}_{it}$ がマイナスの値をとる場合が発生し、処置効果の推定に大きなバイアスを与えてしまう。

Rでの実行

RでStaggered DiDを実行するには、fixest packageを用いる。

いつくかのサンプルデータをモデルに当てはめて、モデルの特定を把握しておきます。

処置効果に異質性がある場合

サンプルデータとして以下のようなデータを作成した。

$N(100,10)$ に従う時系列を200件のサンプルを生成して、4つのグループに分けています。

そして、1~3のグループにおいてT=15, 20, 25で+1.5t, +2.5t, +3.5tとなるような介入効果を設定しています。

group=4が対照群としている。

did_df %>% 
  group_by(group, t) %>% 
  summarise(y = mean(y)) %>% 
  ggplot()+
  aes(x=t, y=y, col= group)+
  geom_line()+
  ylim(80,180)

TWFEモデルを当てはめる場合は以下のようになる。

i()が交互作用項であり、ref=が多重共線性を避けるための基準とするタイミングtを設定するためのパラメータ。

res_twfe = feols(y ~ i(time_to_treatment, ref = c(-1)) | t+id, data = did_df, cluster = "id")

時刻ごとの介入効果を確認。 TWFEモデルに当てはめた場合は、観測初期に介入影響があると算出されてしまう。

iplot(res_twfe)

Staggered DiDモデルを当てはめる場合は以下のようになる。

sunab()がStaggered DiDの効果項を定義するものであり、コホートと介入タイミングを設定できる。 ref.cはダミー変数を作成しないレベルを指定する。(多重共線性を避けるため)

res_stagg = feols(y ~ sunab(treat, t, ref.c = max(treat)) | t+id, data = did_df, cluster = "id")

iplot(res_stagg)

Staggered DiDを当てはめた場合は、観測初期の介入影響はない。

ATTは以下で算出することができる。ここでのATTとは時間方向と介入効果における加重平均を取ったものとなる。

> summary(res_stagg, agg = "att")
OLS estimation, Dep. Var.: y
Observations: 8,200 
Fixed-effects: t: 41,  id: 200
Standard-errors: Clustered (id) 
    Estimate Std. Error t value  Pr(>|t|)    
ATT  31.5197    1.14343  27.566 < 2.2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
RMSE: 9.72441     Adj. R2: 0.822263
                Within R2: 0.594405