【R言語】層別化Cox回帰の完全ガイド！strata()でハザード比が消える理由から論文の記述作法まで

研究者: 査読者から「施設間差を層別化Coxモデルで調整しなさい」と言われ、strata() を使って解析したのですが……指定した変数のハザード比（HR）が出力から消えてしまいました！エラーでしょうか？

統計ER: ご安心ください、それはエラーではなく層別化Cox回帰の正しく正常な動作です！ 指定した変数のリスクは層ごとの基準ハザード関数に吸収されるため、最初からハザード比が計算されない構造になっているんですよ。さっそく仕組みと活用法を解説しましょう！

多変量解析としてCox比例ハザードモデルを実行している際、論文の査読者（Reviewer）や指導医から以下のような指摘を受けたことはないだろうか。

• 「特定変数の比例ハザード性が破綻しているので、層別化Coxモデル（Stratified Cox Model）を適用しなさい」
• 「多施設共同研究における施設間差（Center Effect）は、共変量としてではなく層別化でコントロールしなさい」

指示通りにRの coxph() 関数内で strata() を使って解析を再実行したところ、今度は新たな問題が発生してパニックになる研究者が後を絶たない。

「strata() に指定した変数のハザード比（HR）やP値が出力結果から消えてしまった……。コードを間違えたのだろうか？」

本記事では、研究者が一度は焦る「ハザード比が消える理由」を分かりやすく解説するとともに、層別化Cox回帰が必要な背景、通常の共変量調整との根本的な違い、Rでの書き方、交互作用の評価法、論文掲載時の英文記述作法までを網羅して解説する。

なぜ「層別化Cox回帰」を使うのか？通常の多変量調整との違い

まず、通常の多変量Cox回帰（+ covariate）と、層別化Cox回帰（+ strata(covariate)）の統計学的な決定的な違いを整理しよう。

1. 通常の多変量調整（+ covariate）の仕組みと限界

通常の多変量Cox回帰では、観察対象の全集団に対して「ただ1つの共通な基準ハザード関数 $h_0(t)$」を仮定する。

そして、調整変数（例: 年齢や病期 Stage）は、この共通の生存曲線を「上下に平行移動させる影響（一定のハザード比 $\exp(\beta)$）」を持つと見なす。

このモデルが成立するための絶対条件が、「比例ハザード性の仮定（Proportional Hazards Assumption）」である。時間経過に伴って各群のハザード比が変化しない（生存曲線同士が交差しない）ことが前提となる。

2. 層別化（+ strata()）が解決する2つの課題（PH仮定破綻・施設差）

一方、層別化Cox回帰では、指定した層（Stratum / グループ）ごとに「全く異なる形状の基準ハザード関数 $h_{0g}(t)$」を持つことを許容する。つまり、層ごとに生存曲線の形が異なっていても、途中で曲線同士が交差（Cross）していても構わないのである。

ポイント：層別化Cox回帰が威力を発揮する2大臨床シーン

1. 比例ハザード性が破綻している変数の調整 > Schoenfeld残差プロット等で時間経過に伴うリスク変動（$p < 0.05$）が確認された変数（例: 病期 Stage や 年齢階層）を層化し、バイアスを回避する。
2. 多施設共同研究における施設間差（Center Effect）のコントロール > 施設ごとのベースラインリスク（初期死亡率や診断技術の違い）を層ごとに自由に吸収させつつ、主要な治療効果（HR）を一括で正確に推計する。

【実践】Rでの層別化Cox回帰の実装（coxph + strata）とサマリーの読み解き

Rで層別化Cox回帰を実行するのは非常に簡単である。survival パッケージの coxph() 関数内で、層別化したい変数を strata() で囲むだけである。

初心者が陥りやすい注意点：連続変数の層別化 > strata() に渡す変数は、必ずカテゴリ変数（因子型：factor）でなければならない。連続変数（例: age）をそのまま渡すと、値ごとに膨大な層が作られて解析が破綻するため、あらかじめ cut() 等で年代別にカテゴリ化しておこう。

あらかじめ必要なパッケージを読み込み、動作確認用のサンプルデータを作成しよう。

library(tidyverse)
library(survival)

# サンプルデータの作成
set.seed(123)
n <- 200
df <- tibble(
  time   = rexp(n, rate = 0.05),
  status = sample(c(0, 1), n, replace = TRUE, prob = c(0.3, 0.7)),
  treat  = factor(sample(c("Control", "Drug"), n, replace = TRUE)),
  age    = rnorm(n, mean = 60, sd = 10),
  center = factor(sample(c("Center_A", "Center_B", "Center_C"), n, replace = TRUE)) # 施設間差
)

基本コードの書き方（survival::coxph）

