생존분석 (R)

 

  생존분석은 의학연구에서 굉장히 흔히 쓰이는 분석기법이지만, 꽤 흔한 빈도로 잘 못 사용되는 분석법 이기도 하다. 언뜻 생존분석은 로지스틱 회귀분석과 비슷하다는 생각이 들 수도 있지만, 생존분석에는 시간이라는 개념이 들어가기 때문에, 시간에 따른 변화를 읽어내지 못한 채 분석을 하게 되면 잘못된 결과를 얻을 수 있다. 이번 포스트에서는 생존분석의 실질적인 측면에 대해 살펴보기로 하자.

 

 

1. Kaplan-Meier curve

 

  상당수의 의학논문에서 생존분석을 이용해서 분석을 하는 경우, Kaplan-Meier curve를 제시하고 그 후 cox 분석을 주요분석법으로 사용한다. Kaplan curve를 보여주는 이유는 아무래도 단조로운 논문에 그래도 그림이 하나쯤은 필요하기 때문이기도 하지만 더 큰 이유는 cox 분석을 사용하기 위한 주요 가정 중 하나인 proportional hazard assumption을 만족하는지를 Kaplan curve가 rough하게 보여주기 때문이다. 비례위험가정은 생존분석의 가장 중요한 가정이기 때문에 뒤에서 따로 다룰 예정이다. 우선 Kaplan curve에 중심을 두고 진행해 보자.

 

library(survival)
library(asaur)
library(dplyr)
library(ggplot2)

  우선 생존분석을 다루기 위해서는 시간을 다루는 것에 익숙해져야 한다. 대개 질병과 사망을 다루는 경우, 연구대상자는 일정 시점(time 1)에 특정 질환이 진단되어 연구대상자에 포함되게 된다. 그 후 일정 시간이 흐른 후(time 2) 연구대상자는 사망하거나(event) 다른 이유로 연구대상에서 빠지게 된다(censoring).   

 

  A~D까지 4명의 연구대상자 가 있다고 가정해보자. 대상자 A는 2002년(time1)에 연구에 진입한 후 2008년(time2) 사망(event)하였다. 대상자 B는 2001년(time1) 진단 후 2005년(time2) 연구대상에서 loss (censoring)되었다. 대상자 C는 2002년(time1) 진단 후 2010년(time2) 연구 종료 때까지 생존(censoring)하였다. 이를 R로 표현해보면 (대개 status는 censoring을 0, 사망을 1로 입력하는 것이 일반적이다.) 

 

t1 <- c(2,1,2,4)
t2 <- c(8,5,10,9)
st <- c(1,0,0,1)
Surv(time = t1,time2 = t2, event = st)

  위와 같은 방식으로 표현되며, 이를 Kaplan curve로 그려 보면 아래와 같이 표현할 수 있다.

 

plot(survfit(Surv(time = t1,time2 = t2, event = st)~1))

  

  중간에 알아챘을 수도 있을 것 같은데, 사실 위와 같이 단순한 경우는 시간이 2개가 들어갈 필요가 없다. 만약 진단시점(time1)을 0으로 고정하게 되면

tt <- c(6,4,8,5)
st <- c(1,0,0,1)
Surv(tt,st)

 

훨씬 단순히 표현 가능하다.

plot(survfit(Surv(tt,st)~1))

  언뜻 두 그래프가 다르다는 것을 눈치 챈 사람도 있겠지만 사실 동일한 통계량을 지닌다. 다만 진단 시점을 0으로 고정했을 따름이다. 첫 번째 그림을 다시 보자. 10년 동안 추적한 대상자가 있는가? 그리고 대부분의 경우 후자와 같이 진단 시점을 0으로 고정한 후 분석을 진행한다. time1과 time2가 둘 다 필요한 경우는 나중에 다루겠지만 어쩔 수 없는 경우에 억지로 사용하는 경우가 대부분이다.    

 

 

  자 그럼 본격적으로 그림을 그려보자.

data("pharmacoSmoking")
head(pharmacoSmoking)
m1 <- survfit(Surv(ttr,relapse)~grp, data=pharmacoSmoking)
m1

