Rで行える生存解析を出来る限り初歩から解説してみる~ハザード率、cox比例ハザード、ログランク検定等~

今回は生存解析をやろうと思います。
生存解析に関する解説をしているブログは巷に溢れかえっていますが、他のブログよりもハザード率の解説を丁寧に行うことにより差別化を図っていければいいなと思っています。

ハザード率(Hazard Rates)
- 2種類のハザード率の定義
- ハザード率の仮定及び推定
  - ノンパラメトリック推定
  - パラメトリック推定
打ち切りデータにおけるハザード率と生存確率
cox 比例ハザードモデル
生存解析を学ぶ上で役に立つ参考書

ハザード率(Hazard Rates)

ハザード率という概念は、医療や保険の現場でよく利用されます。

2種類のハザード率の定義

具体例

仮に保険屋さんの例で考えてみます。
保険業界では保険契約者のことをPolicy holder と呼んだりします。下の表はある保険会社Aが抱えているPolicy holderから得られたデータです。

Age : 年齢
num : Policy holder の人数
D : Policy holder の内亡くなった人数

D=sample(1:20,length(30:90),replace=TRUE)
Do=order(D)
D=D[Do]
LifeTable=data.frame(Age=30:90,num=sample(100:500,length(30:90),replace=TRUE),D=D)
> LifeTable
   Age num  D
1   30 495  2
2   31 206  2
3   32 412  3
4   33 482  3
5   34 355  4
6   35 335  4
7   36 353  4
8   37 367  4
9   38 312  5
10  39 114  5
11  40 253  6
12  41 121  7
13  42 354  7
14  43 238  7
15  44 194  8
16  45 437  9
17  46 106  9
18  47 410  9
19  48 216  9
20  49 176  9
21  50 234 10
22  51 411 10
23  52 288 10
24  53 241 11
25  54 227 11
26  55 398 11
27  56 404 11
28  57 444 11
29  58 320 12
30  59 238 12
31  60 451 12
32  61 126 13
33  62 306 13
34  63 234 13
35  64 287 13
36  65 350 14
37  66 340 14
38  67 137 14
39  68 200 15
40  69 234 15
41  70 294 16
42  71 318 16
43  72 459 16
44  73 403 17
45  74 488 17
46  75 144 17
47  76 258 17
48  77 265 17
49  78 499 18
50  79 346 18
51  80 422 18
52  81 209 18
53  82 382 18
54  83 218 19
55  84 191 19
56  85 246 19
57  86 483 20
58  87 310 20
59  88 335 20
60  89 247 20
61  90 188 20

(当然データ生成過程からもわかるように架空のデータです。)

ここで保険屋さんは保険料を設定するためにはPolicy holderの生存確率を考えなければなりません。例えば、今30歳の人たちが70歳まで生きる確率はどのくらいか。などなど

この生存確率を考えたいというのが今回の解析のモチベーションです。

ハザード率とは何か

では、考えてみます。

まず、i歳で亡くなる確率、つまり言い換えると寿命がi年である確率を仮に $Pr(X = i)$ と置くことにします。

このように設定することでi歳以降で亡くなる確率、つまり少なくともi-1歳までは生存している確率 $Pr(X≧i)$ は次のようにわかりますね。

$Pr(X≧i) = \sum_{j≧i} Pr(X=j)$

ということは、ベイズの定理より少なくともi-1歳までは生きているが、i歳で亡くなってしまう確率 $Pr(X=i | X≧i)$ もわかります。

$Pr(X=i | X≧i) =\frac{Pr(X = i)}{Pr(X≧i) }= \frac{Pr(X = i)}{\sum_{j≧i} Pr(X=j)}$

この $Pr(X=i | X≧i)$ がハザード率です。
ちなみに、今は暗黙的に離散時間でハザード率を考えていることに注意してください。

ハザード率から見える生存確率

