【R言語】カルテの「繰り返し測定データ」や「途中の治療変更」をそのままCox回帰に入れていませんか？時間依存共変量とtmerge関数の完全ガイド

研究者: 観察期間中に毎月測定した血圧の変化や、追跡途中で追加した抗がん剤の影響をそのままCox回帰モデルの共変量（1/0）として投入しても大丈夫でしょうか？

統計ER: ちょっと待ってください！何の気なしに扱おうとしているそのデータは、統計学で『時間依存共変量（Time-dependent covariate）』と呼ばれる超高度なデータ構造です。思い込みや直感で普通のCox回帰に入れてしまうと「不死時間バイアス」が生じ、論文が一発で即却下（Reject）されます！正しく勉強して適切なデータ変換（tmerge）を行い、正しく解析しましょう！

臨床研究やリアルワールドデータ（RWD）の解析において、電子カルテやレジストリからデータを抽出し、以下のような因子を生存時間解析（Cox比例ハザードモデル）に組み込もうと考えたことはないだろうか。

• 「追跡期間中に毎月測定した血圧や血清マーカー（CRPやHbA1cなど）の動的な変化を考慮したい」
• 「観察開始から6ヶ月目に、副作用の出現や病勢進行によって抗がん剤を追加・変更した影響を正しく評価したい」

「臨床的に価値のあるデータを余すことなくモデルに入れたい」という着眼点自体は非常に素晴らしい。

しかし、直感や思い込みだけで、これらの「追跡途中に測定・発生したデータ」をそのまま通常のCox回帰の共変量（例: 新薬使用＝1, 未使用＝0）として投入してしまうのは極めて危険である。

本来であれば、こうした複雑な縦持ちデータや追跡途中のイベント処理は、専門の生物統計家（バイオスタティスティシャン）に相談してモデルを構築すべき領域である。

だが、日本の臨床研究の現状では、多忙な医師や医療従事者が自力で解析に取り組まざるを得ないケースが依然として多い。

だからこそ、「ノリや思い込みで解析を進める」のは絶対にやめなければならない。

自分が扱おうとしているデータが統計学的にどのような性質を持っているのかを正しく勉強し、適切なデータ構造へ変換した上で解析に挑む必要がある。

本記事では、初心者が最も陥りやすい「不死時間バイアス（Immortal Time Bias）」の恐怖と、それを解消するためのRの神関数 survival::tmerge() を使った時間依存Cox回帰の実装手順を分かりやすく解説する。

なぜ「通常のCox回帰」にそのまま入れてはいけないのか？（不死時間バイアスの恐怖）

追跡途中に発生するデータを普通の共変量としてモデルに入れてはいけない最大の理由は、統計学における「時間依存共変量」のルールを無視することによって、致命的なバイアスが発生するためである。

1. 時間固定共変量（Time-fixed）と時間依存共変量（Time-dependent）の違い

まずは変数の定義を整理しよう。

• 時間固定共変量（Time-fixed covariate）:観察開始時点（ベースライン: $t = 0$）で一決決定され、追跡期間中に変化しない変数（例: 性別、登録時の年齢、遺伝子変異の有無など）。
• 時間依存共変量（Time-dependent covariate）:観察開始後の時間経過に伴って、値やステータスが動的に変化する変数（例: 追跡中の血圧値、追跡途中での併用療法の開始、有害事象の発症など）。

通常のCox回帰（coxph(Surv(time, status) ~ x)）は、すべての共変量が「ベースライン時点（$t = 0$）で固定されていること（時間固定）」を大前提として構築されている。

2. 論文が即却下される「不死時間バイアス（Immortal Time Bias）」とは？

例えば、「がん患者の生存率に対する新薬Aの効果」を検証する研究を考えてみよう。

患者によっては、診断直後（$t = 0$）から新薬Aを飲む人もいれば、病状が悪化した追跡6ヶ月目に初めて新薬Aを開始する人もいる。

ここで、追跡6ヶ月目に新薬を開始した患者を、何の気なしに「最初から新薬群（新薬フラグ = 1）」として通常のCox回帰に投入してしまったらどうなるだろうか？

注意！不死時間バイアス（Immortal Time Bias）のメカニズム

その患者は「6ヶ月目に新薬を開始できた」ということは、少なくとも最初の6ヶ月間は『絶対に死亡しなかった（生存していた）』という未来の情報を持っている。

本来は未治療であった最初の6ヶ月間までも「新薬群の生存時間」としてカウントしてしまうため、新薬群に「絶対に死亡が起きない無敵の期間（不死時間: Immortal Time）」が不当に加算され、新薬の効果が恐ろしいほど過大評価されてしまうのである。

査読者（Reviewer）はこの「不死時間バイアス」に極めて敏感である。何の気なしに手作業でフラグを立てた解析を提出すると、統計学的な妥当性なしと判断され、論文は一発で却下される。

