Values below the detection limit in R

Substituting half the detection limit for a non-detect is not neutral, and a laboratory that improves over time turns a flat series into a decline. Worked in R.

Corrected 14 October 2026: the threshold of about a quarter of the record below the limit was offered as one worth remembering, and it is now tied to the spread of the concentrations as well, which the sweep holds fixed; the paragraph on uncertain limits no longer calls the fabricated trend a lower bound on what a changing limit can do; and the passage on summaries no longer says that half the limit landing close here is why everyone uses it.

Corrected 1 October 2026: the description of the plotted record said the detected points show no trend, and it now says their lower edge steps down with the limit; the threshold of about a quarter of the record below the limit is now tied to a limit that falls, since a fixed limit makes no trend; and the passage on summaries no longer says that no substituted value gets the mean and the spread both right, since half the limit lands close on both.

Updated 26 September 2026: a new section, When the file says <0.1, measures what as.numeric() and stripping the “<” do to a laboratory export on the improving-laboratory record, and the list of fields in the opening no longer includes qPCR, which a new passage in Honest limits explains.

The laboratory returns “< 0.05”. That is not a missing value, because the sample was measured and the measurement carries information: the concentration is somewhere between zero and 0.05. It is also not 0.05, and not 0.025, and not zero, although the spreadsheet will end up containing one of those three.

Substituting a number for a non-detect is the standard practice in nutrient chemistry, ecotoxicology and trace elements, and it has two costs that are easy to miss. The first is bias in the summary statistics, and it runs in opposite directions for the mean and the spread (Helsel 2006 calls substitution fabricating data). The second is worse and is the subject of most of this post: when the detection limit changes over the life of a monitoring programme, substitution manufactures a trend in a population that never moved.

What substitution does to a summary

A lognormal concentration, one detection limit, four of the usual rules.

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"),
          axis.text        = element_text(colour = te_body))
}

mu_log <- -0.6
sd_log <- 1.0

set.seed(20260806)
truth <- rlnorm(20000, mu_log, sd_log)
dl    <- 0.5

rules <- c(zero = 0, half = dl / 2, root2 = dl / sqrt(2), full = dl)
subst <- sapply(rules, function(v) ifelse(truth < dl, v, truth))

round(rbind(mean = c(true = mean(truth), apply(subst, 2, mean)),
            sd   = c(true = sd(truth),   apply(subst, 2, sd))), 3)

      true  zero  half root2  full
mean 0.895 0.771 0.886 0.934 1.001
sd   1.132 1.203 1.134 1.108 1.075

46 per cent of these samples are non-detects, which is ordinary for a trace analyte. The true mean is 0.895 and the true standard deviation 1.132. Substituting zero pulls the mean down to 0.771 and pushes the spread up to 1.203; substituting the full limit pushes the mean up to 1.001 and squeezes the spread down to 1.075.

The mean and the standard deviation move in opposite directions as the substituted value rises, so a value that gets one of them right gets the other right only by luck. Half the limit happens to land close on both here, and it is close by coincidence of this distribution and this censoring fraction rather than by any property of the rule.

The same programme, an improving laboratory

Now the case that matters. Twenty years of monitoring, thirty samples a year, and a concentration distribution that does not change at all. The only thing that changes is the laboratory: the detection limit drops from 1.0 to 0.5 to 0.2 to 0.1 across four five-year eras, because the instruments got better.

n_year  <- 30
years   <- 1:20
dl_era  <- c(rep(1.0, 5), rep(0.5, 5), rep(0.2, 5), rep(0.1, 5))

one_series <- function() {
  d <- data.frame(year = rep(years, each = n_year),
                  dl   = rep(dl_era, each = n_year))
  d$conc     <- rlnorm(nrow(d), mu_log, sd_log)     # never changes
  d$censored <- d$conc < d$dl
  d$reported <- ifelse(d$censored, d$dl / 2, d$conc)
  d
}

set.seed(14)
series <- one_series()
c(censored_overall = mean(series$censored),
  censored_era1    = mean(series$censored[series$year <= 5]),
  censored_era4    = mean(series$censored[series$year > 15]))

censored_overall    censored_era1    censored_era4 
      0.37000000       0.72000000       0.05333333 

