Dose-response curves and the LC50

R
ecotoxicology
GLM
ecology tutorial
ggplot2
Fit a dose-response curve in R with three binomial links and measure what changes: the LC50 barely moves, the LC10 moves, and AIC cannot tell the links apart.
Author

Tidy Ecology

Published

2026-07-22

A laboratory runs a 48 hour acute test on Daphnia magna. Seven concentrations in a doubling series, twenty-five animals at each, and at the end somebody counts how many animals are immobilised. The count sheet has seven rows on it. What leaves the building is a single number, the LC50, and increasingly a second one: a low-effect concentration such as the LC10, which is what a risk assessment divides by an assessment factor to get a threshold. Neither number is on the count sheet. Both are read off a curve that was fitted to it, at a place on that curve chosen by a regulation.

That is the whole subject of this tutorial. An effective dose is a quantile of a fitted curve, so it inherits everything the curve family assumes, and it inherits more of it the further from the middle of the data you read. This post fits the same seven counts with a probit, a logit and a complementary log-log link and measures how far the answer moves. The short version of the measurements below: at the LC50 the two symmetric links agree to within a ratio of 1.0022 and the asymmetric one does not agree at all; at the LC1 the three answers span more than a factor of two; and over 300 simulated experiments from each of three generating links AIC recovers the right link in 0.5778 of them, picking the wrong symmetric link more often than the right one.

The fitting machinery is the binomial GLM of Logistic regression for presence-absence data, and nothing here changes it. What changes is the question put to the fit. There, the object of interest is a coefficient, an estimate of how strongly a predictor moves the probability. Here the coefficients are a nuisance and the object of interest is a ratio of them, the concentration at which the curve crosses a stated level. Ratios of estimates are worse behaved than the estimates they are built from, which is why this post keeps returning to the same measurement. The log-logistic distribution also appears in Parametric survival and the AFT model, where it runs along a time axis and the quantity read off it is a survival time; here the same algebra runs along a concentration axis at a fixed observation time, and the quantity read off it is a dose. The R package written for this work is drc. Everything below uses glm instead, so that each step is visible.

library(ggplot2)

te_pal <- list(forest = "#275139", green = "#2f8f63", sage = "#93a87f",
               clay = "#b5534e", gold = "#cda23f", line = "#dad9ca",
               ink = "#16241d", paper = "#f5f4ee")

theme_te <- function() {
  theme_minimal(base_size = 12) +
    theme(panel.grid.minor = element_blank(),
          panel.grid.major = element_line(colour = "#e7e6dc"),
          plot.background = element_rect(fill = "#f5f4ee", colour = NA),
          panel.background = element_rect(fill = "#f5f4ee", colour = NA),
          plot.title = element_text(face = "bold", colour = te_pal$ink),
          axis.title = element_text(colour = "#2c3a31"))
}

A doubling series, and three curves that agree at the LC50

The design is the standard one: concentrations of 0.8, 1.6, 3.2, 6.4, 12.8, 25.6 and 51.2 micrograms per litre, twenty-five animals at each, 175 animals in total, scored at 48 hours. A range-finding test has already put the LC50 somewhere near 10, which is why the series is centred where it is.

The generating curve needs a link, and any choice of link is a choice about the shape of the tails, so all three are set up here on equal terms. Each truth is a straight line on log concentration through its own link, each passes through a response of one half at exactly 9 micrograms per litre, and each has the same steepness there: a slope of 0.65 in response per natural log unit of concentration, at the point where the response is one half. The three generating curves are therefore indistinguishable at the LC50 and as close as three different functional forms can be made in the middle of the design.

conc  <- c(0.8, 1.6, 3.2, 6.4, 12.8, 25.6, 51.2)
n_per <- 25
lc50_true  <- 9
slope_mid  <- 0.65
links <- c("logit", "probit", "cloglog")
lev   <- c(0.5, 0.2, 0.1, 0.05, 0.01)

