---
title: "個人レベルデータで実施する Event DID:震災の Treatable Death への効果"
subtitle: "福島・岩手・宮城 vs 44都道府県 — 「個人データそのまま」と「集計データに加工」の比較(再現可能チュートリアル / デモデータ)"
author: "Event DID 学習用チュートリアル"
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
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}
::: {.callout-important}
## 目的と前提
- 震災(東日本大震災・2011年3月)が **全年齢の treatable death(extended-treatable mortality)** に与えた効果を、**個人レベルデータ**で推定する **Event DID(event-study difference-in-differences)** を、コードを動かしながら学びます。
- **「個人データをそのまま使う解析」と「集計データに加工して使う解析」を並べて比較**することが主眼です。
- データは **統計法に基づき政府統計の個票(人口動態統計 死亡票)+ 国勢調査/推計人口の分母を申請取得した**という現実的シナリオを想定し、その形の **デモ(シミュレーション)データ**を本ノート内で生成します。
- **本ノートの数値は完全なシミュレーションであり、実際の福島・岩手・宮城の死亡率や震災の効果を主張するものではありません。** 真のデータ生成過程(DGP)を既知にすることで、各推定量の挙動を理解しやすくする教育目的です。
:::
::: {.callout-tip}
## 読み方
- **コード**は折り畳めます(右上の `</>` → `Show All Code` で全展開)。
- **専門用語**は色付き callout box で初出時に説明します。
- 比較の主軸は「**手法A:個人データそのまま(person-time の Poisson / 線形確率モデル)**」 vs 「**手法B:集計データに加工(都道府県×年 ASR の Callaway–Sant'Anna DID)**」です。
:::
## 環境準備:パッケージ
```{r}
#| label: setup-packages
#| code-fold: true
#| code-summary: "パッケージのインストール・ロード(クリックで展開)"
pkgs <- c(
"tidyverse", # データ操作・可視化の基盤
"fixest", # 高速な固定効果推定(個人レベルの event-study に使用)
"did", # Callaway-Sant'Anna 推定量(集計パネルの中核)
"broom", # モデル結果を tidy データへ
"modelsummary", # 結果の整形表
"gt", # 整形テーブル
"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))
# 重み付き(分数)応答の Poisson で出る情報ノートを抑制(エラーではない)
fixest::setFixest_notes(FALSE)
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)
# 全体で固定する乱数シード(再現性)
set.seed(20260604)
```
---
# 推定したいもの(estimand)と研究デザイン
> **この章のねらい** — 解析を始める前に「何を因果効果として推定するのか(estimand)」を定義する。アウトカム(全年齢 treatable death)、比較群(被災3県 vs 44県)、効果の指標 ATT_k、そして DID がなぜ ATT を識別できるのか(前提)を確認する。
## アウトカム:全年齢 treatable death(extended-treatable mortality)
本研究の**主要アウトカム**は、OECD/Eurostat の "treatable causes of death" リスト(= 適時かつ適切な医療介入で回避しうる死亡)に該当する死亡を、**全年齢**に適用したもの。本ノートではこれを **treatable death**(extended-treatable mortality)と呼びます。本研究に固有の定義上の決定は次の3点です(treatable / preventable / treatable share の一般的な定義は下記「用語」ボックスを参照)。
- **年齢範囲**:OECD 標準の **0–74 歳上限を撤廃し全年齢**に適用した独自構築定義(震災関連死が 75 歳以上に集中するため)。
- **数え方**:単純な総死亡数ではなく、各死亡を **treatable share(治療可能割合)で重み付けした合計 `Σ treatable_share`** として定義。本ノートはこの**重み付け定義のみ**を用い(本番パイプラインと同一)、二値化(share>0 を一律1件)は使いません。具体的な実装は後述の「treatable share の重み付け」で扱います。
- **指標**:**年齢標準化死亡率(ASR, age-standardized rate;人口10万人・年あたり)**。標準人口は **OECD 2010 standard population**。
::: {.callout-note}
## 用語:avoidable / preventable / treatable mortality と treatable share
OECD/Eurostat は「**回避可能死亡(avoidable mortality)**」を、介入の種類によって2つに二分し、各死因をどちらか一方、または両方に割り当てます。
- **予防可能(preventable)**:喫煙対策・予防接種・事故防止など**一次予防**で防げる死亡(=「行動・生活習慣・社会経済・環境など広い公衆衛生要因への介入で回避しえた死亡」)。
- **治療可能(treatable)**:発症後の**適時・適切な医療**で防げる死亡(=「最適な質の医療で回避しえた死亡」)。
**treatable share(治療可能割合)** とは「その死因の死亡のうち treatable 側に数える割合」で、標準リストでは **{0, 0.5, 1.0}** の3値をとります。
- **1.0(完全に treatable)**:その死因の死亡を**まるごと1件**として treatable mortality に算入。例 — 肺炎・敗血症・大腸癌・乳癌・急性虫垂炎(医療アクセスがあれば多くが救命可能な死因)。
- **0.5(予防可能・治療可能の折半)**:OECD が予防・治療の両方に該当すると判断した死因。treatable としては**半分(0.5件)**だけ寄与し、残り 0.5 は preventable 側に計上される。例 — 虚血性心疾患(急性心筋梗塞)・脳血管疾患・糖尿病・高血圧性疾患・子宮頸癌。
- **0(treatable ではない)**:純粋に予防可能(例:肺癌・交通事故)、または回避可能リスト外の死因。treatable mortality には**算入しない**。
なお上記 {0, 0.5, 1.0} は OECD リストが**個々の死因**に与える分類であり、**複数の ICD をまとめた死因(レンジ/グループ)では配下コードの加重平均となって連続値**を取ります(本デモの share も本番ブリッジに準じた連続値を割り当て)。
**根拠**:
- OECD/Eurostat. *Avoidable mortality: OECD/Eurostat lists of preventable and treatable causes of death*(2019年版リスト、2022年1月公表). <https://www.oecd.org/content/dam/oecd/en/data/datasets/oecd-health-statistics/avoidable-mortality-2019-joint-oecd-eurostat-list-preventable-treatable-causes-of-death.pdf> — 予防/治療の優劣に明確な根拠がない死因(虚血性心疾患・脳卒中・糖尿病など)は二重計上を避けるため **50%–50% で両カテゴリに按分**し、回避可能死亡の年齢上限は **0–74 歳**(本研究はこの上限を撤廃)。
- Eurostat. *Preventable and treatable mortality statistics — Statistics Explained.* <https://ec.europa.eu/eurostat/statistics-explained/index.php?title=Preventable_and_treatable_mortality_statistics> — treatable / preventable の定義、および上記 50/50 按分ルールの出典。
:::
## 比較群:被災3県 vs 44非被災県(nested natural experiment の外側)
- **処置群(treated cohort)**:福島(07)+ 岩手(03)+ 宮城(04)の **3県**。3県とも **2011年に同時に処置**(= 非段階的 / single-timing)。
- **対照群(never-treated)**:残り **44都道府県**。
- 処置時点が全県同一なので、本研究は **段階的採択(staggered adoption)ではなく 2群 × 多時点**の event-study です。
::: {.callout-note}
## 用語:ATT_k(event-time別の平均処置効果)
イベント(2011年)からの相対年 $k = \text{year} - 2011$ ごとに、処置群が「処置を受けなかった場合の反実仮想」と比べてどれだけ outcome が変化したかを表す量:
$$
\text{ATT}_k = \mathbb{E}\big[\, Y_{i,\,2011+k}(\text{処置}) - Y_{i,\,2011+k}(\text{無処置}) \;\big|\; i \in \text{被災3県} \,\big]
$$
- $k < 0$(処置前)の ATT_k は **平行トレンド(parallel trends)検定**に使う(理想的には 0 近傍)。
- $k = -1$(= 2010年)を **基準期(reference period)**に固定し、すべての ATT_k はこの基準との相対値。
:::
::: {.callout-important}
## 重要:縦断的に個人を追跡したパネルは**不要**(repeated cross-section で DID は成立)
「政府統計は個人を年をまたいで紐づけられないが、Event DID は可能か?」 → **可能です。**
- DID/event-study の識別に必要なのは「**処置割当の単位(=都道府県)が各年観測されている**」ことであり、**同一個人を追跡すること**ではありません。
- 人口動態統計の **死亡票は年ごとの断面イベント**(個人は年をまたいで紐づかない)、分母は国勢調査/推計人口の集計値 — この「**反復横断(repeated cross-section)**」で DID は成立します。
- Callaway–Sant'Anna の `did::att_gt()` も `panel = FALSE` で反復横断を正式サポートします(その場合 `idname` は無視)。本ノートの集計版は**都道府県**を単位とするバランスパネル(`panel = TRUE`)、個人版は**反復横断**です。<https://bcallaway11.github.io/did/reference/att_gt.html>
- 個人が毎年入れ替わることによる**人口構成の変化**は、**個人版=年齢階級・性別の共変量調整**、**集計版=直接法による年齢標準化**で吸収します(両者の点推定が一致する理由でもあります)。
:::
## FAQ
::: {.callout-tip}
## FAQ:個人レベルデータで treatable death は計算できる?
**できます。むしろ個票(個人レベル)の方が自然です。** 「treatable かどうか」は **死因(ICD-10)の属性**なので、1行=1死亡の個票なら各死亡に `treatable_share`(死因ごとの治療可能割合)を **1件ずつ直接付与**できます(集計データは先に率へ丸めるため、この死因単位の情報を後から復元できません)。
- **率にするには 分子 + 分母**:分子 = 個票の treatable 死亡 = `sum(treatable_share)`(各死亡をその share で重み付け;0.5 死因は 0.5 件)、分母 = 国勢調査/推計人口の人年。両者を突き合わせると「**その年の treatable 寄与(0、または死因の share)**」という person-year レコードになる。
- **「個人レベル」≠ 個人を追跡**:死亡票は年ごとの断面(反復横断)。treatable 判定は断面で死因から行うので、同一人物の年次リンクは不要。
- **必要な変数は3つ**:① 死因(ICD-10)→ treatable 判定、② 年齢 → 標準化(本研究は OECD の 0–74 歳上限を撤廃し全年齢に適用)、③ 分母人口(別ソースの集計)。死亡票は ①② を備えるため、分子側は個票だけで完結する。
- **個人展開 ≡ セル Poisson(lossless)**:1行=1人の線形確率モデルと、同一共変量をまとめたセル形式 Poisson は同じ尤度。treatable は稀(≈0.2%)なのでセル Poisson で軽量化しつつ、年齢・性・死因の情報を保持する。
実装は後述の §「データ:政府統計申請で得る『個票』と『分母』を模す」および §「『1行=1人』と『セル集計』は同じ尤度(lossless)」で示す。
:::
::: {.callout-note collapse="true"}
## なぜ「年齢標準化」死亡率か(単純な死亡率ではだめか)/ なぜ OECD 2010 標準人口か
### 単純(粗)死亡率ではなぜだめか
- **粗死亡率(crude rate)= 総死亡数 / 総人口** は、その集団の**年齢構成に強く依存**します。treatable death は加齢とともに急増するため、**年齢別リスクが全く同じでも、高齢な集団ほど粗率が高く**出ます。
- 被災3県(特に東北)は全国平均より高齢 → 粗率の差には「震災の効果」と「もともとの年齢構成の差」が**混在**してしまう。
- DID は県固定効果で**時間不変**の交絡(恒常的な年齢構成差)を、年固定効果で全国共通の変化を除きますが、**年齢構成は時間不変ではありません**。高齢化の進み方や、震災後の避難(若年世帯の流出・高齢者の残留)で**処置群と対照群の年齢構成が異なるトレンド**を持つと、粗率の DID は**平行トレンド前提を破り、バイアス**が生じます。
- → **年齢標準化**(または個人モデルでの年齢調整)で「年齢構成の違い・その経時変化」を取り除き、**純粋な年齢別リスクの差**だけを比較します。
- 本ノートの対応関係:**集計版=直接法による年齢標準化(ASR)**、**個人版=年齢階級の固定効果/共変量**。両者は同じ交絡(年齢)を別の方法で除いており、後述で**点推定が一致する理由**になります。
### 直接法標準化と「標準人口」とは
$$
\text{ASR} = \sum_a w_a \, m_a \times 100{,}000
$$
- $m_a$ = 年齢階級 $a$ の死亡率、$w_a$ = **標準人口**の年齢重み($\sum_a w_a = 1$)。標準人口は全県・全年に共通の「ものさし」です。
- 重要:**同一の標準人口をすべての県・年に適用する限り**、どの標準を選んでも**相対比較(% 効果)はほぼ不変**です。標準人口の選択が効くのは ①**ASR の絶対水準**、②**国際比較可能性**、③**高齢層への重み付け**の3点です。
### 標準人口の選択肢
| 標準人口 | 由来・用途 | 高齢層の重み |
|---|---|---|
| **OECD 2010**(本研究で採用) | OECD/Eurostat が avoidable(preventable/treatable)mortality の算出に使用 | 中〜高(現実的) |
| 日本 1985 年(昭和60年)モデル人口 | 厚労省の従来標準。**Fukuma 2017 が使用** | 低(若い構造) |
| 日本 2015 年(平成27年)モデル人口 | 厚労省が 2020 年統計以降に採用 | 高 |
| WHO 世界標準人口 | 世界平均ベースの国際比較用 | 低 |
| 欧州標準人口(ESP 2013) | Eurostat が欧州比較で使用 | 中〜高 |
### なぜ OECD 2010 が良いか(本研究の文脈)
- **outcome 定義との一貫性**:treatable mortality の死因リスト自体が OECD/Eurostat 2022 由来。OECD が同指標を算出する標準人口と揃えることで定義が内的に整合し、**公表されている OECD の avoidable mortality と直接比較**できる。
- **高齢層を過小評価しない**:WHO 世界標準や日本 1985 年は若い構造で 75 歳以上の重みが小さい。本研究は震災関連死・treatable death が集中する **75 歳以上を捉えるため OECD 標準の under-75 上限を撤廃**した設計なので、高齢層に十分な重みを置く OECD 2010 が適切。
- **再現性・監査可能性**:推定を要しない公表固定表であり、解析前に固定(事前指定)できる。誰が計算しても同じ重みになる。
- **頑健性は感度分析で確認**:標準人口選択への依存は、**日本 1985 年(Fukuma 2017 との比較)・under-75 cap・WHO 標準**などを並行報告して点検します(本ノートはコア比較に集中するため割愛)。
根拠:OECD/Eurostat 2022(前掲);直接法標準化率の信頼区間は Fay & Feuer 1997(後掲 References);日本のモデル人口は厚生労働省「人口動態統計」標準。
:::
::: {.callout-tip}
## Q. 個人データで年齢を調整すれば、年齢標準化(ASR)は要らないのでは?
**ほぼその通りです。** 個人レベルで年齢を調整するのは正しい方針で、本ノートの **手法A がまさにそれ**(年齢階級の固定効果/共変量による調整)。**わざわざ標準化率に丸めてから DID する必要はありません。** ただし「個人データでの年齢調整」と「直接法による年齢標準化(ASR)」は**近いが完全に同一ではない**点に注意します。
- **両者の違い = 重みの出どころ**
- **直接法標準化(ASR)**:年齢別の率を **外部の標準人口(OECD 2010)の重み**で合成。「もし人口構成が標準人口だったら」という解釈で、**外部基準に対して固定**され国際・経時比較が可能。
- **個人モデルの年齢 FE**:年齢の主効果を吸収し、ATT は **データ内の人年(= その集団の年齢構成)で重み付けした**年齢別効果の平均。重みは外部標準ではなく**観測データ由来**。
- **いつ一致し、いつズレるか**
- 震災の効果が**年齢でほぼ一様**なら、年齢 FE の ATT と標準化 ATT はほぼ一致(本ノートで 手法A ≈ 手法B となる理由 → 後述「ブリッジ」)。
- 効果が**年齢で異質**(例:高齢ほど大きい)で、かつ**群間で年齢構成が異なる/異なるトレンドを持つ**場合、「年齢 FE だけ」の平均(データ重み)と「標準人口重み」の平均は**別の数値**になりうる。
- **個人モデルで標準化 ASR と同じ estimand を出すには**(= 回帰標準化 / g-computation)
1. `i(event_time, nested_treated)`×**年齢階級**の**交互作用**まで入れて**年齢別 ATT_k**を推定、
2. 予測した年齢別(反実仮想)率を **OECD 2010 標準人口の重み $w_a$ で加重平均** → 直接法標準化 ATT を**厳密に再現**できる。
- つまり標準化は「個人モデルの後段で**外部重みを掛けて周辺化する操作**」とみなせます。`age_band` を主効果(FE)で入れるだけだと、この周辺化が**データ重みで暗黙に**行われている、という関係です。
- **年齢調整にも「分母」は必須**:個人で年齢調整する場合も、率を推定するには**年齢別の人年(分母)**が要ります(死亡個票だけでは率になりません)。本ノートの person-time / 人年セルがその役割を担います。
- **では なぜ ASR も報告するのか**:因果推定に必要だからではなく、**①単一の率として解釈しやすい**、**②外部基準で国際・経時比較できる**(公表 OECD 値と並べられる)、**③伝達しやすい**から。すなわち**報告・比較の都合**であって、妥当性の要件ではありません。
**結論**:年齢は**個人データで調整すれば十分**(手法A)。標準化 ASR(手法B)は「個人モデルを外部標準人口で周辺化した特別な場合」に対応し、両者の点推定が一致することは後述「ブリッジ」で実証します。
:::
::: {.callout-note collapse="true"}
## なぜ DID 推定量は ATT を識別するのか(導出)
**結論を先に**:DID が計算する $(\Delta_{treated}) - (\Delta_{control})$ は、**平行トレンド + 非予期 + 波及なし**の仮定の下で **ATT(処置群における平均処置効果)に一致**します。ATE(全員に対する平均効果)ではない点が肝心です。
### 1. 潜在結果と「観測できない反実仮想」
各単位 $i$・時点 $t$ に 2 つの潜在結果 $Y_{it}(1)$(処置あり)/ $Y_{it}(0)$(処置なし)を考えます。観測できるのは実際に起きた一方だけ。被災3県の post 期について知りたいのは
$$
\text{ATT} = \mathbb{E}\!\left[\,Y_{post}(1) - Y_{post}(0)\mid \text{treated}\,\right]
$$
第1項 $\mathbb{E}[Y_{post}(1)\mid treated]$ は**観測可能**(実際に処置を受けた被災3県の post 値)。問題は第2項 **$\mathbb{E}[Y_{post}(0)\mid treated]$=「被災3県がもし震災を受けなかったら」の post 値** で、これは**決して観測できません**(因果推論の根本問題)。DID の仕事は**この反実仮想を埋める**ことです。
### 2. 平行トレンド仮定(鍵となる識別仮定)
> もし処置がなかったら、**被災3県の outcome の変化量**は、**44対照県の変化量と同じ**だったはず。
$$
\mathbb{E}\!\left[\,Y_{post}(0)-Y_{pre}(0)\mid treated\,\right]
=
\mathbb{E}\!\left[\,Y_{post}(0)-Y_{pre}(0)\mid control\,\right]
$$
加えて **非予期(no anticipation)**:処置は pre 期の outcome に影響しない。これにより $k=-1$(2010年)を**汚染されていない基準**として使えます。
### 3. 反実仮想を「対照群の変化」で代入する
非予期から $\mathbb{E}[Y_{pre}(0)\mid treated]=\mathbb{E}[Y_{pre}\mid treated]$(被災3県の pre は観測値そのもの)。平行トレンドから、欠けている post の反実仮想は
$$
\mathbb{E}[Y_{post}(0)\mid treated]
= \underbrace{\mathbb{E}[Y_{pre}\mid treated]}_{\text{被災3県 pre(観測)}}
+ \underbrace{\big(\mathbb{E}[Y_{post}\mid control]-\mathbb{E}[Y_{pre}\mid control]\big)}_{\text{対照県の変化(観測)}}
$$
### 4. ATT に代入すると DID になる
$$
\begin{aligned}
\text{ATT}
&= \mathbb{E}[Y_{post}\mid treated] - \mathbb{E}[Y_{post}(0)\mid treated]\\[2pt]
&= \big(\mathbb{E}[Y_{post}\mid treated]-\mathbb{E}[Y_{pre}\mid treated]\big)
- \big(\mathbb{E}[Y_{post}\mid control]-\mathbb{E}[Y_{pre}\mid control]\big)\\[2pt]
&= \underbrace{\Delta_{treated} - \Delta_{control}}_{=\ \text{DID}}
\end{aligned}
$$
→ **すべて観測量だけ**で書けました。これが「DID 推定量 = ATT」の正体です。回帰式の $\beta$($Treated\times Post$ の係数)は、ちょうどこの差を推定しています。
### 5. なぜ ATT であって ATE ではないか
平行トレンドが反実仮想を埋めてくれるのは **被災3県についてだけ**(対照県の変化を借りてくる)。対照県が「もし処置されたらどうなるか」については何も言っていないので、効果が群で異質なら **ATT ≠ ATE**。DID が答えるのは常に「**処置を受けた群における**効果」です。
### 6. event-study(ATT_k)への一般化
post を各 event-time $k$、pre を基準期 $k=-1$ に置き換えるだけ:
$$
\text{ATT}_k = (\,\bar Y_{k}-\bar Y_{-1}\,)_{treated} - (\,\bar Y_{k}-\bar Y_{-1}\,)_{control}
$$
- $k<0$ の ATT_k が **0 近傍**なら平行トレンドと整合(**間接的な検証であって証明ではない**)。
- $k\ge 0$ の ATT_k が処置後の動学(即時 → 調整 → 再均衡)。
### 7. 成立に必要な前提(本研究での留意点)
- **平行トレンド**:反実仮想なので**直接は検証不能**(pre-trend で間接点検)。共変量で条件付けて成立を狙うのが CS-DID の回帰調整。
- **非予期**:震災は2011年に外生的に発生 → 妥当。
- **波及なし(SUTVA)**:対照県が処置の影響を受けないこと。**避難者の流入**で近隣県の分母・死亡が動くと違反 → 本研究では大規模避難受入県を比較群(donor)から除外。
- **構成の安定/調整**:反復横断では各年の標本が母集団を代表し、年齢構成の差は年齢調整(標準化/年齢FE)で吸収。
:::
:::: {.callout-note}
## なぜ個人が紐づかなくても推定できるのか
DID が比較するのは、同じ個人の前後変化そのものではなく、**処置群と対照群における集団平均の前後変化の差**です。
$$
\text{DID}
= \left(\bar{Y}_{treated, post} - \bar{Y}_{treated, pre}\right)
- \left(\bar{Y}_{control, post} - \bar{Y}_{control, pre}\right)
$$
この式で必要なのは、各時点の「処置群の平均」と「対照群の平均」を推定できることです。したがって、2010年の個人と2011年の個人が同一人物としてリンクされていなくても、各年の標本がその県・年の母集団を代表していれば、反復横断データとして DID を推定できます。
本ノートでいう **個人レベル DID** は、次のような「個人レコードを使った反復横断 DID」です。
$$
Y_{igt}
= \alpha
+ \beta \left(Treated_g \times Post_t\right)
+ \lambda_g
+ \tau_t
+ X_{igt}'\gamma
+ \varepsilon_{igt}
$$
ここで $i$ は個人、$g$ は都道府県、$t$ は年です。$i$ が年をまたいで同じ人物である必要はありません。推定している $\beta$ は、**同一個人の変化**ではなく、**処置県の集団平均が、対照県と比べてどれだけ変化したか**です。
| できること | 個人が年次リンクされない反復横断データ |
|---|---|
| 集団平均の DID / event-study | 可能 |
| 年齢・性別などの個人共変量調整 | 可能 |
| 年齢層・性別などの異質性分析 | 可能 |
| 同一個人の前後変化 | 不可 |
| 個人固定効果 DID | 不可 |
したがって、「サマリーデータ化しないと DID できない」わけではありません。サマリーデータ化は、県×年などの平均・率を**先に作る**実装です。一方、個人レベル DID は、個人行を残したまま回帰の中で県×年の平均差を取り出す実装です。処置が県単位で、重み付け・標準化・共変量調整が対応していれば、両者は非常に近い estimand を見ています。
::::
---
# データ:政府統計申請で得る「個票」と「分母」を模す
> **この章のねらい** — 以降の手法A・手法B 双方の入力データを用意する。政府統計の申請で実際に得られる形(死亡個票 + 人口分母)を模したデモデータを生成し、真の効果(DGP)を既知にして推定量の挙動を検証できるようにする。
実データでは、(1) **人口動態統計の死亡票(個票)**= 1行1死亡(都道府県・年・年齢・性・死因 ICD-10)と、(2) **国勢調査/推計人口**= 都道府県×年×年齢×性の人口(分母)を別々に取得します。ここでは同じ構造のデモデータを生成します。
## データの流れ(パイプライン概要)
データは**一本道の前方フロー**で、最後に各推定器へ `data =` で**1つのデータフレーム**を渡します(手法A=`cell_an`/`cell_an_adj`、手法B=`did_panel`)。生成順とともに、どのデータセットがどの推定に渡るかを示します。
```text
[生成] cells (県×年×年齢×性 = 33,840 セル;セルごとに Poisson で death 発生)
│
├─→ deaths_micro (個票 1行=1死亡, 569,448 行) ← 実データ: 人口動態統計の死亡票
└─→ population_py (分母 県×年×年齢×性, 33,840) ← 実データ: 国勢調査/推計人口
│
└─(突合 + treatable_share で重み付け)→ cells_w (33,840; 重み付き treatable 死亡)
│
├─→ cell_an (33,840) ───────────────────────→ 手法A 主 : fepois(率) / feols(LPM)
│ └─→ cell_an_adj (+ baseline SES z 4本) → 手法A 感度: fepois(SES調整)
│
└─→ asr_panel (県×年 ASR, 940) ──→ did_panel (940) → 手法B : did::att_gt
└─→ plot_dat (記述トレンド図) / ri_delta (randomization inference)
```
| 段階 | データセット | 単位 / 行数 | 渡す先(推定器・図) |
|---|---|---|---|
| 生データ相当 | `deaths_micro` | 個票 / 569,448 | (集計して `cells_w` へ) |
| 生データ相当 | `population_py` | 県×年×年齢×性 / 33,840 | (分母として `cells_w` へ) |
| 共有ベース | `cells_w` | 県×年×年齢×性 / 33,840 | `cell_an` / `asr_panel` を派生 |
| 解析入力 | `cell_an` | 同上 / 33,840 | **手法A 主**:`fepois`(mA_pois) ＋ `feols`(mA_lpm) |
| 解析入力 | `cell_an_adj` | 同上(+SES z) / 33,840 | **手法A 感度**:`fepois`(mA_pois_adj) |
| 共有中間 | `asr_panel` | 県×年 ASR / 940 | `did_panel`・`plot_dat`・`ri_delta` を派生 |
| 解析入力 | `did_panel` | 県×年 / 940 | **手法B**:`did::att_gt` → `aggte` |
要点:**推定に直接渡すデータセットは 3 つ**(`cell_an`・`cell_an_adj`・`did_panel`)で、系統で見れば**手法A・手法B の 2 系統**。すべては `cells_w`(個人側)と `asr_panel`(集計側)の 2 つの共有データから分岐します。実データでは `deaths_micro` と `population_py` が申請取得する 2 ファイルに対応します。
## 都道府県マスターと2010年共変量
外部参照値として総務省統計局 2010年国勢調査、内閣府県民経済計算など利用します。
```{r}
#| label: prefecture-master
#| code-fold: true
#| code-summary: "47都道府県マスター + 2010年人口・共変量(参照値)"
prefecture_master <- tibble::tribble(
~prefecture_code, ~prefecture_name, ~region,
"01", "Hokkaido", "Hokkaido", "02", "Aomori", "Tohoku", "03", "Iwate", "Tohoku",
"04", "Miyagi", "Tohoku", "05", "Akita", "Tohoku", "06", "Yamagata", "Tohoku",
"07", "Fukushima", "Tohoku", "08", "Ibaraki", "Kanto", "09", "Tochigi", "Kanto",
"10", "Gunma", "Kanto", "11", "Saitama", "Kanto", "12", "Chiba", "Kanto",
"13", "Tokyo", "Kanto", "14", "Kanagawa", "Kanto", "15", "Niigata", "Chubu",
"16", "Toyama", "Chubu", "17", "Ishikawa", "Chubu", "18", "Fukui", "Chubu",
"19", "Yamanashi", "Chubu", "20", "Nagano", "Chubu", "21", "Gifu", "Chubu",
"22", "Shizuoka", "Chubu", "23", "Aichi", "Chubu", "24", "Mie", "Kinki",
"25", "Shiga", "Kinki", "26", "Kyoto", "Kinki", "27", "Osaka", "Kinki",
"28", "Hyogo", "Kinki", "29", "Nara", "Kinki", "30", "Wakayama", "Kinki",
"31", "Tottori", "Chugoku", "32", "Shimane", "Chugoku", "33", "Okayama", "Chugoku",
"34", "Hiroshima", "Chugoku", "35", "Yamaguchi", "Chugoku", "36", "Tokushima", "Shikoku",
"37", "Kagawa", "Shikoku", "38", "Ehime", "Shikoku", "39", "Kochi", "Shikoku",
"40", "Fukuoka", "Kyushu", "41", "Saga", "Kyushu", "42", "Nagasaki", "Kyushu",
"43", "Kumamoto", "Kyushu", "44", "Oita", "Kyushu", "45", "Miyazaki", "Kyushu",
"46", "Kagoshima", "Kyushu", "47", "Okinawa", "Okinawa"
)
# 2010年総人口(人)
population_2010 <- tibble::tribble(
~prefecture_code, ~pop_2010,
"01", 5506419, "02", 1373339, "03", 1330147, "04", 2348165, "05", 1085997,
"06", 1168924, "07", 2029064, "08", 2969770, "09", 2007683, "10", 2008068,
"11", 7194556, "12", 6216289, "13", 13159388, "14", 9048331, "15", 2374450,
"16", 1093247, "17", 1169788, "18", 806314, "19", 863075, "20", 2152449,
"21", 2080773, "22", 3765007, "23", 7410719, "24", 1854724, "25", 1410777,
"26", 2636092, "27", 8865245, "28", 5588133, "29", 1400728, "30", 1002198,
"31", 588667, "32", 717397, "33", 1945276, "34", 2860750, "35", 1451338,
"36", 785491, "37", 995842, "38", 1431493, "39", 764456, "40", 5071968,
"41", 849788, "42", 1426779, "43", 1817426, "44", 1196529, "45", 1135233,
"46", 1706242, "47", 1392818
)
# 2010年 共変量(65歳以上比率, GDP/人[千円], 人口密度[人/km^2])
# gdp_pc(所得)・pop_density(人口密度)は実データの参照値。
# 学歴(edu_low)・失業(unemp_rate)は本デモでは下で決定論的に生成(実データでは
# 国勢調査の学歴別人口・労働力人口から取得)。いずれも 2011年震災「前」の 2010年断面で
# 固定する baseline 共変量であり、因果経路上(mediator)ではない点が重要(過調整回避)。
covar_2010 <- tibble::tribble(
~prefecture_code, ~elderly_share, ~gdp_pc, ~pop_density,
"01", 0.247, 3895, 70, "02", 0.254, 2790, 142, "03", 0.275, 3267, 87,
"04", 0.220, 3530, 322, "05", 0.296, 2918, 93, "06", 0.273, 3061, 125,
"07", 0.252, 3801, 147, "08", 0.222, 4174, 487, "09", 0.227, 4121, 313,
"10", 0.236, 4198, 315, "11", 0.196, 3203, 1894, "12", 0.210, 3324, 1206,
"13", 0.205, 7044, 5995, "14", 0.201, 3514, 3744, "15", 0.262, 3613, 188,
"16", 0.265, 4221, 257, "17", 0.241, 4015, 280, "18", 0.249, 3949, 192,
"19", 0.250, 3633, 193, "20", 0.265, 3551, 159, "21", 0.241, 3589, 196,
"22", 0.241, 4351, 484, "23", 0.205, 5252, 1432, "24", 0.247, 3826, 320,
"25", 0.214, 4221, 351, "26", 0.228, 3681, 571, "27", 0.221, 4119, 4670,
"28", 0.231, 3568, 666, "29", 0.241, 3055, 379, "30", 0.272, 3066, 211,
"31", 0.262, 3127, 168, "32", 0.293, 3151, 107, "33", 0.252, 3661, 273,
"34", 0.235, 4135, 337, "35", 0.279, 4087, 237, "36", 0.265, 3522, 189,
"37", 0.255, 3700, 530, "38", 0.262, 3196, 252, "39", 0.290, 3052, 107,
"40", 0.219, 3676, 1018, "41", 0.249, 3147, 348, "42", 0.262, 2982, 343,
"43", 0.262, 3041, 245, "44", 0.275, 3470, 188, "45", 0.270, 2814, 146,
"46", 0.266, 2854, 185, "47", 0.182, 2483, 610
)
# baseline 2010 SES の追加(文献ベース調整変数)。RNG を使わず決定論的に生成し、既存の
# 所得(gdp_pc)と相関させる(低所得県ほど低学歴・高失業)。県別の滑らかな揺らぎ(sin/cos)で
# 横断的なばらつきを付与。値は全国平均周辺の妥当域(edu_low≈0.55–0.70, unemp_rate≈0.04–0.06)。
# - edu_low : 低学歴割合(高卒以下など)の代理。Fukuda 2007 / Ping&Oshio 2023 ほか
# - unemp_rate: 失業率。Fukuda 2007(地域剥奪指標の失業領域)
# これらは震災の影響を受けない 2010年断面の属性として固定する(post 値は mediator のため不使用)。
covar_2010 <- covar_2010 %>%
dplyr::mutate(
.z_inc = (gdp_pc - mean(gdp_pc)) / stats::sd(gdp_pc),
.pref_int = as.integer(prefecture_code),
edu_low = round(0.62 - 0.030 * .z_inc + 0.012 * sin(.pref_int), 3),
unemp_rate = round(0.050 - 0.008 * .z_inc + 0.004 * cos(.pref_int), 4)
) %>%
dplyr::select(-.z_inc, -.pref_int)
# 被災3県(処置コホート)
treated_codes <- c("07", "03", "04") # 福島・岩手・宮城
prefecture_master <- prefecture_master %>%
mutate(
nested_treated = as.integer(prefecture_code %in% treated_codes),
fukushima = as.integer(prefecture_code == "07")
)
years <- 2005:2024
event_year <- 2011
```
## 年齢構造・標準人口・年齢別 treatable 死亡率
```{r}
#| label: age-structure
#| code-fold: true
#| code-summary: "年齢18階級:人口シェア・OECD2010標準人口重み・年齢別 treatable 率"
age_tbl <- tibble::tribble(
~age_band, ~pop_share, ~oecd_std, ~rate_per100k,
"0-4", 0.040, 5750, 15,
"5-9", 0.045, 5500, 3,
"10-14", 0.045, 5500, 3,
"15-19", 0.050, 5500, 6,
"20-24", 0.055, 6000, 8,
"25-29", 0.060, 6000, 10,
"30-34", 0.065, 6500, 15,
"35-39", 0.070, 7000, 22,
"40-44", 0.070, 7000, 35,
"45-49", 0.065, 7000, 55,
"50-54", 0.060, 7000, 85,
"55-59", 0.060, 6500, 130,
"60-64", 0.065, 6000, 200,
"65-69", 0.060, 5500, 320,
"70-74", 0.050, 5000, 520,
"75-79", 0.040, 4000, 850,
"80-84", 0.025, 2500, 1500,
"85+", 0.025, 1500, 3000
) %>%
mutate(
pop_share = pop_share / sum(pop_share), # 合計1に正規化
std_weight = oecd_std / sum(oecd_std), # OECD2010重み(合計1)
age_band = factor(age_band, levels = age_band)
)
# 性別シェア
sex_tbl <- tibble::tibble(sex = c("male", "female"), sex_share = c(0.49, 0.51))
```
```{r}
#| label: peek-age-structure
#| code-fold: true
#| code-summary: "2.2 で作った age_tbl / sex_tbl を確認"
# 年齢18階級の標準化テーブル(小さい参照表なので全行表示)。rate_per100k が加齢で急上昇する点に注目
age_tbl %>%
knitr::kable(
digits = 4,
caption = "age_tbl — 年齢階級別の人口シェア(pop_share)・OECD2010 標準人口重み(std_weight)・年齢別 treatable 率(rate_per100k)"
)
sex_tbl %>%
knitr::kable(caption = "sex_tbl — 性別シェア")
```
::: {.callout-note collapse="true"}
## DGP(データ生成過程)の設計(クリックで展開)
- **規模**:計算を軽くするため日本の人口を **1/10 スケール**で生成(率は不変、死亡数・人口が1/10)。
- **年齢別 treatable 率** `rate_per100k` は加齢で上昇(若年で低く、85+で高い)。OECD2010重みで標準化すると死亡頭数ベースで ASR ≈ 200/10万。これを `treatable_share` で重み付けすると平均 share(≈0.69)倍となり、**重み付き treatable ASR ≈ 140/10万**(本ノートの主アウトカム)。
- **都道府県差**:被災3県・東北をやや高め、東京・大阪をやや低めに設定。
- **全国トレンド**:treatable 死亡率は年率 −1.5%(医療進歩)で緩やかに低下。
- **都道府県×年の AR(1)残差**:同一県内・同一年で相関(→ 後述の「クラスタリングが必要」の根拠)。
- **災害ショック(処置3県のみ、event-time別)**:想定する4期(即時 / 調整 / 再均衡)に対応。
- 福島:k=0 +30%、k=1–3 +20%→+14%、k≥4 減衰。
- 宮城・岩手:k=0 +12%、k=1–3 +6%→、k≥4 0%。
- **人口(分母)**は2010年値で固定(時間変化は簡略化;率の効果に集中)。
:::
## 死亡票(個票)と分母の生成
```{r}
#| label: generate-data
#| code-fold: true
#| code-summary: "セル(県×年×年齢×性)で death を発生 → 個票へ展開 + 分母テーブル"
# --- 1. 都道府県×年 AR(1) 残差(同一県内で相関を持たせる) ------------------
gen_ar1 <- function(n, rho = 0.6, sd = 0.05) {
e <- numeric(n)
e[1] <- rnorm(1, 0, sd)
for (i in 2:n) e[i] <- rho * e[i - 1] + rnorm(1, 0, sd * sqrt(1 - rho^2))
e
}
pref_year_noise <- expand_grid(
prefecture_code = prefecture_master$prefecture_code,
year = years
) %>%
arrange(prefecture_code, year) %>%
group_by(prefecture_code) %>%
mutate(ar1 = gen_ar1(n())) %>%
ungroup()
# --- 2. 都道府県の baseline 率 乗数 ----------------------------------------
pref_mult <- prefecture_master %>%
mutate(mult = case_when(
prefecture_code == "07" ~ 1.12, # Fukushima
prefecture_code %in% c("03", "04") ~ 1.08, # Iwate, Miyagi
region == "Tohoku" ~ 1.06,
prefecture_code %in% c("13", "27") ~ 0.92, # 大都市
TRUE ~ 1.00
)) %>%
select(prefecture_code, mult)
# --- 3. セル単位の期待死亡数 → Poisson 発生 --------------------------------
SCALE <- 10 # 1/10 スケール
cells <- expand_grid(
prefecture_code = prefecture_master$prefecture_code,
year = years,
age_band = age_tbl$age_band,
sex = sex_tbl$sex
) %>%
left_join(prefecture_master, by = "prefecture_code") %>%
left_join(population_2010, by = "prefecture_code") %>%
left_join(pref_mult, by = "prefecture_code") %>%
left_join(pref_year_noise, by = c("prefecture_code", "year")) %>%
left_join(age_tbl %>% select(age_band, pop_share, rate_per100k),
by = "age_band"
) %>%
left_join(sex_tbl, by = "sex") %>%
mutate(
k = year - event_year,
# person-years(分母):2010人口 × 年齢シェア × 性シェア / スケール
person_years = pop_2010 * pop_share * sex_share / SCALE,
# 全国トレンド
trend = exp(-0.015 * (year - 2010)),
# 災害ショック(処置3県のみ)
shock = case_when(
prefecture_code == "07" & k == 0 ~ 0.30,
prefecture_code == "07" & k %in% 1:3 ~ 0.20 - 0.03 * (k - 1),
prefecture_code == "07" & k >= 4 ~ pmax(0, 0.10 - 0.015 * (k - 4)),
prefecture_code %in% c("03", "04") & k == 0 ~ 0.12,
prefecture_code %in% c("03", "04") & k %in% 1:3 ~ 0.06 - 0.015 * (k - 1),
TRUE ~ 0
),
# 男性をやや高めに(treatable はおおむね男性多め)
sex_mult = if_else(sex == "male", 1.15, 0.85),
# 期待死亡数
lambda = person_years * (rate_per100k / 1e5) * mult * sex_mult *
trend * (1 + shock) * exp(ar1),
deaths = rpois(n(), lambda)
)
# --- 4. 分母テーブル(国勢調査/推計人口に相当) ----------------------------
population_py <- cells %>%
select(prefecture_code, prefecture_name, year, age_band, sex, person_years)
# --- 5. 死亡票(個票):1行=1死亡。死因と treatable share を付与 -----------
# treatable_share は本番 bridge(appendix_A_part_A2_kantan_to_oecd_bridge.csv)に準じた
# 連続値。レンジ/グループ死因は配下 ICD の治療可能性の加重平均となり半端な小数を取る
# (肺炎 J12-J18 ≈ 5/7、敗血症 A40-A41 ≈ 0.85)。一律 50/50 の死因は 0.5 のまま。
cause_tbl <- tibble::tribble(
~icd10_cause, ~cause_label, ~treatable_share, ~prob,
"I21", "急性心筋梗塞", 0.5, 0.18,
"I60-I69", "脳血管疾患", 0.5, 0.20,
"E10-E14", "糖尿病", 0.5, 0.07,
"I10-I13", "高血圧性疾患", 0.5, 0.05,
"J12-J18", "肺炎", 0.71, 0.16,
"A40-A41", "敗血症", 0.85, 0.05,
"C18-C21", "大腸癌", 1.0, 0.10,
"C50", "乳癌", 1.0, 0.05,
"C53", "子宮頸癌", 0.5, 0.02,
"K35", "急性虫垂炎ほか外科治療可", 0.90, 0.12
) %>% mutate(prob = prob / sum(prob))
deaths_micro <- cells %>%
filter(deaths > 0) %>%
select(prefecture_code, prefecture_name, year, age_band, sex, deaths) %>%
tidyr::uncount(deaths) %>% # 1行=1死亡 に展開
mutate(
record_id = row_number(),
draw = sample(seq_len(nrow(cause_tbl)), n(),
replace = TRUE,
prob = cause_tbl$prob
),
icd10_cause = cause_tbl$icd10_cause[draw],
cause_label = cause_tbl$cause_label[draw],
treatable_share = cause_tbl$treatable_share[draw],
event_time = year - event_year,
nested_treated = as.integer(prefecture_code %in% treated_codes)
) %>%
select(
record_id, prefecture_code, prefecture_name, year, age_band, sex,
icd10_cause, cause_label, treatable_share, event_time, nested_treated
)
# --- 6. 正準オブジェクト:重み付き treatable 死亡(= Σ treatable_share) -------
# 各死亡を treatable_share で重み付けして集計。本番 sum(deaths * oecd_treatable_share)
# と同一定義(0.5 死因は 0.5 件)。これが本ノート唯一の treatable アウトカム。
cells_w <- cells %>%
select(prefecture_code, prefecture_name, year, age_band, sex, person_years) %>%
left_join(
deaths_micro %>%
count(prefecture_code, prefecture_name, year, age_band, sex,
wt = treatable_share, name = "treatable_deaths"
),
by = c("prefecture_code", "prefecture_name", "year", "age_band", "sex")
) %>%
mutate(
treatable_deaths = coalesce(treatable_deaths, 0), # 死亡ゼロセルを復元
nested_treated = as.integer(prefecture_code %in% treated_codes),
event_time = year - event_year
)
# デモ平均 treatable share(較正の目安:二値カウントとの差はおおむねこの倍率)
mean_share_demo <- with(cause_tbl, sum(prob * treatable_share))
stopifnot(
all(cells_w$treatable_deaths >= 0),
any(cells_w$treatable_deaths %% 1 != 0) # 重み付きなので分数を含む
)
cat(
"死亡票(個票)の行数:", format(nrow(deaths_micro), big.mark = ","),
" 分母テーブルのセル数:", format(nrow(population_py), big.mark = ","),
" 平均 treatable share:", round(mean_share_demo, 3), "\n"
)
```
```{r}
#| label: peek-cells-cause
#| code-fold: true
#| code-summary: "2.3 で作った cells(DGP核)と cause_tbl を head() で確認"
# (DGP核)cells:県×年×年齢×性のセル。年齢別率 rate_per100k → 期待死亡 lambda → deaths を Poisson 発生
cells %>%
select(prefecture_code, year, age_band, sex, person_years, rate_per100k, lambda, deaths) %>%
head(6) %>%
knitr::kable(digits = 2, caption = "cells(DGP核, 先頭6行)— rate_per100k × person_years から lambda を作り deaths を Poisson 実現")
# 死因 × treatable_share(連続値)の対応表
cause_tbl %>%
knitr::kable(digits = 3, caption = "cause_tbl — 死因ごとの treatable_share(治療可能割合)と出現確率 prob")
```
```{r}
#| label: peek-micro
#| code-fold: true
#| code-summary: "個票の中身(先頭10行)"
deaths_micro %>%
slice_head(n = 10) %>%
knitr::kable(caption = "死亡票(個票)の例 — 1行が1人の死亡(シミュレーション)")
```
```{r}
#| label: peek-cells-w
#| code-fold: true
#| code-summary: "2.3 の正準オブジェクト cells_w(重み付き treatable 死亡)を head() で確認"
# (解析入力)cells_w:treatable_deaths = Σ treatable_share(重み付きなので分数を取りうる)。手法A/手法B の元データ
cells_w %>%
mutate(person_years = round(person_years), treatable_deaths = round(treatable_deaths, 2)) %>%
select(prefecture_code, year, age_band, sex, person_years, treatable_deaths, nested_treated, event_time) %>%
head(6) %>%
knitr::kable(caption = "cells_w(解析入力, 先頭6行)— treatable_deaths は重み付き(Σ share)で分数を含む")
```
::: {.callout-note}
## treatable share の重み付け
個票では各死亡に **treatable_share**(死因ごとの治療可能割合)が付きます。本ノートの**主解析**は、各死亡を treatable_share で**重み付け**して数えます(**0.5 死因は 0.5 件**、1.0 死因は 1 件)。すなわち treatable 死亡数 = `sum(treatable_share)` で、これは本番パイプラインの `sum(deaths × oecd_treatable_share)`(`oecd_treatable_share > 0` で絞った後)と一致します。二値化(share>0 を一律 1 件)は採用しません。share の値は本番ブリッジに準じた連続値(例:肺炎 ≈0.71、敗血症 ≈0.85)を用います。
**なぜ 0.5 でなく 0.71 になるのか** — `0.5` と連続値は**別レベルの数字**です。
- **0.5(死因"内"の折半)**:OECD が**1つの死因**を「予防・治療のどちらも効く」と判断したときの 50/50 配分(例:急性心筋梗塞 I21・脳血管疾患・糖尿病)。死因の中での配分。
- **連続値(死因"群"の加重平均)**:日本の公表死亡は **簡単分類**(複数 ICD コードの束)で集計されるため、カテゴリの share は**配下コードの治療可能割合の加重平均**になり、半端な小数を取る。カテゴリの中での構成比。
肺炎カテゴリ(簡単分類 10200)は `J12–J18` の**7コードの束**で、
| ICD-10 | 内容 | 分類 |
|---|---|---|
| J13 | 肺炎球菌性肺炎 | **Treatable** |
| J14 | インフルエンザ菌性肺炎 | **Treatable** |
| J15 | 細菌性肺炎(他に分類されないもの) | **Treatable** |
| J16 | その他の感染性肺炎(他に分類されないもの) | **Treatable** |
| J18 | 病原体不明の肺炎 | **Treatable** |
| J12 | ウイルス性肺炎(ワクチンで一次予防) | Preventable |
| J17 | 他疾患に分類される肺炎(二次分類コード) | Preventable |
→ **7コード中 5つが treatable** なので `treatable_share = 5/7 ≈ 0.71`(残り 2/7 ≈ 0.29 は preventable)。同様に敗血症(A40–A41)は配下14コード(4桁)の加重平均で ≈0.85。つまり **0.71 は「肺炎が71%治療可能」ではなく、「このカテゴリの死亡のうち治療可能扱いのコードが 5/7 を占める」という構成比**です(分類は本研究ブリッジの符号化に基づく)。
:::
---
# 記述:被災3県と全国の treatable ASR トレンド
> **この章のねらい** — 推定の前にデータを目で点検する。処置前(2010年以前)に被災3県と対照群が平行に動いているか、2011年以降に乖離が生じているかを確認し、DID が想定する構造がデータに入っているかを見る(推定結果ではなく入力点検)。
集計版の入力となる **都道府県×年の年齢標準化 treatable 死亡率(ASR)** を作り、まず目で見ます。
```{r}
#| label: build-asr
#| code-fold: true
#| code-summary: "セル → 都道府県×年 ASR(OECD2010 直接法標準化)"
# 年齢別(性込み)に集計 → 直接法標準化
asr_panel <- cells_w %>%
group_by(prefecture_code, prefecture_name, year, age_band) %>%
summarise(
treatable_deaths = sum(treatable_deaths), person_years = sum(person_years),
.groups = "drop"
) %>%
left_join(age_tbl %>% select(age_band, std_weight), by = "age_band") %>%
mutate(rate = treatable_deaths / person_years * 1e5) %>% # 年齢別 重み付き crude rate
group_by(prefecture_code, prefecture_name, year) %>%
summarise(asr = sum(std_weight * rate), .groups = "drop") %>% # ASR
left_join(prefecture_master %>% select(prefecture_code, nested_treated, fukushima),
by = "prefecture_code"
) %>%
left_join(population_2010, by = "prefecture_code") %>%
mutate(event_time = year - event_year)
# baseline(処置3県・処置前)ASR — 後で「%効果」へ変換するのに使う
baseline_asr_treated <- asr_panel %>%
filter(nested_treated == 1, year %in% 2005:2010) %>%
summarise(m = weighted.mean(asr, pop_2010)) %>%
pull(m)
```
```{r}
#| label: fig-trends
#| fig-cap: "treatable ASR のトレンド(デモデータ)。赤破線=処置(2011)。被災3県は2011年以降に上昇。実データの主張ではない。"
#| code-fold: true
plot_dat <- asr_panel %>%
mutate(grp = case_when(
prefecture_code == "07" ~ "Fukushima",
prefecture_code == "03" ~ "Iwate",
prefecture_code == "04" ~ "Miyagi",
TRUE ~ "44 controls (mean)"
))
ctrl_mean <- plot_dat %>%
filter(nested_treated == 0) %>%
group_by(year) %>%
summarise(asr = weighted.mean(asr, pop_2010), .groups = "drop") %>%
mutate(grp = "44 controls (mean)")
bind_rows(
plot_dat %>% filter(nested_treated == 1) %>% select(year, asr, grp),
ctrl_mean
) %>%
ggplot(aes(year, asr, color = grp)) +
geom_vline(xintercept = 2010.5, linetype = "dashed", color = "red") +
geom_line(linewidth = 1) +
geom_point(size = 1.5) +
scale_color_brewer(palette = "Set1") +
labs(
title = "全年齢 treatable death の年齢標準化死亡率(ASR)",
subtitle = "被災3県 vs 44県平均(シミュレーション)",
x = "年", y = "ASR(人口10万・年あたり)", color = NULL,
caption = "デモデータ。実際の死亡率を示すものではない。"
)
```
**読み方**:処置前(2010年以前)は3県と対照群が**ほぼ平行**に推移し、2011年以降に被災3県(特に福島)が上方に乖離していれば、DID/event-study が想定する構造がデータに入っています。これは推定結果ではなく、入力データの確認です。
---
# 手法A:個人データを「そのまま」使う Event DID
> **この章のねらい** — 個人データ(死亡個票 + 人年)を**集計(率に丸める)せずそのまま**使い、Poisson 率モデル/線形確率モデルで event-study DID を推定する。年齢・性は回帰で調整。後段の手法B(集計)と対比する「個人側」の実装。
## 手法A が使うデータ形式
手法A は **個人粒度を保ったまま**、次の3つの形を使います(率に丸めません)。
- **(分子)死亡個票 `deaths_micro`** — 1行 = 1死亡。死因(ICD-10)と `treatable_share`(治療可能割合)を保持。
- **(分母)人年 `population_py`** — 都道府県 × 年 × 年齢階級 × 性 の **person-years(人年)**。
- **(解析入力)person-time セル `cells_w`** — 上の分子・分母を**同一共変量セル**で突き合わせた表。1行 = 県 × 年 × 年齢階級 × 性、`treatable_deaths`(= Σ treatable_share、**分子**)+ `person_years`(**分母**)。これがそのまま `fepois`/`feols` に渡ります(`offset = log(person_years)`、応答 = `treatable_deaths`)。
```{r}
#| label: peek-method-a-data
#| code-fold: true
#| code-summary: "手法A の入力データを head() で確認"
# (分子)死亡個票:1行 = 1死亡
deaths_micro %>%
select(record_id, prefecture_code, year, age_band, sex,
icd10_cause, cause_label, treatable_share, event_time, nested_treated) %>%
head(6) %>%
knitr::kable(caption = "deaths_micro(分子・個票)— 1行 = 1死亡。死因と treatable_share を保持")
# (分母)人年:都道府県 × 年 × 年齢階級 × 性
population_py %>%
mutate(person_years = round(person_years)) %>%
head(6) %>%
knitr::kable(caption = "population_py(分母)— person-years(人年)")
# (解析入力)person-time セル:分子 + 分母を同一セルで突き合わせ
cells_w %>%
mutate(person_years = round(person_years),
treatable_deaths = round(treatable_deaths, 2)) %>%
select(prefecture_code, year, age_band, sex,
person_years, treatable_deaths, nested_treated, event_time) %>%
head(6) %>%
knitr::kable(caption = "cells_w(解析入力)— 1行 = 県×年×年齢×性。treatable_deaths = Σ share(分子), person_years = 分母")
```
## 個人データそのまま、とは(person-time モデル)
::: {.callout-note}
## 「1行=1人」と「セル集計」は同じ尤度(lossless)
**率(rate;人年あたりの死亡率)**を推定するには **分子(死亡=個票)+ 分母(人年, person-time)** が要ります。死亡票を分母と突き合わせると、各個人の「その年の treatable 寄与(死ななければ 0、treatable 死亡なら死因の **treatable_share**)」という **person-year(人年)レコード**が得られます。
- これを **1行=1人年**に展開した線形確率モデル(LPM)と、
- **同じ共変量を持つ個人をセルにまとめ**、重み付き treatable 死亡数(Σ share)を応答・人年を offset とした Poisson 回帰
は **数学的に同一の尤度**を最適化します(共変量が同一のレコードをまとめても情報は失われない)。死亡が稀(treatable ≈ 0.2%)なため、数百万行に展開する代わりに **セル形式の Poisson** を使うのが実務的で、これは「集計データ(率に丸めたもの)」とは異なり、**個人の年齢・性・死因情報をそのまま保持**しています。
**「率(rate)」か「割合(risk / proportion)」か**:本ノートの outcome は **率(rate)** です — 分母が **人年(person-time)**、Poisson の `offset = log(person_years)` で「人年あたりの死亡」を推定し、単位は「人口10万・年あたり」(ASR も率)。一方 **割合 / リスク** は分母が **人数**の無次元量(0–1)で、本来は別概念です。ただし**観察が1年単位 ＋ イベントが稀**(treatable ≈ 0.2%)なので、**1人年あたりの率 ≈ 1年あたりのリスク**で数値的にほぼ一致し、「1行=1人年・died 0/1 の平均」も近似的に率として読めます。
:::
```{r}
#| label: uncount-proof
#| code-fold: true
#| code-summary: "「1行=1人」展開と セル率が一致することの確認(1セル)"
one_cell <- cells_w %>% filter(prefecture_code == "07", year == 2011, age_band == "70-74", sex == "male")
stopifnot(nrow(one_cell) == 1)
# このセルに属する個票(treatable 死亡)の share を取り出す
cell_deaths <- deaths_micro %>%
filter(prefecture_code == "07", year == 2011, age_band == "70-74", sex == "male")
n_py <- as.integer(round(one_cell$person_years))
stopifnot(n_py >= nrow(cell_deaths))
# person-year 展開:死亡者は重み = treatable_share、非死亡者は 0
person_weight <- c(cell_deaths$treatable_share, rep(0, n_py - nrow(cell_deaths)))
tibble(
`セルの treatable_deaths/person_years` = one_cell$treatable_deaths / one_cell$person_years,
`1行=1人 展開の mean(weight)` = mean(person_weight)
) %>% knitr::kable(
digits = 6,
caption = "個人の treatable 重み(0 または share)の平均 = セルの重み付き率(福島・2011・70-74歳・男性)"
)
```
## 個人レベル event-study の推定(fixest)
ここでは同じ個人 person-time データに、**2つの等価な指定(A1・A2)**を当てはめます。どちらも DID 構造(県 + 年 + 年齢階級 FE、`event-time × 処置` 交互作用、県クラスタ SE)は同じで、**効果を測る「ものさし(スケール)」だけが違います**。
- **A1:Poisson 率モデル**(`fepois`)— **log(率)スケール**。効果は**率比** $\exp(\beta_k)$(乗法的)。本ノートの主推定。
- **A2:線形(率)モデル(LPM)**(`feols`)— **加法スケール**。効果は**率の差**(人年あたり)。A1 の対照として併記。
A1 と A2 は表現が違うだけで同じ DID を見ており、**%効果に直すとほぼ重なります**(後段の図・「ブリッジ」で確認)。以下、共通の設計を整理してから両者の式を示します。
**なぜ Poisson(率モデル)を主推定にするか**
- **outcome が率(count ÷ person-time)**:分子=死亡数、分母=人年。`offset = log(PY)` 付き Poisson が「人年あたりの死亡」を表す自然な尤度。
- **ゼロセルに頑健**:稀イベント(treatable ≈ 0.2%)で死亡 0 のセルが多発。Poisson は応答が 0 でも問題なし。
- **係数がそのまま率比**:$\exp(\beta_k)$ = 率比で解釈が直接的。死亡率は乗法的に動くことが多く log スケールが自然。
- **QMLE で分布仮定に非依存**:`fepois` は真に Poisson でなくても条件付き平均さえ正しければ一致推定(過分散・系列相関は県クラスタ SE で対処)。
**他の選択肢との比較**
- **線形(率)LPM(= A2)**:加法スケール、estimand は「率の差」。対照として併記し後段ブリッジで % 効果はほぼ一致。差の解釈には有用だが、ゼロ近傍で負の予測率を許す。
- **log(率)を OLS**:ゼロセルで $\log 0$ が定義不能。さらに対数変換 → 逆変換のバイアス(retransformation bias)を被る。Poisson なら回避。
- **負の二項(NB)**:過分散に対応する反面、分散構造を誤ると不一致になりうる。Poisson QMLE + 県クラスタ SE の方が分散仮定に頼らず頑健なため不採用。
- **ロジスティック / 二項**:分母が「人数」のリスク(0–1)向け。本ノートは分母が人年の**率**なので不適(稀イベントゆえ率 ≈ リスクで数値は近いが概念が別)。
推定の構成:
- **応答**:個人の死亡(person-time)= 重み付き treatable 死亡(Σ share)
- **固定効果**:都道府県 + 年 + 年齢階級
- **共変量**:性別
- **処置効果**:`i()` による **event-time × 処置**の交互作用
- **標準誤差**:都道府県でクラスタ
**モデル式**(添字:$g$=都道府県, $t$=年, $a$=年齢階級, $s$=性別;$D$=treatable 死亡数, $PY$=人年)
**(A1) Poisson 率モデル(主)** — $\log$(率)に線形:
$$
\log \mathbb{E}[D_{gtas}]
= \underbrace{\log PY_{gtas}}_{\text{offset}}
+ \alpha_g + \tau_t + \gamma_a
+ \delta\,\mathrm{Male}_s
+ \sum_{k \neq -1} \beta_k \,\mathbf{1}\{\,t-2011=k\,\}\cdot \mathrm{Treated}_g
$$
- $\alpha_g, \tau_t, \gamma_a$:都道府県・年・年齢階級の固定効果
- $\mathrm{Treated}_g$:被災3県ダミー / $\mathbf{1}\{\cdot\}$:指示関数
- **event-time $k$ の効果は率比** $\exp(\beta_k)$(%効果 $=(\exp\beta_k-1)\times100$)
- 基準期 $k=-1$ で $\beta_{-1}=0$、$\mathrm{ATT}_k=\beta_k$
**(A2) 線形(率)モデル LPM** — 加法スケール、(A1)と同情報:
$$
\frac{D_{gtas}}{PY_{gtas}}
= \alpha_g + \tau_t + \gamma_a
+ \delta\,\mathrm{Male}_s
+ \sum_{k \neq -1} \beta_k \,\mathbf{1}\{\,t-2011=k\,\}\cdot \mathrm{Treated}_g
+ \varepsilon_{gtas}
$$
- 人年で加重($\text{weights}=PY$)
- $\beta_k$ は**率の加法的変化**(人年あたり)
- SE は県クラスタ
**コードとの対応**:
- `i(event_time, nested_treated, ref = -1)` $= \sum_{k\neq-1}\beta_k\,\mathbf{1}\{k\}\,\mathrm{Treated}_g$
- `| prefecture_code + year + age_band` $= \alpha_g+\tau_t+\gamma_a$
- `+ sex` $= \delta\,\mathrm{Male}_s$
- `offset = ~log(person_years)` $= \log PY_{gtas}$
- `cluster = ~prefecture_code` = 県クラスタ
::: {.callout-note}
## 「Poisson で DID」— 乗法(log率)スケールの event-study DID
A1 は **Poisson 回帰で DID(event-study)** を行っています。仕組みは2元固定効果 DID と同じで、**処置効果は `event-time × 処置` の交互作用 $\beta_k$** から取り出します(県 FE $\alpha_g$ が「群」、年 FE $\tau_t$ が「時間」の主効果を吸収)。Poisson はそれを **log(率)スケール**に乗せているだけです。
- **何が DID か**:$\beta_k$ は処置群と対照群の「前後変化の差」。ただし log スケールなので、実体は **率比の比(ratio of rate ratios)= 乗法的 DID**。$\exp(\beta_k)$ が率比、$(\exp\beta_k-1)\times100$ が %。
- **平行トレンドの前提が違う**:Poisson DID は **log 率での平行トレンド(処置がなければ群間の率比が一定 = 比例的)** を仮定。一方 LPM(A2)や ASR の CS-DID(手法B)は **水準(率の差)での平行トレンド(加法的)** を仮定。両者は別の仮定で、どちらが自然かは outcome 次第(死亡率は比例的に動くことが多く、log スケールがしばしば自然)。
- **この設計では素直**:処置時点が3県とも2011年で共通(非段階)なので、TWFE event-study の負の重み問題(Goodman-Bacon / Sun-Abraham)は生じません。
- **`fepois` は擬似最尤(QMLE)**:分布が真に Poisson でなくても、条件付き平均が正しければ率比を一致推定します(過分散・系列相関は県クラスタ SE で対処)。
- **なぜ手法B と % が近いか**:効果が中程度なら「加法 %(= ATT/baseline)」≈「乗法 %(= $\exp\beta_k-1$)」で数値が接近(後述ブリッジ)。ただし estimand は別物(差の差 vs 比の比)。
:::
::: {.callout-important}
## 固定効果は「都道府県+年+年齢階級」— **個人固定効果は使えない/不要**
同一個人を追跡していない(反復横断)ため個人 FE は定義できませんが、**不要**です。処置は都道府県単位なので都道府県 FE が交絡(地域の恒常差)を、年 FE が全国共通ショックを吸収します。年齢階級 FE と性別共変量が、毎年入れ替わる人口の**構成変化**を調整します(= 回帰版の年齢標準化)。
:::
```{r}
#| label: model-A
#| code-fold: true
#| code-summary: "手法A:個人レベル event-study(Poisson rate モデル + LPM)"
# 解析用セル(重み付き treatable 死亡。nested_treated/event_time は生成時に付与済)
cell_an <- cells_w %>% filter(person_years > 0)
# (A1) Poisson 率モデル:重み付き treatable 死亡(分数)を log-rate で。offset = log(人年)
mA_pois <- fepois(
treatable_deaths ~ i(event_time, nested_treated, ref = -1) + sex |
prefecture_code + year + age_band,
offset = ~ log(person_years),
cluster = ~prefecture_code,
data = cell_an
)
# (A2) 線形確率モデル(LPM)等価:応答=重み付き率, 重み=人年。1行=1人 LPM と同値
cell_an <- cell_an %>% mutate(rate = treatable_deaths / person_years)
mA_lpm <- feols(
rate ~ i(event_time, nested_treated, ref = -1) + sex |
prefecture_code + year + age_band,
weights = ~person_years,
cluster = ~prefecture_code,
data = cell_an
)
# event-study 係数の抽出ヘルパー(fixest の i() 項を (k, est, se) へ)
extract_es <- function(model, label, multiplicative = FALSE, base_rate = NULL) {
out <- broom::tidy(model) %>%
filter(str_detect(term, "^event_time::-?\\d+:nested_treated$")) %>%
mutate(
k = as.integer(str_match(term, "event_time::(-?\\d+):")[, 2]),
model = label
) %>%
{
if (multiplicative) { # Poisson: log-rate → %
mutate(.,
pct = (exp(estimate) - 1) * 100,
lo = (exp(estimate - 1.96 * std.error) - 1) * 100,
hi = (exp(estimate + 1.96 * std.error) - 1) * 100
)
} else { # 加法(率): / baseline → %
mutate(.,
pct = estimate / base_rate * 100,
lo = (estimate - 1.96 * std.error) / base_rate * 100,
hi = (estimate + 1.96 * std.error) / base_rate * 100
)
}
} %>%
select(model, k, estimate, std.error, pct, lo, hi) %>%
arrange(k)
if (nrow(out) == 0) {
stop("Event-study coefficients were not found. Check the fixest i() term names.")
}
out
}
base_rate_treated <- cell_an %>%
filter(nested_treated == 1, year %in% 2005:2010) %>%
summarise(r = sum(treatable_deaths) / sum(person_years) * 1e5) %>%
pull(r) / 1e5
esA_pois <- extract_es(mA_pois, "A1: 個人 Poisson", multiplicative = TRUE)
esA_lpm <- extract_es(mA_lpm, "A2: 個人 LPM",
multiplicative = FALSE,
base_rate = base_rate_treated
)
knitr::kable(esA_pois %>% select(k, `log-rate係数` = estimate, SE = std.error, `効果%` = pct),
digits = 3, caption = "手法A(個人 Poisson)event-study 係数。k=0 で約 +log(1.x)。"
)
```
**読み方(手法A)**:
- Poisson の係数は **log 率比**。`効果% = (exp(係数) − 1) × 100`。DGP の真値は **k=0 で福島+30%・宮城/岩手+12% の3県平均**なので、コホート ATT_0 はおおむね **+20%前後**を回収するはずです。
- $k \le -2$ の係数が **0 近傍**であれば平行トレンド前提と整合(基準は $k=-1$ なので必ず 0)。
- $k=0$ で正に跳ね、$k$ が増えると**減衰**していれば、「即時ショック → 調整 → 再均衡」の動学を捉えています。
最後に、手法A の event-study 係数(A1 Poisson と A2 LPM)を **%効果**で可視化します。
```{r}
#| label: fig-method-a
#| fig-cap: "手法A の個人レベル event-study(%効果)。A1=Poisson(率比ベース)、A2=LPM(率の加法ベース)。点=ATT_k 推定、帯=95%CI、縦破線=基準 k=-1(2010年)、横線=効果ゼロ。デモデータ。"
#| code-fold: true
bind_rows(esA_pois, esA_lpm) %>%
filter(k >= -6, k <= 10) %>%
ggplot(aes(k, pct, color = model, fill = model)) +
geom_hline(yintercept = 0, color = "grey50") +
geom_vline(xintercept = -1, linetype = "dashed", color = "red") +
geom_ribbon(aes(ymin = lo, ymax = hi), alpha = 0.12, color = NA) +
geom_line(linewidth = 0.8) +
geom_point(size = 1.8) +
scale_x_continuous(breaks = seq(-6, 10, 2)) +
scale_color_brewer(palette = "Set1") +
scale_fill_brewer(palette = "Set1") +
labs(
title = "手法A:個人レベル event-study(被災3県 vs 44県)",
subtitle = "全年齢 treatable death への %効果。基準 k = -1(2010年)。",
x = "event time k(= 年 - 2011)", y = "treatable death への効果(%)",
color = NULL, fill = NULL,
caption = "デモデータ。実際の効果ではない。"
) +
theme(legend.position = "top")
```
**図の読み方**:点が ATT_k(%)、帯が 95%CI。基準 $k=-1$ で 0、処置前($k<0$)が **0 近傍**なら平行トレンドと整合、$k=0$ で立ち上がり以降**減衰**していれば、手法A 単独で「即時 → 調整 → 再均衡」の動学を捉えています。A1(Poisson, 率比)と A2(LPM, 率差)は尺度が異なりますが、%効果ではほぼ重なります(後述「ブリッジ」で手法B とも重ね描きします)。
## ベースライン SES 共変量による調整(事前指定の感度)
> **この節のねらい** — 主解析(無調整の手法A)に対し、**2010年(震災前)のベースライン社会経済(SES)共変量**を加えた**事前指定の感度仕様**を示す。ここでは「**どの変数を**」「**どの方法で**」入れるかを、既存文献と本デザインの制約に基づいて明示する。
### 調整変数の選択(なぜこの4変数か)
採用するのは **2010年県レベル baseline の4変数**:**所得**(`gdp_pc`)・**学歴**(`edu_low`)・**失業**(`unemp_rate`)・**人口密度**(`log_density`)。
- **領域マッピング**:所得・学歴・失業は、日本の**地域剥奪指標**の中核領域であり、県・市区町村の **<74歳死亡率**と相関する(Fukuda 2007, *Public Health*; 7領域の z スコア指標で男 r=0.65)。人口密度は剥奪と**独立**に全死亡を予測する地域要因(Nakaya 2014, *PLoS One*; 個人 SES 調整後も有意)。所得は age-adjusted 県死亡率と関連(Fukuda 2007, *Biosci Trends*)。学歴は日本の自己評価健康・健康格差で頑健な決定因(Ping & Oshio 2023; Takahashi 2020; Furuya 2015)。災害下の SES 格差の検討でも教育・所得が指標(Kino 2023)。
- **baseline 限定の理由**:震災**後**の所得・失業・人口密度は、災害→経済混乱・避難の**下流(中間変数/合流点)**である。後値で調整すると **ATT を null 方向に歪める過調整**になる(Schisterman 2009; VanderWeele 2019)。よって**全て 2010年値で固定**する。
- **入れない変数の明示**:(i) 高齢化率(`elderly_share`)は手法A では **`age_band` 固定効果が年齢を完全調整**するため、**二重調整を避けて追加しない**。(ii) アウトカム予測子は分散を下げバイアスを増やさず投入価値がある(Brookhart 2006)ため、この4変数は妥当。
- **parsimony(自由度)**:処置はわずか**3県**。共変量×event-time 項は自由度を食い、pre-trend を**過剰適合**し得る。よって**事前指定の最小限(≤4–5)**に限定する。
### 調整方法の選択(なぜ Poisson回帰×交互作用で、PS ではないか)
- **PS-IPW を使わない理由**:処置は**県単位で決定論的**(47県中3県、県内分散ゼロ)。処置割当の傾向スコアは 0/1 に張り付き、**overlap/正値性が無く退化**する。集計レベルの「重み付け」の正解は PS ではなく **synthetic control**(ドナーをこれら baseline 予測子で釣り合わせる)であり、本研究では手法B/SCM がその役割を担う。
- **回帰調整を選ぶ理由**:年齢・性は `age_band` FE + `sex` で**回帰標準化**済み(estimand=ASR と整合、g-computation で外部標準人口へ周辺化可)。死亡は**全数**で survey weight 不要 → 反復横断サーベイ向けの PS×survey-weight 法(Ye 2025 ほか)の重み半分は**本解析に不適用**。
- **×event-time 交互作用にする技術理由**:時間不変の baseline 変数を `| prefecture_code` 入りモデルに**加法主効果**で入れると、**県 FE に吸収されて消える**。条件付き平行トレンド(トレンドが baseline SES で異なる)を効かせるには、**baseline × event-time** 交互作用 `i(event_time, z_X, ref = -1)` として投入する必要がある。
- **DR/手法B との役割分担**:共変量条件付き平行トレンドの**主たる処理は手法B**(CS-DID, `xformla`, 二重頑健)。手法A の SES 調整は**事前指定の感度**であり primary ではない。
- **PS が正解になる場面の明示**:所得・学歴・失業・独居を**個人属性**として持つ **demand-side 個票**(国民生活基礎調査/クレーム)を反復横断で母集団 ATT 推定するときのみ(群メンバーシップに overlap あり)→ Ye 法 or DR-DiD。**死亡個票には不適用**。
```{r}
#| label: model-A-adj
#| code-fold: true
#| code-summary: "手法A(感度):baseline 2010 SES を event-time 交互作用で調整した Poisson"
# baseline SES を県レベルで z 標準化(係数解釈と fepois 収束の安定化)
ses_z <- covar_2010 %>%
transmute(
prefecture_code,
z_income = as.numeric(scale(gdp_pc)),
z_edu = as.numeric(scale(edu_low)),
z_unemp = as.numeric(scale(unemp_rate)),
z_density = as.numeric(scale(log(pop_density)))
)
cell_an_adj <- cell_an %>% left_join(ses_z, by = "prefecture_code")
# baseline SES × event-time 交互作用(主効果は県 FE が吸収するので交互作用で投入)
mA_pois_adj <- fepois(
treatable_deaths ~ i(event_time, nested_treated, ref = -1) + sex
+ i(event_time, z_income, ref = -1) # baseline 所得 × event-time
+ i(event_time, z_edu, ref = -1) # baseline 学歴 × event-time
+ i(event_time, z_unemp, ref = -1) # baseline 失業 × event-time
+ i(event_time, z_density, ref = -1) | # baseline 人口密度 × event-time
prefecture_code + year + age_band,
offset = ~ log(person_years),
cluster = ~prefecture_code,
data = cell_an_adj
)
# ATT_k は nested_treated 交互作用のみを抽出(extract_es は ^...:nested_treated$ で限定)
esA_pois_adj <- extract_es(mA_pois_adj, "A1-adj: 個人 Poisson + baseline SES",
multiplicative = TRUE
)
# 無調整(primary)と調整(感度)の ATT_k を並べる
compare_adj <- esA_pois %>%
select(k, `無調整 %` = pct) %>%
left_join(esA_pois_adj %>% select(k, `SES調整 %` = pct), by = "k") %>%
mutate(`差(pt)` = `SES調整 %` - `無調整 %`)
knitr::kable(compare_adj,
digits = 2,
caption = "手法A:無調整(primary)vs baseline SES 調整(感度)の ATT_k(%効果)。baseline は処置効果と独立に生成したため、両者は概ね近いはず(大きく乖離する場合は join/標準化のバグを疑う)。"
)
```
::: {.callout-warning}
## 解釈上の注意 — 3処置県では SES×event-time は弱識別
`z_X × event_time` 項は **44 非処置県の横断的な baseline SES 差**から識別されますが、処置が **3県**しかないため、`nested_treated × event_time`(=ATT)と baseline SES 項が**部分的に共線**になりやすく、`fepois` が一部の項を落とす/係数が不安定になることがあります。これは「調整変数の選択」で述べた **parsimony(≤4–5)** の根拠そのものです。識別が弱い場合の**節約的代替**は、`i(event_time, z_X)` を **`z_X : post`(post 期一括の1係数)** に置き換える形です。本ノートでは event-study 準拠の ×event-time を既定とし、SES 調整は**感度**として位置づけます。
:::
## FAQ
::: {.callout-note}
## なぜ「率」ではなくこの形を渡すのか
手法B が渡すのは県×年あたり **1つの率(ASR)** に丸めた集計データですが、手法A は **`treatable_deaths`(分子)と `person_years`(分母)を分けたまま** GLM に渡します。こうすると年齢・性・死因の情報が保持され、年齢調整(年齢階級 FE)や感度分析を回帰の中で柔軟に行えます。死亡が稀なため、同一共変量の個人をセルにまとめても情報は失われません(`1行=1人年` の展開と同じ尤度)。
**(付記)Poisson で ASR を扱えるか**:ASR(直接法標準化率)は外部の標準人口重みで**合成済みの率**なので、そのまま Poisson の応答(カウント)には**できません**。ただし **Poisson に年齢を入れて推定 → 予測した年齢別率を標準人口重み $w_a$ で加重平均(回帰標準化 / g-computation)** すれば、直接法標準化した率・率比(= ASR ベースの効果)を Poisson から**再現**できます。手法A の年齢階級 FE つき Poisson は、年齢別率比を**人年(内部)重み**で平均したもので、効果が年齢で一様なら直接法標準化と一致します(異質なら外部重みでの周辺化が必要;前述「年齢標準化が要らないのでは?」と同じ論点)。なお Poisson + offset + 年齢調整は **間接標準化(SMR)** の標準モデルでもあります。
:::
::: {.callout-important}
## Fukushima Study では ASR と「死亡数 + 年齢調整」をどう選ぶか
Fukushima Study の主結果として読者に示す estimand は、**全年齢 extended-treatable mortality の OECD 2010 ASR** がもっとも明確です。ASR は「もし全県・全年が同じ標準人口構成だったら、treatable death は人口10万人・年あたりどれだけ増減したか」という単一の率なので、県間・年次・既存の OECD avoidable mortality 指標との比較がしやすいからです。
一方で、実装上は「死亡数そのもの」を outcome にするのではなく、**treatable 死亡数(分子) + 人年(分母) + 年齢調整**として扱います。つまり本セクションの A1 は、`treatable_deaths` を応答、`offset = log(person_years)` を分母、`age_band` を固定効果として入れる **死亡率モデル**です。これは粗死亡数の DID ではなく、年齢構成の変化を回帰内で調整した person-time DID です。
選び方は次のように整理します。
| 目的 | 主に使う表現 |
|---|---|
| 論文の主結果として「何/10万人・年」の効果を示す | **ASR** |
| 個票情報を保持し、年齢・性別・死因別の調整や異質性を扱う | **死亡数 + 人年 offset の Poisson** |
| ASR と完全に同じ標準人口重みの estimand を個人モデルから出す | 年齢別 ATT を推定し、OECD 2010 重みで **g-computation** |
| 実際の追加死亡数(excess deaths)を述べる | 率モデルから予測死亡数へ変換する副次解析 |
したがって本ノートの実務的な位置づけは、**主 estimand は ASR、推定の強い実装・検証は死亡数 + 人年 + 年齢調整モデル**です。年齢効果がほぼ一様なら A1/A2 と ASR ベースの手法Bは近い値になります。年齢別効果が大きく異質なら、年齢階級 FE だけでは観測人年で平均した効果になり、ASR と同じ外部標準人口重みにするには、年齢別予測率を OECD 2010 重みで周辺化します。
:::
::: {.callout-note collapse="true"}
## なぜ DID 専用パッケージ(`did` 等)でなく `fixest` か
手法A で `fixest`(`fepois`/`feols`)、手法B で `did`(Callaway–Sant'Anna)と**使い分けて**います。手法A を fixest にする理由:
- **処置が非段階(2011年で共通)**:Callaway–Sant'Anna / Sun–Abraham / Borusyak / de Chaisemartin–D'Haultfœuille らの DID 専用推定量は、**処置時点がバラバラ(staggered)で効果が異質なときに TWFE が負の重みで歪む問題**を直すためのもの。本研究は3県とも2011年処置 + クリーンな never-treated(44県)なので、**素の TWFE event-study `i(event_time, treated, ref=-1)` が既に偏りのない2群 DID に一致**します(専用パッケージが解く問題が起きない)。
- **個人 person-time の率モデルが必要**:手法A は **Poisson 率モデル(`offset = log(人年)`)+ 年齢・性の調整 + 県クラスタ SE** を、個人/セル粒度・高次元 FE で回します。これは `fixest::fepois` がネイティブ&高速に扱える一方、`did::att_gt` は**集計パネルの連続アウトカム向けの 2×2 二重頑健/IPW/回帰推定**が中心で、率の offset・個人レベル GLM を自然には扱いません。→ 個人側=fixest、集計 ASR 側=`did` の役割分担。
- **計算量**:大規模な個票/セル + 高次元 FE では fixest が圧倒的に速い。
- **fixest にも DID-robust はある**:必要なら fixest 内の `sunab()`(Sun–Abraham 2021)で staggered にも対応可能。「fixest を使う = 現代 DID を無視」ではなく、**今回は非段階なので不要**なだけ。
**逆に専用パッケージが要る場面**:処置が staggered で効果が異質/共変量条件付き平行トレンドを二重頑健で/uniform(sup-t)band が欲しい — これらの理由で **手法B は `did::att_gt` を採用**しています。
:::
---
# 手法B:集計データに「加工」して使う(Callaway–Sant'Anna DID)
> **この章のねらい** — 同じデータを**都道府県×年の年齢標準化率(ASR)に集計**し、Callaway–Sant'Anna の DID(`did::att_gt`)で ATT_k を推定する。事前標準化 + 共変量の回帰調整 + 少数クラスタ band を備えた「集計側」の実装で、手法A と突き合わせる。
## 主解析(primary)の推定量
::: {.callout-note}
## 用語:Callaway–Sant'Anna (CS) DID と回帰調整(RA)
- `did::att_gt()` は **group-time ATT** を推定し、`aggte(type="dynamic")` で event-time 別 ATT_k に集約します。
- **never-treated(44県)**を比較群に、**回帰調整(`est_method = "reg"`;Sant'Anna–Zhao 2020)**で **2010年(震災前)の文献ベース baseline SES(所得・学歴・失業・人口密度)+ baseline ASR** を調整します。共変量セットは手法A(感度)と**同一**に揃えています。高齢化率は ASR の年齢標準化(および手法A の `age_band` FE)と重複するため**除外**(意図的選択)。
- 単位は**都道府県のバランスパネル**(`panel = TRUE`)。応答は**事前に年齢標準化した ASR**(=「集計データに加工」)。
:::
```{r}
#| label: model-B
#| code-fold: true
#| code-summary: "手法B:都道府県×年 ASR に CS-DID(att_gt → aggte dynamic)"
did_panel <- asr_panel %>%
left_join(covar_2010, by = "prefecture_code") %>%
left_join(
asr_panel %>% filter(year %in% 2005:2010) %>%
group_by(prefecture_code) %>%
summarise(baseline_asr = mean(asr), .groups = "drop"),
by = "prefecture_code"
) %>%
mutate(
unit_id = as.integer(prefecture_code),
first_treat = if_else(nested_treated == 1, event_year, 0L),
log_density = log(pop_density)
)
set.seed(20260604)
att_obj <- did::att_gt(
yname = "asr",
tname = "year",
idname = "unit_id",
gname = "first_treat",
xformla = ~ gdp_pc + edu_low + unemp_rate + log_density + baseline_asr,
data = as.data.frame(did_panel),
panel = TRUE,
control_group = "nevertreated",
est_method = "reg",
weightsname = "pop_2010",
base_period = "universal",
bstrap = TRUE, cband = TRUE, biters = 499
)
agg_dyn <- did::aggte(att_obj, type = "dynamic", na.rm = TRUE)
esB <- tibble(
model = "B: 集計 CS-DID",
k = agg_dyn$egt,
estimate = agg_dyn$att.egt,
std.error = agg_dyn$se.egt,
pct = agg_dyn$att.egt / baseline_asr_treated * 100,
lo = (agg_dyn$att.egt - 1.96 * agg_dyn$se.egt) / baseline_asr_treated * 100,
hi = (agg_dyn$att.egt + 1.96 * agg_dyn$se.egt) / baseline_asr_treated * 100
)
knitr::kable(esB %>% select(k, `ATT(ASR)` = estimate, SE = std.error, `効果%` = pct),
digits = 2, caption = "手法B(集計 CS-DID)の event-time 別 ATT_k(ASR単位)。"
)
```
**読み方(手法B)**:ATT_k は **ASR(人口10万あたり)の加法的変化**。`効果% = ATT_k / baseline_ASR × 100` で手法Aと同じ%尺度に並べられます。CS-DID は共変量の動学差を回帰調整し、never-treated を比較群に固定する点が単純 TWFE と異なります。
---
# ブリッジ:個人 ⇔ 集計はなぜ(ほぼ)一致するか
> **この章のねらい** — 手法A(個人)と手法B(集計)の ATT_k を同じ%尺度で重ね描きし、**なぜ点推定がほぼ一致するのか**(年齢調整の流儀が違うだけで同じ県×年の変動を見ている)を確認する。
```{r}
#| label: fig-bridge
#| fig-cap: "event-study の重ね描き(%尺度)。手法A(個人 Poisson/LPM)と手法B(集計 CS-DID)は同じ動学を回収する。点は推定、帯は95%CI。"
#| code-fold: true
es_all <- bind_rows(esA_pois, esA_lpm, esB) %>%
filter(k >= -6, k <= 10)
ggplot(es_all, aes(k, pct, color = model, fill = model)) +
geom_hline(yintercept = 0, color = "grey50") +
geom_vline(xintercept = -0.5, linetype = "dashed", color = "red") +
geom_ribbon(aes(ymin = lo, ymax = hi), alpha = 0.12, color = NA) +
geom_line(linewidth = 0.8) +
geom_point(size = 1.8) +
scale_color_brewer(palette = "Dark2") +
scale_fill_brewer(palette = "Dark2") +
labs(
title = "個人レベル event-study vs 集計 CS-DID(%効果)",
subtitle = "被災3県コホート vs 44県。基準 k = -1(2010年)。",
x = "event time k(= 年 - 2011)", y = "treatable death への効果(%)",
color = NULL, fill = NULL,
caption = "デモデータ。3手法が同じ ATT_k 動学を回収する。"
) +
theme(legend.position = "top")
```
```{r}
#| label: tbl-bridge
#| code-fold: true
#| code-summary: "k=0 と post 平均(k≥0)の3手法比較表"
post_summary <- es_all %>%
filter(k >= 0) %>%
group_by(model) %>%
summarise(`post平均 効果%` = mean(pct), .groups = "drop")
k0 <- es_all %>%
filter(k == 0) %>%
select(model, `k=0 効果%` = pct)
left_join(k0, post_summary, by = "model") %>%
knitr::kable(digits = 1, caption = "3手法の効果(%)。点推定はおおむね一致する。")
```
::: {.callout-tip}
## 個人 vs 集計の相違点(要点)
| 論点 | 手法A:個人データそのまま | 手法B:集計データに加工 |
|---|---|---|
| 年齢調整 | 年齢階級 FE / 共変量(**回帰版**標準化) | 事前に**直接法**年齢標準化(ASR) |
| 重み | データ件数(人年)で自然加重 | `weightsname = "pop_2010"` で人口加重 |
| 死因・年齢の保持 | **個票の情報をそのまま保持**(感度分析が容易) | 率に丸めると失われる |
| 推論の単位 | **県クラスタ**(後述) | 県(パネル単位) |
| 計算 | セル形式なら軽量・lossless | 最軽量 |
点推定が一致するのは、両者とも **「県×年」の処置×時間の変動**から同じ ATT を取り出しているからです(年齢調整の流儀が違うだけ)。
:::
---
# 推論(最重要):個人データでも「精度」は増えない
> **この章のねらい** — 点推定ではなく**不確実性(SE・p 値)**に注目する。処置が県単位=有効クラスタは3つなので、個人 n を増やしても精度は増えない。素朴 SE → 県クラスタ SE → randomization inference の順に、不確実性を正しく評価する。
::: {.callout-important}
## 処置は県単位 → 有効クラスタ数 = treated 3
個人を数百万行に増やしても、**処置は都道府県単位で割り当て**られているため、独立な「処置の繰り返し」は **3県**しかありません。個人を独立扱いした素朴(iid)な標準誤差は**過小**で、p 値を誤って小さく見せます。
:::
```{r}
#| label: se-comparison
#| code-fold: true
#| code-summary: "同じ点推定・3種類の標準誤差(iid vs 県クラスタ)"
# k=0 の係数で比較(手法A Poisson)
b0 <- coef(mA_pois)["event_time::0:nested_treated"]
se_iid <- summary(mA_pois, vcov = "iid")$se["event_time::0:nested_treated"]
se_cl <- summary(mA_pois, vcov = ~prefecture_code)$se["event_time::0:nested_treated"]
tibble(
推論 = c("素朴 iid(個人を独立扱い・誤り)", "都道府県クラスタ(正しい方向)"),
`k=0 係数(log率比)` = c(b0, b0),
SE = c(se_iid, se_cl),
`素朴比` = c(1, round(se_cl / se_iid, 1))
) %>% knitr::kable(
digits = 4,
caption = "点推定は同じでも、クラスタSEは iid より大きい(= 個人nは精度を生まない)。"
)
```
## 少数クラスタ向け randomization inference(RI)
クラスタSEも treated が3つだと不正確になりがちです。**処置の割当を入れ替えるプラセボ検定(randomization inference)**で、観測された効果が偶然どれだけ起きにくいかを評価します(本研究の主推論である studentized randomization inference を簡略化したもの)。
```{r}
#| label: randomization-inference
#| code-fold: true
#| code-summary: "47県から3県をランダムに「処置」にするプラセボ検定"
all_codes <- prefecture_master$prefecture_code
# 統計量:post-pre の単純 DID。バランスパネルなので feols(asr ~ tr:post | 県 + 年)
# の post×treated 係数と同じ情報を高速に計算できる。
ri_delta <- asr_panel %>%
group_by(prefecture_code) %>%
summarise(
delta = mean(asr[year >= event_year]) - mean(asr[year < event_year]),
.groups = "drop"
)
ri_stat <- function(treat_set) {
tr <- ri_delta$prefecture_code %in% treat_set
mean(ri_delta$delta[tr]) - mean(ri_delta$delta[!tr])
}
obs_stat <- ri_stat(treated_codes)
set.seed(99)
n_perm <- 1000
placebo <- replicate(n_perm, ri_stat(sample(all_codes, length(treated_codes))))
ri_p <- (1 + sum(abs(placebo) >= abs(obs_stat))) / (1 + n_perm) # Phipson-Smyth 補正
cat(sprintf("観測 DID(post×treated, ASR): %.2f /10万\n", obs_stat))
cat(sprintf("RI 両側 p 値(プラセボ %d 回): %.3f\n", n_perm, ri_p))
```
```{r}
#| label: fig-ri
#| fig-cap: "プラセボ分布(47県から3県を処置に再割当)。赤線=実際の被災3県の推定。実推定が分布の外側にあれば「偶然では起きにくい」。"
#| code-fold: true
tibble(placebo = placebo) %>%
ggplot(aes(placebo)) +
geom_histogram(bins = 40, fill = "grey70", color = "white") +
geom_vline(xintercept = obs_stat, color = "red", linewidth = 1) +
geom_vline(xintercept = -obs_stat, color = "red", linetype = "dotted") +
labs(
title = "Randomization inference:プラセボ DID 分布",
subtitle = sprintf("両側 p = %.3f(デモ)", ri_p),
x = "プラセボ3県の post×treated DID(ASR/10万)", y = "頻度"
)
```
**読み方(推論)**:
- iid SE → 県クラスタ SE → RI、と進むほど不確実性は**正直に広がります**。個人 n の多さは見かけの精度を生むだけで、**真の情報量はクラスタ(県)数で決まる**。
- RI p 値は「47県から無作為に選んだ3県が、被災3県と同等以上の効果を**偶然**示す確率」。小さいほど効果は頑健。
- 根拠:Conley & Taber 2011(少数政策変更の推論)、MacKinnon & Webb 2018/2020(few/wild bootstrap・RI)、Phipson & Smyth 2010(p 値に +1 補正)。
---
# 期間集約(4期の解釈)
> **この章のねらい** — event-time 別の ATT_k を、即時(k=0)/調整(k=1–3)/再均衡(k≥4)の3期にまとめ、ショック後の動学を要約して読みやすくする。
```{r}
#| label: period-aggregation
#| code-fold: true
#| code-summary: "即時 / 調整 / 再均衡 の期間別平均効果"
period_of <- function(k) {
dplyr::case_when(
k == 0 ~ "即時 (k=0, 2011)",
k %in% 1:3 ~ "調整 (k=1-3, 2012-14)",
k >= 4 ~ "再均衡 (k>=4, 2015-)",
TRUE ~ "処置前 (k<0)"
)
}
es_all %>%
mutate(period = period_of(k)) %>%
filter(k >= 0) %>%
group_by(model, period) %>%
summarise(`平均効果%` = mean(pct), .groups = "drop") %>%
pivot_wider(names_from = period, values_from = `平均効果%`) %>%
knitr::kable(digits = 1, caption = "期間別の平均効果(%)。即時に大きく、その後減衰(デモ)。")
```
**読み方**:正の効果は「treatable death(医療で回避しえた死亡)の**増加**」=**医療システムの応答性低下**を示唆します。即時(k=0)に大きく、調整・再均衡で縮小していれば「**ショック → 回復**」の解釈と整合します。
---
# 先行研究の解析モデル
> **この章のねらい** — 本研究が踏まえる先行研究(Fukuma 2017 ほか)の解析手法と、本ノートで用いた方法論(DID/合成対照/少数クラスタ推論)の出典を整理し、本研究の位置づけを示す。
## 直接の先行研究:Fukuma ら 2017
- Fukuma S, Ahmed S, Goto R, Inui TS, Atun R, Fukuhara S. *Fukushima after the Great East Japan Earthquake: lessons for developing responsive and resilient health systems.* J Glob Health. 2017;7(1):010501. **PMID 28400956**. DOI 10.7189/jogh.07.010501. <https://pubmed.ncbi.nlm.nih.gov/28400956/>
- **アウトカム**:全年齢の **age-adjusted all-cause mortality** と **死因別**(急性心筋梗塞・脳血管疾患・癌・肺疾患・自殺)。標準人口は **日本 1985 年モデル人口**。
- **解析**:2005–2014年の**記述的トレンド比較**(日本平均・福島・宮城・岩手の前後比較)。**因果推論デザイン(DID/合成対照)は不使用**。
- **本研究の発展**:本研究は (1) アウトカムを **extended-treatable mortality** に絞り、(2) **被災3県 vs 44県の event-study DID / 合成対照**という反実仮想ベースの因果デザインに拡張し、(3) **少数クラスタ推論**を導入します。
## 関連する災害・経済ショックと死亡
- Maruthappu M, Watkins J, Noor AM, Williams C, Ali R, Sullivan R, Zeltner T, Atun R. *Economic downturns, universal health coverage, and cancer mortality in high-income and middle-income countries, 1990–2010: a longitudinal analysis.* Lancet. 2016;388(10045):684–695. **PMID 27236345**. <https://pubmed.ncbi.nlm.nih.gov/27236345/>
- Santos-Burgoa C, Sandberg J, Suárez E, et al. *Differential and persistent risk of excess mortality from Hurricane Maria in Puerto Rico.* Lancet Planet Health. 2018;2(11):e478–e488. **PMID 30318387**. <https://pubmed.ncbi.nlm.nih.gov/30318387/>
## 共変量選択と地域剥奪・社会経済要因(調整変数の根拠)
手法A(感度)・手法B で調整する **baseline 2010 SES(所得・学歴・失業・人口密度)** の選択根拠。
- Fukuda Y, Nakamura K, Takano T. *Higher mortality in areas of lower socioeconomic position measured by a single index of deprivation in Japan.* Public Health. 2007;121(3):163–173. **PMID 17222876**. DOI 10.1016/j.puhe.2006.10.015. <https://pubmed.ncbi.nlm.nih.gov/17222876/>
- 失業・低学歴・低所得を含む**7領域の地域剥奪指標**が県・市区町村の**<74歳死亡率**と相関(男 r=0.65)。所得・学歴・失業を選んだ直接の根拠。
- Nakaya T, Honjo K, Hanibuchi T, et al. *Associations of all-cause mortality with census-based neighbourhood deprivation and population density in Japan: a multilevel survival analysis.* PLoS One. 2014;9(6):e97802. **PMID 24905731**. DOI 10.1371/journal.pone.0097802. <https://pubmed.ncbi.nlm.nih.gov/24905731/>
- 地域剥奪と**人口密度**が個人 SES 調整後も全死亡を予測。人口密度を別個に入れる根拠。
- Fukuda Y, Nakao H, Yahata Y, Imai H. *Are health inequalities increasing in Japan? The trends of 1955 to 2000.* Biosci Trends. 2007;1(1):38–42. **PMID 20103865**. <https://pubmed.ncbi.nlm.nih.gov/20103865/>
- **県民所得 per capita** と age-adjusted 死亡率の関連(Poisson)。所得を入れる根拠。
- Ping R, Oshio T. *Educational inequalities in self-rated health and their mediators in late adulthood: comparison of China and Japan.* PLoS One. 2023;18(9):e0291661. **PMID 37713366**. DOI 10.1371/journal.pone.0291661. <https://pubmed.ncbi.nlm.nih.gov/37713366/>
- Takahashi S, Jang SN, Kino S, Kawachi I. *Gender inequalities in poor self-rated health: cross-national comparison of South Korea and Japan.* Soc Sci Med. 2020;252:112919. **PMID 32224365**. DOI 10.1016/j.socscimed.2020.112919. <https://pubmed.ncbi.nlm.nih.gov/32224365/>
- Furuya Y, Kondo N, Yamagata Z, Hashimoto H. *Health literacy, socioeconomic status and self-rated health in Japan.* Health Promot Int. 2015;30(3):505–513. **PMID 24131729**. DOI 10.1093/heapro/dat071. <https://pubmed.ncbi.nlm.nih.gov/24131729/>
- 学歴(および所得・就業)が日本の自己評価健康の頑健な決定因/交絡。学歴を入れる根拠。
- Kino S, Aida J, Kondo K, Kawachi I. *Do disasters exacerbate socioeconomic inequalities in health among older people?* Int J Disaster Risk Reduct. 2023;98:104071. **PMID 37982017**. DOI 10.1016/j.ijdrr.2023.104071. <https://pubmed.ncbi.nlm.nih.gov/37982017/>
- 東日本大震災の被災地 vs 非被災地を **DID** で比較し、教育・所得を SES 指標に用いた直接の災害文脈の先行例。
**共変量選択・過調整の方法論**(疫学方法論のため PMID + DOI):
- Schisterman EF, Cole SR, Platt RW. *Overadjustment bias and unnecessary adjustment in epidemiologic studies.* Epidemiology. 2009;20(4):488–495. **PMID 19525685**. DOI 10.1097/EDE.0b013e3181a819a1. <https://pubmed.ncbi.nlm.nih.gov/19525685/>
- 中間変数(やその下流)の調整=**過調整バイアス**。震災後 SES を入れない(baseline 限定)根拠。
- VanderWeele TJ. *Principles of confounder selection.* Eur J Epidemiol. 2019;34(3):211–219. **PMID 30840181**. DOI 10.1007/s10654-019-00494-6. <https://pubmed.ncbi.nlm.nih.gov/30840181/>
- Brookhart MA, Schneeweiss S, Rothman KJ, Glynn RJ, Avorn J, Stürmer T. *Variable selection for propensity score models.* Am J Epidemiol. 2006;163(12):1149–1156. **PMID 16624967**. DOI 10.1093/aje/kwj149. <https://pubmed.ncbi.nlm.nih.gov/16624967/>
- アウトカム関連変数は投入価値あり(分散減・バイアス増なし)。SES4変数の妥当性の根拠。
## 方法論(Event DID / 合成対照 / 少数クラスタ推論)
これらは PubMed 非収載の計量経済・統計手法論文のため **DOI + URL** で示します。
- Callaway B, Sant'Anna PHC. *Difference-in-differences with multiple time periods.* J Econometrics. 2021;225(2):200–230. DOI 10.1016/j.jeconom.2020.12.001. <https://doi.org/10.1016/j.jeconom.2020.12.001>
- Sun L, Abraham S. *Estimating dynamic treatment effects in event studies with heterogeneous treatment effects.* J Econometrics. 2021;225(2):175–199. DOI 10.1016/j.jeconom.2020.09.006. <https://doi.org/10.1016/j.jeconom.2020.09.006>
- Goodman-Bacon A. *Difference-in-differences with variation in treatment timing.* J Econometrics. 2021;225(2):254–277. DOI 10.1016/j.jeconom.2021.03.014. <https://doi.org/10.1016/j.jeconom.2021.03.014>
- Borusyak K, Jaravel X, Spiess J. *Revisiting event study designs: robust and efficient estimation.* Rev Econ Stud. 2024;91(6):3253–3285. DOI 10.1093/restud/rdae007. <https://doi.org/10.1093/restud/rdae007>
- Roth J. *Pretest with caution: event-study estimates after testing for parallel trends.* Am Econ Rev: Insights. 2022;4(3):305–322. DOI 10.1257/aeri.20210236. <https://doi.org/10.1257/aeri.20210236>
- Rambachan A, Roth J. *A more credible approach to parallel trends.* Rev Econ Stud. 2023;90(5):2555–2591. DOI 10.1093/restud/rdad018. <https://doi.org/10.1093/restud/rdad018>
- Sant'Anna PHC, Zhao J. *Doubly robust difference-in-differences estimators.* J Econometrics. 2020;219(1):101–122. DOI 10.1016/j.jeconom.2020.06.003. <https://doi.org/10.1016/j.jeconom.2020.06.003>
- Abadie A, Diamond A, Hainmueller J. *Synthetic control methods for comparative case studies.* J Am Stat Assoc. 2010;105(490):493–505. DOI 10.1198/jasa.2009.ap08746. <https://doi.org/10.1198/jasa.2009.ap08746>
- Abadie A. *Using synthetic controls: feasibility, data requirements, and methodological aspects.* J Econ Lit. 2021;59(2):391–425. DOI 10.1257/jel.20191450. <https://doi.org/10.1257/jel.20191450>
- Ben-Michael E, Feller A, Rothstein J. *The augmented synthetic control method.* J Am Stat Assoc. 2021;116(536):1789–1803. DOI 10.1080/01621459.2021.1929245. <https://doi.org/10.1080/01621459.2021.1929245>
- Conley TG, Taber CR. *Inference with "difference in differences" with a small number of policy changes.* Rev Econ Stat. 2011;93(1):113–125. DOI 10.1162/REST_a_00049. <https://doi.org/10.1162/REST_a_00049>
- MacKinnon JG, Webb MD. *The wild bootstrap for few (treated) clusters.* Econometrics J. 2018;21(2):114–135. DOI 10.1111/ectj.12107. <https://doi.org/10.1111/ectj.12107>
- MacKinnon JG, Webb MD. *Randomization inference for difference-in-differences with few treated clusters.* J Econometrics. 2020;218(2):435–450. DOI 10.1016/j.jeconom.2020.04.019. <https://doi.org/10.1016/j.jeconom.2020.04.019>
- Athey S, Imbens GW. *Design-based analysis in difference-in-differences settings with staggered adoption.* J Econometrics. 2022;226(1):62–79. DOI 10.1016/j.jeconom.2020.10.012. <https://doi.org/10.1016/j.jeconom.2020.10.012>
- Fay MP, Feuer EJ. *Confidence intervals for directly standardized rates: a method based on the gamma distribution.* Stat Med. 1997;16(7):791–801. **PMID 9131766**. <https://pubmed.ncbi.nlm.nih.gov/9131766/>
- Phipson B, Smyth GK. *Permutation P-values should never be zero.* Stat Appl Genet Mol Biol. 2010;9(1):Article 39. **PMID 21044043**. DOI 10.2202/1544-6115.1585. <https://pubmed.ncbi.nlm.nih.gov/21044043/>
- OECD/Eurostat. *Avoidable mortality: OECD/Eurostat lists of preventable and treatable causes of death.* 2019年版リスト(2022年1月公表). <https://www.oecd.org/content/dam/oecd/en/data/datasets/oecd-health-statistics/avoidable-mortality-2019-joint-oecd-eurostat-list-preventable-treatable-causes-of-death.pdf>
- Eurostat. *Preventable and treatable mortality statistics — Statistics Explained.*(treatable/preventable の定義と 50/50 按分ルールの一次解説). <https://ec.europa.eu/eurostat/statistics-explained/index.php?title=Preventable_and_treatable_mortality_statistics>
---
# まとめ・限界
> **この章のねらい** — 個人(手法A)と集計(手法B)の使い分けを整理し、本ノートの限界(シミュレーション前提など)を明示する。
## 個人 vs 集計の使い分け
- **個人データそのまま(手法A)**:treatable share による重み付けを主アウトカムとし(本番と同一)、年齢・性・死因の情報を保持するため、share の細分化や年齢×効果の異質性などの追加解析が柔軟。死亡が稀ならセル形式 Poisson が lossless かつ軽量。
- **集計データに加工(手法B)**:主解析(primary)の CS-DID。事前標準化した ASR で解釈が明快、共変量の回帰調整・少数クラスタ band(`cband`)を実装済み。
- **点推定は一致**するが、**推論はどちらでも県クラスタ数(=3)で決まる**。個人 n の多さは精度を生まない → RI / wild bootstrap が必須。
- **共変量調整**:年齢・性は両手法で回帰標準化(手法A=`age_band` FE、手法B=ASR)。**baseline 2010 SES(所得・学歴・失業・人口密度)**は、手法A では `i(event_time, z_X)` の交互作用、手法B では `xformla` で**同一セット**を投入。SES 調整は **primary ではなく事前指定の感度**で、PS-IPW は処置決定論(3県)で退化するため不使用。
## 限界(このノート)
- すべて**シミュレーション**。人口は2010年で固定、年齢構成は全国共通、treatable share はデモ用に少数死因へ簡略化(値は本番準拠の連続値、数え方=重み付けは本番と同一)。実データの主張ではない。
- 合成対照(`augsynth::multisynth`)・三重差分(DDD)・studentized randomization inference の完全実装は本ノートの範囲外とし、コア比較に集中。
- **baseline SES の生成**:本ノートの学歴・失業は所得と相関させた決定論的デモ値(実データでは国勢調査・労働力調査から取得)。処置効果と独立に生成しているため、SES 調整は ATT を実質的に動かさない設計。
- **3処置県での弱識別**:`SES × event-time` 交互作用は 44 非処置県の横断差から識別され、処置3県では ATT 項と部分的に共線になりやすい。SES 調整は感度に留め、節約版(`z_X : post`)も選択肢。実データ解析では事前登録 SAP の共変量セットに従う。
- 実データでの本解析は、事前登録した解析計画と、盲検化手順を含む再現可能なパイプラインに従う。
## 引用
本文中の各文献は **PMID(医学・疫学)/ DOI(手法)+ URL** を直接付しています(`.bib` 不使用)。一次情報は上記リンクから確認してください。