Over the whole record 37 per cent of samples are non-detects, but the fraction falls from 72 per cent in the first era to 5 per cent in the last. That decline is real and it is a fact about the laboratory. The substituted value that stands in for those non-detects falls with it, from 0.5 to 0.05, and that is where the trend comes from.

naive <- lm(log(reported) ~ year, data = series)
round(coef(summary(naive)), 5)

            Estimate Std. Error  t value Pr(>|t|)
(Intercept) -0.35842    0.08188 -4.37747  0.00001
year        -0.02091    0.00683 -3.05949  0.00232

The slope is -0.0209 per year on the log scale, with p = 0.002. Over the twenty years that is a 32.8 per cent decline in a quantity that was drawn from the same distribution every single year.

series$kind <- ifelse(series$censored, "non-detect, substituted", "detected")

ggplot(series, aes(x = year, y = reported)) +
  geom_point(aes(colour = kind, shape = kind), size = 1.5, alpha = 0.6) +
  geom_step(aes(y = dl), direction = "hv", linewidth = 0.9, colour = te_ink) +
  geom_smooth(method = "lm", formula = y ~ x, se = FALSE,
              linewidth = 1.1, colour = te_rust) +
  annotate("text", x = 1, y = 1.35, hjust = 0, size = 3.4, colour = te_ink,
           label = "detection limit") +
  scale_y_log10() +
  scale_colour_manual(values = c("detected" = te_forest,
                                 "non-detect, substituted" = te_gold)) +
  scale_shape_manual(values = c(16, 17)) +
  labs(x = "year", y = "reported concentration, log scale",
       colour = NULL, shape = NULL,
       title = "A flat population with a downward trend") +
  theme_datasheet() +
  theme(legend.position = "top")

Figure 1: One simulated twenty-year record. Detected values are shown as points, non-detects at the substituted half-limit, with the stepped detection limit and the fitted trend through the reported values.

The likelihood that does not need a substituted value

A non-detect is an interval observation: the concentration is below the limit, and the probability of seeing that is the distribution function evaluated at the limit. Detected values contribute their density. Adding the two gives a likelihood that uses every sample without inventing a number for any of them, and it handles a limit that changes between samples without any extra machinery, because each observation carries its own limit.

The mechanics of interval-censored likelihoods are set out in the post on rounded and coarsened measurements; this is the one-sided case, where the lower edge of the interval is zero.

nll <- function(par, d) {
  m <- par[1] + par[2] * d$year
  s <- exp(par[3])
  ll <- ifelse(d$censored,
               plnorm(d$dl, m, s, log.p = TRUE),
               dlnorm(d$conc_obs, m, s, log = TRUE))
  -sum(ll)
}

series$conc_obs <- ifelse(series$censored, NA, series$conc)
fit_ml <- optim(c(mu_log, 0, log(sd_log)), nll, d = series,
                method = "BFGS", hessian = TRUE)
se_ml  <- sqrt(diag(solve(fit_ml$hessian)))[2]

round(c(slope = fit_ml$par[2], se = se_ml,
        lower = fit_ml$par[2] - 1.96 * se_ml,
        upper = fit_ml$par[2] + 1.96 * se_ml,
        sd = exp(fit_ml$par[3]), converged = fit_ml$convergence), 5)

    slope        se     lower     upper        sd converged 
  0.00053   0.00882  -0.01674   0.01781   1.08503   0.00000 

The maximum likelihood slope is 0.00053 per year with a standard error of 0.00882, so the interval runs from -0.0167 to 0.0178 and sits squarely on the truth of zero. The recovered standard deviation on the log scale is 1.085 against a truth of 1.0. The same data that gave a 33 per cent decline gives no detectable trend once the non-detects are treated as what they are.

The two estimates use exactly the same observations. The difference is that one of them was told which values are bounds.

How much censoring it takes

With the limit falling in the same pattern, scaled up or down, the size of the fabricated trend is a function of how much of the record is below it, so it is worth seeing the whole curve rather than one point on it.

