Confidence Intervals for Lethal Dose (LD) for Logistic Regression in R

Question

I want to find Lethal Dose (LD50) with its confidence interval in R. Other softwares line Minitab, SPSS, SAS provide three different versions of such confidence intervals. I could not find such intervals in any package in R (I also used findFn function from sos package).

How can I find such intervals? I coded for one type of intervals based on Delta method (as not sure about it correctness) but would like to use any established function from R package. Thanks

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

Answer 1

From the package drc, you can get the ED50 (same calculation), along with confidence intervals.

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

Which matches the manual output.

The "finney71" data is included in this package, and your calculation of confidence intervals exactly matches the example given by the drc folks, down to the "# from MASS" comment. You should give credit to them, rather than claiming you wrote the code.

---

There's a few other ways to figure this out. One is using parametric bootstrap, which is conveniently available through the boot package.

First, we'll refit the model.

library(boot)

finney71 <- finney71[finney71$dose != 0,] # pre-clean data 
fm1 <- glm(cbind(affected, total-affected) ~ log(dose),
                 family=binomial(link = logit), 
                 data=finney71)

And for illustration, we can figure out the LD50 and 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 function can work out a variety of confidence intervals for us, using this object.

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 )  

The confidence intervals using the normal approximation are thrown off quite a bit by a few extreme values, which the basic and percentile-based intervals are more robust to.

One interesting thing to note: if the sign of the slope is sufficiently unclear, we can get some rather extreme values (simulated as in this answer, and discussed more thoroughly in this blog post by 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")

The delta method above gives us mean = 6.448, lower ci = -36.22, and upper ci = 49.12, and all of the bootstrap CIs give us similarly extreme estimates.

##Level      Normal              Basic              Percentile     
##95%   (-232.19,  247.76 )   ( -20.17,   45.13 )   ( -32.23,   33.06 )