3. 「時間依存共変量」と「時間依存係数（比例ハザード性破綻）」の区別

ここで臨床医が最も混乱しやすい2つの用語の違いを整理しておこう。

• 時間依存共変量 (Time-dependent covariate):変数 $X(t)$ 自体が時間とともに変化する（例: 検査値の変化、途中の治療開始）。➔ 本記事で解説する tmerge() で対処
• 時間依存係数 (Time-varying coefficient):効果の大きさ（ハザード比 $\beta(t)$）が時間とともに変化する（例: 治療初期はハザード比が良いが、後期になると効果が薄れて生存曲線が交差する＝比例ハザード性の破綻）。

ハザード比自体が時間によって変化してしまう「時間依存係数」や比例ハザード性の破綻（交差する生存曲線）への対処法（RMSTや層別化、期間分割モデルなど）については、以下の記事で詳しく解説しているため、ぜひ併せて参照してほしい。

➔ Cox回帰で比例ハザード性が怪しい時の対処法！RMST（制限付き平均生存時間）の使い方と層別化・期間分割の使い分け

Rでの解決策：Counting Process 形式 (start, stop, event) へのデータ変換

不死時間バイアスを回避し、時間依存共変量をCox回帰で正しく扱うための数学的アプローチが「Counting Process 形式（カウント過程形式）」である。

これは、1人の患者（1行データ）を「ステータスが一定であった時間区間（Interval）」ごとに複数行に切り分ける手法である。

例えば、「6ヶ月目に新薬を開始し、18ヶ月目に死亡した患者」の場合、データ行を以下のように2行に分割する。

患者ID	開始時間 (tstart)	終了時間 (tstop)	新薬ステータス	イベント (status)	解釈
ID: 1	0 ヶ月	6 ヶ月	0 (未服用)	0 (生存)	新薬を飲んでいない期間としてカウント
ID: 1	6 ヶ月	18 ヶ月	1 (服用中)	1 (死亡)	新薬を飲んでいる期間としてカウント

このようにデータを区間分割すれば、0〜6ヶ月目は「コントロール群のリスク」、6〜18ヶ月目は「新薬群のリスク」として正しくハザード計算に組み込まれ、不死時間が完全に解消される。

【実践】神関数 survival::tmerge() でデータを一瞬で整形する

データ行を自分で手動計算して複数行に切り分けるのは、計算ミスやバグの温床となる。

Rの標準パッケージ survival に含まれる tmerge() 関数の仕組みを使えば、ベースラインデータに時間依存のイベントや繰り返し測定データを安全かつ自動で統合・区間分割できる。

手元のR環境で1ステップずつコードを実行しながら、データの構造変化を確認していこう。

Step 1: 必要なパッケージの読み込みとベースラインデータの作成

まずは、患者1人につき1行の基本情報（ID、年齢、追跡期間、転帰）を持つベースラインデータを作成する。

library(tidyverse)
library(survival)

# 1. ベースラインデータの作成（1人1行）
set.seed(123)
n <- 200
base_df <- tibble(
  id      = 1:n,
  age     = rnorm(n, mean = 60, sd = 10),
  futime  = runif(n, min = 1, max = 36), # 追跡期間（ヶ月）
  death   = sample(c(0, 1), n, replace = TRUE, prob = c(0.7, 0.3)) # 転帰
)

# 作成したベースラインデータの確認（1人1行になっている）
head(base_df)

Step 2: 追跡途中の「時間依存イベントデータ」の準備

次に、追跡途中で発生したイベント（例: 新薬の服用開始）とその発生タイミングを持つ別データフレームを用意する。

# 2. 追跡途中の「治療開始イベント」データの作成
# 一部の患者のみ、追跡途中で新薬を開始したとする
set.seed(456)
drug_df <- tibble(
  id        = sample(1:n, 80),
  drug_time = runif(80, min = 0.5, max = 20) # 新薬を開始したタイミング（ヶ月）
) %>% 
  # 追跡期間より後に新薬開始となっている不整合データを防ぐフィルター
  inner_join(base_df, by = "id") %>% 
  filter(drug_time < futime) %>% 
  select(id, drug_time)

# イベントデータの確認（誰が何ヶ月目に新薬を開始したか）
head(drug_df)

Step 3: tmerge() の1回目：ベースラインデータと最終転帰の登録

tmerge() 関数の第1段階として、ベースラインデータに全追跡期間（futime）と最終イベント（death）を割り当てる。

# 1回目：ベースラインデータと追跡期間・最終イベントを設定
df_split <- tmerge(
  data1 = base_df,
  data2 = base_df,
  id    = id,
  event = event(futime, death) # 追跡時間と最終イベント
)

# 1回目実行後のデータ構造確認（tstart, tstop が自動生成される）
head(df_split)

Step 4: tmerge() の2回目：tdc() による時間依存共変量の結合と区間分割