link_q <- function(p, link) switch(link, logit = log(p / (1 - p)),
                                   probit = qnorm(p), cloglog = log(-log(1 - p)))
link_p <- function(eta, link) switch(link, logit = plogis(eta),
                                     probit = pnorm(eta), cloglog = 1 - exp(-exp(eta)))
dens_half <- function(link) switch(link, logit = 0.25, probit = dnorm(0),
                                   cloglog = log(2) / 2)
slope_of <- function(link) slope_mid / dens_half(link)
true_p  <- function(x, link)
  link_p(link_q(0.5, link) + slope_of(link) * (log(x) - log(lc50_true)), link)
true_lc <- function(p, link)
  exp(log(lc50_true) + (link_q(p, link) - link_q(0.5, link)) / slope_of(link))

c(concentrations = length(conc), animals_each = n_per,
  animals_total = length(conc) * n_per, test_duration_hours = 48,
  true_lc50 = lc50_true, slope_at_lc50 = slope_mid)
     concentrations        animals_each       animals_total test_duration_hours 
               7.00               25.00              175.00               48.00 
          true_lc50       slope_at_lc50 
               9.00                0.65 
print(c(response_levels_percent = 100 * lev))
response_levels_percent1 response_levels_percent2 response_levels_percent3 
                      50                       20                       10 
response_levels_percent4 response_levels_percent5 
                       5                        1 
print(round(rbind(concentration = conc, true_response = true_p(conc, "logit")), 4))
                [,1]   [,2]   [,3]   [,4]    [,5]    [,6]    [,7]
concentration 0.8000 1.6000 3.2000 6.4000 12.8000 25.6000 51.2000
true_response 0.0018 0.0111 0.0636 0.2919  0.7142  0.9381  0.9892
print(round(sapply(links, slope_of), 4))
  logit  probit cloglog 
 2.6000  1.6293  1.8755 
truth_tab <- sapply(links, function(l) sapply(lev, function(p) true_lc(p, l)))
rownames(truth_tab) <- paste0("LC", 100 * lev)
print(round(cbind(truth_tab,
                  spread = apply(truth_tab, 1, function(z) max(z) / min(z))), 4))
      logit probit cloglog spread
LC50 9.0000 9.0000  9.0000 1.0000
LC20 5.2806 5.3692  4.9180 1.0917
LC10 3.8657 4.0987  3.2962 1.2435
LC5  2.9001 3.2795  2.2456 1.4604
LC1  1.5371 2.1585  0.9417 2.2922

Under the logit truth the seven concentrations produce expected responses of 0.0018, 0.0111, 0.0636, 0.2919, 0.7142, 0.9381 and 0.9892. That is a well-shaped experiment by any standard: the LC50 sits between the fourth and fifth concentration, and there is a real gradient across the whole series rather than a jump from nothing to everything.

The last table is worth a pause, because it is written before any data exist. Three curves that are identical at the LC50 and equally steep there have LC10 values of 3.8657, 4.0987 and 3.2962. By the LC1 they read 1.5371, 2.1585 and 0.9417, a spread of 2.2922 between the largest and the smallest. No estimation has happened yet. The disagreement is a property of the three functional forms, sitting there before the first animal is counted, and estimation can only add to it.

Recovering an LC50 that we already know

One dataset, drawn from the logit truth. The counts are what a technician would write down.

set.seed(20260729)
dead <- rbinom(length(conc), n_per, true_p(conc, "logit"))
print(rbind(concentration = conc, immobilised = dead,
            proportion = round(dead / n_per, 4)))
              [,1] [,2] [,3] [,4]  [,5]  [,6] [,7]
concentration  0.8 1.60  3.2 6.40 12.80 25.60 51.2
immobilised    0.0 1.00  0.0 8.00 19.00 23.00 25.0
proportion     0.0 0.04  0.0 0.32  0.76  0.92  1.0
c(immobilised_total = sum(dead), animals_total = length(conc) * n_per)
immobilised_total     animals_total 
               76               175

