Introductory statistics with R: chapter 14

---

Survival analysis

 

14.1 Essential concepts

 

생존분석의 큰 특징은 censored data를 포함한다는 것이다. 일정 기간동안 환자군을 관찰한다고 했을 때, 특정 event(사망 혹은 재발 등)가 발생하는 경우와 발생하지 않는 경우 이외에도 환자가 단순히 follow up loss가 되거나 다른 원인으로 사망하는 등의 경우가 있을 수 있고 이를 censored data 라고 함.

 

Survival function: S(t) = 1-F(t)

특정 시간(기간)에 생존해있을 확률이며, event 발생에 대한 누적확률함수 F(t)를 1에서 빼준 값이 된다.

 

Hazard function: h(t) = f(t)/S(t)

해당 t 시간에 생존한 환자가 아주 짧은 기간 내에 사망할 확률. 

 

14.2 Survival obejcts

생존분석에 사용하는 함수들은 survival 이라는 library에 포함되어 있다.

 

* Surv object

생존분석의 기본이 되는 object로, f/u 기간과 censoring 관련 정보를 혼합한 object이다. ISwR 패키지의 melanom 이라는 데이셋을 예로 들어 사용하는데,

 

melanom은 악성 흑색종(malignant melanoma) 수술 후 환자 생존을 조사한 데이터이다.

no는 환자 식별번호, status는 end of follow up 시 환자의 상태(1: melanoma로 사망, 2: 생존, 3: 다른 원인으로 사망), 그리고 days는 f/u days를 말하며 ulc는 melanoma에 ulceration이 있었는지 여부 (1: 있음, 2: 없음), thick는 tumor의 thickness(1/100 mm), 그리고 sex는 환자 성별(1: 여자, 2:남자) 이다.

 

attach(melanom) 을 한 뒤, Surv object를 다음과 같이 만들 수 있다.

Surv(days, status==1)

argument로 days of follow up과 end of f/u 시 환자 상태(특히, 'uncensored' 상태가 어떤 것인지 지정)를 넣어준다.

 

그럼 아래와 같은 결과가 출력되는데, +가 붙은 것은 censored data(즉, alive 혹은 died from other caues)를 의미한다. 이 자체만으로는 알아낼 수 있는 정보가 많지 않고, KM 생존곡선을 그려서 시각화하는 것이 유용하다.

 

14.3 Kaplan-Meier estimates

 

survfit(Surv object) 함수를 이용해 KM estimator를 적합할 수 있다.

 

위 결과를 보면, survfit 객체는 event 발생을 시간순으로 배열하고 각각의 시점에서 survival과 95% CI를 계산해 가지고 있다. plot을 통해 이를 시각화해볼 수 있다.

 

plot(surv.all)

이는 전체 데이터에 대한 생존곡선이며, conf.int=FALSE 를 사용하면 CI 곡선을 제거 가능함.

lm과 마찬가지로, survfit 함수는 model formula를 사용할 수 있다.

 

surv.bysex <- survfit(Surv(days, status==1)~sex)

plot(surv.bysex, col=c('black','grey'))

legend('topright', legend=c('female','male'), col=c('black','grey'),lty='solid')

14.4 The log-rank test

KM을 시각화해서 female과 male이 다른 생존곡선을 보이는 것 같다는 clue를 얻었다. 이 둘의 차이가 통계적으로 유의미한 수준인지는 log-rank test를 통해 확인할 수 있다. Log-rank test는 간략히 말하면, 각 event 발생 시점에서 기대빈도 (관찰대상수에 비례)를 구해서 모두 합하고, 실제 관찰빈도와의 차의 제곱을 구하여 기대빈도로 나눠주면 chi-square 통계랑이 구해지고, 이를 통해 p-value를 구할 수 있다. 마치 각 칸의 기대도수와 실제 도수의 차의 제곱을 통해 chi-square 검정을 하는 것과 유사하다.

아래 함수를 이용해 test를 시행할 수 있다.

survdiff(Surv(days,status==1)~sex)

층화 분석(stratified analysis): 층화된 로그순위검정

결과에 영향을 줄 수 있는 "층화된" 범주로 묶인 경우, 이를 보정해 분석을 시행할 수 있다.

이 경우, 아래를 보면 p값의 유의성이 사라졌는데 이는 "남자는 질병이 진행된 상태일 수록(ulceration의 단계가 높을 수록) 더 적극적인 치료를 받는 경향이 있다" 정도로 설명될 수 있을 것이다.

14.5 The Cox proportional hazards model(phm)

로그순위법은 단변수에 대한 분석만 가능하다. 따라서 생존곡선에 대한 다변수의 영향을 파악하기 위해서는 콕스의 비례위험모형을 사용한다.

 

비례위험모형은 위험함수비율(hazard ratio)를 통해 독립변수가 위험률(hazard rate)에 미치는 영향을 보인다.

 

위험률(hazard rate): t 시점까지 생존한 사람이 그 직후 사망한 확률, 즉 순간 사망률의 개념.

앞에서 h(t) = f(t)/S(t)로 정의를 했으며, h(t)를 시간에 영향을 받는 부분과 covariate에 영향을 받는 부분으로 나누어 아래와 같이 바꿔쓸 수 있다. 

위험함수비율(HR)는 h(t)의 비율로 정의된다. Cox의 비례위험모형에서는 "비례위험 가정"을 전제로 하는데, 이는 쉽게 말하면 covariate에 대해 영향을 받는 요소가 아닌, 시간에 따라 영향을 받는 부분은 동일하다는 것이다. 이 가정을 만족하게 된다면 서로 다른 변수에 대해 각각 위험률을 구할 때, 그 비율로 정의되는 위험함수비율(HR)은 동일한 시간 요소(h0(t))가 상쇄되며 시간에 영향을 받지 않는, covariate에 영향을 받는 부분의 비율로만 남는다.

비례위험 가정: HR이 시간에 관계 없이 일정하다는 가정을 만족해야 한다.

만약 비례위험가정이 만족되지 않는다면 층화된 비례위험모형을 사용해야 한다.

 

coxph(Surv object에 대한 model formula) 함수를 이용해 cox object를 만들 수 있다.

 

lm과 마찬가지로, 어떤 모형을 선택할지는 step() 함수 등을 이용하여 설명력, 모형의 복잡도 등을 고려해 최종 선택을 할 수 있다. 그리고 아래 예시는 ulc에 대해 비례위험가정을 만족하지 못해 층화 비례위험모형을 만들었다.

 

비례위험모형을 만족하는지는 coxph object에 대해 log-log plot을 그리거나, cox.zph() 함수 등을 이용하면 알 수 있다.