更にハザード率がわかると、少なくともi-1歳まで生きて、i歳時点では亡くならない確率 $1-Pr(X=i | X≧i)$ もわかるわけです。

よって、少なくともi-1歳まで生きて、更にそこからj歳まで生きる確率 $Pr(X>j|X≧i)$ は次のような確率の掛け算で得られます。

$Pr(X>j|X≧i) = \prod_{k=i}^j (1-Pr(X=k | X≧k))$

ハザード率を連続時間に拡張する

知っているハザード率の定義と違う！と思った方少々お待ちください。
今の話をそのまま連続時間に拡張すると一般的なハザード率の定義になります。

正の連続値確率変数 $T$ を考えます。
この時、 $T$ が $t$ 以上となる確率 $Pr(T≧t)$ は確率密度関数 $f(t)$ を用いて次のように表せるはずです。

$Pr(T≧t)=\int_{t}^\infty f(x) dx$

これがさっき説明したところの $Pr(X≧i)$ です。

そうするとハザード率 $hazard(t)$ は
$hazard(t)=\frac{Pr(T=t)}{Pr(T≧t)}=\frac{f(t)}{\int_{t}^\infty f(x) dx}$

これもさっきの理屈と同じです。解釈もさっきの理屈と同じです。
一点違うのは時間が離散から連続になっているという点ですね。

これは見方を変えると次のように見ることも出来ますね。

$dt→0$ のもとで
$hazard(t)dt = Pr(t≦T≦t+dt|T≧t)$

そのため、よくハザード率の説明では次のような定義が成されます。

$hazard(t) =\displaystyle \lim_{dt \to 0} \frac{Pr(t≦T≦t+dt|T≧t)}{dt}$

これがよく見かける連続時間におけるハザード率の定義です。

ハザード率の仮定及び推定

離散時間、連続時間の2種類の観点からハザード率について説明しました。
しかし、話が一般的すぎます。
離散時間で考えると $Pr(X=i | X≧i)$ がハザード率だ！
連続時間で考えると $\displaystyle \lim_{dt \to 0} \frac{Pr(t≦T≦t+dt|T≧t)}{dt}$ がハザード率だ！

と言われても、だからどうした、つまりそれをどう扱えばよいのだとなりかねません。
そこで確率分布を何らかの方法で特定しましょう。

ノンパラメトリック推定

離散時間でものを考えるのならば、話は中学生の時の確率レベルで構いません。
要は $i$ 歳まで生きていた人の人数 $n_i$ と、i歳時点で亡くなってしまった人の人数 $D_i$ が分かれば良いわけです。
よって、ハザード率の推定量を次のように置くことにしましょう。

$\hat{hazard(i)}=\frac{D_i}{n_i}$

このように単純な相対度数を取った推定方法を経験分布による推定と呼びます。
確率論の定義に則れば、サンプルサイズ無限大で相対度数は真の分布と定義上等しいので、真っ当な推定方法です。
これをもう少し大きな枠組みでノンパラメトリックな推定と呼んだりします。

ただし、一般的にノンパラメトリックな推定はパラメトリックな推定と比べるとefficientではないという欠点があります。

パラメトリック推定

ノンパラメトリック推定は分布の仮定をかなり広く取った(正確にはパラメータ空間を広く取った)推定方法でしたが、もう少し分布を限定してしまうことも出来ます。

例えば連続時間で考えた上で、分布 $f(t)=\frac{1}{c}exp(-\frac{t}{c})$ を考えます。
これを片側指数分布と呼びます。

この場合 $\int_{t}^\infty f(x) dx=exp(-\frac{t}{c})$ のため、
$hazard(t) = \frac{1}{c}$
になります。

つまり、時間に依存せず、ハザード率は等しいというような分布の仮定になります。
これは生存解析の文脈でハザード率の説明をする際には具体例として必ず取り上げられる特徴的な仮定です。

打ち切りデータにおけるハザード率と生存確率

