無作為化比較試験の主解析で、治療のハザード比が0.65だったとします。層別因子で調整した多変量Coxモデルにかけ直すと0.59になりました。共著者から「割付は無作為なのに、なぜ調整で値が動くんですか」と聞かれます。
交絡が残っていたのではありません。ハザード比という指標自体が、調整する変数を変えると値が変わる性質を持っています。
ここではCoxモデルが基準ハザードを消す仕組みを確認したうえで、coxphの出力をどう読むか、そして「調整済み」という言葉が何を意味しているのかを見ていきます。
基準ハザードを消すという発想
Coxモデルは、共変量ベクトル$\mathbf{x}$を持つ個体のハザードを次の形に置きます。
$$h(t \mid \mathbf{x}) = h_0(t) \exp(\boldsymbol{\beta}^\top \mathbf{x}) \tag{1}$$
式(1)は、時間の関数である$h_0(t)$と、共変量の効果である$\exp(\boldsymbol{\beta}^\top\mathbf{x})$の掛け算という構造です。時間を含むのは前の因子だけで、共変量の効果に$t$は入っていません。だから2個体のハザード比をとると$h_0(t)$が約分され、時間に依存しない定数が残ります。これが比例ハザードという仮定の中身です。
この構造が効いてくるのは推定のときです。イベントが起きた時刻$t_{(j)}$で、そのときのリスク集合$\mathcal{R}(t_{(j)})$の中から実際にイベントを起こしたのがその個体だった条件付き確率を書くと、分子と分母の両方に$h_0(t_{(j)})$が現れて消えます。
$$L_j(\boldsymbol{\beta}) = \frac{\exp(\boldsymbol{\beta}^\top \mathbf{x}_{(j)})}{\sum_{l \in \mathcal{R}(t_{(j)})} \exp(\boldsymbol{\beta}^\top \mathbf{x}_l)} \tag{2}$$
式(2)には$h_0(t)$がどこにも残っていません。基準ハザードの形を何も決めないまま、回帰係数だけを推定できます。これをすべてのイベント時点にわたって掛け合わせたものが部分尤度で、その理論的な正当化は部分尤度で扱います。
実務上の含意は単純です。ハザードの時間的な形について何も仮定しなくてよい代わりに、比例ハザードという仮定を1つ背負います。パラメトリック生存モデルとは、どの仮定を負うかの選択が違うだけです。
coxphの出力を読む
lungデータで、性別、年齢、ECOG全身状態指標を入れた多変量モデルを当てはめます。欠測を除いた227例、死亡164件です。
library(survival)
analysis <- na.omit(lung[, c("time", "status", "sex", "age", "ph.ecog")])
fit <- coxph(Surv(time, status) ~ factor(sex) + age + factor(ph.ecog),
data = analysis)
summary(fit)#> exp(coef) lower .95 upper .95 Pr(>|z|)
#> factor(sex)2 0.5794 0.4166 0.8056 0.00118
#> age 1.0109 0.9926 1.0295 0.24637
#> factor(ph.ecog)1 1.5069 1.0190 2.2284 0.03995
#> factor(ph.ecog)2 2.4677 1.5782 3.8587 7.48e-05
#> factor(ph.ecog)3 7.0607 0.9383 53.1289 0.05767
#>
#> Likelihood ratio test = 30.87 on 5 df, p = 1e-05
読み方を順に確認します。性別の0.579は、年齢とECOG指標を固定したときの女性対男性のハザード比です。年齢の1.011は1歳あたりの比なので、10歳差なら$1.011^{10} = 1.11$と計算します。ECOG指標はダミー変数なので、それぞれECOG0を基準とした比です。
ここで注意が要るのはECOG3の行です。ハザード比7.06、信頼区間0.94から53.13。区間の幅が50倍以上あります。このデータでECOG3の患者は1名しかいません。点推定値だけを表に載せると、読者は7倍という数字を受け取ってしまいます。私なら、この行はECOG2と統合するか、統合しない理由を本文に書きます。
共変量をいくつまで入れてよいかの目安として、イベント数が共変量あたり10以上という基準がよく参照されます。今回は164件に対してパラメータ5つなので32.8で、余裕があります。ただしこの基準はダミー変数の1つ1つを数えるので、カテゴリ変数を細かく分けるとすぐに枯渇します。ECOG3のような症例数の少ないカテゴリで区間が爆発するのは、この基準を満たしていても起こります。
調整しても交絡が消えるだけではない
冒頭の疑問に戻ります。無作為化されていれば共変量は群間で釣り合っているので、調整しても値は動かないはずだ。ロジスティック回帰でも同じことが起きますが、線形回帰の感覚で考えるとこうなります。
実際に確かめてみます。治療を無作為に割り付け、それとは完全に独立な予後因子を1つ置いた架空データを作ります。真の条件付きハザード比は0.60です。以下は説明のための架空データです。
set.seed(42)
n <- 4000
treatment <- rbinom(n, 1, 0.5) # 無作為割付
prognostic <- rnorm(n) # 割付と独立な予後因子
risk <- log(0.60) * treatment + log(2.5) * prognostic
event_time <- rexp(n, rate = 0.05 * exp(risk))
cens_time <- runif(n, 0, 40)
obs <- pmin(event_time, cens_time)
died <- as.numeric(event_time <= cens_time)
exp(coef(coxph(Surv(obs, died) ~ treatment)))
exp(coef(coxph(Surv(obs, died) ~ treatment + prognostic))[1])#> 調整なし 0.646
#> 調整あり 0.586群間の予後因子の平均は0.013と-0.022で、釣り合っています。交絡はありません。それでも調整なしの推定値は0.646、調整ありは0.586で、真値の0.60に近いのは調整したほうです。
何が起きているか。予後の悪い患者から先にイベントを起こして脱落するので、時間が経つほど両群のリスク集合には予後のよい患者が残ります。この選択は予後の悪い側の群でより速く進むため、群間のハザードの差が見かけ上縮まります。調整しない推定値は、この時間を通じた集団構成の変化を平均した値になります。
ハザード比が持つこの性質を非可縮性（non-collapsibility）と呼びます。リスク差やリスク比にはない性質で、オッズ比とハザード比だけが持ちます。実務上の含意は3つに分かれます。
第一に、粗ハザード比と調整済みハザード比が違っていても、それだけでは交絡の証拠になりません。「調整で値が動いたので交絡があった」と考察に書くのは誤りです。第二に、どの共変量を入れたかによって推定対象そのものが変わるので、調整した変数の一覧を必ず明記します。同じデータでも変数の組が違えば比較できません。第三に、無作為化試験で事前に規定した調整モデルがあるなら、そちらが主解析です。粗解析との差は交絡の指標ではなく、条件付き効果と周辺効果の違いとして説明します。
モデルから生存曲線を出す
ハザード比は比の指標なので、それだけでは患者の生存率が分かりません。共変量の値を決めれば、Coxモデルから生存曲線を描けます。
profile <- data.frame(sex = c(1, 2, 1, 2),
age = median(analysis$age),
ph.ecog = c(0, 0, 2, 2))
plot(survfit(fit, newdata = profile), xlab = "Days", ylab = "S(t)")
ここで区別が要るのは、この曲線が「特定の共変量値を持つ患者の生存曲線」だということです。集団全体の生存曲線ではありません。集団の効果を知りたい場合は、共変量の分布全体にわたって平均した周辺生存曲線を計算します。論文でCoxモデルの生存曲線を出すときは、どちらなのかを図の説明に書きます。年齢を中央値に固定した曲線を「治療群の生存曲線」と書くと、読者は集団全体の曲線だと受け取ります。
結果を書くときに決めること
調整する共変量は、データを見る前に決めます。単変量でp値が小さかった変数を集めて多変量モデルを作る手順は勧めません。選択の過程が推定値と信頼区間に反映されず、区間が実際より狭く出るためです。臨床的に交絡が疑われる変数と、事前に規定した層別因子を入れる。これが出発点になります。
報告には、調整済みハザード比と95%信頼区間、そして調整した共変量の一覧を書きます。イベント数とパラメータ数も添えておくと、読者が推定の安定性を判断できます。比例ハザード性の確認結果も本文か補足資料に置きます。
その比例ハザード性が崩れていた場合にどうするかは比例ハザード仮定の検証、変数選択と外れ値の扱いはCoxモデルの変数選択と診断、時間とともに値が変わる共変量を入れる場合は時変共変量で扱います。
無作為化試験で事前規定の調整モデルをどう組むか、観察研究で交絡変数をどこまで入れるか。このあたりは研究デザインと臨床的な文脈の両方を見ないと決まりません。Dr.データサイエンスでは、こうした解析計画のご相談を承っています。
参考文献
Cox DR. Regression models and life-tables. Journal of the Royal Statistical Society: Series B. 1972;34(2):187-202.
Martinussen T, Vansteelandt S. On collapsibility and confounding bias in Cox and Aalen regression models. Lifetime Data Analysis. 2013;19(3):279-296.
Hernán MA. The hazards of hazard ratios. Epidemiology. 2010;21(1):13-15.
Peduzzi P, Concato J, Feinstein AR, Holford TR. Importance of events per independent variable in proportional hazards regression analysis. Journal of Clinical Epidemiology. 1995;48(12):1503-1510.
Therneau TM, Grambsch PM. Modeling Survival Data: Extending the Cox Model. Springer; 2000.
