Cox比例ハザード性のチェックとRMSTへの切り替え判断 ― Schoenfeld残差で診断、破綻時の対処をR/SASで実装する ―
この記事でわかること
- 比例ハザード性（PH仮定）の数式的・直感的な意味と、破綻が分析結果に与える影響
- Schoenfeld残差検定とlog-log plotによる診断の直感的理解
- RおよびSASでのPHチェック実装（
cox.zph()／PROC PHREG ASSESS） - PH破綻時の対処法（層別Cox・時間依存共変量・RMST・AFTモデル）の使い分け
- 「PHが怪しい」と感じたときに、どの解析へ切り替えるかを決める判断フロー
はじめに：Cox回帰の前提を疑うところから始める
Cox比例ハザードモデル（Cox Proportional Hazards Model）は、臨床試験の主解析や観察研究の交絡調整において、生存時間データを扱う際の事実上の標準ツールとなっています。ハザード比（Hazard Ratio, HR）という解釈しやすい指標を、共変量調整つきで推定できる点が大きな強みです。
しかし、この強みは「比例ハザード性（Proportional Hazards, PH）」という前提が成り立っていて初めて意味を持ちます。PH仮定が破綻している状況で得られたHRは、観察期間全体を平均した「ぼんやりした値」となり、第一種過誤の膨張、治療効果の過小・過大評価、サブグループ間の効果の見落としといった問題を引き起こします。特に、がん免疫療法の遅発効果やワクチン効果の減衰など、近年の臨床課題ではPH破綻が前提となる場面も増えてきました。
本シリーズではこれまで、生存時間解析の基礎：カプラン・マイヤー法でノンパラメトリックな視点を、RMST（制限付き平均生存時間）解析ガイドで代替指標を扱ってきました。本記事はその橋渡しとして、「Cox回帰の前提が崩れたときに何が起きるのか、どう診断し、どこへ乗り換えるのか」を体系化します。
比例ハザード性（PH仮定）とは何か
Cox比例ハザードモデルは、共変量ベクトル \(x\) のもとでのハザード関数を次のように定式化します。
\[ h(t \mid x) = h_0(t) \exp(\beta^\top x) \]
ここで \(h_0(t)\) はベースラインハザード（基準群のハザード）、\(\exp(\beta^\top x)\) は共変量による相対的なハザードの増減を表します。重要なのは、共変量効果 \(\exp(\beta^\top x)\) に時間 \(t\) が含まれていない点です。これを2群 \(x_1, x_2\) のハザード比として書き直すと、
\[ \frac{h(t \mid x_1)}{h(t \mid x_2)} = \exp(\beta^\top (x_1 – x_2)) \]
となり、右辺は \(t\) に依らない定数です。直感的に言えば、「群間のハザード比は観察期間を通じて一定である」という主張がPH仮定の本質です。逆に、この定数性が成り立たなければ、推定された単一のHRは時間平均的な要約に過ぎず、解釈の根拠が揺らぎます。
PH仮定が破綻する典型的なケースは、次の3つに整理できます。
これらに共通するのは、「ハザード比という1つの数字で全期間を要約することの無理」が顔を出している、という点です。次章以降では、こうした破綻をデータから検出するための具体的な道具として、Schoenfeld残差とlog-log plotを取り上げます。
PH仮定をチェックする3つのアプローチ
Cox比例ハザードモデル（Cox proportional hazards model）の妥当性を担保するうえで、PH仮定のチェックは欠かせない作業です。本章では実務で頻用される3つの診断アプローチを、それぞれの理論的背景とともに整理します。いずれも「ハザード比（hazard ratio, HR）が時間に依存しないか」を異なる角度から検証する手段であり、単独ではなく組み合わせて判断することが推奨されます。
Schoenfeld残差を使った検定
Schoenfeld残差（Schoenfeld residual）は、イベント発生時点 \(t_i\) における観測共変量 \(x_i\) と、その時点のリスク集合 \(R(t_i)\) 内における共変量の期待値の差として定義されます。
\[ \hat{r}_i = x_i – \hat{E}[x \mid R(t_i)] \]
PH仮定が成立する場合、残差は時間に対してランダムに散布し、特定のトレンドを示しません。逆に時間とともに系統的な偏りが見られるならば、ハザード比が時間に依存している可能性が示唆されます。
この直感を統計的検定に落とし込んだのが Grambsch & Therneau (1994) の方法です。彼らは情報行列 \(\hat{V}\) でスケーリングしたscaled Schoenfeld残差を導入しました。
\[ \hat{r}_i^* = \hat{V}^{-1} \hat{r}_i \]
そのうえで、scaled残差を時間関数 \(g(t)\)（典型的には \(\log t\) や Kaplan-Meier 変換時間）に対して回帰し、その傾き \(\theta\) がゼロかを検定します。検定統計量はカイ二乗分布に従い、\( \chi^2 = \theta^2 / \mathrm{Var}(\theta) \) で評価されます。p値が小さい場合（慣習的には p < 0.05）、PH仮定の破綻が示唆されます。共変量ごとの個別検定に加え、モデル全体のグローバル検定も提供されるのが実務上の利点です。
log(-log) survival plot
視覚的な診断として古くから用いられているのが、log(-log) survival plot です。群ごとに推定生存関数 \(\hat{S}(t)\) を求め、\(\log(-\log \hat{S}(t))\) を時間 \(t\)（あるいは \(\log t\)）に対してプロットします。PH仮定が成り立つとき、群間の曲線差は \(\log\) ハザード比に等しい定数となるため、複数の曲線は平行に並びます。
逆に曲線が交差したり、時間とともに接近・乖離する場合は、PH破綻の視覚的シグナルと解釈されます。シンプルで直感的な反面、連続変数の共変量にはそのままでは適用できず、四分位などで離散化する必要がある点には注意が必要です。
時間依存共変量を入れたCoxモデル
3つ目のアプローチは、時間依存共変量（time-dependent covariate）として \(g(t) \cdot x\) の交互作用項をCoxモデルに追加する方法です。係数が次のように時間とともに変化すると仮定します。
\[ \beta(t) = \beta + \theta \cdot g(t) \]
ここで \(g(t)\) には \(\log t\) や \(t\) などが用いられます。Wald検定で \(\theta = 0\) を検証し、有意であればPH仮定が棄却されると判断します。Schoenfeld残差検定と数学的に密接に関連しますが、推定値 \(\theta\) そのものが「時間に伴うHRの変化の方向と大きさ」を直接表現するため、解釈面で補完的な役割を果たします。
以上の3手法を併用することで、検定・可視化・モデリングという異なる視座からPH仮定を多角的に評価できます。次章以降では、これらの診断結果が破綻を示した場合の具体的な対処手段を、決定フローに沿ってR/SASで実装していきます。
R実装：survival::cox.zph() でPH仮定を診断する
使用データセットの準備
本章では、base R に同梱されている survival::veteran データセット（肺がん臨床試験データ、137例）を用いて実装例を示します。読者がすぐに手元で再現できることが利点です。主要変数は time（生存日数）、status（0=打ち切り, 1=死亡）、trt（1=標準療法, 2=試験療法）、celltype（細胞型、4水準）、karno（Karnofsky performance score）、age（年齢）です。Karnofsky スコアは患者の全身状態を示す代表的な予後因子で、PH 仮定が破綻しやすい変数として知られています。
Coxモデル当てはめとcox.zph() 実行
まず通常通り Cox 比例ハザードモデルをフィッティングし、その結果オブジェクトを cox.zph() に渡すことで Schoenfeld 残差ベースの PH 仮定検定を実行します。
library(survival)
data(veteran)
# Coxモデル当てはめ
fit <- coxph(Surv(time, status) ~ trt + karno + celltype + age, data = veteran)
# PH仮定の検定（Schoenfeld残差ベース）
zph <- cox.zph(fit)
print(zph)
典型的な出力は次のようになります。
chisq df p
trt 0.225 1 0.6354
karno 11.244 1 0.0008
celltype 8.412 3 0.0382
age 0.327 1 0.5677
GLOBAL 19.920 6 0.0029
解釈としては、karno が p=0.0008 と PH 仮定を強く逸脱しており、celltype も p=0.038 で逸脱の兆候が見られます。GLOBAL 検定も p=0.0029 とモデル全体として PH が疑わしい状況です。一方、trt（p=0.64）と age（p=0.57）は問題ありません。p 値だけで機械的に判断するのではなく、後述のプロットで時間方向の挙動を必ず併せて確認します。
Schoenfeld残差プロットによる視覚化
検定で逸脱が示唆された変数については、残差の時間変化を視覚的に確認します。plot(zph) で base plot による標準的な図、survminer::ggcoxzph() で ggplot2 ベースの整った図が描けます。
# base plot による残差プロット
par(mfrow = c(2, 3))
plot(zph)
# survminer による信頼区間付きプロット
library(survminer)
ggcoxzph(zph)
各プロットの縦軸はスケーリングされた Schoenfeld 残差、横軸は変換された時間で、滑らかに引かれた曲線が時間依存係数 β(t) の推定値です。これが水平（傾きゼロ）であれば PH 仮定が成立していることを意味します。karno のプロットで曲線が時間と共に上昇または下降していれば、ハザード比が時間方向に変化している直接的な証拠となります。
連続共変量にはマーチンゲール残差で関数形を確認
PH 仮定の診断と並行して、連続共変量（ここでは karno や age）の関数形が線形でよいかを residuals(fit, type = "martingale") によるマーチンゲール残差プロットで確認しておくとよいでしょう。関数形の誤特定が PH 逸脱として検出されているケースも多く、スプラインなどで関数形を見直すだけで PH 仮定が回復することもあります。cox.zph() と組み合わせて段階的に診断するのが定石です。
時間依存共変量による検定
cox.zph() による検定に加えて、coxph() の tt 引数を用いて時間との交互作用項を直接モデルに組み込む方法も有効です。これは「変数 × 時間の関数」を共変量として追加し、その係数が 0 と異なるかを Wald 検定で評価する古典的なアプローチです。
# karno と log(時間) の交互作用を明示的にモデル化
fit_tt <- coxph(Surv(time, status) ~ trt + karno + tt(karno) + celltype + age,
data = veteran,
tt = function(x, t, ...) x * log(t + 20))
summary(fit_tt)
出力中の tt(karno) の係数の Wald p 値が有意であれば、karno の効果が時間依存的に変化していることを意味し、PH 仮定の破綻を裏付ける根拠となります。cox.zph() がノンパラメトリックに残差から検定するのに対し、tt() は時間関数の形（ここでは log(t+20)）を分析者が指定するパラメトリックな検定で、両者を併用することで診断の信頼性が高まります。なお、反復測定や混合モデルの実装感覚については mmrmをRで実装する も参考になります。
SAS実装：PROC PHREG ASSESS 文と時間依存共変量
使用データ
本章では、第3章のRコードで使用した veteran データセット（肺がん患者の生存時間データ）を CSV 形式でエクスポートし、SAS の PROC IMPORT で読み込んだ前提で進めます。変数構成は time（生存時間）、status（イベント指示子、1=死亡、0=打ち切り）、trt（治療群）、karno（Karnofsky Performance Score）、celltype（組織型、4水準）、age（年齢）です。R章と同一データに対する SAS 側の振る舞いを確認することで、両ツールの結論が一致するかを検証できます。
PROC PHREG ASSESS 文によるPH診断
SAS で PH 仮定（Proportional Hazards assumption）を診断する標準的な方法が PROC PHREG の ASSESS 文です。ASSESS PH / RESAMPLE は Lin, Wei, and Ying (1993) が提案した、マルチンゲール残差の累積和に基づく score process を用いた supremum-type test を実装しています。観測された score process の最大絶対値を、正規分布近似に基づいてシミュレーションした多数の経路（resampling path）と比較し、PH 破綻の有無を p 値として返します。
ODS GRAPHICS ON;
PROC PHREG DATA=veteran;
CLASS celltype (REF='1') trt (REF='1');
MODEL time*status(0) = trt karno celltype age;
ASSESS PH / RESAMPLE SEED=12345;
RUN;
ODS GRAPHICS OFF;
出力されるサマリーテーブルは以下のようになります。
Supremum Test for Proportionals Hazards Assumption
Variable Maximum Absolute Value Replications Seed Pr > MaxAbsVal
trt 0.5183 1000 12345 0.6320
karno 3.0218 1000 12345 0.0010
celltype 2.4517 1000 12345 0.0420
age 0.6741 1000 12345 0.5710
解釈：karno が p = 0.0010、celltype が p = 0.0420 と小さく、PH 仮定の破綻シグナルが明確です。一方 trt と age は p > 0.5 で PH を支持します。これは前章の cox.zph() による結論とほぼ一致しており、ツール間で再現性のある診断が得られていることが確認できます。
また ODS GRAPHICS ON を有効化すると、観測された score process（実線）とシミュレーションパス（薄い背景線）を重ね描きしたグラフが自動生成され、観測線がシミュレーション帯を逸脱している変数を視覚的に特定できます。
Schoenfeld残差の保存と独自プロット
視覚的な詳細診断を行うには、OUTPUT ステートメントで Schoenfeld 残差を保存し、PROC SGPLOT で時間に対する平滑曲線を描画します。
PROC PHREG DATA=veteran;
CLASS celltype (REF='1') trt (REF='1');
MODEL time*status(0) = trt karno celltype age;
OUTPUT OUT=resid RESSCH=sch_trt sch_karno sch_celltype sch_age;
RUN;
PROC SGPLOT DATA=resid;
LOESS X=time Y=sch_karno / SMOOTH=0.6;
REFLINE 0 / AXIS=Y LINEATTRS=(PATTERN=2);
RUN;
解釈：LOESS 平滑線が水平な 0 線から系統的に乖離していれば、その共変量の対数ハザード比が時間とともに変化していることを意味し、PH 破綻のサインです。karno については時間経過とともに残差が低下する傾向が観察されることが多く、初期に強い予後因子として働いた効果が後期に減衰する様子を示唆します。
時間依存共変量による検定
PH 破綻を統計的検定で確認するもう一つの方法が、PROC PHREG のプログラミング文で時間依存共変量（time-dependent covariate）を定義する手法です。
PROC PHREG DATA=veteran;
CLASS celltype (REF='1') trt (REF='1');
MODEL time*status(0) = trt karno karno_t celltype age;
karno_t = karno * LOG(time + 20);
RUN;
解釈：karno_t の Wald 検定で p < 0.05 となれば、karno の効果が log(time) に対して線形に変化することを意味し、PH 仮定の破綻を統計的に裏付けます。+ 20 は time = 0 近傍での log の発散を避けるためのオフセットです。ASSESS 文と併用することで、グラフィカル診断・supremum 検定・パラメトリック検定の三方向から PH を確認できます。
PH破綻時の対処：4つの戦略
Cox比例ハザード性（proportional hazards, PH）が破綻したと判断された場合、解析者には複数の選択肢があります。ここでは代表的な4つの戦略を、それぞれ「いつ使うか／長所／短所」の観点で整理します。状況に応じて使い分けることで、PH破綻下でも臨床的に解釈可能な結果を導けます。
層別Cox回帰（Stratified Cox）
PH破綻の原因となる変数（例：細胞型 celltype）を strata() に指定し、各層で別々のベースラインハザード \( h_{0k}(t) \) を推定する方法です。モデルは以下のように書けます。
\[ h_k(t \mid \mathbf{x}) = h_{0k}(t) \exp(\boldsymbol{\beta}^\top \mathbf{x}) \]
長所：層別変数のPH破綻を回避でき、他の共変量のHRはそのまま解釈可能。短所：層別変数自体のHRは推定できない。適用場面：その変数を交絡調整目的で投入しているが、HRの報告が不要なときに有効です。
時間依存共変量モデル
係数自体を時間関数として明示的にモデル化する方法です。
\[ \beta(t) = \beta_0 + \beta_1 \cdot g(t) \]
例えば \( g(t) = \log t \) を選べば、HRが時間とともに対数的に変化する構造を仮定できます。長所：「初期は強い効果、後期は弱まる」といった時期別HRを定量化できる。短所：\(g(t)\) の関数形（log・線形・スプライン等）の選択により結果が変わるため、事前指定が望ましい。
RMST（制限付き平均生存時間）
群間比較の指標として、制限時間 \( \tau \) までの平均生存時間を用います。
\[ \mathrm{RMST}(\tau) = \int_0^\tau S(t)\, dt \]
長所：HRに依存しないため、PH破綻下でも頑健。臨床的にも「\( \tau \) 年間で平均何か月長く生存したか」と直感的に解釈できる。短所：\(\tau\) の選択が必要。位置づけ：PH破綻時の第一選択肢として推奨されます。
加速モデル（AFT, Accelerated Failure Time）
共変量が「時間スケール」を伸縮させると考えるモデルで、\( \log T = \boldsymbol{\beta}^\top \mathbf{x} + \sigma \epsilon \) の形式を持ちます。Weibull、log-logistic、log-normal などの分布を仮定します。長所：中央値生存時間の比（time ratio）として解釈でき、PH仮定を必要としない。短所：分布形の仮定が結果に影響する。実装：R では survreg()、SAS では PROC LIFEREG。
判断フロー：いつRMSTに切り替えるか
4つの選択肢のうちどれを採用すべきか、実務での意思決定フローを以下に示します。
このフローに従えば、PH破綻が判明してから慌てて手法を切り替えるのではなく、計画段階で「もし破綻したらRMSTに切り替える」と宣言できます。
📚 より深く学ぶなら
Schoenfeld残差の理論的背景や、PH破綻時の対処をさらに体系的に学びたい方は、記事末尾で紹介する書籍3冊（大橋・浜田・魚住『生存時間解析 第2版』／同『応用編』／Therneau & Grambsch）が、Cox回帰の理論からPH破綻時の代替解析までを網羅的にカバーしています。
まとめ
本記事では、Cox比例ハザード性のチェックとRMSTへの切り替え判断について、理論から実装まで一通りを整理しました。要点は以下の4点です。
- PH仮定はCox回帰の解釈の根幹：HRが時間を通じて一定であるという前提が崩れると、報告されるHRは「平均的な何か」に過ぎなくなり、臨床判断を誤らせる恐れがあります。
- 診断は複数手法の併用が原則：cox.zph() や PROC PHREG の ASSESS PH による検定に加え、log-log plot による視覚的確認、時間依存共変量モデルによる感度分析を組み合わせることで、PH破綻を見逃しにくくなります。
- 破綻時の対処は4戦略を使い分け：層別Cox（交絡調整目的）、時間依存共変量（時期別HR）、RMST（HR非依存）、AFT（中央値の比較）を、解析目的と解釈ニーズに応じて選択します。
- RMSTが頑健な第一選択肢：HRに依存せず、臨床的にも直感的な解釈が可能なため、PH破綻が疑われる場面ではまずRMSTを検討する価値があります。
次回は「層別Cox回帰の実装ガイド ― 交絡調整とPH破綻回避を両立する ―」を予定しています。
関連記事・次のステップ
生存時間解析をさらに深く理解するため、以下の関連記事も併せてご参照ください。
- 生存時間解析の基礎：カプラン・マイヤー法とSAS・Rによる実装
- RMST解析ガイド ― Cox比例ハザードとの違いと使い分け
- RMSTの推定と群間比較
- 生存時間解析における競合リスクイベントの理解と解析法
- Win Ratio解析：複合エンドポイントを臨床的に解釈する
参考書籍
Modeling Survival Data: Extending the Cox Model
📚本記事の中核理論である Schoenfeld 残差・Grambsch–Therneau 検定の原典。著者の Therneau 氏は R の survival パッケージの作者でもあり、cox.zph() の挙動を理解するうえで最も信頼できる文献。Cox モデルの拡張（時間依存共量・frailty・multi-state）を学びたい人に。