In total 76 of the 175 animals were immobilised. Note the shape of the low end: none at 0.8, one at 1.6, none at 3.2. That non-monotone run is ordinary sampling noise, and it is also the entire empirical content of the region where a low-effect concentration will later be read off.

Fitting is one line. The estimate of a level-p effective dose is the concentration at which the linear predictor equals the link function evaluated at p, which is a ratio of two estimated coefficients. Two intervals are worth having. The delta method gives a standard error on the log concentration scale by linearising that ratio. Fieller’s theorem does better: it inverts the test that the linear predictor equals its target value at a given concentration, which handles the ratio exactly rather than by first-order approximation, and it can return an unbounded interval when the slope is poorly determined, which is honest rather than a defect.

n_warn <- 0
fit_link <- function(y, link)
  withCallingHandlers(
    glm(cbind(y, n_per - y) ~ log(conc), family = binomial(link = link)),
    warning = function(w) { n_warn <<- n_warn + 1; invokeRestart("muffleWarning") })

lc_hat <- function(fit, p, link) {
  b <- coef(fit)
  unname(exp((link_q(p, link) - b[1]) / b[2]))
}
lc_se_log <- function(fit, p, link) {
  b <- coef(fit); V <- vcov(fit); q <- link_q(p, link)
  g <- c(-1 / b[2], -(q - b[1]) / b[2]^2)
  sqrt(as.numeric(t(g) %*% V %*% g))
}
fieller <- function(fit, p, link, conf = 0.95) {
  b <- coef(fit); V <- vcov(fit); q <- link_q(p, link)
  z <- qnorm(1 - (1 - conf) / 2)
  aa <- b[2]^2 - z^2 * V[2, 2]
  bb <- 2 * (b[2] * (b[1] - q) - z^2 * V[1, 2])
  cc <- (b[1] - q)^2 - z^2 * V[1, 1]
  disc <- bb^2 - 4 * aa * cc
  if (aa <= 0 || disc < 0) return(c(NA_real_, NA_real_))
  unname(exp(sort(c((-bb - sqrt(disc)) / (2 * aa), (-bb + sqrt(disc)) / (2 * aa)))))
}

f_logit <- fit_link(dead, "logit")
print(round(summary(f_logit)$coefficients, 4))
            Estimate Std. Error z value Pr(>|z|)
(Intercept)  -5.9063     0.9554 -6.1822        0
log(conc)     2.7143     0.4217  6.4374        0
round(c(residual_deviance = deviance(f_logit),
        residual_df = df.residual(f_logit),
        lack_of_fit_p = pchisq(deviance(f_logit), df.residual(f_logit),
                               lower.tail = FALSE),
        warnings_so_far = n_warn), 4)
residual_deviance       residual_df     lack_of_fit_p   warnings_so_far 
           5.4364            5.0000            0.3650            0.0000 
one_row <- function(p) {
  ci <- fieller(f_logit, p, "logit")
  c(estimate = lc_hat(f_logit, p, "logit"), se_log = lc_se_log(f_logit, p, "logit"),
    fieller_lower = ci[1], fieller_upper = ci[2], truth = true_lc(p, "logit"))
}
print(round(rbind(LC50 = one_row(0.5), LC10 = one_row(0.1)), 4))
     estimate se_log fieller_lower fieller_upper  truth
LC50   8.8108 0.1013        7.1575       10.8594 9.0000
LC10   3.9216 0.1607        2.5692        5.0794 3.8657

The fitted slope is 2.7143 against a true 2.6, the residual deviance is 5.4364 on 5 degrees of freedom, and the lack-of-fit probability is 0.3650. The LC50 comes out at 8.8108 against a truth of 9, with a standard error of 0.1013 on the log scale and a Fieller interval running from 7.1575 to 10.8594. The LC10 comes out at 3.9216 against a truth of 3.8657, with an interval from 2.5692 to 5.0794. Both are correct, and the second interval is already noticeably wider relative to its centre than the first.