plot(m1, col=c("red","blue"))
legend("topright", lty=1, col=c("red","blue"),legend = c("combination","patch only"))

  금연 후 실패에 관한 데이터이다. 금연 연구에 왠 생존분석? 이라고 생각할 수도 있지만 금연 후 실패를 사망으로 간주한다면 생존분석으로 분석 못 할 이유는 없다. 그런 경우 생존분석(survival analysis)라는 용어대신 event history analysis라는 용어를 쓰기도 하지만 의학연구에서는 상관없이 생존분석 용어를 사용하는 경우가 더 많다. conf.time 명령과 mark.time 명령을 통해 신뢰구간과 타임마크를 더해줄 수도 있다. 

plot(m1, col=c("red","blue"), conf.int = T, mark.time = T)
legend("topright", lty=1, col=c("red","blue"),legend = c("combination","patch only"))

  상위급 저널에서는 Kaplan curve를 제시하면서  아래와 같이 No. at risk를 제시하기를 원하는 경우도 있다 이런 경우 

summary(m1, times = c(0,50,100,150))

원하는 시간을 입력하면 그 시간에서의 no. at risk를 쉽게 구할 수 있다(median 값은 설명 도입부 첫 그림과 같이 "m1" 만 입력할 경우 구할 수 있음). 

 

  마지막으로 proportional hazard에 대한 언급을 해 보자. cox 분석은 여러 분석법 중 신경쓸 일 많지 않은 쉬운 분석법 중에 하나이다. 하지만 비례위험 가정은 꼭!!! 어기면 안된다.  비례위험 가정에 대해 예를 들어보면, 암환자의 경우 진단 후 초/중반에는 사망위험이 높지만 시간이 지날수록 사망 위험도는 점차 감소하다,  5년 이상 생존하는 경우 사망위험이 일반인보다 약간 높을 정도까지 낮아진다. 만약 새로운 항암제의 효과에 대해 연구하고 싶다면 신/구 항암제의 위험도 차이가 초반에나 후반에나 동일해야 한다는 것이다.     

 왼쪽 그림과 같이 시간에 상관없이 위험도는 동일한 간격을 가져야 한다. 만약 오른쪽 그림과 같이 위험도가 일정한 간격을 지니지 않을 경우 시간에 대해 별도의 배려를 해야만 한다.  

data("pancreatic2")
m2 <- survfit(Surv(os,status)~stage,data=pancreatic2)
plot(m2, col=c("red","blue"))
legend("topright",lty=1, col=c("red","blue"),legend = c("Locally Advanced","Metastatic"))

  위험도 그래프를 직접 그리는 것은 상당히 귀찮은 일이라서... Kaplan curve를 통해서 보면 두 선이 만나는 부분이 보인다(400부근 & 600부근). 만약 위와 같이 두 선이 만나는 접점이 있을 경우 비례위험 가정은 만족하지 못하는 것으로 간주되며, 따라서 시간에 대해 추가적인 배려를 하지 않은 연구결과는 인정받을 수 없다. 

 

중간 정리를 해보면

 

1. 접점 유무를 보기 위해 Kaplan curve를 요구하는 것이다.

2. 만약 접점이 있는 경우에도 무시하고 일반적인 방법으로 분석하면... 결과는 무조건 reject 이다.

 

 

 

 

2. Cox 분석

 

  Cox 분석의 정식 명칭은 proportional hazards model regression analysis이다. 정식 명칭을 보면 앞에서 왜 그토록  proportional hazard assumption을 강조했는지와, Cox 분석도 일반화선형회귀분석의 한 분류라는 것을 알 수 있다. 전체적인 분위기는 로지스틱회귀분석과 놀라울 정도로 닮아 있다. 잔차분석이 그리 중요하지 않다던가... 연속형 변수 사용에 유의해야 한다던가... 하여간 앞의 금연 데이터를 이용해서 Cox 분석을 맛봐보자.

