如何生成，使用R

我想從包含隨時間變化的協變量Cox比例風險模型產生的存活時間與時間相關的協變量生存數據。該模型是如何生成，使用R

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

其中從Binomial(1,0.5)產生Xi和mi(t)是time-dependent covariate。

對於時間無關的協我產生如下

#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)

任何人都可以請幫我隨時間變化的協變量產生生存數據。

來源

2014-07-23 Sheikh

晚回覆，但你可以嘗試從CoxFlexBoost包封裝genSurv或rSurvTime功能genTDCM。此外，還有這個帖子有類似的問題http://stats.stackexchange.com/questions/109237/how-to-generate-survival-data-with-time-dependent-covariates-using-r –

時間相關協變量通過在協變量在時間0或在檢查時重新進入隊列時檢查觀察值而被輸入到Cox模型中。因此，協變量矩陣是根據事件觀察的前/後期以間隔生成的，並且針對非事件和針對事件的多對一進行合併。您可以在事件後刪除協變量修改。協變量值和失敗的併發變化不是由Cox模型處理的，所以我們必須排除這種可能性。

對於模擬生存結果，可以生成協變量值和切換點。然後根據基線協變量值模擬生存結果。如果第一個協變化的時間超過了故障時間，保持故障時間和參與者沒有改變協，否則審查那個故障時間的觀察，他們在終檢的時間與新協值重新進入隊列。模擬第二個協變量值的變化（如果存在的話）和第二個可能的失效時間，迭代評估它們，只有在失效時間先於協變量變化值時才結束。

鋪設了這一點，有人可以提供更有效的代碼比我，但一個簡單的方法來做到這一點是與遞歸。我會暫時假設存在一個持續的基準危險函數（指數生存），但根據任意基線危險函數模擬生存結果的原理已在其他地方描述，我將其留給您。爲了簡單起見，我還會假設m是二進制的。再次，這是你打包的基礎。

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

比較-1.917係數爲-2爲指數存活結果日誌危害比。

來源

2018-01-18 16:31:58 AdamO

如何生成，使用R

回答

相關問題