생존분석 - 모수적 모형 (Weibull model)

생존분석 - 모수적 모형

Cox의 비례위험 모형

일반적으로 생존분석에 사용되는 모델인 Cox의 비례위험모형의 경우 준모수적 (Semi-parametric) 방법입니다. 이 방법은 생존함수 (Survival function)을 추정하기 위해 Hazard function에 대한 모델링을 합니다. 시간에 따라 변하는 Baseline hazard가 존재하고, Covariate 들의 effect에 따라 hazard가 달라진다고 모델링합니다. 그리고 hazard의 effect인 Hazard ratio가 시간에 따라 동일하다는 비례위험을 가정하여 Hazard function과 Survival function 을 추정하게 됩니다. 이 방법의 장점은 Survival function에 분포에 대해서 가정을하지 않기 때문에 강건하며, Covariate를 보정할 수 있다는 것입니다. 또한 Baseline hazard function이 어떤 함수인지를 알 필요가 없이 계수 추정을 통해 Hazard ratio만 구할 수 있다면 Covariate의 effect를 파악할 수 있습니다. 

모수적 모형

반면 모수적 모형은 Hazard function과 Survival function의 모양에 대해 Cox 보다 더 강한 '가정'을 합니다. 이 가정은 Hazard function과 Survival function이 어떤 분포를 따를 것이라는 믿음에서 출발합니다. 가장 많이 사용되는 모수적 모델은 Weibull 모형입니다. 

Weibull 모형의 Hazard function과 Survival function 은 다음과 같이 정의됩니다.

$$ h(t) = \lambda p t^{p-1} $$

$$ S(t) = exp(-\lambda t^p) $$ 

h(t)와 S(t)에 정의에 따라 h(t)가 저렇게 정의되면 S(t)가 저렇게 계산되는 것을 확인할 수 있습니다.

$$ \lambda = exp(\beta_0 + \beta_1 TRT) $$ 

만약 Treatment 변수에 대해 Weibull model을 만들 경우, 이 때 위와 같이 식을 놓고, beta0, beta1, p 3개의 계수를 추정하게 됩니다. 

또는 Hazard를 기준으로 계수를 추정하는 것이 아니라 Survival time을 기준으로 추정할 수 있는데, 이 경우가 모수적 방법에서 많이 사용되는 방법입니다. 이를 Acceleration Failure Time Model이라고 합니다. 

$$ S(t) = exp(-\lambda t^p) $$

$$ t = [-lnS(t)]^{1/p} \times \frac{1}{\lambda^{\frac{1}{p}}} $$

$$ \frac{1}{\lambda^{\frac{1}{p}}} = exp(\alpha_0 + \alpha_1 TRT) $$

AFT 모델에서는 위와 같이 놓고 alpha0, alpha1, p를 추정하게 됩니다. 

이렇게 두 가지 방법으로 구한 계수 beta와 alpha의 의미는 다른데, 만약 p가 고정되있다고 한다면, beta는 exponential을 취하면 Hazard ratio를 뜻하며, alpha는 exponential을 취하면 Accerleration factor가 됩니다. Acceleration factor에 대해서는 본 포스팅에서는 다루지 않을 예정이므로 이 자료를 참고 바랍니다. 