head(pharmacoSmoking)

   데이터를 보면 대충 감이 오겠지만, 이 데이터는 금연치료에 패치만 쓴 그룹과 패치와 약물치료를 같이 한 그룹으로 나누어 금연치료의 효율성을 본 연구이다. 따라서 이 연구의 결과변수(Y)는 시간(ttr)과 재발여부(relapse)가 결합된 결과이며, 원인변수(X)는 치료그룹(grp)이 된다. 그 외 나이라던가 흡연년수...은 covariate로 보면 될 것이다.

with(pharmacoSmoking,hist(ttr))
with(pharmacoSmoking,table(relapse))
with(pharmacoSmoking,table(relapse,ttr))

 

  시간의 분포를 보면 연구의 대략적인 결과를 파악할 수 있다. 생존분석에서 시간은 사망하거나(event), 연구에서 중도탈락(censoring)한 두 가지 경우 중 하나이다. 시간의 histogram을 보면 전체적으로 right shifted 된 느낌인 가운데 0과 200근처가 우뚝 서있는 것을 볼 수 있다. 0에 가까운 봉우리는 아마도 작심삼일 들일 것이며, 200에 가까운 봉우리는 연구 종료에 따른 censoring일 것이다. 아래 table을 보면 그 정황을 확인 가능하다. 그리고 아래 table을 보면 연구 종료를 제외하면 censoring이 없음을 알 수 있다. 이는 이 연구가 굉장히 잘 관리된 연구 결과임을 보여준다.       

 

  대부분의 연구에서 무시되고 있는 내용인데, 생존분석에서 중도탈락은 무작위(random)하게 이루어지는 것으로 가정된다. 하지만 실제 연구에서는 중도탈락이 의미를 지니는 경우가 많다. 치료약의 부작용이 심하거나, 환자가 스스로 판단하기에 호전이 없을 경우, 무작위 배정 결과가 중간에 환자에게 노출되어 위약군 환자가 연구에서 빠져나가는 경우... 어찌 보면 중도탈락은 사망결과만큼이나 중요한 결과이지만 일반적인 cox분석 결과에서는 중도탈락을 결과물로 제시하지 않는다. 하지만 연구 결과를 판단할 때 신약의 치료효과가 아무리 좋더라도 이해할 수 없는 이유로 치료군 혹은 placebo군의 중도탈락이 반대편에 비해 의미가 있을 정도로 많다면 그 결과는 신뢰하지 않는 것이 바람직하다. 

 

m1 <- coxph(Surv(ttr,relapse)~grp,data=pharmacoSmoking)
summary(m1)

  Cox분석 결과를 보면 피치만 사용한 경우 약과 같이 치료한 경우보다 금연에 실패할 가능성(상대위험도)이 1.83배(95% CI: 1.20~2.80)높다고 나온다(로지스틱회귀분석에서는 exp(β)가 오즈비가 되지만 생존분석에서는 상대위험도로 불린다.). 그리고 원하는 결과의 기준이 복합치료라면 exp(-β)값인 0.55를 읽는 방법도 있고 아래와 같이 reference 값을 바꿔주면 원하는 결과를 얻을 수 있다.

pharmacoSmoking$grp <- relevel(pharmacoSmoking$grp, ref = "patchOnly")
m1 <- coxph(Surv(ttr,relapse)~grp,data=pharmacoSmoking)
summary(m1)

 

 

 

 

3. Martingale residual -  Cox 분석에서 연속변수의 처리법

  로지스틱회귀분석에서와 같이 생존분석에서도 연속변수를 분석하는데에는 주의가 따른다. 따라서 적절한 값을 기준으로 잘라서 범주형변수로 분석하는 것이 합리적이다. 하지만 대상자의 수가 많지 않아서 범주형 변수로 잘랐을 때 통계적으로 유의해지지 않는다던가, 범주형 변수로 변환하더라도 어느 값을 기준으로 몇 개의 범주로 나누는 것이 합리적인지 판단할 때에는 연속형 변수를 이용한 분석이 도움이 될 때도 있다. 이 때 사용되는 것이 martingale residual 이다.

with(pharmacoSmoking, hist(age))

   age변수의 histogram을 보면 정확히 정규분포를 만족해 보이지는 않지만 이 정도는 특별히 변환 필요 없이 분석해도 될 것 같은 느낌이 든다 (이런 느낌적인 느낌이 싫으면 정규분포 확인해봐도 된다.) 