ただ、現実はそうそう甘くない。先ほどの保険屋さんのデータ例は非常に綺麗なデータでしたが、Policy holder が途中で離脱してしまうということも考えられます。つまり契約の解約です。この場合、少なくとも〇〇歳まで生きていたことはわかるが、それ以上はよくわかりませんというデータになるわけです。

カプランマイヤー推定量

具体例

それでは、打ち切りデータの具体例を見てみます。
RのパッケージsurvBootOutliersにはget.whas100.dataset()という関数でwhas100という打ち切りデータを取ってくることが出来るので、それを使います。

> library(survival)
> library(survBootOutliers)
> Data=get.whas100.dataset()
> head(Data)
  X los age gender      bmi times status
1 1   4  65      0 31.38134     6      1
2 2   5  88      1 22.65790   374      1
3 3   5  77      0 27.87892  2421      1
4 4   9  81      1 21.47878    98      1
5 5   4  78      0 30.70601  1205      1
6 6   7  82      1 26.45294  2065      1

病気の患者のデータですね。
statusが0の場合は生きているか、観察打ち切り(退院、転院等)のどちらかの打ち切りデータです。
statusが1の場合はもう亡くなってしまっていると考えます。
timesは入院から(死亡or離脱or生存)の期間を表します。

survivalパッケージのsurvを利用すると、timesとstatusの関係を整理することが出来ます。

> Surv(Data$times,Data$status)
  [1]    6   374  2421    98  1205  2065  1002  2201   189  2719+ 2638+  492 
 [13]  302  2574+ 2610+ 2641+ 1669  2624  2578+ 2595+  123  2613+  774  2012 
 [25] 2573+ 1874  2631+ 1907   538   104     6  1401  2710   841   148  2137+
 [37] 2190+ 2173+  461  2114+ 2157+ 2054+ 2124+ 2137+ 2031  2003+ 2074+  274 
 [49] 1984+ 1993+ 1939+ 1172    89   128  1939+   14  1011  1497  1929+ 2084+
 [61]  107   451  2183+ 1876+  936   363  1048  1889+ 2072+ 1879+ 1870+ 1859+
 [73] 2052+ 1846+ 2061+ 1912+ 1836+  114  1557  1278  1836+ 1916+ 1934+ 1923+
 [85]   44  1922+  274  1860+ 1806  2145+  182  2013+ 2174+ 1624   187  1883+
 [97] 1577    62  1969+ 1054

数字はtimesです。そして、status=0の場合は数字の横に+がついています。

このような打ち切りデータの場合、先ほどのようなハザード率・生存確率の推定では、打ち切りを考慮できていません。

カプランマイヤー推定量

先ほどのデータを時間順に並び替えてみることにします

> S=Surv(Data$times,Data$status)
> S[order(S)]
  [1]    6     6    14    44    62    89    98   104   107   114   123   128 
 [13]  148   182   187   189   274   274   302   363   374   451   461   492 
 [25]  538   774   841   936  1002  1011  1048  1054  1172  1205  1278  1401 
 [37] 1497  1557  1577  1624  1669  1806  1836+ 1836+ 1846+ 1859+ 1860+ 1870+
 [49] 1874  1876+ 1879+ 1883+ 1889+ 1907  1912+ 1916+ 1922+ 1923+ 1929+ 1934+
 [61] 1939+ 1939+ 1969+ 1984+ 1993+ 2003+ 2012  2013+ 2031  2052+ 2054+ 2061+
 [73] 2065  2072+ 2074+ 2084+ 2114+ 2124+ 2137+ 2137+ 2145+ 2157+ 2173+ 2174+
 [85] 2183+ 2190+ 2201  2421  2573+ 2574+ 2578+ 2595+ 2610+ 2613+ 2624  2631+
 [97] 2638+ 2641+ 2710  2719+

