Faecal egg count reduction tests at low counts

parasitology

disease ecology

bootstrap

simulation

ecology tutorial

At McMaster egg counts a drug that works 99 per cent usually reads as 100 and its interval collapses. Coverage and resistance calls are measured in R.

Author

Tidy Ecology

Published

2026-09-08

A smallholder with ten ewes wants to know whether the white drench still works. The vet takes a faecal sample from each ewe on the day of treatment and again two weeks later, and the laboratory counts eggs on a McMaster slide. One egg seen under the grid is recorded as fifty eggs per gram of faeces, because the chamber holds a fixed fraction of a diluted sample. The pre-treatment counts come back at a few hundred eggs per gram, the post-treatment counts are all zero, and the report says the drug reduced egg output by 100 per cent. A bootstrap interval, if anyone computes one, runs from 100 to 100.

The faecal egg count reduction test is the standard field test for anthelmintic resistance in sheep, goats, cattle and horses. In its paired form it divides the mean post-treatment count by the mean pre-treatment count and compares the reduction with a threshold. The problem this post measures is well known in veterinary parasitology, and nothing here is a new result. Levecke and colleagues (2012) showed by simulation that the outcome of the test depends on the detection limit of the counting method, the number of animals and the level and aggregation of egg output; Torgerson, Paul and Furrer (2014) built the eggCounts package around a model that treats the slide count, not the eggs per gram, as the observation; and the 2023 guideline of the World Association for the Advancement of Veterinary Parasitology (Kaplan and colleagues 2023) replaced the 1992 rule (Coles and colleagues 1992) with one that sizes the test by the number of eggs counted and asks for positive evidence before calling a drug either susceptible or resistant. What the post adds is a small, runnable demonstration: how often the report reads 100, how often its interval covers the truth, and how often a resistant drug is called, under both rules.

Three posts on this site deal with measurements that stop resolving at the low end, and none of them is this problem. Values below the detection limit treats a continuous concentration that is censored below a limit and fits the censored likelihood; its test is a trend regression, and nothing in it is a ratio or a treatment decision. Rounded and coarsened measurements works with a known grid that adds a fixed amount of variance to a continuous value. Parasite burdens truncated by host death is about worm counts and the aggregation parameter. Here the limit is a counting multiplier applied to aggregated counts, and it acts on a ratio at its top end, which is exactly where the treatment decision sits. The failure of the bootstrap that follows is a cousin of the one in checking a bootstrap, where resampling cannot produce a value beyond the sample maximum: here it cannot produce a post-treatment egg when none was counted.

library(ggplot2)

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),
          strip.text       = element_text(colour = te_ink))
}

The generating model and three intervals

Each animal carries its own true egg output, drawn from a gamma distribution with mean 150 eggs per gram and shape 0.7, so that egg output is strongly aggregated: a few animals carry most of the eggs. Egg output also varies from day to day within an animal, with a coefficient of variation of 0.3, drawn independently for the two sampling days. The treatment multiplies the animal’s output by one minus the efficacy. The slide sees a Poisson number of eggs with mean equal to the output divided by the multiplication factor, and the laboratory reports that count times the factor. All of these constants were fixed before anything was run.

cv_day   <- 0.3                 # day-to-day excretion CV within an animal
shp_day  <- 1 / cv_day^2
k_set    <- 0.7                 # gamma shape of egg output among animals
mu_set   <- 150                 # mean eggs per gram before treatment
n_set    <- 10                  # animals in the treatment group
mf_set   <- 50                  # McMaster multiplication factor
eff_set  <- 0.99                # true efficacy
n_rep    <- 600                 # simulated tests per design cell
n_boot   <- 999                 # paired bootstrap resamples per test
mcse_max <- sqrt(0.25 / n_rep)

sim_test <- function(R, n, mu, k, mf, eff) {
  epg  <- matrix(rgamma(R * n, shape = k, rate = k / mu), R)
  pre  <- matrix(rpois(R * n, epg * rgamma(R * n, shp_day, shp_day) / mf), R)
  post <- matrix(rpois(R * n, epg * (1 - eff) * rgamma(R * n, shp_day, shp_day) / mf), R)
  list(pre = pre, post = post, mf = mf)   # raw eggs seen on the slide
}

fecr_of <- function(tst) 100 * (1 - rowSums(tst$post) / rowSums(tst$pre))

The reduction is computed from group arithmetic means, as in the 1992 guideline; the 2023 guideline moves the decision onto the interval instead of a prescribed mean. The multiplication factor cancels from that ratio, so the point estimate can be computed from raw slide counts. The intervals are where the choices differ, and three are compared.

The first is a paired bootstrap: resample animals with replacement, keeping each animal’s two counts together, and take percentiles of the resampled reduction. Resampling animals is the same as drawing multinomial weights for them, which lets all resamples of all simulated tests be computed with two matrix products.

boot_lim <- function(tst, probs, B = n_boot) {
  n_an  <- ncol(tst$pre)
  wts   <- rmultinom(B, n_an, rep(1 / n_an, n_an))
  bpre  <- tst$pre %*% wts
  bpost <- tst$post %*% wts
  bf <- 100 * (1 - bpost / bpre)
  bf[!is.finite(bf)] <- NA          # resamples with no pre-treatment eggs
  t(apply(bf, 1, quantile, probs = probs, na.rm = TRUE, names = FALSE))
}

The second uses the counts. Under the model in the eggCounts package, an animal’s pre- and post-treatment slide counts are Poisson with means proportional to the same gamma-distributed output, the second multiplied by one minus the efficacy. Conditional on the animal’s total count, its post-treatment count is then binomial with probability r over one plus r, where r is one minus the efficacy, and the animal’s own output drops out. The profile likelihood for the efficacy is therefore a binomial likelihood in the share of all counted eggs that were counted after treatment, whatever the aggregation. The interval below inverts that likelihood ratio. When no post-treatment egg is seen its lower limit has a closed form that depends only on the number of eggs counted, which is the quantity the bootstrap never sees.

The day-to-day variation breaks the binomial step, because the two days no longer share one mean. The third interval is the quasi-likelihood version of the second: it estimates a dispersion from the Pearson residuals across animals, divides the likelihood ratio by it when it exceeds one, and uses an F critical value.

to_eff <- function(p_share) 100 * (1 - p_share / (1 - p_share))

lik_lim <- function(tst, level, quasi = FALSE) {
  z_post <- rowSums(tst$post)
  n_all  <- rowSums(tst$pre) + z_post
  p_hat  <- z_post / n_all
  loglik <- function(p) ifelse(z_post > 0, z_post * log(p), 0) +
    ifelse(n_all - z_post > 0, (n_all - z_post) * log1p(-p), 0)
  l_max <- loglik(p_hat)
  phi  <- rep(1, length(z_post))
  crit <- rep(qchisq(level, 1), length(z_post))
  if (quasi) {
    tot_an <- tst$pre + tst$post
    pear <- ifelse(tot_an > 0, (tst$post - tot_an * p_hat)^2 /
                     (tot_an * p_hat * (1 - p_hat)), 0)
    df_an <- pmax(rowSums(tot_an > 0) - 1, 1)
    phi_raw <- rowSums(pear) / df_an
    over <- z_post > 0 & z_post < n_all & phi_raw > 1
    phi[over]  <- phi_raw[over]
    crit[over] <- qf(level, 1, df_an[over])
  }
  gap <- function(p) 2 * (l_max - loglik(p)) / phi - crit
  bisect <- function(lo, hi, upper_side) {
    for (it in 1:60) {
      mid <- (lo + hi) / 2
      inside <- gap(mid) < 0
      if (upper_side) {
        lo <- ifelse(inside, mid, lo); hi <- ifelse(inside, hi, mid)
      } else {
        hi <- ifelse(inside, mid, hi); lo <- ifelse(inside, lo, mid)
      }
    }
    (lo + hi) / 2
  }
  p_up <- bisect(p_hat, rep(1 - 1e-12, length(p_hat)), TRUE)
  p_lo <- ifelse(z_post == 0, 0, bisect(rep(1e-300, length(p_hat)), p_hat, FALSE))
  cbind(lo = to_eff(p_up), hi = to_eff(p_lo))
}

Ten ewes and a factor of fifty

set.seed(4417)
ex_draws <- sim_test(20, n_set, mu_set, k_set, mf_set, eff_set)
all_zero <- which(rowSums(ex_draws$post) == 0 & rowSums(ex_draws$pre) > 0)
i_ex <- all_zero[1]
ex <- list(pre  = ex_draws$pre[i_ex, , drop = FALSE],
           post = ex_draws$post[i_ex, , drop = FALSE], mf = mf_set)
