Force of infection from age prevalence

disease ecology

epidemiology

GLM

ecology tutorial

A cloglog binomial GLM with a log age offset returns the force of infection directly. Fitting catalytic models to wildlife serology data in R, from scratch.

Author

Tidy Ecology

Published

2026-08-09

Modified

2026-10-15

Corrected 15 October 2026: the post said that a single cross-section cannot separate disease induced mortality from antibody waning, although its own simulation picks the right mechanism in most surveys; it now says that a cross-section of ordinary size cannot separate them reliably.

Corrected 7 October 2026: the limits section said the catalytic model is more sensitive to age error than a logistic regression because its offset ties the curve to the age scale rather than to the ordering of ages; a logistic regression on age uses the age values as well, so the section now says only that the catalytic model is not exempt and that the size of the effect was not simulated.

Corrected 30 September 2026: the catalytic-model section said the residual deviance shows that a constant rate describes the data; it now says the deviance check finds no evidence against a constant rate, which does not by itself establish one. The band-model section also said that no amount of extra sampling in old age classes would narrow the interval on the oldest band; it now says that prevalence there has not saturated and that sampling more old animals is the direct way to narrow it.

Updated 25 September 2026: a new section, Maternal antibody in the youngest class, measures what the usual repair of dropping the youngest animals leaves behind, and What to report no longer says that the fit tests none of the host assumptions.

A serological survey gives one number per animal: seropositive or not, plus an age. The result is usually summarised as a prevalence, and prevalence is a stock. It says what fraction of the population currently carries the mark of infection, not how fast animals are acquiring it.

The rate is what the ecology usually turns on, and a single cross-section can supply it, because age is time. An animal aged five has been exposed for five years, and the fraction of five year olds carrying antibodies is the probability of having been infected at least once in five years of exposure. That is a survival problem wearing different clothes, and its answer is a hazard: the force of infection.

The model that makes that translation is nearly a century old, and in R it is one call to glm() with a link most people never use and an offset. The purpose of this post is to write that line, show why the logistic regression people reach for instead is a different model, extend it when the rate is not constant, and then describe the one thing a cross-section cannot do reliably.

A cross-section of a wildlife host

Twelve age classes, one serological result per animal, no repeat sampling. The generating process is an irreversible infection acquired at a constant rate, so the truth is known.

library(ggplot2)
library(patchwork)

te_paper  <- "#f5f4ee"
te_ink    <- "#16241d"
te_body   <- "#2c3a31"
te_forest <- "#275139"
te_rust   <- "#b5534e"
te_gold   <- "#c9b458"
te_line   <- "#dad9ca"

theme_datasheet <- function() {
  theme_minimal(base_size = 12) +
    theme(plot.background  = element_rect(fill = te_paper, colour = NA),
          panel.background = element_rect(fill = te_paper, colour = NA),
          panel.grid.major = element_line(colour = te_line, linewidth = 0.3),
          panel.grid.minor = element_blank(),
          text             = element_text(colour = te_body),
          plot.title       = element_text(colour = te_ink, face = "bold"),
          plot.subtitle    = element_text(colour = te_body),
          axis.text        = element_text(colour = te_body))
}

foi_true <- 0.22
n_class  <- 12
sampled  <- c(180, 160, 140, 120, 100, 85, 70, 55, 45, 35, 25, 20)

set.seed(17)
sero <- data.frame(age = seq_len(n_class), n = sampled)
sero$pos  <- rbinom(n_class, sero$n, 1 - exp(-foi_true * sero$age))
sero$neg  <- sero$n - sero$pos
sero$prev <- sero$pos / sero$n

wilson <- function(x, n, z = 1.96) {
  ph <- x / n
  mid <- (ph + z^2 / (2 * n)) / (1 + z^2 / n)
  hw  <- z * sqrt(ph * (1 - ph) / n + z^2 / (4 * n^2)) / (1 + z^2 / n)
  cbind(lo = pmax(0, mid - hw), hi = pmin(1, mid + hw))
}
sero[, c("lo", "hi")] <- wilson(sero$pos, sero$n)
n_total <- sum(sero$n)

1035 animals, 561 of them seropositive, 54.2 per cent overall. That headline figure is the stock. It mixes a lightly exposed cohort of yearlings with animals that have had a decade to meet the pathogen, and on its own it says nothing about transmission.

The test is treated as perfect here. Correcting an observed prevalence for imperfect sensitivity and specificity is a separate step covered elsewhere on the site, and it should be done before anything below.

The catalytic model is one line of glm

If infection arrives as a Poisson process of constant rate lambda and never reverses, then the probability that an animal aged a has been infected at least once is one minus the probability of no events, which is 1 - exp(-lambda * a). Rearranging gives the whole model:

log(-log(1 - p)) = log(lambda) + log(a)

The left side is the complementary log-log of the prevalence. The right side is an intercept plus a term in log(a) whose coefficient is fixed at one. So a binomial GLM with a cloglog link and log(age) as an offset returns log(lambda) as its intercept, with a standard error and everything else a GLM gives.

m_cat <- glm(cbind(pos, neg) ~ 1 + offset(log(age)), data = sero,
             family = binomial(link = "cloglog"))

foi_hat <- unname(exp(coef(m_cat)[1]))
foi_ci  <- unname(exp(confint.default(m_cat)[1, ]))
se_log  <- summary(m_cat)$coefficients[1, 2]

nll_const <- function(log_l) {
  p <- 1 - exp(-exp(log_l) * sero$age)
  -sum(sero$pos * log(p) + sero$neg * log(1 - p))
}
foi_optim <- exp(optimize(nll_const, c(-5, 2))$minimum)
dev_cat   <- deviance(m_cat)
df_cat    <- m_cat$df.residual
ci_pct    <- 95
mean_age_inf <- 1 / foi_hat

The estimate is 0.211 infections per animal per year, with a 95 per cent interval of 0.193 to 0.231 against a true value of 0.22. Writing the likelihood out by hand and minimising it returns 0.211, which is the same number: the offset trick is not an approximation to the catalytic model, it is the catalytic model.

The intercept is a log rate, not a log odds. Its exponent has units of one over time, its reciprocal 4.7 years is the mean age at first infection, and that quantity connects to the reproduction number by an identity the site derives in the context of recurrent epidemics. The residual deviance is 17.7 on 11 degrees of freedom, so the deviance check finds no evidence against a constant rate; as the section on maternal antibody below shows, passing that check does not by itself establish that the rate is constant.

A logistic regression on age is a different model

The reflex is glm(cbind(pos, neg) ~ age, family = binomial). It has one more parameter than the catalytic model, it will produce a rising curve, and it will look fine.

m_logit <- glm(cbind(pos, neg) ~ age, data = sero, family = binomial)

dev_logit  <- deviance(m_logit)
df_logit   <- m_logit$df.residual
aic_gap    <- AIC(m_logit) - AIC(m_cat)
p_at_birth <- unname(plogis(coef(m_logit)[1]))
logit_slope <- unname(coef(m_logit)[2])
slope_ratio <- logit_slope / foi_hat

It fits worse while spending more: a residual deviance of 37.2 on 10 degrees of freedom against 17.7 on 11, and an AIC 21.5 higher. That gap is not a small preference. It is the model saying that the shape imposed by the logistic link is the wrong shape for a process of accumulating exposure.

The clearest symptom is at the left edge. The logistic fit puts prevalence at birth at 20.4 per cent, which for an infection acquired after birth is an impossible claim, and the catalytic model puts it at zero because the algebra forces it to. The slope, 0.366 per year, is a change in log odds and reading it as a rate would overstate the force of infection by a factor of 1.7.

grid_age <- data.frame(age = seq(0, n_class, length.out = 200))
grid_age$catalytic <- 1 - exp(-foi_hat * grid_age$age)
grid_age$logistic  <- predict(m_logit, grid_age, type = "response")
fit_long <- data.frame(
  age = rep(grid_age$age, 2),
  p   = c(grid_age$catalytic, grid_age$logistic),
  fit = rep(c("catalytic (cloglog, log age offset)", "logistic in age"),
            each = nrow(grid_age)))

ggplot(sero, aes(age, prev)) +
  geom_errorbar(aes(ymin = lo, ymax = hi), width = 0.18,
                colour = te_body, linewidth = 0.4) +
  geom_point(size = 2.4, colour = te_ink) +
  geom_line(data = fit_long, aes(age, p, colour = fit), linewidth = 0.9) +
  scale_colour_manual(values = c(te_forest, te_rust), name = NULL) +
  scale_y_continuous(limits = c(0, 1)) +
  labs(x = "age (years)", y = "seroprevalence",
       title = "Two curves through the same points",
       subtitle = "one of them says a fifth of newborns are already positive") +
  theme_datasheet() +
  theme(legend.position = "bottom")