One dataset proves nothing about an estimator, so the same experiment is run 2000 times.

set.seed(20260729)
n_sim <- 2000
cal <- t(sapply(seq_len(n_sim), function(i) {
  y <- rbinom(length(conc), n_per, true_p(conc, "logit"))
  f <- fit_link(y, "logit")
  ci <- lapply(c(0.5, 0.1, 0.01), function(p) fieller(f, p, "logit"))
  hit <- function(k, p) ci[[k]][1] < true_lc(p, "logit") &&
    true_lc(p, "logit") < ci[[k]][2]
  c(lc50 = lc_hat(f, 0.5, "logit"), lc10 = lc_hat(f, 0.1, "logit"),
    cov50 = hit(1, 0.5), cov10 = hit(2, 0.1), cov1 = hit(3, 0.01),
    wid50 = ci[[1]][2] / ci[[1]][1], wid10 = ci[[2]][2] / ci[[2]][1],
    wid1 = ci[[3]][2] / ci[[3]][1])
}))
c(datasets = n_sim, unbounded_intervals = sum(is.na(cal[, "cov50"])),
  warnings_so_far = n_warn)
           datasets unbounded_intervals     warnings_so_far 
               2000                   0                   0 
round(c(median_lc50 = median(cal[, "lc50"]), true_lc50 = lc50_true,
        median_lc10 = median(cal[, "lc10"]), true_lc10 = true_lc(0.1, "logit"),
        coverage_lc50 = mean(cal[, "cov50"]), coverage_lc10 = mean(cal[, "cov10"]),
        coverage_lc1 = mean(cal[, "cov1"]),
        median_width_ratio_lc50 = median(cal[, "wid50"]),
        median_width_ratio_lc10 = median(cal[, "wid10"]),
        median_width_ratio_lc1 = median(cal[, "wid1"])), 4)
            median_lc50               true_lc50             median_lc10 
                 9.0546                  9.0000                  3.9083 
              true_lc10           coverage_lc50           coverage_lc10 
                 3.8657                  0.9575                  0.9560 
           coverage_lc1 median_width_ratio_lc50 median_width_ratio_lc10 
                 0.9540                  1.5229                  1.9901 
 median_width_ratio_lc1 
                 3.3890

The machinery is calibrated. The median estimate over 2000 experiments is 9.0546 for an LC50 of 9 and 3.9083 for an LC10 of 3.8657. The Fieller intervals cover at 0.9575, 0.9560 and 0.9540 for the LC50, the LC10 and the LC1, which is what a nominal 0.95 procedure should do, and not one of the 2000 returned an unbounded interval on this design.

The interval widths are the interesting column. Expressed as the ratio of the upper limit to the lower, the median LC50 interval spans a factor of 1.5229, the LC10 interval 1.9901, and the LC1 interval 3.3890. Same animals, same fit, same nominal confidence: the further down the curve the question is asked, the wider the answer, and that is with the model form known to be exactly right. It is not right in any real analysis, which is the rest of this post.

Where the three curves part company

The dataset above carries sampling noise, so it cannot separate what the links do from what the draw did. The clean way to ask is to hand each link an infinite dataset from the same generating curve. That is what a misspecified maximum likelihood estimator converges to: the parameter values that minimise the Kullback-Leibler divergence from the truth, which for grouped binomial data is the expected deviance. Minimising it directly with optim gives the large-sample fit of each link with no noise in it at all.

kl_fit <- function(true_link, link) {
  pt <- true_p(conc, true_link)
  start <- c(link_q(0.5, link) - slope_of(link) * log(lc50_true), slope_of(link))
  obj <- function(b) {
    ph <- link_p(b[1] + b[2] * log(conc), link)
    ph <- pmin(pmax(ph, 1e-12), 1 - 1e-12)
    -sum(n_per * (pt * log(ph) + (1 - pt) * log(1 - ph)))
  }
  optim(start, obj, method = "BFGS", control = list(reltol = 1e-12))$par
}
asym_b <- lapply(links, function(l) kl_fit("logit", l))
names(asym_b) <- links
asym_tab <- sapply(links, function(l)
  sapply(lev, function(p) exp((link_q(p, l) - asym_b[[l]][1]) / asym_b[[l]][2])))
