[R] Cox 분석을 위한 생존 시간 데이터 생성 (시뮬레이션 코드)

본 게시글을 Cox PH 모형 분석을 위한 특정 분포를 따르는 Survival Time 데이터를 생성하는 방법에 대한 내용이다.

참고에 Stackoveflow와 해당 게시글에서 언급된 Bender et al. (2005)를 참고하여 작성하였다.

---

데이터 생성하는 법

정규분포, 이항분포, 포아송분포 등의 분포를 따르는 시뮬레이션 데이터는 R의 기본 내장함수를 통해 쉽게 사용할 수 있다.

하지만 어떠한 Cox 모형의 확장된 분석방법이 시뮬레이션 데이터에서 잘 작동하는지 확인하고자 할 때,

Cox 모형은 관찰된 생존 시간에 대한 어떠한 가정을 하는 것이 아닌, Hazard Function에 대한 가정을 하고 있기 때문에 어떻게 생성할지 막막해진다.

 

2005년도에 Stat Med에서 Cox 모형을 위한 생존시간 데이터를 생성하는 법에 대한 내용이 발표되었다,

논문의 Table 2. 는 Baseline Hazard에 대한 모수적인 가정을 하였을 때,

공변량 x에 따라 크기게 비례하는 형태를 갖는 Hazard Function과 그 Hazard Function으로부터 유도한 Survival Time T에 대한 식을 나타낸다. 이때 U는 Uniform 분포를 따른다.

 

현재 아래 Table엔 Baseline Hazard가 각각 $\lambda$, $\lambda vt^{v-1}$,... 이렇게 모수 가정하여 표현한 것이고,

원래 Cox 비례모형은 Baseline Hazard에 대한 모수적인 가정을 하지 않는다.

 

 

일반적으로 Cox 모형에서의 관심사는 $\beta$이고, $\lambda_{0}(t)$가 비모수적 추론되는 부분이다.

해당 값에 어떠한 가정을 하지 않으면 데이터 자체를 생성할 수 없기에, 각각의 분포 가정을 통해 아래 표처럼 작성했다고 생각하면 이해가 쉽다.

 

Bender et al. (2005)

---

코드

위의 표에서 언급된 수식과 일치시켜서 작성된 코드는 다음과 같다.

nn=500; p=3;
beta= c(-1, 1, 2); # interest
lambda <- v <- 1

# covariate
x <- matrix(rnorm(nn*p), ncol=p, nrow=nn)

U <- runif(n=nn)
T <- (- log(U) / (lambda * exp(x %*% beta)))^(1 / v) # Survival Time
C <- runif(nn, 0, 5) # Time for Censoring

time <- pmin(T, C) # observed time
status <- as.numeric(T <= C)

 

해당 데이터를 이용하여 Cox 모형의 분석 결과, 사전에 사용하였던 -1, 1, 2 값을 그대로 얻을 수 있었다.

생성한 데이터 확인

 

---

참고

• Bender, R., Augustin, T., & Blettner, M. (2005). Generating survival times to simulate Cox proportional hazards models. Statistics in medicine, 24(11), 1713-1723.
• Stack Overflow: https://stats.stackexchange.com/questions/135124/how-to-create-a-toy-survival-time-to-event-data-with-right-censoring