다음 과정은 martingale residual 계산이다.

m1 <- coxph(Surv(ttr,relapse)~1,data = pharmacoSmoking)
pharmacoSmoking$resid <- residuals(m1, type = "martingale")
pharmacoSmoking %>% ggplot(aes(age,resid))+geom_point()+geom_smooth()+theme_classic()

  이전 선형회귀분석 포스트에서 다뤘던 내용이기는 한데, 잔차는 회귀분석에 투입된 변수를 보정한 후 남은 값이다. 생존분석에서는 이 역할을 martingale 잔차가 맡고 있다. 따라서 martingale 잔차는 투입된 변수들을 보정한 후 hazard값이 된다. 위의 cox분석에서는 보정할 변수를 넣지 않았으므로 resid=0는 평균 hazard이며, +값은 hazard가+, -값은 반대를 말한다. 위 그림의 파란색 선은 loess 분석선인데 이를 보면 40대 이전은 hazard 값이 높으며, 그 후 어느 정도 일정한 값으로 유지되다, 50대 이후 감소되는 모습을 보인다. 이런 양상을 이해하고 있으면 본 data에서는 age를 agegroup2와 agegroup4 2종류 변수로 범주화시켜 놓았는데 그 이유를 이해할 수 있을 것이다.       

 

with(pharmacoSmoking, hist(longestNoSmoke))
with(pharmacoSmoking, hist(log(longestNoSmoke+1)))
pharmacoSmoking %>% ggplot(aes(log(longestNoSmoke+1),resid))+geom_point()+geom_smooth()+theme_classic()

  이전 최대금연일 같은 경우는 log 씌워주면 상태가 좋아져서 이대로 확인해 보았다. loess 그래프상 완전한 직선은 아니지만 중후반부는 완만한 직선을 그리고 있으며 그 선을 연장해보면 초반에도 90% 신뢰구간을 완전히 벗어나지는 않는다. 따라서 억지로 우긴다면 연속형 변수 형태로 투입 가능하다(만약 이 변수가 주요 변수라면 이런 짓은 해서는 안된다. 하지만 어차피 보정을 위한 covariate이기 때문에... 다른 말로 하면 엑스트라이기 때문에 이런 짓이 용납될 가능성이 있는 것이다. 만약 대상자 수가 충분해서 모든 covariate를 범주화해도 통계적 유의성을 유지하는 경우는 당연히 이런 짓은 필요 없다.).  

m1 <- coxph(Surv(ttr,relapse)~log(longestNoSmoke+1),data = pharmacoSmoking)
summary(m1)

m1 <- coxph(Surv(ttr,relapse)~cut_number(longestNoSmoke,n = 3),data = pharmacoSmoking)
summary(m1)

연속형변수로 처리한 결과와 3개의 범주를 지니는 범주형변수로 변환한 후의 결과이다. 어떤 결과가 더 어울릴지는 상황에 따라 판단해야 한다. 

 

 

 

4. Scheonfeld residual - proportional hazard assumption

  Kaplan curve 설명할 때 proportional hazard assumption의 중요성에 대해서 이미 언급한 바 있다. 생존분석을 제대로 하느냐 못하느냐는 이 가정을 얼마나 이해하고 있느냐에 달려있을만큼 중요한 내용이다. 

 

  proportional hazard 가정을 확인하는 방법으로는

 

1) Kaplan curve 이용한 육안적 확인

2) Schenofeld residual 그래프 이용한 육안적 확인

3) Scheonfeld method 이용한 통계적 확인

4) time dependent covariate 투입 후 통계적 유의성 확인 

 

  4가지 방법이 있다. 네 번째 방법은 진단법 이라기보다 치료법에 가까우니 우선 1~3까지 방법을 알아보자.

 

 

 

4-1) Kaplan curve 이용한 육안적 확인

 

m1 <- survfit(Surv(ttr,relapse)~ageGroup4,data=pharmacoSmoking)
plot(m1, col=c("red","blue","green","violet"))

ageGroup4 변수를 이용한 Kaplan curve이다. 어째 시작부터 암울한 기운이 드는 것이... 자세히 보면 군데군데 교차점들이 보인다. 상당히 위험한 상황이며 만약 주요변수를 이용한 Kaplan curve가 이모양이면...