A plot of seroprevalence against age from one to twelve years, with points and vertical error bars. A green catalytic curve rises from the origin and bends towards a plateau near ninety per cent. A red logistic curve starts at about twenty per cent at age zero, crosses the points, and rises more steeply through the middle ages. — Figure 1: Observed seroprevalence by age with Wilson intervals, and two fitted curves.

When the rate is not constant

A constant force of infection is a hypothesis, and it is often wrong: juveniles may be protected by maternal antibody, or exposure may peak when animals disperse. The offset trick does not extend directly, because a piecewise constant rate enters through a cumulative hazard that is a sum rather than a product, so the model is no longer log linear. Writing the likelihood out takes six lines.

band_edge <- c(0, 3, 7, 14)
band_foi  <- c(0.07, 0.30, 0.14)

cum_hazard <- function(a, edge, lam) {
  vapply(a, function(x)
    sum(lam * pmax(pmin(x, edge[-1]) - edge[-length(edge)], 0)), 0)
}

n_class2 <- 13
sampled2 <- c(200, 180, 165, 150, 135, 120, 105, 95, 85, 78, 70, 62, 55)
set.seed(23)
vary <- data.frame(age = seq_len(n_class2), n = sampled2)
vary$pos  <- rbinom(n_class2, vary$n, 1 - exp(-cum_hazard(vary$age, band_edge, band_foi)))
vary$neg  <- vary$n - vary$pos
vary$prev <- vary$pos / vary$n
vary[, c("lo", "hi")] <- wilson(vary$pos, vary$n)

m_flat <- glm(cbind(pos, neg) ~ 1 + offset(log(age)), data = vary,
              family = binomial(link = "cloglog"))
foi_flat <- unname(exp(coef(m_flat)[1]))

nll_bands <- function(par) {
  p <- 1 - exp(-cum_hazard(vary$age, band_edge, exp(par)))
  p <- pmin(pmax(p, 1e-12), 1 - 1e-12)
  -sum(vary$pos * log(p) + vary$neg * log(1 - p))
}
fit_bands <- optim(rep(log(0.2), length(band_foi)), nll_bands,
                   method = "BFGS", hessian = TRUE)
stopifnot(fit_bands$convergence == 0)   # check before inverting the Hessian
foi_bands <- exp(fit_bands$par)
se_bands  <- sqrt(diag(solve(fit_bands$hessian)))
band_lo   <- exp(fit_bands$par - 1.96 * se_bands)
band_hi   <- exp(fit_bands$par + 1.96 * se_bands)

deviance_of <- function(p, obs) {
  p <- pmin(pmax(p, 1e-12), 1 - 1e-12)
  2 * sum(ifelse(obs$pos > 0, obs$pos * log(obs$prev / p), 0) +
          ifelse(obs$neg > 0, obs$neg * log((1 - obs$prev) / (1 - p)), 0))
}
dev_bands <- deviance_of(1 - exp(-cum_hazard(vary$age, band_edge, foi_bands)), vary)
df_bands  <- n_class2 - length(band_foi)
lrt_stat  <- deviance(m_flat) - dev_bands
lrt_df    <- length(band_foi) - 1
lrt_p     <- pchisq(lrt_stat, lrt_df, lower.tail = FALSE)
band_width <- (band_hi - band_lo) / foi_bands
dev_flat   <- deviance(m_flat)
df_flat    <- m_flat$df.residual

Forcing a constant rate on these data returns 0.156 with a residual deviance of 94.7 on 12 degrees of freedom, which is a decisive rejection. The three band model recovers 0.070, 0.323 and 0.128 against generating values of 0.07, 0.30 and 0.14, with a deviance of 16.5 on 10 degrees of freedom. The likelihood ratio against the constant model is 78.2 on 2 degrees of freedom.

grid2 <- data.frame(age = seq(0, n_class2, length.out = 300))
grid2$band <- 1 - exp(-cum_hazard(grid2$age, band_edge, foi_bands))
grid2$flat <- 1 - exp(-foi_flat * grid2$age)
fit2 <- data.frame(age = rep(grid2$age, 2), p = c(grid2$band, grid2$flat),
                   fit = rep(c("three bands", "constant rate"), each = nrow(grid2)))

p_curve <- ggplot(vary, aes(age, prev)) +
  geom_errorbar(aes(ymin = lo, ymax = hi), width = 0.2,
                colour = te_body, linewidth = 0.4) +
  geom_point(size = 2.2, colour = te_ink) +
  geom_line(data = fit2, aes(age, p, colour = fit), linewidth = 0.9) +
  scale_colour_manual(values = c(te_rust, te_forest), name = NULL) +
  scale_y_continuous(limits = c(0, 1)) +
  labs(x = "age (years)", y = "seroprevalence", title = "Prevalence") +
  theme_datasheet() + theme(legend.position = "bottom")

band_tab <- data.frame(
  band = factor(sprintf("%g to %g", band_edge[-length(band_edge)], band_edge[-1]),
                levels = sprintf("%g to %g", band_edge[-length(band_edge)], band_edge[-1])),
  est = foi_bands, lo = band_lo, hi = band_hi, truth = band_foi)

p_band <- ggplot(band_tab, aes(band, est)) +
  geom_errorbar(aes(ymin = lo, ymax = hi), width = 0.14,
                colour = te_forest, linewidth = 0.6) +
  geom_point(size = 3, colour = te_forest) +
  geom_point(aes(y = truth), size = 2.6, shape = 4, colour = te_rust, stroke = 1.1) +
  labs(x = "age band (years)", y = "force of infection per year",
       title = "Rate by band", subtitle = "red crosses: generating values") +
  theme_datasheet()

p_curve + p_band + plot_annotation(theme = theme_datasheet())

Two panels. The left panel shows seroprevalence points with error bars, a green stepped curve from the three band model tracking them closely, and a red constant rate curve that sits far above the youngest ages and below the middle ones. The right panel shows the three fitted band rates as points with vertical intervals and red crosses at the generating values, the middle band highest and the oldest band carrying by far the tallest interval. — Figure 2: An age varying force of infection, and what each band is worth.

The right panel carries the warning. The interval on the oldest band spans 1.4 times its own estimate, against 0.5 for the youngest, and part of the reason is visible in the left panel: by age seven about three quarters of the animals are already seropositive, so comparatively few in the oldest band remain to seroconvert, and the design samples the fewest animals at those ages. A cross-section carries information about the rate only where seronegative animals remain. That makes old animals expensive to learn from, not useless: prevalence here has not saturated, and per animal the oldest classes carry the most information about the oldest band, so sampling more of them is the direct way to narrow that interval. Only as prevalence approaches one does an extra old animal add almost nothing.

There is a second constraint worth naming. In this model prevalence can only rise with age, because the cumulative hazard can only accumulate. A survey that shows prevalence falling in old animals cannot be explained by any non negative force of infection, and the explanation has to come from somewhere else.

Maternal antibody in the youngest class

A young animal can test positive because its mother was infected. Peel and colleagues (2018) fitted bat serology with a model that keeps maternally derived and infection acquired antibody apart, and a simpler joint model of the same kind is one of the three fits compared here. The host is the one from the first section, with a constant rate of 0.22, plus a class of 345 juveniles aged half a year. A newborn carries maternal antibody with the probability that its mother is seropositive, taken as the prevalence at age four, and that antibody stops being detectable with a half-life h. Prevalence at age a is then

p(a) = 1 - exp(-lambda * a) * (1 - pi_m * 2^(-a / h))

The three fits are the catalytic GLM on every class, the same GLM after the usual repair of dropping animals under a year, and a likelihood with the maternal term that estimates lambda, pi_m and h together (from three starting points, keeping the best); the same GLM without animals under two is reported alongside as a stricter version of the repair. The bias of the catalytic fits is arithmetic: fitting them to the expected counts, with no sampling noise at all, gives the value each one converges on as the survey grows. The simulation is needed only for interval coverage and for how often the deviance check notices.

mat_age <- c(0.5, seq_len(n_class))
mat_n   <- c(345, sampled)
mat_pim <- 1 - exp(-foi_true * 4)
mat_h   <- c(0.5, 1)

p_maternal <- function(a, lam, pim, h) 1 - exp(-lam * a) * (1 - pim * 2^(-a / h))

