我想知道如何模拟包含 I 型右删失观察的 n 个 Weibull 分布寿命的样本。例如,让 n = 3,shape = 3,scale = 1,审查率 = .15,审查时间 = .88。我知道如何生成 Weibull 样本,但我不知道如何生成 R 中类型 I 右删失的删失数据。
T = rweibull(3, shape=.5, scale=1)
我想知道如何模拟包含 I 型右删失观察的 n 个 Weibull 分布寿命的样本。例如,让 n = 3,shape = 3,scale = 1,审查率 = .15,审查时间 = .88。我知道如何生成 Weibull 样本,但我不知道如何生成 R 中类型 I 右删失的删失数据。
T = rweibull(3, shape=.5, scale=1)
(作为 R
编码风格问题,最好不要 T
用作变量名,因为它是 的别名 TRUE
,那样做难免会出问题。)
你的问题有点模棱两可;有几种解释它的方法。让我们来看看它们:
您规定要模拟类型 1 审查。这通常被认为意味着实验运行了一段时间,并且到那时任何没有发生过事件的研究单位都会被审查。如果这就是您的意思,那么(必然)不可能同时规定形状和比例参数,以及审查时间和速率。已经规定了任何三个,最后一个必然是固定的。
(尝试)求解形状参数:
失败;无论形状参数是什么,似乎不可能在 0.88 的审查时间与韦布尔分布的尺度参数保持为 1 的情况下获得 15% 的审查率。
optim(.5, fn=function(shp){(pweibull(.88, shape=shp, scale=1, lower.tail=F)-.15)^2})
# $par
# [1] 4.768372e-08
# ...
# There were 46 warnings (use warnings() to see them)
pweibull(.88, shape=4.768372e-08, scale=1, lower.tail=F)
# [1] 0.3678794
optim(.5, fn=function(shp){(pweibull(.88, shape=shp, scale=1, lower.tail=F)-.15)^2},
control=list(reltol=1e-16))
# $par
# [1] 9.769963e-16
# ...
# There were 50 or more warnings (use warnings() to see the first 50)
pweibull(.88, shape=9.769963e-16, scale=1, lower.tail=F)
# [1] 0.3678794
求解尺度参数:
optim(1, fn=function(scl){(pweibull(.88, shape=.5, scale=scl, lower.tail=F)-.15)^2})
# $par
# [1] 0.2445312
# ...
pweibull(.88, shape=.5, scale=0.2445312, lower.tail=F)
# [1] 0.1500135
求解审查时间:
qweibull(.15, shape=.5, scale=1, lower.tail=F)
# [1] 3.599064
求解审查率:
pweibull(.88, shape=.5, scale=1, lower.tail=F)
# [1] 0.3913773
另一方面,我们可以将审查视为在整个研究过程中随机(通常是独立地)发生的,例如由于辍学。在这种情况下,程序是模拟两组 Weibull 变量。然后您只需注意哪个先出现:您使用较小的值作为端点,如果较小的值是审查时间,则将该单位称为审查单位。例如:
set.seed(0775)
t = rweibull(3, shape=.5, scale=1)
t # [1] 0.7433678 1.1325749 0.2784812
c = rweibull(3, shape=.5, scale=1.5)
c # [1] 3.3242417 2.8866217 0.9779436
time = pmin(t, c)
time # [1] 0.7433678 1.1325749 0.2784812
cens = ifelse(c<t, 1, 0)
cens # [1] 0 0 0
只是为了确保我们谈论的是同一件事,I 型审查是什么时候
...一个实验有一定数量的受试者或项目,并在预定时间停止实验,此时剩余的任何受试者都将被右删失。
要使用censoring time = 0.88生成正确的删失数据,您只需使用以下min
函数:
T <- rweibull(3, shape=.5, scale=1)
censoring_time <- 0.88
T_censored <- min(censoring_time, T)
但是,我不完全确定您所说的“审查率 = 0.15 ”是什么意思……您的意思是说您的科目中有 15% 是正确审查的吗?这些关于审查的说明似乎表明I 型审查所需的唯一参数是审查时间,所以我不确定这个比率是如何影响的。