sweep_dl <- function(scale_up, reps = 60) {
  out <- data.frame()
  for (k in scale_up) {
    res <- replicate(reps, {
      d <- one_series()
      d$dl <- d$dl * k
      d$censored <- d$conc < d$dl
      d$reported <- ifelse(d$censored, d$dl / 2, d$conc)
      d$conc_obs <- ifelse(d$censored, NA, d$conc)
      b_naive <- coef(lm(log(reported) ~ year, data = d))[["year"]]
      p_naive <- coef(summary(lm(log(reported) ~ year, data = d)))["year", "Pr(>|t|)"]
      f <- optim(c(mu_log, 0, log(sd_log)), nll, d = d, method = "BFGS")
      c(cens = mean(d$censored), naive = b_naive, ml = f$par[2],
        sig = p_naive < 0.05)
    })
    out <- rbind(out, data.frame(
      censored = mean(res["cens", ]),
      naive    = mean(res["naive", ]),
      ml       = mean(res["ml", ]),
      sig      = mean(res["sig", ]),
      decline  = 100 * (1 - exp(mean(res["naive", ]) * 19))))
  }
  out
}

set.seed(3)
sw <- sweep_dl(c(0.4, 0.7, 1.0, 1.5, 2.2))
round(sw, 4)

  censored   naive      ml    sig decline
1   0.1431  0.0000  0.0000 0.0167  0.0002
2   0.2565 -0.0080 -0.0008 0.2167 14.0431
3   0.3450 -0.0184 -0.0019 0.8833 29.5440
4   0.4612 -0.0326  0.0009 0.9833 46.1770
5   0.5699 -0.0528 -0.0004 1.0000 63.3380

At 14 per cent censoring there is nothing to see: the fabricated decline is 0.0 per cent and the naive test calls it significant in 2 per cent of runs, which is the nominal five per cent to within the noise of sixty replicates. By 26 per cent it is a 14.0 per cent decline, significant 22 per cent of the time; at 34 per cent, which is roughly the record plotted above, 29.5 per cent and 88 per cent; and at 57 per cent, 63.3 per cent in every single run. The maximum likelihood slope stays within 0.0019 of zero across the whole range.

The threshold in this sweep, for a limit that falls the way this one does, is the one between the first two rows. Somewhere around a quarter of the record being non-detect, substitution stops being a rounding convention and starts being the thing that produces the result. The sweep moves only the limit and holds the spread of the concentrations fixed, so that share belongs to this spread as well as to this limit. What makes the trend is a limit that moves during the record, and the share below it only sets the size: with one limit held throughout, every year has the same share replaced however large it is, and the last row of the table in the next section shows substitution then making no trend.

long <- rbind(
  data.frame(censored = sw$censored, decline = sw$decline,
             method = "half-limit substitution"),
  data.frame(censored = sw$censored,
             decline = 100 * (1 - exp(sw$ml * 19)),
             method = "censored maximum likelihood"))

ggplot(long, aes(x = 100 * censored, y = decline, colour = method,
                 shape = method)) +
  geom_hline(yintercept = 0, colour = te_ink, linewidth = 0.5) +
  geom_line(linewidth = 0.9) +
  geom_point(size = 2.9) +
  scale_colour_manual(values = c("half-limit substitution" = te_rust,
                                 "censored maximum likelihood" = te_forest)) +
  scale_shape_manual(values = c(16, 17)) +
  labs(x = "per cent of the record below the detection limit",
       y = "apparent decline over twenty years, per cent",
       colour = NULL, shape = NULL,
       title = "The trend is bought entirely with substituted values") +
  theme_datasheet() +
  theme(legend.position = "top")

Figure 2: Fabricated twenty-year decline from half-limit substitution against the fraction of the record below the limit, with the censored maximum likelihood estimate for comparison. Sixty simulated programmes per point.

When the file says <0.1

Every chunk so far started from a column of numbers. A laboratory export is a column of text, with each non-detect written as "<0.1" beside the measured values, and one of two things usually happens to it. as.numeric() cannot read "<0.1", so it returns NA for every non-detect and warns once, “NAs introduced by coercion”, however many values it removed; lm() then drops those rows as it drops any missing value. The usual cure for the warning, stripping the < with sub(), puts the full limit in place of each non-detect with nothing left to say it was a bound. One is deletion, the other substitution of the full limit, and neither was chosen.

lab_file <- function(conc, dl)          # what the laboratory export holds
  ifelse(conc < dl, paste0("<", dl), as.character(signif(conc, 3)))
