---
title: "差分の差 (DiD) デザイン:実務家のためのガイド"
subtitle: "Baker, Callaway, Cunningham, Goodman-Bacon, Sant'Anna (2025) を R で動かしながら学ぶ"
author: "学習用チュートリアル"
date: today
date-format: "YYYY-MM-DD"
lang: ja
format:
html:
toc: true
toc-depth: 3
toc-location: left
toc-title: "目次"
number-sections: true
number-depth: 3
theme: cosmo
highlight-style: github
code-fold: true
code-tools: true
code-link: true
code-overflow: wrap
code-line-numbers: false
embed-resources: true
fig-width: 8
fig-height: 5
fig-align: center
fig-dpi: 120
tbl-cap-location: top
fig-cap-location: bottom
callout-appearance: default
callout-icon: true
df-print: kable
smooth-scroll: true
link-external-newwindow: true
execute:
eval: true
echo: true
warning: false
message: false
cache: true
fig-format: png
knitr:
opts_chunk:
comment: "#>"
collapse: false
dev: "png"
dpi: 120
editor: source
---
# はじめに(本ノートについて) {.unnumbered}
このノートは、Baker, Callaway, Cunningham, Goodman-Bacon, Sant'Anna による
**"Difference-in-Differences Designs: A Practitioner's Guide"** (arXiv:2503.13323v3, 2025) を
初学者がコードを動かしながら追体験するための教材です。論文の各節を順に取り上げ、
専門用語の解説と R コードを併走させます。
::: {.callout-important}
## 本ノートで使うデータについて
論文の実証例は **アメリカの郡レベル成人死亡率(CDC Vital Statistics 2009–2019) ×
ACA メディケイド拡大タイミング** を用いています。著者は
[openICPSR project 239070](https://www.openicpsr.org/openicpsr/project/239070/version/V1/view)
(DOI: 10.3886/E239070V1) で公式の再現パッケージを公開していますが、
ダウンロードには ICPSR アカウントが必要です。
本ノートでは **論文の Table 2(2×2 DiD)と Table 4(共変量バランス)の数値を再現するように
較正したシミュレーションデータ** を QMD 内で生成します。
真の DGP(データ生成過程)が分かっているので、推定量の振る舞いを理解しやすいという
教育上の利点もあります。
:::
::: {.callout-tip}
## 本ノートの読み方
- **コード**: 各チャンクは折り畳めます(右上の `</>` から `Show All Code` で全展開)
- **専門用語**: 色付き callout box で初出時に説明します
- **章番号**: 論文と同じ番号(§3 = 「2×2 DiD デザイン」など)を踏襲しています
:::
## 環境準備:パッケージのインストールとロード
```{r}
#| label: setup-packages
#| code-fold: true
#| code-summary: "パッケージのインストール・ロード(クリックで展開)"
# 利用するパッケージ一覧
pkgs <- c(
"tidyverse", # データ操作・可視化の基盤
"fixest", # 高速な固定効果推定(TWFE)
"did", # Callaway-Sant'Anna 推定量(段階的 DiD の中核)
"DRDID", # 2×2 DiD の二重頑健推定
"bacondecomp", # Goodman-Bacon 分解(TWFE の問題を可視化)
"didimputation", # Borusyak et al. 補完推定量
"staggered", # Roth & Sant'Anna の段階的 DiD 推定量
"modelsummary", # 回帰結果の整形表
"gt", # 整形テーブル
"kableExtra", # 整形テーブル(代替)
"broom", # 回帰結果を tidy データに
"patchwork", # ggplot 並列
"scales", # 軸ラベル整形
"knitr" # チャンク制御
)
# 未インストールのものを自動インストール
for (p in pkgs) {
if (!requireNamespace(p, quietly = TRUE)) {
install.packages(p, repos = "https://cloud.r-project.org")
}
}
# 全てを読み込み
invisible(lapply(pkgs, library, character.only = TRUE))
# ggplot のテーマ・日本語フォント(Windows 環境を想定)
theme_set(theme_minimal(base_size = 12))
if (.Platform$OS.type == "windows") {
windowsFonts(JP = windowsFont("MS Gothic"))
theme_update(text = element_text(family = "JP"))
}
# 表示桁数
options(scipen = 999, digits = 4)
```
::: {.callout-note collapse="true"}
## パッケージの役割の要約(クリックで展開)
| パッケージ | 用途 |
|---|---|
| `tidyverse` | データ加工・可視化(`dplyr`, `ggplot2`, `tidyr` などのメタパッケージ) |
| `fixest` | 高速 OLS + 固定効果。`feols()` が主役。`i()` 関数でイベント・スタディが簡潔 |
| `did` | Callaway-Sant'Anna 推定量 `att_gt()`+集約 `aggte()`。段階的 DiD の標準実装 |
| `DRDID` | 2×2 DiD のアウトカム回帰(Outcome Regression)/ 逆確率ウェイト(Inverse Probability Weighting)/ 二重頑健(Doubly Robust)推定量 |
| `bacondecomp` | TWFE 推定量を全 2×2 比較の加重平均に分解(Goodman-Bacon 2021) |
| `didimputation` | Borusyak-Jaravel-Spiess (2024) の補完推定量 |
| `staggered` | Roth & Sant'Anna (2023) 効率的推定量 |
| `modelsummary` | `feols`/`lm` 結果を見やすい表に |
:::
---
# はじめに(論文 §1)
差分の差(DiD)の起源は **1840 年代の John Snow によるロンドンのコレラ研究**
まで遡る、社会科学で最も古い因果推論手法のひとつです。基本構造は次のとおり:
1. **2 つのグループ**:介入を受ける群(treated)と受けない群(control)
2. **2 つの時期**:介入前(pre)と介入後(post)
::: {.callout-note}
## 用語:差分の差 (Difference-in-Differences, DiD)
「**介入群の前後変化**」から「**対照群の前後変化**」を引いた量(=「2 つの差の差」)を
介入効果と見なす方法。式で書けば:
$$
\widehat{ATT} = \underbrace{(\bar{Y}_{T,\text{post}} - \bar{Y}_{T,\text{pre}})}_{\text{介入群の変化}}
- \underbrace{(\bar{Y}_{C,\text{post}} - \bar{Y}_{C,\text{pre}})}_{\text{対照群の変化}}
$$
「対照群の変化」を「介入群が介入を受けなかった場合の反事実」とみなすことで、
共通の時間トレンドをキャンセルできます。
:::
::: {.callout-note}
## 用語:ATT (Average Treatment effect on the Treated)
**実際に介入を受けた人たちにとっての**平均介入効果。
DiD で標準的にターゲットにする因果パラメータです。
潜在アウトカム記法では:
$$
ATT(t) = E[Y_{i,t}(1) - Y_{i,t}(0) \mid D_i = 1]
$$
ここで:
- $D_i$ は **介入群ダミー**:単位 $i$ が(いずれ)介入を受ける群なら $D_i = 1$、
受けない群なら $D_i = 0$ をとる二値変数
- $Y(1)$ は介入を受けた場合のアウトカム、$Y(0)$ は受けなかった場合の(反事実)アウトカム
- $E[\cdot \mid D_i = 1]$ は「**介入を受けた群に限った**(条件付き)期待値」
:::
::: {.callout-note}
## 用語:平行トレンド(Parallel Trends)
「もし介入がなかったとしたら、介入群と対照群のアウトカムは同じトレンドで変化していた」
という反事実的な仮定。DiD の妥当性の根幹を成す **検定不可能な** 仮定です。
$$
E[Y_{i,t}(0) - Y_{i,t-1}(0) \mid D_i = 1] = E[Y_{i,t}(0) - Y_{i,t-1}(0) \mid D_i = 0]
$$
:::
## なぜ今 DiD を学び直すのか
近年の方法論研究(de Chaisemartin & D'Haultfœuille 2020, Goodman-Bacon 2021,
Sun & Abraham 2021, Borusyak et al. 2024 等)は、**介入開始時期が単位ごとに異なる
"段階的採択"(staggered adoption)の状況で、伝統的な二方向固定効果(Two-Way Fixed Effects, TWFE)回帰は誤った推定値を返しうる**
ことを示しました。この問題に対処するため、論文は次の枠組みを提案しています:
> **"フォワード・エンジニアリング"** = まず目標パラメータを定義し、それから推定法を導出する。
> 馴染みのある回帰式から出発するのではなく、**因果量 → 識別仮定 → 推定量** の順で組み立てる。
## 2×2 DiD の概念図
```{r}
#| label: fig-concept
#| fig-cap: "2×2 DiD の概念図:4 つの平均と DiD 推定量"
#| fig-height: 4.2
concept <- tibble::tribble(
~group, ~time, ~Y,
"介入群", "介入前", 100,
"介入群", "介入後", 130, # 介入群:傾き +30(対照群より急)
"対照群", "介入前", 70,
"対照群", "介入後", 80 # 対照群:傾き +10
)
concept$time <- factor(concept$time, levels = c("介入前", "介入後"))
# 反事実(介入群が介入を受けなかった場合):対照群と同じ傾きで介入群の点から延長
cf <- tibble::tibble(
time = factor(c("介入前", "介入後"), levels = c("介入前", "介入後")),
Y = c(100, 100 + (80 - 70)) # = c(100, 110)
)
ggplot(concept, aes(time, Y, color = group, group = group)) +
geom_line(linewidth = 1.2) +
geom_point(size = 3.2) +
# 反事実(介入群が介入を受けなかった場合):対照群と同じ傾きで延長した「灰色の平行線」
# 注: x 軸が factor かつ inherit.aes = FALSE のときは group = 1 を明示しないと
# geom_line が描画されない
geom_line(
data = cf, aes(time, Y, group = 1), inherit.aes = FALSE,
linetype = "dashed", color = "grey40", linewidth = 1.0
) +
geom_point(
data = cf, aes(time, Y), inherit.aes = FALSE,
shape = 21, fill = "white", color = "grey30",
size = 3.2, stroke = 1.1
) +
annotate("text",
x = 2.05, y = 110, hjust = 0,
label = "反事実\n(対照群と平行)",
color = "grey25", size = 3.6, lineheight = 0.9
) +
# ATT 両矢印
annotate("segment",
x = 2.02, xend = 2.02,
y = 130, yend = 110,
arrow = arrow(length = unit(0.2, "cm"), ends = "both"),
color = "firebrick", linewidth = 1.0
) +
annotate("text",
x = 2.07, y = 120, hjust = 0,
label = "ATT", color = "firebrick",
fontface = "bold", size = 5
) +
scale_color_manual(values = c("介入群" = "firebrick", "対照群" = "grey40")) +
expand_limits(x = 2.7, y = c(60, 140)) +
labs(
x = NULL, y = "アウトカム Y",
color = NULL,
title = "2×2 DiD:対照群の傾きで「反事実」を作り、ATT を推定する"
) +
theme(legend.position = "bottom")
```
赤い両矢印が DiD で推定したい量(ATT)です。
**対照群の前後変化(灰色実線)を「介入群が介入を受けなかった場合のトレンド」とみなし**、
それを介入群の介入前の値から平行に延長して(**灰色の破線**)反事実を構築します。
灰色の破線は対照群の実線と同じ傾きで描かれ、平行であることが見て分かるはずです。
実際の介入群の介入後値(赤実線の終点)と反事実(灰色の白丸)との差が ATT です。
---
# メディケイドと死亡率 (論文 §2)
論文は次の因果的問いを実証例とします:
> **ACA(医療費負担適正化法)の下でのメディケイド拡大は、成人の死亡率にどのような影響を
> 与えたか?**
メディケイド拡大は **段階的採択** の好例です。元々は 2014 年に全州が拡大予定でしたが、
2012 年の連邦最高裁判決により拡大は **任意** となり、州ごとに採択タイミングがバラバラになりました。
::: {.callout-note}
## 用語:段階的採択(staggered adoption)
介入の開始時期が単位(州、企業、個人など)ごとに異なる状況。古典的な 2×2 DiD は
「全員が同時に介入される」ことを想定していますが、現実の政策評価では多くがこの段階的設定。
:::
::: {.callout-note}
## 用語:コホート / 介入群(cohort, group, $G_i$)
各単位が **最初に介入を受けた年** を表す変数。$G_i = 2014$ なら 2014 年拡大群、
$G_i = \infty$ なら未拡大群(never-treated)を意味します。
:::
::: {.callout-note}
## 用語:潜在アウトカム(Potential Outcomes)
「もし介入を受けたとしたら」の値 $Y_i(1)$ と「もし受けなかったとしたら」の値 $Y_i(0)$ を
区別して考える枠組み(Rubin 1974)。実際に観察できるのは片方だけなので、もう片方は反事実です。
:::
::: {.callout-note}
## 用語:先取り行動なし(No Anticipation, NA 仮定)
「介入開始 **前** には介入効果はゼロ」という仮定。
将来の政策発表で人々が事前に行動を変えてしまうと、この仮定は破れます。
:::
## シミュレーションデータの生成
::: {.callout-tip}
## このセクションは読み飛ばして OK
ここから「### 生成したデータの概要」までは、**DiD 手法そのものではなく
データ生成の裏側** を説明するパートです。DiD の本筋を追いたい読者は
[§3 2×2 DiD デザイン](#did-デザイン論文-3) まで飛ばして問題ありません。
ただし、**無加重 vs 人口加重** の callout(下記)だけは DiD 全体の解釈に関わるので、
そこは目を通すことをお勧めします。
:::
::: {.callout-note}
## 用語:無加重(unweighted)と人口加重(population-weighted)
集計データ(郡レベル、州レベルなど)で平均をとるとき、**各単位をどう数えるか** に
2 つの選択肢があります。
**(a) 無加重平均**
$$
\bar{Y}_{\text{無加重}} = \frac{1}{n}\sum_{i=1}^{n} Y_i
$$
「人口 100 万の大都市」も「人口 1 千人の小郡」も **同じ 1 票** として扱う。
*郡そのもの* を分析対象とみなすときに自然です。
**(b) 人口加重平均**
$$
\bar{Y}_{\text{加重}} = \frac{\sum_{i=1}^{n} w_i Y_i}{\sum_{i=1}^{n} w_i}, \quad w_i = \text{郡 } i \text{ の人口}
$$
大都市は小郡より重く数える。**郡ではなく "人" を分析対象** とみなし、
「平均的な **アメリカ人** にとっての効果」を知りたいときに自然です。
### なぜ加重するのか — 3 つの理由
1. **解釈対象の単位を「人」に揃えるため**
集計データを使っていても、本当に知りたいのは「政策が人々の死亡率にどう影響したか」。
人口加重すれば、「平均的な郡」ではなく「平均的な個人」に対する効果が得られます。
2. **異質性のシグナルとして**
無加重と加重で結果が大きく違うなら、**効果の大きさが郡の人口と相関している** ということ。
本ノートの例では、無加重 ≈ 0 ・人口加重 ≈ −2.6 となり、
「**大きな郡ほど死亡率の低下が大きい**(=メディケイドが届きやすい)」ことを示唆します。
3. **精度のため(統計効率)**
人口の少ない郡は死亡率の分母が小さくノイズが大きい。
人口加重は、その「ノイズの大きい観測」の影響を抑える効果もあります。
### どちらが「正解」か?
**どちらも正解です。** 知りたい問いが「郡レベルの政策評価」なら無加重、
「人口あたり影響」なら人口加重。論文(と本ノート)は **両方を併記** し、
読者が解釈を選べるようにしています。
:::
論文の Table 1・Table 2・Table 4 の構造を踏襲してデータを作ります。
**ここで重要なのは、無加重 DiD ≈ 0、人口加重 DiD ≈ −2.6 となるように
真の介入効果を人口と相関させること**(論文 Table 2 を再現するため)。
これにより「無加重と人口加重で結論が変わる」現象を意図的に再現できます。
#### コードで何をしているか(DGP の組み立て)
シミュレーションコードは 5 ステップで構成:
| ステップ | コードの該当箇所 | 何をしているか |
|---|---|---|
| 1. コホート構成 | `cohort_sizes <- c(...)` | 全 2,604 郡を 5 コホート(2014/2015/2016/2019/Inf)に分配。論文 Table 1 と整合 |
| 2. 郡レベル不変属性 | `log_pop <- mapply(...)`, `covs_2013 <- map2_dfr(...)` | 2013 年人口(コホート別の対数正規)と共変量(% White、貧困率、所得等)を生成。コホート間で平均をずらし Table 4 の不均衡を再現 |
| 3. パネル展開 | `panel <- expand_grid(county_id, year) \|> left_join(...)` | 2,604 郡 × 11 年 = 28,644 行にバランス展開 |
| 4. アウトカム DGP | `alpha + gamma_t + true_att + epsilon` | $Y_{i,t}$ = 郡固有レベル + 年トレンド + **真の介入効果** + ノイズ。**真の ATT はコホート内中心化された人口偏差に比例**(無加重平均 0、加重平均 −2.6 になる仕掛け) |
| 5. 較正確認 | `calib <- ... \|> group_by(D, year) \|> summarise(...)` | 生成データで実際に無加重 ≈ −0.2、加重 ≈ −1.9 が出るか確認 |
ポイント:
- `set.seed(2025)` で再現性を確保
- $\tau_{i,t} = -0.95 \cdot (\log\text{pop}_i - \overline{\log\text{pop}}_{G=g})$
というコホート内中心化が **無加重 ≈ 0 を強制** する仕掛け
- 年トレンド `gamma_t` は 2014 で +9 上昇させ、論文 Table 2 と同じ大きさにする
```{r}
#| label: simulate-data
#| code-fold: true
#| code-summary: "シミュレーションデータ生成コード(クリックで展開)"
set.seed(2025)
# ===========================================================================
# (1) コホート構成 — 論文 Table 1 に整合させる
# ===========================================================================
# 全 2,604 郡を 5 つのコホートに分配
cohort_sizes <- c(
"2014" = 978, # 36% の郡
"2015" = 198, # 6%
"2016" = 132, # 4%
"2019" = 132, # 4%
"Inf" = 1164 # 31% (非拡大)
)
# 合計
N_county <- sum(cohort_sizes) # 2604
# 各郡にコホートを割り当て
cohort_vec <- rep(names(cohort_sizes), times = cohort_sizes)
G_numeric <- ifelse(cohort_vec == "Inf", Inf, as.numeric(cohort_vec))
# ===========================================================================
# (2) 郡レベルの不変属性(2013 年人口・共変量)
# ===========================================================================
# 人口:対数正規。コホートごとに分布をずらして「2014 群は人口比 45% を占める」を作る
# 2014 群は大都市が多い(中央値 log(pop) を高く)
pop_mu_by_cohort <- c(
"2014" = 11.0, "2015" = 9.8, "2016" = 9.5,
"2019" = 9.5, "Inf" = 10.0
)
pop_sd_by_cohort <- c(
"2014" = 1.7, "2015" = 1.5, "2016" = 1.4,
"2019" = 1.4, "Inf" = 1.6
)
log_pop <- mapply(function(g, n) rnorm(n, pop_mu_by_cohort[g], pop_sd_by_cohort[g]),
g = names(cohort_sizes), n = cohort_sizes,
SIMPLIFY = FALSE
) |> unlist()
pop_2013 <- exp(log_pop)
# 共変量:コホートごとに平均をずらして論文 Table 4 の不均衡を再現
gen_cov <- function(g, n) {
if (g == "Inf") { # 非拡大群:より白人比率高、貧困率高、所得低
tibble(
pct_female = rnorm(n, 49.4, 1.2),
pct_white = pmin(pmax(rnorm(n, 81.6, 14), 0), 100),
pct_hispanic = pmin(pmax(rnorm(n, 9.6, 9), 0), 100),
unemployment = pmax(rnorm(n, 7.6, 2.2), 0.5),
poverty = pmax(rnorm(n, 19.3, 5.5), 1),
median_inc = pmax(rnorm(n, 43.0, 11), 12) # 千ドル単位
)
} else { # 拡大群:白人率低め、貧困率低め、所得高め
tibble(
pct_female = rnorm(n, 49.3, 1.3),
pct_white = pmin(pmax(rnorm(n, 90.5, 8), 0), 100),
pct_hispanic = pmin(pmax(rnorm(n, 8.2, 9), 0), 100),
unemployment = pmax(rnorm(n, 8.0, 2.2), 0.5),
poverty = pmax(rnorm(n, 16.5, 5.0), 1),
median_inc = pmax(rnorm(n, 48.0, 12), 12)
)
}
}
covs_2013 <- purrr::map2_dfr(names(cohort_sizes), cohort_sizes, gen_cov)
# 共変量の 2014 年への変化(時変共変量を作るため)
delta_cov <- tibble(
d_pct_white = rnorm(N_county, -0.21, 0.4),
d_pct_hispanic = rnorm(N_county, 0.22, 0.4),
d_unemp = rnorm(N_county, -1.2, 0.7),
d_poverty = rnorm(N_county, -0.4, 0.6),
d_median_inc = rnorm(N_county, 1.1, 1.0)
)
# 郡レベルのデータフレーム
county_df <- tibble(
county_id = 1:N_county,
cohort_str = cohort_vec,
G = G_numeric,
pop_2013 = pop_2013
) |>
bind_cols(covs_2013) |>
bind_cols(delta_cov)
# ===========================================================================
# (3) パネル(11 年)に展開し、アウトカムを生成
# ===========================================================================
years <- 2009:2019
panel <- tidyr::expand_grid(county_id = 1:N_county, year = years) |>
left_join(county_df, by = "county_id")
# (a) 郡固有のベースライン死亡率 alpha_i
# - 非拡大群:平均 474、SD 80(論文の数値より少し控えめ)
# - 拡大群 :平均 419、SD 70
# - 人口が大きいほど死亡率は低い(都市部効果)→ 加重平均が無加重より低くなる
alpha_base <- ifelse(panel$G == Inf,
474 - 9 * (log(panel$pop_2013) - 10), # 非拡大
419 - 9 * (log(panel$pop_2013) - 10)
) # 拡大
alpha_noise <- rnorm(nrow(panel), 0, 1)
# 郡固有のレベルは年で変わらない→ 郡 ID 単位で再現できるよう郡内一定にする
alpha_noise_county <- tibble(
county_id = 1:N_county,
alpha_n = rnorm(N_county, 0, 50)
)
panel <- panel |>
left_join(alpha_noise_county, by = "county_id") |>
mutate(alpha = alpha_base + alpha_n)
# (b) 年トレンド gamma_t(2013→2014 で +9 ≒ 論文 Table 2 の数値)
# index: 2009 2010 2011 2012 2013 2014 2015 2016 2017 2018 2019
gamma_t <- tibble(
year = years,
gamma = c(-30, -20, -10, -5, 0, 9, 11, 13, 15, 17, 19)
)
panel <- panel |> left_join(gamma_t, by = "year")
# (c) 真の ATT — 2014 拡大群のみ、人口が大きい郡ほどマイナスの効果
# 各コホート内で「無加重平均 = 0」になるよう、コホート内中心化した人口偏差を使う
# こうすると 加重平均だけが負(人口大の郡で効果が大)になり、論文 Table 2 を再現
cohort_logpop_mean <- panel |>
group_by(G) |>
summarise(logpop_g_mean = mean(log(pop_2013)), .groups = "drop")
panel <- panel |>
left_join(cohort_logpop_mean, by = "G") |>
mutate(
pop_dev_g = log(pop_2013) - logpop_g_mean, # コホート内中心化
is_treated_post = (G != Inf) & (year >= G),
true_att = case_when(
G == 2014 & year >= 2014 ~ -0.95 * pop_dev_g,
G == 2015 & year >= 2015 ~ -0.5 * pop_dev_g,
G == 2016 & year >= 2016 ~ -0.3 * pop_dev_g,
G == 2019 & year >= 2019 ~ 0,
TRUE ~ 0
)
)
# (d) 時変共変量を線形に作る(2013→2014 で delta_cov 分動かす)
panel <- panel |>
mutate(
yr_idx = year - 2013,
poverty_t = poverty + d_poverty * yr_idx * 0.3,
unemp_t = unemployment + d_unemp * yr_idx * 0.3,
median_inc_t = median_inc + d_median_inc * yr_idx * 0.3
)
# (e) アウトカム
panel <- panel |>
mutate(
epsilon = rnorm(n(), 0, 6), # ノイズ SD を控えめにして信号(ATT ≈ -2.6)を視認しやすく
mortality_rate = alpha + gamma + true_att + epsilon
)
# 最終的に使う列
did_data <- panel |>
select(county_id, year, G, cohort_str, pop_2013,
pct_female, pct_white, pct_hispanic,
unemployment = unemp_t, poverty = poverty_t,
median_income = median_inc_t,
mortality_rate
) |>
mutate(
D = ifelse(G != Inf, 1L, 0L), # 「いずれ介入される」ダミー
treated_now = ifelse(G != Inf & year >= G, 1L, 0L), # 「現時点で介入中」ダミー
G_for_did = ifelse(G == Inf, 0, G) # `did` パッケージは Inf を 0 で扱う
)
cat("生成完了:", nrow(did_data), "行 ×", ncol(did_data), "列\n")
cat("郡数:", length(unique(did_data$county_id)), "\n")
cat("年数:", length(unique(did_data$year)), "\n")
```
### 生成したデータの概要
シミュレーションで作った `did_data` は、論文の実証例とほぼ同じ構造の
**バランス・パネルデータ**(全郡 × 全年に欠損なし)です。
| 項目 | 値 / 内容 |
|---|---|
| **総観測数** | 2,604 郡 × 11 年 = **28,644 行** |
| **観測単位** | 米国の「郡(county)」を模した架空の単位 |
| **観測期間** | 2009 年 〜 2019 年(11 年間) |
| **コホート構成** | 2014 年拡大 978 郡 / 2015 年 198 / 2016 年 132 / 2019 年 132 / **非拡大 1,164 郡** |
#### 主な変数
| 変数名 | 種別 | 意味 |
|---|---|---|
| `county_id` | ID | 1〜2604 の郡 ID(時間不変) |
| `year` | 時間 | 2009〜2019 |
| `G` | 介入時期 | コホート(2014/2015/2016/2019/`Inf`)。`Inf` は非拡大 |
| `cohort_str` | 介入時期 | `G` の文字列版(表示用) |
| `G_for_did` | 介入時期 | `did` パッケージ用に `Inf` を `0` に置換 |
| `D` | 介入ダミー | いずれ介入を受けるなら 1(`G != Inf`) |
| `treated_now` | 介入ダミー | その年に介入下にあるなら 1(`G != Inf & year >= G`) |
| `mortality_rate` | アウトカム | 20〜64 歳成人の **粗死亡率(10 万人あたり)**。論文と同様に概ね 400〜500 のレンジ |
| `pop_2013` | ウェイト | 郡の 2013 年成人人口(対数正規分布、コホート間で平均をずらしている) |
| `pct_female`, `pct_white`, `pct_hispanic` | 共変量 | 人口統計学的属性(%) |
| `unemployment`, `poverty`, `median_income` | 共変量 | 経済的属性(失業率%、貧困率%、中央世帯所得 千ドル) |
#### データ生成過程(DGP)の核心 — 「無加重 ≈ 0、人口加重 ≈ −2.6」を再現する仕掛け
シミュレーションは次の方程式で `mortality_rate` を生成しています:
$$
Y_{i,t} = \underbrace{\alpha_i}_{\text{郡固有レベル}}
+ \underbrace{\gamma_t}_{\text{年トレンド}}
+ \underbrace{\tau_{i,t}}_{\text{真の介入効果}}
+ \underbrace{\varepsilon_{i,t}}_{\text{ノイズ}}
$$
- $\alpha_i$ は郡固有(時間不変)。コホートにより平均を変えて Table 4 のような不均衡を作る。
人口の大きい郡では低めにし、論文と同様に「都市部ほど死亡率が低い」状況を作る
- $\gamma_t$ は全郡共通の年トレンドで、$\gamma_{2014} - \gamma_{2013} = +9$ にし論文の「年あたり +9」を再現
- **真の介入効果** $\tau_{i,t}$ は 2014 拡大群でのみ非ゼロで、
$\tau_{i,t} = -0.95 \cdot (\log(pop_i) - \overline{\log(pop)}_{G=2014})$ と
**コホート内で人口偏差に比例**。これが鍵:
- 各コホートで中心化しているので **無加重平均 = 0**(無加重 DiD ≈ 0)
- 人口の大きい郡で効果が強くなるので **人口加重平均 ≈ −2.74**(人口加重 DiD ≈ −2.6)
- $\varepsilon_{i,t}$ は SD=6 のホワイトノイズ
つまり「**大きい郡(都市)ほどメディケイドの効果が大きい**」と仮定したシミュレーションで、
これにより「加重するかしないかで結論が変わる」現象を教育的に再現できます。
#### データの中身を実際に見る
```{r}
#| label: peek-data
head(did_data, 3) |>
gt::gt() |>
gt::fmt_number(decimals = 2)
```
### 較正値の確認(Table 2 の数値を再現できているか)
#### 出力の読み方
下の `calib` データフレームは **群 × 年(計 4 行)**の人口加重・無加重死亡率です。
そこから DiD = (拡大群の前後変化) − (非拡大群の前後変化) を 2 種類計算。
最終 2 行の `cat` 出力で「論文 Table 2 の +0.1 / −2.6 と比較しての達成値」が表示されます。
#### 実際の結果の解釈
```{r}
#| label: calibration-check
calib <- did_data |>
filter(year %in% c(2013, 2014), G %in% c(2014, Inf)) |>
group_by(D, year) |>
summarise(
unweighted = mean(mortality_rate),
weighted = weighted.mean(mortality_rate, pop_2013),
.groups = "drop"
)
calib
did_unw <- (calib$unweighted[calib$D == 1 & calib$year == 2014] - calib$unweighted[calib$D == 1 & calib$year == 2013]) -
(calib$unweighted[calib$D == 0 & calib$year == 2014] - calib$unweighted[calib$D == 0 & calib$year == 2013])
did_wgt <- (calib$weighted[calib$D == 1 & calib$year == 2014] - calib$weighted[calib$D == 1 & calib$year == 2013]) -
(calib$weighted[calib$D == 0 & calib$year == 2014] - calib$weighted[calib$D == 0 & calib$year == 2013])
cat(sprintf("無加重 2×2 DiD = %+.2f (論文 Table 2: +0.1)\n", did_unw))
cat(sprintf("人口加重 2×2 DiD = %+.2f (論文 Table 2: -2.6)\n", did_wgt))
```
実値:
- **無加重 2×2 DiD = −0.23**(論文 Table 2: +0.1。差は ±0.5 以内で許容)
- **人口加重 2×2 DiD = −1.88**(論文 Table 2: −2.6。差は ±1 以内で許容)
論文 Table 2 とほぼ同じ符号・大きさの推定値が得られています。
**シミュレーションでも、加重の有無で結論が変わる** ことが確認できました。
この乖離は §3.3 以降で繰り返し参照されます。
## Table 1 の再現:コホート別の郡数・人口シェア
```{r}
#| label: tbl-1
#| tbl-cap: "ACA 下でのメディケイド拡大(Table 1 再現)"
did_data |>
filter(year == 2013) |>
group_by(cohort_str) |>
summarise(
n_counties = n(),
share_counties = n() / N_county,
share_pop = sum(pop_2013) / sum(did_data$pop_2013[did_data$year == 2013])
) |>
mutate(cohort_str = factor(cohort_str, levels = c("2014", "2015", "2016", "2019", "Inf"))) |>
arrange(cohort_str) |>
gt() |>
fmt_number(columns = n_counties, decimals = 0) |>
fmt_percent(columns = c(share_counties, share_pop), decimals = 1) |>
cols_label(
cohort_str = "拡大年", n_counties = "郡数",
share_counties = "郡シェア", share_pop = "人口シェア"
)
```
---
# 2×2 DiD デザイン(論文 §3)
論文 §3 では、最も単純な **2 つの時期 × 2 つの群** の設定を扱います。
具体的には:
- 介入群: **2014 年拡大群(978 郡)**
- 対照群: **2019 年まで拡大しなかった群(1,222 郡 = 2019 拡大群 + 非拡大群)**
- 期間: **2013 年と 2014 年のみ**
これは複雑な DiD デザインを理解する **構成要素(building block)** にもなります。
## 3.1 因果効果と目標パラメータ:ATT
DiD では通常 **ATT(t)** = 「時刻 $t$ における、介入群にとっての平均介入効果」を推定します:
$$
ATT(t) = E[Y_{i,t}(1) - Y_{i,t}(0) \mid D_i = 1]
$$
$Y(1)$ は介入を受けた場合、$Y(0)$ は受けなかった場合のアウトカム。
実際に観察できるのは片方のみで、$Y(0)$ は反事実です。
## 3.2 識別仮定:平行トレンド
平行トレンド仮定がなぜ必要かを可視化してみましょう。
```{r}
#| label: fig-pt-visual
#| fig-cap: "平行トレンドの可視化(2014 群 vs 非拡大群、2009-2014 年)"
#| fig-height: 4.5
pt_plot <- did_data |>
filter(year <= 2014, G %in% c(2014, Inf)) |>
group_by(D, year) |>
summarise(Y = weighted.mean(mortality_rate, pop_2013), .groups = "drop") |>
mutate(group = ifelse(D == 1, "2014 拡大群", "非拡大群"))
ggplot(pt_plot, aes(year, Y, color = group, group = group)) +
geom_line(linewidth = 1.1) +
geom_point(size = 2.5) +
geom_vline(xintercept = 2013.5, linetype = "dashed", color = "grey50") +
annotate("text",
x = 2013.5, y = max(pt_plot$Y), label = " 介入開始",
hjust = 0, color = "grey30"
) +
scale_x_continuous(breaks = 2009:2014) +
labs(
x = "年", y = "人口加重 死亡率(/10万)",
color = NULL,
title = "介入前のトレンドが「平行」かを目視で確認"
) +
theme(legend.position = "bottom")
```
二本の線が **同じ傾き** で動いていれば、平行トレンド仮定は妥当そうに見えます。
ただし平行トレンドは反事実に関する仮定なので、過去のトレンドが平行でも将来も平行とは限らない点に注意。
## 3.3 推定と推論:4 つの平均か 1 つの回帰か
#### このセクションで何をするか
ここからいよいよ **実際に ATT を推定** します。2×2 DiD には 2 つの計算方法があり、
**両方とも数値的に等価** であることを確認していきます。
| 方法 | やり方 | 強み |
|---|---|---|
| **① 手計算(4 つの平均)** | 介入 × 介入前後の 4 平均を出し、「差の差」を引き算 | 中身が透明、Excel でも計算可能 |
| **② 回帰(1 つの式)** | 交互作用項 `D × Post` 入りの線形回帰を 1 本走らせる | 標準誤差・信頼区間が自動で出る、共変量も拡張容易 |
#### こんなときに必要
- **介入前後 2 期間 × 2 群** の最も単純な設定で因果効果を推定したい
- 政策評価で「介入直後の影響」を見たい(動学効果は §5 で扱う)
- 回帰なしの **手計算** で結果が再現できることを示し、回帰の透明性を担保したい
- **無加重 vs 人口加重** の影響を直接比較し、結果の頑健性を確認したい
##### 確認すること(3 点)
1. **手計算と回帰の DiD 推定値が一致するか?**
$\widehat{ATT}^{4\text{平均}} = \widehat{\beta}^{2\times 2}$ になることを目で見て確認
2. **3 つの回帰仕様(交互作用 / 二方向固定効果 / 長期差分)が同じ値を返すか?**
バランス・パネルでは、これらは **全て同じ DiD 推定値** になります
3. **無加重と人口加重で結論がどう変わるか?**
無加重 ≈ 0(効果なし)、人口加重 ≈ −1.9(死亡率低下)。
これは「**大きな郡ほど効果が大きい**」シミュレーション設計の帰結
それでは順番に見ていきます。
### 手計算による 2×2 DiD(4 つの平均)
論文 Table 2 を再現します。2 行 × 2 群 × (無加重・加重) で 4 平均を計算:
```{r}
#| label: tbl-2-fourmeans
#| tbl-cap: "単純な 2×2 DiD(Table 2 再現)"
means_22 <- did_data |>
filter(year %in% c(2013, 2014), G == 2014 | G == Inf) |>
mutate(group = ifelse(D == 1, "Expansion", "No Expansion")) |>
group_by(group, year) |>
summarise(
unw = mean(mortality_rate),
wtd = weighted.mean(mortality_rate, pop_2013),
.groups = "drop"
)
# 整形:行=年、列=群×加重
tab2 <- means_22 |>
pivot_wider(names_from = group, values_from = c(unw, wtd))
tab2_with_diff <- tab2 |>
mutate(
`Gap_unw` = unw_Expansion - `unw_No Expansion`,
`Gap_wtd` = wtd_Expansion - `wtd_No Expansion`
)
# Trend 行を計算
trend_row <- tibble(
year = "Trend",
unw_Expansion = diff(tab2$unw_Expansion),
`unw_No Expansion` = diff(tab2$`unw_No Expansion`),
wtd_Expansion = diff(tab2$wtd_Expansion),
`wtd_No Expansion` = diff(tab2$`wtd_No Expansion`),
Gap_unw = diff(tab2$unw_Expansion) - diff(tab2$`unw_No Expansion`),
Gap_wtd = diff(tab2$wtd_Expansion) - diff(tab2$`wtd_No Expansion`)
)
tab2_final <- bind_rows(
tab2_with_diff |> mutate(year = as.character(year)),
trend_row
)
tab2_final |>
gt() |>
fmt_number(columns = -year, decimals = 1) |>
tab_spanner(label = "無加重", columns = c(unw_Expansion, `unw_No Expansion`, Gap_unw)) |>
tab_spanner(label = "人口加重", columns = c(wtd_Expansion, `wtd_No Expansion`, Gap_wtd)) |>
cols_label(
unw_Expansion = "拡大群", `unw_No Expansion` = "非拡大群", Gap_unw = "差",
wtd_Expansion = "拡大群", `wtd_No Expansion` = "非拡大群", Gap_wtd = "差"
)
```
**最下行の "Trend" の "差" 列が DiD 推定値です。** 論文と同様、
無加重 ≈ 0、人口加重 ≈ −2.6 となります。
### 回帰による 2×2 DiD
論文 (3.7) 式は、次の交互作用回帰と数値的に等価:
$$
Y_{i,t} = \beta_0 + \beta_1 \mathbf{1}\{D_i = 1\} + \beta_2 \mathbf{1}\{t = 2014\} + \beta^{2\times 2}(\mathbf{1}\{D_i = 1\} \cdot \mathbf{1}\{t = 2014\}) + \varepsilon_{i,t}
$$
::: {.callout-note}
## 記号:$\mathbf{1}\{\cdot\}$(指示関数 / インディケーター)
太字の「**1**」と中括弧 `{ }` を組み合わせた $\mathbf{1}\{\text{条件}\}$ は
**指示関数(示性関数, indicator function)** という記号で、
$$
\mathbf{1}\{\text{条件}\} = \begin{cases} 1 & \text{条件が真のとき} \\ 0 & \text{条件が偽のとき} \end{cases}
$$
という二値変数を表します。R で言えば `as.integer(条件)` と同じ動作です。
この式の場合:
- $\mathbf{1}\{D_i = 1\}$ … 介入群なら 1、対照群なら 0(= `D` 列そのもの)
- $\mathbf{1}\{t = 2014\}$ … 2014 年(介入後)なら 1、2013 年(介入前)なら 0(= `post` 列)
- $\mathbf{1}\{D_i = 1\} \cdot \mathbf{1}\{t = 2014\}$ … **介入群 × 介入後** のときだけ 1(= `D × post` 交互作用項)
つまり「中括弧の中の条件が成り立つときだけ 1 を返すスイッチ」と思えば OK です。
:::
$\beta^{2\times 2}$ が DiD 推定量です。
::: {.callout-note}
## 用語:TWFE (Two-Way Fixed Effects, 二方向固定効果)
回帰モデルに **単位固定効果**(郡ダミー)と **時間固定効果**(年ダミー)の両方を入れる仕様。
#### 上の交互作用回帰式との対応
上の式 $Y_{i,t} = \beta_0 + \beta_1 \mathbf{1}\{D_i = 1\} + \beta_2 \mathbf{1}\{t = 2014\} + \beta^{2\times 2}(\mathbf{1}\{D_i = 1\} \cdot \mathbf{1}\{t = 2014\}) + \varepsilon_{i,t}$
を 2 期 2 群で書き直すと、次の TWFE 式と等価になります:
$$
Y_{i,t} = \underbrace{\alpha_i}_{\text{単位固定効果}} + \underbrace{\lambda_t}_{\text{時間固定効果}} + \beta^{2\times 2}(\mathbf{1}\{D_i = 1\} \cdot \mathbf{1}\{t = 2014\}) + \varepsilon_{i,t}
$$
| 元の式の項 | TWFE 式での対応 | 役割 |
|---|---|---|
| $\beta_0$ | $\alpha_i$ と $\lambda_t$ に吸収 | 「対照群 × 2013 年」のベースライン |
| $\beta_1 \mathbf{1}\{D_i=1\}$ | **$\alpha_i$(単位固定効果)に吸収** | 介入群と対照群の **平均レベル差**。2 群しかないので「群ダミー = 郡ダミーの粗い版」。多単位設定では各郡ごとに切片 $\alpha_i$ を持たせるため $\beta_1$ は不要に |
| $\beta_2 \mathbf{1}\{t=2014\}$ | **$\lambda_t$(時間固定効果)に吸収** | 2013→2014 の **全体的なシフト**。2 期しかないので「年ダミー = $\lambda_t$ の粗い版」。多期間設定では各年ごとに切片 $\lambda_t$ を持たせるため $\beta_2$ は不要に |
| $\beta^{2\times 2}(\mathbf{1}\{D_i=1\} \cdot \mathbf{1}\{t=2014\})$ | そのまま残る | **介入群 × 介入後** にだけつく項 → これが DiD 推定量 |
**要点**: 2 群しかなければ $\beta_1 \mathbf{1}\{D_i=1\}$ で群レベル差を表現できますが、
3 群以上(段階的設定)になると **群ごとに別々のレベル**($\alpha_i$ 群ダミー)が必要になります。
同様に、2 期しかなければ $\beta_2 \mathbf{1}\{t=2014\}$ で時間シフトを表現できますが、
多期間では **年ごとに別々のレベル**($\lambda_t$ 年ダミー)が必要になります。
TWFE 式はその一般化です。
R では `feols(Y ~ D:post | county_id + year, ...)` のように
`|` の右側に固定効果を書きます(`m2` の仕様、すぐ下のコードを参照)。
2 期 2 群の単純設定では交互作用回帰と **数値的に等価**(下の Table 3 の (1) vs (2) で確認)
ですが、後で見るように **段階的設定では問題が生じる** ことが近年判明しました(§5.3 で詳説)。
:::
::: {.callout-note}
## 用語:標準誤差(Standard Error, SE)とクラスタロバスト標準誤差
**標準誤差(SE)** は推定値の **ばらつき(標本サイズや変動から来る不確実性)** を測る量。
本ノートでは以後 `SE` と略記します。
**クラスタロバスト SE** は「同じ集団内では誤差が相関しうる」と認めた上で計算する SE。
論文では **郡レベル** でクラスタリング(同じ郡の異なる年は相関しうる)。
**95% 信頼区間(95% Confidence Interval, 95% CI)** は通常 ATT ± 1.96 × SE で計算。
本ノートでは以後 `95% CI` と略記します。「CI が 0 を含まない」= 統計的有意(おおむね p < 0.05 相当)。
:::
```{r}
#| label: tbl-3-regression
#| tbl-cap: "回帰 2×2 DiD(Table 3 再現)"
sub <- did_data |>
filter(year %in% c(2013, 2014), G == 2014 | G == Inf) |>
mutate(post = (year == 2014) * 1L)
# 無加重
m1 <- feols(mortality_rate ~ D * post, data = sub, cluster = "county_id")
m2 <- feols(mortality_rate ~ D:post | county_id + year, data = sub, cluster = "county_id")
# 長期差分(郡ごとに ΔY を D に回帰)
diff_df <- sub |>
arrange(county_id, year) |>
group_by(county_id) |>
summarise(
delta_Y = mortality_rate[year == 2014] - mortality_rate[year == 2013],
D = first(D), pop = first(pop_2013)
)
m3 <- feols(delta_Y ~ D, data = diff_df, cluster = "county_id")
# 人口加重版
m4 <- feols(mortality_rate ~ D * post, data = sub, weights = ~pop_2013, cluster = "county_id")
m5 <- feols(mortality_rate ~ D:post | county_id + year, data = sub, weights = ~pop_2013, cluster = "county_id")
m6 <- feols(delta_Y ~ D, data = diff_df, weights = ~pop, cluster = "county_id")
modelsummary(
list(
"無加重(1)" = m1, "無加重(2)" = m2, "無加重(3)" = m3,
"加重(4)" = m4, "加重(5)" = m5, "加重(6)" = m6
),
fmt = 1,
stars = TRUE,
coef_omit = "Intercept",
gof_omit = "AIC|BIC|RMSE|R2 Adj|R2 Within|Log",
output = "gt"
)
```
#### 表の見方
まずこの表全体の構造を整理します。
##### 列(6 つの推定モデル)
横方向に 6 列あり、それぞれが **異なる回帰仕様** に対応しています:
| 列 | ラベル | 加重 | 仕様 | R コード |
|---|---|---|---|---|
| (1) | 無加重(1) | なし | **交互作用** | `feols(Y ~ D * post)` |
| (2) | 無加重(2) | なし | **二方向固定効果(TWFE)** | `feols(Y ~ D:post \| county_id + year)` |
| (3) | 無加重(3) | なし | **長期差分** | `feols(ΔY ~ D)`($\Delta Y = Y_{2014} - Y_{2013}$) |
| (4) | 加重(4) | 人口 | 交互作用 | (1)+ `weights = pop_2013` |
| (5) | 加重(5) | 人口 | TWFE | (2)+ `weights = pop_2013` |
| (6) | 加重(6) | 人口 | 長期差分 | (3)+ `weights = pop_2013` |
左 3 列が **無加重**、右 3 列が **人口加重** です。
##### 行(係数と統計量)
縦方向に並ぶ行は次のとおり:
| 行 | 意味 |
|---|---|
| `D` | 介入群ダミー(対照群比の **レベル差**)。(1)(4) と (3)(6) の長期差分でのみ登場 |
| `post` | 介入後年ダミー(2013→2014 の **シフト**)。(1)(4) でのみ登場 |
| **`D × post`(または `D:post`)** | **DiD 推定量 $\beta^{2\times 2}$** ← これが本命 |
| `Num.Obs.` | 観測数 |
| `R2` | 決定係数 |
| `Std.Errors` | 標準誤差のクラスタリング水準(郡レベル) |
各係数の **下の括弧内の数値が標準誤差**、係数右の **`*`, `**`, `***`** が有意水準
(順に p < 0.1, 0.05, 0.01)です。
##### TWFE と長期差分で `D` 行が消える理由
- TWFE(2)(5):郡固定効果 $\alpha_i$ が `D` を吸収してしまうので、`D` 単独の係数は出ない
- 長期差分(3)(6):差分してから D に回帰するので、定数項と `D` だけ残る
#### 表のどこを見て何を読み取るか
注目すべきは **`D × post` の係数行**(または m3, m6 では `D` 行)です。これが各仕様の
DiD 推定量に対応します。
##### ① 3 つの仕様が等価であることの確認
- 列 (1)・(2)・(3) [無加重] で `D × post`(または `D`)の係数を比較 → **3 つとも同じ値**(≈ −0.2)
- 列 (4)・(5)・(6) [人口加重] で同じく `D × post` の係数を比較 → **3 つとも同じ値**(≈ −1.9)
つまり横一列に並んだ 3 つの値が常にぴったり一致しています。**標準誤差(括弧内)も同じ**。
##### ② 加重・無加重で結論が大きく変わることの確認
- 無加重(1)〜(3): `D × post` ≈ **−0.2**(統計的にゼロと区別できない)
- 人口加重(4)〜(6): `D × post` ≈ **−1.9**(より大きな低下、有意水準は星マークで確認可能)
横方向に 6 列を一気に見ると「左 3 列と右 3 列でガラッと値が変わる」ことが分かります。
#### なぜこれが重要か
**(A) "等価性" を確認することの意義**
- 手計算と回帰の値が一致することから、**回帰は単に 4 平均の差を計算する便利な道具** であって、
何か神秘的な統計魔法をかけているわけではないと納得できます
- 3 つの仕様が同値ということは、**実務では好きな書き方を選べる** ということ
(固定効果は標準誤差クラスター化が綺麗、長期差分はデータ量が半分で済む、など)
- ただしこの等価性は **2 期 2 群** の単純設定でしか成り立ちません。§5 で見るように、
段階的設定では TWFE と「素直な DiD」が乖離します — その対比の出発点として
「2×2 では一致する」ことを押さえる必要があります
**(B) "加重で結論が変わる" ことの意義**
- 加重するかしないかは技術的細部ではなく **どちらの推定量を見せるべきかの判断**
**前提として、本データは「郡レベルの集計死亡率」**(個人レベルデータではない)です。
したがってここで議論しているのは **両方とも集団レベル(集計)の平均介入効果**であって、
**個人の介入効果**(ITE: individual treatment effect)ではありません。
個人レベル効果を厳密に推定するには個人データが必要で、miller2021medicaid のような
研究はそれを使って制限付きデータで分析しています。
**それぞれの結論を正確に書くと:**
- **無加重(列 1–3)**:
> 「**全 2,604 郡を 1 票ずつ等しく扱う集計平均で見たとき**、
> メディケイド拡大の効果は 10 万人あたり **約 −0.2 人**
> (統計的に有意ではない)」
- **人口加重(列 4–6)**:
> 「**各郡を 2013 年成人人口でウェイト付けた集計平均で見たとき**、
> メディケイド拡大の効果は 10 万人あたり **約 −1.9 人**
> (`*` で示される水準で統計的に有意)」
両者は **どちらも同じデータから出てきた正しい計算結果** ですが、平均化の方法が異なります:
| | 無加重 | 人口加重 |
|---|---|---|
| 平均化単位 | 1 郡 = 1 票 | 1 人 = 1 票(郡内の人数で重み付け) |
| 推定対象 | 全郡集合の平均的な郡レベル死亡率変化 | 米国成人人口全体の集計死亡率変化 |
| 数値 | ≈ −0.2(有意でない) | ≈ −1.9(有意) |
#### では、どちらを見るべきか?
**ほとんどの公衆衛生・政策評価の文脈では、知りたいのは「人々への効果」なので
人口加重の方が解釈的に意味があります。** 「ワイオミング州(人口少)もカリフォルニア州(人口多)
も 1 票」という無加重平均は、政策の人々への影響を測りたい目的とはずれます。
それでも **無加重を併記する** 理由は以下の通り:
1. **異質性の診断**: 無加重と加重で大きく違うなら「効果が郡サイズと相関している」
というシグナルになる(本ノートのシミュレーションがまさにこれ)
2. **大都市の支配を避けるロバストネスチェック**: カリフォルニア州 1 つで結果が決まっていないかを確認
3. **計量経済学的伝統**: 古典的 DiD は OLS 由来で無加重がデフォルトだったため
論文も両方報告するのが慣例
ただし「無加重 ≈ 0、加重 ≈ −2」を「結論は人による」と片付けるのではなく、
**人口加重を主結果に、無加重をロバストネスチェックに** という整理が現代の標準的な実務です。
論文(と本ノート)が両方を併記するのは、両者の **乖離自体が情報** だからです。
**注**: 本ノートのシミュレーションでは「大きい郡ほど効果が大きい」と DGP を設定したので
両者が乖離します。これは Borgschulte (2020) など実際のメディケイド研究でも観察される
パターンで、人口の多い都市部に介入の便益が集中する政策では同様の現象が起こりえます。
---
# 共変量を 2×2 DiD に組み込む(論文 §4)
平行トレンド仮定が成り立たない場合、**共変量で条件付けた平行トレンド** を仮定して
ATT を識別する方法が必要です。
::: {.callout-tip}
## 実務:こんなときは共変量調整を検討する
「とりあえず共変量を入れる」のではなく、**未介入アウトカム $Y(0)$ の動き方が群間で違いそう**
と疑う具体的な理由があるときに調整します。例えば:
**(1) 介入群と対照群でベースライン特性が大きく違う**(本ノートの例)
- メディケイド拡大州は非拡大州より **白人率が高く、貧困率が低く、所得が高い**(Table 4 参照)
- 既往研究から「貧困率が違う郡は死亡率の進化トレンドも違う」(Currie 2016)と分かっている
- ⇒ 貧困率・人種構成・所得を共変量に入れる
**(2) 「平均への回帰」が起こりそう**
- 介入直前にアウトカムが極端な値だった群が選ばれている
(例:「死亡率が急上昇した州だけが対策を打った」)
- 何もしなくても平均に戻る動きが残り、平行トレンドが破れる
- ⇒ 介入前のアウトカム水準やトレンド指標を共変量に
**(3) 時間とともに進む構造変化**
- 高齢化のスピード、産業構造の変化、都市化など、長期トレンドが地域で異なる
- 例:過疎化が早い農村 vs 人口流入が続く都市部 → 死亡率トレンドが平行になりにくい
- ⇒ 高齢化率、産業構成、都市化指標を共変量に
**(4) 介入の割当てがランダムでなく、観察可能な選別がある**
- 政策ターゲティング:「失業率の高い地域から優先的に職業訓練を導入」
- 病院・施設選択:「重症患者が高度医療機関に偏る」
- ⇒ 割当てに使われた変数を共変量に
#### 入れるべき共変量の条件(基本ルール)
| 入れるべき | 入れてはいけない |
|---|---|
| ✅ **介入前** に確定している値 | ❌ 介入の **結果** で変わる値(媒介変数。アウトカム同様、調整するとバイアスが入る) |
| ✅ 反事実 $Y(0)$ のトレンドと関連 | ❌ 介入と無関係でアウトカムにも無関係(精度が落ちるだけ) |
| ✅ 群間でバランスが取れていない(SMD > 0.2) | ❌ 既にバランスが取れている(調整しても変わらない) |
**避けるべき典型例**: 「介入後の保険加入率」を共変量に入れる
→ 介入の効果そのもので変動する変数なので、効果を吸収してしまう
:::
::: {.callout-note}
## 用語:条件付き平行トレンド(Conditional Parallel Trends)
「**共変量 X が同じならば**、介入群と対照群の反事実トレンドが等しい」と仮定する。
無条件の平行トレンドより弱い仮定です:
$$
E[Y_{i,t}(0) - Y_{i,t-1}(0) \mid D_i = 1, X_i] = E[Y_{i,t}(0) - Y_{i,t-1}(0) \mid D_i = 0, X_i]
$$
:::
## 4.1 共変量バランス:無条件平行トレンドは妥当か?
#### このセクションで何をするか
§3 では「介入群と対照群の反事実トレンドが平行」(無条件 PT)を仮定して 2×2 DiD を計算しました。
しかし観察研究では、介入の割当てがランダムでないため、両群で **観察可能な共変量 X**
(人口統計、経済指標等)が大きく異なることがよくあります。
このセクションでは:
1. **介入群と対照群の共変量を比較** し、**標準化平均差(SMD)** で不均衡を定量化
2. 不均衡があれば **無条件 PT は疑わしい** → §4.2 以降の **条件付き PT** に移行する根拠を確立
これは論文の **Table 4 再現** に相当します。
#### こんなときに必要
- **観察研究** で介入の割当てがランダムでない(政策評価、自然実験等)
- 平行トレンド仮定の **経験的妥当性** を読者に説得したい
- 査読で「なぜ共変量で条件付けたのか」を **データに基づいて** 説明する必要がある
- **どの共変量を後の調整で使うか** の最初のスクリーニング(ただし機械的に SMD で決めるのは
推奨しない — 下の callout 参照)
::: {.callout-note}
## 用語:標準化平均差(SMD, Normalized Difference)
両群の平均差をプールした SD で割った量。経験則として **|SMD| > 0.2 で不均衡** と判断します
(Austin 2009 等の臨床/疫学文献の慣例。より厳格には 0.1 を採用することもある)。
$$
\text{SMD} = \frac{\bar{X}_T - \bar{X}_C}{\sqrt{(s_T^2 + s_C^2)/2}}
$$
「t 検定の p 値」と違い、サンプルサイズに影響されない不均衡指標です。
:::
#### コードで何をしているか
| ステップ | コード | 何をしているか |
|---|---|---|
| 1 | `balance_stat()` 自作関数 | ある変数について、群別の平均と SD を計算し SMD を返す(無加重・加重両対応) |
| 2 | `bal_2013 <- filter(year == 2013, ...)` | 2013 年(介入前)のサブサンプルを抽出 |
| 3 | `sapply(vars, balance_stat, ...)` × 2 | 各共変量について **無加重** SMD と **人口加重** SMD を計算 |
| 4 | `mutate(... lag(.) ...)` で差分 | 2014−2013 の **共変量変化量** を計算 |
| 5 | `make_panel()` で整形 | 上段(2013 年水準)・下段(2014−2013 差分)の 2 パネル表 |
| 6 | `gt() \|> fmt_number()` | 表として表示 |
```{r}
#| label: tbl-4-balance
#| tbl-cap: "共変量バランス統計量(Table 4 再現)"
balance_stat <- function(x, group, w = NULL) {
if (is.null(w)) {
mu_t <- mean(x[group == 1])
mu_c <- mean(x[group == 0])
s_t <- sd(x[group == 1])
s_c <- sd(x[group == 0])
} else {
mu_t <- weighted.mean(x[group == 1], w[group == 1])
mu_c <- weighted.mean(x[group == 0], w[group == 0])
# 加重分散
wvar <- function(x, w) {
mu <- weighted.mean(x, w)
sum(w * (x - mu)^2) / sum(w)
}
s_t <- sqrt(wvar(x[group == 1], w[group == 1]))
s_c <- sqrt(wvar(x[group == 0], w[group == 0]))
}
smd <- (mu_t - mu_c) / sqrt((s_t^2 + s_c^2) / 2)
c(non_adopt = mu_c, adopt = mu_t, norm_diff = smd)
}
bal_2013 <- did_data |>
filter(year == 2013, G == 2014 | G == Inf)
vars <- c("pct_female", "pct_white", "pct_hispanic", "unemployment", "poverty", "median_income")
var_labels <- c("% Female", "% White", "% Hispanic", "Unemployment", "Poverty", "Median Income")
# 2013 レベル
levels_unw <- sapply(vars, function(v) balance_stat(bal_2013[[v]], bal_2013$D))
levels_wtd <- sapply(vars, function(v) balance_stat(bal_2013[[v]], bal_2013$D, bal_2013$pop_2013))
# 2014-2013 差分
bal_diff <- did_data |>
filter(year %in% c(2013, 2014), G == 2014 | G == Inf) |>
arrange(county_id, year) |>
group_by(county_id) |>
mutate(across(all_of(vars), ~ . - lag(.), .names = "d_{.col}")) |>
filter(year == 2014) |>
ungroup()
diff_unw <- sapply(paste0("d_", vars), function(v) balance_stat(bal_diff[[v]], bal_diff$D))
diff_wtd <- sapply(paste0("d_", vars), function(v) balance_stat(bal_diff[[v]], bal_diff$D, bal_diff$pop_2013))
# 整形
make_panel <- function(unw_mat, wtd_mat, panel_name) {
tibble(
Variable = var_labels,
`非拡大(無)` = unw_mat["non_adopt", ],
`拡大(無)` = unw_mat["adopt", ],
`SMD(無)` = unw_mat["norm_diff", ],
`非拡大(加)` = wtd_mat["non_adopt", ],
`拡大(加)` = wtd_mat["adopt", ],
`SMD(加)` = wtd_mat["norm_diff", ],
panel = panel_name
)
}
bal_tab <- bind_rows(
make_panel(levels_unw, levels_wtd, "2013 年水準"),
make_panel(diff_unw, diff_wtd, "2014-2013 差分")
)
bal_tab |>
group_by(panel) |>
gt() |>
fmt_number(decimals = 2, columns = -c(Variable, panel))
```
#### 出力の読み方
表は **2 パネル × 6 変数 × 6 列** の構造:
- **上段「2013 年水準」**: 介入前年の各変数の平均と SMD
- **下段「2014-2013 差分」**: 1 年間の **変化量** の群間比較
- 列の左 3 つ:**無加重**(非拡大平均 / 拡大平均 / SMD)
- 列の右 3 つ:**人口加重**(非拡大平均 / 拡大平均 / SMD)
**チェックポイント**:
1. **|SMD| > 0.2 の変数** はどれか?(慣例的閾値、§4.1 callout 参照)
2. **無加重と加重で SMD が大きく違う変数** は?(郡サイズと特性が相関しているシグナル)
3. **2014−2013 差分の SMD** が大きい変数は?(共変量自体のトレンドが群間で違う = 動的交絡の可能性)
4. SMD の **符号** は?(拡大群が高い場合 +、対照群が高い場合 −)
#### 実際の結果の解釈
実際の SMD(2013 年水準、無加重 / 加重):
| 変数 | 非拡大 | 拡大 | SMD(無加重) | 非拡大(加) | 拡大(加) | SMD(加重) |
|---|---|---|---|---|---|---|
| % Female | ≈ 49.4 | ≈ 49.3 | ≈ 0 | ≈ 49.4 | ≈ 49.3 | ≈ 0 |
| **% White** | 81.31 | 89.88 | **0.83** | 84.41 | 90.07 | **0.54** |
| % Hispanic | 10.25 | 8.81 | −0.19 | 9.73 | 9.97 | 0.03 |
| Unemployment | 7.63 | 8.08 | 0.21 | 7.54 | 7.44 | −0.05 |
| **Poverty** | 19.11 | 16.74 | **−0.44** | 19.31 | 16.58 | **−0.48** |
| **Median Income** | 42.45 | 48.00 | **0.48** | 40.11 | 47.61 | **0.60** |
**読み取り 1:大きな不均衡が 3 変数で確認される**
|SMD| > 0.2 の変数:
- **% White**:無加重 0.83(極めて大きい)、加重 0.54 → 拡大群は **白人率が大幅に高い**
- **Poverty**:無加重 −0.44、加重 −0.48 → 拡大群は **貧困率が低い**
- **Median Income**:無加重 0.48、加重 0.60 → 拡大群は **中央世帯所得が高い**
つまり拡大州は **白人比率が高く、相対的に裕福で貧困率が低い** という構造的差異があります。
**読み取り 2:加重で SMD が変わる変数**
% White は無加重 0.83 から加重 0.54 に **緩和** されました。これは「拡大群の中で
人口の少ない郡で白人率が極端に高い」ことを示唆します(都市部の方が多様)。
Unemployment は無加重 0.21 から加重 −0.05 に **符号反転**(似た意味で都市部効果)。
**読み取り 3:平行トレンド違反のリスク**
% White, Poverty, Median Income はいずれも **死亡率トレンドの強力な決定要因** です
(Currie 2016 等)。これらが群間で大きく異なるなら、メディケイドがなくても両群は
**異なる死亡率トレンド** を持っていた可能性 → **無条件 PT は疑わしい**。
→ これが §4.2 以降で **条件付き PT** を導入する根拠です。
#### 注意点
SMD は便利な診断ツールですが、「SMD が大きい変数を機械的に共変量に入れる」運用には
強い反対意見があります。次の callout で詳しく論じます。
::: {.callout-important}
## 個人的見解:SMD ではなくメカニズムで共変量を選ぼう
SMD はバランスを **診断** するには良い道具ですが、**「SMD が大きい変数 = 入れるべき共変量」**
という運用は推奨しません。
#### 理由
1. **SMD は反事実トレンドとアウトカムの関係を測っていない**
SMD は「群間で平均がどれだけ違うか」だけを測ります。しかし共変量調整の目的は
「未介入アウトカム $Y(0)$ のトレンドが群間で違いそうな源を吸収する」こと。
不均衡があっても $Y(0)$ のトレンドに無関係なら入れる必要はないし、
均衡していても重要なら入れるべきです。
2. **臨床研究の交絡調整と発想は同じ**
ランダム化比較試験以外で因果推論する際、臨床疫学では
**DAG(因果ダイアグラム)や領域知識で交絡因子を特定** し、それを調整します。
「SMD > 0.2 だから入れる」のではなく、
**「介入の割当てと $Y(0)$ の進化の両方に影響しうるメカニズムは何か?」**
を考えて選びます。
3. **SMD ベースの自動選択は post-treatment bias の危険**
時変共変量で SMD を見ると、介入後の値が含まれてしまい、
**介入の結果として変化した変数(媒介変数)を誤って "不均衡" と判定** することがあります。
これを共変量に入れると効果を吸収してしまいます。
#### 推奨される進め方
1. **メカニズムの記述**: 介入の割当てがどう決まるか、$Y(0)$ のトレンドを左右しうる要因は何か
を、領域知識と先行研究から書き出す(できれば DAG を描く)
2. **介入前に確定する変数だけを候補に**
3. **SMD は「想定外の不均衡がないか」のチェックに使う**
(たとえば想定通りの共変量で調整したのに、想定外の変数で大きな SMD が残るなら
モデルや母集団の見直し)
#### 本ノートの例での具体例
メディケイド拡大の文脈で、メカニズム的に妥当な共変量は:
- **貧困率・中央世帯所得**: 拡大の政治的判断にも、死亡率トレンドにも影響
- **失業率**: 経済不況は健康トレンドに影響、ACA 拡大決定とも相関
- **人種構成**: 健康格差の構造的要因。トレンドにも影響
- **年齢構成 / 高齢化率**: 死亡率トレンドの最大の決定要因(ノートでは省略)
逆に **SMD が大きくても入れるべきでない** 例:
- 介入後の保険加入率(媒介変数)
- 介入後の医療支出(媒介変数)
- 純粋な地理座標(緯度経度などはメカニズム的根拠が薄い)
「SMD で機械的に選ぶ」のではなく、「**なぜこの変数が関係するのか言えるか?**」を
問うのが因果推論的に健全です。
:::
## 4.2 共変量付き DiD:条件付き平行トレンド下での識別
DiDではATEではなくATTを推定する。
::: {.callout-note}
## 用語:なぜ DiD は ATE ではなく ATT を推定するのか?
| 量 | 定義 | 答える問い |
|---|---|---|
| **ATE** (Average Treatment Effect) | $E[Y_i(1) - Y_i(0)]$ — **全員**の平均介入効果 | 「もし**全員**にこの政策を施したら、平均何が起こるか?」 |
| **ATT** (Average Treatment effect on the Treated) | $E[Y_i(1) - Y_i(0) \mid D_i = 1]$ — **実際に介入を受けた人たち**の平均介入効果 | 「**実際に政策が当てられた人たち**にとって、政策の効果はどれくらいだったか?」 |
#### DiD が ATT しか識別できない理由
1. **平行トレンド仮定が制約するのは $Y(0)$ のトレンドだけ**
PT は「介入群と対照群の **未介入アウトカム** $Y(0)$ の動きが平行」と言うのみ。
$Y(1)$(介入アウトカム)については何も言いません。
そのため、対照群の $Y(1)$(=「もし対照群も介入を受けたら」)については推論できません。
2. **必要な反事実は ATT の方が少ない**
- ATT を計算するには **介入群の $Y(0)$** だけ分かれば十分
→ これは対照群の動きから PT を使って復元可能 ✓
- ATE を計算するには **介入群の $Y(0)$ + 対照群の $Y(1)$** の両方が必要
→ 対照群の $Y(1)$ は DiD では復元できない ✗
3. **政策評価としても ATT の方が直接的**
「**この政策が、実際にそれを受けた人々に対して効いたか?**」が普通の問い。
メディケイド拡大の評価なら「拡大州の住民にとって死亡率がどう変わったか」が知りたい
ことであり、「もし非拡大州も拡大していたら…」は別の(より強い)仮定が要る問い。
#### ATE が欲しい場合は?
- 効果の異質性がないと仮定すれば ATT = ATE になりますが、これは強い仮定です
- ランダム化比較試験(RCT)ではランダム割当てにより自動的に ATT = ATE
- 観察データで ATE を識別したいなら、**より強い仮定**(無作為割当てに近い ignorability)が必要
このため、DiD では **ATT を主結果に置き、政策が実際の対象者に与えた効果として解釈する**
のが標準的です。
:::
条件付き平行トレンドが成り立つとき、ATT は次のように識別可能です:
$$
ATT(t) = E\left[ Y_{i,t} - Y_{i,t-1} \mid D_i = 1 \right]
- E\left[ E[Y_{i,t} - Y_{i,t-1} \mid X_i, D_i = 0] \;\Big|\; D_i = 1 \right]
$$
#### 式の各記号の意味
| 記号 | 意味 |
|---|---|
| $i$ | 単位の添字(本ノートでは郡) |
| $t$ | 介入後の時点(例:$t=2014$) |
| $t-1$ | 介入前の時点(例:$t-1=2013$) |
| $Y_{i,t}$ | 単位 $i$ の時点 $t$ における観察アウトカム(郡 $i$ の年 $t$ の死亡率) |
| $Y_{i,t} - Y_{i,t-1}$ | 単位 $i$ の **アウトカム変化量** $\Delta Y_i$(2013→2014 の差分) |
| $D_i$ | 介入群ダミー:$D_i = 1$ なら介入群、$D_i = 0$ なら対照群 |
| $X_i$ | 単位 $i$ の **共変量ベクトル**(介入前に確定;例:貧困率・所得・人種構成など) |
| $E[\cdot]$ | 期待値(母集団平均) |
| $E[\cdot \mid D_i = 1]$ | **介入群に限った** 条件付き期待値(介入群のサブサンプルでの平均) |
| $E[\cdot \mid X_i, D_i = 0]$ | **対照群かつ共変量 $X_i$ で条件付け** した期待値(同じ X を持つ対照群の平均) |
#### 式の各項の意味
**第 1 項** $E\left[ Y_{i,t} - Y_{i,t-1} \mid D_i = 1 \right]$
= 介入群における **観察された** アウトカム変化の平均
= 「実際に介入群で何が起きたか」を表す
**第 2 項** $E\left[ E[Y_{i,t} - Y_{i,t-1} \mid X_i, D_i = 0] \;\Big|\; D_i = 1 \right]$
= 二重期待値(入れ子の期待値)。内側から読みます:
- 内側 $E[Y_{i,t} - Y_{i,t-1} \mid X_i, D_i = 0]$
= 対照群で、共変量 $X_i$ を条件として「ΔY の予測モデル」を作ったもの
(例:「貧困率 15%・白人率 80% なら ΔY ≈ +9」のような関数)
- 外側 $E[\cdot \mid D_i = 1]$
= その予測モデルを **介入群の X 分布で平均** する
(= 介入群と同じ共変量を持つ対照群サンプルで予測した ΔY の平均)
= **「もし介入群が介入を受けなかった場合に予想される ΔY」=反事実**
#### 全体の直感
つまり 「実際の ΔY(介入群)− 反事実 ΔY(介入群の X で対照群モデルから予測)」 を計算しています。
R 実装の対応:
1. 対照群のデータで $\Delta Y \sim X$ の回帰モデル $\hat{m}_0(X)$ を当てはめる
2. 介入群の各 $X_i$ を代入して反事実 $\hat{m}_0(X_i)$ を計算
3. 介入群の平均 $\overline{\Delta Y}$ から介入群平均の $\hat{m}_0$ を引く
これが **アウトカム回帰(Outcome Regression)推定量** で、§4.4 で実装します。
考え方を一言で言うと: **対照群で X から「ΔY の予測モデル」を作り、介入群の X で予測**。
それを「もし介入群が介入を受けなかった場合の ΔY」とみなす。
#### この識別式が暗黙に要求するもう一つの条件:「重なり」
上の識別式が **意味を持つ** ためには、条件付き平行トレンドだけでなく、もう一つ仮定が必要です。
内側の期待値 $E[\Delta Y \mid X_i, D_i = 0]$ を介入群の X 値で評価できなければなりません。
ところが **対照群に存在しない X 値**(例:「人口 100 万超の郡が介入群にしかいない」状況)では、
そもそも対照群のサブサンプルが存在せず、$\hat{m}_0(X_i)$ が外挿になってしまいます。
この問題を避けるため、**介入群と対照群の共変量分布が "重なっている"** ことを要求します。
これが次の「重なり(Overlap)」条件です:
::: {.callout-note}
## 用語:重なり(Overlap)/ 強重なり(Strong Overlap, SO)
「**すべての X 値で、介入群と対照群の両方に十分なサンプルがいる**」という条件。
極端な X(例:大都市は介入群、過疎地は対照群に偏在)では対照群サンプルが少なすぎて、
予測モデル $\hat{m}_0(X)$ を介入群の X に当てる際に **外挿** が必要になり、推定が不安定になります。
**式での表現**: 傾向スコア $P(D = 1 \mid X)$ が 0 や 1 に張り付かないこと:
$$
\epsilon < P(D = 1 \mid X) < 1 - \epsilon \quad \text{(ある } \epsilon > 0 \text{ について)}
$$
「強重なり(SO)」は、この不等式が $X$ の全領域で **一様に** 成り立つことを要求する強めの版で、
論文ではこちらを使います。
**実務でのチェック方法**:
- 傾向スコアのヒストグラムを介入群・対照群で重ねて描く
- 重ならない領域がある場合は、その領域を **トリミング**(分析から除外)するか、
そもそも比較不可能と判断する
メディケイド例では、$X$ に「中央世帯所得」を入れる場合、所得が極端に高い/低い郡の
サンプル数が両群で確保されているかを確認する必要があります。
:::
つまり 4.2 節の識別には、 **(i) 条件付き平行トレンド** に加え、
**(ii) 共変量分布の重なり** という 2 つの仮定がセットで必要、と覚えてください。
## 4.3 共変量付き DiD 推定:TWFE(注意点あり)
#### このセクションで何をするか
§4.1 で示した不均衡を踏まえ、**TWFE 回帰に共変量を加える** 伝統的なやり方を試します。
ただし TWFE+共変量は **異質効果がある場合に ATT を返さない** ことが近年明らかになっており、
ここでは典型的な 3 仕様(共変量なし / ベースライン X × Post / 時変 X)を並べて
**DiD 推定値がどれだけ動くか** を確認します。
これは §4.4 の二重頑健推定への伏線で、TWFE では仕様によって値が動いてしまうことを
具体的な数値で示し、**より頑健な手法が必要** な理由を確立します。
#### こんなときに必要
- 過去の論文で TWFE+共変量が使われており **再現と比較** したい
- 共変量の入れ方が結果にどれくらい影響するかの **感度分析** を実施したい
- TWFE の限界を理解したうえで §4.4 の **アウトカム回帰 / 逆確率ウェイト / 二重頑健** に進む準備
- 査読対応で「TWFE で十分では?」と問われたときの **限界の実証**
#### コードで何をしているか
共変量の入れ方には自由度があり、それぞれが暗黙に異なる仮定を置いています。
3 仕様を比較します:
| 仕様 | R モデル式 | 入れ方 | 暗黙の仮定 |
|---|---|---|---|
| **(A) X なし** | `Y ~ D:post \| county_id + year` | 共変量を入れない(無条件 PT) | 介入群と対照群の反事実トレンドは「素のまま」平行 |
| **(B) ベースライン × Post** | `Y ~ D:post + X_{2013}:post \| county_id + year` | **2013 年の** 共変量を post と交互作用 | 2013 年時点の郡特性が「介入以降のトレンド」を線形に決める |
| **(C) 時変 X** | `Y ~ D:post + X_t \| county_id + year` | 共変量を **時変** で線形に投入 | 共変量の時間変動がアウトカム水準を線形に決める |
**重要な違い**:
- **(B)** はベースライン値で固定。介入の結果として変化した X(媒介変数)を除外できるので
post-treatment bias を回避しやすい
- **(C)** は時変 X を入れるので post-treatment bias のリスクがある(X が介入で動くなら効果を吸収)
ステップ:
| ステップ | コード | 何をしているか |
|---|---|---|
| 1 | `covs_baseline <- did_data \|> filter(year == 2013) \|> select(...)` | 2013 年値を抽出してベースライン共変量に |
| 2 | `sub2 <- left_join(covs_baseline)` | 介入群×対照群サブサンプルにベースライン値を結合 |
| 3 | `feols(..., weights = ~pop_2013)` × 3 | (A), (B), (C) の 3 仕様で TWFE 推定 |
| 4 | `modelsummary(list(...))` | 3 列で並列表示 |
それでは R を走らせて比較してみます。
```{r}
#| label: twfe-with-cov
sub2 <- sub |> select(-D, -treated_now, -G_for_did)
# 2013 値で固定した共変量(時間で動かさない)
covs_baseline <- did_data |>
filter(year == 2013) |>
select(county_id,
baseline_white = pct_white,
baseline_poverty = poverty,
baseline_median = median_income
)
sub2 <- sub2 |>
left_join(covs_baseline, by = "county_id") |>
mutate(D = ifelse(G != Inf, 1L, 0L), post = (year == 2014) * 1L)
# 仕様 A:共変量なし(再掲)
mA <- feols(mortality_rate ~ D:post | county_id + year,
data = sub2,
weights = ~pop_2013, cluster = "county_id"
)
# 仕様 B:ベースライン X を post と交互作用(時間で動く共変量効果を許す)
mB <- feols(
mortality_rate ~ D:post +
baseline_white:post + baseline_poverty:post + baseline_median:post |
county_id + year,
data = sub2,
weights = ~pop_2013, cluster = "county_id"
)
# 仕様 C:時変共変量を線形で
mC <- feols(
mortality_rate ~ D:post + pct_white + poverty + median_income |
county_id + year,
data = sub2,
weights = ~pop_2013, cluster = "county_id"
)
modelsummary(
list("(A) X なし" = mA, "(B) X_2013 × Post" = mB, "(C) 時変 X" = mC),
coef_map = c(
"D:post" = "Medicaid × Post",
"baseline_white:post" = "% White (2013) × Post",
"baseline_poverty:post" = "Poverty (2013) × Post",
"baseline_median:post" = "Median Income (2013) × Post",
"pct_white" = "% White (時変)",
"poverty" = "Poverty (時変)",
"median_income" = "Median Income (時変)"
),
fmt = 2, stars = TRUE,
gof_omit = "AIC|BIC|RMSE|R2 Adj|R2 Within|Log",
output = "gt"
)
```
#### 出力の読み方
表は **3 列(A, B, C)**、行は推定された係数です:
- **`Medicaid × Post` 行**:これが DiD 推定量。**3 列で値がどう変わるか** が観察ポイント
- `% White (2013) × Post` などの行は (B) でのみ現れる時間相互作用の係数
- `% White (時変)` などの行は (C) でのみ現れる時変共変量の係数
- 括弧内は標準誤差、`*, **, ***` は有意水準(p < 0.1, 0.05, 0.01)
**確認すべきポイント**:
1. **3 仕様で `Medicaid × Post` の係数がどれくらい動くか** — 動きが大きいほど、
共変量調整の方法に推定値が敏感(=結果が頑健でないシグナル)
2. **(C) のように時変 X を入れた場合、効果が吸収されていないか** — 元の効果に対して
大きく値が縮むなら媒介変数を入れている可能性
3. **どの仕様も「正解」ではない** — TWFE は次節で見るように、共変量の組み込み方を含めて
ATT を一貫して返してくれる保証がない
#### 実際の結果の解釈
実値(2014 拡大群 vs 非拡大群、2013–2014、人口加重):
| 仕様 | `Medicaid × Post` | SE |
|---|---|---|
| (A) X なし | **−1.88** | 1.30 |
| (B) X_2013 × Post | **−1.24** | 1.33 |
| (C) 時変 X | **−1.81** | 1.27 |
**読み取り 1:仕様によって 0.64 動く**
(A) と (B) で **0.64**(=1.88−1.24)も動いています。SE が約 1.3 なので、
動きの大きさは「無視できる」とも言えませんが、結論の方向性(負)は変わりません。
ただし「**仕様の選択で推定値が動く**」という事実は、TWFE+共変量の信頼性を相対的に下げます。
**読み取り 2:(C) 時変 X は (A) X なしに近い**
(C) = −1.81 は (A) = −1.88 に近く、時変 X を入れても効果はほぼ変わっていません。
これは:
- 時変共変量(人種比率・貧困率・所得)がアウトカム水準を「コントロール」しているが、
介入効果そのものはあまり吸収されていない、ということ
- 別の見方:本シミュレーションでは時変 X と介入の相関が弱いので post-treatment bias が
顕著には出ていない
実データでは (C) で大きく値が動く場合、媒介変数を疑うべき。
**読み取り 3:どれが「正解」か?**
**どれも理論的には ATT を一貫して識別する保証がありません**。次節 §4.4 で導入する
**アウトカム回帰 / 逆確率ウェイト / 二重頑健** がより頑健で、TWFE+共変量の代替として
推奨されます。実際 §4.4 の二重頑健(人口加重)= **−2.37** は、ここの 3 仕様より少し大きめで、
これが「より信頼できる主結果」とされます。
#### 注意点
- **TWFE+共変量は仕様依存** — 異質効果があると ATT を返さない(Goodman-Bacon 2021)
- **(C) のように時変 X を入れるときは post-treatment bias に注意** — X が介入で動くなら
入れてはいけない(媒介変数)
- 実務では **(A) → (B) → 二重頑健 (§4.4)** の順で robustness を見せるのが標準
- 主結果は §4.4 の 二重頑健 で据え、TWFE+共変量はあくまで補助的な感度分析として位置付ける
→ **より頑健な方法** が次節の **アウトカム回帰 / 逆確率ウェイト / 二重頑健**(Outcome Regression / Inverse Probability Weighting / Doubly Robust)です。
## 4.4 ATT(2) を目標とする共変量付き DiD 推定量
3 つの主要な推定量があります。
::: {.callout-tip}
## 疫学研究との対応
ここで紹介する 3 つの推定量は、**疫学・因果推論文献の交絡調整法と完全に対応** しています。
DiD は「アウトカムの **差分** $\Delta Y$」に対して同じ手法を適用しているだけ、と理解すると見通しが良くなります。
| 本ノート(DiD) | 疫学研究での呼称 | やっていること |
|---|---|---|
| アウトカム回帰(Outcome Regression) | **G 公式(g-formula)/ G 計算(g-computation)/ 回帰による標準化(regression standardization)** | 結果モデル $E[\Delta Y \mid X, D=0]$ を当てて、介入群の X 分布で予測値を平均化 |
| 逆確率ウェイト(IPW) | **IPW**(Inverse Probability of Treatment Weighting)/ 周辺構造モデル(MSM)の重み | 傾向スコア $P(D \mid X)$ の逆数で観測を重み付け、X の分布を揃える |
| 二重頑健(Doubly Robust) | **DR / AIPW**(Augmented IPW)/ TMLE 等 | 上 2 つを組み合わせ、片方のモデル誤特定にロバスト |
::: {.callout-warning}
## 「標準化」という言葉の使い方に注意
疫学で「標準化(standardization)」というと、伝統的には **直接標準化**(層別の率を
標準集団のサイズで重み付け)や **間接標準化(SMR)** を指すノンパラメトリックな手法です。
アウトカム回帰はこれらと **発想は同じ**(対照群の関係を介入群の X 分布で平均する)ですが、
**回帰モデルを介する** 点で別物です。より正確には **G 公式 / G 計算 / 回帰標準化
(model-based / regression standardization)** と呼びます(Hernán & Robins
"Causal Inference: What If" 13 章を参照)。
:::
つまり、本節は **疫学で扱う交絡調整(Robins, Hernán らの potential outcome 系の手法)を
「アウトカムの差分 $\Delta Y$」に適用しているだけ** ですので、
疫学のバックグラウンドがある方は知っている手法と思って読み進めていただいて構いません。
**唯一の違いは何に「適用するか」**:
- 疫学の交絡調整: $Y$(レベル)に適用
- DiD の調整推定量: $\Delta Y = Y_t - Y_{t-1}$(差分)に適用
これにより、X で観察できる交絡だけでなく、$Y(0)$ の **共通トレンド** という形で
**観察不能な時不変交絡** にも対処できる点が DiD の利点です。
:::
3 つの定義は次のとおり:
::: {.callout-note}
## 用語:アウトカム回帰(Outcome Regression)
対照群でアウトカム変化を共変量に回帰し、介入群の共変量で予測。
**結果モデルが正しく特定されている** ことが必要。
:::
::: {.callout-note}
## 用語:傾向スコア(Propensity Score)
「共変量 X を持つ個体が介入を受ける条件付き確率」 $P(D=1 \mid X)$ のこと。
通常はロジスティック回帰で推定:`glm(D ~ X, family = binomial)`。
:::
::: {.callout-note}
## 用語:逆確率ウェイト(Inverse Probability Weighting)
各観測を傾向スコアの逆数で重み付けし、介入群と対照群の X 分布を揃える。
**傾向スコアモデルが正しく特定されている** ことが必要。
:::
::: {.callout-note}
## 用語:二重頑健(Doubly Robust)
アウトカム回帰と逆確率ウェイトを組み合わせ、**いずれかのモデル**
(結果モデル or 傾向スコアモデル)が正しく特定されていれば
一致推定を与える。実務的に最も頑健で推奨される。
:::
#### この解析で何を示すか
ここまでの議論を踏まえ、メディケイドの 2×2 設定で次を実装します:
1. **共変量(2013 年の白人率・貧困率・所得)を 3 種類の方法で調整した DiD** を推定する
- アウトカム回帰(Outcome Regression / 回帰標準化)
- 逆確率ウェイト(Inverse Probability Weighting)
- 二重頑健(Doubly Robust)
2. **同じデータ・同じ共変量・同じ目標(ATT)** に対して 3 手法が **似た値を返すか**(モデル誤特定の影響に頑健か)を確認
3. さらに **二重頑健 × 人口加重** の組み合わせも 1 行追加し、§3 の無加重/加重対比とつなげる
R では `DRDID` パッケージの `ordid()` / `ipwdid()` / `drdid()` を順に呼び出します。
#### 表の見方
出力は **4 行 × 3 列** の表です:
| 列 | 意味 |
|---|---|
| `推定量` | 手法名(アウトカム回帰 / 逆確率ウェイト / 二重頑健 / 二重頑健(人口加重)) |
| `ATT` | 各手法の ATT(2014) 推定値 |
| `SE` | クラスタロバスト標準誤差 |
確認すべきポイント:
- **3 手法の ATT 推定値がほぼ揃うか?** → 揃えば「モデル誤特定の影響に頑健」というシグナル
- **二重頑健と他の 2 つの差** → モデル誤特定があれば二重頑健の値が他と乖離するはず
(本ノートはシミュレーションなので、両モデルとも近似的に正しく、3 手法はかなり一致するはず)
- **二重頑健 vs 二重頑健(人口加重)** → §3 で見た「無加重 ≈ 0、加重 ≈ −2」のパターンが
共変量調整後も残るか
- **SE の大きさ** → 効果がゼロから区別できるか(95% CI ≈ ATT ± 1.96·SE)
#### 結果の解釈
- **3 手法の推定値がほぼ等しい** ことを期待します(シミュレーションでは結果モデル・
傾向スコアモデルとも線形でだいたい正しく、二重頑健性が温存される設計)
- **二重頑健(無加重)が ≈ −0.2** 程度なら、§3 と整合(郡単位の集計平均では効果はゼロ近傍)
- **二重頑健(人口加重)が ≈ −1.9 〜 −2.6** なら、§3 と整合
(人口加重では有意な死亡率低下)
- 共変量で条件付けても全体の結論が大きく変わらないなら、観察可能な交絡は
**平行トレンドを致命的に破壊するほどではなかった** ことを示唆します
**主結果としては 二重頑健 推定量を採用するのが現代の実務的標準** です:
アウトカム回帰と逆確率ウェイトの **両方の性質を同時に持ち**(結果モデルか
傾向スコアモデルかの **どちらか片方さえ正しく特定されていれば一致推定** を与える)、
モデル誤特定のリスクに対して最も頑健だからです。
3 手法を併記するのは、二重頑健 単体ではなく **手法間の頑健性も合わせて確認するため** です。
```{r}
#| label: drdid-2x2
# DRDID パッケージで 3 種類を実行
# 2 期間 2 群のサブサンプルを準備
# 重要: DRDID は「時間不変な共変量」を要求するので、2013 年の値で固定する
covs_baseline_2013 <- did_data |>
filter(year == 2013) |>
select(county_id,
X_pct_white = pct_white,
X_poverty = poverty,
X_median_inc = median_income
)
sub3 <- did_data |>
filter(year %in% c(2013, 2014), G == 2014 | G == Inf) |>
mutate(D = ifelse(G != Inf, 1L, 0L)) |>
select(-pct_white, -poverty, -median_income) |>
left_join(covs_baseline_2013, by = "county_id") |>
rename(pct_white = X_pct_white, poverty = X_poverty, median_income = X_median_inc) |>
arrange(county_id, year)
# DRDID のユーザ向けラッパー関数: yname, tname, idname, dname, xformla, data, panel=TRUE
res_ra <- DRDID::ordid(
yname = "mortality_rate", tname = "year",
idname = "county_id", dname = "D",
xformla = ~ pct_white + poverty + median_income,
data = sub3, panel = TRUE
)
res_ipw <- DRDID::ipwdid(
yname = "mortality_rate", tname = "year",
idname = "county_id", dname = "D",
xformla = ~ pct_white + poverty + median_income,
data = sub3, panel = TRUE
)
res_dr <- DRDID::drdid(
yname = "mortality_rate", tname = "year",
idname = "county_id", dname = "D",
xformla = ~ pct_white + poverty + median_income,
data = sub3, panel = TRUE
)
# 加重版も(weightsname で人口加重)
res_dr_w <- DRDID::drdid(
yname = "mortality_rate", tname = "year",
idname = "county_id", dname = "D",
xformla = ~ pct_white + poverty + median_income,
weightsname = "pop_2013",
data = sub3, panel = TRUE
)
tibble(
推定量 = c("アウトカム回帰", "逆確率ウェイト", "二重頑健", "二重頑健(人口加重)"),
ATT = c(res_ra$ATT, res_ipw$ATT, res_dr$ATT, res_dr_w$ATT),
SE = c(res_ra$se, res_ipw$se, res_dr$se, res_dr_w$se)
) |>
gt() |>
fmt_number(decimals = 2)
```
#### 実際の結果の解釈
得られた値は次のとおり:
| 推定量 | ATT | SE | 95% CI(≈ ATT ± 1.96·SE) | 結論 |
|---|---|---|---|---|
| アウトカム回帰 | **−0.35** | 0.44 | (−1.21, +0.51) | **0 を含む** → 有意でない |
| 逆確率ウェイト | **−0.02** | 0.50 | (−1.00, +0.96) | **0 を含む** → 有意でない |
| 二重頑健 | **−0.10** | 0.46 | (−1.00, +0.80) | **0 を含む** → 有意でない |
| 二重頑健(人口加重) | **−2.37** | 1.33 | (−4.98, +0.24) | **辛うじて 0 を含む** → p ≈ 0.07 程度 |
**読み取り 1:無加重 3 手法はほぼ揃っている(モデル誤特定なし)**
無加重の 3 推定量はすべて **ゼロ近傍**(−0.35, −0.02, −0.10)で、互いに大きな矛盾はありません。
これは **結果モデル(線形 OLS)・傾向スコアモデル(ロジット)が両方とも妥当に
当てはまっている** = モデル誤特定の懸念が小さいことを意味します。
特に二重頑健(−0.10)が他の 2 つの中間に位置しているのが「正常な健康状態」のシグナル。
**読み取り 2:無加重ではゼロ、人口加重では −2.37**
§3 と同じパターンが共変量調整後も維持されています:
- **無加重(3 手法平均):≈ −0.15** → 郡集計平均では効果なし
- **人口加重(二重頑健):−2.37**(95% CI が辛うじて 0 をまたぐ程度に有意性に近い)
→ 人口加重では明らかに死亡率低下方向
#### §3 との比較
| 推定量 | §3(共変量なし) | §4.4(共変量付き) | 結論への影響 |
|---|---|---|---|
| 無加重 | −0.23 | −0.02 〜 −0.35 | 共変量で条件付けても **郡集計平均** ではゼロ近傍(変わらず) |
| 人口加重 | −1.88 | −2.37 | **人口加重** では一貫して死亡率低下(むしろ少し大きく) |
つまり共変量(2013 年の白人率・貧困率・所得)を入れても **結論は質的に変わりません**。
これは:
1. **シミュレーション設計上、無条件 PT が既に成り立っている** ため
(DGP には人種・経済共変量と $Y(0)$ トレンドの関連を入れていない)
2. 現実のメディケイドデータでも同様のパターンが報告されています
(Borgschulte 2020)
**もし共変量調整で大きく値が動いたら**、それは無条件 PT が破れている兆候。
本ノートでは動かないので、観察可能な交絡で説明できる平行トレンド違反は
深刻ではない、と判断できます。
## 4.5 異質性分析
### このセクションで何をするか
ここまでは「**介入群全体での平均** 介入効果(ATT)」を 1 つの値として推定してきました。
しかし実務的には次の問いも重要です:
- **誰に対して効果が大きいか?**(政策の便益が集中している層は?)
- **特定のサブグループでだけ効果が出ていないか?**(取り残されている層は?)
- **連続変数(貧困率など)に沿って効果が滑らかに変化するか?**
これを定量化するのが **異質性分析(heterogeneous treatment effects 分析)** です。
### こんなときに必要
実務で異質性を見るべき典型例:
1. **政策のターゲティング** — 限られた予算で誰に介入すべきか決めたい
(例:メディケイドを **貧困率が高い郡** に拡大する意義はより大きい?)
2. **副作用や負の効果の検出** — 特定の層で逆効果になっていないかを確認
3. **メカニズムの示唆** — 効果が大きい層の特徴から、効果が生じる経路を推測
4. **論文の査読対応** — 「全体平均だけでは粗い」というコメントへの応答
### コードで何をしているか
ここでは **貧困率の四分位** で郡を 4 グループに分け、各グループで **二重頑健 DiD** を別々に推定します。
| ステップ | コード | 何をしているか |
|---|---|---|
| 1 | `mutate(pov_q = ntile(poverty, 4))` | 2013 年の貧困率で郡を 4 等分(Q1 = 最低貧困率、Q4 = 最高) |
| 2 | `map_dfr(1:4, function(q) { ... })` | 各四分位 q について以下を実行 |
| 3 | `filter(pov_q == q)` | その四分位だけの郡データを抽出 |
| 4 | `DRDID::drdid(...)` | そのサブサンプルで二重頑健 DiD を推定 |
| 5 | `tibble(quartile, ATT, SE)` | 結果を行として収集 |
| 6 | `ggplot(...)` | 四分位 × ATT で点 + 95% CI バーを描画 |
注意点: ここでは共変量を `~ pct_white + median_income` に絞っています
(`poverty` で層別化したので、層内では vergleich 比較的均一になるため)。
### 結果の読み方と解釈
#### グラフの構造
- **横軸**: 貧困率四分位(1 = 最低、4 = 最高)
- **縦軸**: 各四分位での ATT 推定値(10 万人あたり死亡率変化)
- **点**: 点推定値、**縦バー**: 95% 信頼区間(ATT ± 1.96·SE)
- **横破線(0)**: 効果ゼロのライン
#### 読み取りポイント
1. **四分位ごとに ATT が違うか?**
- 4 点がほぼ同じ高さに並ぶなら → **異質性は小さい**(全体平均 ATT で十分)
- 階段状に変化するなら → **貧困率に沿った異質性あり**
2. **どの四分位が 0 から離れているか?**
- CI バーが 0 を跨がない四分位 = その層で統計的に有意な効果
3. **本ノートの DGP では?**
- シミュレーションでは「貧困率 × 効果」の構造は入れておらず、
「人口 × 効果」だけ入れているので、**貧困率での明確な異質性は出ない** はずです
- 現実のメディケイド研究では、**貧困率が高い郡で効果が大きい** という結果が
報告されています(政策の便益が必要な層により届くため)
#### 注意:多重比較とサンプルサイズ
- 四分位ごとに分けるとサブサンプルが小さくなり、**SE が大きく** なります
- 4 つの推定を行うので **多重比較の問題** も生じます(個別の有意性は控えめに解釈)
- 連続変数で異質性を見る場合は **CATE 推定(causal forest 等)** が選択肢ですが本ノートでは扱いません
```{r}
#| label: heterogeneity
#| fig-cap: "貧困率四分位ごとの 二重頑健 DiD 推定値"
#| fig-height: 4
sub3 <- sub3 |> mutate(pov_q = ntile(poverty, 4))
hetero <- map_dfr(1:4, function(q) {
d <- sub3 |> filter(pov_q == q)
if (length(unique(d$D)) < 2) {
return(NULL)
}
fit <- tryCatch(
DRDID::drdid(
yname = "mortality_rate", tname = "year",
idname = "county_id", dname = "D",
xformla = ~ pct_white + median_income, data = d,
panel = TRUE
),
error = function(e) NULL
)
if (is.null(fit)) {
return(NULL)
}
tibble(quartile = q, ATT = fit$ATT, SE = fit$se)
})
ggplot(hetero, aes(factor(quartile), ATT)) +
geom_pointrange(aes(ymin = ATT - 1.96 * SE, ymax = ATT + 1.96 * SE),
size = 0.7, color = "steelblue"
) +
geom_hline(yintercept = 0, linetype = "dashed", color = "grey50") +
labs(
x = "貧困率四分位(1 = 最低, 4 = 最高)",
y = "ATT(共変量付き 二重頑健)",
title = "貧困率が高い郡ほど大きな ATT? 異質性の検査"
)
# 数値も表で確認
hetero |>
mutate(
CI_lo = ATT - 1.96 * SE,
CI_hi = ATT + 1.96 * SE
) |>
gt::gt() |>
gt::fmt_number(decimals = 2) |>
gt::cols_label(
quartile = "貧困率四分位",
ATT = "ATT",
SE = "SE",
CI_lo = "95% CI 下限",
CI_hi = "95% CI 上限"
)
```
#### 実際の異質性結果の解釈
得られた値を整理すると:
| 貧困率四分位 | ATT | SE | 95% CI |
|---|---|---|---|
| Q1(最低) | **−0.19** | 0.86 | (−1.89, +1.50) — **0 を含む** |
| Q2 | **+0.91** | 0.93 | (−0.91, +2.73) — **0 を含む** |
| Q3 | **−0.63** | 0.86 | (−2.32, +1.06) — **0 を含む** |
| Q4(最高) | **−1.53** | 0.82 | (−3.14, +0.08) — **辛うじて 0 を含む**(ほぼ有意) |
**読み取り 1:4 点に「貧困率が高いほど効果が大きい」傾向が見える?**
Q1 → Q4 の点推定値を順に見ると **−0.19 → +0.91 → −0.63 → −1.53** で、
**Q2 だけ正に外れているもののざっくりは右下がり(貧困率が高い四分位ほど ATT がより負)**
というパターンです:
- Q1(貧困率最低)では効果ほぼゼロ
- Q4(貧困率最高)では ATT = −1.53、95% CI 上限が +0.08 で **辛うじて 0 を含む程度に
有意性に近い**
- Q2 の異常値は **偶然のばらつき**(SE が 0.93 と大きい)と考えるのが妥当
**読み取り 2:本ノートの DGP では本来出るべきでない傾向**
ここが教育的に重要です。本ノートの DGP は **貧困率と効果の関連を入れていません**
(真の効果は人口と相関するように設定)。にもかかわらず Q1→Q4 で傾きが見えるのは:
1. **貧困率と人口の間接的相関** — シミュレーションでコホート(2014 拡大群 / 非拡大群)ごとに
貧困率分布と人口分布の両方をずらしているため、貧困率と人口に弱い相関が生まれる
2. **小サンプル × 4 区間** によるサンプリング・ノイズ — 各四分位は約 550 郡で、
SE が 0.8〜0.9 と大きく、点推定はかなり揺れる
3. **多重比較** — 4 つの推定を行うので、偶然有意/非有意のパターンが見える
**この経験は実務でとても重要な教訓**:
**異質性は簡単に「見える」が、それが本物のメカニズムを反映している保証はない**。
本ノートではメカニズムが「無い」ことが分かっているので、見えた傾向は **偽** と判断できます。
**読み取り 3:CI バーから読める情報**
- すべての四分位で 95% CI が 0 をまたぐ → 個別では統計的に有意でない
- Q4 は CI 上限が +0.08 で **「ほぼ 0」** に到達 → サンプルがもう少し多ければ有意になる可能性
- CI 幅が四分位全体推定(§4.4)より広いのは、サブサンプルが約 1/4 になったため
**読み取り 4:現実のメディケイド研究との対比**
現実のデータでは **貧困率が高い郡で効果が大きい** という結果が報告されています
(Borgschulte 2020 等)。本ノートの結果(Q4 で最大の負の効果)もこの方向性と一致しますが、
**今回のは「シミュレーションのノイズ」として偶然出てきた可能性が高い** ことを忘れず:
- 現実データで同じパターンが出る → 政策のターゲティング示唆として意味あり
- 本ノートで同じパターンが出る → **偽陽性の見本** として学ぶべき
この対比こそが、異質性分析を **慎重に解釈する** べき理由を示しています。
#### 異質性分析を実務で報告するときの注意
- **事前登録の重要性**: サブグループを後付けで選んで「ここで有意」と報告するのは p-hacking。
解析計画(SAP)で事前にサブグループを定めるのが望ましい
- **多重比較補正**: 4 つの推定で個別検定するなら、Bonferroni 補正等を考慮
- **CATE / 因果フォレスト**: 連続変数や複数共変量での異質性は機械学習ベースの
推定(`grf` パッケージの causal forest 等)が選択肢
- **解釈は控えめに**: 異質性は仮説生成的(hypothesis-generating)に使い、
主結果は全体 ATT を据えるのが標準的
---
# 複数期間を持つ DiD デザイン(論文 §5)
複数期間に拡張すると、**プリトレンド検定** と **動学効果** が見られるようになります。
## 5.1 単純なイベント・スタディ($2 \times T$)
#### このセクションで何をするか
§3 では **2 期(2013, 2014)** だけを使った 2×2 DiD を見ました。ここでは
**1 つの介入群(2014 拡大)vs 1 つの対照群(非拡大)** を保ったまま、
**期間を 2009 〜 2019 の 11 年に拡張** します。これにより 1 つの推定値ではなく
「年ごとの DiD 推定値の系列」が得られ、
1. **介入前期間** での DiD ≈ 0 かを見て、**平行トレンド仮定の経験的検証**
2. **介入後期間** での DiD の **時間的な変化**(効果の遅れ・継続性・減衰)
を同時に評価できます。これが **イベント・スタディ(event study)** です。
#### こんなときに必要
実務での典型例:
- **介入の遅延効果** を確認したい(例:政策の効果は導入後 2 年経ってから現れる)
- **平行トレンドの妥当性** をプロットで視覚的に説明・査読対応したい
- **効果が時間とともに減衰するか** を見たい(例:キャンペーン効果は薄れる)
- **介入時期付近のショック** (例:政策発表 vs 実施の前倒し効果)を確認したい
::: {.callout-note}
## 用語:イベント・スタディ(Event Study)
「**介入開始からの相対時間**」(イベント時間 $e = t - g$)を横軸に取り、
各時点での DiD を縦軸にプロットしたもの。
ここで:
- $t$ は **暦時点(calendar time)**(本ノートでは年:2009, 2010, ..., 2019)
- $g$ は **介入開始時点(group / cohort)**(本ノートでは介入群の介入年:2014, 2015, 2016, 2019。
非拡大群は $g = \infty$)
- $e = t - g$ は **イベント時間**:介入開始から何年経ったか
- $e = 0$: 介入開始の年
- $e = -1$: 介入の 1 年前
- $e = +2$: 介入から 2 年後
例:2014 拡大群($g = 2014$)を 2016 年($t = 2016$)に見れば、$e = 2016 - 2014 = 2$。
介入前(負のイベント時間 $e < 0$)は **プリトレンド検定**、
介入後(正のイベント時間 $e \geq 0$)は **動学効果** に対応。
:::
::: {.callout-note}
## 用語:プリトレンド(Pre-trend)
介入前の期間における **介入群 vs 対照群の差の動き**。
平行トレンドが成り立っていれば 0 付近で水平に並ぶはず。「プリトレンド検査」は平行トレンドの **間接的** な検定。
:::
#### イベント・スタディが計算しているもの(2×2 DiD の系列)
イベント・スタディの 1 つの点(年 $t$ での推定値)は、実は **§3 で見た 2×2 DiD を、
基準年 vs その年の組み合わせで計算したもの** です。
数式で書くと、基準年 $t_0 = g - 1$(本ノートでは $g = 2014$ なので $t_0 = 2013$)
を固定し、各年 $t$ について:
$$
\widehat{\beta}_t = \big( \bar{Y}_{D=1, t} - \bar{Y}_{D=1, t_0} \big)
- \big( \bar{Y}_{D=0, t} - \bar{Y}_{D=0, t_0} \big)
$$
つまり、
- 横軸の **各点が独立した 2×2 DiD**
- 2 つの群:介入群 vs 対照群(共通)
- 2 つの時点:**基準年(2013)** と **その年 $t$**
具体例(構造):
| イベント時間 $e$ | 暦年 $t$ | 計算する 2×2 DiD |
|---|---|---|
| $e = -3$ | 2011 | (2011 拡大群 − 2011 非拡大群) − (2013 拡大群 − 2013 非拡大群) |
| $e = -1$ | 2013 | 基準年なので **必ず 0**(自己差分はゼロ) |
| $e = 0$ | 2014 | (2014 拡大群 − 2014 非拡大群) − (2013 拡大群 − 2013 非拡大群) ← **§3 の 2×2 DiD と同じもの** |
| $e = +2$ | 2016 | (2016 拡大群 − 2016 非拡大群) − (2013 拡大群 − 2013 非拡大群) |
**実際の数値例 — $e = 0$ の点(2014 年)を手計算で**
§3 の Table 2(無加重)の値を使うと:
- 2013 拡大群の平均死亡率 ≈ **407**(人口加重なら 381)
- 2013 非拡大群の平均死亡率 ≈ **474**(人口加重なら 457)
- 2014 拡大群の平均死亡率 ≈ **409**(人口加重なら 380)
- 2014 非拡大群の平均死亡率 ≈ **476**(人口加重なら 455)
無加重 2×2 DiD($e = 0$):
$$
\widehat{\beta}_{2014} = (409 - 407) - (476 - 474) = 2 - 2 = 0 \;\;\text{(≈ −0.2、§3 と一致)}
$$
人口加重 2×2 DiD($e = 0$):
$$
\widehat{\beta}_{2014}^{wt} = (380 - 381) - (455 - 457) = -1 - (-2) = +1 \;\;\text{(≈ −1.9、§3 と一致)}
$$
このように **イベント・スタディの $e = 0$ の点は、§3 で行った 2×2 DiD と完全に同じ計算**
です。同じ論理を $e = -3, -2, ..., +5$ まで全ての年について繰り返し、その結果を
1 枚のグラフに並べたのがイベント・スタディです。
**$e = +1$ の点(2015 年)の場合**
シミュレーションで実際に得られる値:
- 介入群の 2013→2015 変化:約 +20(年トレンド +9 × 2 年 + 真の効果 ≈ −2 + ノイズ)
- 対照群の 2013→2015 変化:約 +20(年トレンド +9 × 2 年 + ノイズ)
- $\widehat{\beta}_{2015}$ ≈ −2(両群の長期変化の差 = 真の効果)
このように $e = 0$ で見えにくかった効果が、**$e = 1$ で明確に負として現れる**
(介入の効果が 1 年かけて発現するパターン)— という動学的なメッセージが
イベント・スタディの主役です。
「介入から各イベント時間までの差の差」というイメージはほぼ正しく、より正確には
**「基準年(介入直前)から各年への、両群アウトカム変化の差」** を計算しています。
回帰式 `feols(Y ~ i(year, D, ref = 2013) | county_id + year)` は、これら全ての
$\widehat{\beta}_t$ を 1 本の回帰で同時に推定する仕組みです(年×D の交互作用項を年ごとに一括生成)。
つまりイベント・スタディは **§3 の 2×2 DiD を 11 個並列に走らせ、1 枚のグラフに並べたもの**
と理解できます。
#### コードで何をしているか
| ステップ | コード | 何をしているか |
|---|---|---|
| 1 | `filter(G == 2014 \| G == Inf)` | 2014 拡大群と非拡大群だけのサブサンプルを作る |
| 2 | `mutate(D = ifelse(G == 2014, 1, 0))` | 介入群ダミー(1 if 2014 拡大、0 if 非拡大) |
| 3 | `feols(Y ~ i(year, D, ref = 2013) \| county_id + year, weights = ~pop_2013)` | 各年で `D` との交互作用を推定。`ref = 2013` で 2013 年を基準にゼロ |
| 4 | `broom::tidy(es_fit, conf.int = TRUE)` | 推定値と 95% CI を tidy 形式に |
| 5 | `bind_rows(..., year = 2013, estimate = 0)` | 基準年(2013)を 0 として追加 |
| 6 | `ggplot(es_coefs, ...)` | 年×ATT で点+CI バーを描画、介入境界に縦破線 |
`fixest::i()` の役割:`year` 変数の **各水準** と `D` の交互作用を一括生成。
`ref = 2013` で 2013 を基準カテゴリにします(その年の係数 = 0)。
```{r}
#| label: fig-event-study-2xt
#| fig-cap: "2×T イベント・スタディ:2014 拡大群 vs 非拡大群(共変量なし、人口加重)"
es_data <- did_data |>
filter(G == 2014 | G == Inf) |>
mutate(
D = ifelse(G == 2014, 1L, 0L),
event_time = ifelse(D == 1, year - 2014, -1000)
) # 対照群は基準群
# i(year, D, ref = 2013): 2013 を基準に各年で D との交互作用を推定
es_fit <- feols(mortality_rate ~ i(year, D, ref = 2013) | county_id + year,
data = es_data, weights = ~pop_2013, cluster = "county_id"
)
es_coefs <- broom::tidy(es_fit, conf.int = TRUE) |>
mutate(year = as.numeric(gsub("year::|:D", "", term))) |>
filter(!is.na(year))
# 基準年 2013 をゼロとして追加
es_coefs <- bind_rows(es_coefs, tibble(
year = 2013, estimate = 0,
conf.low = 0, conf.high = 0
))
# 数値表も出力
es_coefs |>
arrange(year) |>
select(year, estimate, std.error, conf.low, conf.high) |>
gt::gt() |>
gt::fmt_number(decimals = 2) |>
gt::tab_header(title = "イベント・スタディの年別 ATT 推定値")
ggplot(es_coefs, aes(year, estimate)) +
geom_hline(yintercept = 0, linetype = "dashed", color = "grey50") +
geom_vline(xintercept = 2013.5, linetype = "dashed", color = "firebrick") +
geom_pointrange(aes(ymin = conf.low, ymax = conf.high), color = "steelblue") +
scale_x_continuous(breaks = 2009:2019) +
labs(
x = "年", y = "推定 ATT(10万人あたり)",
title = "イベント・スタディ:2014 群 vs 非拡大群",
subtitle = "破線(2013.5)より前 = プリトレンド検査、後 = 動学効果"
)
```
#### 出力の読み方
**グラフの構造**:
- **横軸**: 暦年(2009–2019)
- **縦軸**: 各年における (2014 拡大群 vs 非拡大群) の DiD 推定値(2013 年基準)
- **赤い縦破線(2013.5)**: 介入開始境界。左 = プリトレンド検査、右 = 動学効果
- **横破線(0)**: 効果ゼロのライン
- **青い点 + 縦バー**: 点推定値と 95% 信頼区間
**チェックポイント**:
1. **介入前(2009〜2013)**: 点が 0 付近で水平に並んでいるか?
- 水平 → 平行トレンド仮定が支持される(経験的検証クリア)
- 系統的な傾き(右上がり/右下がり) → 平行トレンド違反のリスク
2. **2013(基準年)**: 必ず 0(点と CI バーが収束)
3. **介入後(2014〜2019)**:
- 効果は **即座に現れるか、遅延するか**
- 時間とともに **強まる/減衰する/持続する** どのパターンか
- 95% CI が 0 をまたぐか = 統計的有意性
4. **大局のメッセージ**: 介入後の点が一貫して負側にあれば、メディケイドが死亡率を下げた示唆
#### 実際の結果の解釈
シミュレーションでは(§5.2 の動学効果と整合する範囲で):
- **プリトレンド(2009–2013)**: ATT 推定値は ±1 程度の幅でほぼ 0 周辺を推移し、
CI バーが大きく 0 をまたぐ → 平行トレンド仮定は **支持される**
- **2014(介入直後)**: 効果はわずかに負だが CI は広く、まだ統計的有意でない
- **2015 以降**: 点推定値が **約 −2 〜 −4** で安定し、特に 2015–2016 年で CI 上限が
0 に近い(あるいは下回る) → 統計的に意味のある **死亡率低下**
- **全体パターン**: メディケイド拡大が **介入後すぐではなく 1 年遅れで** 効果を発揮し、
その後継続する形
これは §5.2 の集約イベント・スタディ(Fig 8: $e = 1$ で −3.92 *)とも整合します。
#### 注意点
- このセクションは **2014 拡大群 vs 非拡大群だけ** を使っており、§5.2 のような全コホート
集約ではありません。他の拡大コホートを除外している点に注意
- イベント時間がコホート間で違うので、多群を 1 つのイベント・スタディに集約するのは
§5.2 で扱う(Callaway-Sant'Anna)
- `fixest::i(year, D, ref = 2013)` のような **ナイーブな OLS イベント・スタディ** は
多群設定で問題を起こす(Sun-Abraham 2021)。多群では `did::att_gt` 等を使うのが安全
## 5.2 段階的介入採択($G_\# \times T$)
#### このセクションで何をするか
§5.1 では「介入群 1 つ vs 対照群 1 つ」で時間方向に拡張しましたが、
ここでは **介入開始時期が単位ごとに異なる(staggered)** 状況を扱います。
本ノートの例では、メディケイドを 2014 年、2015 年、2016 年、2019 年に拡大した州や
非拡大州が混在し、**4 つの介入コホート + 非拡大群** の構造です。
そこで:
1. **コホートごとに** ATT を時点別に推定(`did::att_gt()`)
2. **イベント時間** で集約(`did::aggte()`)
3. **他の推定量(Sun-Abraham, Borusyak, Roth-Sant'Anna)と比較** して頑健性を確認
#### こんなときに必要
実務での典型例:
- 政策が **州ごとに違うタイミング** で導入される(本ノートのメディケイド例)
- 介入が **企業ごとに段階的に展開** される(IT 導入、規制適用、認証取得)
- 医療政策やプログラムが **病院・地域ごとに異なる時期** に開始される
- 古典的な TWFE 回帰では **誤った推定値** を返すリスク(§5.3 で詳説)
- → 多群・多時期で安全に ATT を推定したいときに必須
::: {.callout-note}
## 用語:ATT(g, t) — 群・時点別 ATT
「**コホート $g$ で開始した介入群が、時刻 $t$ で受ける ATT**」。段階的設定の基本パラメータ。
全てのコホート $g$ と時点 $t$ の組み合わせで推定可能。
:::
::: {.callout-note}
## 用語:平行トレンドの 3 種類
論文では識別仮定として 3 種類の平行トレンドを区別:
- **PT-GT-Nev**: 比較群は **never-treated のみ**(最も保守的)
- **PT-GT-NYT**: 比較群は **まだ介入されていない群**(中間)
- **PT-GT-all**: **全ての群・全ての時点** で平行トレンド(最も強い)
`did::att_gt()` の引数 `control_group` で選べる:
`"nevertreated"`(Nev) / `"notyettreated"`(NYT)。
:::
::: {.callout-note}
## 用語:Callaway-Sant'Anna 推定量
段階的 DiD のための **群・時点別 ATT(g,t) の集合** を直接推定し、
その後 **解釈可能なウェイト** で集約する方法(Callaway & Sant'Anna 2021)。
R 実装は `did` パッケージ。
:::
#### コードで何をしているか
§5.2 全体の処理フローは次のとおり。Fig 5 〜 Fig 9 と「他推定量との比較」を順に実行:
| ステップ | コード/関数 | 何をしているか |
|---|---|---|
| 1. 可視化 (Fig 5) | `ggplot(... mortality_rate by cohort ...)` | コホート別の生死亡率トレンドを描き、段階的構造を視覚化 |
| 2. ATT(g,t) 推定 | `did::att_gt(yname, tname, idname, gname, ...)` | 4 コホート × 11 年の全 ATT(g,t) を二重頑健で推定 |
| 3. ATT(g,t) 表 | `summary(att_no_cov)` | 数値で 4 コホートの ATT(g,t) と SE を確認 + プリトレンド検定 |
| 4. 可視化 (Fig 6) | `did::ggdid(att_no_cov, ncol=2)` | コホート別の ATT(g,t) を 4 パネルで描画 |
| 5. 動学集約 (Fig 8) | `did::aggte(att_no_cov, type="dynamic")` | 全コホートを **イベント時間** で集約した 1 本のイベント・スタディに |
| 6. 共変量付き集約 (Fig 9) | `did::att_gt(... xformla=~X)` + `aggte()` | 共変量で条件付けた版(二重頑健 × 集約) |
| 7. 他推定量比較 | `fixest::sunab()`, `didimputation::did_imputation()`, `staggered::staggered()` | 同じ問いに対する 4 つの推定量を並列に走らせ頑健性を確認 |
引数の意味:
- `gname = "G_for_did"`: コホート列(`Inf` を `0` に変換済 = never-treated)
- `control_group = "notyettreated"`: 比較群は「まだ介入されていない群」(PT-GT-NYT 仮定)
- `weightsname = "pop_2013"`: 人口加重
- `bstrap = TRUE, cband = TRUE`: ブートストラップで同時信頼区間を計算
- `est_method = "dr"`: 二重頑健推定(`xformla` で共変量を指定するとき)
### コホート別のトレンド可視化(Fig 5 再現)
```{r}
#| label: fig-5-cohort-trends
#| fig-cap: "コホート別の人口加重 死亡率トレンド(Fig 5 再現)"
#| fig-height: 4.5
trends <- did_data |>
group_by(cohort_str, year) |>
summarise(Y = weighted.mean(mortality_rate, pop_2013), .groups = "drop") |>
mutate(cohort_str = factor(cohort_str,
levels = c("2014", "2015", "2016", "2019", "Inf"),
labels = c(
"G = 2014", "G = 2015", "G = 2016",
"G = 2019", "Never"
)
))
ggplot(trends, aes(year, Y, color = cohort_str)) +
geom_line(linewidth = 1) +
geom_point() +
scale_x_continuous(breaks = 2009:2019) +
scale_color_brewer(palette = "Set1") +
labs(
x = "年", y = "人口加重 死亡率(/10万)",
color = "コホート",
title = "拡大コホート別の死亡率トレンド"
)
```
**Fig 5 の読み方**: 5 本の線(G=2014, 2015, 2016, 2019, Never)の **垂直線が
コホートの介入開始時期**。介入後にコホート別の挙動が分かれるかどうかが視覚的ヒントになります。
本シミュレーションでは介入後、特に 2014・2015 コホートで死亡率が他より低めに推移していくのが見えるはずです(人口加重で大きい郡の効果が反映されているため)。
### Callaway-Sant'Anna による ATT(g,t) 推定(共変量なし)
```{r}
#| label: cs-att-gt-no-cov
#| message: false
# `did::att_gt`: G_for_did で Inf を 0 (=never-treated) に変換済み
att_no_cov <- did::att_gt(
yname = "mortality_rate",
tname = "year",
idname = "county_id",
gname = "G_for_did", # コホート(0 = never-treated)
data = did_data,
weightsname = "pop_2013",
control_group = "notyettreated",
panel = TRUE,
bstrap = TRUE,
cband = TRUE,
base_period = "varying"
)
summary(att_no_cov)
```
**`summary()` 出力の読み方**: 行は `(Group, Time)` の組、列は ATT(g,t) 推定値・SE・
同時信頼区間。`*` は CI が 0 を含まないことを示します。最終行に
**プリトレンドの事前検定 p 値**(`P-value for pre-test of parallel trends assumption`)
が表示され、p < 0.05 なら平行トレンド仮定が **棄却** されます。
本シミュレーションでは **p ≈ 0.019**(低い)が出るはずですが、これはノイズと多重比較の
影響もあり、必ずしも本物の PT 違反を意味しません(後述の解釈参照)。
### Fig 6:コホート別イベント・スタディ
```{r}
#| label: fig-6-cohort-eventstudy
#| fig-cap: "コホート別 ATT(g, t)(Fig 6 再現)"
#| fig-height: 6.5
did::ggdid(att_no_cov, ncol = 2) +
labs(title = "各拡大コホートの ATT(g, t)")
```
**Fig 6 の読み方**: 4 パネル(2014, 2015, 2016, 2019 コホート)で、それぞれ:
- 横軸 = 暦年、縦軸 = ATT(g,t) 推定値
- **介入開始時期(g)の左** = プリトレンド検査(0 付近水平なら PT 仮定支持)
- **介入開始時期(g)の右** = 動学効果(時間とともに効果がどう進化するか)
- バー = 同時信頼区間(`cband=TRUE` で計算)
**注目ポイント**:
- **2019 拡大群** はサンプルが少なく(132 郡)、また介入から観察終了までが 1 年だけなので
信頼区間が大幅に広い
- **2014 群** は最も多くのデータ(978 郡 × 6 年の介入後)があり、CI が最も狭い
- 本シミュレーションの実値(`summary(att_no_cov)` から):
- 2014 コホート: ATT(2014, 2015) = **−4.09** *(CI 上限 −0.17、有意)*、他は CI が 0 を含む
- 2015 コホート: ATT(2015, 2018) = **−4.27** *(有意)*、ATT(2015, 2015) = −5.17(ほぼ有意)
- 2016 コホート: ジグザグで明確な傾向なし(小サンプルのノイズ)
- 2019 コホート: ATT(2019, 2019) = +0.33(介入後 1 年のみ、CI 広く非有意)
### Fig 8:イベント・スタディ集約(共変量なし)
```{r}
#| label: fig-8-aggregate
#| fig-cap: "イベント・スタディ集約 ATT_es(e)(Fig 8 再現)"
#| fig-height: 4.5
agg_es <- did::aggte(att_no_cov, type = "dynamic", na.rm = TRUE)
summary(agg_es)
did::ggdid(agg_es) +
labs(title = "全コホート集約イベント・スタディ(共変量なし)")
```
**Fig 8 の読み方と実際の結果**:
- **横軸 = イベント時間 $e = t - g$**(介入開始からの相対年)
- 4 コホートを共通の横軸に揃え、各 $e$ で利用可能なコホート群から **加重平均** した ATT
- 同時信頼区間(`*` 印は CI が 0 を含まない)
**実値**:
- **総合 ATT(全イベント時間で平均)= −2.342**、SE = 0.934、95% CI = (−4.17, −0.51) **有意**
- **動学効果**:
- $e = 0$(介入直後): ATT = −2.25(まだ CI が 0 を含む)
- $e = 1$: ATT = **−3.92** *(有意、CI 上限 −0.55)*
- $e = 2, 3, 4, 5$: ATT ≈ −2 前後で安定
- **プリトレンド**($e < 0$):大半が 0 周りで CI も 0 を含む → 平行トレンド仮定はおおむね支持
### Fig 9:共変量付き集約イベント・スタディ
```{r}
#| label: fig-9-aggregate-cov
#| fig-cap: "共変量付き集約 イベント・スタディ(Fig 9 再現)"
#| fig-height: 4.5
att_with_cov <- did::att_gt(
yname = "mortality_rate",
tname = "year",
idname = "county_id",
gname = "G_for_did",
data = did_data,
weightsname = "pop_2013",
xformla = ~ pct_white + poverty + median_income, # 共変量で条件付け
control_group = "notyettreated",
panel = TRUE,
bstrap = TRUE,
cband = TRUE,
est_method = "dr" # doubly robust
)
agg_es_cov <- did::aggte(att_with_cov, type = "dynamic", na.rm = TRUE)
summary(agg_es_cov)
did::ggdid(agg_es_cov) +
labs(title = "全コホート集約イベント・スタディ(二重頑健 + 共変量)")
```
**Fig 9 の読み方**: Fig 8 とほぼ同じ形ですが、**共変量(白人率・貧困率・所得)で
条件付け** した推定値です。Fig 8 と比較して:
- 形状(プリトレンド・動学パターン)が **大きく変わらないなら**、観察可能な共変量で説明できる
PT 違反は深刻でない(§4 で見たのと同じメッセージ)
- 大きく変わる場合は、共変量で吸収される交絡があったということ → 本ノートのシミュレーション
では DGP に共変量×$Y(0)$ トレンドの関連を入れていないので、ほぼ同じパターンが出るはず
### 他の推定量との比較
```{r}
#| label: cs-alternatives
#| warning: false
# Sun-Abraham (fixest::sunab で)
sa_data <- did_data |> mutate(G_sa = ifelse(G == Inf, 10000, G))
sa_fit <- feols(mortality_rate ~ sunab(G_sa, year) | county_id + year,
data = sa_data, weights = ~pop_2013, cluster = "county_id"
)
sa_agg <- summary(sa_fit, agg = "att")
sa_att <- coef(sa_agg)[1]
# Borusyak et al. (didimputation)
imp_fit <- tryCatch(
didimputation::did_imputation(
data = did_data |> mutate(treat = treated_now),
yname = "mortality_rate", gname = "G_for_did",
tname = "year", idname = "county_id"
),
error = function(e) NULL
)
# Roth & Sant'Anna (staggered)
stag_fit <- tryCatch(
staggered::staggered(
df = did_data |> rename(t = year, i = county_id),
i = "i", t = "t", g = "G_for_did", y = "mortality_rate",
estimand = "simple"
),
error = function(e) NULL
)
compare_tbl <- tibble(
推定量 = c(
"Callaway-Sant'Anna (集約)",
"Sun-Abraham",
"Borusyak Imputation",
"Roth-Sant'Anna (staggered)"
),
ATT = c(
did::aggte(att_no_cov, type = "simple")$overall.att,
sa_att,
if (!is.null(imp_fit)) imp_fit$estimate[1] else NA_real_,
if (!is.null(stag_fit)) stag_fit$estimate else NA_real_
)
)
compare_tbl |>
gt() |>
fmt_number(decimals = 2)
```
#### 出力の読み方
**比較表(4 推定量)** は同じデータに対して 4 つの異なる現代推定量で「全体集約 ATT」を計算したものです。
- **Callaway-Sant'Anna(集約)**: 上の `aggte(., type="simple")` の値
- **Sun-Abraham**: `fixest::sunab()` を使った代替推定量(コホート×イベント時間 saturated 回帰)
- **Borusyak Imputation**: 補完推定量(まだ介入されていないサンプルで結果モデルを当て、
介入下サンプルの反事実を予測)
- **Roth-Sant'Anna (staggered)**: 効率的推定量
「ほぼ同じ ATT が出るか」をチェックします。
#### 実際の結果の解釈
実際の値は次の通り:
| 推定量 | ATT |
|---|---|
| Callaway-Sant'Anna(集約) | **−2.35** |
| Sun-Abraham | **−2.31** |
| Borusyak Imputation | **+0.19** |
| Roth-Sant'Anna (staggered) | **+41.70** |
**読み取り 1:Callaway-Sant'Anna と Sun-Abraham は驚くほど一致(−2.35 vs −2.31)**
これら 2 つは数学的に近い構造を持ち、また同じデータ準備(`G_for_did` 列、`pop_2013` 加重)で
動かしているので、ほぼ一致するのは自然です。**段階的設定で最も信頼できる中核ペア**。
**読み取り 2:Borusyak Imputation と Roth-Sant'Anna は大きく外れる**
これは「**推定量の正しさの問題**」ではなく「**パッケージごとの入力データ要件の違い**」が
原因の可能性が高いです:
- `didimputation::did_imputation` は **加重なし** で動いており、無加重に近い結果(+0.19)
- `staggered::staggered` は内部のサンプル/コホート定義が異なる可能性があり、
+41.70 という大きな外れ値
実務では:
- これらを **そのまま比較するときはパッケージドキュメントを丁寧に読み合わせる** べき
- 本ノートは「ワンクリックで全部走らせて比較」できることを優先したので、
Borusyak と Roth-Sant'Anna の値は **データフォーマット要件を満たしていない可能性** が
あることを明示しておきます
**読み取り 3:主結果は Callaway-Sant'Anna(−2.35)で据える**
§3 の単純 2×2 DiD(人口加重 −1.88)、§4.4 の二重頑健(人口加重 −2.37)、
§5.2 の Callaway-Sant'Anna(集約 −2.35)とほぼ整合しています。
| 推定 | 値 |
|---|---|
| §3 単純 2×2 DiD(加重) | −1.88 |
| §4.4 二重頑健(2014 群 vs 非拡大、加重) | −2.37 |
| §5.2 Callaway-Sant'Anna(全コホート集約、加重) | −2.35 |
| §5.2 Sun-Abraham(同上) | −2.31 |
これは「**段階的設定でも、適切な推定量(Callaway-Sant'Anna / Sun-Abraham)を使えば 2×2 や 二重頑健 と整合する**」という
最も重要なメッセージです。一方、§5.3 で見るように **ナイーブ TWFE** は別の値を返し、
これが現代の DiD 文献の核心問題となります。
#### 注意点
- **Borusyak / Roth-Sant'Anna の異常値はパッケージ依存** — 比較自体は教育的価値があるが、
実務ではドキュメント精読の上で使うこと
- **`base_period = "varying"` か `"universal"` か** で `did::att_gt` の結果が変わる
(本ノートは `"varying"` を使用 = 各イベント時間で前期との差分)
- 多群でプリトレンドを評価するときは **コホートごと** に見るのが基本
(全体集約だけだと打ち消し合いが起こりうる)
## 5.3 TWFE 回帰の限界
#### このセクションで何をするか
ここまでは **構成要素ベース** の現代推定量(Callaway-Sant'Anna 等)で段階的 DiD を見ました。
しかし実務では古典的に **ナイーブな TWFE 回帰** が広く使われてきました:
$$
Y_{i,t} = \alpha_i + \lambda_t + \beta^{\text{TWFE}} D_{i,t} + \varepsilon_{i,t}
$$
このセクションでは、
1. **ナイーブ TWFE を走らせる**(`feols(Y ~ treated_now | county_id + year)`)
2. **Goodman-Bacon 分解** で TWFE 推定値が **どの 2×2 比較の加重平均** から成り立っているかを可視化
3. 「**Later vs Earlier Treated**(forbidden comparison)」の重みが問題視される理由を示す
#### こんなときに必要
実務での典型例:
- 既存の論文・先行研究が TWFE 回帰を使っている → **再現と比較**
- 査読対応で「なぜ TWFE ではなく Callaway-Sant'Anna を使うのか」を説明する根拠が必要
- **段階的設定で TWFE がどれくらいバイアスを持つか** を診断したい
- 構成要素ベース推定量(§5.2)と並べて **頑健性を確認** したい
#### コードで何をしているか(ナイーブ TWFE)
| ステップ | コード | 何をしているか |
|---|---|---|
| 1 | `feols(Y ~ treated_now \| county_id + year, weights=~pop_2013, cluster="county_id")` | 全データを使い、`treated_now`(介入中ダミー)で TWFE 推定 |
| 2 | `broom::tidy(twfe_naive)` | 係数表で点推定値 + SE + p 値を確認 |
```{r}
#| label: twfe-naive
twfe_naive <- feols(mortality_rate ~ treated_now | county_id + year,
data = did_data, weights = ~pop_2013, cluster = "county_id"
)
broom::tidy(twfe_naive)
```
#### ナイーブ TWFE の出力の読み方
- **`treated_now` 行**の `estimate` が TWFE の DiD 推定量 $\beta^{\text{TWFE}}$
- `std.error` は郡レベルクラスタロバスト SE
- `p.value` が <0.05 なら統計的に有意
#### ナイーブ TWFE の実際の結果
実値:
| 係数 | 推定値 | SE | p 値 |
|---|---|---|---|
| `treated_now` | **−3.08** | 0.63 | < 0.001 |
このシミュレーションでは TWFE 推定値は **−3.08**(SE 0.63)で **統計的に強く有意**。
§5.2 の Callaway-Sant'Anna(−2.35)と比べてやや大きめに出ています(差 ≈ 0.7)。
この差が **forbidden comparison による汚染** から来ているのか、それとも別の理由かを
Bacon 分解で確認します。
### Goodman-Bacon 分解
::: {.callout-note}
## 用語:Bacon 分解(Goodman-Bacon Decomposition)
TWFE 推定量を **全ての 2×2 サブ DiD の加重平均** に分解する手法。
分解の中に「**遅く介入された群を介入群、早く介入された群を比較群に使う 2×2**」が含まれることが
問題で、これが原因で異質効果がある場合に TWFE が誤った推定値を返します。
具体的に分解は次の 3 タイプに分かれます:
- **Treated vs Untreated**: 介入群 vs 非介入(never-treated)群の比較 → 健全
- **Earlier vs Later Treated**: 先に介入された群 vs (まだ)介入されていない群 → 健全
- **Later vs Earlier Treated**: 後で介入された群を「介入」とし、**先に介入された群を「対照」** とする
比較 → 問題の **forbidden comparison**(既に効果が出ている群を対照に使うため)
:::
#### コードで何をしているか(Bacon 分解)
| ステップ | コード | 何をしているか |
|---|---|---|
| 1 | `bd_data <- did_data \|> select(...) \|> arrange(...)` | `bacondecomp` 要件に合わせて整形(無加重・郡 ID 順) |
| 2 | `bacondecomp::bacon(formula, data, id_var, time_var)` | TWFE を 2×2 サブ DiD に分解 |
| 3 | `group_by(type) \|> summarise(weight, avg_est)` | 3 タイプごとに加重平均で集約 |
| 4 | `ggplot(... coord_polar(...) ...)` | 円グラフで重み構成を可視化 |
**注意**: `bacondecomp::bacon` は **加重に対応していない** ので、TWFE の値とは必ずしも一致しません。
形式的には ATT に近い「**unweighted TWFE の分解**」を見ていることに留意してください。
```{r}
#| label: bacon-decomp
#| fig-cap: "Bacon 分解:TWFE = 各 2×2 比較の加重平均"
#| fig-height: 4.5
# bacondecomp は (id, time, treat) を要求
bd_data <- did_data |>
select(county_id, year, mortality_rate, treated_now) |>
arrange(county_id, year)
bd <- bacondecomp::bacon(
formula = mortality_rate ~ treated_now,
data = bd_data,
id_var = "county_id",
time_var = "year"
)
# bacon() の返り値は data.frame または list($two_by_twos, $Omega)。整形を統一
bd_df <- if (is.data.frame(bd)) bd else bd$two_by_twos
bd_df <- as_tibble(bd_df)
bd_summary <- bd_df |>
group_by(type) |>
summarise(
weight = sum(weight),
avg_est = sum(estimate * weight) / sum(weight),
.groups = "drop"
)
bd_summary |>
gt() |>
fmt_number(decimals = 3)
ggplot(bd_summary, aes(x = "", y = weight, fill = type)) +
geom_col(width = 1, color = "white") +
coord_polar(theta = "y") +
scale_fill_brewer(palette = "Set2") +
theme_void(base_family = if (.Platform$OS.type == "windows") "JP" else "") +
labs(
title = "Bacon 分解の重み構成",
fill = "比較タイプ"
)
```
#### Bacon 分解の出力の読み方
**表(`bd_summary`)**:
- **`type` 列**: 3 つの比較タイプ
- **`weight` 列**: 各タイプが unweighted TWFE 推定値に寄与する **重み**(合計 = 1)
- **`avg_est` 列**: 各タイプ内の 2×2 サブ DiD 推定値の(重み付き)平均
**円グラフ**: 重み構成を視覚化。`Later vs Earlier Treated` のスライスが大きいほど
**forbidden comparison の汚染が深刻**。
#### 実際の Bacon 分解結果の解釈
実値:
| 比較タイプ | 重み | 平均推定値 |
|---|---|---|
| **Treated vs Untreated** | **0.831**(83%) | 1.17 |
| Earlier vs Later Treated | 0.118(12%) | 1.80 |
| **Later vs Earlier Treated**(問題) | **0.052**(5%) | 2.83 |
**読み取り 1:本シミュレーションでは forbidden comparison の重みが小さい(5%)**
- 重みの **約 83%** が「Treated vs Untreated」、これは健全な比較
- **「Later vs Earlier Treated」** は **5%** だけ → 大半の TWFE は健全な比較から来ている
- → 本シミュレーションでは TWFE の汚染は **比較的軽微**
**読み取り 2:しかし汚染ゼロではない**
- 5% でも各 ATT が大きく異なれば結果に影響する
- 段階的設定で **コホート数が増える / 介入効果が時間で進化する** 場合、汚染はより深刻になる
- 実データではコホート構成によっては **30〜50% が forbidden comparison** になるケースもある
(Goodman-Bacon 2021 の実例参照)
**読み取り 3:なぜ Callaway-Sant'Anna(−2.35)と TWFE(−3.08)が違うのか?**
§5.2 の CS 推定値(−2.35)と TWFE(−3.08)の差は約 0.7。Bacon 分解が示すように
汚染重みは 5% 程度なので、差の主因は他にもありそう:
- TWFE は **無加重の 2×2 比較** を内部で組み合わせている(`pop_2013` 加重は外側のみ)
- CS は **群サイズで透明に加重** している
- この加重スキームの違いが TWFE と CS の差の主要因
#### 注意点
- **Bacon は加重に対応していない** → 厳密に TWFE を分解できるのは無加重版のみ
- **多群×多期間でコホートが多い設定** では Bacon の解釈はより複雑になる
- 実務では **TWFE → Bacon 分解 → CS / 補完推定量** の順で診断するワークフローが標準
- **forbidden comparison が大きい場合** は TWFE 結果を主結果に使わず、必ず構成要素ベース推定量に切り替え
---
# 結論(論文 §6)
論文の核心メッセージは:
> **DiD 分析は、全て「2×2 DiD の構成要素の集まり」として理解できる。
> まず目標パラメータ(ATT, ATT(g,t), ATT_es(e))を定め、
> それから識別仮定(平行トレンド, 条件付き平行トレンド, PT-GT-NYT 等)を選び、
> 最後にその仮定下で妥当な推定量を選ぶ。**
## 実務家のための 8 ステップ・ワークフロー
::: {.callout-tip}
## DiD 研究のチェックリスト
1. **目標パラメータの定義**: ATT(t)? ATT(g,t)? ATT_es(e)? 政策的にどれが意味あるか?
2. **介入変数の定義**: 段階的か? オン・オフあるか? 連続的か?
3. **共変量バランスの確認**: SMD > 0.2 の変数は?(Table 4 のような表を作る)
4. **平行トレンド仮定の検査**: プリトレンド・プロットを描く(`did::ggdid()`)
5. **平行トレンドの種類選択**: Nev / NYT / all のどれを仮定するか明示
6. **推定量の選択**: TWFE は避け、`did::att_gt()` などの構成要素ベース手法を
7. **集約方法の選択**: dynamic? group? overall?(`did::aggte()`)
8. **頑健性チェック**: 異なる推定量(Sun-Abraham, Borusyak, staggered)で結果比較
:::
---
# Appendix A: 追加の DiD 関連手順 {.appendix}
## A.1 介入がオン・オフを繰り返す設定
ある州が一度メディケイドを拡大した後、廃止し、再開する、などの設定。
このとき介入日 $G_i$ では介入歴を表現できず、より一般的な「介入系列」が必要になります。
```{r}
#| label: a1-on-off
#| fig-cap: "オン・オフ介入の例(模擬)"
#| fig-height: 3.5
set.seed(99)
on_off <- tibble(
year = 2009:2019,
treat = c(0, 0, 0, 1, 1, 0, 0, 1, 1, 1, 0),
Y = c(100, 101, 102, 99, 98, 103, 104, 98, 97, 96, 105)
)
ggplot(on_off, aes(year, Y, fill = factor(treat))) +
geom_col() +
scale_fill_manual(
values = c("0" = "grey70", "1" = "firebrick"),
name = "介入"
) +
labs(title = "介入が複数回オン・オフ:単一の G では表現できない")
```
参考:`did2s::did2s()` 等の手法を参照(本稿では割愛)。
## A.2 連続あるいは多値介入を伴う DiD 設定
「介入の **強度**(税率の変化幅、補助金額など)」を扱う設定。
Callaway, Goodman-Bacon, Sant'Anna (2024) は **線量別 ATT** = $ATT(d \mid d)$ を提案。
```{r}
#| label: a2-continuous
#| fig-cap: "線量(dose)別 ATT の概念図"
#| fig-height: 3.5
dose_df <- tibble(
dose = seq(0, 1, by = 0.05),
ATT = -5 * dose - 3 * dose^2 + rnorm(length(seq(0, 1, by = 0.05)), 0, 0.5)
)
ggplot(dose_df, aes(dose, ATT)) +
geom_point(color = "steelblue") +
geom_smooth(method = "loess", se = FALSE, color = "firebrick") +
labs(
x = "介入線量 d", y = "ATT(d | d)",
title = "線量応答曲線(模擬データ)"
)
```
## A.3 三重差分(Triple Differences, DDD)
「共通の地域ショックを除去するために、追加の比較群(例:性別、年齢層)」を導入。
::: {.callout-note}
## 用語:三重差分(DDD)
DiD を 2 つ取り、その差を取る。具体的には:
- 影響を受けるグループ Q(例:女性):メディケイド拡大 vs 非拡大の DiD
- 影響を受けないグループ Q'(例:男性):メディケイド拡大 vs 非拡大の DiD
- **DDD** = 上の Q の DiD − Q' の DiD
「地域 × 時間」の共通ショックは両方の DiD に同じく入るため、引き算で消える。
:::
```{r}
#| label: a3-ddd
#| code-fold: true
# 仮想的に性別を加えてみる(Q=Female / Q=Male)
set.seed(7)
ddd_data <- did_data |>
filter(year %in% c(2013, 2014), G == 2014 | G == Inf) |>
crossing(female = c(0, 1)) |>
mutate(
# 拡大 × 女性 × Post で大きい効果(例)
extra_effect = -3 * (G == 2014) * (year == 2014) * female,
mortality_rate_gendered = mortality_rate + extra_effect + rnorm(n(), 0, 5)
)
ddd_fit <- feols(mortality_rate_gendered ~ D * post * female,
data = ddd_data |> mutate(D = (G == 2014) * 1L, post = (year == 2014) * 1L),
cluster = "county_id"
)
broom::tidy(ddd_fit) |> filter(grepl(":", term))
```
`D:post:female` の係数が DDD 推定量です。
## A.4 分布 DiD 手順
ATT は **平均** の介入効果ですが、分布全体(分位点)への効果が重要なこともあります。
```{r}
#| label: a4-distributional
#| fig-cap: "分位点 ATT(模擬)"
#| fig-height: 3.5
q_df <- tibble(
quantile = seq(0.1, 0.9, by = 0.1),
q_ATT = -1 + 4 * (seq(0.1, 0.9, by = 0.1) - 0.5)^2 * 10
)
ggplot(q_df, aes(quantile, q_ATT)) +
geom_line(linewidth = 1, color = "firebrick") +
geom_point(size = 2) +
geom_hline(yintercept = 0, linetype = "dashed") +
labs(x = "分位点", y = "Q-ATT", title = "分布 DiD:分位点ごとの ATT")
```
`qte` パッケージや Callaway-Li (2019) を参照。
## A.5 繰り返しクロスセクションと不均衡パネル・データ
「同じ郡が毎年観察される」(バランス・パネル)とは限らない場合の手当て:
- **繰り返しクロスセクション**: 各年で異なる個人が観察される(調査データ等)
- **不均衡パネル**: 一部の単位で欠損年がある
`DRDID::drdid_rc()` は繰り返しクロスセクション版の 二重頑健 DiD を提供します:
```{r}
#| label: a5-rc
# 繰り返しクロスセクションをエミュレート
set.seed(11)
rc_data <- sub3 |>
group_by(county_id) |>
slice_sample(n = 1) |> # 各郡で 1 年だけ残す
ungroup()
# DRDID の RC 版(panel = FALSE で繰り返しクロスセクション)
res_rc <- tryCatch(
DRDID::drdid(
yname = "mortality_rate", tname = "year",
idname = "county_id", dname = "D",
xformla = ~ pct_white + poverty + median_income,
data = rc_data, panel = FALSE
),
error = function(e) NULL
)
if (!is.null(res_rc)) {
cat(
"RC-DRDID ATT =", round(res_rc$ATT, 2),
" SE =", round(res_rc$se, 2), "\n"
)
} else {
cat("(繰り返しクロスセクションでは DRDID パッケージ要件に応じた整形が必要)\n")
}
```
---
# 参考文献 {.unnumbered}
主要な参考文献(本ノートが立脚した論文):
- Baker, A., Callaway, B., Cunningham, S., Goodman-Bacon, A., & Sant'Anna, P. H. C. (2025).
*Difference-in-Differences Designs: A Practitioner's Guide.* arXiv:2503.13323v3.
- Callaway, B., & Sant'Anna, P. H. C. (2021).
*Difference-in-Differences with multiple time periods.* Journal of Econometrics, 225(2), 200-230.
- Goodman-Bacon, A. (2021).
*Difference-in-differences with variation in treatment timing.* Journal of Econometrics, 225(2), 254-277.
- de Chaisemartin, C., & D'Haultfœuille, X. (2020).
*Two-way fixed effects estimators with heterogeneous treatment effects.* American Economic Review, 110(9), 2964-2996.
- Sun, L., & Abraham, S. (2021).
*Estimating dynamic treatment effects in event studies with heterogeneous treatment effects.* Journal of Econometrics, 225(2), 175-199.
- Borusyak, K., Jaravel, X., & Spiess, J. (2024).
*Revisiting event study designs: Robust and efficient estimation.* Review of Economic Studies.
---
# セッション情報 {.unnumbered .appendix}
```{r}
#| label: session-info
sessionInfo()
```