rownames(asym_tab) <- paste0("LC", 100 * lev)
print(round(cbind(asym_tab,
                  spread = apply(asym_tab, 1, function(z) max(z) / min(z)),
                  logit_vs_probit = asym_tab[, "logit"] / asym_tab[, "probit"],
                  truth = sapply(lev, function(p) true_lc(p, "logit"))), 4))
      logit probit cloglog spread logit_vs_probit  truth
LC50 9.0000 8.9805 10.4376 1.1623          1.0022 9.0000
LC20 5.2806 4.9716  5.0593 1.0621          1.0621 5.2806
LC10 3.8657 3.6497  3.1323 1.2342          1.0592 3.8657
LC5  2.9001 2.8275  1.9775 1.4666          1.0257 2.9001
LC1  1.5371 1.7517  0.6979 2.5098          0.8775 1.5371
print(round(rbind(true_response = true_p(conc, "logit"),
                  sapply(conc, function(x)
                    sapply(links, function(l)
                      link_p(asym_b[[l]][1] + asym_b[[l]][2] * log(x), l)))), 4))
                [,1]   [,2]   [,3]   [,4]   [,5]   [,6]   [,7]
true_response 0.0018 0.0111 0.0636 0.2919 0.7142 0.9381 0.9892
logit         0.0018 0.0111 0.0636 0.2919 0.7142 0.9381 0.9892
probit        0.0003 0.0070 0.0710 0.3148 0.6930 0.9320 0.9934
cloglog       0.0124 0.0362 0.1032 0.2756 0.6148 0.9405 0.9998

With the noise removed, the logit recovers its own truth exactly, as it must. The probit lands on 8.9805 at the LC50, 3.6497 at the LC10 and 1.7517 at the LC1, so it tracks the logit closely in the middle and drifts away from it downwards. The complementary log-log fit lands on 10.4376, 3.1323 and 0.6979, and its LC1 is less than half the logit value. The three-way spread is 1.1623 at the LC50, dips to 1.0621 at the LC20 where the complementary log-log curve crosses the other two, then climbs to 1.2342 at the LC10, 1.4666 at the LC5 and 2.5098 at the LC1.

The second table says where that comes from. At 0.8 micrograms per litre the true response is 0.0018, the large-sample complementary log-log fit predicts 0.0124 and the probit predicts 0.0003. At 3.2 the truth is 0.0636 and the complementary log-log fit says 0.1032. None of that could be caught by eye on a plot of the fitted curves, because 0.0124 and 0.0018 are both a flat line along the bottom of the panel. It is exactly the error that matters when the reported quantity is the concentration at which the response reaches one tenth.

Now count the evidence available in that region.

lc10_logit <- lc_hat(fits[["logit"]], 0.1, "logit")
below <- conc < lc10_logit
c(concentrations_below_estimated_lc10 = sum(below),
  animals_below = sum(below) * n_per, immobilised_below = sum(dead[below]),
  lowest_tested = min(conc))
concentrations_below_estimated_lc10                       animals_below 
                                3.0                                75.0 
                  immobilised_below                       lowest_tested 
                                1.0                                 0.8 
grid_lo <- exp(seq(log(0.05), log(60), length.out = 2000))
pred_m <- sapply(links, function(l)
  link_p(coef(fits[[l]])[1] + coef(fits[[l]])[2] * log(grid_lo), l))