lab <- lab_file(series$conc, series$dl)
msgs <- character(0)
as_number <- withCallingHandlers(as.numeric(lab), warning = function(w) {
  msgs <<- c(msgs, conditionMessage(w)); invokeRestart("muffleWarning") })
stopifnot(length(msgs) == 1, msgs == "NAs introduced by coercion",
          nobs(lm(log(as_number) ~ series$year)) == sum(!series$censored))
lab       <- trimws(lab)                       # the parse that keeps the flag
lab_flag  <- grepl("^<", lab)
lab_value <- as.numeric(sub("^<", "", lab))
stopifnot(!anyNA(lab_value))                   # every value parsed
parsed <- data.frame(year = series$year, dl = series$dl, censored = lab_flag,
                     conc_obs = ifelse(lab_flag, NA, lab_value))
fit_parsed <- optim(c(mu_log, 0, log(sd_log)), nll, d = parsed, method = "BFGS")
head(lab, 8)

[1] "<1"   "3.06" "4.58" "2.45" "<1"   "1.88" "<1"   "1.6" 

In the plotted record as.numeric() turned 222 of the 600 values into NA with a single warning. The lines from the comment down drop stray spaces, keep the flag and the limit, and check that every value parsed; a code such as ND or <LOD needs its own rule. The censored likelihood on the parsed file gives a slope of 0.00052 per year against 0.00053 from the exact values, the difference being the export’s three significant figures. Reading field data into R has the same type rule turning a count column into text over one n/a, and Missing values in summaries has na.rm = TRUE fabricating a decline from months nobody walked. Here the rows that go are the bottom of every year, and a smaller slice of it in the later eras.

decline   <- function(b) 100 * (1 - exp(19 * b))
fit_slope <- function(y, yr) coef(summary(lm(log(y) ~ yr)))[2, c(1, 4)]
lab_rules <- c("deleted", "stripped", "half")
rule_run <- function(L, reps = 200) replicate(reps, {   # L: limit in each year
  d <- transform(one_series(), dl = rep(L, each = n_year))
  s <- lab_file(d$conc, d$dl)
  d$censored <- grepl("^<", s)
  value <- as.numeric(sub("^<", "", s))
  d$conc_obs <- ifelse(d$censored, NA, value)
  f <- optim(c(mu_log, 0, log(sd_log)), nll, d = d, method = "BFGS")
  c(cens = mean(d$censored), ml = f$par[2],
    deleted  = fit_slope(suppressWarnings(as.numeric(s)), d$year),
    stripped = fit_slope(value, d$year),
    half     = fit_slope(ifelse(d$censored, d$dl / 2, value), d$year),
    mean_true = mean(d$conc), mean_deleted = mean(value[!d$censored]))
})
expected_decline <- function(L, rule) {        # closed form, no simulation
  a <- (log(L) - mu_log) / sd_log
  p <- pnorm(a)                                   # share below the limit
  above <- mu_log + sd_log * dnorm(a) / (1 - p)   # inverse Mills ratio
  era_mean <- switch(rule, deleted = above, stripped = p * log(L) + (1 - p) * above,
                     half = p * log(L / 2) + (1 - p) * above)
  kept <- if (rule == "deleted") 1 - p else rep(1, length(L))
  decline(coef(lm(era_mean ~ years, weights = kept))[[2]])
}
set.seed(344)
arms <- list(x0.7 = dl_era * 0.7, x1.0 = dl_era, x1.5 = dl_era * 1.5,
             fixed0.5 = rep(0.5, 20))
rank_tab <- t(sapply(arms, function(L) {
  r <- rule_run(L)
  b <- sapply(lab_rules, function(v) r[paste0(v, ".Estimate"), ])
  c(censored = 100 * mean(r["cens", ]), apply(b, 2, function(x) decline(mean(x))),
    ml = decline(mean(r["ml", ])),
    cf = sapply(lab_rules, function(v) expected_decline(L, v)),
    sig_deleted = mean(r["deleted.Pr(>|t|)", ] < 0.05),
    mc_se = max(1900 * exp(19 * colMeans(b)) * apply(b, 2, sd) / sqrt(nrow(b))),
    mean_true = mean(r["mean_true", ]), mean_deleted = mean(r["mean_deleted", ]))
}))
round(rank_tab[, 1:9], 2)

         censored deleted stripped  half    ml cf.deleted cf.stripped cf.half
