12  비례위험모형

12.1 필요한 패키지

library(tidyverse)
library(kableExtra)
library(ggplot2)
library(here)

library(survival)
library(KMsurv)
library(asaur)
library(survminer)

12.2 비례위험모형의 개요

설명변수 또는 공변량 \(x\) 의 값에 따라서 반응변수 \(y\)의 평균이 어떻게 다른지 나타내는 모형을 회귀모형(regression model) 이라고 한다.

일반화 선형모형(generalized linear model)은 반응변수의 평균 \(E(y|x)\)과 선형예측식 \(\eta =\beta_0+\beta_1 x\) 의 관계를 다음과 같이 연결함수 \(g\) 를 사용하여 모형을 만든다.

\[ g[E(y|x)]= \eta = \beta_0+\beta_1 x \tag{12.1}\]

보기 12.1 (정규분포 - 직선 회귀모형) 반응변수 \(y\) 가 정규분포를 따르는 경우 일반화 선형모형 12.1 의 연결함수는 항등함수 \(g(x) = x\) 이며 일반적인 직선 회귀모형으로 나타난다.

\[ E(y|x) = \beta_0+\beta_1 x \]

\(\blacksquare\)

보기 12.2 (이항분포 - 로지스틱 회귀모형) 반응변수 \(y\) 가 이항분포를 따르는 경우 일반화 선형모형 12.1 의 연결함수는 로짓함수이며 로지스틱 회귀모형으로 나타난다.

\[ \log \left( \frac{E(y|x)}{1-E(y|x)} \right) = \beta_0+\beta_1 x \] 또는

\[ E(y|x) = P(y=1|x) = \frac{\exp(\beta_0+\beta_1 x)}{1+\exp(\beta_0+\beta_1 x)} \] \(\blacksquare\)

생존시간을 다루는 생존분석에서도 공변량이 생존시간 \(T\) 에 영향을 미치는 모형을 위와 같은 일반화 선형모형으로 고려하는 경우도 있다. 하지만 생존분석의 회귀모형에서는 생존시간의 평균보다 위험함수를 더 중요한 요소로 보기 때문에 공변량과 위험함수의 관계에 설정하는 회귀모형을 더 많이 사용한다.

생존분석에서 회귀모형을 고려하는 경우 주의할 점은 위험함수 \(h(t)\) 는 하나의 모수가 아닌 시간 \(t\)의 함수라는 것이다.

생존분석에는 공변량의 값에 따라서 위험함수가 변하게 되며 공변량의 값이 0인 경우 나타나는 기저 위험함수(baseline hazard function) 은 평균과 같이 하나의 값이 아닌 시간의 함수로 주어진다는 것을 유의하자.

이제 공변량 \(x\)과 위험함수 \(h(t)\)의 관계를 다음과 같은 비례위험 모형으로 표현할 수 있다. \[ h(t) = h_0(t) \gamma(x) \tag{12.2}\]

식 12.2 를 비례위험모형(proportional hazard model) 이라고 부르며 \(h_0(t)\) 는 기저 위험함수(baseline hazrard function)라고 부른다.

일반적으로 식 12.2 에서 함수 \(\gamma(x)\)를 선택하는 경우 공변량 \(x=0\) 이면 함수의 값이 1이 되도록 설정한다 (\(\gamma(0)=1\)). \(\gamma(x)=1\) 인 경우, 즉 공변량이 0인 경우의 위험함수가 기저 위험함수 \(h_0(t)\) 가 된다.

식 12.2 을 비례위험모형이라고 부르는 공변량이 위험함수에 미치는 영향이 가법적인 아닌 비례적으로 나타나며 시간에 따라 변하지 않는다는 의미이다. 예를 들어 두 개의 서로 다음 공변량 \(x_1\) 과 \(x_2\) 를 고려하면 모형 12.2 에 의하여 두 위험함수의 비는 다음과 같이 나타난다.