Strata 별로 p가 동일하다고 가정한 경우, Weibull model은 Proportional Hazard 가정과 Acceleration Failure Time 가정을 만족합니다. (https://en.wikipedia.org/wiki/Accelerated_failure_time_model)

PBC 데이터를 통한 Weibull model의 구축

pbc.dat

library(survival)

data <- read.table("pbc.dat")

# Column rename
colnames(data) <- c('id', 'futime', 'status', 'drug', 'age', 'sex', 'ascites', 'hepato', 'spiders', 'edema', 'bili', 'chol', 'albumin', 'copper', 'alk_phos', 'sgot', 'trig', 
                   'platelet', 'protime', 'stage')

head(data)

data <- data[data$drug != '.', ]
data$drug <- as.character(data$drug)
data$drug[(data$drug == 2)] = 0
data$drug <- factor(data$drug)

Drug = 0 플라시보, Drug = 1 페니실린에 대해 Death에 대한 Time to event 분석을 수행합니다.

Weibull Model Fitting

library(survival)
library(eha)

# AFT  
swei = survreg(s~drug, dist = "weibull",data = data) 

# PH  
sweiph = phreg(s~drug, dist = "weibull",data = data)

summary(swei) 
summary(sweiph)

Call:
survreg(formula = s ~ drug, data = data, dist = "weibull")
              Value Std. Error     z      p
(Intercept)  8.3167     0.1094 75.99 <2e-16
drug1       -0.0438     0.1431 -0.31  0.759
Log(scale)  -0.1532     0.0729 -2.10  0.036

Scale= 0.858 

Weibull distribution
Loglik(model)= -1348.2   Loglik(intercept only)= -1348.3
	Chisq= 0.09 on 1 degrees of freedom, p= 0.76 
Number of Newton-Raphson Iterations: 5 
n= 312 

Call:
phreg(formula = s ~ drug, data = data, dist = "weibull")

Covariate          W.mean      Coef Exp(Coef)  se(Coef)    Wald p
drug 
               0    0.491     0         1           (reference)
               1    0.509     0.051     1.052     0.167     0.759 

log(scale)                    8.317               0.109     0.000 
log(shape)                    0.153               0.073     0.036 

Events                    144 
Total time at risk        625985 
Max. log. likelihood      -1348.2 
LR test statistic         0.09 
Degrees of freedom        1 
Overall p-value           0.759387

Median Survival Time Comparison

predict(swei,newdata=data.frame(drug='0'),type=c("quantile"),se.fit=FALSE,p=0.5)
predict(swei,newdata=data.frame(drug='1'),type=c("quantile"),se.fit=FALSE,p=0.5)

1: 2987.74956079141

1: 2859.63735172644

Weibull Model을 통해 구한 Survival Function과 Kaplen-Meier estimation의 비교

data$event_all=ifelse(data$status==0,0,1)
data$futime = ifelse(data$futime==0,0.1,data$futime)
s = Surv(data$futime,data$event_all)

sKM <- survfit(s ~ drug, data=data)
plot(sKM)
lines(predict(swei,newdata=list(drug='0'),type="quantile", p=seq(.01,.99,by=.01)),seq(.99,.01,by=-.01),col='red')
lines(predict(swei,newdata=list(drug='1'),type="quantile", p=seq(.01,.99,by=.01)),seq(.99,.01,by=-.01),col='red')

위 그림에서 검은색 직선이 Kaplen-Meier Survival function estimation이며, 붉은색과 파란색이 Weibull 모델로 추정한 Survival function 입니다. Weibull 모델의 S(t)의 모양에 대한 가정이 맞다면 두 직선은 가까울 것입니다. 그리고 위 직선들은 p가 같다는 가정하에 구해진 직선들입니다. 즉, PH, AFT 가정을 한 후에 그려진 직선입니다. 그렇기 때문에 Weibull 모형의 Survival function에 대한 가정, 그리고, PH, AFT 가정을 해도 될지에 대해서 체크를 해보아야합니다.

Weibull 가정의 체크 

(https://stat.ethz.ch/education/semesters/ss2011/seminar/contents/handout_9.pdf)

Weibull 모형에서 가정을 체크하는 방법은 X 축을 log(t), Y 축을 log(-log(S(t)) 로 그래프를 그려보는 것입니다. 이 때, S(t)는 Kaplen-Meier 방법으로 구한 Survival function 입니다. 이 때, 그룹간의 그래프의 모양을 보고 Weibull 생존모형 가정이 맞는지, 또한 두 그룹에서 PH, AFT 가정이 맞는지를 체크할 수 있습니다. 

먼저, 두 그룹의 직선이 평행하게 나오면 Weibull 가정을 만족합니다. 그렇지 않으면 Weibull 가정이 만족되지 못하며 AFT 가정도 만족되지 못합니다. 또한 두 그룹의 기울기의 비가 일정하면 PH 가정을 만족합니다. 

sKM <- survfit(s~ drug, data=data)
kmfit.TRT1 = summary(sKM)
kmfit.TRT2 = data.frame(kmfit.TRT1$strata, kmfit.TRT1$time, kmfit.TRT1$surv)
names(kmfit.TRT2) = c("TRT","time","survival")
TRT0=kmfit.TRT2[kmfit.TRT2$TRT=="drug=0", ]
TRT1=kmfit.TRT2[kmfit.TRT2$TRT=="drug=1", ]
plot(log(TRT0$time), log(-log(TRT0$survival)), pch=16, xlab="time in days using log-arithmic scale", ylab="log-log survival",type="p", xlim=c(3,10), ylim=c(-6,1), font=2, lwd=2, font.lab=2, col="Black",main="log-log curves by treatment")
par(new=T)
plot(log(TRT1$time), log(-log(TRT1$survival)), pch=3, col="grey50", ylab="",xlab="",  xlim=c(3,10), ylim=c(-6,1), font=2, lwd=2, font.lab=2)

이 그래프의 경우 우선 그룹에서 거의 직선이라는 것을 확인할 수 있지만, 두 그래프가 평행하지는 않습니다. (교차가 일어납니다.) 그렇기 때문에 Weibull 가정을 만족하지만 PH 가정과 AFT 가정을 만족하지 않으며, 두 그룹에서 p가 다르다는 알 수 있습니다.