同じデータセットを2人が別々に解析して、ハザード比が0.59と0.60になったとします。データも変数も同じ、追跡期間も同じ。片方はR、もう片方はSASを使っていました。
どちらも間違っていません。同着、つまり同じ日にイベントが起きた患者が複数いる場合の処理方法が、既定値として違うだけです。RのcoxphはEfron法、SASのPROC PHREGはBreslow法を既定にしています。
生存時間が日単位や月単位で記録される臨床データでは、同着は例外ではありません。ここでは部分尤度が基準ハザードを消す仕組みを確認したうえで、同着の処理で推定値がどこまで動くか、そして推定がうまくいかないときに何が起きているかを見ていきます。
リスク集合の中での競争として書く
イベントが起きた時刻$t_{(j)}$で、そのとき追跡中だった患者の集合をリスク集合$R(t_{(j)})$と呼びます。部分尤度は、この集合の中で誰かにイベントが起きたという条件のもとで、実際にそれが個体$j$だった確率を書き下したものです。
$$L_j(\boldsymbol{\beta}) = \frac{\exp(\boldsymbol{\beta}^\top \mathbf{x}_j)}{\sum_{i \in R(t_{(j)})} \exp(\boldsymbol{\beta}^\top \mathbf{x}_i)} \tag{1}$$
式(1)の分子と分母の両方に基準ハザード$h_0(t_{(j)})$が掛かるので、約分されて消えます。時刻そのものは、リスク集合を決めるためだけに使われます。イベントが何日目に起きたかという情報は捨てて、誰が先に起きたかという順序だけを使う。これが部分尤度の構造です。
対数をとって微分すると、スコア関数の形が推定の中身を教えてくれます。
$$U(\boldsymbol{\beta}) = \sum_{j=1}^{D} \left[ \mathbf{x}_j – \bar{\mathbf{x}}(t_{(j)}, \boldsymbol{\beta}) \right] \tag{2}$$
式(2)は、実際にイベントを起こした患者の共変量と、そのときリスク集合にいた人たちの重み付き平均との差を、全イベント時点で足し合わせたものです。これがゼロになる$\boldsymbol{\beta}$を探します。イベントを起こした人が平均的な人と変わらなければ差はゼロに近く、係数もゼロに近づきます。
同着が増えると推定値が動く
式(1)は「その時点でイベントを起こしたのは1人」という前提で書かれています。同じ日に3人が亡くなったとき、この式をそのまま使うことはできません。
Breslow法は、3人それぞれについて同じ分母を3回使います。実際には1人目が起きた直後のリスク集合には2人目と3人目が残っていて、2人目の時点では3人目だけが残っているのですが、その減り方を無視します。分母を大きめに見積もることになるので、係数は1に向かって縮みます。
Efron法は、この減り方を平均的に補正します。同着した3人がどの順で起きたか分からないので、リスク集合から少しずつ抜いていく形で分母を調整する。計算量はBreslow法とほとんど変わりません。Exact法は3人の順序の全通りを数え上げるもので、正確ですが同着数が増えると計算が破綻します。
どのくらい違うのかを実データで確かめます。lungデータの観測時間を、日単位から週、月、四半期、半年へと丸めていき、同着の割合を変えながら3つの方法で推定します。
library(survival)
analysis <- na.omit(lung[, c("time", "status", "sex", "age")])
for (unit in c(1, 30, 90, 180)) {
analysis$rounded <- ceiling(analysis$time / unit) * unit
hr <- sapply(c("breslow", "efron", "exact"), function(m)
exp(coef(coxph(Surv(rounded, status) ~ factor(sex) + age,
data = analysis, ties = m))[1]))
cat(unit, round(hr, 4), "\n")
}#> 記録単位 breslow efron exact
#> 1日 0.5990 0.5986 0.5982
#> 30日 0.5975 0.5907 0.5752
#> 90日 0.6088 0.5834 0.5431
#> 180日 0.6209 0.5663 0.4943
日単位で記録されていれば、どの方法でも実質的に同じ値になります。問題になるのは、外来受診の間隔でしかイベント時点が分からない研究です。年2回の定期検診で判定するアウトカムなら、同着はほぼ全部の時点で起きます。
Efron法を使うことを勧めます。計算コストがBreslow法とほぼ同じで、Exact法に近い値を返すからです。Rの既定はこれなので、Rを使っているなら何もしなくてかまいません。SASで解析する場合はties=EFRONを明示的に指定します。そして論文には、どの方法を使ったかを書きます。同着が少ないデータでは書いても書かなくても数値は変わりませんが、書いておけば読者が再現できます。
Wald検定が効果を見落とすとき
coxphの出力には3つの検定が並びます。Wald検定、尤度比検定、スコア検定。標本が大きければ3つとも同じ結論になるので、ふだんは意識しません。
食い違うのは、効果が非常に強いときです。極端な例として、片方の群でイベントが1件も起きなかったデータを作ります。以下は説明のための架空データです。
set.seed(23)
n <- 40
group <- rep(0:1, each = n / 2)
event_time <- rexp(n, rate = 0.05 * exp(log(0.02) * group))
cens_time <- runif(n, 0, 50)
obs <- pmin(event_time, cens_time)
died <- as.numeric(event_time <= cens_time)
table(group, died)
summary(coxph(Surv(obs, died) ~ group))#> group died=0 died=1
#> 0 6 14
#> 1 20 0
#>
#> Likelihood ratio test = 30.22 on 1 df, p=4e-08
#> Wald test = 0.00 on 1 df, p=1
#> Score (logrank) test = 26.07 on 1 df, p=3e-07
#>
#> Warning: Loglik converged before variable 1; coefficient may be infinite.Wald検定のp値は1です。係数の表だけを見た読者は「差はなかった」と受け取ります。実際には介入群で1件もイベントが起きておらず、これ以上ないほど強い効果が出ています。尤度比検定のp値は4e-08で、Rも係数が無限大になりうるという警告を出しています。
Wald統計量は、最尤点まわりで対数尤度が2次関数に近いことを前提にしています。図3のように片側へ平らに伸びる形では、この近似が崩れます。効果が強いほどWald統計量が小さくなるこの現象はHauck-Donner効果と呼ばれ、ロジスティック回帰でも同じことが起きます。
完全分離までいかなくても、イベント数が少なく効果が大きい状況では同じ方向のずれが出ます。私なら、Wald検定と尤度比検定のp値が桁違いに離れていたら、Wald側を疑います。信頼区間も、Wald型ではなく尤度比に基づくプロファイル区間を出したほうが安全です。Rではconfint関数がこれを計算します。
ただし、完全分離が起きている状況で信頼区間を出しても、片側は無限大に伸びます。この場合に本当に必要なのは統計的な工夫ではなく、なぜ片群でイベントがゼロなのかを確認することです。追跡期間が短すぎないか、群の定義がアウトカムと重複していないか。数値が壊れているときは、たいていデータか設計の側に理由があります。
推定できる範囲を先に見積もる
部分尤度が使う情報はイベント数で決まります。追跡している人数ではありません。1000例を5年追跡してもイベントが12件なら、推定に使える情報は12件ぶんです。共変量1つあたりのイベント数が10を下回ると係数が不安定になり、5を下回ると過学習が目に見えてきます。
症例数設計の段階でイベント数を見積もっておくと、あとで困りません。予定登録数と想定イベント率から追跡期間内のイベント数を出し、それを入れたい共変量の数で割る。この値が10を切るなら、共変量を減らすか追跡を延ばすかを先に決めます。解析段階で気づいても打つ手はほとんどありません。
共変量の数が症例数を超えるような設定では、通常の部分尤度は解を持ちません。罰則付きの推定に切り替えることになりますが、その場合ハザード比の信頼区間の解釈が変わるので、予測モデルとして使うのか効果推定として使うのかを先に決めておく必要があります。
報告に書くこと
同着の処理方法、イベント数と共変量の数、そして検定の種類。この3つを書いておけば、読者が数値を再現できます。ソフトウェアとバージョンも添えると確実です。
基準ハザードは部分尤度では推定されないので、生存曲線を出すには別途Breslow推定量を使います。この推定量の形はNelson-Aalen推定量の分母をリスクスコアの和に置き換えたものです。
モデルの定式化と調整済みハザード比の解釈はCox比例ハザードモデル、推定の前提である比例ハザード性の確認は比例ハザード仮定の検証、当てはまりと外れ値の診断はCoxモデルの変数選択と診断で扱います。
イベント数が足りない解析でどこまで共変量を入れるか、同着の多いデータでどの方法を選ぶか。このあたりはデータの記録粒度と研究の目的を見ないと決まりません。Dr.データサイエンスでは、こうした解析方針のご相談を承っています。
参考文献
Cox DR. Partial likelihood. Biometrika. 1975;62(2):269-276.
Efron B. The efficiency of Cox’s likelihood function for censored data. Journal of the American Statistical Association. 1977;72(359):557-565.
Breslow N. Covariance analysis of censored survival data. Biometrics. 1974;30(1):89-99.
Hauck WW, Donner A. Wald’s test as applied to hypotheses in logit analysis. Journal of the American Statistical Association. 1977;72(360):851-853.
Andersen PK, Gill RD. Cox’s regression model for counting processes: a large sample study. The Annals of Statistics. 1982;10(4):1100-1120.
