研究者: 治療法が3群（A薬・B薬・C薬）ある生存時間解析を行っています。全体のログランク検定で有意差が出たのですが、どの群同士に差があるのかの調べ方や、Cox回帰の基準群の設定方法、グラフの綺麗な描き方が分かりません……。
統計ER: 3群以上の比較では「多重性の補正」と「基準群（Reference）の制御」が重要です！
survminer::pairwise_survdiff()を使えばHolm法等で補正した事後比較が一瞬で出力でき、ggsurvplot()で論文用のKM曲線も綺麗に描けます。比例ハザード性は検定に頼らず、Schoenfeld残差プロットで視覚的に確認しましょう！
臨床研究や観察研究において、治療アプローチが3つ以上に分かれるケース（例：A薬 vs B薬 vs C薬、または外科手術 vs 化学療法 vs 放射線療法）は非常に多い。しかし、生存時間解析において比較対象が3群以上になった途端、以下のような統計学・R実装上の壁にぶち当たる研究者が急増する。
- 「全体のログランク検定で有意差（$p < 0.05$）が出たが、どの群とどの群の間に有意差があるのか分からない」
- 「2群ずつのログランク検定を総当たりで繰り返しても統計学的に許されるのか？」
- 「Cox比例ハザードモデルで、自分が対照群（Reference）にしたい群を基準にしてハザード比（HR）を算出する方法が分からない」
- 「3群の生存曲線（カプランマイヤー曲線）を描くと、線や Risk Table がごちゃごちゃして論文用に綺麗に整えられない」
生存時間解析の3群比較には、多重性の調整や因子水準の制御といった特有の「作法」が存在する。これを誤ると、偽陽性を増やしたり、査読者（Reviewer）から厳しく指摘されたりする原因となる。
本記事では、全体のログランク検定から事後ペアワイズ多重比較、Cox回帰の基準群の考え方、論文掲載レベルのKM曲線描画、グラフィカルな比例ハザード性の確認までをRで一気通貫に完結させる手順を解説する。
ログランク検定で差が出たら！pairwise_survdiff() による事後多重比較
3群以上の生存曲線に差があるかを調べる際、まずは全体に対して一括でログランク検定（Global Log-rank Test）を行う。
Rでは survival::survdiff() を用いるが、ここで $p < 0.05$ となった場合、分かるのは「3群のどこか少なくとも1箇所に分布の差がある」ということだけである。具体的にどの群間に差があるのか（事後比較 / Post-hoc Comparison）を突き止めたい場合、多重比較をする必要がある。
このとき、「調整なしのログランク検定を2群ずつ3回繰り返す」のは厳禁である。検定を繰り返すほど、全体としての偽陽性率（タイプIエラー）が膨らんでしまう（多重性の問題）。
これを解決するのが、survminer パッケージに用意されている pairwise_survdiff() 関数である。この関数を使えば、Bonferroni（ボンフェローニ）法や Holm（ホルム）法といったP値の補正を自動で行い、適切な事後多重比較を一発で実行できる。
あらかじめ必要なパッケージを読み込んでおこう。
library(tidyverse)
library(survival)
library(survminer)
事後多重比較の実践コード
# サンプルデータの作成（3群の生存時間データ）
set.seed(456)
n <- 150
df_surv <- tibble(
group = factor(sample(c("A", "B", "C"), n, replace = TRUE)),
time = rexp(n, rate = case_when(group == "A" ~ 0.05, group == "B" ~ 0.08, group == "C" ~ 0.12)),
status = sample(c(0, 1), n, replace = TRUE, prob = c(0.2, 0.8))
)
# 1. 全体のログランク検定（Global Log-rank Test）
global_diff <- survdiff(Surv(time, status) ~ group, data = df_surv)
print(global_diff) # ここで有意差が出た場合、事後多重比較に進む
# 2. pairwise_survdiff による事後多重比較（Holm法でP値を補正）
pairwise_diff <- pairwise_survdiff(
Surv(time, status) ~ group,
data = df_surv,
p.adjust.method = "holm" # "bonferroni" なども指定可能
)
print(pairwise_diff)
pairwise_survdiff() の出力はマトリクス形式で表示され、「A vs B」「A vs C」「B vs C」それぞれの補正済みP値が一目で確認できる。論文に結果を記載する際は、全体のログランク検定のP値とともに、この補正済みP値を明記すれば統計的な厳密性を担保できる。
Cox比例ハザードモデルでの「基準群（Reference）」の設定とハザード比の読み方
多変量解析としてCox比例ハザードモデル（coxph）に3群の変数を投入する場合、Rは内部的に自動でダミー変数を作成して解析を行う。
この際、デフォルトでは因子の文字コード順（アルファベット順や五十音順）で一番若い群（上記の例では「A群」）が基準群（Reference）として固定される。出力されるハザード比（HR）は、「A群に対するB群のハザード比」「A群に対するC群のハザード比」の2つとなる。
しかし、実務においては「標準治療であるC群を基準（Reference）にして、新規治療であるA群やB群の優越性をハザード比で示したい」というケースが多々ある。
意図した群を基準に設定するためには、モデルに投入する前に因子の水準（Factor Level）を適切に並べ替えておく必要がある。
因子水準の変更方法や、回帰モデルにおける基準群の切り替えコード、さらには作図時の軸順序の制御といった詳細なロジックについては、以前に執筆した解説記事に網羅されているため、実装の際は併せて参照してほしい。
➔ Rでの因子水準（factor）の並べ替え・基準群（Reference）の変更・ggplot2の軸順序の変更
意図通りに因子水準を再定義（Relevel）した上で coxph() を実行すれば、望み通りのハザード比と95%信頼区間がスッキリと得られる。
論文掲載レベル！3群の生存曲線（KM曲線）を ggsurvplot で美しく整える
3群以上の生存時間解析では、カプランマイヤー（KM）曲線の視覚的な分かりやすさも極めて重要である。しかし、3本の曲線と3段の Risk Table（各時点での生存人数一覧）をデフォルトのまま並べると、文字や線が重なって非常に見づらくなる。
survminer::ggsurvplot() の実務的なオプションを駆使し、論文にそのまま載せられるクオリティまでグラフを整形するテンプレートコードを以下に示す。
# 生存曲線オブジェクトの作成
fit <- survfit(Surv(time, status) ~ group, data = df_surv)
# 論文用カスタムKM曲線の描画
ggsurvplot(
fit,
data = df_surv,
pval = TRUE, # 全体のログランク検定のP値を表示
pval.coord = c(0, 0.1), # P値の表示位置を調整（x, y）
conf.int = FALSE, # 3群で信頼区間を描くと重なって見づらいためFALSE推奨
palette = c("#E41A1C", "#377EB8", "#4DAF4A"), # 識別しやすい配色
censor.shape = 124, # センサリングマークを縦線「|」にしてすっきりさせる
censor.size = 3,
# Risk Table（生存人数表）のカスタマイズ
risk.table = TRUE, # Risk Tableを表示
risk.table.col = "strata", # 文字色を群のカラーに合わせる
risk.table.y.text = FALSE, # 縦軸の群名を非表示にしてスペースを確保
tables.theme = theme_cleantable(),# テーブルの背景や枠線をクリーンに
# 全体レイアウト
ggtheme = theme_bw(), # 白背景の標準テーマ
xlab = "Time (Months)", # 横軸ラベル
ylab = "Overall Survival", # 縦軸ラベル
legend.labs = c("Group A", "Group B", "Group C"), # 凡例名
legend.title = "Treatment"
)
この設定を適用することで、3つの群の情報が整理され、査読時にも視認性が高く洗練された印象を与える図が完成する。
比例ハザード性は「グラフ（残差プロット）」で確認するのが鉄則
Cox比例ハザード回帰を行う際、前提条件となる「比例ハザード性（時間の経過によってハザード比が変化しないこと）」の確認は必須である。
なお、比例ハザード性の確認において「仮定の仮説検定（P値の判定）」に頼りすぎるのは好ましくない。帰無仮説（比例ハザードである）が棄却されなかった（有意にならなかった）からといって、「比例ハザード性が証明された」わけではなく単なる判断保留に過ぎないためである（サンプルサイズが小さいと不適格でも有意になりにくく、逆に大サンプルではわずかなズレでも有意になってしまう問題がある）。
したがって、比例ハザード性は「Schoenfeld残差プロット」および「KM曲線」のグラフで視覚的に確認するのが最も理にかなったアプローチである。
# 1. Cox回帰モデルの構築
fit_cox <- coxph(Surv(time, status) ~ group, data = df_surv)
# 2. Schoenfeld残差の計算
cz <- cox.zph(fit_cox)
# 3. Schoenfeld残差プロットの描画（グラフで視覚的に確認）
ggcoxzph(cz)
グラフによる視覚的判定のポイントと対処法
- 判定ポイント: 描画された平滑化ライン（スプライン曲線）が水平（y = 0 のラインに平行）に近く、時間経過とともに直線から大きく傾いていなければ、比例ハザード性が概ね保持されていると判断してよい。赤い点は残差そのもの。縦軸は時間におけるハザード比で、これが一定（横一線＝比例ハザード性）であることが求められている。
- 対処法: もし残差プロットのラインが右肩上がり／右肩下がりに大きく傾いている場合や、KM曲線のライン同士が途中で交差（Cross）している場合は、比例ハザード性の破綻が疑われる。その場合は、観察期間を前期・後期に分割して解析するランドマーク解析（Landmark Analysis）や、時間依存性共変量をモデルに組み込むアプローチへの変更を検討しよう。
【まとめ】コピペで動く！3群生存時間解析の一括Rコードテンプレート
本記事で解説した「データ作成 ➔ 全体検定 ➔ 事後多重比較 ➔ 基準群設定付きCox回帰 ➔ 美しいKM曲線描画 ➔ グラフによる仮定確認」の一連のパイプラインコードを以下にまとめた。自身の手元のデータセット名・変数名に置き換えて活用してほしい。
# === 3群以上の生存時間解析一括テンプレート ===
library(tidyverse)
library(survival)
library(survminer)
# 1. 全体のログランク検定
global_diff <- survdiff(Surv(time, status) ~ group, data = df_surv)
print(global_diff)
# 2. 事後ペアワイズ多重比較（Holm法で補正）
pairwise_diff <- pairwise_survdiff(Surv(time, status) ~ group, data = df_surv, p.adjust.method = "holm")
print(pairwise_diff)
# 3. 基準群の並べ替え（※詳細は過去記事を参照）
# df_surv <- df_surv %>% mutate(group = fct_relevel(group, "C"))
# 4. Cox比例ハザード回帰
fit_cox <- coxph(Surv(time, status) ~ group, data = df_surv)
summary(fit_cox)
# 5. 比例ハザード性の視覚的確認（Schoenfeld残差プロット）
cz <- cox.zph(fit_cox)
ggcoxzph(cz)
# 6. 論文用カプランマイヤー曲線の描画
fit_km <- survfit(Surv(time, status) ~ group, data = df_surv)
ggsurvplot(
fit_km, data = df_surv, pval = TRUE, palette = "Set1",
risk.table = TRUE, risk.table.y.text = FALSE, tables.theme = theme_cleantable(),
ggtheme = theme_minimal()
)
2群比較で培った生存時間解析の知識を3群以上の多群比較へと正しく拡張し、多重性の罠や基準群の迷子を回避し、グラフによる適切な仮定確認を行うことで、査読に耐えうる頑健で美しい臨床研究論文を完成させよう。
コメント