ratio_lo <- apply(pred_m, 1, function(z) max(z) / min(z))
sep_dose <- function(thr) max(grid_lo[ratio_lo > thr])
round(c(spread_exceeds_1.5_below = sep_dose(1.5),
        spread_exceeds_2_below = sep_dose(2),
        spread_exceeds_5_below = sep_dose(5),
        logit_response_at_that_dose = unname(link_p(coef(fits[["logit"]])[1] +
          coef(fits[["logit"]])[2] * log(sep_dose(2)), "logit"))), 4)
   spread_exceeds_1.5_below      spread_exceeds_2_below 
                     3.3324                      2.4738 
     spread_exceeds_5_below logit_response_at_that_dose 
                     1.5435                      0.0308

Three of the seven concentrations sit below the estimated LC10, so 75 animals were tested there. One of them was immobilised. That single animal is the entire observational basis for the part of the curve a low-effect concentration is read from; everything else the fit knows about that region it knows by assuming a shape and continuing it downwards from where the responses were common.

The predicted responses of the three fits differ by more than a factor of 1.5 below 3.3324 micrograms per litre, by more than a factor of 2 below 2.4738, and by more than a factor of 5 below 1.5435. The lowest concentration actually tested is 0.8, so the separation does not begin outside the design at all: it begins inside it, at a concentration where the logit fit predicts a response of 0.0308, which no group of twenty-five animals can tell apart from a different curve shape.

tail_df <- do.call(rbind, lapply(links, function(l)
  data.frame(x = grid_lo, p = link_p(coef(fits[[l]])[1] +
                                       coef(fits[[l]])[2] * log(grid_lo), l), link = l)))
tail_df <- tail_df[tail_df$x >= 0.35 & tail_df$x <= 14, ]
tail_df$link <- factor(tail_df$link, levels = links)
lc10_df <- data.frame(x = sapply(links, function(l) lc_hat(fits[[l]], 0.1, l)),
                      link = factor(links, levels = links))
floor_p <- 5e-4
seen <- obs_df[obs_df$x <= 14, ]
zero_df <- seen[seen$p == 0, ]
zero_df$p <- floor_p

ggplot(tail_df, aes(x, p, colour = link)) +
  geom_hline(yintercept = 0.1, linetype = "dashed", colour = "#5d6b61") +
  geom_line(linewidth = 0.9) +
  geom_segment(data = lc10_df, aes(x = x, xend = x, y = 0.05, yend = 0.2),
               linewidth = 1.1) +
  geom_point(data = seen[seen$p > 0, ], aes(x, p), inherit.aes = FALSE,
             size = 2.6, colour = te_pal$ink) +
  geom_point(data = zero_df, aes(x, p), inherit.aes = FALSE, size = 2.6,
             shape = 6, colour = te_pal$ink) +
  scale_x_log10(breaks = c(0.4, 0.8, 1.6, 3.2, 6.4, 12.8)) +
  scale_y_log10(breaks = c(0.001, 0.003, 0.01, 0.03, 0.1, 0.3),
                labels = c("0.001", "0.003", "0.01", "0.03", "0.1", "0.3")) +
  scale_colour_manual(values = link_cols, name = NULL) +
  coord_cartesian(ylim = c(4e-4, 0.45)) +
  labs(x = "Concentration (ug per litre)", y = "Proportion immobilised",
       title = "Below a ten per cent response the three fits stop agreeing",
       subtitle = paste("Open triangles: doses at which no animal was immobilised",
                        "(zero has no place on a log axis)")) +
  theme_te() +
  theme(legend.position = "top",
        plot.subtitle = element_text(size = 9.5, colour = "#2c3a31"))
Three curves rising steeply from the bottom left on a panel with both axes logarithmic. At the right of the panel the curves lie on top of one another; to the left they fan apart, the complementary log-log curve staying well above the other two and the probit curve dropping fastest. Three short vertical ticks mark where each curve crosses a dashed line at a response of one tenth. Two open triangles sit along the floor of the panel, at the two tested doses where no animal was immobilised, and a line of subtitle text below the title says so.
Figure 2: The lower tail of the same three fitted curves, with the response on a logarithmic axis so that ratios between the curves are visible. The vertical ticks mark the LC10 of each link and the dashed line is a response of one tenth. The two concentrations at which no animal was immobilised are drawn as open triangles along the floor of the panel, since zero has no place on a log axis.

