r - 模拟竞争风险数据

我的目标是模拟可用于测试竞争风险的数据集
模型。我只是尝试使用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分布发生。我希望能够在数据中指定分层和聚类。审查可以遵循Weibull或Bernouli分布。

最佳答案

要恢复指定的估计值，可以将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原因1将接近beta0.ev
m2$coefficients原因2将接近beta0.ev

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

m1$scale原因1将接近beta0.ev
m2$scale原因1将接近beta0.ev

> 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

您可以使用特定原因的Cox回归来恢复Beta参数。
在这里，我们添加了一个对原因1没有影响的二进制协变量，对原因2没有对数的logHR为0.69（HR为2）（请注意使用负beta参数，因为survsim包使用加速的故障时间参数化，而）

set.seed(105)
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,
                                 beta = list(c(0,-0.69)),
                                  x = list(c("bern",0.3)),
                                     nsit = 2)

cox_cause1 <- coxph(Surv(time, cause==1)~factor(x), data = stack_data)
cox_cause2 <- coxph(Surv(time, cause==2)~factor(x), data = stack_data)

cox_cause1$coefficients
cox_cause2$coefficients

 factor(x)1
-0.02391504 (crisk.sim beta = 0)

factor(x)1
 0.7273799 (crisk.sim beta = -0.69)