\[ \frac{ h(t|x_2)}{h(t|x_1)} = \frac{ h_0(t) \gamma(x_1) }{h_0(t) \gamma(x_2)} = \frac{ \gamma(x_1) }{\gamma(x_2)} \tag{12.3}\]

위의 모형 12.3 에서 유의할 점은 시간에 따라서 위험함수의 상대적 비율 \(\gamma(x_1)/\gamma(x_2)\) 은 변하지 않는다. 즉, 공변량의 영향이 모든 시간에서 모두 동일하다는 것이다.

만약에 공변량의 영향이 시간에 따라서 변하면 함수 \(\gamma\)도 시간에 따라서 변화하는 함수 \(\gamma(x,t)\)를 고려해야 한다. 시간에 따라서 공변량의 형태나 영향이 바뀌는 것을 시간의존 공변량(time-varying covariate) 라고 한다.

식 12.2 에 나타난 함수 \(\gamma(x)\) 를 다음과 같이 회귀 계수 \(\beta\) 포함된 지수함수의 형태로 나타낸 모형을 Cox 의 비례위험모형이라고 부른다.

\[ h(t) = h_0(t) \exp (\beta x) \tag{12.4}\]

유의할 점은 Cox 의 모형에서는 일반적인 회귀모형에서 사용하는 절편을 포함시키지 않는다는 것이다. 즉, 일반적인 회귀모형에서 사용되는 절편의 역활을 Cox 모형에서는 기저 위험함수 \(h_0(t)\) 가 하는 것이다. 따라서 Cox 모형에서는 공변량의 값이 0인 경우의 위험함수가 기저 위험함수 \(h_0(t)\) 가 된다.

Cox 비례위험모형에서는 위험함수의 비가 다음과 같이 나타난다. 아래에서 볼 수 있듯이 Cox 비례위험모형에서는 공변량이 1 단위 증가하면 위험함수가 \(\exp(\beta)\) 에 비례하게 증가한다.

