생존분석 - Cox 비례위험 모형의 가정 체크

Cox 비례위험 모형의 가정 체크

개념 : Cox 비례위험 모형은 "비례위험 가정" 을한다. 

비례위험가정은 시간에 따라 어떤 변수에서의 Hazard ratio가 일정해야한다는 것이다.

비례위험가정이 어긋나면, Stratified cox regression 또는 time-dependant variable 정의를 통한 extended cox regression을 수행할 수 있다. 

비례위험 가정의 체크는 아래방식으로 할 수 있다. 

1. Log-log plot

2. Observed-expected plot

3. Goodness of fit test

본 포스팅에서는 R의 survival, survminer 패키지를 통한 실습을 통해 cox 비례 위험 모형의 가정을 체크하는 방법을 알아본다. 

addict.csv

library(data.table)
library(survival)
data <- read.csv("addict.csv")

head(data)

ID	clinic	event	time	prison	dose
1	1	1	428	0	50
2	1	1	275	1	55
3	1	1	262	0	55
4	1	1	183	0	30
5	1	1	259	1	65

먼저 데이터의 형태는 위와 같이 생겼다. 

관심변수 : clinic 에 대하여 Kaplan-Meier 플롯 그리기

# Treatment 별로 그룹지어서 KM Estimation
fit = survfit(y~clinic, data=data)
par(mai=c(1,1,1,1))
plot(fit,main="KM Curve", xlab="Time(Week)", ylab=expression(paste(hat(S),"(t)")), col=c("blue", "red"))

우리는 이것이 비례위험 가정을 만족하는지를 확인하고 싶다. 비례위험 가정을 만족하려면 모든 timepoints에서 Hazard ratio가 같아야한다. hazard란 f(t)/S(t) (참고 : https://3months.tistory.com/349?category=743476) 로 위 그림만 봐서는 가정을 만족하는지 한 눈에 알 수 없다.  

로그로그 플롯

개념 : 로그로그 플롯을 그리고 그룹간에 평행하게 나오면 비례위험가정을 만족한다. 

plot(fit, fun="cloglog", lty=1:2, col=c("Black", "Grey50"), lwd=2, font=2, ylim=c(-3,2), font.lab=2, main="log-log KM curves by Rx", ylab="log-log survival", xlab="time in weeks using logarithmic scale")
legend("bottomright",lty=1:2,legend=c("Group1","Group2"),bty="n", text.font=2, lwd=2, col=c("Black", "Grey50"))

관찰값-기댓갓 플롯

개념 : Kaplen-Meier 생존 함수와 Cox regression을 통해 구한 생존함수를 동시에 그려본다.

이 때, Expected plot은 Cox 비례위험 모형에서 나오는 생존함수로 관찰한 생존함수와 모양의 차이가 많이나면 비례위험 가정을 만족하지 않는 것이다. 

kmfit <- survfit(y~clinic, data=data)
#observed
plot(kmfit, lty="dashed", col=c("Black", "Grey50"), lwd=2, font=2, font.lab=2, main="Observed Versus Expected Plots by Clinic", ylab="Survival probability", xlab="Survival time")
par(new=T)
#expected
exp <- coxph(y~clinic, data=data)
new_df <- data.frame(clinic=c(1,2))
kmfit.exp <- survfit(exp, newdata=new_df)
plot(kmfit.exp, lty="solid", col=c("Black", "Grey50"), lwd=2, font=2, font.lab=2)

Goodness of fit test

개념 : Goodness of fit test(GoF Test) 는 관찰값-기댓값의 차이를 통계적으로 검정하는 방법이다. (카이스퀘어 분포 이용)

# Gooodness of fit test
mod1 <- coxph(y~clinic, data=data)
cox.zph(mod1, transform=rank)

plot(cox.zph(mod1, transform=rank), se=F, var="clinic")
abline(h=0, lty=3)
plot(cox.zph(mod1, transform=rank), se=F, var="prison")
abline(h=0, lty=3)

        rho chisq       p
clinic -0.255  10.5 0.00118

p-value가 0.001이므로 clinic 변수는 비례위험 가정을 만족하지 않는다. 이는 그림을 통해서도 확인할 수 있다. 위 그림에서 직선이 평행하지를 확인하면되는데, 직선은 추정된 beta값을 의미하는데 직선이 평행하지 않고 어떠한 경향성이 보이면 beta 값이 시간에 따라 바뀐다는 의미로 비례위험 가정에 어긋난다. 

참고 : https://www.r-bloggers.com/cox-model-assumptions/