# Grouped mortality: number killed out of n tested at each dose
dose   <- c(1.6907, 1.7242, 1.7552, 1.7842, 1.8113, 1.8369, 1.8610, 1.8839)
n      <- c(59, 60, 62, 56, 63, 59, 62, 60)
killed <- c(6, 13, 18, 28, 52, 53, 61, 60)

model <- glm(cbind(killed, n - killed) ~ dose,
             family = binomial(link = "probit"))

sm <- summary(model)$coefficients   # estimate, std. error, z value
wald <- sm[, "z value"]^2           # Wald chi-square on 1 df
cbind(sm, Wald = wald,
      p_Wald = pchisq(wald, df = 1, lower.tail = FALSE))

lev <- c(0.50, 0.90)                      # LD50 and LD90 (any level in 0-1)
ed  <- MASS::dose.p(model, p = lev)       # estimate on the fitted dose scale
est <- as.numeric(ed)
se  <- as.numeric(attr(ed, "SE"))         # delta-method standard error

z <- qnorm(0.975)                         # 95% CI (alpha = 0.05), delta method
data.frame(level    = paste0("LD", lev * 100),
           estimate = est,
           SE       = se,
           lower    = est - z * se,
           upper    = est + z * se)

# If a log10 dose scale is used, back-transform the estimate and interval
# (est_orig = 10^est) and report the SE on the same scale via the delta method:
#   se_orig = se * 10^est * log(10)

# Likelihood-ratio test of the overall dose effect
lr    <- model$null.deviance - model$deviance
lr_df <- model$df.null - model$df.residual
pchisq(lr, lr_df, lower.tail = FALSE)

# Pseudo R-squared at the individual (Bernoulli) level.
# N = total subjects; the null log-likelihood uses the overall response rate.
N   <- sum(n); r <- sum(killed); p0 <- r / N
ll0 <- r * log(p0) + (N - r) * log(1 - p0)   # individual-level null log-likelihood
R2_CS <- 1 - exp(-lr / N)                     # Cox & Snell
c(CoxSnell   = R2_CS,
  Nagelkerke = R2_CS / (1 - exp(2 * ll0 / N)))

# Valid lack-of-fit test for GROUPED count data only.
pearson  <- sum(residuals(model, type = "pearson")^2)
deviance <- sum(residuals(model, type = "deviance")^2)
df       <- model$df.residual
data.frame(statistic = c("Pearson", "Deviance"),
           chisq = c(pearson, deviance),
           df    = df,
           p     = pchisq(c(pearson, deviance), df, lower.tail = FALSE))
# A large p-value indicates no evidence of lack of fit / overdispersion.

xs <- seq(min(dose), max(dose), length.out = 200)
pr <- predict(model, newdata = data.frame(dose = xs),
              type = "link", se.fit = TRUE)
z  <- qnorm(0.975)                        # 95% band
curve <- data.frame(dose = xs,
                    fitted = pnorm(pr$fit),
                    lower  = pnorm(pr$fit - z * pr$se.fit),
                    upper  = pnorm(pr$fit + z * pr$se.fit))
# Observed proportions killed/n are overlaid at each dose.

dat <- data.frame(dose = dose, dead = killed, alive = n - killed)

set.seed(42)                              # reproducible resampling
n_boot <- 1000
boot_ld50 <- replicate(n_boot, {
  idx <- sample(nrow(dat), replace = TRUE)               # case resampling
  bf  <- glm(cbind(dead, alive) ~ dose,
             family = binomial(link = "probit"), data = dat[idx, ])
  as.numeric(MASS::dose.p(bf, p = 0.5))
})
quantile(boot_ld50, c(0.025, 0.975), names = FALSE)      # percentile interval