How much of the disagreement is model form and how much is the particular draw? Repeating the experiment 400 times and recording the spread of the three estimates at every level answers it.

set.seed(20260729)
n_spread <- 400
lev_grid <- sort(unique(c(exp(seq(log(0.005), log(0.5), length.out = 40)),
                          0.5, 0.1, 0.05, 0.01)))
sp <- vapply(seq_len(n_spread), function(i) {
  y <- rbinom(length(conc), n_per, true_p(conc, "logit"))
  bs <- lapply(links, function(l) coef(fit_link(y, l)))
  est <- sapply(seq_along(links), function(j)
    exp((link_q(lev_grid, links[j]) - bs[[j]][1]) / bs[[j]][2]))
  apply(est, 1, function(z) max(z) / min(z))
}, numeric(length(lev_grid)))
sp_med <- apply(sp, 1, median)

asym_curve <- sapply(links, function(l)
  exp((link_q(lev_grid, l) - asym_b[[l]][1]) / asym_b[[l]][2]))
sys_all <- apply(asym_curve, 1, function(z) max(z) / min(z))
sys_sym <- apply(asym_curve[, c("logit", "probit")], 1,
                 function(z) max(z) / min(z))

c(datasets = n_spread, warnings_so_far = n_warn)
       datasets warnings_so_far 
            400              23 
sel_lev <- c(0.5, 0.1, 0.05, 0.01)
idx <- match(sel_lev, lev_grid)
print(round(rbind(response_level = lev_grid[idx], systematic_all = sys_all[idx],
                  systematic_symmetric_pair = sys_sym[idx],
                  median_sampled_all = sp_med[idx]), 4))
                            [,1]   [,2]   [,3]   [,4]
response_level            0.5000 0.1000 0.0500 0.0100
systematic_all            1.1623 1.2342 1.4666 2.5098
systematic_symmetric_pair 1.0022 1.0592 1.0257 1.1396
median_sampled_all        1.1566 1.1573 1.3491 2.1955

The median spread over 400 fresh datasets is 1.1566 at the LC50 and 2.1955 at the LC1, against large-sample values of 1.1623 and 2.5098. Sampling noise is not what makes the links disagree. The disagreement is systematic, it is already there in an infinitely large experiment, and running more animals through the same design does not touch it.

Keeping the two symmetric links apart from the asymmetric one sharpens the picture further. Probit against logit alone spreads 1.0022 at the LC50 and 1.1396 at the LC1, so the symmetric pair never separates far, at any level. The 2.5098 at the LC1 is entirely the price of the tail shape that the complementary log-log link imposes, and that link is the only one of the three that is asymmetric about a response of one half.

lad <- rbind(
  data.frame(p = lev_grid, ratio = sys_all, what = "All three links, large sample"),
  data.frame(p = lev_grid, ratio = sys_sym, what = "Probit against logit, large sample"),
  data.frame(p = lev_grid, ratio = sp_med, what = "All three links, median of 400 datasets"))
lad$what <- factor(lad$what, levels = unique(lad$what))

ggplot(lad, aes(100 * p, ratio, colour = what, linetype = what)) +
  geom_hline(yintercept = 1, colour = te_pal$line, linewidth = 0.8) +
  geom_line(linewidth = 0.9) +
  scale_x_log10(breaks = c(0.5, 1, 2, 5, 10, 20, 50)) +
  scale_colour_manual(values = c(te_pal$clay, te_pal$forest, te_pal$gold), name = NULL) +
  scale_linetype_manual(values = c("solid", "solid", "22"), name = NULL) +
  labs(x = "Response level read off the curve (per cent)",
       y = "Largest estimate divided by smallest",
       title = "The links disagree more the further down the curve you read") +
  theme_te() +
  theme(legend.position = "top", legend.text = element_text(size = 8.5))
