模拟竞争风险数据答案

【问题标题】：Simulate competing risk data模拟竞争风险数据
【发布时间】：2015-12-27 16:04:12
【问题描述】：

我的目标是模拟一个可用于测试竞争风险的数据集模型。我只是在尝试使用 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 的示例并确保模型 survreg 拟合的参数化与 crisk.sim 相同。
survreg 的帮助页面特别警告我们有关 Weibull 参数的参数化。 # survreg's scale = 1/(rweibull shape) # survreg's intercept = log(rweibull scale)

标签： r simulation survival-analysis weibull

【解决方案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 原因 2 将接近 beta0.ev

> 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

【讨论】：