ex_epg_pre  <- mf_set * ex$pre
ex_mean_pre <- mean(ex_epg_pre)
ex_eggs     <- sum(ex$pre)
ex_zero_pre <- sum(ex$pre == 0)
ex_max_epg  <- max(ex_epg_pre)
ex_fecr     <- fecr_of(ex)
ex_boot     <- boot_lim(ex, c(0.025, 0.975))
ex_lik      <- lik_lim(ex, 0.95)
ex_lo_closed <- to_eff(1 - exp(-qchisq(0.95, 1) / (2 * ex_eggs)))
ex_epg_pre

     [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
[1,]  100   50  250  450  100   50    0  150    0    50

This is the first of twenty simulated tests from the design above in which no egg was seen after treatment, and it was taken as the first such draw rather than chosen. The ten ewes had a mean pre-treatment count of 120 eggs per gram, the largest 450, and 2 of them showed no egg at all before treatment. Behind that mean are 24 eggs actually counted under the microscope. The reported reduction is 100 per cent against a true 99.

The paired bootstrap interval runs from 100.0 to 100.0. Every resample contains only zeros after treatment, so every resampled reduction is 100, and the interval has no width at all. The likelihood interval runs from 91.7 to 100.0. Its lower limit has a closed form. With no post-treatment egg among 24 counted, the likelihood ratio reaches the chi-squared critical value at a post-treatment share of one minus the exponential of minus 1.92 over 24, which is an efficacy of 91.7, and the interval says what the slide supports: a drug somewhere above about nine tenths effective, with no evidence either way about 95 or 99 per cent.

How often a 99 per cent drug reads as 100

The report of exactly 100 per cent happens when every post-treatment slide is empty, and that probability has a closed form. If every animal excreted the same number of eggs per gram, the post-treatment total would be Poisson with mean n times the output times one minus the efficacy, divided by the factor, and the chance of seeing nothing is the exponential of minus that. With gamma-distributed output among animals each slide is empty with the negative binomial zero probability, and the group probability is its n-th power. The day-to-day variation is not in either formula, so the simulation checks how much it matters.

mf_grid  <- c(50, 25, 15, 1)
n_grid   <- c(10, 20, 40)
eff_grid <- c(0.99, 0.95, 0.90)
n_zero_rep <- 20000

p_zero_fixed <- function(n, mu, mf, eff) exp(-n * mu * (1 - eff) / mf)
p_zero_gamma <- function(n, mu, k, mf, eff) (1 + mu * (1 - eff) / (k * mf))^(-k * n)
pz_fixed_set <- p_zero_fixed(n_set, mu_set, mf_set, eff_set)
pz_gamma_set <- p_zero_gamma(n_set, mu_set, k_set, mf_set, eff_set)

set.seed(2023)
zero_tab <- do.call(rbind, lapply(n_grid, function(n_an) {
  do.call(rbind, lapply(mf_grid, function(mf) {
    tst <- sim_test(n_zero_rep, n_an, mu_set, k_set, mf, eff_set)
    ok  <- rowSums(tst$pre) > 0
    data.frame(n = n_an, mf = mf,
               sim   = mean(rowSums(tst$post)[ok] == 0),
               fixed = p_zero_fixed(n_an, mu_set, mf, eff_set),
               gamma = p_zero_gamma(n_an, mu_set, k_set, mf, eff_set))
  }))
}))
zero_gap_gamma <- max(abs(zero_tab$sim - zero_tab$gamma))
zero_gap_fixed <- max(abs(zero_tab$sim - zero_tab$fixed))
mcse_zero <- sqrt(0.25 / n_zero_rep)
zt <- function(n_an, mf, col) zero_tab[zero_tab$n == n_an & zero_tab$mf == mf, col]

For ten animals at 150 eggs per gram and a factor of 50, the fixed-output formula gives 0.741 and the gamma formula 0.745. Simulating 20000 tests per design gives 0.749. Across all twelve designs the gamma formula stays within 0.0048 of the simulation, against a Monte Carlo standard error of at most 0.0035, so day-to-day variation barely moves this probability; the fixed-output formula is off by up to 0.030, because aggregation concentrates the eggs in fewer animals and leaves more slides empty.

The exponent is the whole story: n times the mean count times one minus the efficacy, over the factor, is the expected number of eggs on all post-treatment slides together. At a factor of 50 that is 0.3 eggs for ten animals. Doubling the animals and switching to a factor of 15 raises it to 2.0, and the share of empty post-treatment totals falls to 0.159. Only the method that counts every egg (factor 1) removes the problem outright.

zero_tab$group <- factor(sprintf("%d animals", zero_tab$n),
                         levels = sprintf("%d animals", n_grid))
mf_fine <- exp(seq(log(1), log(50), length.out = 120))
zero_curve <- do.call(rbind, lapply(n_grid, function(n_an)
  data.frame(mf = mf_fine, gamma = p_zero_gamma(n_an, mu_set, k_set, mf_fine, eff_set),
             group = factor(sprintf("%d animals", n_an), levels = levels(zero_tab$group)))))
ggplot(zero_tab, aes(mf, sim, colour = group)) +
  geom_line(data = zero_curve, aes(y = gamma), linewidth = 0.9) +
  geom_point(size = 2.4) +
  scale_x_log10(breaks = c(1, 2, 5, 15, 25, 50)) +
  scale_colour_manual(values = c(te_rust, te_gold, te_forest), name = NULL) +
  scale_y_continuous(limits = c(0, 1)) +
  labs(x = "multiplication factor (eggs per gram per egg seen)",
       y = "share reporting exactly 100 per cent",
       title = "A perfect result is mostly a counting result",
       subtitle = "true efficacy 99 per cent, mean 150 eggs per gram, shape 0.7") +
  theme_datasheet() + theme(legend.position = "bottom")

Three rising curves on warm off-white paper, one per group size, plotting the share of tests that report exactly 100 per cent against the multiplication factor on a log axis from 1 to 50. All three start at zero at a factor of 1 and stay near zero until a factor of about 3 for ten animals, 5 for twenty and 10 for forty. At a factor of 50 the red ten-animal curve reaches about 0.75, the gold twenty-animal curve about 0.56 and the green forty-animal curve about 0.31. Simulated points sit on the closed-form curves at factors 1, 15, 25 and 50. — Figure 1: Share of simulated tests of a 99 per cent effective drug that report exactly 100 per cent reduction, by multiplication factor and group size. Lines: gamma closed form; points: simulation.

The bootstrap interval collapses where the counts do

The share reporting 100 is arithmetic. The interval is not, so it is measured. Each design cell below simulates 600 tests, a number fixed before running, which puts the Monte Carlo standard error of any coverage at 0.020 or less. Tests with no egg at all before treatment are dropped, as a laboratory would drop them.

run_cell <- function(n, mu, k, mf, eff, R = n_rep) {
  tst <- sim_test(R, n, mu, k, mf, eff)
  ok  <- rowSums(tst$pre) > 0
  tst$pre  <- tst$pre[ok, , drop = FALSE]
  tst$post <- tst$post[ok, , drop = FALSE]
  fecr <- fecr_of(tst)
  bq  <- boot_lim(tst, c(0.025, 0.05, 0.95, 0.975))
  l95 <- lik_lim(tst, 0.95);               l90 <- lik_lim(tst, 0.90)
  q95 <- lik_lim(tst, 0.95, quasi = TRUE); q90 <- lik_lim(tst, 0.90, quasi = TRUE)
  te <- 100 * eff
  covers  <- function(lo, hi) mean(lo <= te & hi >= te)
  class23 <- function(lo, hi) c(res = mean(hi < 99), sus = mean(lo >= 95 & hi >= 99))
  c23b <- class23(bq[, 2], bq[, 3])
  c23l <- class23(l90[, 1], l90[, 2])
  c23q <- class23(q90[, 1], q90[, 2])
  data.frame(n = n, mu = mu, k = k, mf = mf, eff = eff, used = sum(ok),
             eggs = n * mu / mf, rep100 = mean(fecr == 100),
             cov_boot = covers(bq[, 1], bq[, 4]), cov_lik = covers(l95[, 1], l95[, 2]),
             cov_quasi = covers(q95[, 1], q95[, 2]),
             r92_boot = mean(fecr < 95 & bq[, 1] < 90),
             r92_lik  = mean(fecr < 95 & l95[, 1] < 90),
             res_boot = c23b[["res"]], sus_boot = c23b[["sus"]],
             res_lik = c23l[["res"]], sus_lik = c23l[["sus"]],
             res_quasi = c23q[["res"]], sus_quasi = c23q[["sus"]])
}

set.seed(1992)
grid_tab <- do.call(rbind, lapply(eff_grid, function(e) do.call(rbind,
  lapply(n_grid, function(n_an) do.call(rbind,
    lapply(mf_grid, function(mf) run_cell(n_an, mu_set, k_set, mf, e)))))))

gv <- function(n_an, mf, e, col) grid_tab[grid_tab$n == n_an & grid_tab$mf == mf &
                                             abs(grid_tab$eff - e) < 1e-9, col]
n_dropped <- sum(n_rep - grid_tab$used)
top <- grid_tab[abs(grid_tab$eff - 0.99) < 1e-9, ]
boot_low_eggs <- range(top$cov_boot[top$eggs <= 60])
lik_mf1  <- range(grid_tab$cov_lik[grid_tab$mf == 1])
quasi_all <- range(grid_tab$cov_quasi)
quasi_mf1 <- range(grid_tab$cov_quasi[grid_tab$mf == 1])
quasi_rest_min <- min(grid_tab$cov_quasi[grid_tab$mf != 1])
boot_90 <- range(grid_tab$cov_boot[abs(grid_tab$eff - 0.90) < 1e-9])

At ten animals and a factor of 50 the 99 per cent drug reports exactly 100 in 0.758 of tests, and the paired bootstrap interval covers the true 99 in 0.242, against a nominal 0.95 (0 of the 21600 simulated tests in the grid were dropped for an empty pre-treatment slide). Where the grid expects 60 or fewer eggs to be counted before treatment, its coverage for this drug is between 0.242 and 0.457. Twenty animals at a factor of 15 lift the coverage to 0.813, and forty animals at a factor of 1 give 0.933.

The likelihood interval, on exactly the same simulated counts, covers the true 99 in 0.982 of tests at ten animals and a factor of 50. It does not need more counting to be honest at the top end, because a count of zero out of a known number of eggs is information it can use. It fails in the opposite corner. When every egg is counted, day-to-day variation in excretion is larger than the Poisson counting noise the binomial assumes, and the plain likelihood interval covers between 0.543 and 0.892 across the factor-of-1 designs, worst for the less effective drugs. The quasi-likelihood interval, which estimates that extra variation from the animals, covers between 0.842 and 0.988 across all thirty-six cells; every cell in which it falls below 0.913 is a factor-of-1 design (those range from 0.842 to 0.948), where a dispersion estimated from ten to forty aggregated animals is itself noisy.

cov_long <- do.call(rbind, lapply(
  list(c("cov_boot", "paired bootstrap"), c("cov_lik", "binomial likelihood"),
       c("cov_quasi", "quasi-likelihood")), function(m)
    data.frame(eggs = grid_tab$eggs, eff = grid_tab$eff, cover = grid_tab[[m[1]]],
               interval = m[2])))
cov_long$interval <- factor(cov_long$interval,
                            levels = c("paired bootstrap", "binomial likelihood", "quasi-likelihood"))
cov_long$panel <- factor(sprintf("true efficacy %.0f per cent", 100 * cov_long$eff),
                         levels = sprintf("true efficacy %.0f per cent", 100 * eff_grid))
cov_mean <- aggregate(cover ~ eggs + panel + interval, data = cov_long, FUN = mean)

ggplot(cov_long, aes(eggs, cover, colour = interval)) +
  geom_hline(yintercept = 0.95, linetype = "dashed", colour = te_body, linewidth = 0.5) +
  geom_line(data = cov_mean, linewidth = 0.8) +
  geom_point(size = 2, alpha = 0.9) +
  facet_wrap(~ panel, ncol = 3) +
  scale_x_log10(breaks = c(30, 100, 300, 1000, 3000)) +
  scale_colour_manual(values = c(te_rust, te_gold, te_forest), name = NULL) +
  scale_y_continuous(limits = c(0, 1)) +
  labs(x = "eggs expected on the pre-treatment slides (log scale)",
       y = "coverage of the true efficacy",
       title = "Each interval fails at a different end",
       subtitle = "dashed line: nominal 0.95; lines join cell means at equal egg totals") +
  theme_datasheet() + theme(legend.position = "bottom", panel.spacing.x = unit(1.2, "lines"))

Three panels for true efficacies of 99, 95 and 90 per cent, each plotting interval coverage against the expected number of pre-treatment eggs on a log axis from 30 to 6000, with a dashed line at 0.95. In the 99 per cent panel the red paired bootstrap line starts at about 0.24 at 30 eggs and climbs steeply to about 0.9 by 400 eggs, while the gold binomial likelihood and green quasi-likelihood lines stay near 0.95 to 0.99 at low counts. In the 95 and 90 per cent panels the bootstrap starts at about 0.73 and 0.86 and stays below the dashed line; the gold likelihood line holds near 0.92 up to 400 eggs and then drops to about 0.7 and about 0.55 above 1000 eggs; the green quasi-likelihood line stays between about 0.84 and 0.97 throughout. — Figure 2: Coverage of three 95 per cent intervals against the number of eggs expected on the pre-treatment slides (animals times mean count over the factor), for three true efficacies. Each point is one design cell of 600 simulated tests.

For the effective drug the horizontal axis carries most of the pattern. Ten animals at a factor of 25 and twenty animals at a factor of 50 both expect 60 eggs on the pre-treatment slides, and their bootstrap coverages for the 99 per cent drug are 0.418 and 0.457. For the 90 per cent drug the bootstrap is at or below nominal across the whole grid, between 0.838 and 0.945, but never collapses, because post-treatment eggs are seen in most tests and the resamples have something to vary.

Calling resistance under the 1992 and 2023 rules

The 1992 rule calls resistance when the reduction is below 95 per cent and the lower limit of the 95 per cent interval is below 90. The original rule was written for a treated group against an untreated control, with a variance formula for the interval; here it is applied to the paired design with the bootstrap and likelihood intervals. The 2023 guideline recommends the paired design, uses a 90 per cent interval, and compares it with two numbers for the host and drug: an expected efficacy, 99 per cent in the guideline’s sheep example, and a lower efficacy threshold, taken here as 95 per cent (the research protocol’s grey zone of 95 to 99; the clinical protocol uses 90 to 99). A test is resistant when the upper limit is below the expected efficacy, susceptible when the lower limit reaches the threshold and the upper limit reaches the expected efficacy, and inconclusive otherwise. The guideline suggests a Bayesian approach built for these data and recommends one of two web tools for the analysis (the Bayesian eggCounts interface, or fecrt.com with a hybrid frequentist and Bayesian method) unless a qualified statistician does it, and warns that only a few formulae give an interval when the observed reduction is 100 per cent; the bootstrap below is what a laboratory gets if it ignores that advice, and the likelihood intervals are a frequentist stand-in for the model-based analysis, not the analysis itself. The guideline also sizes the treatment group by the total number of eggs counted before the multiplication factor is applied, not by a minimum count per animal.

c92_boot_10 <- gv(10, 50, 0.90, "r92_boot"); c92_boot_20 <- gv(20, 15, 0.90, "r92_boot")
c92_lik_10  <- gv(10, 50, 0.90, "r92_lik")
fp92_max    <- max(grid_tab$r92_boot[abs(grid_tab$eff - 0.99) < 1e-9])

# the headline cell gets its own, larger run so its rate carries a small MC error
n_rep_esc <- 5000
set.seed(1992 + 10)
esc_tst <- sim_test(n_rep_esc, n_set, mu_set, k_set, mf_set, 0.90)
esc_ok  <- rowSums(esc_tst$pre) > 0
esc_tst$pre <- esc_tst$pre[esc_ok, ]; esc_tst$post <- esc_tst$post[esc_ok, ]
esc_fecr <- fecr_of(esc_tst)
esc_call_boot <- mean(esc_fecr < 95 & boot_lim(esc_tst, c(0.025, 0.975))[, 1] < 90)
esc_call_lik  <- mean(esc_fecr < 95 & lik_lim(esc_tst, 0.95)[, "lo"] < 90)
escape_10   <- 1 - esc_call_boot
mcse_grid_esc <- sqrt(c92_boot_10 * (1 - c92_boot_10) / gv(10, 50, 0.90, "used"))
mcse_escape <- sqrt(escape_10 * (1 - escape_10) / sum(esc_ok))
esc_hi_fecr <- mean(esc_fecr >= 95)

designs <- data.frame(n = c(10, 20, 40), mf = c(50, 15, 1),
                      lab = c("10 animals, factor 50", "20 animals, factor 15",
                              "40 animals, factor 1"))
set.seed(1203)
miss_tst <- sim_test(n_rep, n_set, mu_set, k_set, mf_set, 0.90)
miss_ok  <- rowSums(miss_tst$pre) > 0
miss_tst$pre <- miss_tst$pre[miss_ok, ]; miss_tst$post <- miss_tst$post[miss_ok, ]
miss_up  <- boot_lim(miss_tst, c(0.05, 0.95))[, 2]
missed   <- miss_up >= 99
miss_at100 <- mean(miss_up[missed] == 100)
miss_few   <- mean(rowSums(miss_tst$post[missed, ] > 0) <= 2)

class_long <- do.call(rbind, lapply(seq_len(nrow(designs)), function(i) {
  do.call(rbind, lapply(eff_grid, function(e) {
    do.call(rbind, lapply(list(c("boot", "bootstrap"), c("lik", "likelihood"),
                               c("quasi", "quasi")), function(m) {
      res <- gv(designs$n[i], designs$mf[i], e, paste0("res_", m[1]))
      sus <- gv(designs$n[i], designs$mf[i], e, paste0("sus_", m[1]))
      data.frame(design = designs$lab[i], eff = e, interval = m[2],
                 class = c("resistant", "inconclusive", "susceptible"),
                 share = c(res, 1 - res - sus, sus))
    }))
  }))
}))

With ten animals and a factor of 50, a drug that removes only 90 per cent of egg output is called resistant by the 1992 rule with the bootstrap in 0.805 of the grid’s 600 tests. Because this is the headline cell, it was rerun on its own with 5000 tests: the drug escapes a resistance call in 0.236 of them (Monte Carlo standard error 0.006), a little under one test in four (the grid’s run put the escape at 0.195, a low draw given its own standard error of 0.016). The likelihood interval barely changes this, calling resistance in 0.773 against 0.764 for the bootstrap on the same counts, because the escapes come from the point estimate: the rule needs the reduction itself below 95, and with 30 eggs expected on the pre-treatment slides the 90 per cent drug shows a reduction of 95 or more in 0.227 of tests. Twenty animals at a factor of 15 raise the call rate to 0.937. The 1992 rule never calls a 99 per cent drug resistant in more than 0.060 of tests in this grid.

The 2023 rule fed a bootstrap interval fails at low counts in both directions. For ten animals at a factor of 50, it calls the 99 per cent drug susceptible in 0.793 of tests, on the strength of an interval that has collapsed onto 100, and calls the 90 per cent drug resistant in only 0.385. It even calls a 95 per cent drug, which is below the expected efficacy, susceptible in 0.273. The likelihood interval on the same counts calls the 90 per cent drug resistant in 0.778, the 99 per cent drug susceptible in only 0.377, and the 95 per cent drug susceptible in 0.082; its most common answer for the effective drug is inconclusive, in 0.587 of tests, which is what 30 expected eggs can support.

The missed resistance calls for the 90 per cent drug have a plain cause. In a separate run of 600 such tests, 380 were not called resistant; among those, the share with an upper bootstrap limit of exactly 100 was 1.000, and the share whose post-treatment eggs came from two animals or fewer was 1.000. Resamples that leave those animals out read 100, and there are enough of them to hold the upper limit there.

class_long$class <- factor(class_long$class, levels = c("resistant", "inconclusive", "susceptible"))
class_long$interval <- factor(class_long$interval, levels = c("bootstrap", "likelihood", "quasi"))
class_long$design <- factor(class_long$design, levels = designs$lab)
class_long$effl <- factor(sprintf("true %.0f%%", 100 * class_long$eff),
                          levels = sprintf("true %.0f%%", 100 * eff_grid))
ggplot(class_long, aes(interval, share, fill = class)) +
  geom_col(width = 0.75, colour = te_paper, linewidth = 0.3) +
  facet_grid(effl ~ design) +
  scale_fill_manual(values = c(te_rust, te_line, te_forest), name = NULL) +
  scale_y_continuous(breaks = c(0, 0.5, 1)) +
  labs(x = NULL, y = "share of simulated tests",
       title = "A bootstrap interval misleads the 2023 rule",
       subtitle = "rows: true efficacy; columns: design") +
  theme_datasheet() +
  theme(legend.position = "bottom", axis.text.x = element_text(size = 9))

A grid of stacked bars, rows for true efficacies of 99, 95 and 90 per cent and columns for ten animals at factor 50, twenty at factor 15 and forty at factor 1, with one bar each for bootstrap, likelihood and quasi intervals split into red resistant, grey inconclusive and green susceptible shares. In the ten-animal column the bootstrap bar for the 99 per cent drug is about four fifths green, while the likelihood and quasi bars are under two fifths green and mostly grey; for the 90 per cent drug the bootstrap bar is under two fifths red and mostly grey, against about four fifths and two thirds red for likelihood and quasi. In the other two columns the 99 per cent row is mostly green and the 95 and 90 per cent rows are mostly or entirely red for all three intervals. — Figure 3: Classification under the 2023 rule (90 per cent interval, expected efficacy 99, threshold 95) for three designs, three true efficacies and three intervals.

With forty animals counted at a factor of 1 the three intervals agree on the classification almost everywhere. The one visible difference there is that the plain likelihood interval calls the 99 per cent drug resistant in 0.095 of tests against 0.045 for the bootstrap, which is its undercoverage at high counts turning into false alarms.

Aggregation matters less than the eggs counted

The negative binomial shape and the mean count were both held fixed above. The grid below varies them at twenty animals and a factor of 50, for the 99 and 90 per cent drugs.

k_grid  <- c(0.3, 0.7, 2)
mu_grid <- c(50, 150, 600)
set.seed(3009)
agg_tab <- do.call(rbind, lapply(c(0.99, 0.90), function(e) do.call(rbind,
  lapply(k_grid, function(kk) do.call(rbind,
    lapply(mu_grid, function(mu) run_cell(20, mu, kk, mf_set, e)))))))
av <- function(kk, mu, e, col) agg_tab[agg_tab$k == kk & agg_tab$mu == mu &
                                         abs(agg_tab$eff - e) < 1e-9, col]
agg_top <- agg_tab[abs(agg_tab$eff - 0.99) < 1e-9, ]
spread_k_top <- max(tapply(agg_top$cov_boot, agg_top$mu, function(v) diff(range(v))))
spread_mu_top <- max(tapply(agg_top$cov_boot, agg_top$k, function(v) diff(range(v))))
agg_low <- agg_tab[abs(agg_tab$eff - 0.90) < 1e-9, ]
spread_k_low  <- max(tapply(agg_low$cov_boot, agg_low$mu, function(v) diff(range(v))))
spread_mu_low <- max(tapply(agg_low$cov_boot, agg_low$k, function(v) diff(range(v))))

For the 99 per cent drug, changing the shape from 0.3 to 2 at a fixed mean moves the bootstrap coverage by at most 0.118, while changing the mean from 50 to 600 at a fixed shape moves it by up to 0.697. At a mean of 50 eggs per gram, twenty animals expect twenty eggs on their slides, and the bootstrap covers the true 99 in 0.205 of tests at shape 0.7. The likelihood interval covers in 0.982 there.

For the 90 per cent drug the two effects are of similar size, because eggs are seen after treatment and the bootstrap is resampling a skewed ratio rather than a column of zeros: the shape moves its coverage by up to 0.107 and the mean by up to 0.142. At a mean of 50, its coverage is 0.747 at shape 0.3 and 0.853 at shape 2; the quasi-likelihood interval covers 0.942 and 0.943.

agg_long <- rbind(
  data.frame(agg_tab[, c("k", "mu", "eff")], cover = agg_tab$cov_boot, interval = "paired bootstrap"),
  data.frame(agg_tab[, c("k", "mu", "eff")], cover = agg_tab$cov_quasi, interval = "quasi-likelihood"))
agg_long$shape <- factor(sprintf("shape %.1f", agg_long$k), levels = sprintf("shape %.1f", k_grid))
agg_long$panel <- factor(sprintf("true efficacy %.0f per cent", 100 * agg_long$eff),
                         levels = sprintf("true efficacy %.0f per cent", c(99, 90)))
ggplot(agg_long, aes(mu, cover, colour = shape, linetype = interval)) +
  geom_hline(yintercept = 0.95, linetype = "dashed", colour = te_body, linewidth = 0.5) +
  geom_line(linewidth = 0.8) + geom_point(size = 2) +
  facet_wrap(~ panel) +
  scale_x_log10(breaks = mu_grid) +
  scale_colour_manual(values = c(te_rust, te_gold, te_forest), name = NULL) +
  scale_linetype_manual(values = c("solid", "longdash"), name = NULL) +
  scale_y_continuous(limits = c(0, 1)) +
  guides(colour = guide_legend(order = 1), linetype = guide_legend(order = 2, nrow = 2)) +
  labs(x = "mean eggs per gram before treatment (log scale)", y = "coverage of the true efficacy",
       title = "The mean count sets the failure",
       subtitle = "twenty animals, factor 50; dashed line: nominal 0.95") +
  theme_datasheet() + theme(legend.position = "bottom", legend.key.width = unit(2.2, "lines"),
                          panel.spacing.x = unit(2, "lines"))

Two panels for true efficacies of 99 and 90 per cent plotting coverage against mean eggs per gram at 50, 150 and 600, with red, gold and green lines for shapes 0.3, 0.7 and 2, solid for the paired bootstrap and long-dashed for the quasi-likelihood interval, and a dashed reference at 0.95. In the 99 per cent panel the three solid bootstrap lines rise together from about 0.16 to 0.2 at a mean of 50, through about 0.4 to 0.45 at 150, to between about 0.77 and 0.89 at 600, while the dashed lines sit near 0.99 falling to about 0.92. In the 90 per cent panel the solid lines start between about 0.75 and 0.86 at 50 and flatten near 0.9, and the dashed lines stay between about 0.92 and 0.97. — Figure 4: Coverage of the paired bootstrap and quasi-likelihood intervals at twenty animals and a factor of 50, against the mean pre-treatment count, for three levels of aggregation.

What to report

Report the raw number of eggs counted before and after treatment, as well as the multiplication factor. Those two totals are the information in the test; the eggs per gram figures are a rescaling of them, and a reduction of 100 per cent from 30 eggs means something quite different from 100 per cent from 3000.

Do not report a bootstrap interval of zero width. When every post-treatment count is zero the paired bootstrap has nothing to resample and returns the point estimate twice; quoting it as an interval invites the reader to treat a counting result as a precise one. A likelihood interval on the counts, whose lower limit depends on the eggs counted, is a defensible replacement at low counts, and the model-based Bayesian intervals in eggCounts are built on the same observation; analysing the counts with a method built for these data is also the guideline’s own advice.

Say which rule was applied and with which interval. Under the 2023 rule the interval is the whole decision, and the same counts produced a susceptible call for the effective drug in 0.793 of tests with one interval and 0.377 with another. An inconclusive result from a small, coarse test is a correct result.

Plan the test on eggs counted rather than on animals. For the effective drug, the two designs here that expect the same 60 pre-treatment eggs from different group sizes and factors gave bootstrap coverages of 0.418 and 0.457, and the mean count moved coverage far more than aggregation did; that is the same logic as the 2023 guideline sizing groups by counted eggs.

Honest limits

The simulation uses a single gamma-Poisson generator with independent day-to-day variation of fixed size. Real post-treatment counts also carry zero inflation (animals with no patent infection), a change in the proportion of worm species that survive treatment, and the gap between egg output and worm burden that a reduction in egg output does not measure. The binomial likelihood interval is exact only under the model without day-to-day variation, and its failure at high counts was shown here with one value of that variation; a smaller one would shrink the failure and a larger one would widen it.

The Bayesian models in eggCounts, and other model-based intervals a laboratory might use, were not fitted. Their behaviour at low counts depends on the priors, and nothing above says how close the likelihood intervals here come to them. The comparison is between the paired bootstrap, which is easy to compute, and two frequentist intervals built on the counts.

The 1992 rule was applied to a paired design with a resampling interval rather than to the treated-versus-control design and variance formula it was written for, and the 2023 thresholds were taken from the guideline’s sheep example with the research protocol; values for other hosts and drugs may differ. The clinical protocol’s lower threshold of 90 per cent would call more drugs susceptible; the resistant call, which depends only on the upper limit and the expected efficacy of 99, would not change. The low resistance subclass of the 2023 guideline was not separated out.

The scene is a small flock of sheep, and the numbers transfer to cattle, goats and horses only through the eggs counted. Wildlife egg count surveys share the low counts and the coarse factors but rarely have a paired treatment design, and the reduction test does not apply to them as it stands.

References

Coles GC, Bauer C, Borgsteede FHM, Geerts S, Klei TR, Taylor MA, Waller PJ 1992 Veterinary Parasitology 44(1-2):35-44 (10.1016/0304-4017(92)90141-U)

Levecke B, Dobson RJ, Speybroeck N, Vercruysse J, Charlier J 2012 Veterinary Parasitology 188(3-4):391-396 (10.1016/j.vetpar.2012.03.020)

Torgerson PR, Paul M, Furrer R 2014 International Journal for Parasitology 44(5):299-303 (10.1016/j.ijpara.2014.01.005)

Kaplan RM, Denwood MJ, Nielsen MK et al. 2023 Veterinary Parasitology 318:109936 (10.1016/j.vetpar.2023.109936)

--- title: "Faecal egg count reduction tests at low counts" description: "At McMaster egg counts a drug that works 99 per cent usually reads as 100 and its interval collapses. Coverage and resistance calls are measured in R." date: "2026-09-08 09:00" categories: [R, parasitology, disease ecology, bootstrap, simulation, ecology tutorial] image: thumbnail.png image-alt: "A line chart on warm off-white paper of interval coverage for a 99 per cent effective drug against the number of eggs expected on the pre-treatment slides, on a log axis from 30 to 6000, with a dashed line at 0.95. A red paired bootstrap line starts at about 0.24 at 30 eggs and climbs steeply to about 0.9 by 400 eggs, then stays just below the dashed line. A gold binomial likelihood line and a green quasi-likelihood line start near 0.98 at low counts; above 1000 eggs the gold line sags to under 0.9 while the green line stays close to 0.93 to 0.95." --- A smallholder with ten ewes wants to know whether the white drench still works. The vet takes a faecal sample from each ewe on the day of treatment and again two weeks later, and the laboratory counts eggs on a McMaster slide. One egg seen under the grid is recorded as fifty eggs per gram of faeces, because the chamber holds a fixed fraction of a diluted sample. The pre-treatment counts come back at a few hundred eggs per gram, the post-treatment counts are all zero, and the report says the drug reduced egg output by 100 per cent. A bootstrap interval, if anyone computes one, runs from 100 to 100. The faecal egg count reduction test is the standard field test for anthelmintic resistance in sheep, goats, cattle and horses. In its paired form it divides the mean post-treatment count by the mean pre-treatment count and compares the reduction with a threshold. The problem this post measures is well known in veterinary parasitology, and nothing here is a new result. Levecke and colleagues (2012) showed by simulation that the outcome of the test depends on the detection limit of the counting method, the number of animals and the level and aggregation of egg output; Torgerson, Paul and Furrer (2014) built the `eggCounts` package around a model that treats the slide count, not the eggs per gram, as the observation; and the 2023 guideline of the World Association for the Advancement of Veterinary Parasitology (Kaplan and colleagues 2023) replaced the 1992 rule (Coles and colleagues 1992) with one that sizes the test by the number of eggs counted and asks for positive evidence before calling a drug either susceptible or resistant. What the post adds is a small, runnable demonstration: how often the report reads 100, how often its interval covers the truth, and how often a resistant drug is called, under both rules. Three posts on this site deal with measurements that stop resolving at the low end, and none of them is this problem. [Values below the detection limit](../values-below-the-detection-limit/) treats a continuous concentration that is censored below a limit and fits the censored likelihood; its test is a trend regression, and nothing in it is a ratio or a treatment decision. [Rounded and coarsened measurements](../rounded-and-coarsened-measurements/) works with a known grid that adds a fixed amount of variance to a continuous value. [Parasite burdens truncated by host death](../truncated-parasite-burden/) is about worm counts and the aggregation parameter. Here the limit is a counting multiplier applied to aggregated counts, and it acts on a ratio at its top end, which is exactly where the treatment decision sits. The failure of the bootstrap that follows is a cousin of the one in [checking a bootstrap](../checking-a-bootstrap/), where resampling cannot produce a value beyond the sample maximum: here it cannot produce a post-treatment egg when none was counted. ```{r setup} #| message: false #| warning: false library(ggplot2) 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), strip.text = element_text(colour = te_ink)) } ``` ## The generating model and three intervals Each animal carries its own true egg output, drawn from a gamma distribution with mean 150 eggs per gram and shape 0.7, so that egg output is strongly aggregated: a few animals carry most of the eggs. Egg output also varies from day to day within an animal, with a coefficient of variation of 0.3, drawn independently for the two sampling days. The treatment multiplies the animal's output by one minus the efficacy. The slide sees a Poisson number of eggs with mean equal to the output divided by the multiplication factor, and the laboratory reports that count times the factor. All of these constants were fixed before anything was run. ```{r design} cv_day <- 0.3 # day-to-day excretion CV within an animal shp_day <- 1 / cv_day^2 k_set <- 0.7 # gamma shape of egg output among animals mu_set <- 150 # mean eggs per gram before treatment n_set <- 10 # animals in the treatment group mf_set <- 50 # McMaster multiplication factor eff_set <- 0.99 # true efficacy n_rep <- 600 # simulated tests per design cell n_boot <- 999 # paired bootstrap resamples per test mcse_max <- sqrt(0.25 / n_rep) sim_test <- function(R, n, mu, k, mf, eff) { epg <- matrix(rgamma(R * n, shape = k, rate = k / mu), R) pre <- matrix(rpois(R * n, epg * rgamma(R * n, shp_day, shp_day) / mf), R) post <- matrix(rpois(R * n, epg * (1 - eff) * rgamma(R * n, shp_day, shp_day) / mf), R) list(pre = pre, post = post, mf = mf) # raw eggs seen on the slide } fecr_of <- function(tst) 100 * (1 - rowSums(tst$post) / rowSums(tst$pre)) ``` The reduction is computed from group arithmetic means, as in the 1992 guideline; the 2023 guideline moves the decision onto the interval instead of a prescribed mean. The multiplication factor cancels from that ratio, so the point estimate can be computed from raw slide counts. The intervals are where the choices differ, and three are compared. The first is a paired bootstrap: resample animals with replacement, keeping each animal's two counts together, and take percentiles of the resampled reduction. Resampling animals is the same as drawing multinomial weights for them, which lets all resamples of all simulated tests be computed with two matrix products. ```{r bootstrap} boot_lim <- function(tst, probs, B = n_boot) { n_an <- ncol(tst$pre) wts <- rmultinom(B, n_an, rep(1 / n_an, n_an)) bpre <- tst$pre %*% wts bpost <- tst$post %*% wts bf <- 100 * (1 - bpost / bpre) bf[!is.finite(bf)] <- NA # resamples with no pre-treatment eggs t(apply(bf, 1, quantile, probs = probs, na.rm = TRUE, names = FALSE)) } ``` The second uses the counts. Under the model in the `eggCounts` package, an animal's pre- and post-treatment slide counts are Poisson with means proportional to the same gamma-distributed output, the second multiplied by one minus the efficacy. Conditional on the animal's total count, its post-treatment count is then binomial with probability r over one plus r, where r is one minus the efficacy, and the animal's own output drops out. The profile likelihood for the efficacy is therefore a binomial likelihood in the share of all counted eggs that were counted after treatment, whatever the aggregation. The interval below inverts that likelihood ratio. When no post-treatment egg is seen its lower limit has a closed form that depends only on the number of eggs counted, which is the quantity the bootstrap never sees. The day-to-day variation breaks the binomial step, because the two days no longer share one mean. The third interval is the quasi-likelihood version of the second: it estimates a dispersion from the Pearson residuals across animals, divides the likelihood ratio by it when it exceeds one, and uses an F critical value. ```{r likelihood} to_eff <- function(p_share) 100 * (1 - p_share / (1 - p_share)) lik_lim <- function(tst, level, quasi = FALSE) { z_post <- rowSums(tst$post) n_all <- rowSums(tst$pre) + z_post p_hat <- z_post / n_all loglik <- function(p) ifelse(z_post > 0, z_post * log(p), 0) + ifelse(n_all - z_post > 0, (n_all - z_post) * log1p(-p), 0) l_max <- loglik(p_hat) phi <- rep(1, length(z_post)) crit <- rep(qchisq(level, 1), length(z_post)) if (quasi) { tot_an <- tst$pre + tst$post pear <- ifelse(tot_an > 0, (tst$post - tot_an * p_hat)^2 / (tot_an * p_hat * (1 - p_hat)), 0) df_an <- pmax(rowSums(tot_an > 0) - 1, 1) phi_raw <- rowSums(pear) / df_an over <- z_post > 0 & z_post < n_all & phi_raw > 1 phi[over] <- phi_raw[over] crit[over] <- qf(level, 1, df_an[over]) } gap <- function(p) 2 * (l_max - loglik(p)) / phi - crit bisect <- function(lo, hi, upper_side) { for (it in 1:60) { mid <- (lo + hi) / 2 inside <- gap(mid) < 0 if (upper_side) { lo <- ifelse(inside, mid, lo); hi <- ifelse(inside, hi, mid) } else { hi <- ifelse(inside, mid, hi); lo <- ifelse(inside, lo, mid) } } (lo + hi) / 2 } p_up <- bisect(p_hat, rep(1 - 1e-12, length(p_hat)), TRUE) p_lo <- ifelse(z_post == 0, 0, bisect(rep(1e-300, length(p_hat)), p_hat, FALSE)) cbind(lo = to_eff(p_up), hi = to_eff(p_lo)) } ``` ## Ten ewes and a factor of fifty ```{r example} set.seed(4417) ex_draws <- sim_test(20, n_set, mu_set, k_set, mf_set, eff_set) all_zero <- which(rowSums(ex_draws$post) == 0 & rowSums(ex_draws$pre) > 0) i_ex <- all_zero[1] ex <- list(pre = ex_draws$pre[i_ex, , drop = FALSE], post = ex_draws$post[i_ex, , drop = FALSE], mf = mf_set) ex_epg_pre <- mf_set * ex$pre ex_mean_pre <- mean(ex_epg_pre) ex_eggs <- sum(ex$pre) ex_zero_pre <- sum(ex$pre == 0) ex_max_epg <- max(ex_epg_pre) ex_fecr <- fecr_of(ex) ex_boot <- boot_lim(ex, c(0.025, 0.975)) ex_lik <- lik_lim(ex, 0.95) ex_lo_closed <- to_eff(1 - exp(-qchisq(0.95, 1) / (2 * ex_eggs))) ex_epg_pre ``` This is the first of twenty simulated tests from the design above in which no egg was seen after treatment, and it was taken as the first such draw rather than chosen. The ten ewes had a mean pre-treatment count of `r sprintf("%.0f", ex_mean_pre)` eggs per gram, the largest `r sprintf("%.0f", ex_max_epg)`, and `r sprintf("%d", ex_zero_pre)` of them showed no egg at all before treatment. Behind that mean are `r sprintf("%d", ex_eggs)` eggs actually counted under the microscope. The reported reduction is `r sprintf("%.0f", ex_fecr)` per cent against a true `r sprintf("%.0f", 100 * eff_set)`. The paired bootstrap interval runs from `r sprintf("%.1f", ex_boot[1])` to `r sprintf("%.1f", ex_boot[2])`. Every resample contains only zeros after treatment, so every resampled reduction is 100, and the interval has no width at all. The likelihood interval runs from `r sprintf("%.1f", ex_lik[1, "lo"])` to `r sprintf("%.1f", ex_lik[1, "hi"])`. Its lower limit has a closed form. With no post-treatment egg among `r sprintf("%d", ex_eggs)` counted, the likelihood ratio reaches the chi-squared critical value at a post-treatment share of one minus the exponential of minus `r sprintf("%.2f", qchisq(0.95, 1) / 2)` over `r sprintf("%d", ex_eggs)`, which is an efficacy of `r sprintf("%.1f", ex_lo_closed)`, and the interval says what the slide supports: a drug somewhere above about nine tenths effective, with no evidence either way about 95 or 99 per cent. ## How often a 99 per cent drug reads as 100 The report of exactly 100 per cent happens when every post-treatment slide is empty, and that probability has a closed form. If every animal excreted the same number of eggs per gram, the post-treatment total would be Poisson with mean n times the output times one minus the efficacy, divided by the factor, and the chance of seeing nothing is the exponential of minus that. With gamma-distributed output among animals each slide is empty with the negative binomial zero probability, and the group probability is its n-th power. The day-to-day variation is not in either formula, so the simulation checks how much it matters. ```{r zeros} mf_grid <- c(50, 25, 15, 1) n_grid <- c(10, 20, 40) eff_grid <- c(0.99, 0.95, 0.90) n_zero_rep <- 20000 p_zero_fixed <- function(n, mu, mf, eff) exp(-n * mu * (1 - eff) / mf) p_zero_gamma <- function(n, mu, k, mf, eff) (1 + mu * (1 - eff) / (k * mf))^(-k * n) pz_fixed_set <- p_zero_fixed(n_set, mu_set, mf_set, eff_set) pz_gamma_set <- p_zero_gamma(n_set, mu_set, k_set, mf_set, eff_set) set.seed(2023) zero_tab <- do.call(rbind, lapply(n_grid, function(n_an) { do.call(rbind, lapply(mf_grid, function(mf) { tst <- sim_test(n_zero_rep, n_an, mu_set, k_set, mf, eff_set) ok <- rowSums(tst$pre) > 0 data.frame(n = n_an, mf = mf, sim = mean(rowSums(tst$post)[ok] == 0), fixed = p_zero_fixed(n_an, mu_set, mf, eff_set), gamma = p_zero_gamma(n_an, mu_set, k_set, mf, eff_set)) })) })) zero_gap_gamma <- max(abs(zero_tab$sim - zero_tab$gamma)) zero_gap_fixed <- max(abs(zero_tab$sim - zero_tab$fixed)) mcse_zero <- sqrt(0.25 / n_zero_rep) zt <- function(n_an, mf, col) zero_tab[zero_tab$n == n_an & zero_tab$mf == mf, col] ``` For ten animals at 150 eggs per gram and a factor of 50, the fixed-output formula gives `r sprintf("%.3f", pz_fixed_set)` and the gamma formula `r sprintf("%.3f", pz_gamma_set)`. Simulating `r sprintf("%d", n_zero_rep)` tests per design gives `r sprintf("%.3f", zt(10, 50, "sim"))`. Across all twelve designs the gamma formula stays within `r sprintf("%.4f", zero_gap_gamma)` of the simulation, against a Monte Carlo standard error of at most `r sprintf("%.4f", mcse_zero)`, so day-to-day variation barely moves this probability; the fixed-output formula is off by up to `r sprintf("%.3f", zero_gap_fixed)`, because aggregation concentrates the eggs in fewer animals and leaves more slides empty. The exponent is the whole story: n times the mean count times one minus the efficacy, over the factor, is the expected number of eggs on all post-treatment slides together. At a factor of 50 that is `r sprintf("%.1f", n_set * mu_set * (1 - eff_set) / mf_set)` eggs for ten animals. Doubling the animals and switching to a factor of 15 raises it to `r sprintf("%.1f", 20 * mu_set * (1 - eff_set) / 15)`, and the share of empty post-treatment totals falls to `r sprintf("%.3f", zt(20, 15, "sim"))`. Only the method that counts every egg (factor 1) removes the problem outright. ```{r fig-zeros} #| fig-cap: "Share of simulated tests of a 99 per cent effective drug that report exactly 100 per cent reduction, by multiplication factor and group size. Lines: gamma closed form; points: simulation." #| fig-alt: "Three rising curves on warm off-white paper, one per group size, plotting the share of tests that report exactly 100 per cent against the multiplication factor on a log axis from 1 to 50. All three start at zero at a factor of 1 and stay near zero until a factor of about 3 for ten animals, 5 for twenty and 10 for forty. At a factor of 50 the red ten-animal curve reaches about 0.75, the gold twenty-animal curve about 0.56 and the green forty-animal curve about 0.31. Simulated points sit on the closed-form curves at factors 1, 15, 25 and 50." zero_tab$group <- factor(sprintf("%d animals", zero_tab$n), levels = sprintf("%d animals", n_grid)) mf_fine <- exp(seq(log(1), log(50), length.out = 120)) zero_curve <- do.call(rbind, lapply(n_grid, function(n_an) data.frame(mf = mf_fine, gamma = p_zero_gamma(n_an, mu_set, k_set, mf_fine, eff_set), group = factor(sprintf("%d animals", n_an), levels = levels(zero_tab$group))))) ggplot(zero_tab, aes(mf, sim, colour = group)) + geom_line(data = zero_curve, aes(y = gamma), linewidth = 0.9) + geom_point(size = 2.4) + scale_x_log10(breaks = c(1, 2, 5, 15, 25, 50)) + scale_colour_manual(values = c(te_rust, te_gold, te_forest), name = NULL) + scale_y_continuous(limits = c(0, 1)) + labs(x = "multiplication factor (eggs per gram per egg seen)", y = "share reporting exactly 100 per cent", title = "A perfect result is mostly a counting result", subtitle = "true efficacy 99 per cent, mean 150 eggs per gram, shape 0.7") + theme_datasheet() + theme(legend.position = "bottom") ``` ## The bootstrap interval collapses where the counts do The share reporting 100 is arithmetic. The interval is not, so it is measured. Each design cell below simulates `r sprintf("%d", n_rep)` tests, a number fixed before running, which puts the Monte Carlo standard error of any coverage at `r sprintf("%.3f", mcse_max)` or less. Tests with no egg at all before treatment are dropped, as a laboratory would drop them. ```{r grid} run_cell <- function(n, mu, k, mf, eff, R = n_rep) { tst <- sim_test(R, n, mu, k, mf, eff) ok <- rowSums(tst$pre) > 0 tst$pre <- tst$pre[ok, , drop = FALSE] tst$post <- tst$post[ok, , drop = FALSE] fecr <- fecr_of(tst) bq <- boot_lim(tst, c(0.025, 0.05, 0.95, 0.975)) l95 <- lik_lim(tst, 0.95); l90 <- lik_lim(tst, 0.90) q95 <- lik_lim(tst, 0.95, quasi = TRUE); q90 <- lik_lim(tst, 0.90, quasi = TRUE) te <- 100 * eff covers <- function(lo, hi) mean(lo <= te & hi >= te) class23 <- function(lo, hi) c(res = mean(hi < 99), sus = mean(lo >= 95 & hi >= 99)) c23b <- class23(bq[, 2], bq[, 3]) c23l <- class23(l90[, 1], l90[, 2]) c23q <- class23(q90[, 1], q90[, 2]) data.frame(n = n, mu = mu, k = k, mf = mf, eff = eff, used = sum(ok), eggs = n * mu / mf, rep100 = mean(fecr == 100), cov_boot = covers(bq[, 1], bq[, 4]), cov_lik = covers(l95[, 1], l95[, 2]), cov_quasi = covers(q95[, 1], q95[, 2]), r92_boot = mean(fecr < 95 & bq[, 1] < 90), r92_lik = mean(fecr < 95 & l95[, 1] < 90), res_boot = c23b[["res"]], sus_boot = c23b[["sus"]], res_lik = c23l[["res"]], sus_lik = c23l[["sus"]], res_quasi = c23q[["res"]], sus_quasi = c23q[["sus"]]) } set.seed(1992) grid_tab <- do.call(rbind, lapply(eff_grid, function(e) do.call(rbind, lapply(n_grid, function(n_an) do.call(rbind, lapply(mf_grid, function(mf) run_cell(n_an, mu_set, k_set, mf, e))))))) gv <- function(n_an, mf, e, col) grid_tab[grid_tab$n == n_an & grid_tab$mf == mf & abs(grid_tab$eff - e) < 1e-9, col] n_dropped <- sum(n_rep - grid_tab$used) top <- grid_tab[abs(grid_tab$eff - 0.99) < 1e-9, ] boot_low_eggs <- range(top$cov_boot[top$eggs <= 60]) lik_mf1 <- range(grid_tab$cov_lik[grid_tab$mf == 1]) quasi_all <- range(grid_tab$cov_quasi) quasi_mf1 <- range(grid_tab$cov_quasi[grid_tab$mf == 1]) quasi_rest_min <- min(grid_tab$cov_quasi[grid_tab$mf != 1]) boot_90 <- range(grid_tab$cov_boot[abs(grid_tab$eff - 0.90) < 1e-9]) ``` At ten animals and a factor of 50 the 99 per cent drug reports exactly 100 in `r sprintf("%.3f", gv(10, 50, 0.99, "rep100"))` of tests, and the paired bootstrap interval covers the true 99 in `r sprintf("%.3f", gv(10, 50, 0.99, "cov_boot"))`, against a nominal 0.95 (`r sprintf("%d", n_dropped)` of the `r sprintf("%d", n_rep * nrow(grid_tab))` simulated tests in the grid were dropped for an empty pre-treatment slide). Where the grid expects 60 or fewer eggs to be counted before treatment, its coverage for this drug is between `r sprintf("%.3f", boot_low_eggs[1])` and `r sprintf("%.3f", boot_low_eggs[2])`. Twenty animals at a factor of 15 lift the coverage to `r sprintf("%.3f", gv(20, 15, 0.99, "cov_boot"))`, and forty animals at a factor of 1 give `r sprintf("%.3f", gv(40, 1, 0.99, "cov_boot"))`. The likelihood interval, on exactly the same simulated counts, covers the true 99 in `r sprintf("%.3f", gv(10, 50, 0.99, "cov_lik"))` of tests at ten animals and a factor of 50. It does not need more counting to be honest at the top end, because a count of zero out of a known number of eggs is information it can use. It fails in the opposite corner. When every egg is counted, day-to-day variation in excretion is larger than the Poisson counting noise the binomial assumes, and the plain likelihood interval covers between `r sprintf("%.3f", lik_mf1[1])` and `r sprintf("%.3f", lik_mf1[2])` across the factor-of-1 designs, worst for the less effective drugs. The quasi-likelihood interval, which estimates that extra variation from the animals, covers between `r sprintf("%.3f", quasi_all[1])` and `r sprintf("%.3f", quasi_all[2])` across all thirty-six cells; every cell in which it falls below `r sprintf("%.3f", quasi_rest_min)` is a factor-of-1 design (those range from `r sprintf("%.3f", quasi_mf1[1])` to `r sprintf("%.3f", quasi_mf1[2])`), where a dispersion estimated from ten to forty aggregated animals is itself noisy. ```{r fig-coverage} #| fig-cap: "Coverage of three 95 per cent intervals against the number of eggs expected on the pre-treatment slides (animals times mean count over the factor), for three true efficacies. Each point is one design cell of 600 simulated tests." #| fig-alt: "Three panels for true efficacies of 99, 95 and 90 per cent, each plotting interval coverage against the expected number of pre-treatment eggs on a log axis from 30 to 6000, with a dashed line at 0.95. In the 99 per cent panel the red paired bootstrap line starts at about 0.24 at 30 eggs and climbs steeply to about 0.9 by 400 eggs, while the gold binomial likelihood and green quasi-likelihood lines stay near 0.95 to 0.99 at low counts. In the 95 and 90 per cent panels the bootstrap starts at about 0.73 and 0.86 and stays below the dashed line; the gold likelihood line holds near 0.92 up to 400 eggs and then drops to about 0.7 and about 0.55 above 1000 eggs; the green quasi-likelihood line stays between about 0.84 and 0.97 throughout." cov_long <- do.call(rbind, lapply( list(c("cov_boot", "paired bootstrap"), c("cov_lik", "binomial likelihood"), c("cov_quasi", "quasi-likelihood")), function(m) data.frame(eggs = grid_tab$eggs, eff = grid_tab$eff, cover = grid_tab[[m[1]]], interval = m[2]))) cov_long$interval <- factor(cov_long$interval, levels = c("paired bootstrap", "binomial likelihood", "quasi-likelihood")) cov_long$panel <- factor(sprintf("true efficacy %.0f per cent", 100 * cov_long$eff), levels = sprintf("true efficacy %.0f per cent", 100 * eff_grid)) cov_mean <- aggregate(cover ~ eggs + panel + interval, data = cov_long, FUN = mean) ggplot(cov_long, aes(eggs, cover, colour = interval)) + geom_hline(yintercept = 0.95, linetype = "dashed", colour = te_body, linewidth = 0.5) + geom_line(data = cov_mean, linewidth = 0.8) + geom_point(size = 2, alpha = 0.9) + facet_wrap(~ panel, ncol = 3) + scale_x_log10(breaks = c(30, 100, 300, 1000, 3000)) + scale_colour_manual(values = c(te_rust, te_gold, te_forest), name = NULL) + scale_y_continuous(limits = c(0, 1)) + labs(x = "eggs expected on the pre-treatment slides (log scale)", y = "coverage of the true efficacy", title = "Each interval fails at a different end", subtitle = "dashed line: nominal 0.95; lines join cell means at equal egg totals") + theme_datasheet() + theme(legend.position = "bottom", panel.spacing.x = unit(1.2, "lines")) ``` For the effective drug the horizontal axis carries most of the pattern. Ten animals at a factor of 25 and twenty animals at a factor of 50 both expect 60 eggs on the pre-treatment slides, and their bootstrap coverages for the 99 per cent drug are `r sprintf("%.3f", gv(10, 25, 0.99, "cov_boot"))` and `r sprintf("%.3f", gv(20, 50, 0.99, "cov_boot"))`. For the 90 per cent drug the bootstrap is at or below nominal across the whole grid, between `r sprintf("%.3f", boot_90[1])` and `r sprintf("%.3f", boot_90[2])`, but never collapses, because post-treatment eggs are seen in most tests and the resamples have something to vary. ## Calling resistance under the 1992 and 2023 rules The 1992 rule calls resistance when the reduction is below 95 per cent and the lower limit of the 95 per cent interval is below 90. The original rule was written for a treated group against an untreated control, with a variance formula for the interval; here it is applied to the paired design with the bootstrap and likelihood intervals. The 2023 guideline recommends the paired design, uses a 90 per cent interval, and compares it with two numbers for the host and drug: an expected efficacy, 99 per cent in the guideline's sheep example, and a lower efficacy threshold, taken here as 95 per cent (the research protocol's grey zone of 95 to 99; the clinical protocol uses 90 to 99). A test is resistant when the upper limit is below the expected efficacy, susceptible when the lower limit reaches the threshold and the upper limit reaches the expected efficacy, and inconclusive otherwise. The guideline suggests a Bayesian approach built for these data and recommends one of two web tools for the analysis (the Bayesian `eggCounts` interface, or fecrt.com with a hybrid frequentist and Bayesian method) unless a qualified statistician does it, and warns that only a few formulae give an interval when the observed reduction is 100 per cent; the bootstrap below is what a laboratory gets if it ignores that advice, and the likelihood intervals are a frequentist stand-in for the model-based analysis, not the analysis itself. The guideline also sizes the treatment group by the total number of eggs counted before the multiplication factor is applied, not by a minimum count per animal. ```{r calls} c92_boot_10 <- gv(10, 50, 0.90, "r92_boot"); c92_boot_20 <- gv(20, 15, 0.90, "r92_boot") c92_lik_10 <- gv(10, 50, 0.90, "r92_lik") fp92_max <- max(grid_tab$r92_boot[abs(grid_tab$eff - 0.99) < 1e-9]) # the headline cell gets its own, larger run so its rate carries a small MC error n_rep_esc <- 5000 set.seed(1992 + 10) esc_tst <- sim_test(n_rep_esc, n_set, mu_set, k_set, mf_set, 0.90) esc_ok <- rowSums(esc_tst$pre) > 0 esc_tst$pre <- esc_tst$pre[esc_ok, ]; esc_tst$post <- esc_tst$post[esc_ok, ] esc_fecr <- fecr_of(esc_tst) esc_call_boot <- mean(esc_fecr < 95 & boot_lim(esc_tst, c(0.025, 0.975))[, 1] < 90) esc_call_lik <- mean(esc_fecr < 95 & lik_lim(esc_tst, 0.95)[, "lo"] < 90) escape_10 <- 1 - esc_call_boot mcse_grid_esc <- sqrt(c92_boot_10 * (1 - c92_boot_10) / gv(10, 50, 0.90, "used")) mcse_escape <- sqrt(escape_10 * (1 - escape_10) / sum(esc_ok)) esc_hi_fecr <- mean(esc_fecr >= 95) designs <- data.frame(n = c(10, 20, 40), mf = c(50, 15, 1), lab = c("10 animals, factor 50", "20 animals, factor 15", "40 animals, factor 1")) set.seed(1203) miss_tst <- sim_test(n_rep, n_set, mu_set, k_set, mf_set, 0.90) miss_ok <- rowSums(miss_tst$pre) > 0 miss_tst$pre <- miss_tst$pre[miss_ok, ]; miss_tst$post <- miss_tst$post[miss_ok, ] miss_up <- boot_lim(miss_tst, c(0.05, 0.95))[, 2] missed <- miss_up >= 99 miss_at100 <- mean(miss_up[missed] == 100) miss_few <- mean(rowSums(miss_tst$post[missed, ] > 0) <= 2) class_long <- do.call(rbind, lapply(seq_len(nrow(designs)), function(i) { do.call(rbind, lapply(eff_grid, function(e) { do.call(rbind, lapply(list(c("boot", "bootstrap"), c("lik", "likelihood"), c("quasi", "quasi")), function(m) { res <- gv(designs$n[i], designs$mf[i], e, paste0("res_", m[1])) sus <- gv(designs$n[i], designs$mf[i], e, paste0("sus_", m[1])) data.frame(design = designs$lab[i], eff = e, interval = m[2], class = c("resistant", "inconclusive", "susceptible"), share = c(res, 1 - res - sus, sus)) })) })) })) ``` With ten animals and a factor of 50, a drug that removes only 90 per cent of egg output is called resistant by the 1992 rule with the bootstrap in `r sprintf("%.3f", c92_boot_10)` of the grid's `r sprintf("%d", n_rep)` tests. Because this is the headline cell, it was rerun on its own with `r sprintf("%d", n_rep_esc)` tests: the drug escapes a resistance call in `r sprintf("%.3f", escape_10)` of them (Monte Carlo standard error `r sprintf("%.3f", mcse_escape)`), a little under one test in four (the grid's run put the escape at `r sprintf("%.3f", 1 - c92_boot_10)`, a low draw given its own standard error of `r sprintf("%.3f", mcse_grid_esc)`). The likelihood interval barely changes this, calling resistance in `r sprintf("%.3f", esc_call_lik)` against `r sprintf("%.3f", esc_call_boot)` for the bootstrap on the same counts, because the escapes come from the point estimate: the rule needs the reduction itself below 95, and with `r sprintf("%.0f", n_set * mu_set / mf_set)` eggs expected on the pre-treatment slides the 90 per cent drug shows a reduction of 95 or more in `r sprintf("%.3f", esc_hi_fecr)` of tests. Twenty animals at a factor of 15 raise the call rate to `r sprintf("%.3f", c92_boot_20)`. The 1992 rule never calls a 99 per cent drug resistant in more than `r sprintf("%.3f", fp92_max)` of tests in this grid. The 2023 rule fed a bootstrap interval fails at low counts in both directions. For ten animals at a factor of 50, it calls the 99 per cent drug susceptible in `r sprintf("%.3f", gv(10, 50, 0.99, "sus_boot"))` of tests, on the strength of an interval that has collapsed onto 100, and calls the 90 per cent drug resistant in only `r sprintf("%.3f", gv(10, 50, 0.90, "res_boot"))`. It even calls a 95 per cent drug, which is below the expected efficacy, susceptible in `r sprintf("%.3f", gv(10, 50, 0.95, "sus_boot"))`. The likelihood interval on the same counts calls the 90 per cent drug resistant in `r sprintf("%.3f", gv(10, 50, 0.90, "res_lik"))`, the 99 per cent drug susceptible in only `r sprintf("%.3f", gv(10, 50, 0.99, "sus_lik"))`, and the 95 per cent drug susceptible in `r sprintf("%.3f", gv(10, 50, 0.95, "sus_lik"))`; its most common answer for the effective drug is inconclusive, in `r sprintf("%.3f", 1 - gv(10, 50, 0.99, "res_lik") - gv(10, 50, 0.99, "sus_lik"))` of tests, which is what `r sprintf("%.0f", n_set * mu_set / mf_set)` expected eggs can support. The missed resistance calls for the 90 per cent drug have a plain cause. In a separate run of `r sprintf("%d", n_rep)` such tests, `r sprintf("%d", sum(missed))` were not called resistant; among those, the share with an upper bootstrap limit of exactly 100 was `r sprintf("%.3f", miss_at100)`, and the share whose post-treatment eggs came from two animals or fewer was `r sprintf("%.3f", miss_few)`. Resamples that leave those animals out read 100, and there are enough of them to hold the upper limit there. ```{r fig-calls} #| fig-cap: "Classification under the 2023 rule (90 per cent interval, expected efficacy 99, threshold 95) for three designs, three true efficacies and three intervals." #| fig-alt: "A grid of stacked bars, rows for true efficacies of 99, 95 and 90 per cent and columns for ten animals at factor 50, twenty at factor 15 and forty at factor 1, with one bar each for bootstrap, likelihood and quasi intervals split into red resistant, grey inconclusive and green susceptible shares. In the ten-animal column the bootstrap bar for the 99 per cent drug is about four fifths green, while the likelihood and quasi bars are under two fifths green and mostly grey; for the 90 per cent drug the bootstrap bar is under two fifths red and mostly grey, against about four fifths and two thirds red for likelihood and quasi. In the other two columns the 99 per cent row is mostly green and the 95 and 90 per cent rows are mostly or entirely red for all three intervals." class_long$class <- factor(class_long$class, levels = c("resistant", "inconclusive", "susceptible")) class_long$interval <- factor(class_long$interval, levels = c("bootstrap", "likelihood", "quasi")) class_long$design <- factor(class_long$design, levels = designs$lab) class_long$effl <- factor(sprintf("true %.0f%%", 100 * class_long$eff), levels = sprintf("true %.0f%%", 100 * eff_grid)) ggplot(class_long, aes(interval, share, fill = class)) + geom_col(width = 0.75, colour = te_paper, linewidth = 0.3) + facet_grid(effl ~ design) + scale_fill_manual(values = c(te_rust, te_line, te_forest), name = NULL) + scale_y_continuous(breaks = c(0, 0.5, 1)) + labs(x = NULL, y = "share of simulated tests", title = "A bootstrap interval misleads the 2023 rule", subtitle = "rows: true efficacy; columns: design") + theme_datasheet() + theme(legend.position = "bottom", axis.text.x = element_text(size = 9)) ``` With forty animals counted at a factor of 1 the three intervals agree on the classification almost everywhere. The one visible difference there is that the plain likelihood interval calls the 99 per cent drug resistant in `r sprintf("%.3f", gv(40, 1, 0.99, "res_lik"))` of tests against `r sprintf("%.3f", gv(40, 1, 0.99, "res_boot"))` for the bootstrap, which is its undercoverage at high counts turning into false alarms. ## Aggregation matters less than the eggs counted The negative binomial shape and the mean count were both held fixed above. The grid below varies them at twenty animals and a factor of 50, for the 99 and 90 per cent drugs. ```{r aggregation} k_grid <- c(0.3, 0.7, 2) mu_grid <- c(50, 150, 600) set.seed(3009) agg_tab <- do.call(rbind, lapply(c(0.99, 0.90), function(e) do.call(rbind, lapply(k_grid, function(kk) do.call(rbind, lapply(mu_grid, function(mu) run_cell(20, mu, kk, mf_set, e))))))) av <- function(kk, mu, e, col) agg_tab[agg_tab$k == kk & agg_tab$mu == mu & abs(agg_tab$eff - e) < 1e-9, col] agg_top <- agg_tab[abs(agg_tab$eff - 0.99) < 1e-9, ] spread_k_top <- max(tapply(agg_top$cov_boot, agg_top$mu, function(v) diff(range(v)))) spread_mu_top <- max(tapply(agg_top$cov_boot, agg_top$k, function(v) diff(range(v)))) agg_low <- agg_tab[abs(agg_tab$eff - 0.90) < 1e-9, ] spread_k_low <- max(tapply(agg_low$cov_boot, agg_low$mu, function(v) diff(range(v)))) spread_mu_low <- max(tapply(agg_low$cov_boot, agg_low$k, function(v) diff(range(v)))) ``` For the 99 per cent drug, changing the shape from 0.3 to 2 at a fixed mean moves the bootstrap coverage by at most `r sprintf("%.3f", spread_k_top)`, while changing the mean from 50 to 600 at a fixed shape moves it by up to `r sprintf("%.3f", spread_mu_top)`. At a mean of 50 eggs per gram, twenty animals expect twenty eggs on their slides, and the bootstrap covers the true 99 in `r sprintf("%.3f", av(0.7, 50, 0.99, "cov_boot"))` of tests at shape 0.7. The likelihood interval covers in `r sprintf("%.3f", av(0.7, 50, 0.99, "cov_lik"))` there. For the 90 per cent drug the two effects are of similar size, because eggs are seen after treatment and the bootstrap is resampling a skewed ratio rather than a column of zeros: the shape moves its coverage by up to `r sprintf("%.3f", spread_k_low)` and the mean by up to `r sprintf("%.3f", spread_mu_low)`. At a mean of 50, its coverage is `r sprintf("%.3f", av(0.3, 50, 0.90, "cov_boot"))` at shape 0.3 and `r sprintf("%.3f", av(2, 50, 0.90, "cov_boot"))` at shape 2; the quasi-likelihood interval covers `r sprintf("%.3f", av(0.3, 50, 0.90, "cov_quasi"))` and `r sprintf("%.3f", av(2, 50, 0.90, "cov_quasi"))`. ```{r fig-aggregation} #| fig-cap: "Coverage of the paired bootstrap and quasi-likelihood intervals at twenty animals and a factor of 50, against the mean pre-treatment count, for three levels of aggregation." #| fig-alt: "Two panels for true efficacies of 99 and 90 per cent plotting coverage against mean eggs per gram at 50, 150 and 600, with red, gold and green lines for shapes 0.3, 0.7 and 2, solid for the paired bootstrap and long-dashed for the quasi-likelihood interval, and a dashed reference at 0.95. In the 99 per cent panel the three solid bootstrap lines rise together from about 0.16 to 0.2 at a mean of 50, through about 0.4 to 0.45 at 150, to between about 0.77 and 0.89 at 600, while the dashed lines sit near 0.99 falling to about 0.92. In the 90 per cent panel the solid lines start between about 0.75 and 0.86 at 50 and flatten near 0.9, and the dashed lines stay between about 0.92 and 0.97." agg_long <- rbind( data.frame(agg_tab[, c("k", "mu", "eff")], cover = agg_tab$cov_boot, interval = "paired bootstrap"), data.frame(agg_tab[, c("k", "mu", "eff")], cover = agg_tab$cov_quasi, interval = "quasi-likelihood")) agg_long$shape <- factor(sprintf("shape %.1f", agg_long$k), levels = sprintf("shape %.1f", k_grid)) agg_long$panel <- factor(sprintf("true efficacy %.0f per cent", 100 * agg_long$eff), levels = sprintf("true efficacy %.0f per cent", c(99, 90))) ggplot(agg_long, aes(mu, cover, colour = shape, linetype = interval)) + geom_hline(yintercept = 0.95, linetype = "dashed", colour = te_body, linewidth = 0.5) + geom_line(linewidth = 0.8) + geom_point(size = 2) + facet_wrap(~ panel) + scale_x_log10(breaks = mu_grid) + scale_colour_manual(values = c(te_rust, te_gold, te_forest), name = NULL) + scale_linetype_manual(values = c("solid", "longdash"), name = NULL) + scale_y_continuous(limits = c(0, 1)) + guides(colour = guide_legend(order = 1), linetype = guide_legend(order = 2, nrow = 2)) + labs(x = "mean eggs per gram before treatment (log scale)", y = "coverage of the true efficacy", title = "The mean count sets the failure", subtitle = "twenty animals, factor 50; dashed line: nominal 0.95") + theme_datasheet() + theme(legend.position = "bottom", legend.key.width = unit(2.2, "lines"), panel.spacing.x = unit(2, "lines")) ``` ## What to report Report the raw number of eggs counted before and after treatment, as well as the multiplication factor. Those two totals are the information in the test; the eggs per gram figures are a rescaling of them, and a reduction of 100 per cent from 30 eggs means something quite different from 100 per cent from 3000. Do not report a bootstrap interval of zero width. When every post-treatment count is zero the paired bootstrap has nothing to resample and returns the point estimate twice; quoting it as an interval invites the reader to treat a counting result as a precise one. A likelihood interval on the counts, whose lower limit depends on the eggs counted, is a defensible replacement at low counts, and the model-based Bayesian intervals in `eggCounts` are built on the same observation; analysing the counts with a method built for these data is also the guideline's own advice. Say which rule was applied and with which interval. Under the 2023 rule the interval is the whole decision, and the same counts produced a susceptible call for the effective drug in `r sprintf("%.3f", gv(10, 50, 0.99, "sus_boot"))` of tests with one interval and `r sprintf("%.3f", gv(10, 50, 0.99, "sus_lik"))` with another. An inconclusive result from a small, coarse test is a correct result. Plan the test on eggs counted rather than on animals. For the effective drug, the two designs here that expect the same 60 pre-treatment eggs from different group sizes and factors gave bootstrap coverages of `r sprintf("%.3f", gv(10, 25, 0.99, "cov_boot"))` and `r sprintf("%.3f", gv(20, 50, 0.99, "cov_boot"))`, and the mean count moved coverage far more than aggregation did; that is the same logic as the 2023 guideline sizing groups by counted eggs. ## Honest limits The simulation uses a single gamma-Poisson generator with independent day-to-day variation of fixed size. Real post-treatment counts also carry zero inflation (animals with no patent infection), a change in the proportion of worm species that survive treatment, and the gap between egg output and worm burden that a reduction in egg output does not measure. The binomial likelihood interval is exact only under the model without day-to-day variation, and its failure at high counts was shown here with one value of that variation; a smaller one would shrink the failure and a larger one would widen it. The Bayesian models in `eggCounts`, and other model-based intervals a laboratory might use, were not fitted. Their behaviour at low counts depends on the priors, and nothing above says how close the likelihood intervals here come to them. The comparison is between the paired bootstrap, which is easy to compute, and two frequentist intervals built on the counts. The 1992 rule was applied to a paired design with a resampling interval rather than to the treated-versus-control design and variance formula it was written for, and the 2023 thresholds were taken from the guideline's sheep example with the research protocol; values for other hosts and drugs may differ. The clinical protocol's lower threshold of 90 per cent would call more drugs susceptible; the resistant call, which depends only on the upper limit and the expected efficacy of 99, would not change. The low resistance subclass of the 2023 guideline was not separated out. The scene is a small flock of sheep, and the numbers transfer to cattle, goats and horses only through the eggs counted. Wildlife egg count surveys share the low counts and the coarse factors but rarely have a paired treatment design, and the reduction test does not apply to them as it stands. ## References Coles GC, Bauer C, Borgsteede FHM, Geerts S, Klei TR, Taylor MA, Waller PJ 1992 Veterinary Parasitology 44(1-2):35-44 (10.1016/0304-4017(92)90141-U) Levecke B, Dobson RJ, Speybroeck N, Vercruysse J, Charlier J 2012 Veterinary Parasitology 188(3-4):391-396 (10.1016/j.vetpar.2012.03.020) Torgerson PR, Paul M, Furrer R 2014 International Journal for Parasitology 44(5):299-303 (10.1016/j.ijpara.2014.01.005) Kaplan RM, Denwood MJ, Nielsen MK et al. 2023 Veterinary Parasitology 318:109936 (10.1016/j.vetpar.2023.109936) ## Related tutorials - [Values below the detection limit in R](../values-below-the-detection-limit/) - [Parasite burdens truncated by host death](../truncated-parasite-burden/) - [Checking a bootstrap: when resampling fails](../checking-a-bootstrap/) - [Prevalence from an imperfect test](../prevalence-from-an-imperfect-test/)