Survival analysis, lab 3
yongchan Kwon, ykwon.0407@gmail.com
Summary : 다음을 R을 통하여 실습합니다. (1) Review (2) Two-sample testing을 해봅니다.
Review
Lecture note chapter 4에서 Likelihood Inference를 배웠습니다. Right censoring이 포함된 data에서 likelihood function을 다음과 유도할 수 있었습니다. \[ L(\theta; x, \delta) \propto \prod_{i=1}^n h(x; \theta)^{\delta_i}S(x; \theta) \] Likelihood fucntion이 well-defined됬으니 \(\theta\)를 추정할 수 있었습니다. (어떻게 할 수 있을까요? 생각해봅시다!)
Chapter 5 에서는 chapter 3에 이어 counting process를 배웠습니다. Data with right censoring (\(\{ (T_i, \delta_i) \}_{i=1, 2, \ldots, n}\))에서 Process들을 다음과 같이 정의했습니다. \[ Y(t) := \sum_{i=1}^{n} I(T_i \geq t), \quad N(t) := \sum_{i=1}^{n} I(T_i \leq t, \delta_i = 1) \\ \Lambda (t) := \int _0 ^t \lambda(s) ds, \quad M(t) := N(t) - \Lambda (t) \] 위의 주요 네가지 process 정의들로 부터 Martingale 이론(Central Limit Theorem)을 생각할 수 있었고 \(\sqrt{n} [\hat{H}(t)-H^*(t)]\) 분포를 추정할 수 있었습니다.
오늘 실습시간에 할 내용은 Lecture note chapter 5.2: Hypothesis testing입니다. 가설검정을 하기위해서는 적절한 가설이 필요합니다. 우리는 수업시간에 가설을 세워서 test statistic을 공부했습니다. 한 번더 살펴봅시다.
가설 : \(H_0 : h_1(t) = h_2(t)\) for all \(t > 0\) versus \(H_1 : h_1(t) \neq h_2(t)\) for some \(t > 0\).
아래 모든 통계량들은 \(H_0\) 하에서 \(\chi^2 (1)\) 분포를 따른다. Reject \(H_0\) for large values.
Test statistics.
Right-censored data에 그룹 정보\(Z_i\)를 추가하여 \(\{(T_i, \delta_i, Z_i)\}_{i=1}^{n}\) (\(Z_i=1\) or \(2\)) 을 관측하였다고 가정하고, 이로부터 (at-)risk process와 event \(Y_1(t)\), \(Y_2(t)\), \(Y(t)\), \(dN_1(t)\), \(dN_2(t)\), \(dN(t)\) 등을 정의합시다. 그리고 \(t_1 < t_2 < \ldots < t_D\)는 \(n\)개의 관측치를 중에 uncensored observation 이 있는 (i.e., event 혹은 failure가 있는) 시점을 가리킨다 하자.
Log-rank test : \[ Z^2 := \frac{ \left[ \sum_{k=1}^{D} \! \left\{ dN_2(t_k) - Y_2(t_k)\frac{dN(t_k)}{Y(t_k)} \right\} \right]^2 } { \sum_{k=1}^{D}\! Y_2(t_k) \frac{dN(t_k)}{Y(t_k)} \frac{Y(t_k)-dN(t_k)}{Y(t_k)} \frac{Y_1(t_k)}{Y(t_k)-1} } \] Weighted log-rank test : \[ Z_W^2 := \frac{ \left[ \sum_{k=1}^{D} W(t_k) \left\{ dN_2(t_k) - Y_2(t_k)\frac{dN(t_k)}{Y(t_k)} \right\} \right]^2 } { \sum_{k=1}^{D} W(t_k)^2 Y_2(t_k) \frac{dN(t_k)}{Y(t_k)} \frac{Y(t_k)-dN(t_k)}{Y(t_k)} \frac{Y_1(t_k)}{Y(t_k)-1} } \] 여기서 weight \(W(t)\)에 따라 test의 이름이 달라지는데, \(W(t) \equiv 1\) (log-rank), \(W(t) = Y(t)\) (Gehan), \(W(t) = Y(t)^{1/2}\) (Tarone-Ware), 그리고 \(W(t) = \hat{S}(t^-)^p (1-\hat{S}(t^-))^q \) (Fleming-Harrington) 등이 있습니다.
Some remarks
위의 통계량들을 자세히 살펴보면 \[ \frac{\{\sum({\rm Observation} - \rm Expectation)\}^2}{\sum \rm Variance} ~~{\rm or}~~\frac{({\rm Observation} - {\rm Expectation})^2}{\rm Variance}\] 이러한 형태임을 알 수 있습니다.
당연한 얘기지만 위에서 그룹 1과 2의 역할을 뒤집어도 equivalent한 통계량을 얻습니다.
위의 two-sample testing 통계량이 \(Z^2 / V\)의 꼴로 보일 수도 있다. \(K\)-sample testing에서는 \(K-1\) dimension의 벡터 \(\bf Z\)를 quadratic form \({\bf Z}^{\top} {\bf V}^{-1} {\bf Z}\) 으로 변형한 형태가 test statistic이 되며 null distribution은 \(\chi^2(K-1)\)이 된다.
R의 survival 패키지에서 함수survdiff()는 \(W(t) = Y(t)^{\rho}\) 꼴의 test를 지원합니다. 그럼 실습을 시작해봅시다.
Two-sample testing
Two-sample testing은 survival 패키지의 survdiff() 함수를 통해서 할 수 있습니다. survival 패키지를 설치하고 로드합시다.
#install.packages("survival") ;
library(survival) ;## Loading required package: splinesworking directory를 설정하고 eTL에 있는 kidney dataset을 다운 받아봅시다. Textbook chapter 1.4에서 데이터에 대한 설명과 데이터를 볼 수 있습니다. (txt 파일을 확인해봅시다!) 저번 실습시간에 배운대로 Surv object를 만들어봅시다.
setwd("C:\\Users\\snu\\Downloads\\"); # working directory setting
kidney = read.table("kidney.txt", header=T, sep="\t") ;
survobj.kidney = Surv(time = kidney$survtime, event = kidney$censoring) ;
survobj.kidney;## [1] 1.5 3.5 4.5 4.5 5.5 8.5 8.5 9.5 10.5 11.5 15.5
## [12] 16.5 18.5 23.5 26.5 2.5+ 2.5+ 3.5+ 3.5+ 3.5+ 4.5+ 5.5+
## [23] 6.5+ 6.5+ 7.5+ 7.5+ 7.5+ 7.5+ 8.5+ 9.5+ 10.5+ 11.5+ 12.5+
## [34] 12.5+ 13.5+ 14.5+ 14.5+ 21.5+ 21.5+ 22.5+ 22.5+ 25.5+ 27.5+ 0.5
## [45] 0.5 0.5 0.5 0.5 0.5 2.5 2.5 3.5 6.5 15.5 0.5+
## [56] 0.5+ 0.5+ 0.5+ 0.5+ 0.5+ 0.5+ 0.5+ 0.5+ 0.5+ 1.5+ 1.5+
## [67] 1.5+ 1.5+ 2.5+ 2.5+ 2.5+ 2.5+ 2.5+ 3.5+ 3.5+ 3.5+ 3.5+
## [78] 3.5+ 4.5+ 4.5+ 4.5+ 5.5+ 5.5+ 5.5+ 5.5+ 5.5+ 6.5+ 7.5+
## [89] 7.5+ 7.5+ 8.5+ 8.5+ 8.5+ 9.5+ 9.5+ 10.5+ 10.5+ 10.5+ 11.5+
## [100] 11.5+ 12.5+ 12.5+ 12.5+ 12.5+ 14.5+ 14.5+ 16.5+ 16.5+ 18.5+ 19.5+
## [111] 19.5+ 19.5+ 20.5+ 22.5+ 24.5+ 25.5+ 26.5+ 26.5+ 28.5+Testing : survdiff(formula, rho=0)
이번 실습ㅇ베서는 two-sample testing만 다루겠습니다만 역시 같은 문법으로 \(K\)-sample testing에도 적용할 수 있습니다.
formula에는 “(survival object) \(\sim\) (group-indicating vector)” 꼴을 입력하면 됩니다. rho는 weight 함수 \(W(t) = Y(t)^{\rho}\)에서의 \(\rho\)에 대응됩니다. 이를테면 아래 셋은 rho=0을 지정하였으므로 log-rank test의 결과를 보여 줍니다. 문법상에 미묘한 차이가 있으나 모두 같은 결과를 저장합니다.
test.kidney = survdiff(survobj.kidney ~ as.factor(kidney$group), rho=0) ;
test.kidney = survdiff(survobj.kidney ~ kidney$group, rho=0) ;
test.kidney = survdiff(survobj.kidney ~ group, rho=0, data=kidney) ;그럼 결과를 살펴봅시다. 주요 결과가 (Review에서 보았던) Observation, Expectation, Variance의 형태로 확인할 수 있습니다.
test.kidney## Call:
## survdiff(formula = survobj.kidney ~ group, data = kidney, rho = 0)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## group=1 43 15 11 1.42 2.53
## group=2 76 11 15 1.05 2.53
##
## Chisq= 2.5 on 1 degrees of freedom, p= 0.112names(test.kidney);## [1] "n" "obs" "exp" "var" "chisq" "call"test.kidney$obs;test.kidney$exp;test.kidney$var;test.kidney$chisq## [1] 15 11## [1] 11.04 14.96## [,1] [,2]
## [1,] 6.211 -6.211
## [2,] -6.211 6.211## [1] 2.53[생각해보기] Review에서 보았던 검정통계량은 1차원인데 왜 위 결과의 variance matrix는 2 by 2일까요?
survfit() revisited
survfit() 또한 input formula에 “(survival object) \(\sim\) (group-indicating vector)”의 형태를 사용할 수 있습니다.
우선 fit.kidney에 survfit object를 저장하고 다양한 결과를 확인해봅시다.
fit.kidney = survfit( survobj.kidney ~ group, data=kidney ) ;
fit.kidney # 그룹별로 survival function이 나뉘어 적합되어 있습니다.## Call: survfit(formula = survobj.kidney ~ group, data = kidney)
##
## records n.max n.start events median 0.95LCL 0.95UCL
## group=1 43 43 43 15 18.5 15.5 NA
## group=2 76 76 76 11 NA NA NAnames(fit.kidney)## [1] "n" "time" "n.risk" "n.event" "n.censor"
## [6] "surv" "type" "strata" "std.err" "upper"
## [11] "lower" "conf.type" "conf.int" "call"fit.kidney$surv## [1] 0.9767 0.9767 0.9523 0.8994 0.8722 0.8722 0.8722 0.8024 0.7659 0.7276
## [11] 0.6872 0.6872 0.6872 0.6872 0.6247 0.5623 0.4998 0.4998 0.4998 0.3748
## [21] 0.3748 0.1874 0.1874 0.9211 0.9211 0.8882 0.8700 0.8700 0.8700 0.8452
## [31] 0.8452 0.8452 0.8452 0.8452 0.8452 0.8452 0.8452 0.7848 0.7848 0.7848
## [41] 0.7848 0.7848 0.7848 0.7848 0.7848 0.7848 0.7848summ.kidney = summary(fit.kidney) ;
names(summ.kidney)## [1] "n" "time" "n.risk" "n.event" "n.censor"
## [6] "surv" "type" "strata" "std.err" "upper"
## [11] "lower" "conf.type" "conf.int" "call" "table"summ.kidney$table## records n.max n.start events median 0.95LCL 0.95UCL
## group=1 43 43 43 15 18.5 15.5 NA
## group=2 76 76 76 11 NA NA NA# try various expressions끝으로, plot(fit.kidney)를 입력하여 보면 두 그룹별로 나뉘어 그려진 곡선을 볼 수 있습니다.
plot(fit.kidney) ;
plot(fit.kidney, mark.time=FALSE) ; # censored observation의 표시여부를 결정
# 사실 censored observation mark의 모양 (기본은 십자(+)) 및 크기/색상도 조절가능
plot(fit.kidney, col=c("red", "blue")) ; # 그룹별로 생존함수 곡선의 색을 변경
plot(fit.kidney, conf.int=TRUE, col=c("red", "blue")) ;
# 그룹별 (pointwise) CI를 추가
plot(fit.kidney, mark.time=FALSE, col=c("red", "blue"),
xlab="t", ylab="S(t)", main="Estimated (infection-free) survival function") ;
# 그래프 제목, x축 이름, y축 이름 등 plot()의 기본 지원 옵션을 추가할 수 있습니다.
legend("bottomleft", col=c("red", "blue"), lty=c(1,1),
c("Surgical (GP 1)", "Percutaneous (GP 2)") ) ;
# 그려진 그래프 위에 범례를 추가하는 건 R 그래픽의 기본구조상 자명하게 가능합니다.Disclaimer : 이 자료는 최영근 조교님의 2013년 가을 Survival anlysis 실습 수업자료를 바탕으로 만들어졌습니다.