pseudo_true <- function(keep, q) {
  nll <- function(log_l) {
    p <- 1 - exp(-exp(log_l) * mat_age[keep])
    -sum(mat_n[keep] * (q[keep] * log(p) + (1 - q[keep]) * log(1 - p)))
  }
  exp(optimize(nll, c(-6, 2))$minimum) / foi_true
}
mat_ratio <- sapply(mat_h, function(h) {
  q <- p_maternal(mat_age, foi_true, mat_pim, h)
  c(naive = pseudo_true(mat_age > 0, q), drop1 = pseudo_true(mat_age >= 1, q),
    drop2 = pseudo_true(mat_age >= 2, q))
})

fit_maternal <- function(pos) {
  obs <- data.frame(age = mat_age, pos = pos, neg = mat_n - pos)
  wald <- function(m) {
    s <- summary(m)$coefficients[1, 1:2]
    c(exp(s[1] + c(0, -1.96, 1.96) * s[2]),
      pchisq(deviance(m), m$df.residual, lower.tail = FALSE))
  }
  m_naive <- glm(cbind(pos, neg) ~ 1 + offset(log(age)), data = obs,
                 family = binomial(link = "cloglog"))
  m_drop  <- glm(cbind(pos, neg) ~ 1 + offset(log(age)), data = obs[obs$age >= 1, ],
                 family = binomial(link = "cloglog"))
  m_drop2 <- glm(cbind(pos, neg) ~ 1 + offset(log(age)), data = obs[obs$age >= 2, ],
                 family = binomial(link = "cloglog"))
  nll <- function(par) {
    p <- p_maternal(mat_age, exp(par[1]), plogis(par[2]), exp(par[3]))
    p <- pmin(pmax(p, 1e-12), 1 - 1e-12)
    -sum(pos * log(p) + (mat_n - pos) * log(1 - p))
  }
  starts <- list(c(log(0.2), 0, log(0.5)), c(log(0.2), 1, 0), c(log(0.3), -1, log(0.25)))
  fits <- lapply(starts, function(s) optim(s, nll, method = "BFGS", hessian = TRUE))
  best <- fits[[which.min(vapply(fits, function(f) f$value, 0))]]
  se_ml <- sqrt(solve(best$hessian)[1, 1])
  c(wald(m_naive), wald(m_drop), exp(best$par[1] + c(0, -1.96, 1.96) * se_ml),
    wald(m_drop2))
}

set.seed(2018)
mat_demo <- fit_maternal(rbinom(length(mat_age), mat_n, p_maternal(mat_age, foi_true, mat_pim, 1)))

n_mat_rep <- 200
mat_sim <- lapply(mat_h, function(h)
  replicate(n_mat_rep, fit_maternal(rbinom(length(mat_age), mat_n,
                                           p_maternal(mat_age, foi_true, mat_pim, h)))))
covers <- function(res, i) mean(res[i + 1, ] <= foi_true & res[i + 2, ] >= foi_true)
mat_tab <- sapply(mat_sim, function(res) c(
  cover_naive = covers(res, 1), cover_drop = covers(res, 5), cover_ml = covers(res, 9),
  cover_drop2 = covers(res, 12),
  flag_naive  = mean(res[4, ] < 0.05), flag_drop = mean(res[8, ] < 0.05),
  pass_miss   = mean(res[8, ] >= 0.05 & (res[6, ] > foi_true | res[7, ] < foi_true)),
  width_drop  = median(res[7, ] - res[6, ]), width_ml = median(res[11, ] - res[10, ])))
mat_p1       <- 1 - exp(-foi_true)
mat_excess1  <- sapply(mat_h, function(h) p_maternal(1, foi_true, mat_pim, h)) - mat_p1
mat_widen    <- mat_tab["width_ml", ] / mat_tab["width_drop", ] - 1
mat_share    <- mat_n[1] / sum(mat_n)
mat_duration <- mat_h / log(2)
mat_cover_se <- sqrt(0.25 / n_mat_rep)

The juveniles make up 25 per cent of the 1380 animals, and the prevalence in mothers is 0.585. With a half-life of half a year, the catalytic fit on every class converges on 1.248 times the true rate and the fit without juveniles on 1.061 times; with a half-life of one year the two limits are 1.442 and 1.165. Dropping the juveniles does not remove the problem, because the next class is aged exactly one, when a quarter of the maternal antibody is still detectable at h = 0.5 and half of it at h = 1. That adds 11.7 and 23.5 percentage points to a catalytic prevalence of 19.7 per cent at age one, and a single rate has to absorb the excess. Dropping everything under two years brings the limits to 1.013 and 1.068.

Over 200 simulated surveys at each half-life (Monte Carlo standard error at most 3.5 points), the interval from the fit on every class covers the true rate in 0.0 and 0.0 per cent of surveys, and the deviance check of the catalytic fit above flags 100.0 and 100.0 per cent of those fits: ignoring maternal antibody altogether is caught. The fit without juveniles covers in 77.5 per cent at h = 0.5 and 7.5 per cent at h = 1. At h = 1 its own deviance check flags 99.5 per cent of surveys, so in a survey of this size the problem announces itself. At h = 0.5 the check flags 65.5 per cent, and in 8.5 per cent of all surveys the reduced fit passes its own check while its interval misses the true rate. Dropping everything under two years covers in 95.0 and 65.0 per cent. The joint model covers in 96.5 and 95.0 per cent, and its median interval is 7 and 14 per cent wider than the interval of the fit without juveniles. In a single survey drawn at h = 1 the three estimates are 0.316 (0.292 to 0.343), 0.255 (0.233 to 0.278) and 0.213 (0.189 to 0.240) against 0.22, and both catalytic fits fail the deviance check with p below 0.001.

How long maternal antibody stays detectable decides which of these cases applies. Epstein and colleagues (2013) followed pteropid bat pups and found antibody against canine distemper virus and Hendra virus waning over a mean of about 229 and 255 days. Under the exponential loss assumed here, a half-life of half a year means a mean detectable duration of 263 days, so h = 0.5 is the realistic case for such hosts and h = 1 (1.44 years) is the long tail. Even the realistic case leaves the fit without juveniles biased in a survey of this size, because the yearlings still carry maternal antibody. Either drop the yearlings as well, which is enough here at h = 0.5 and not at h = 1, or fit the maternal term, which costs some interval width and needs no guess at the half-life.

Two mechanisms, one curve

What the catalytic model leaves out is the honest limit of the whole exercise. Two processes it leaves out flatten an age prevalence curve below one, and a single cross-section of ordinary size cannot separate them reliably. Infection may kill: seropositive animals are removed faster, so the older age classes are enriched with survivors who were never infected. Or antibodies may wane: animals revert to seronegative and are counted as never infected.

p_mortality <- function(a, lam, alpha) {
  inf <- lam * (exp(-alpha * a) - exp(-lam * a)) / (lam - alpha)
  sus <- exp(-lam * a)
  inf / (sus + inf)
}
p_reversion <- function(a, lam, rho) lam / (lam + rho) * (1 - exp(-(lam + rho) * a))

age_grid  <- seq_len(15)
foi_mort  <- 0.30
alpha_mort <- 0.45
curve_mort <- p_mortality(age_grid, foi_mort, alpha_mort)

match_rev <- optim(log(c(0.3, 0.2)), function(par)
  sum((p_reversion(age_grid, exp(par[1]), exp(par[2])) - curve_mort)^2),
  method = "BFGS")
rev_par    <- exp(match_rev$par)
curve_rev  <- p_reversion(age_grid, rev_par[1], rev_par[2])
max_gap    <- max(abs(curve_rev - curve_mort))
half_life  <- log(2) / rev_par[2]
years_left <- 1 / alpha_mort

The best matching seroreversion model has a force of infection of 0.228 and a reversion rate of 0.135, and its curve differs from the mortality curve by at most 2.7 percentage points anywhere between ages one and 15. The two describe entirely different animals. One says infection removes a host in 2.2 years on average; the other says infection is harmless and antibodies have a half life of 5.1 years.

n_rep    <- 200
sampled3 <- c(200, 180, 165, 150, 135, 120, 105, 95, 85, 78, 70, 62, 55, 48, 42)

fit_two <- function(pos) {
  nll <- function(par, fun) {
    p <- fun(age_grid, exp(par[1]), exp(par[2]))
    p <- pmin(pmax(p, 1e-9), 1 - 1e-9)
    -sum(pos * log(p) + (sampled3 - pos) * log(1 - p))
  }
  c(2 * optim(log(c(0.3, 0.4)), nll, fun = p_mortality,  method = "BFGS")$value,
    2 * optim(log(c(0.3, 0.2)), nll, fun = p_reversion, method = "BFGS")$value)
}

