解析の途中経過を共著者に見せたら、こう聞かれたとします。「生存曲線があるのに、なぜ累積ハザードのグラフも作っているんですか」。生存関数と累積ハザードは$S(t) = \exp(-H(t))$で一対一に対応しているので、片方があればもう片方は決まります。情報としては同じものです。
それでも累積ハザードを別に推定する場面があります。比例ハザード性を目で確かめるとき、再発するイベントを数えるとき、そしてCox回帰から個々の患者の生存曲線を出すときです。
ここではNelson-Aalen推定量の定義を確認したうえで、この3つの場面で累積ハザードのスケールが何をしてくれるのかを見ていきます。
累積ハザードは何を数えているか
イベントが観測された時刻$t_{(j)}$で、直前のリスク集合の大きさを$n_j$、イベント数を$d_j$とします。Nelson-Aalen推定量は各時点の$d_j/n_j$をそのまま足していきます。
$$\hat{H}(t) = \sum_{t_{(j)} \le t} \frac{d_j}{n_j} \tag{1}$$
式(1)は、各時点で「今いる人のうち何割にイベントが起きたか」を積み上げた量です。Kaplan-Meier推定量が同じ$d_j/n_j$を$1 – d_j/n_j$の形にして掛け算していくのに対し、こちらは引き算も掛け算もせずに足すだけです。
足すだけなので、$\hat{H}(t)$は1を超えます。生存確率と違って上限がありません。この上限のなさが、あとで出てくる3つの使い道のもとになります。
library(survival)
months <- c(3, 4, 5, 5, 8, 8, 9, 11, 12, 14)
event <- c(1, 0, 1, 0, 1, 1, 0, 1, 0, 0)
km <- survfit(Surv(months, event) ~ 1, stype = 1) # 積の形
fh <- survfit(Surv(months, event) ~ 1, stype = 2) # exp(-H)の形
cbind(time = km$time, H = km$cumhaz, KM = km$surv, FH = fh$surv)#> time H KM FH
#> [1,] 3 0.1000 0.9000 0.9048
#> [2,] 5 0.2250 0.7875 0.7985
#> [3,] 8 0.5583 0.5250 0.5722
#> [4,] 11 0.8917 0.3500 0.4100
累積ハザードから生存関数を戻したものをFleming-Harrington推定量と呼びます。$\log(1-x) < -x$が常に成り立つので、この推定量はKaplan-Meier推定量より必ず大きくなります。標本が大きくなれば差は消えますが、10例のこのデータでは6ポイントありました。論文でどちらを使ったかを書いておかないと、他の人が再現したときに数値が合いません。
比例ハザード性を目で確かめる
Cox回帰にかける前に、群間のハザード比が追跡期間を通じて一定と見なせるかを確認します。ここで累積ハザードのスケールが効いてきます。
比例ハザードが成り立つなら、2群の累積ハザードは定数倍の関係になります。両辺の対数をとると、定数倍が上下の平行移動に変わります。
$$\log H_1(t) = \log \mathrm{HR} + \log H_0(t) \tag{2}$$
式(2)は、$\log \hat{H}(t)$を時間に対して描いたとき、比例ハザードなら2本の曲線が一定の縦幅を保つことを表しています。$H(t) = -\log S(t)$なので、この縦軸は$\log(-\log \hat{S}(t))$と同じものです。Rではfun引数にcloglogを渡すだけで描けます。
plot(survfit(Surv(time, status) ~ sex, data = lung),
fun = "cloglog", xlab = "Time (log scale)", ylab = "log(-log S(t))")
veteran$kg <- ifelse(veteran$karno >= 70, "high", "low")
plot(survfit(Surv(time, status) ~ kg, data = veteran), fun = "cloglog")
この2つを検定にかけると、lungの性別はp=0.09、veteranのKarnofsky指標はp=0.0003でした。図の見た目と検定結果が一致しています。
ただし、私ならこのプロットを検定の代わりには使いません。曲線の左端は追跡開始直後でイベントが数件しかなく、$\log(-\log \hat{S}(t))$が大きく振れます。図2でも左端の階段が粗いのはそのためで、ここを見て非平行と判断すると誤ります。プロットで形を把握し、検定で判断する。順序としてはこちらを勧めます。形を見ておくと、検定が有意になったときに「どういう非比例性なのか」を答えられます。この続きは比例ハザード仮定の検証で扱います。
再発するイベントを数える
COPDの増悪、心不全の再入院、てんかん発作、痛風発作。同じ患者に何度も起きるイベントを主要評価項目にする研究では、生存関数という枠組み自体が合いません。1回目の増悪で追跡が終わるわけではないからです。
最初の1回だけを見てKaplan-Meier曲線を描く方法はよく使われますが、2回目以降の情報を捨てています。年に4回増悪する患者と1回の患者が、初回までの時間が同じなら同じ扱いになります。
累積ハザードのスケールなら、この制約がありません。$\hat{H}(t)$は上限を持たないので、「1人あたり平均何回起きたか」をそのまま表せます。この形で使うときは平均累積関数と呼ばれます。
mcf <- survfit(Surv(tstart, tstop, status) ~ arm, data = episodes, id = id)
summary(mcf, times = c(90, 180, 365))#> arm time MCF
#> Placebo 90 0.44
#> Placebo 180 0.96
#> Placebo 365 1.98
#> Treatment 90 0.31
#> Treatment 180 0.56
#> Treatment 365 1.21
縦軸の単位が回数なので、群間差をそのまま引き算できます。上の例なら1年あたり0.77回の減少です。生存確率の差では、こういう読み方はできません。
ここで注意が要るのは、死亡が競合する場合です。心不全の再入院を数えている途中で患者が亡くなると、その後は再入院が起こりえません。死亡を単なる打ち切りとして扱うと、平均累積回数を過大に評価します。この構造は競合イベントの扱いと同じで、死亡を別の事象として組み込む定式化が必要になります。予後の悪い群ほど早く亡くなり、そのぶん再入院の機会も失う。この効き方があるので、群間差が縮まる方向にも広がる方向にも動きます。
累積ハザードの数値をどう書くか
3つ目の使い道はCox回帰の内側にあります。部分尤度で回帰係数を求めたあと、ベースラインの累積ハザードをデータから推定する必要があります。ここで使われるBreslow推定量は、式(1)の分母をリスクスコアの合計に置き換えた形をしています。共変量の値を決めれば、その患者の生存曲線が出てきます。予測モデルの較正を見るときも、この推定値が土台になります。
報告のときに困るのは、累積ハザードの値そのものが直感的でないことです。$\hat{H}(3\,\text{年}) = 0.35$と書かれても、臨床的に多いのか少ないのか分かりません。
対応としては、生存確率に換算した値を併記します。$\hat{H} = 0.35$なら$\exp(-0.35) = 0.70$で、3年生存率70%です。$\hat{H} = 1.0$は$\exp(-1.0) = 0.37$にあたります。
もう1つ、100人あたりの期待イベント数として$100 \times \hat{H}(t)$を出す書き方があります。ただしこの読み替えが使えるのは$\hat{H}(t)$が小さいときだけです。$\hat{H} = 0.05$なら100人あたり約5件でよいのですが、$\hat{H} = 1.0$を100人あたり100件と書くと、実際には63人にしか起きていないので誤りになります。イベントが1回しか起きない設定では、$\hat{H}$が0.2を超えたあたりからこの近似は使わないほうがよいでしょう。再発イベントの平均累積関数として使う場合は、$\hat{H}$がそのまま回数なので、この制約はかかりません。
どちらを図に出すか
読者に生存確率を見せたいなら生存曲線です。論文の主要な図はほぼこちらになります。累積ハザードの図を本文に載せるのは、再発イベントの平均累積回数を報告する場合と、比例ハザード性の議論を本文でする場合。それ以外は補足資料に置けば足ります。
分布の形を仮定して当てはめる方向に進むならパラメトリック生存モデル、共変量を入れた比較に進むならCox比例ハザードモデルが続きです。
再発イベントを主要評価項目にした研究で、死亡をどう組み込むか。予測モデルのベースライン生存関数をどう報告するか。このあたりは研究デザインとエンドポイントの定義を見ないと決まりません。Dr.データサイエンスでは、こうした解析方針のご相談を承っています。
参考文献
Nelson W. Theory and applications of hazard plotting for censored failure data. Technometrics. 1972;14(4):945-966.
Aalen O. Nonparametric inference for a family of counting processes. The Annals of Statistics. 1978;6(4):701-726.
Fleming TR, Harrington DP. Counting Processes and Survival Analysis. Wiley; 1991.
Rogers JK, Pocock SJ, McMurray JJV, et al. Analysing recurrent hospitalizations in heart failure. European Journal of Heart Failure. 2014;16(1):33-40.
Therneau TM, Grambsch PM. Modeling Survival Data: Extending the Cox Model. Springer; 2000.
