从包装中drc http://cran.r-project.org/web/packages/drc/index.html，你可以得到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”数据，以及您的置信区间计算exactly与给出的示例匹配drc伙计们，请注意“# from MASS”评论。您应该赞扬他们，而不是声称您编写了代码。
还有其他一些方法可以解决这个问题。一种是使用参数引导程序，可以通过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)
The 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 )
使用正态近似的置信区间会受到一些极值的影响，而基本区间和基于百分位数的区间对此更加稳健。
需要注意的一件有趣的事情是：如果斜率的符号足够不清楚，我们可以获得一些相当极端的值（模拟为这个答案 https://stats.stackexchange.com/questions/46523/how-to-simulate-artificial-data-for-logistic-regression，并在这篇博文 http://andrewgelman.com/2011/06/21/inference_for_a/作者：安德鲁·格尔曼）。
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 )