\[ \frac{ h(t|x_1 +1)}{h(t|x_1)} = \frac{\exp(\beta[ x_1+1])}{\exp(\beta(x_1)} =\exp[\beta (x_1+1-x_1)] = \exp(\beta) \tag{12.5}\]

만약 공변량이 하나가 아니라 \(p\) 개를 가지는 경우 다음과 같이 모형을 확장할 수 있다.

\[ h(t) = h_0(t) \exp (\beta_1 x_1 + \beta_2 x_2 + \cdots + \beta_p x_p) \]

Cox 비례위험모형 12.2 에서 양변에 로그를 하면 다음과 같은 모형이 된다.

\[ \log h(t) = \alpha(t) + \beta x \tag{12.6}\]

모형 12.6 에서 \(\alpha(t) = \log h_0(t)\) 는 시간에 따라 변하는 기저 위험함수를 나타낸다.

모형 12.6 에 \(\alpha(t) =\alpha_0\) 로 상수이면 생존시간이 지수분포를 떠른다고 가정한 모형이 되고 \(\alpha(t) = \alpha_0 + \alpha_1 \log(t)\) 이면 생존시간의 분포가 와이블(weibull) 분포를 가정한 모형이 된다.

또한 생존시간과 위험함수의 관계 \(S(t) = \exp(-H(t))\)를 이용하면 Cox 비례위험모형에서는 다음과 같이 나타낼 수 있다.

\[ S(t) = S_0(t)^{\exp(\beta x)} \]

위의 식에서 \(S_0(t)\) 는 기저 위험함수 \(h_0(t)\) 를 가지는 생존함수이다.

12.3 부분 가능도 함수

Cox 비례위험모형은 특정한 분포를 가정하지 않고 공변량과 위험함수의 관계만을 설정하고 있다. 또한 기저 위험함수의 특정한 형태도 가정하지 않는다. 따라서 Cox 비례위험모형의 가정 하에서는 완전한 가능도 함수(full likelihood function)를 구할 수 없다. 이러한 경우 부분 가능도함수(partial likelihood function)를 이용하면 특정한 분포의 가정없이 중도절단이 포함된 생존자료에서 회귀계수를 추정할 수 있다.

이제 독립적으로 관측된 \(n\) 개의 생존 시간 \(t_i\) 와 중도절단을 나타내는 함수 \(\delta_i\), 그리고 공변량 \(x_i\)가 주어졌다고 하자.

\[ (t_1, \delta_1, x_1), (t_2, \delta_2, x_2), \cdots , (t_n, \delta_n, x_n) \] 또한 위험함수는 Cox 비례위험모형 12.4 을 따른다고 가정하자. 즉 \(i\) 번째 관측값에 대한 위험함수 \(h_i(t)\)를 다음과 같이 나타낼 수 있다.

\[ h_i(t) = h_0 (t) \exp(\beta x_i) \tag{12.7}\]

만약 실제 생존 시간이 관측된 경우 \((\delta_i=1)\) 관측값 \(t_i\) 의 가능도 함수 \(L_i(\beta)\) 는 다음과 같다.

\[ L_i(\beta) = P(T=t_i) = f(t_i) = h_i(t_i) S_i(t_i) \tag{12.8}\]

위의 식에서 마지막 결과는 위험함수의 정의 \(h(t) = f(t)/S(t)\)를 이용한 것이다.

만약 중도절단이 일어난 경우 \((\delta_i=0)\) 관측값 \(t_i\) 의 가능도 함수 \(L_i(\beta)\) 는 다음과 같다.

\[ L_i(\beta) = P(T > t_i) = S_i(t_i) \tag{12.9}\]

이제 \(R(t)\) 를 시간 \(t\) 에서의 위험집단(risk set), 즉 시간 \(t\) 까지 살아있는 개체들의 집합이라고 하자.

식 12.8 과 12.9 를 합쳐서 모든 자료에 대한 가능도 함수 \(L(\beta)\) 를 다음과 같이 나타낼 수 있다.

\[ \begin{aligned} L(\beta) & = \prod_{i=1}^n L_i(\beta) \\ & = \prod_{i=1}^n h_i (t)^{\delta_i} S_i (t_i) \\ & = \prod_{i=1}^n \left [ \frac{h_i (t_i)} { \sum_{j \in R(t_i)} h_j (t_i)} \right ]^{\delta_i} \left [ \sum_{j \in R(t_i)} h_j (t_i) \right ]^{\delta_i} S(t_i) \\ & = \underbrace{ \prod_{i=1}^n \left [ \frac{h_i (t_i)} { \sum_{j \in R(t_i)} h_j (t_i)} \right ]^{\delta_i} }_{\text{부분 가능도함수}} \underbrace{ \prod_{i=1}^n \left [ \sum_{j \in R(t_i)} h_j (t_i) \right ]^{\delta_i} S(t_i) }_{\text{기저 위험함수를 포함}} \end{aligned} \tag{12.10}\]

이제 가능도 함수 \(L(\beta)\) 12.10 의 첫 번째 부분을 먼저 자세히 살펴보자. 아래에서는

\[ \begin{aligned} \prod_{i=1}^n \left [ \frac{h_i (t_i)} { \sum_{j \in R(t_i)} h_j (t_i)} \right ]^{\delta_i} & = \prod_{i=1}^n \left [ \frac{h_0(t_i) \exp(\beta x_i)} { \sum_{j \in R(t_i)} h_0(t_i) \exp(\beta x_j)} \right ]^{\delta_i} \\ & = \prod_{i=1}^n \left [ \frac{\exp(\beta x_i)} { \sum_{j \in R(t_i)} \exp(\beta x_j)} \right ]^{\delta_i} \end{aligned} \tag{12.11}\]

가능도 함수 \(L(\beta)\)의 앞 부분은 위의 식 12.11 에서 나타났듯이 기저 위험함수 \(h_0(t)\) 가 나타나지 않고 공변량과 회귀계수 \(\beta\) 만의 함수로 구성되어 있다. 또한 중도절단이 일어난 관측값(\(\delta_i=0\))은 식 12.11 의 계산에서 제외된다. 식 12.11 에 주어진 함수를 부분 가능도함수라고 부른다.

가능도 함수 \(L(\beta)\) 12.10 의 두 번째 부분은 기저 위험함수 \(h_0(t)\) 가 포함된 부분으로서 부분 가능도함수를 이용하여 회귀계수를 추정하는 경우는 이 부분을 제외하고 추정하는 것이다.

Cox 비례위험모형에서 최대 가능도 방법은 아래 주어진 부분 가능도함수 \(L_p(\beta)\) 를 최대화하는 회귀계수를 구하는 것이다.

\[ \begin{aligned} L_p(\beta) & = \prod_{i=1}^n \left [ \frac{h_i (t_i)} { \sum_{j \in R(t_i)} h_j (t_i)} \right ]^{\delta_i} \\ & = \prod_{i=1}^n \left [ \frac{ h_0(t) \exp(\beta x_i)} { \sum_{j \in R(t_i)} h_0(t) \exp(\beta x_j)} \right ]^{\delta_i} \\ & = \prod_{i=1}^n \left [ \frac{ \exp(\beta x_i)} { \sum_{j \in R(t_i)} \exp(\beta x_j)} \right ]^{\delta_i} \end{aligned} \]

보기 12.3 이제 Moore (2016) 에 주어진 간단한 예제를 통하여 Cox 비례위험모형을 적합하는 방법을 살펴보자. 다음과 같은 자료가 주어졌다고 하자.

Patient	Survtime	Censor	Group	x_i
1	6	1	C	0
2	7	0	C	0
3	10	1	T	1
4	15	1	C	0
5	19	0	T	1
6	25	1	T	1

위의 자료에서 생존시간(Survtime)과 중도절단여부(Censor)를 이용하여 Cox 비례위험모형을 적합해보자. 설명변수는 다음과 같이 그룹 \(C\) 이면 \(x=0\) 이고 그룹 \(T\) 이면 \(x=1\) 이다.

\[ x_i = \begin{cases} 0 & \text{if group is C} \\ 1 & \text{if group is T} \end{cases} \tag{12.12}\]

먼저 첫 번째 생존 시간 \(t_1=6\) 에 대한 부분 가능도함수를 구해보자. 이 경우 위험그룹 \(R(6)\) 는 모든 환자이며 부분 가능도함수는 다음과 같다.

\[ L_1(\beta) = \frac{ \exp(\beta x_1)} { \sum_{j \in R(6)} \exp(\beta x_j) } = \frac{ \exp(\beta x_1)} { \sum_{j=1}^6 \exp(\beta x_j) } \]

두 번째 생존 시간 \(t_2=7\) 은 중도절단이므로 부분 가능도함수에 포함되지 않는다.

세 번째 생존 시간 \(t_3=10\) 에 대한 부분 가능도함수는 다음과 같다. 이때 위험그룹 \(R(10)\) 는 3, 4, 5, 6 번 환자이다.

\[ L_3(\beta) = \frac{ \exp(\beta x_3)} {\sum_{j \in R(10)} \exp(\beta x_j) } = \frac{ \exp(\beta x_3)} { \sum_{j=3}^6 \exp(\beta x_j) } \] 나머지 생존 시간에 대한 부분 가능도함수를 구하면 다음과 같다.

\[ L_4(\beta) = \frac{ \exp(\beta x_4)} { \sum_{j=4}^6 \exp(\beta x_j) }, \quad L_6(\beta) = \frac{ \exp(\beta x_6)} { \sum_{j=6}^6 \exp(\beta x_j) } = 1 \] 이제 위의 결과를 모두 곱하면 부분 가능도함수 \(L_p(\beta)\) 를 구할 수 있다. 아래 부분 가능도함수를 유도하는 경우 12.12 에서 정의된 설명변수를 사용하였다.

\[ \begin{aligned} L_p(\beta) & = L_1(\beta) L_3(\beta) L_4(\beta) L_6(\beta) \\ & = \left [ \frac{ \exp(\beta x_1)} { \sum_{j=1}^6 \exp(\beta x_j) } \right ] \left [ \frac{ \exp(\beta x_3)} { \sum_{j=3}^6 \exp(\beta x_j) } \right ] \left [ \frac{ \exp(\beta x_4)} { \sum_{j=4}^6 \exp(\beta x_j) } \right ] \left [ \frac{ \exp(\beta x_6)} { \sum_{j=6}^6 \exp(\beta x_j) } \right ] \\ & = \left [ \frac{ \exp(\beta x_1)} { \sum_{i=j}^6 \exp(\beta x_j) } \right ] \left [ \frac{ \exp(\beta x_3)} { \sum_{j=3}^6 \exp(\beta x_j) } \right ] \left [ \frac{ \exp(\beta x_4)} { \sum_{j=4}^6 \exp(\beta x_j) } \right ] (1) \\ & = \ \frac{ \exp(\beta x_1 + \beta x_3 + \beta x_4)} { [\sum_{j=1}^6 \exp(\beta x_j) ] [\sum_{j=3}^6 \exp(\beta x_j)] [\sum_{j=4}^6 \exp(\beta x_j)] } \\ & = \ \frac{ \exp(\beta) }{[3+3\exp(\beta)] [1+3\exp(\beta)] [1+2\exp(\beta)]} \end{aligned} \] 이제 위의 부분 가능도함수 \(L_p(\beta)\) 에 로그를 취한 로그 부분 가능도함수 \(\log L_p(\beta)\) 를 구하면 다음과 같다.

\[ \log L_p(\beta) = \beta - \log [3+3\exp(\beta)] - \log [1+3\exp(\beta)] - \log [1+2\exp(\beta)] \]

위의 로그 부분 가능도함수를 최대화하는 회귀계수 \(\beta\) 를 R 에서 제공하는 optim 함수를 사용하여 직접 구해보자.

logL <- function(beta) {
  psi <- exp(beta)
  result <- log(psi) - log(3*psi + 3) -
  log(3*psi + 1) - log(2*psi + 1)
  result }

result <- optim(par=0, fn = logL, method = "L-BFGS-B",  control=list(fnscale = -1),  lower = -3, upper = 1)

result$par

[1] -1.326129

위의 결과를 보먄 최대 가능도 추정량은 \(\hat{\beta} = -0.6931\) 이다. 회귀계수 \(\beta\) 의 의미는 두 그룹 간의 위험비를 나타낸다. 두 그룹 간의 위험비는 \(\exp(\beta)\) 이므로 그룹 \(T\) 의 위험은 그룹 \(C\) 의 26.5% 로 작게 나타난다.

\[ \frac{h(t | T)}{h(t |C)} = \exp(\hat{\beta}) = \exp(-1.326) = 0.265 \]

위에서 계산한 부분 가능도함수를 이용하여 Cox 비례위험모형 방법은 R 의 coxph() 함수를 사용하여 쉽게 구할 수 있다.

df <- data.frame(time = c(6, 7, 10, 15, 19, 25), 
                 censor = c(1, 0, 1, 1, 0, 1), 
                 treat = c("C", "C", "T", "C", "T", "T"))
df$treat <- factor(df$treat)
coxph(Surv(time, censor) ~ treat, data = df)

Call:
coxph(formula = Surv(time, censor) ~ treat, data = df)

          coef exp(coef) se(coef)     z     p
treatT -1.3261    0.2655   1.2509 -1.06 0.289

Likelihood ratio test=1.21  on 1 df, p=0.2715
n= 6, number of events= 4 

\(\blacksquare\)

보기 12.4 교과서 Jaewon Lee (2005) 의 예제 7.5은 HSV_2 를 가진 48명의 환자들이 gD2(antigen glycoprotein)을 기초로 하는 새로운 백신연구의 연구 결과 자료이다. 생존시간은 처음 재발할 때 까지의 시간(time) 리고 연구 참여 전 발생한 에피소드의 횟수(x)를 공변량으로 하여 Cox 비례위험모형을 적합해 보자. 처리 그룹(treat)은 gD2 로 치료받은 그룹(gD2)과 위약으로 치료반은 제어군(pbo)이다.

자료는 example75.txt에 저장되어 있다.

df75 <- read.csv(here("data","example75.txt"), header = T, sep="")
df75$treat <- factor(df75$treat)
head(df75)

  treat  x time censor
1   gd2 12    8      1
2   gd2  6   44      1
3   gd2 11   35      1
4   gd2  9   52      0
5   gd2 10    9      1
6   gd2 13   13      1

함수 coxph() 는 Cox 비례위험모형을 적함하는 함수이다.

먼저 공변량을 포함하지 않고 처리효과의 차이 만을 고려한 Cox 비례위험모형을 적합해보자.

ex75logrank <- coxph(Surv(time, censor) ~ treat , data= df75)
summary(ex75logrank)

Call:
coxph(formula = Surv(time, censor) ~ treat, data = df75)

  n= 48, number of events= 31 

           coef exp(coef) se(coef)     z Pr(>|z|)  
treatpbo 0.7138    2.0417   0.3678 1.941   0.0523 .
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

         exp(coef) exp(-coef) lower .95 upper .95
treatpbo     2.042     0.4898     0.993     4.198

Concordance= 0.595  (se = 0.048 )
Likelihood ratio test= 3.78  on 1 df,   p=0.05
Wald test            = 3.77  on 1 df,   p=0.05
Score (logrank) test = 3.91  on 1 df,   p=0.05

위의 결과를 보면 제어군의 재발 위험이 gD2 를 이용한 치료군에 비해 약 2배 높다는 것을 알 수 있다. 하지만 회귀계수는 통계적으로 유의하지 않다(p-value = 0.0523).

\[ \frac{h(t | {\textbf pbo} )}{h(t | {\textbf gD2})} = \exp(\hat{\beta}) = \exp(0.7138) = 2.0417 \]

다음으로 공변량 \(x\) 을 포함하는 Cox 비례위험모형을 적합해보자.

ex75cox <- coxph(Surv(time, censor) ~ treat + x, data= df75)
summary(ex75cox)

Call:
coxph(formula = Surv(time, censor) ~ treat + x, data = df75)

  n= 48, number of events= 31 

            coef exp(coef) se(coef)     z Pr(>|z|)   
treatpbo 0.90157   2.46346  0.37476 2.406  0.01614 * 
x        0.17577   1.19217  0.06561 2.679  0.00738 **
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

         exp(coef) exp(-coef) lower .95 upper .95
treatpbo     2.463     0.4059     1.182     5.135
x            1.192     0.8388     1.048     1.356

Concordance= 0.645  (se = 0.063 )
Likelihood ratio test= 10.81  on 2 df,   p=0.004
Wald test            = 10.73  on 2 df,   p=0.005
Score (logrank) test = 11  on 2 df,   p=0.004

위의 결과를 보면 연구 참여 전 발생한 에피소드의 횟수 \(x\) 가 1 증가할 때마다 재발 위험이 약 19% 증가한다는 것을 알 수 있다. 회귀계수는 통계적으로 유의하다(p-value = 0.01614).

\[ \frac{h(t | x+1 )}{h(t | x )} = \exp(\hat{\beta}) = \exp(0.17577) = 1.19217 \]

또한 공변량 \(x\)을 포함함 경우, 제어군의 재발 위험이 gD2 를 이용한 치료군에 비해 약 2.5 배 높다는 것을 알 수 있다. 회귀계수는 통계적으로 유의하다(p-value = 0.00738 ).

\[ \frac{h(t | {\textbf pbo, x} )}{h(t | {\textbf gD2,x})} = \exp(\hat{\beta}) = \exp(0.90157) = 2.463 \]

\(\blacksquare\)

Jaewon Lee, Hanna Yu, Mira Park. 2005. 생명과학연구를 위한 통계적 방법. 1st ed. 자유아카데미.

Moore, Dirk F. 2016. Applied survival analysis using R. Springer.