模拟竞争风险数据

Question

我的目标是模拟可用于测试竞争性风险模型的数据集。我只是用survsim::crisk.sim函数尝试一个简单的例子，但是它并没有带来我期望的结果。

 require(survival)
 simulated_data <- survsim::crisk.sim(n = 100,
                                      foltime = 200,
                                      dist.ev = rep("weibull", 2),
                                      anc.ev = c(0.8, 0.9),
                                      beta0.ev = c(2, 4),
                                      anc.cens = 1,
                                      beta0.cens = 5,
                                      nsit = 2)

 model <- survreg(Surv(time, status) ~ 1 + strata(cause), data = simulated_data)

 exp(model$scale)

 ## cause=1  cause=2 
 ## 4.407839 2.576357 

我希望这些数字与beta0.ev相同。关于我可能做错什么的任何指示，或者其他如何模拟竞争风险数据的建议。

完成:我希望模拟数据中的事件发生在威布尔分布之后，而威布尔分布对于每个风险都是不同的。我希望能够在数据中指定一个阶层和一个集群。审查可以遵循Weibull或Bernouli分布。

Answer 1

要恢复指定的估计数，可以使用survreg和特定于原因的表示法。

这个例子使用您的参数，但是对于更多的病人，需要更精确的估计：

set.seed(101)
stack_data <- survsim::crisk.sim(n = 2000,
                                     foltime = 200,
                                     dist.ev = rep("weibull", 2),
                                     anc.ev = c(0.8, 0.9),
                                     beta0.ev = c(2, 4),
                                     anc.cens = 1,
                                     beta0.cens = 5,
                                     nsit = 2)

m1 <- survreg(Surv(time, cause==1) ~ 1, data =stack_data, dist = "weibull")
m2 <- survreg(Surv(time, cause==2) ~ 1, data = stack_data, dist = "weibull")

m1$coefficients这将接近beta0.ev的原因1

m2$coefficients这将接近beta0.ev的原因2

> m1$coefficients
(Intercept) 
   1.976449 
> m2$coefficients
(Intercept) 
   3.995716 

m1$scale这将接近anc.ev的原因1

m2$scale这将接近anc.ev的原因2

> m1$scale
[1] 0.8088574
> m2$scale
[1] 0.8923334

不幸的是，这只适用于统一审查，或低非均匀审查(如在您的例子)。

如果我们增加了审查的风险，那么拦截器就不代表beta0.ev参数

set.seed(101)
stack_data <- survsim::crisk.sim(n = 2000,
                                     foltime = 200,
                                     dist.ev = rep("weibull", 2),
                                     anc.ev = c(0.8, 0.9),
                                     beta0.ev = c(2, 4),
                                     anc.cens = 1,
                                     beta0.cens = 2, #reduced from 5, increasing the hazard function for censoring rate
                                     nsit = 2)

m1 <- survreg(Surv(time, cause==1) ~ 1, data =stack_data, dist = "weibull")
m2 <- survreg(Surv(time, cause==2) ~ 1, data = stack_data, dist = "weibull")

> m1$coefficients
(Intercept) 
   1.531818 
> m2$coefficients
(Intercept) 
   3.553687 
> 
> m1$scale
[1] 0.8139497
> m2$scale
[1] 0.8910465