x0.7        25.73   64.23    46.49 12.65 -0.16      64.26       46.47   12.83
x1.0        34.65   70.78    59.05 26.86 -0.56      71.29       59.40   27.34
x1.5        46.04   76.99    71.68 46.26 -0.84      77.05       71.69   46.30
fixed0.5    46.18    0.52     0.12 -0.01  0.01       0.00        0.00    0.00
         sig_deleted
x0.7            1.00
x1.0            1.00
x1.5            1.00
fixed0.5        0.04

The first three rows are the improving laboratory with its limits scaled by 0.7, 1 and 1.5, 200 simulated programmes per row. Deletion turns the flat population into a decline of 64.2 to 77.0 per cent over the twenty years, significant in every run of every row, and stripping the < gives 46.5 to 71.7 per cent. Both are larger than the half-limit substitution the rest of the post warns about, 12.6 to 46.3 per cent (the sweep above rerun with 200 programmes per row instead of 60), and the censored likelihood stays within 0.9 points of zero. The Monte Carlo standard error of each deleted, stripped and half-limit decline is at most 1.0 points. None of the deleted, stripped and half-limit columns is a finding of the simulation. The values that survive in an era are a lognormal cut at that era’s limit L, and their mean log is mu + sigma * dnorm(a) / (1 - pnorm(a)) with a = (log(L) - mu) / sigma (mu and sigma are mu_log and sd_log): sigma times the inverse Mills ratio above mu. A line through those era means, weighted by the share kept, gives cf.deleted; mixing the limit or half of it back in with weight pnorm(a) gives cf.stripped and cf.half. The simulation reproduces the closed forms to within 0.6 points.

The ranking belongs to this design, a limit that falls while the population stands still. With one limit of 0.5 throughout (the last row, 46 per cent non-detects) every year loses the same slice, and deletion makes no trend: 0.5 per cent, significant in 4 per cent of runs. It moves the level instead: the mean of what is left is 1.44 against a true 0.90. Parse the file into a flag and a value whatever the design; what deletion costs depends on how the limit moved.

What to record and what to report

Keep the flag and the limit, not the substituted number. Two columns, one logical and one numeric, and the analysis can do anything it likes afterwards. A file that contains only 0.25 has destroyed the distinction between a measurement and a bound, and no later analysis can recover it.

Report the censoring fraction alongside any summary from censored data, because a mean from a record that is 14 per cent non-detect and one that is 57 per cent non-detect are not comparable quantities even when they carry the same units. And if the limit changed during the programme, say when and by how much, because that is a covariate whether or not anyone modelled it.

Honest limits

The censored likelihood needs a distributional assumption where substitution needs none, and that is a real trade (Shoari and Dube 2018 argue for analysing concentration data with the nondetects kept in). Lognormal is conventional for concentrations and is what generated the data here, so this comparison is friendly to the method. A misspecified distribution biases the estimate in its own way, and the check is the usual one: fit an alternative family and see whether the conclusion moves. Non-parametric alternatives exist, Kaplan-Meier run on the reversed scale being the common one (Helsel 2012 sets out the methods), and they carry less structure at the cost of less precision.

Very heavy censoring defeats everything. Past roughly three quarters non-detect there is little left to identify the lower tail with, and the honest report is the detection frequency itself rather than a mean nobody can estimate.

qPCR looks like the same problem and is not, because at low copy number a qPCR non-detect is mostly a reaction that received no template molecule: a Poisson zero, with probability exp(-lambda) for a mean of lambda copies per reaction, not a concentration hiding below a limit. The chunk below simulates studies of 200 samples with a median of 0.3 copies per reaction and three replicate reactions each, reads the copies off a standard curve with some Cq noise, and fits this post’s likelihood (nll() with the slope held at zero) with each non-detect censored at half a copy, below almost every reading a single template molecule gives.

