我想在 LD50 中找到具有置信区间的致命剂量( R )。 Minitab、SPSS、SAS 等其他软件提供了这种置信区间的三种不同版本。我在 R 的任何包中都找不到这样的间隔(我也使用了 findFn 包中的 sos 函数)。
我怎样才能找到这样的间隔？我基于 Delta 方法编码了一种类型的间隔(因为不确定它的正确性)，但想使用 R 包中的任何已建立的函数。谢谢
MWE:
dose <- c(10.2, 7.7, 5.1, 3.8, 2.6, 0)
total <- c(50, 49, 46, 48, 50, 49)
affected <- c(44, 42, 24, 16, 6, 0)
finney71 <- data.frame(dose, total, affected)
fm1 <- glm(cbind(affected, total-affected) ~ log(dose),
family=binomial(link = logit), data=finney71[finney71$dose != 0, ])
summary(fm1)$coef
Estimate Std. Error z value Pr(>|z|)
(Intercept) -4.886912 0.6429272 -7.601035 2.937717e-14
log(dose) 3.103545 0.3877178 8.004650 1.198070e-15
library(MASS)
xp <- dose.p(fm1, p=c(0.50, 0.90, 0.95)) # from MASS
xp.ci <- xp + attr(xp, "SE") %*% matrix(qnorm(1 - 0.05/2)*c(-1,1), nrow=1)
zp.est <- exp(cbind(xp, attr(xp, "SE"), xp.ci[,1], xp.ci[,2]))
dimnames(zp.est)[[2]] <- c("LD", "SE", "LCL","UCL")
zp.est
LD SE LCL UCL
p = 0.50: 4.828918 1.053044 4.363708 5.343724
p = 0.90: 9.802082 1.104050 8.073495 11.900771
p = 0.95: 12.470382 1.133880 9.748334 15.952512
最佳答案
从包 drc 中，您可以获得 ED50(相同的计算)以及置信区间。
library(drc) # Directly borrowed from the drc manual
mod <- drm(affected/total ~ dose, weights = total,
data = finney71[finney71$dose != 0, ], fct = LL2.2(), type = "binomial")
#intervals on log scale
ED(mod, c(50, 90, 95), interval = "fls", reference = "control")
Estimated effective doses
(Back-transformed from log scale-based confidence interval(s))
Estimate Lower Upper
1:50 4.8289 4.3637 5.3437
1:90 9.8021 8.0735 11.9008
1:95 12.4704 9.7483 15.9525
与手动输出匹配。
“finney71”数据包含在此包中，您对 置信区间的计算与 完全匹配
drc 人员给出的示例，直到“# from MASS”注释。您应该赞扬他们，而不是声称您编写了代码。还有其他一些方法可以解决这个问题。一种是使用参数 bootstrap ，它可以通过
boot 包方便地获得。首先，我们将重新拟合模型。
library(boot)
finney71 <- finney71[finney71$dose != 0,] # pre-clean data
fm1 <- glm(cbind(affected, total-affected) ~ log(dose),
family=binomial(link = logit),
data=finney71)
为了说明，我们可以找出 LD50 和 LD75。
statfun <- function(dat, ind) {
mod <- update(fm1, data = dat[ind,])
coefs <- coef(mod)
c(exp(-coefs[1]/coefs[2]),
exp((log(0.75/0.25) - coefs[2])/coefs[1]))
}
boot_out <- boot(data = finney71, statistic = statfun, R = 1000)
boot.ci 函数可以使用这个对象为我们计算出各种置信区间。boot.ci(boot_out, index = 1, type = c('basic', 'perc', 'norm'))
##BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
##Based on 999 bootstrap replicates
##
##CALL :
##boot.ci(boot.out = boot_out, type = c("basic", "perc", "norm"),
## index = 1)
##Intervals :
##Level Normal Basic Percentile
##95% ( 3.976, 5.764 ) ( 4.593, 5.051 ) ( 4.607, 5.065 )
使用正态近似的置信区间被一些极值甩掉了很多，基本和基于百分位的区间对这些极值更加稳健。
一件有趣的事情要注意:如果斜率的符号不够清晰，我们可以得到一些相当极端的值(在 this answer 中模拟，并在 this blog post 中由 Andrew Gelman 进行了更彻底的讨论)。
set.seed(1)
x <- rnorm(100)
z = 0.05 + 0.1*x*rnorm(100, 0, 0.05) # small slope and more noise
pr = 1/(1+exp(-z))
y = rbinom(1000, 1, pr)
sim_dat <- data.frame(x, y)
sim_mod <- glm(y ~ x, data = sim_dat, family = 'binomial')
statfun <- function(dat, ind) {
mod <- update(sim_mod, data = dat[ind,])
-coef(mod)[1]/coef(mod)[2]
}
sim_boot <- boot(data = sim_dat, statistic = statfun, R = 1000)
hist(sim_boot$t[,1], breaks = 100,
main = "Bootstrap of simulated model")
上面的 delta 方法为我们提供了均值 = 6.448，较低的 ci = -36.22 和较高的 ci = 49.12，并且所有引导 CI 都为我们提供了类似的极端估计。
##Level Normal Basic Percentile
##95% (-232.19, 247.76 ) ( -20.17, 45.13 ) ( -32.23, 33.06 )
关于r - R 中 Logistic 回归的致死剂量 (LD) 置信区间，我们在Stack Overflow上找到一个类似的问题：https://stackoverflow.com/questions/28053542/