# 1. 通常の多変量Cox回帰（施設間差 center をそのまま投入）
fit_normal <- coxph(Surv(time, status) ~ treat + age + center, data = df)
summary(fit_normal)

# 2. 層別化Cox回帰（施設間差 center を層別化）
fit_strat <- coxph(Surv(time, status) ~ treat + age + strata(center), data = df)
summary(fit_strat)

なぜ層化変数のハザード比（HR）が出力から消えるのか？

実行結果を summary() で確認してみよう。

summary(fit_strat)

出力結果を見てみると、treat（治療群）や age（年齢）のハザード比（exp(coef)）やP値は表示されるが、strata(center) として指定した center のハザード比の行が完全に消えている。

消える理由の数学的ロジック > 層別化Cox回帰のモデル式は、層$g$（各施設）において以下のように表される。$$h_g(t) = h_{0g}(t) \exp(\beta_1 \cdot \text{treat} + \beta_2 \cdot \text{age})$$

施設ごとのリスクの違いは、すべて層ごとの基準ハザード関数$h_{0g}(t)$の内部に完全に吸収される。

つまり、層別化変数（center）に対しては最初から回帰係数$\beta$を定義・推定しない構造になっている。そのため、「ハザード比が表示されない」のが統計学的に正しく、正常な出力結果なのである。

【発展】層ごとに治療効果が異なるか？（交互作用・効果修飾の評価）

デフォルトの層別化Cox回帰（strata()）は、「基準ハザードは層ごとに異なるが、治療（treat）のハザード比 $\beta_1$ は全層で共通（等しい）」という仮定を置いている。

しかし、臨床研究では「施設ごと、あるいは病期（Stage）ごとに治療の効果（ハザード比）そのものが異なるのではないか？（サブグループ解析 / 交互作用）」を検証したい場面がある。

層ごとに個別のハザード比を推定したい場合は、モデル内で交差項（treat * strata(...)）を指定する。

# 層ごとに個別治療効果（ハザード比）を求める交差項モデル
fit_inter <- coxph(Surv(time, status) ~ age + treat * strata(center), data = df)
summary(fit_inter)

このコードを実行すると、Center_A を基準とした各施設における治療のハザード比や、施設ごとの効果の違い（交互作用のP値）を評価することができる。

査読者（Reviewer）を納得させる論文・学会発表での記述作法

層別化Cox回帰を用いて得られた結果を論文に掲載する際は、査読者に「正しく層別化解析が行われたこと」を伝える標準的な英文表現を用いる必要がある。

解析方法（Methods）の英文テンプレート

• 比例ハザード仮定の破綻に対応した場合:“To account for the non-proportional hazards of disease stage, a stratified Cox proportional hazards model was fitted by stratifying on disease stage.”
• 多施設共同研究の施設間差をコントロールした場合:“To control for potential clustering and center-specific baseline risks, a stratified Cox proportional hazards model was utilized with study center as a stratifying factor.”

結果表（Table 2など）の脚注（Footnote）記載例

層別化を行った変数は表の行から消えるため、必ず脚注に「どの変数で層別化したか」を明記しておくのがマナーである。

“Hazard ratios (HRs) and 95% confidence intervals (CIs) were estimated using a Cox proportional hazards model stratified by study center and adjusted for age and sex.”

【まとめ】コピペで動く！層別化Cox回帰の一括Rコードテンプレート

本記事で解説した「通常モデル ➔ 層別化モデル ➔ サマリー確認 ➔ 交互作用モデル」の一連の流れをまとめた一括コードを以下に示す。自身の手元のデータセット名・変数名に置き換えて活用してほしい。

# === 層別化Cox回帰 解析テンプレート ===
library(tidyverse)
library(survival)

# 1. 通常の多変量Cox回帰
fit_normal <- coxph(Surv(time, status) ~ treat + age + sex + stage, data = df)
summary(fit_normal)

# 2. 層別化Cox回帰（PH仮定が破綻した stage や施設差を strata に投入）
# ※ stage は事前に因子型（factor）になっていることを確認！
fit_strat <- coxph(Surv(time, status) ~ treat + age + sex + strata(stage), data = df)

# サマリー確認（※ stage のHRが出力されないのは正常動作！）
summary(fit_strat)

# 3. （発展）層ごとの治療効果の違い（交互作用）を評価したい場合
fit_inter <- coxph(Surv(time, status) ~ age + sex + treat * strata(stage), data = df)
summary(fit_inter)

「strata() を指定するとその変数のハザード比が出なくなる」という仕様を正しく理解しておけば、もう出力結果を見て焦る必要はない。

比例ハザード性の破綻や多施設データの交絡をスマートにクリアし、査読に通る頑健で信頼性の高い研究論文を完成させよう。

おすすめ書籍

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