하여간 Kaplan curve에서 선끼리 교차가 이루어지면 위험신호라고 받아들이는 게 옳다.

 

 

 

4-2)  Scheonfeld method

 

Kaplan curve만 가지고 단정 짓는 일은 약간 빠르니 다른 방법을 통해 확인해 보자.

m2 <- coxph(Surv(ttr,relapse)~ageGroup4,data = pharmacoSmoking)
summary(m2)
cox.zph(m2)

Scheonfeld method는 Scheonfeld residual을 구해서 그 통계량을 확인하는 방법이다. cox.zph 명령을 이용해서 통계량을 얻을 수 있다. 다행히 ageGroup4 변수는 proportional hazard 가정에 어긋나지 않는다 (p=0.351).

 

 

 

4-3) Scheonfeld residual을 이용한 시각적 방법

 사실 시각적으로 확인한 후 통계량을 봐야 하는데, 순서가 바뀐 듯도 하지만...  R에서 구하는 순서에 따르기로 했다.

plot(cox.zph(m2))

martingale 잔차는 시간의 개념을 포함하고 잇는데 비해서 Sceonfeld 잔차는 시간의 개념을 포함하지 않는다. 따라서 Sceonfeld 잔차의 시간에 따른 기울기가 0이라면 hazard 또한 시간의 변화에 따라 바뀌지 않는다 말할 수 있다. Sceonfeld method는 위의 그래프들의 기울기가 0인지 아닌지를 알려주는 통계량이다. 

 

data("lung")
m3 <- coxph(Surv(time,status)~ph.karno,data=lung)
m4 <- cox.zph(m3)
summary(m3)
m4
plot(m4)
abline(h=0, col="red", lty=2)
abline(h=-0.016,col="blue")

 

다른 예를 들어보자. physician이 측정한 initial Karnofsky score에 관한 내용이다. Karnofsky score는 폐암환자의 생존에 영향을 미치는데 Scheonfeld method에서의 통계량(p=0.0048)은 시간에 따라 hazard가 변함을 알려주고 있다. 그래서 확인한 아래의 그림에서는 대략 time 200까지는 높은 Karnofsky score가 낮은 위험도를 보이지만 (exp(-0.016)보다 낮다.),  time 270 이후에는 위험도의 95% 신뢰구간이 1 (exp(0))에 걸치게 된다. 다른 말로 표현하면 Karnofsky score가 높은 경우 대략 time 270까지는 낮은 위험도를 보이지만 그 이후에는 유의하지 않게 된다고 설명 가능하다.    

 

 

 

4) time dependent covariate 투입 후 통계적 유의성 확인 

m5 <- coxph(Surv(time,status)~ph.karno+tt(ph.karno),data=lung)
m5

  마지막 방법은 time dependent covariate를 투입하는 방법이다.  만약 시간의 변화에 따라 hazard가 변화한다면 time*covariate의 상호작용은 통계적으로 유의할 것이다. 우선 이 방법의 긍정적인 측면은 비례위험가정의 진단 및 치료라는 점이다. 만약 ph.karno의 time transfer function (tt)가 유의하다면 비례위험가정이 어긋난다는 의미이며, 그 해결방안은 tt(ph.karno) 항목을 남겨놓는 것이다. 다만 문제는... Scheonfeld residual을 이용한 방법에 비해 신뢰도가 떨어지며, 좀 더 정확하게는 time, log(time), rank(time)... time을 변형한 수치도 투입해서 확인해야 한다는 점이다. (위의 예제는 time transfer function에 아무 옵션을 주지 않았기 때문에 rank(time)이 기본치이다.)   

 

  가장 현실적인 방법은 Scheonfeld residual graph와 통계량을 이용하는 방법이다. 그럼에도 불구하고 네 가지 방법을 모두 언급한 이유는 1) Kaplan curve의 경우 선끼리 교차가 있으면 다른 방법으로 확인을 했다고 해도 쓸데없는 오해를 남길 수 있다는 점, 2) 일부 통계 패키지의 경우 Scheonfeld residual을 구할 수 없다는 이유이다.  

 

 

 

 