ちなみにこのデータはデータ数100です。
timesに重複があるのがちょっとややこしいですが、一旦説明の都合上完全に順序があると考えてください。
つまり、6日で2人亡くなっているわけですが、時間軸上先に一人目が亡くなって、次に二人目が亡くなったと考えてください。
つまり
$t_{1}＜t_{2}＜,\cdots,＜t_{100}$

だと考えているものとします。
更に、 $t_{i}$ 時の切断の有無を $d_{i}$ と表し、切断が成された場合(+がついている場合)は $d_{i}=0$ とします。

この時、 $t_{j}$ 時点での生存確率はノンパラメトリックに(経験分布的に相対度数で)考えると次のようになりますね。

$\prod_{k≦j} (\frac{n-k}{n-k+1})^{d{k}}$

この生存確率の推定量をカプランマイヤー推定量と呼びます。

生存曲線

このカプランマイヤー推定量は次のようにplotして確認してみるのが一般的です。
ちなみに下のプロットを生存曲線と呼びます。

plot(survfit(Surv(Data$times,Data$status)~Data$gender),lty=c(1,2))
legend(locator(1),c("male","female"), lty=c(1,2))

f:id:bananarian:20190214024928p:plain
ちなみに、普通生存曲線は比較する二つの群をプロットして、どちらの生存確率が高いかといったことを確認するものなので、survfit関数では比較する群を識別できるようformulaで記述します。
今回は男女で分けてみました。理由は特にありません。

信頼区間

カプランマイヤー推定量の信頼区間もsurvfitでは求めることができます。

plot(survfit(Surv(Data$times,Data$status)~Data$gender),lty=c(1,2))
lines(survfit(Surv(Data$times,Data$status)~Data$gender,conf.type="log"),mark.T=F,conf.int=TRUE,col=c("black","red"))
legend(locator(1),c("male","female"),lty=1,col=c("black","red"))

f:id:bananarian:20190214121046p:plain

ただし、この信頼区間の推定、恐らくグリーンウッドの公式を用いた方法が取られていると思われますが、この公式、生存確率が0か1に近い場合過少推定が起こると指摘されています。
以下参照
https://www.jstage.jst.go.jp/article/dds/30/5/30_474/_pdf

パラメトリックアプローチ

カプランマイヤーはノンパラメトリックな推定方法でしたが、パラメトリックに考えることも当然可能です。

Rではsurvreg関数が用意されており、distに分布を指定すればそれに合わせてパラメトリックに生存解析を行うことが可能です。

> survreg(Surv(Data$times,Data$status)~Data$gender,dist="logistic")
Call:
survreg(formula = Surv(Data$times, Data$status) ~ Data$gender, 
    dist = "logistic")

Coefficients:
(Intercept) Data$gender 
  2293.8562   -676.2864 

Scale= 871.8654 

Loglik(model)= -478.8   Loglik(intercept only)= -480.8
        Chisq= 3.94 on 1 degrees of freedom, p= 0.0472 
n= 100

2群の検定(ログランク検定)

2群のハザード率・生存確率を比較するにあたり、２群のハザード率が等しいか否かを検定する方法にログランク検定があります。ログランク検定の帰無仮説は次のようになります。

帰無仮説: 　全ての時間において2群のハザード率が等しい or (生存確率が等しい)

Rではsurvdiff関数によって行えます。

> survdiff(Surv(Data$times,Data$status)~Data$gender)
Call:
survdiff(formula = Surv(Data$times, Data$status) ~ Data$gender)

               N Observed Expected (O-E)^2/E (O-E)^2/V
Data$gender=0 65       28     34.6      1.27      3.97
Data$gender=1 35       23     16.4      2.68      3.97

 Chisq= 4  on 1 degrees of freedom, p= 0.05

結局ハザード率って何なの

ここで一度整理しておきます。ハザード率の確率論的解釈についてはこの記事の初めでお話しましたが、結局のところノンパラメトリックな観点から見たハザード率は、確率密度の代わりであると考えることも可能です。