Three lines against response level on a logarithmic axis running from half a per cent on the left to fifty per cent on the right. At the right the lines sit between one and about one point two; moving left towards one per cent the line for all three links climbs past two and a half, the line for the symmetric pair alone stays close to one, and the dashed simulated line runs just below the line for all three links.
Figure 3: The largest of the three effective dose estimates divided by the smallest, as a function of the response level read off the curve. The two large-sample lines carry no sampling noise: they are what the links do to an infinitely large experiment on the same generating curve.

The honest limit

Two of the numbers above deserve to be put side by side. At the LC10 the sampling interval spans a factor of 1.9901 while the three links disagree by 1.2342, so on this design the animals are still the dominant source of uncertainty. At the LC1 the interval spans 3.3890 and the links disagree by 2.5098, which is much closer. The difference between the two sources is that one of them shrinks when you buy more animals and the other does not.

n_now <- length(conc) * n_per
round(c(animals_now = n_now,
        animals_for_lc10_parity =
          n_now * (log(median(cal[, "wid10"])) / log(sys_all[idx[2]]))^2,
        animals_for_lc1_parity =
          n_now * (log(median(cal[, "wid1"])) / log(sys_all[idx[4]]))^2), 1)
            animals_now animals_for_lc10_parity  animals_for_lc1_parity 
                  175.0                  1872.6                   307.9

Treating the interval width as scaling with the square root of the number of animals, the two sources reach parity at about 1872.6 animals for the LC10 and about 307.9 animals for the LC1. The comparison is crude, since an interval width and a spread between three point estimates are not the same object, but the scale is the point: at the LC1 an experiment less than twice the size of this one is already limited by a choice the analyst makes rather than by the animals. Reporting a confidence interval from one link, at that level, understates the uncertainty by a factor that no amount of extra testing will remove.

Three further limits sit outside what was measured here. The generating curve was itself one of the three links under test, which is generous: a real dose-response surface has no reason to be exactly probit, logit or complementary log-log, so the spread measured above is a lower bound on model-form uncertainty rather than an estimate of it. The animals were treated as independent binomial trials, which they are not when they share a vessel, and vessel-level clustering makes every interval here too narrow. And the concentrations were assumed to be known exactly and the control mortality to be zero. Those three assumptions are the subject of Checking a dose-response analysis, and each of them moves the low-effect end of the curve more than the middle, in the same direction as everything measured here.

Where to go next

The intervals in this post were built with Fieller’s theorem and checked once for coverage, which is enough to establish that the estimates are calibrated but not enough to choose an interval method. Confidence intervals for effective doses takes that apart properly, including the design question this post left alone: where to place the concentrations if the low-effect end is what the report will quote.

Two other directions lead out of the same fit. Everything above assumes the response rises monotonically with concentration, and Hormesis and non-monotonic responses measures what a monotone model does to an effective dose when the data carry a low-dose stimulation instead. And the curves here were fitted as generalised linear models because a link function made that possible; when the shape needs an asymptote or an extra parameter, the fit becomes a nonlinear least squares or likelihood problem of the kind in Nonlinear regression with nls.

References

Ritz C 2010 Environmental Toxicology and Chemistry 29(1):220-229 (10.1002/etc.7)

Ritz C, Baty F, Streibig JC, Gerhard D 2015 PLOS ONE 10(12):e0146021 (10.1371/journal.pone.0146021)

Fieller EC 1954 Journal of the Royal Statistical Society Series B 16(2):175-185 (10.1111/j.2517-6161.1954.tb00159.x)

Finney DJ 1971 Probit Analysis. 3rd edition. Cambridge University Press, ISBN 978-0-521-08041-5

Collett D 2003 Modelling Binary Data. 2nd edition. Chapman and Hall/CRC, ISBN 978-1-58488-324-1

Newsletter

Get new tutorials by email

New R and QGIS tutorials for ecologists, straight to your inbox. No spam; unsubscribe anytime.

By subscribing you agree to receive these emails and confirm your address once. See the privacy policy.