5. non-proportional hazard의 교정

 

   비례위험 가정이 깨진 경우 이를 해결하기 위한 방법은 

 

1) 연구기간을 줄일 수 있는지 확인한다.

2) step function을 이용

3) time dependent covariate 투입

 

  의 3가지 방법이 현실적이다.

 

4) 만약 주요변수가 아니라면 stratification method 쓰면 된다.  

 

 

1) 연구기간 확인

 

  아무래도 연구기간이 길어지면 비례위험가정이 깨지는 경우가 흔하다. 당장 위의 폐암 예를 봐도 진단 당시 Karnofsky score는 초기 300일 까지 정도까지만 유효하다. 이때 만약 연구기간을 300일 까지로 변경한다면 비례위험가정은 만족되는 상태가 된다. 특히 암환자 연구에서 5년 생존을 보게 되면 변수 중 몇 개 이상은 비례위험 가정에 어긋나는 경우가 흔하다. 필요 없이 긴 추적관찰 기간이 연구에 오히려 해가 되는 경우가  있을 수 있다. 아래 2가지 방법이 어려운 경우 연구기간을 줄이는 방법이 가장 현실적이다.

 

 

 

2) Step function

 

  연구기간을 줄일 수 없는 경우 기간을 나누어서 위험도를 계산하는 방법이 있다.

  위의 예제를 다시 보면 대략 300일 이전과 이후의 양상이 다름을 알 수 있다. 이를 어떻게 나눌지는 전적으로 연구자의 역할이기는 한데...  270일(95% 신뢰구간의 윗 쪽이 0에 걸치는...)이전과 이후로 나누는 법, 아님 180/360/그 이후로 나누는 법이 있을 듯하다. 아님 365일 이전과 이후로 나누는 법도 있고... 기존 연구들의 결과와 달력(년, 4분기, 월 단위로 나누는 것이 합리적이다.)에 따라 나누자.

  편의상 180/360/그 이후로 나누어 보자.    

lung2 <- survSplit(Surv(time,status)~.,data=lung,cut=c(180,360),episode = "tgroup")
head(lung)
head(lung2)
table(lung2$tgroup)

lung2 데이터의 1~2 번째 줄은 lung 데이터의 1번 째  환자의 자료임을 알 수 있다. 306일 후 사망한 환자의 자료를 0~180일 censor & 180~306일 사망의 2가지 자료로 분리했음을 알 수 있다 (lung 데이터 두 번째 68세 환자의 자료는 lung2 데이터의 3~5줄의 자료로 분리되었음을 보자.). 

이 자료를 분석하자면 2개의 시간을 다루는 것이 필수이다. 그리고 tgroup으로 묶인 부분도 확인하자.

m1 <- coxph(Surv(time,status)~ph.karno,data=lung)
m1
cox.zph(m1)

m2 <- coxph(Surv(time = tstart,time2 = time,event = status)~ph.karno:strata(tgroup),data=lung2)
m2
cox.zph(m2)

 

 

model2의 결과를 보면 0~180일까지만 유의하고 그 이후는 유의하지 않음을 알 수 있고, Sceonfeld method에서도 비례위험 가정을 만족하는 모습을 확인할 수 있다(예를 들기 위해 3군으로 나눴지만 데이터의 실제 양상은 0~270일/그 이후 2군으로 나누는 것이 더 적절하다.).

lung3 <- survSplit(Surv(time,status)~.,data=lung,cut=270,episode = "tgroup")
m3 <- coxph(Surv(time = tstart,time2 = time,event = status)~ph.karno:strata(tgroup),data=lung3)
summary(m3)
cox.zph(m3)

 

 

3) time dependent covariate 투입

 연속형 time dependent covariate를 투입하는 방법도 있다. SPSS의 생존분석은 시간을 하나밖에 다룰 수 없으며, Schenofeld residual도 계산할 수 없기 때문에 time dependent covariate 투입이 유일한 비례위험가정의 진단 및 치료법이다. 때문에 쓸데없이 유명해졌지만 실은 사용이 굉장히 어려운 방법이다. 

 위쪽 LOESS 곡선을 마음의 눈으로 보면 β(t)=a*Time+b 의 선형회귀 방식으로 표현이 가능함을 알 수 있다.

