Summary : 다음을 R을 통하여 실습합니다. (1) Review (2) Left-truncation이 있는 right-censoring 생존 자료 (LTRC data)를 분석해봅시다.

Review

먼저 left truncation과 left censoring의 차이를 생각해봅시다. 편의상 no censoring을 가정하고 time-to-event data를 \(\{ X_i \}_i\)라 합시다. 개체 \(i\)에 대한 또다른 관측시점 \(L_i\)에 대하여 \(X_i \leq L_i\)이면,

Left censoring :

e.g. 처음 마리화나를 피었을 때. 다음과 같은 대답이 가능합니다 : 17세 때부터 (observation), 피운 적 없다 (right censoring), 피운적은 있지만 기억나지 없다 (left censoring)

Left truncation :

e.g. Fyn county diabetes. 1973년 7월 1일에 연구가 시작될 당시 이 지역에는 1499명의 당뇨병 환자가 있었다. 이들은 1983년 1월 1일까지 정기적으로 검진을 받거나, 그 전에 죽거나, 혹은 이사를 갔다. 각 환자에 대해 당뇨병의 합병증이 발병한 때는 알려져 있었다. 다음과 같이 데이터를 생각해 봅시다. :

\(O\) : time of onset of diabetes. 당뇨병 합병증 발생 시점.

\(D\) : time of death. 사망 시점.

\(T\) : survival time from onset to death. 합병증 발병 이후 사망까지의 생존 시간. (\(T=D-O\))

\(L\) : length of time from onset to the start of study. 합병증 발병 이후 연구 시작까지 걸린 시간.

이 연구는 \(L_i \leq T_i\)인 사람들만 관측되었습니다. 직관적으로는 합병증에 걸려도 잘 안 죽는 사람들만 관측되었을 가능성이 있습니다.

[생각해보기] 동일한 자료라도, 보는 관점에 따라 right-censored data가 될 수도 있고 left-truncated and right-censored data가 될 수도 있습니다.

Left-truncated and right-censored (LTRC) data에 대한 통계적 추론은 counting process theory에 의존합니다. 강의노트의 at-risk process \(Y(t)\)에 기존의 \(Y(t)=I(t \leq T)=I(0 < t \leq T)\) 대신 \(Y(t)=I(L < t \leq T)\)를 대입하므로서 대부분의 이론 전개를 할 수 있습니다. 오늘은 LTRC 자료에 대한 R 코드 실습을 하도록 하겠습니다.

Left-truncation right-censoring DATA

R을 이용하여 지수분포에서 LTRC data를 임의로 생성하고, 이를 다시 survfit()으로 적합하여 과연 true curve (지수분포)에 가까운 모양이 얻어지는지 관찰해봅시다.

library(survival) ;
## Loading required package: splines

left-truncation time을 생성하기 전에, 먼저 보통의 right-censored data \(\{(T_i, \delta_i)\}_{i=1}^{n}\)를 만들 것이다. \((T, \delta)\)는 event time \(X\)와 censoring time \(C\)로부터 구성할 수 있습니다. 생존자료의 분석을 위한 기본적인 가정들(\(T:=\min(X,C)\), \(\delta:=I(X \leq C)\), \(X \perp C\))에서 시작합니다. 여기서는 \(n=200\), \(X_i \sim \exp(0.02)\), \(C_i \sim \exp(0.02/3)\), 그리고 각 outcome이 모두 독립이도록 자료를 생성하였습니다.

[생각해보기] \(P(X>C)\)의 값은 얼마가 될까요?

num = 200 ;
rate.event = 0.02 ;   # lambda of X_i
rate.cens = rate.event / 3 ;  # lambda of C_i. want censoring rate to be 1/4
# generate (x_i, c_i)
set.seed(100) ;     # fix the first term of pseudo-random number generator
time.event = rexp(n=num, rate=rate.event) ;  # X_i
time.cens = rexp(n=num, rate=rate.cens) ; # C_i

\(T_i:=\min(X_i,C_i)\)와 \(\delta_i:=I(X_i \leq C_i)\)를 간단(?)하게 계산하려면 함수 apply(matrix, row or column, vector-to-vector function object)와 논리 연산을 활용해봅시다.

event.obs = (time.event <= time.cens) ; # delta_i
event.obs = as.numeric(event.obs) ; # true/false -> 1/0
mean(event.obs)
## [1] 0.805

이제 left-truncation time \(\{L_i\}\)를 생성해봅시다. 구간 \((0,50)\) 위의 균등분포로부터 \(\{L_i\}\)를 뽑고, 이후 \(L_i \geq T_i\)인 자료는 모두 버리고 남은 자료를 data.trunc에 저장합시다.

temp = cbind(time.event, time.cens) ;
time.obs = apply(temp, 1, min) ;  # T_i
time.trunc = runif(n=num)*50 ; # l_i
#time.trunc = runif(n=num,min=40,max=50) ; # 추천추천추천
data.full = cbind(time.trunc, time.obs, event.obs) ; # (L, T, \delta) ;
dim(data.full)
## [1] 200   3
colnames(data.full) = c("trunc", "time", "event") ;
data.trunc = data.full[time.trunc < time.obs, ] ;
dim(data.trunc)
## [1] 118   3

이제 우리는 \(\{(L_i, T_i, \delta_i)\}_i\) 꼴의 LTRC data matrix를 얻었습니다!! 이로부터 기본적인 K-M 추정곡선과 그림을 그려봅시다. 달라지는 부분은 survival object를 만들 때이다. Surv(time=Lvec, time2=Tvec, event=deltavec, type=‘counting’)을 사용한다. at-risk process의 양 끝을 입력한다고 생각하면 매우 자연스러워 보입니다.

survobj.trunc = Surv(time=data.trunc[ ,"trunc"], time2=data.trunc[ ,"time"],
                     event=data.trunc[ ,"event"], type="counting") ;
fit.trunc = survfit(survobj.trunc ~ 1) ;
fit.trunc
## Call: survfit(formula = survobj.trunc ~ 1)
## 
## records   n.max n.start  events  median 0.95LCL 0.95UCL 
##   118.0    57.0     0.0    92.0    25.2    19.0    40.6
#summary(fit.trunc) Hidden
plot(fit.trunc)
curve( eval( exp(-rate.event*x) ), 0, 200, add=T, col="red") ;

plot of chunk unnamed-chunk-5 [생각해보기] summary 함수를 이용해서 survival 추정 값이 어떻게 계산되었는지 생각해봅시다. (n.risk가 계속 바뀌는 것을 확인해보세요.)

[생각해보기] 그래프를 그린 결과를 볼 때, true 곡선과 큰 차이가 나는가? 실제 데이터도 이런 경향을 따를 것인가? left-truncation time의 분포를 다른 방식으로 바꾸어봅시다 (이를테면 \({\rm Unif}(40,50)\)).

Disclaimer : 이 자료는 최영근 조교님의 2013년 가을 Survival anlysis 실습 수업자료를 바탕으로 만들어졌습니다.