2014-07-23 5 views
3

나는 시간 의존 공변량을 포함하는 Cox proportional hazards 모델로부터 생존 시간을 생성하려고합니다. 이 모델은 XiBinomial(1,0.5)에서 생성 mi(t)time-dependent covariate이다시간 의존 공변량으로 생존 데이터를 생성하는 방법 R

h(t|Xi) =h_0(t) exp(gamma*Xi + alpha*mi(t)) 

이다. 누군가가 시간에 따라 변화하는 공변량과 생존 데이터를 생성하는 데 도움 주실 래요

#For time independent case 
# h_0(t) = 1 
gamma <- -1 
u <- runif(n=100,min=0,max=1) 
Xi <- rbinom(n=100,size=1,prob=0.5) 
T <- -log(u)/exp(gamma*Xi) 

을 다음과 같이

시간에 독립적 공변량을 위해 내가 생성.

+0

답장을 보내지 않았지만 CoxFlexBoost 패키지의 genSurv 또는 rSurvTime 패키지에서 genTDCM 기능을 사용해 볼 수 있습니다. 또한 유사한 질문이있는 게시물이 있습니다. http://stats.stackexchange.com/questions/109237/how-to-generate-survival-data-with-time-dependent-covariates-using-r –

답변

0

시간 의존 공변량은 공변량이 변경 될 때 관찰을 검열하고 시간 0 또는 검열 시간에 코호트로 다시 입력하여 Cox 모델에 입력됩니다. 따라서 공변량 행렬은 간격으로 생성되고 이벤트 관찰 전/후 기간을 기반으로하는 이벤트에 대해 다기관과 다 대다가 병합됩니다. 이벤트 다음에 공변량 수정을 삭제할 수 있습니다. 공변량 값과 실패의 동시 변경은 Cox 모델에 의해 처리되지 않으므로이 가능성을 배제해야합니다.

생존 결과를 시뮬레이션하기 위해 공변량 값과 전환점을 생성합니다. 그런 다음 기준 공변량 값에 따라 생존 결과를 시뮬레이션합니다. 첫 번째 공변량 변화 시간이 실패 시간을 초과하면, 실패 시간을 유지하고 그 참가자는 공변량 변화가 없으며, 그렇지 않으면 그 실패 시간에 관찰을 검열하고 새로운 공변량 값으로 검열 시간에 코호트에 다시 입력하십시오. 두 번째 공변산 값 변경 (존재하는 경우)과 두 번째 가능한 실패 시간을 시뮬레이션하고, 반복적으로 평가하고 실패 시간이 공변량 변경 값보다 먼저 나오면 종료하십시오.

이 코드를 작성하면 누군가가 나보다 효율적인 코드를 제공 할 수 있지만,이를 수행하는 간단한 방법은 재귀를 사용하는 것입니다. 당분간은 기저선 위험 함수 (지수 생존)가 있다고 가정 할 것입니다. 그러나 임의의 기저선 위험 함수에 따라 생존 결과를 시뮬레이션하는 원칙은 다른 곳에서 설명되어 있으며 저는이를 여러분에게 남겨 둡니다. 나는 또한 간단하게 m이 2 진이라고 가정 할 것이다. 다시 말하지만, 이것은 당신이 마무리하는 기초입니다. 다음 출력주기

set.seed(123) 
times <- sim(id=1:1000, t0= rep(0, 1000), m0=rep(0, 1000), bfail=c(-1, -2), rchange=0.4) 
times <- as.data.frame(do.call(rbind, times)) 
coxph(Surv(start, stop, event) ~ m, data=times) 

:

sim <- function(id=1:100, t0= rep(0, 100), m0=rep(0, 100), bfail=c(0,0), rchange=1) { 
    tfail <- rexp(length(id), exp(bfail[1] + bfail[2]*m0)) 
    tchange <- rexp(length(id), rchange) 
    tevent <- pmin(tfail, tchange) 
    fevent <- tfail == tevent 
    if (all(fevent)) 
    return(list(cbind('start'=t0, 'stop'=t0+tevent, 'event'=fevent, 'id'=id, 'm'=m0))) 
    c(
    list(cbind('start'=t0, 'stop'=t0+tevent, 'event'=fevent, 'id'=id, 'm'=m0)), 
    sim(id = id[!fevent], t0=(t0+tevent)[!fevent], m0=1-m0[!fevent], bfail, rchange) 
) 

} 

함께 예를 실행할 수

> coxph(Surv(start, stop, event) ~ m, data=times) 
Call: 
coxph(formula = Surv(start, stop, event) ~ m, data = times) 

    coef exp(coef) se(coef)  z  p 
m -1.917  0.147 0.100 -19.1 <2e-16 

Likelihood ratio test=533 on 1 df, p=0 
n= 2852, number of events= 1000 

는 기하 급수적 생존 결과에 대한 로그 위험 비율 -2 -1.917 계수를 비교.

관련 문제