set.seed(5)
two_fit <- replicate(n_rep, fit_two(rbinom(length(age_grid), sampled3, curve_mort)))
aic_diff <- two_fit[2, ] - two_fit[1, ]
wrong_wins <- mean(aic_diff < 0)
n_survey   <- sum(sampled3)
aic_mean   <- mean(aic_diff)
aic_sd     <- sd(aic_diff)
bern_var   <- 0.25
share_se   <- 100 * sqrt(bern_var / n_rep)

With 1590 animals in the survey, generated from the mortality model, the mean AIC advantage of the correct model is 1.4 points with a standard deviation of 2.4, and the wrong model has the lower AIC in 30 per cent of the surveys. This is not a power problem that a bigger sample fixes cheaply; it is two mechanisms writing nearly the same curve.

conf_long <- data.frame(
  age = rep(age_grid, 2), p = c(curve_mort, curve_rev),
  mechanism = rep(c("infection kills", "antibodies wane"), each = length(age_grid)))

p_conf <- ggplot(conf_long, aes(age, p, colour = mechanism)) +
  geom_line(linewidth = 1) + geom_point(size = 2) +
  scale_colour_manual(values = c(te_gold, te_forest), name = NULL) +
  scale_y_continuous(limits = c(0, 1)) +
  labs(x = "age (years)", y = "seroprevalence",
       title = "Same curve, opposite biology") +
  theme_datasheet() + theme(legend.position = "bottom")

p_aic <- ggplot(data.frame(d = aic_diff), aes(d)) +
  geom_histogram(bins = 30, fill = te_forest, colour = te_paper, linewidth = 0.2) +
  geom_vline(xintercept = 0, colour = te_rust, linetype = "dashed", linewidth = 0.8) +
  labs(x = "AIC of waning minus AIC of mortality", y = "surveys",
       title = "And the data choose poorly",
       subtitle = "dashed red: no preference") +
  theme_datasheet()

p_conf + p_aic + plot_annotation(theme = theme_datasheet())

Two panels. The left panel plots two nearly overlapping curves of prevalence against age, one for extra mortality of the infected and one for antibody waning, rising to a plateau just above sixty per cent. The right panel is a histogram of the AIC difference between the two models across two hundred simulated surveys, centred just above zero with a substantial part of its mass below zero. — Figure 3: Disease induced mortality and antibody waning, fitted to produce the same age prevalence curve.

What to report

Give the age distribution of the sample, not only the prevalence. The number of seronegative animals per age class is what carries the information, and a survey that is mostly old animals can have a large sample size and almost no ability to estimate a rate.

Fit the catalytic model rather than a logistic regression on age. It is family = binomial(link = "cloglog") with offset(log(age)), its intercept exponentiates to a rate with units, and it cannot produce a positive prevalence at birth.

Report the force of infection with an interval, and report the age bands if the rate was allowed to vary. Band edges are a modelling choice, and the estimate in a band where prevalence has saturated should be quoted with its width, which will be embarrassing and should be.

Say which biological assumptions the estimate rests on: no seroreversion, no differential mortality, no maternal antibody, and a population at steady state. Every one of them is a claim about the host. The deviance check can show that the curve is not catalytic, as it did for the fit that ignored maternal antibody in 100.0 per cent of the simulated surveys above, but it cannot say which assumption failed, and a fit that passes it is not cleared: at h = 0.5 the fit without juveniles passed in 34.5 per cent of surveys, although it converges on 1.061 times the true rate. If the sample includes animals young enough to carry maternal antibody, say how they were handled.

If mean age at infection or a reproduction number is derived from the rate, give the demographic quantities that went into the derivation as well. Those identities are covered separately on the site and they are only as good as the mortality schedule they assume.

Honest limits

The largest limitation is the one measured above: a cross-section of ordinary size cannot reliably separate disease induced mortality from antibody waning, and both are common in wildlife. Longitudinal recaptures, or a test that distinguishes recent from historical infection, or independent mortality data, are the ways out, and all of them are more expensive than a survey.

Every model here assumes the population is at steady state and that the force of infection has been what it is now for the whole lifespan of the oldest animal sampled. A cross-section reads age as calendar time backwards, so a rate that changed five years ago appears as an age effect in exactly the same place as a genuine age effect. The two are not separable from one survey either, and the standard fix, fitting age and cohort together, needs surveys from more than one year.

Maternal antibody is measured above for one constant rate, one sample design and an exponential loss of detectability. A seropositive call is a threshold on a decaying titre, so for one animal the loss is a step when its titre crosses the threshold, and across animals it follows the spread of those crossing times, which need not be exponential. The yearlings here are sampled at exactly one year, the youngest age the usual repair keeps and so the age at which the kept animals carry the most maternal antibody; a birth pulse sampled once a year puts the classes at 0.5, 1.5, 2.5 years, where the fit without juveniles is less biased, and neither that design nor the joint model’s behaviour in it was simulated. In the band model the maternal positives would land in the first band, and dropping the youngest class there removes data from the band that is otherwise best determined; that version was not simulated.

Ages are treated as exact. Field ages usually come from tooth wear or horn rings and carry error that a regression on age propagates in the familiar way, which the site covers under measurement error. The catalytic model is not exempt: its offset takes the recorded ages as they stand, and how far age error moves the estimated rate was not simulated here.

The band model uses fixed edges chosen with knowledge of the truth. Choosing edges from the same data that fits the rates makes the intervals optimistic, and the honest alternatives are a smooth on the cloglog scale or edges chosen from the host’s biology before the data are seen.

The comparison of the two mechanisms uses 200 simulated surveys, so the reported share of surveys favouring the wrong model has a standard error of about 3.5 percentage points. The conclusion that the two models are nearly indistinguishable does not rest on that share; it rests on the 2.7 percentage point maximum gap between the curves, which involves no sampling at all.

References

Heisey DM, Joly DO, Messier F 2006 Ecology 87(9):2356-2365 (10.1890/0012-9658(2006)87[2356:TFOGFM]2.0.CO;2)

Hens N, Aerts M, Faes C, Shkedy Z, Lejeune O, Van Damme P, Beutels P 2009 Epidemiology and Infection 138(6):802-812 (10.1017/S0950268809990781)

Grenfell BT, Anderson RM 1985 Journal of Hygiene 95(2):419-436 (10.1017/S0022172400062859)

Peel AJ, Baker KS, Hayman DTS, Broder CC, Cunningham AA, et al. 2018 Scientific Reports 8:3859 (10.1038/s41598-018-22236-6)

Epstein JH, Baker ML, Zambrana-Torrelio C, et al. 2013 PLoS ONE 8(6):e67584 (10.1371/journal.pone.0067584)