ここが最も重要なステップである。2回目の tmerge() で drug_df を渡し、tdc(drug_time) 関数を使用する。

tdc() は Time-Dependent Covariate の略であり、「指定した時間（drug_time）に到達するとフラグが 0 から 1 へ変化する」という指示を意味する。

# 2回目：時間依存イベント（新薬開始）を統合して区間分割
df_split <- tmerge(
  data1 = df_split,
  data2 = drug_df,
  id    = id,
  drug  = tdc(drug_time) # tdc() により新薬開始タイミングで区間分割！
)

# 複数行に分割された Counting Process 形式データの確認
head(df_split %>% select(id, tstart, tstop, event, drug, age), 10)

tdc(drug_time) を指定することで、drug_time より前の区間では drug = 0、drug_time 以降の区間では drug = 1 に自動で切り替わるフラグがセットされ、データが複数行に完璧に分割される。

時間依存Cox回帰の実行（coxph）とハザード比の正しい読み解き

データが Counting Process 形式（tstart, tstop, event）に整形されたら、Cox比例ハザード回帰を実行する。

1. coxph(Surv(tstart, tstop, event) ~ ..., cluster(id)) の書き方

Surv() 関数の第一引数・第二引数に tstart と tstop を指定する。

また、1人の患者が複数行にまたがっているため、標準誤差の過小評価を防ぐ目的で cluster(id) オプション（堅牢な標準誤差: Robust Standard Error）を必ず指定する。

# 時間依存Cox比例ハザードモデルの実行
fit_tdc <- coxph(
  Surv(tstart, tstop, event) ~ drug + age + cluster(id),
  data = df_split
)

# 解析結果の出力
summary(fit_tdc)

2. 算出されたハザード比（HR）の臨床的解釈

得られた drug（新薬）のハザード比（exp(coef)）は、以下のように正しく臨床解釈する必要がある。

ハザード比の臨床的解釈

「全観察期間を通じた固定的なリスク比」ではなく、『各時点において、その瞬間に新薬を服用している患者群は、服用していない患者群と比較して、イベント発生リスクがどれほど変化するか』を表す瞬時リスク比である。

査読者（Reviewer）を納得させる論文記述法（Methods ＆ Results）

時間依存共変量を用いた解析を論文へ報告する際の標準的な英文テンプレートを以下に示す。

解析方法（Methods）の英文例:

“To evaluate the association between in-study treatment changes [または repeating biomarker measurements] and overall survival without immortal time bias, a extended Cox proportional hazards model with time-dependent covariates was utilized. The dataset was converted into counting process format (start, stop, event) using the tmerge function in R. Robust standard errors were calculated using cluster-robust variance estimators to account for multiple observations per patient.”

結果欄（Results / Table 脚注）の英文例:

“After adjusting for baseline age and time-dependent treatment status, the initiation of [新薬名] was significantly associated with a reduced risk of mortality (time-dependent hazard ratio, 0.62; 95% CI, 0.45 to 0.86; p = 0.004).”

【まとめ】コピペで動く！時間依存Cox回帰の一括Rコードテンプレート

本記事で解説した「データ作成 ➔ tmerge() によるデータ構造化 ➔ 時間依存Cox回帰の実行」を一貫して試せる一括テンプレートコードを以下にまとめた。自身の手元のデータフレーム名・変数名に置き換えて活用してほしい。

# === 時間依存共変量（Time-Dependent Covariate）Cox回帰 一括テンプレート ===
library(tidyverse)
library(survival)

# Step 1: データの準備（base_data: ベースライン, time_dep_data: 時間依存イベント）
# ※ ご自身の実データに置き換えてください

# Step 2: tmerge() によるデータ分割とカウント過程形式化
df_tdc <- tmerge(
  data1 = base_data,
  data2 = base_data,
  id    = id,
  event = event(futime, status) # 追跡期間とイベント
)

df_tdc <- tmerge(
  data1 = df_tdc,
  data2 = time_dep_data,
  id    = id,
  treat_tdc = tdc(event_time) # 時間依存変数の開始タイミング
)

# Step 3: 時間依存Cox比例ハザードモデルの実行
fit_tdc <- coxph(
  Surv(tstart, tstop, event) ~ treat_tdc + age + sex + cluster(id),
  data = df_tdc
)

# Step 4: 解析結果表示（時間依存ハザード比の確認）
summary(fit_tdc)

臨床現場のカルテデータに溢れている「繰り返し測定値」や「途中の治療変更」を、何の気なしにそのまま回帰モデルに投入してはいけない。

統計的な仕組みを正しく勉強し、tmerge() を用いてデータ構造を適切に整えることで、不死時間バイアスを排した真に信頼性の高い研究成果を論文として発信していこう。

おすすめ書籍

誰も教えてくれなかった　医療統計の使い分け〜迷いやすい解析手法の選び方が，Rで実感しながらわかる！