「先生！心不全の患者さんの再入院データを解析したくて、2回目や3回目の再入院データもたくさんあったのですが、解析方法がよく分からなかったので、とりあえず『初回再入院までの時間（Time-to-first event）』だけでCox比例ハザードモデルを実行しました！無事にハザード比が出せました！」
……ちょっと待ってほしい。それは臨床研究の論文査読（Peer Review）において、統計査読者から「初回イベント日以降の膨大な臨床データを切り捨てており、治療の全疾病負荷（Cumulative Disease Burden）に対する効果を正しく評価できていない」と強く指摘される典型的なパターンである。
かといって、「1人の患者が3回再入院したから、3人分のデータ（3行）として通常のCoxモデルにそのまま放り込む」という解析を行ってしまうと、生存時間解析の根本前提である「観測値の独立性（Independence of observations）」が破破綻し、査読で一発却下されてしまう。
本記事では、なぜ「初回イベントのみ」や「全イベントの単純投入」が致命的な失敗を招くのかというロジックから、繰り返し発生するイベントを正しく解析する標準手法 Andersen-Gill（AG）モデル と cluster(id)（ロバスト標準誤差） の概念、Counting Process形式（tstart, tstop, status）へのデータ構造、Rによる実践コード、そして論文でそのまま使える英文テンプレートまでを徹底的に解説する。
1. なぜ「初回イベントのみ」や「全イベントの単純投入」は査読で即却下されるのか？
罠1：初回イベント解析（Time-to-first event）による情報の大量廃棄
心不全の再入院、がんの複数回再発、喘息やCOPDの増悪など、臨床現場では同じ患者に「同じイベントが何度も繰り返し発生する（Recurrent Events）」ケースが頻繁に存在する。
従来の生存時間解析（Kaplan-Meier法や標準的なCox比例ハザードモデル）で最もよく用いられるのが「初回イベント発生までの時間（Time-to-first event）」を評価するアプローチである。
しかし、このアプローチには致命的な臨床的弱点が存在する。初回イベントが発生した瞬間にその患者の追跡を打ち切り（Censoring）、2回目以降のイベントデータをすべてゴミ箱に捨ててしまう点である。
- 治療の本当の疾病抑止効果が見えなくなる: 新薬の介入によって「初回再入院」の時期は変わらなくても、「2回目・3回目の再入院率」が劇的に抑制されていた場合、初回イベント解析ではその薬の真の価値（全疾病負荷の軽減効果）を見落としてしまう。
- 検出力（Statistical Power）の著しい低下: せっかく長期間追跡して集めた貴重な「2回目以降のイベント情報」を一切活用しないため、サンプルサイズが無駄になり、統計的検出力が大幅に低下する。
罠2：全イベントを単一Cox回帰へ突っ込む「個人内相関（Within-subject correlation）」の無視
「初回イベントで捨ててしまうのが勿体ないなら、1人の患者が3回再入院したデータを3行の縦持ちデータにして、そのままCox回帰に入れればいいのでは？」と考えたくなるところだが、これはさらに危険な致命的誤用である。
標準的なCox比例ハザードモデルは、すべてのデータ行（観測値）が互いに独立していること（Independence of observations）を絶対的な前提条件としている。
しかし、「同じ患者 A さんの1回目の再入院」と「同じ患者 A さんの2回目の再入院」は独立ではない。 そこには患者ごとのクセ（遺伝的素因、生活習慣、基礎疾患の重症度、アドヒアランスなど）が存在し、同一患者内のイベント同士には強い 個人内相関（Within-subject correlation） が働く。
個人内相関を無視して標準的なCox回帰を実行すると、モデルは「独立した新しい患者がこんなにたくさんイベントを起こした！」と錯覚してしまう。その結果、標準誤差（Standard Error: SE）が不当に小さく見積もられ、$P$ 値が人工的に小さくなり（偽陽性 / 第一種の誤り）、間違った有意差を作り出してしまう。これが、査読者から即座に却下される理由である。
2. 全疾病負荷を正しく評価する「Andersen-Gill（AG）モデル」と Robust SE
Andersen-Gill（AG）モデルとは？
再発イベント解析において、世界的に最も標準的かつ汎用性の高いアプローチが Andersen-Gill（AG）モデル である。
AGモデルは、通常のCox比例ハザードモデルを「Counting Process（計数過程）」の枠組みへと拡張した手法である。その基本コンセプトは非常に明快で、「各イベントが発生した後も、患者は直ちに次のイベントのリスク集合（At-risk set）に戻り、全イベントを通して共通の基底ハザード（Baseline Hazard）を持つ」 と仮定する。
（※なお、イベントの回数（1回目・2回目…）によって基準リスクが変化すると仮定したい場合は PWP モデル、患者ごとのランダム効果を組み込みたい場合は Frailty モデルなどが用いられるが、全疾病負荷の全体像を評価する第一選択として国際的に支持されているのがAGモデルである。）
個人内相関を撃退する cluster(id)（ロバスト標準誤差）
AGモデルを実行する際、前述した「個人内相関（同一患者内の相関）」を正しく処理するために不可欠となるのが ロバスト標準誤差（Robust Standard Error / Sandwich Estimator） である。
Rの coxph 関数では、説明変数側に cluster(id) を追記するだけでこの補正が完了する。
- 点推定値（ハザード比：HR）: 全追跡期間における総イベント発生率の比（Overall Event Rate Ratio）として計算される。
- 標準誤差（SE）と $P$ 値:
cluster(id)を指定することで、同一患者（ID）のデータを「1つのクラスタ（固まり）」として扱い、個人内相関によって小さくなりすぎていた標準誤差を適切に膨らませて、正当な $P$ 値と95%信頼区間を算出する。
Andersen-Gillモデルのハザード式
被験者 $i$ の時刻 $t$ における AG モデルのハザード式は以下のように表される。$$\lambda_i(t) = \lambda_0(t) \exp(\boldsymbol{\beta}^T \mathbf{X}_i(t))$$
ここで重要なのは、基底ハザード $\lambda_0(t)$ がすべてのイベント回数（1回目、2回目、3回目…）を通じて共通であるという点である。イベントが発生しても時間はリセットされず、共通のタイムライン $t$ の上でリスク集合へと即座に復帰し、追跡が継続される。
3. Counting Process形式（tstart, tstop, status）へのデータ整形
1患者1行データから Counting Process 形式への変換
再発イベント解析を実行するには、データを標準的な「1患者1行（Time-to-first）」形式から、時間を区切った Counting Process 形式（1患者複数行データ） へと変換する必要がある。
Counting Process 形式では、各行に以下の3つの必須カラムを保持させる。
tstart: その観察区間の開始時間tstop: その観察区間の終了時間status: 区間終了時（tstop）におけるイベント発生の有無（1 = イベント発生, 0 = 打切り/観察継続）
データ構造の具体例
例えば、ID = 1 の患者が「10ヶ月目に1回目の再入院」「25ヶ月目に2回目の再入院」を起こし、「60ヶ月目で追跡終了（生存打切り）」となった場合、データは以下のように整形される。
| id | tstart | tstop | status | group |
| 1 | 0 | 10 | 1 | New_Drug |
| 1 | 10 | 25 | 1 | New_Drug |
| 1 | 25 | 60 | 0 | New_Drug |
- 1行目: 観察開始（0ヶ月）から10ヶ月目まで追跡され、1回目の再入院が発生（
status = 1）。 - 2行目: 10ヶ月目から即座に次のリスク集合に復帰し、25ヶ月目まで追跡されて2回目の再入院が発生（
status = 1）。 - 3行目: 25ヶ月目からリスク集合に復帰し、60ヶ月目まで追跡されてイベントなく観察終了（
status = 0）。
このようにデータを構築することで、全追跡期間における累積イベント発生状況を正確にモデルに組み込むことが可能となる。
※補足（関連記事）
- Counting Process形式へのデータ変換手順や
tmerge関数の詳細については 時間依存性共変量を含むCox回帰 を参照されたい。
4. 【実践Rコード】coxph ＋ cluster(id) による再発イベント解析（Step by Step）
ここからは、Rの survival パッケージを用いて、疑似的な再発イベントデータを作成し、従来の「初回イベント解析」と「Andersen-Gillモデル」を比較しながらStep by Stepで実行する手順を解説する。
Step 1: 疑似再発イベントデータの作成と確認
データのシミュレーション手順を4つのサブステップに分けて解説する。
Step 1-1: ライブラリの読み込みと基本患者データの設定
まずは必要なパッケージ（tidyverse, survival）を読み込み、患者レベルの背景情報（ID、割り付け群、年齢）を生成する。
library(tidyverse)
library(survival)
# 乱数シードの設定
set.seed(123)
n_patients <- 200
# 疑似患者データの生成（ID, グループ, 年齢）
patients <- tibble(
id = 1:n_patients,
group = factor(rep(c("Control", "Treatment"), each = n_patients / 2), levels = c("Control", "Treatment")),
age = round(rnorm(n_patients, mean = 65, sd = 10))
)
Step 1-2: 各患者の累積イベント時間のシミュレーション
各患者に対して、時間の経過に伴い繰り返し発生する再発イベントの時間間隔をシミュレーションする。
# 各患者の複数イベント発生間隔を生成する関数処理
df_sim_list <- map(patients$id, function(pid) {
p_info <- patients %>% filter(id == pid)
# イベント発生率の設定（Treatment群でイベント発生率が低い）
rate_param <- if_else(p_info$group == "Treatment", 0.03, 0.05)
# 指数分布に従うイベント発生間隔のシミュレーション（最大10回分）
intervals <- rexp(10, rate = rate_param)
cum_times <- cumsum(intervals)
# 最長追跡期間（36ヶ月＝3年）
max_followup <- 36
# 最長追跡期間までに発生したイベント時刻を抽出
valid_times <- cum_times[cum_times < max_followup]
all_times <- c(0, valid_times, max_followup)
tibble(
id = pid,
group = p_info$group,
age = p_info$age,
times = list(all_times),
n_events = length(valid_times)
)
})
Step 1-3: Counting Process形式（tstart, tstop, status）へのデータ変換
シミュレーションしたイベント発生時刻を、AGモデルに必要な Counting Process 形式（1患者複数行）へと整形する。追跡時間が0の不整合行が生じないよう、tstart < tstop のバリデーション処理を入れている。
# タイムラインを (tstart, tstop, status) 形式へ変換
df_recurrent <- map_dfr(df_sim_list, function(item) {
all_times <- item$times[[1]]
n_ev <- item$n_events
tibble(
id = item$id,
group = item$group,
age = item$age,
tstart = round(all_times[-length(all_times)], 1),
tstop = round(all_times[-1], 1),
status = c(rep(1, n_ev), 0) # 最後の区間は観察終了（status = 0）
)
}) %>%
filter(tstart < tstop) # 追跡時間が0の不整合行を除外
Step 1-4: 生成された縦持ちデータの構造確認
作成された Counting Process 形式のデータの先頭行を確認し、1人の患者が複数の観察区間を持っていることを確認する。
# データの先頭（ID=1の複数区間データなど）を確認
head(df_recurrent, 10)
Step 2: 従来の「初回イベントのみ」のCox回帰（比較用）
比較のために、初回イベント（各患者の最初の区間 tstart == 0）のみを切り出した従来通りのCox比例ハザードモデルを実行する。
# 初回イベントデータのみの抽出（各患者の1行目データ）
df_first_event <- df_recurrent %>%
group_by(id) %>%
slice(1) %>%
ungroup()
# 初回イベントのみのCox比例ハザードモデル
fit_first <- coxph(Surv(tstop, status) ~ group + age, data = df_first_event)
summary(fit_first)
Step 3: Andersen-Gill（AG）モデルの実行（coxph ＋ cluster(id)）
全再発イベントデータを用い、Counting Process形式の Surv(tstart, tstop, status) と個人内相関を補正する cluster(id) を組み込んでAGモデルを実行する。
# Andersen-Gill モデルの実行（Robust SE 適用）
fit_ag <- coxph(
Surv(tstart, tstop, status) ~ group + age + cluster(id),
data = df_recurrent
)
# モデル結果のサマリー表示
summary(fit_ag)
Step 4: 初回イベント解析 vs AGモデルの結果比較と解釈
出力結果から、治療群（Treatment vs Control）のハザード比（HR）と $P$ 値を比較する。
- 初回イベントモデル（
fit_first）: 初回イベント発生までの遅延効果のみを評価する。2回目以降のイベントが捨てられているため、信頼区間が広く（SEが大きく）なる傾向がある。 - Andersen-Gillモデル（
fit_ag）: 全追跡期間中の「総イベント発生率（Overall Event Rate）」を評価する。すべてのイベント情報を活用するため、イベント数が多くなり検出力が向上すると同時に、cluster(id)によって正確な Robust SE が算出される。
5. 論文でそのまま使える英文テンプレート（Methods & Results）
論文（MethodsおよびResultsセクション）にそのまま記述できる標準的な英文表現テンプレートを示す。
Methods セクション英文テンプレート
“To evaluate the overall disease burden and the effect of intervention on recurrent [events, e.g., heart failure hospitalizations], the Andersen-Gill (AG) extension of the Cox proportional hazards model was performed. Data were formatted into the counting process structure with time-interval variables (
tstart,tstop) and an event indicator (status).To account for within-subject correlation among recurrent events within the same patient, robust standard errors were calculated using a sandwich estimator clustered on patient identity (
cluster(id)). The treatment effect was presented as an overall event rate ratio (hazard ratio [HR]) with its corresponding 95% confidence interval (CI). For comparison, a conventional time-to-first event Cox regression analysis was also conducted. All analyses were performed using R (version 4.x.x) with thesurvivalpackage.”
Results セクション英文テンプレート
“Over a median follow-up of [36] months, a total of [245] recurrent [heart failure hospitalizations] occurred among [200] patients ([140] events in the Control group and [105] events in the Treatment group).
In the primary Andersen-Gill model evaluating all recurrent events with robust standard errors, [Treatment] was associated with a significant reduction in the overall risk of [hospitalization] (Hazard Ratio [HR], [0.75]; 95% CI, [0.60] to [0.94];$P = 0.013$). In contrast, the conventional time-to-first event analysis demonstrated a HR of [0.78] (95% CI, [0.56] to [1.08];$P = 0.14$), underestimating the cumulative clinical benefit of the treatment by ignoring subsequent recurrent events.”
6. まとめ＆一括コピペ用Rスクリプト
本記事のポイントおさらい
- 初回イベント（Time-to-first event）だけの解析は勿体ない: 2回目以降の再発イベントを切り捨ててしまうため、疾患の全負担（Cumulative Disease Burden）を見落とし、検出力を低下させる。
- 全イベントの単純投入は絶対NG: 個人内相関（Within-subject correlation）を無視すると標準誤差が不当に小さくなり、偽陽性（$P < 0.05$）を起こして査読で即却下される。
- 解決策は Andersen-Gill（AG）モデル ＋
cluster(id): Counting Process形式（tstart,tstop,status）にデータを整形し、cluster(id)でロバスト標準誤差を計算することで、個人内相関を正当に補正した全イベント発生率の評価が可能となる。
一括実行用Rスクリプト
以下のコードをコピー＆ペーストすることで、疑似データの作成、Counting Process形式の確認、初回イベント解析、AGモデルの実行、および論文表示用テーブル（gtsummary）の出力までを一気通貫で実行できる。
# ==============================================================================
# 再発イベント生存時間解析（Andersen-Gill Model & Robust SE）一括実行スクリプト
# パッケージ: survival, tidyverse, gtsummary
# ==============================================================================
if (!requireNamespace("survival", quietly = TRUE)) install.packages("survival")
if (!requireNamespace("tidyverse", quietly = TRUE)) install.packages("tidyverse")
if (!requireNamespace("gtsummary", quietly = TRUE)) install.packages("gtsummary")
library(survival)
library(tidyverse)
library(gtsummary)
# 1. 疑似データの作成（Counting Process形式）
set.seed(123)
n_patients <- 200
patients <- tibble(
id = 1:n_patients,
group = factor(rep(c("Control", "Treatment"), each = n_patients / 2), levels = c("Control", "Treatment")),
age = round(rnorm(n_patients, mean = 65, sd = 10))
)
df_recurrent <- map_dframe(patients$id, function(pid) {
p_info <- patients %>% filter(id == pid)
rate_param <- if_else(p_info$group == "Treatment", 0.03, 0.05)
intervals <- rexp(10, rate = rate_param)
cum_times <- cumsum(intervals)
max_followup <- 36
valid_times <- cum_times[cum_times < max_followup]
all_times <- c(0, valid_times, max_followup)
tibble(
id = pid,
group = p_info$group,
age = p_info$age,
tstart = round(all_times[-length(all_times)], 1),
tstop = round(all_times[-1], 1),
status = c(rep(1, length(valid_times)), 0)
)
}) %>% filter(tstart < tstop)
# 2. 初回イベント解析（Time-to-first event）
df_first <- df_recurrent %>% group_by(id) %>% slice(1) %>% ungroup()
fit_first <- coxph(Surv(tstop, status) ~ group + age, data = df_first)
# 3. Andersen-Gill モデル解析（Recurrent Events + Robust SE）
fit_ag <- coxph(Surv(tstart, tstop, status) ~ group + age + cluster(id), data = df_recurrent)
# 4. 結果の出力
cat("\n==============================================================================\n")
cat(" 1. 初回イベントのみの解析 (Time-to-First) \n")
cat("==============================================================================\n")
print(summary(fit_first))
cat("\n==============================================================================\n")
cat(" 2. Andersen-Gill モデル解析 (Recurrent Events + Robust SE) \n")
cat("==============================================================================\n")
print(summary(fit_ag))
# 5. gtsummary による論文用テーブルの作成
tbl_ag <- fit_ag %>%
tbl_regression(exponentiate = TRUE) %>%
bold_p() %>%
modify_caption("**Table. Multivariable Andersen-Gill Model for Recurrent Events**")
print(tbl_ag)
コメント