m1 <- coxph(Surv(time,status)~ph.karno+tt(ph.karno),data=lung,tt=function(x,t, ...) x*t)
summary(m1)

m2 <- coxph(Surv(time,status)~ph.karno+tt(ph.karno),data=lung,tt=function(x,t, ...) x*log(t))

이를 표현하면 β(t)=7.22^-5*Time-3.86^-2 로 표현 가능하다. 이는 시간이 0일때의 상대위험도는 exp(-3.86^-2)=0.962 이지만 시간이 100일때의 상대위험도는 exp(100*7.22^-5-3.86^-2)=0.969로 변한다는 의미이다.

 

이 방법을 사용하기 위해서 몇 가지 주의할 점이 있다. 첫 번째로 β(t)를 a*Time+b 로만 가정하지 말고 최소한 a*log(Time)+b와 a*rank(Time)+b 3가지 방법을 시도해 보아야 한다. 두 번째는 time dependent covariate를 투입하게 되면 Schenofeld residual을 믿을 수 없는 경우가 많아진다. 케이스 바이 케이스로 대응하자.

 

 

 

 

6. Stanford Heart Transplant Data

 

  비례위험가정 이외에도 생존분석에 주의할 부분이 있다. jasa 데이터는 이전에 연구되었던 심장이식 연구의 관찰연구자료이다.

head(jasa)
m1 <- survfit(Surv(futime,fustat)~transplant,data=jasa)
plot(m1,col = c("red","blue"))

 심장이식 그룹의 생존이 훨씬 나은 것을 알 수 있다. 실제로 이대로 이 결과는 저명 학술지에 발표되었고, RCT까지 진행한 후에 효과가 없음이 밝혀졌다. 이는 심장이식 순서를 기다릴 수 있을 정도로 건강한 사람만이 심장이식군으로 포함되었기 때문이다 (심장이식 기다리다 죽은 사람은 control group). 심장수술을 기다리는 기간을 따로 분리하면 다음과 같아진다.

m2 <- survfit(Surv(time = start,time2 = stop, event = event)~transplant,data=jasa1)
plot(m2,col = c("red","blue"))

jasa1 데이터를 보면  심장이식 환자의 경우 수술 전 대기시간을 따로 분리한 자료이다. 

이러한 상황은 관찰연구에서 굉장히 흔한 상황이며, 아직까지도 수많은 연구에서 동일한 오류가 반복되고 있다. 비록 이 포스트에서는 이런 자료를 어떻게 처리할지 다루지 않겠지만, 관찰연구에서는 치료까지의 대기시간의 유무가 연구결과에 영향을 끼칠 수 있음을 알고 있어야 한다.     

 

 

 

 

 

7. Cox 그래프

 

  대개의 경우 논문에 넣을 그림으로 Kaplan curve를 사용하지만 Cox curve를 넣는 경우도 종종 보인다. 앞에서도 언급했다시피 Kaplan curve는 비례위험가정 만족 여부를 확인시켜주기 때문에 흔히 사용되지만, expert 들이 big journal에 싣는 경우 비례위험가정 확인이 그다지 필요 없기 때문에 결과를 더 잘 보여주는 Cox graph를 선호하기도 한다.

 

m1 <- coxph(Surv(ttr,relapse)~grp+age+employment,data=pharmacoSmoking)
p1 <- survfit(m1,newdata = data.frame(grp=c("combination","patchOnly"),age=40,employment="ft"))
plot(p1,lty=1:2)
legend("topright",legend = c("Combination","Patch Only"), lty=1:2)

 

8. 마무리하며

 

  생존분석은 상대적으로 신경 쓸 부분이 적기는 하지만 비례위험가정 만족 여부와 관찰 연구의 경우 치료까지의 대기시간 같은 결과에 영향을 미칠 수 있는 부분은 신중하게 접근해야 한다.  그 외에는 로지스틱회귀분석과 비슷한 방법으로 접근하면 큰 문제는 없을 것이다.