ちなみにパラメトリックな解析においては、パラメータの推定値を利用して逆算的にハザード率を求めることも可能です。

cox 比例ハザードモデル

最後に、２群じゃなくてもっと多くの条件からハザードを確認し、多くの共変量の中から特に影響の大きいものを探したいといったモチベーションがあった場合、今までの手法ではうまくいきません。そこでcox比例ハザードモデルというセミパラメトリックな方法を取ります。

比例ハザード

まず、比例ハザードについて説明します。
ある個人iの情報(共変量)として、 $X_i$ が与えられているとします。
この時、 $i$ の $t$ 時点でのハザード率は次のようになると仮定します。

$hazard_i(t) = hazard_0(t) G(X_i)$

$G(X_i)$ は $X_i$ の関数です。
$hazard_0(t)$ はベースラインハザードと呼ばれます。共変量に全く依存していないことに注意してください。
要は、ハザードは比較実験のための情報が同じであれば一定であるという仮定をおいているわけです。
医療実験の枠組みで解釈するならベースラインハザードはプラセボ群ということになりますね。

更に共変量に添え字 $t$ が無いことにも注意です。このモデルでは、時間によって共変量の影響が変化しないという仮定も入っています。

cox比例ハザードモデル

そして、 $G(X_i)=exp(X_i \beta)$ と置くことで、これをcox比例ハザードモデルと呼びます。

このモデルの尤度を解くことは難しいので、部分尤度から近似推定量を求めます。
これを部分尤度最大化問題と呼びます。

で、この部分尤度最大化問題まで落とし込めれば基本的には解析的に求めることが出来るわけですが、近似解を求める手法も提案されています。

coxph

survivalパッケージではcoxph関数でcox比例ハザードモデルを推定することが出来ます。
methodに対してexactを与えると、部分尤度最大化問題を解析的に求めます。
efronやbreslowを与えると近似解を求めます。近似解はあくまで近似なので、計算量的な問題を気にしないのであればexactの方が良いと思います。デフォルトはefronです。

> summary(coxph(Surv(Data$times,Data$status)~Data$gender+Data$bmi+Data$age, method="exact"))
Call:
coxph(formula = Surv(Data$times, Data$status) ~ Data$gender + 
    Data$bmi + Data$age, method = "exact")

  n= 100, number of events= 51 

                coef exp(coef) se(coef)      z Pr(>|z|)   
Data$gender  0.14484   1.15585  0.30622  0.473  0.63622   
Data$bmi    -0.07093   0.93153  0.03610 -1.965  0.04945 * 
Data$age     0.03713   1.03783  0.01273  2.917  0.00354 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

            exp(coef) exp(-coef) lower .95 upper .95
Data$gender    1.1559     0.8652    0.6342    2.1065
Data$bmi       0.9315     1.0735    0.8679    0.9998
Data$age       1.0378     0.9636    1.0123    1.0640

Concordance= 0.683  (se = 0.037 )
Rsquare= 0.194   (max possible= 0.984 )
Likelihood ratio test= 21.55  on 3 df,   p=8e-05
Wald test            = 19.47  on 3 df,   p=2e-04
Score (logrank) test = 20.83  on 3 df,   p=1e-04

bmiや年齢が影響を及ぼしていることが確認出来ました。

検定によるチェック

で、cox比例ハザードモデル、かなり仮定が強い印象があります。
この仮定が妥当かどうかを検定でチェックしておく必要があります。
survivalではcox.zphによって行えます。
この検定については私は何をやっているのか知らないのですが(リサーチ不足)
とりあえず先人の知恵を盲信して試しておきます。

> cox.zph(coxph(Surv(Data$times,Data$status)~Data$gender+Data$bmi+Data$age))
                rho  chisq     p
Data$gender -0.0325 0.0671 0.796
Data$bmi     0.1420 1.5793 0.209
Data$age     0.0536 0.1986 0.656
GLOBAL           NA 1.6627 0.645