--- title: "Force of infection from age prevalence" description: "A cloglog binomial GLM with a log age offset returns the force of infection directly. Fitting catalytic models to wildlife serology data in R, from scratch." date: "2026-08-09 10:00" date-modified: "2026-10-15" categories: [R, disease ecology, epidemiology, GLM, ecology tutorial] image: thumbnail.png image-alt: "A curve rising from the origin towards a plateau on a warm off-white background, with observed seroprevalence points and error bars by age class, and a straighter logistic curve that starts well above zero at age zero." --- *Corrected 15 October 2026: the post said that a single cross-section cannot separate disease induced mortality from antibody waning, although its own simulation picks the right mechanism in most surveys; it now says that a cross-section of ordinary size cannot separate them reliably.* *Corrected 7 October 2026: the limits section said the catalytic model is more sensitive to age error than a logistic regression because its offset ties the curve to the age scale rather than to the ordering of ages; a logistic regression on age uses the age values as well, so the section now says only that the catalytic model is not exempt and that the size of the effect was not simulated.* *Corrected 30 September 2026: the catalytic-model section said the residual deviance shows that a constant rate describes the data; it now says the deviance check finds no evidence against a constant rate, which does not by itself establish one. The band-model section also said that no amount of extra sampling in old age classes would narrow the interval on the oldest band; it now says that prevalence there has not saturated and that sampling more old animals is the direct way to narrow it.* *Updated 25 September 2026: a new section, Maternal antibody in the youngest class, measures what the usual repair of dropping the youngest animals leaves behind, and What to report no longer says that the fit tests none of the host assumptions.* A serological survey gives one number per animal: seropositive or not, plus an age. The result is usually summarised as a prevalence, and prevalence is a stock. It says what fraction of the population currently carries the mark of infection, not how fast animals are acquiring it. The rate is what the ecology usually turns on, and a single cross-section can supply it, because age is time. An animal aged five has been exposed for five years, and the fraction of five year olds carrying antibodies is the probability of having been infected at least once in five years of exposure. That is a survival problem wearing different clothes, and its answer is a hazard: the force of infection. The model that makes that translation is nearly a century old, and in R it is one call to `glm()` with a link most people never use and an offset. The purpose of this post is to write that line, show why the logistic regression people reach for instead is a different model, extend it when the rate is not constant, and then describe the one thing a cross-section cannot do reliably. ## A cross-section of a wildlife host Twelve age classes, one serological result per animal, no repeat sampling. The generating process is an irreversible infection acquired at a constant rate, so the truth is known. ```{r setup} #| message: false #| warning: false library(ggplot2) library(patchwork) te_paper <- "#f5f4ee" te_ink <- "#16241d" te_body <- "#2c3a31" te_forest <- "#275139" te_rust <- "#b5534e" te_gold <- "#c9b458" te_line <- "#dad9ca" theme_datasheet <- function() { theme_minimal(base_size = 12) + theme(plot.background = element_rect(fill = te_paper, colour = NA), panel.background = element_rect(fill = te_paper, colour = NA), panel.grid.major = element_line(colour = te_line, linewidth = 0.3), panel.grid.minor = element_blank(), text = element_text(colour = te_body), plot.title = element_text(colour = te_ink, face = "bold"), plot.subtitle = element_text(colour = te_body), axis.text = element_text(colour = te_body)) } foi_true <- 0.22 n_class <- 12 sampled <- c(180, 160, 140, 120, 100, 85, 70, 55, 45, 35, 25, 20) set.seed(17) sero <- data.frame(age = seq_len(n_class), n = sampled) sero$pos <- rbinom(n_class, sero$n, 1 - exp(-foi_true * sero$age)) sero$neg <- sero$n - sero$pos sero$prev <- sero$pos / sero$n wilson <- function(x, n, z = 1.96) { ph <- x / n mid <- (ph + z^2 / (2 * n)) / (1 + z^2 / n) hw <- z * sqrt(ph * (1 - ph) / n + z^2 / (4 * n^2)) / (1 + z^2 / n) cbind(lo = pmax(0, mid - hw), hi = pmin(1, mid + hw)) } sero[, c("lo", "hi")] <- wilson(sero$pos, sero$n) n_total <- sum(sero$n) ``` `r sprintf("%d", n_total)` animals, `r sprintf("%d", sum(sero$pos))` of them seropositive, `r sprintf("%.1f", 100 * sum(sero$pos) / n_total)` per cent overall. That headline figure is the stock. It mixes a lightly exposed cohort of yearlings with animals that have had a decade to meet the pathogen, and on its own it says nothing about transmission. The test is treated as perfect here. Correcting an observed prevalence for imperfect sensitivity and specificity is a separate step covered elsewhere on the site, and it should be done before anything below. ## The catalytic model is one line of glm If infection arrives as a Poisson process of constant rate `lambda` and never reverses, then the probability that an animal aged `a` has been infected at least once is one minus the probability of no events, which is `1 - exp(-lambda * a)`. Rearranging gives the whole model: `log(-log(1 - p)) = log(lambda) + log(a)` The left side is the complementary log-log of the prevalence. The right side is an intercept plus a term in `log(a)` whose coefficient is fixed at one. So a binomial GLM with a cloglog link and `log(age)` as an offset returns `log(lambda)` as its intercept, with a standard error and everything else a GLM gives. ```{r catalytic} m_cat <- glm(cbind(pos, neg) ~ 1 + offset(log(age)), data = sero, family = binomial(link = "cloglog")) foi_hat <- unname(exp(coef(m_cat)[1])) foi_ci <- unname(exp(confint.default(m_cat)[1, ])) se_log <- summary(m_cat)$coefficients[1, 2] nll_const <- function(log_l) { p <- 1 - exp(-exp(log_l) * sero$age) -sum(sero$pos * log(p) + sero$neg * log(1 - p)) } foi_optim <- exp(optimize(nll_const, c(-5, 2))$minimum) dev_cat <- deviance(m_cat) df_cat <- m_cat$df.residual ci_pct <- 95 mean_age_inf <- 1 / foi_hat ``` The estimate is `r sprintf("%.3f", foi_hat)` infections per animal per year, with a `r sprintf("%d", ci_pct)` per cent interval of `r sprintf("%.3f", foi_ci[1])` to `r sprintf("%.3f", foi_ci[2])` against a true value of `r sprintf("%.2f", foi_true)`. Writing the likelihood out by hand and minimising it returns `r sprintf("%.3f", foi_optim)`, which is the same number: the offset trick is not an approximation to the catalytic model, it is the catalytic model. The intercept is a log rate, not a log odds. Its exponent has units of one over time, its reciprocal `r sprintf("%.1f", mean_age_inf)` years is the mean age at first infection, and that quantity connects to the reproduction number by an identity the site derives in the context of recurrent epidemics. The residual deviance is `r sprintf("%.1f", dev_cat)` on `r sprintf("%d", df_cat)` degrees of freedom, so the deviance check finds no evidence against a constant rate; as the section on maternal antibody below shows, passing that check does not by itself establish that the rate is constant. ## A logistic regression on age is a different model The reflex is `glm(cbind(pos, neg) ~ age, family = binomial)`. It has one more parameter than the catalytic model, it will produce a rising curve, and it will look fine. ```{r logistic} m_logit <- glm(cbind(pos, neg) ~ age, data = sero, family = binomial) dev_logit <- deviance(m_logit) df_logit <- m_logit$df.residual aic_gap <- AIC(m_logit) - AIC(m_cat) p_at_birth <- unname(plogis(coef(m_logit)[1])) logit_slope <- unname(coef(m_logit)[2]) slope_ratio <- logit_slope / foi_hat ``` It fits worse while spending more: a residual deviance of `r sprintf("%.1f", dev_logit)` on `r sprintf("%d", df_logit)` degrees of freedom against `r sprintf("%.1f", dev_cat)` on `r sprintf("%d", df_cat)`, and an AIC `r sprintf("%.1f", aic_gap)` higher. That gap is not a small preference. It is the model saying that the shape imposed by the logistic link is the wrong shape for a process of accumulating exposure. The clearest symptom is at the left edge. The logistic fit puts prevalence at birth at `r sprintf("%.1f", 100 * p_at_birth)` per cent, which for an infection acquired after birth is an impossible claim, and the catalytic model puts it at zero because the algebra forces it to. The slope, `r sprintf("%.3f", logit_slope)` per year, is a change in log odds and reading it as a rate would overstate the force of infection by a factor of `r sprintf("%.1f", slope_ratio)`. ```{r fig-catalytic} #| fig-cap: "Observed seroprevalence by age with Wilson intervals, and two fitted curves." #| fig-alt: "A plot of seroprevalence against age from one to twelve years, with points and vertical error bars. A green catalytic curve rises from the origin and bends towards a plateau near ninety per cent. A red logistic curve starts at about twenty per cent at age zero, crosses the points, and rises more steeply through the middle ages." grid_age <- data.frame(age = seq(0, n_class, length.out = 200)) grid_age$catalytic <- 1 - exp(-foi_hat * grid_age$age) grid_age$logistic <- predict(m_logit, grid_age, type = "response") fit_long <- data.frame( age = rep(grid_age$age, 2), p = c(grid_age$catalytic, grid_age$logistic), fit = rep(c("catalytic (cloglog, log age offset)", "logistic in age"), each = nrow(grid_age))) ggplot(sero, aes(age, prev)) + geom_errorbar(aes(ymin = lo, ymax = hi), width = 0.18, colour = te_body, linewidth = 0.4) + geom_point(size = 2.4, colour = te_ink) + geom_line(data = fit_long, aes(age, p, colour = fit), linewidth = 0.9) + scale_colour_manual(values = c(te_forest, te_rust), name = NULL) + scale_y_continuous(limits = c(0, 1)) + labs(x = "age (years)", y = "seroprevalence", title = "Two curves through the same points", subtitle = "one of them says a fifth of newborns are already positive") + theme_datasheet() + theme(legend.position = "bottom") ``` ## When the rate is not constant A constant force of infection is a hypothesis, and it is often wrong: juveniles may be protected by maternal antibody, or exposure may peak when animals disperse. The offset trick does not extend directly, because a piecewise constant rate enters through a cumulative hazard that is a sum rather than a product, so the model is no longer log linear. Writing the likelihood out takes six lines. ```{r bands} band_edge <- c(0, 3, 7, 14) band_foi <- c(0.07, 0.30, 0.14) cum_hazard <- function(a, edge, lam) { vapply(a, function(x) sum(lam * pmax(pmin(x, edge[-1]) - edge[-length(edge)], 0)), 0) } n_class2 <- 13 sampled2 <- c(200, 180, 165, 150, 135, 120, 105, 95, 85, 78, 70, 62, 55) set.seed(23) vary <- data.frame(age = seq_len(n_class2), n = sampled2) vary$pos <- rbinom(n_class2, vary$n, 1 - exp(-cum_hazard(vary$age, band_edge, band_foi))) vary$neg <- vary$n - vary$pos vary$prev <- vary$pos / vary$n vary[, c("lo", "hi")] <- wilson(vary$pos, vary$n) m_flat <- glm(cbind(pos, neg) ~ 1 + offset(log(age)), data = vary, family = binomial(link = "cloglog")) foi_flat <- unname(exp(coef(m_flat)[1])) nll_bands <- function(par) { p <- 1 - exp(-cum_hazard(vary$age, band_edge, exp(par))) p <- pmin(pmax(p, 1e-12), 1 - 1e-12) -sum(vary$pos * log(p) + vary$neg * log(1 - p)) } fit_bands <- optim(rep(log(0.2), length(band_foi)), nll_bands, method = "BFGS", hessian = TRUE) stopifnot(fit_bands$convergence == 0) # check before inverting the Hessian foi_bands <- exp(fit_bands$par) se_bands <- sqrt(diag(solve(fit_bands$hessian))) band_lo <- exp(fit_bands$par - 1.96 * se_bands) band_hi <- exp(fit_bands$par + 1.96 * se_bands) deviance_of <- function(p, obs) { p <- pmin(pmax(p, 1e-12), 1 - 1e-12) 2 * sum(ifelse(obs$pos > 0, obs$pos * log(obs$prev / p), 0) + ifelse(obs$neg > 0, obs$neg * log((1 - obs$prev) / (1 - p)), 0)) } dev_bands <- deviance_of(1 - exp(-cum_hazard(vary$age, band_edge, foi_bands)), vary) df_bands <- n_class2 - length(band_foi) lrt_stat <- deviance(m_flat) - dev_bands lrt_df <- length(band_foi) - 1 lrt_p <- pchisq(lrt_stat, lrt_df, lower.tail = FALSE) band_width <- (band_hi - band_lo) / foi_bands dev_flat <- deviance(m_flat) df_flat <- m_flat$df.residual ``` Forcing a constant rate on these data returns `r sprintf("%.3f", foi_flat)` with a residual deviance of `r sprintf("%.1f", dev_flat)` on `r sprintf("%d", df_flat)` degrees of freedom, which is a decisive rejection. The three band model recovers `r sprintf("%.3f", foi_bands[1])`, `r sprintf("%.3f", foi_bands[2])` and `r sprintf("%.3f", foi_bands[3])` against generating values of `r sprintf("%.2f", band_foi[1])`, `r sprintf("%.2f", band_foi[2])` and `r sprintf("%.2f", band_foi[3])`, with a deviance of `r sprintf("%.1f", dev_bands)` on `r sprintf("%d", df_bands)` degrees of freedom. The likelihood ratio against the constant model is `r sprintf("%.1f", lrt_stat)` on `r sprintf("%d", lrt_df)` degrees of freedom. ```{r fig-bands} #| fig-cap: "An age varying force of infection, and what each band is worth." #| fig-alt: "Two panels. The left panel shows seroprevalence points with error bars, a green stepped curve from the three band model tracking them closely, and a red constant rate curve that sits far above the youngest ages and below the middle ones. The right panel shows the three fitted band rates as points with vertical intervals and red crosses at the generating values, the middle band highest and the oldest band carrying by far the tallest interval." grid2 <- data.frame(age = seq(0, n_class2, length.out = 300)) grid2$band <- 1 - exp(-cum_hazard(grid2$age, band_edge, foi_bands)) grid2$flat <- 1 - exp(-foi_flat * grid2$age) fit2 <- data.frame(age = rep(grid2$age, 2), p = c(grid2$band, grid2$flat), fit = rep(c("three bands", "constant rate"), each = nrow(grid2))) p_curve <- ggplot(vary, aes(age, prev)) + geom_errorbar(aes(ymin = lo, ymax = hi), width = 0.2, colour = te_body, linewidth = 0.4) + geom_point(size = 2.2, colour = te_ink) + geom_line(data = fit2, aes(age, p, colour = fit), linewidth = 0.9) + scale_colour_manual(values = c(te_rust, te_forest), name = NULL) + scale_y_continuous(limits = c(0, 1)) + labs(x = "age (years)", y = "seroprevalence", title = "Prevalence") + theme_datasheet() + theme(legend.position = "bottom") band_tab <- data.frame( band = factor(sprintf("%g to %g", band_edge[-length(band_edge)], band_edge[-1]), levels = sprintf("%g to %g", band_edge[-length(band_edge)], band_edge[-1])), est = foi_bands, lo = band_lo, hi = band_hi, truth = band_foi) p_band <- ggplot(band_tab, aes(band, est)) + geom_errorbar(aes(ymin = lo, ymax = hi), width = 0.14, colour = te_forest, linewidth = 0.6) + geom_point(size = 3, colour = te_forest) + geom_point(aes(y = truth), size = 2.6, shape = 4, colour = te_rust, stroke = 1.1) + labs(x = "age band (years)", y = "force of infection per year", title = "Rate by band", subtitle = "red crosses: generating values") + theme_datasheet() p_curve + p_band + plot_annotation(theme = theme_datasheet()) ``` The right panel carries the warning. The interval on the oldest band spans `r sprintf("%.1f", band_width[3])` times its own estimate, against `r sprintf("%.1f", band_width[1])` for the youngest, and part of the reason is visible in the left panel: by age seven about three quarters of the animals are already seropositive, so comparatively few in the oldest band remain to seroconvert, and the design samples the fewest animals at those ages. A cross-section carries information about the rate only where seronegative animals remain. That makes old animals expensive to learn from, not useless: prevalence here has not saturated, and per animal the oldest classes carry the most information about the oldest band, so sampling more of them is the direct way to narrow that interval. Only as prevalence approaches one does an extra old animal add almost nothing. There is a second constraint worth naming. In this model prevalence can only rise with age, because the cumulative hazard can only accumulate. A survey that shows prevalence falling in old animals cannot be explained by any non negative force of infection, and the explanation has to come from somewhere else. ## Maternal antibody in the youngest class A young animal can test positive because its mother was infected. Peel and colleagues (2018) fitted bat serology with a model that keeps maternally derived and infection acquired antibody apart, and a simpler joint model of the same kind is one of the three fits compared here. The host is the one from the first section, with a constant rate of `r sprintf("%.2f", foi_true)`, plus a class of 345 juveniles aged half a year. A newborn carries maternal antibody with the probability that its mother is seropositive, taken as the prevalence at age four, and that antibody stops being detectable with a half-life `h`. Prevalence at age `a` is then `p(a) = 1 - exp(-lambda * a) * (1 - pi_m * 2^(-a / h))` The three fits are the catalytic GLM on every class, the same GLM after the usual repair of dropping animals under a year, and a likelihood with the maternal term that estimates `lambda`, `pi_m` and `h` together (from three starting points, keeping the best); the same GLM without animals under two is reported alongside as a stricter version of the repair. The bias of the catalytic fits is arithmetic: fitting them to the expected counts, with no sampling noise at all, gives the value each one converges on as the survey grows. The simulation is needed only for interval coverage and for how often the deviance check notices. ```{r maternal} mat_age <- c(0.5, seq_len(n_class)) mat_n <- c(345, sampled) mat_pim <- 1 - exp(-foi_true * 4) mat_h <- c(0.5, 1) p_maternal <- function(a, lam, pim, h) 1 - exp(-lam * a) * (1 - pim * 2^(-a / h)) pseudo_true <- function(keep, q) { nll <- function(log_l) { p <- 1 - exp(-exp(log_l) * mat_age[keep]) -sum(mat_n[keep] * (q[keep] * log(p) + (1 - q[keep]) * log(1 - p))) } exp(optimize(nll, c(-6, 2))$minimum) / foi_true } mat_ratio <- sapply(mat_h, function(h) { q <- p_maternal(mat_age, foi_true, mat_pim, h) c(naive = pseudo_true(mat_age > 0, q), drop1 = pseudo_true(mat_age >= 1, q), drop2 = pseudo_true(mat_age >= 2, q)) }) fit_maternal <- function(pos) { obs <- data.frame(age = mat_age, pos = pos, neg = mat_n - pos) wald <- function(m) { s <- summary(m)$coefficients[1, 1:2] c(exp(s[1] + c(0, -1.96, 1.96) * s[2]), pchisq(deviance(m), m$df.residual, lower.tail = FALSE)) } m_naive <- glm(cbind(pos, neg) ~ 1 + offset(log(age)), data = obs, family = binomial(link = "cloglog")) m_drop <- glm(cbind(pos, neg) ~ 1 + offset(log(age)), data = obs[obs$age >= 1, ], family = binomial(link = "cloglog")) m_drop2 <- glm(cbind(pos, neg) ~ 1 + offset(log(age)), data = obs[obs$age >= 2, ], family = binomial(link = "cloglog")) nll <- function(par) { p <- p_maternal(mat_age, exp(par[1]), plogis(par[2]), exp(par[3])) p <- pmin(pmax(p, 1e-12), 1 - 1e-12) -sum(pos * log(p) + (mat_n - pos) * log(1 - p)) } starts <- list(c(log(0.2), 0, log(0.5)), c(log(0.2), 1, 0), c(log(0.3), -1, log(0.25))) fits <- lapply(starts, function(s) optim(s, nll, method = "BFGS", hessian = TRUE)) best <- fits[[which.min(vapply(fits, function(f) f$value, 0))]] se_ml <- sqrt(solve(best$hessian)[1, 1]) c(wald(m_naive), wald(m_drop), exp(best$par[1] + c(0, -1.96, 1.96) * se_ml), wald(m_drop2)) } set.seed(2018) mat_demo <- fit_maternal(rbinom(length(mat_age), mat_n, p_maternal(mat_age, foi_true, mat_pim, 1))) n_mat_rep <- 200 mat_sim <- lapply(mat_h, function(h) replicate(n_mat_rep, fit_maternal(rbinom(length(mat_age), mat_n, p_maternal(mat_age, foi_true, mat_pim, h))))) covers <- function(res, i) mean(res[i + 1, ] <= foi_true & res[i + 2, ] >= foi_true) mat_tab <- sapply(mat_sim, function(res) c( cover_naive = covers(res, 1), cover_drop = covers(res, 5), cover_ml = covers(res, 9), cover_drop2 = covers(res, 12), flag_naive = mean(res[4, ] < 0.05), flag_drop = mean(res[8, ] < 0.05), pass_miss = mean(res[8, ] >= 0.05 & (res[6, ] > foi_true | res[7, ] < foi_true)), width_drop = median(res[7, ] - res[6, ]), width_ml = median(res[11, ] - res[10, ]))) mat_p1 <- 1 - exp(-foi_true) mat_excess1 <- sapply(mat_h, function(h) p_maternal(1, foi_true, mat_pim, h)) - mat_p1 mat_widen <- mat_tab["width_ml", ] / mat_tab["width_drop", ] - 1 mat_share <- mat_n[1] / sum(mat_n) mat_duration <- mat_h / log(2) mat_cover_se <- sqrt(0.25 / n_mat_rep) ``` The juveniles make up `r sprintf("%.0f", 100 * mat_share)` per cent of the `r sprintf("%d", sum(mat_n))` animals, and the prevalence in mothers is `r sprintf("%.3f", mat_pim)`. With a half-life of half a year, the catalytic fit on every class converges on `r sprintf("%.3f", mat_ratio["naive", 1])` times the true rate and the fit without juveniles on `r sprintf("%.3f", mat_ratio["drop1", 1])` times; with a half-life of one year the two limits are `r sprintf("%.3f", mat_ratio["naive", 2])` and `r sprintf("%.3f", mat_ratio["drop1", 2])`. Dropping the juveniles does not remove the problem, because the next class is aged exactly one, when a quarter of the maternal antibody is still detectable at `h = 0.5` and half of it at `h = 1`. That adds `r sprintf("%.1f", 100 * mat_excess1[1])` and `r sprintf("%.1f", 100 * mat_excess1[2])` percentage points to a catalytic prevalence of `r sprintf("%.1f", 100 * mat_p1)` per cent at age one, and a single rate has to absorb the excess. Dropping everything under two years brings the limits to `r sprintf("%.3f", mat_ratio["drop2", 1])` and `r sprintf("%.3f", mat_ratio["drop2", 2])`. Over `r sprintf("%d", n_mat_rep)` simulated surveys at each half-life (Monte Carlo standard error at most `r sprintf("%.1f", 100 * mat_cover_se)` points), the interval from the fit on every class covers the true rate in `r sprintf("%.1f", 100 * mat_tab["cover_naive", 1])` and `r sprintf("%.1f", 100 * mat_tab["cover_naive", 2])` per cent of surveys, and the deviance check of the catalytic fit above flags `r sprintf("%.1f", 100 * mat_tab["flag_naive", 1])` and `r sprintf("%.1f", 100 * mat_tab["flag_naive", 2])` per cent of those fits: ignoring maternal antibody altogether is caught. The fit without juveniles covers in `r sprintf("%.1f", 100 * mat_tab["cover_drop", 1])` per cent at `h = 0.5` and `r sprintf("%.1f", 100 * mat_tab["cover_drop", 2])` per cent at `h = 1`. At `h = 1` its own deviance check flags `r sprintf("%.1f", 100 * mat_tab["flag_drop", 2])` per cent of surveys, so in a survey of this size the problem announces itself. At `h = 0.5` the check flags `r sprintf("%.1f", 100 * mat_tab["flag_drop", 1])` per cent, and in `r sprintf("%.1f", 100 * mat_tab["pass_miss", 1])` per cent of all surveys the reduced fit passes its own check while its interval misses the true rate. Dropping everything under two years covers in `r sprintf("%.1f", 100 * mat_tab["cover_drop2", 1])` and `r sprintf("%.1f", 100 * mat_tab["cover_drop2", 2])` per cent. The joint model covers in `r sprintf("%.1f", 100 * mat_tab["cover_ml", 1])` and `r sprintf("%.1f", 100 * mat_tab["cover_ml", 2])` per cent, and its median interval is `r sprintf("%.0f", 100 * mat_widen[1])` and `r sprintf("%.0f", 100 * mat_widen[2])` per cent wider than the interval of the fit without juveniles. In a single survey drawn at `h = 1` the three estimates are `r sprintf("%.3f", mat_demo[1])` (`r sprintf("%.3f", mat_demo[2])` to `r sprintf("%.3f", mat_demo[3])`), `r sprintf("%.3f", mat_demo[5])` (`r sprintf("%.3f", mat_demo[6])` to `r sprintf("%.3f", mat_demo[7])`) and `r sprintf("%.3f", mat_demo[9])` (`r sprintf("%.3f", mat_demo[10])` to `r sprintf("%.3f", mat_demo[11])`) against `r sprintf("%.2f", foi_true)`, and both catalytic fits fail the deviance check with p `r if (max(mat_demo[c(4, 8)]) < 0.001) "below 0.001" else sprintf("up to %.3f", max(mat_demo[c(4, 8)]))`. How long maternal antibody stays detectable decides which of these cases applies. Epstein and colleagues (2013) followed pteropid bat pups and found antibody against canine distemper virus and Hendra virus waning over a mean of about 229 and 255 days. Under the exponential loss assumed here, a half-life of half a year means a mean detectable duration of `r sprintf("%.0f", 365.25 * mat_duration[1])` days, so `h = 0.5` is the realistic case for such hosts and `h = 1` (`r sprintf("%.2f", mat_duration[2])` years) is the long tail. Even the realistic case leaves the fit without juveniles biased in a survey of this size, because the yearlings still carry maternal antibody. Either drop the yearlings as well, which is enough here at `h = 0.5` and not at `h = 1`, or fit the maternal term, which costs some interval width and needs no guess at the half-life. ## Two mechanisms, one curve What the catalytic model leaves out is the honest limit of the whole exercise. Two processes it leaves out flatten an age prevalence curve below one, and a single cross-section of ordinary size cannot separate them reliably. Infection may kill: seropositive animals are removed faster, so the older age classes are enriched with survivors who were never infected. Or antibodies may wane: animals revert to seronegative and are counted as never infected. ```{r confound} p_mortality <- function(a, lam, alpha) { inf <- lam * (exp(-alpha * a) - exp(-lam * a)) / (lam - alpha) sus <- exp(-lam * a) inf / (sus + inf) } p_reversion <- function(a, lam, rho) lam / (lam + rho) * (1 - exp(-(lam + rho) * a)) age_grid <- seq_len(15) foi_mort <- 0.30 alpha_mort <- 0.45 curve_mort <- p_mortality(age_grid, foi_mort, alpha_mort) match_rev <- optim(log(c(0.3, 0.2)), function(par) sum((p_reversion(age_grid, exp(par[1]), exp(par[2])) - curve_mort)^2), method = "BFGS") rev_par <- exp(match_rev$par) curve_rev <- p_reversion(age_grid, rev_par[1], rev_par[2]) max_gap <- max(abs(curve_rev - curve_mort)) half_life <- log(2) / rev_par[2] years_left <- 1 / alpha_mort ``` The best matching seroreversion model has a force of infection of `r sprintf("%.3f", rev_par[1])` and a reversion rate of `r sprintf("%.3f", rev_par[2])`, and its curve differs from the mortality curve by at most `r sprintf("%.1f", 100 * max_gap)` percentage points anywhere between ages one and `r sprintf("%d", max(age_grid))`. The two describe entirely different animals. One says infection removes a host in `r sprintf("%.1f", years_left)` years on average; the other says infection is harmless and antibodies have a half life of `r sprintf("%.1f", half_life)` years. ```{r confoundpower} n_rep <- 200 sampled3 <- c(200, 180, 165, 150, 135, 120, 105, 95, 85, 78, 70, 62, 55, 48, 42) fit_two <- function(pos) { nll <- function(par, fun) { p <- fun(age_grid, exp(par[1]), exp(par[2])) p <- pmin(pmax(p, 1e-9), 1 - 1e-9) -sum(pos * log(p) + (sampled3 - pos) * log(1 - p)) } c(2 * optim(log(c(0.3, 0.4)), nll, fun = p_mortality, method = "BFGS")$value, 2 * optim(log(c(0.3, 0.2)), nll, fun = p_reversion, method = "BFGS")$value) } set.seed(5) two_fit <- replicate(n_rep, fit_two(rbinom(length(age_grid), sampled3, curve_mort))) aic_diff <- two_fit[2, ] - two_fit[1, ] wrong_wins <- mean(aic_diff < 0) n_survey <- sum(sampled3) aic_mean <- mean(aic_diff) aic_sd <- sd(aic_diff) bern_var <- 0.25 share_se <- 100 * sqrt(bern_var / n_rep) ``` With `r sprintf("%d", n_survey)` animals in the survey, generated from the mortality model, the mean AIC advantage of the correct model is `r sprintf("%.1f", aic_mean)` points with a standard deviation of `r sprintf("%.1f", aic_sd)`, and the wrong model has the lower AIC in `r sprintf("%.0f", 100 * wrong_wins)` per cent of the surveys. This is not a power problem that a bigger sample fixes cheaply; it is two mechanisms writing nearly the same curve. ```{r fig-confound} #| fig-cap: "Disease induced mortality and antibody waning, fitted to produce the same age prevalence curve." #| fig-alt: "Two panels. The left panel plots two nearly overlapping curves of prevalence against age, one for extra mortality of the infected and one for antibody waning, rising to a plateau just above sixty per cent. The right panel is a histogram of the AIC difference between the two models across two hundred simulated surveys, centred just above zero with a substantial part of its mass below zero." conf_long <- data.frame( age = rep(age_grid, 2), p = c(curve_mort, curve_rev), mechanism = rep(c("infection kills", "antibodies wane"), each = length(age_grid))) p_conf <- ggplot(conf_long, aes(age, p, colour = mechanism)) + geom_line(linewidth = 1) + geom_point(size = 2) + scale_colour_manual(values = c(te_gold, te_forest), name = NULL) + scale_y_continuous(limits = c(0, 1)) + labs(x = "age (years)", y = "seroprevalence", title = "Same curve, opposite biology") + theme_datasheet() + theme(legend.position = "bottom") p_aic <- ggplot(data.frame(d = aic_diff), aes(d)) + geom_histogram(bins = 30, fill = te_forest, colour = te_paper, linewidth = 0.2) + geom_vline(xintercept = 0, colour = te_rust, linetype = "dashed", linewidth = 0.8) + labs(x = "AIC of waning minus AIC of mortality", y = "surveys", title = "And the data choose poorly", subtitle = "dashed red: no preference") + theme_datasheet() p_conf + p_aic + plot_annotation(theme = theme_datasheet()) ``` ## What to report Give the age distribution of the sample, not only the prevalence. The number of seronegative animals per age class is what carries the information, and a survey that is mostly old animals can have a large sample size and almost no ability to estimate a rate. Fit the catalytic model rather than a logistic regression on age. It is `family = binomial(link = "cloglog")` with `offset(log(age))`, its intercept exponentiates to a rate with units, and it cannot produce a positive prevalence at birth. Report the force of infection with an interval, and report the age bands if the rate was allowed to vary. Band edges are a modelling choice, and the estimate in a band where prevalence has saturated should be quoted with its width, which will be embarrassing and should be. Say which biological assumptions the estimate rests on: no seroreversion, no differential mortality, no maternal antibody, and a population at steady state. Every one of them is a claim about the host. The deviance check can show that the curve is not catalytic, as it did for the fit that ignored maternal antibody in `r sprintf("%.1f", 100 * min(mat_tab["flag_naive", ]))` per cent of the simulated surveys above, but it cannot say which assumption failed, and a fit that passes it is not cleared: at `h = 0.5` the fit without juveniles passed in `r sprintf("%.1f", 100 * (1 - mat_tab["flag_drop", 1]))` per cent of surveys, although it converges on `r sprintf("%.3f", mat_ratio["drop1", 1])` times the true rate. If the sample includes animals young enough to carry maternal antibody, say how they were handled. If mean age at infection or a reproduction number is derived from the rate, give the demographic quantities that went into the derivation as well. Those identities are covered separately on the site and they are only as good as the mortality schedule they assume. ## Honest limits The largest limitation is the one measured above: a cross-section of ordinary size cannot reliably separate disease induced mortality from antibody waning, and both are common in wildlife. Longitudinal recaptures, or a test that distinguishes recent from historical infection, or independent mortality data, are the ways out, and all of them are more expensive than a survey. Every model here assumes the population is at steady state and that the force of infection has been what it is now for the whole lifespan of the oldest animal sampled. A cross-section reads age as calendar time backwards, so a rate that changed five years ago appears as an age effect in exactly the same place as a genuine age effect. The two are not separable from one survey either, and the standard fix, fitting age and cohort together, needs surveys from more than one year. Maternal antibody is measured above for one constant rate, one sample design and an exponential loss of detectability. A seropositive call is a threshold on a decaying titre, so for one animal the loss is a step when its titre crosses the threshold, and across animals it follows the spread of those crossing times, which need not be exponential. The yearlings here are sampled at exactly one year, the youngest age the usual repair keeps and so the age at which the kept animals carry the most maternal antibody; a birth pulse sampled once a year puts the classes at 0.5, 1.5, 2.5 years, where the fit without juveniles is less biased, and neither that design nor the joint model's behaviour in it was simulated. In the band model the maternal positives would land in the first band, and dropping the youngest class there removes data from the band that is otherwise best determined; that version was not simulated. Ages are treated as exact. Field ages usually come from tooth wear or horn rings and carry error that a regression on age propagates in the familiar way, which the site covers under measurement error. The catalytic model is not exempt: its offset takes the recorded ages as they stand, and how far age error moves the estimated rate was not simulated here. The band model uses fixed edges chosen with knowledge of the truth. Choosing edges from the same data that fits the rates makes the intervals optimistic, and the honest alternatives are a smooth on the cloglog scale or edges chosen from the host's biology before the data are seen. The comparison of the two mechanisms uses `r sprintf("%d", n_rep)` simulated surveys, so the reported share of surveys favouring the wrong model has a standard error of about `r sprintf("%.1f", share_se)` percentage points. The conclusion that the two models are nearly indistinguishable does not rest on that share; it rests on the `r sprintf("%.1f", 100 * max_gap)` percentage point maximum gap between the curves, which involves no sampling at all. ## References Heisey DM, Joly DO, Messier F 2006 Ecology 87(9):2356-2365 (10.1890/0012-9658(2006)87[2356:TFOGFM]2.0.CO;2) Hens N, Aerts M, Faes C, Shkedy Z, Lejeune O, Van Damme P, Beutels P 2009 Epidemiology and Infection 138(6):802-812 (10.1017/S0950268809990781) Grenfell BT, Anderson RM 1985 Journal of Hygiene 95(2):419-436 (10.1017/S0022172400062859) Peel AJ, Baker KS, Hayman DTS, Broder CC, Cunningham AA, et al. 2018 Scientific Reports 8:3859 (10.1038/s41598-018-22236-6) Epstein JH, Baker ML, Zambrana-Torrelio C, et al. 2013 PLoS ONE 8(6):e67584 (10.1371/journal.pone.0067584) ## Related tutorials - [Prevalence from an imperfect test](../prevalence-from-an-imperfect-test/) - [Snapshot prevalence and the phase of an epidemic](../snapshot-prevalence-and-epidemic-phase/) - [Density or frequency transmission in enclosures](../density-or-frequency-dependent-transmission/) - [Host density and endemic disease](../host-density-and-endemic-disease/)