我想从包含时间相关协变量的Cox比例风险模型生成生存时间。该模型是
h(t|Xi) =h_0(t) exp(gamma*Xi + alpha*mi(t))
其中
Xi是从Binomial(1,0.5)生成的，而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)
谁能帮我生成时变协变量的生存数据。
最佳答案
通过检查协变量发生变化时的观察值并在时间0或检查器时将其重新输入到队列中，将时间相关协变量输入Cox模型。因此，将基于间隔生成协变量矩阵，并根据事件观察的前/后周期将非事件多对一合并为事件合并多对二合并。您可以在事件后删除协变量修改。 Cox模型无法处理协变量值和失败的并发变化，因此我们必须排除这种可能性。
为了模拟生存结果，您可以生成协变量值和转换点。然后根据基线协变量值模拟生存结果。如果第一个协变量更改时间超过故障时间，则保留故障时间并且该参与者没有协变量变化，否则检查该故障时间的观察值，并在检查时使用新的协变量值将其重新输入到同类群组中。模拟第二个协变量值更改（如果存在）和第二个可能的故障时间，进行迭代评估，直到故障时间早于协变量更改值时才结束。
对此进行说明，某人可能会提供比我更有效的代码，但是一种简单的方法是递归。我暂时会假设基线危险函数（指数生存）是恒定的，但是根据任意基线危险函数模拟生存结果的原理已在其他地方介绍过，我留给您。为了简单起见，我还将假定m为二进制。同样，这是您总结的基础。
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)
)
}
可以运行以下示例：
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)
提供以下输出：
> 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的指数存活结果的对数危险比进行比较。