qpcr_study <- function(n_sample = 200, n_rep = 3, median_copies = 0.3) {
  lambda <- rlnorm(n_sample, log(median_copies), sd_log)  # copies per reaction
  copies <- rpois(n_sample * n_rep, rep(lambda, n_rep))    # template per well
  cq  <- 38 - 3.32 * log10(copies) + rnorm(length(copies), 0, 0.3)
  est <- ifelse(copies > 0, 10^((38 - cq) / 3.32), NA)    # standard curve
  d <- data.frame(year = 0, dl = 0.5, censored = is.na(est) | est < 0.5)
  d$conc_obs <- ifelse(d$censored, NA, est)
  f <- optim(c(0, 0), function(p) nll(c(p[1], 0, p[2]), d), method = "BFGS")
  c(nd = mean(d$censored), below_limit = mean(est < 0.5, na.rm = TRUE),
    empty = mean(is.na(est[d$censored])), nd_at_3 = mean(is.na(est) | est < 3),
    censored_fit = exp(f$par[1] + exp(f$par[2])^2 / 2),
    nd_as_zero = mean(ifelse(d$censored, 0, est)))
}
set.seed(2026)
qpcr <- replicate(200, qpcr_study())
true_copies <- 0.3 * exp(sd_log^2 / 2)
qpcr_ratio <- qpcr[c("censored_fit", "nd_as_zero"), ] / true_copies
round(c(rowMeans(qpcr[c("nd", "below_limit", "empty", "nd_at_3"), ]),
        rowMeans(qpcr_ratio), quantile(qpcr_ratio["censored_fit", ], c(0.1, 0.9)),
        se_nd_as_zero = sd(qpcr_ratio["nd_as_zero", ]) / sqrt(ncol(qpcr)),
        cq_noise = exp((0.3 * log(10) / 3.32)^2 / 2)), 4)

           nd   below_limit         empty       nd_at_3  censored_fit 
       0.6810        0.0004        0.9998        0.9722        1.4096 
   nd_as_zero           10%           90% se_nd_as_zero      cq_noise 
       1.0233        1.2267        1.5668        0.0092        1.0219 

With 68 per cent non-detects, short of the roughly three quarters past which very heavy censoring defeats everything (a limit of detection of three copies would make it 97 per cent), the censored fit puts the mean at 1.41 times the true mean copies per reaction, averaged over 200 simulated studies. It is above the truth in every one of them, between 1.23 and 1.57 times from the 10th to the 90th percentile over the studies, and only 0.04 per cent of the wells that amplified read below the half-copy limit. In this simulation a non-detect is an empty reaction (99.98 per cent of them are), so entering it as zero adds no bias by construction: the plain mean gives 1.023, reproducing the 1.022 that the Cq noise alone predicts once back-transformed. That repair is for the mean only: the distribution of sample concentrations needs a Poisson-lognormal model, and a non-detect is a Poisson zero only when inhibition and amplification failure are ruled out and the standard curve is taken as known. For concentrations measured against a limit, the advice above stands.

This is censoring, not truncation. A non-detect is present in the record with partial information; a truncated observation is absent altogether, which is a different problem with a different fix, treated in the post on truncated parasite burdens.

The detection limit is itself an estimate, it varies between batches and matrices, and laboratories do not always report the value they actually used. Where the limit is uncertain the interval is wider than modelled here, and the fabricated trend above is one example of what a changing limit can do.

Finally, the mechanism is not confined to detection limits. Any bound that moves during a monitoring programme does the same thing: a reporting threshold, a minimum countable size, a changed survey protocol. The post on splicing a monitoring series treats the version where a method change puts a step in the mean; here an unchanged population acquires a slope.

References

Helsel DR 2006 Chemosphere 65(11):2434-2439 (10.1016/j.chemosphere.2006.04.051)

Shoari N, Dube JS 2018 Environmental Toxicology and Chemistry 37(3):643-656 (10.1002/etc.4046)

Helsel DR 2012 Statistics for Censored Environmental Data Using Minitab and R, 2nd edition, Wiley, ISBN 978-0-470-47988-9

Related tutorials

• Rounded and coarsened measurements
• Parasite burdens truncated by host death
• Faecal egg count reduction tests at low counts
• Compositional analysis and the zero problem

Reproducibility

• Code last run: 2026-11-06
• R version 4.5.3 (2026-03-11), aarch64-apple-darwin20
• Packages: ggplot2 4.0.3