── Row 3: Normal Q–Q on log scale ──────────────────────────────────────────
for (k in 1:3) { lv <- list(log_true, log_zero, log_ros)[[k]] qq <- qqnorm(lv, plot.it = FALSE) plot(qq\(x, qq\)y, pch = 19, cex = 0.85, col = adjustcolor(col3[k], 0.8), xlab = “Theoretical quantiles”, ylab = “Sample quantiles”, main = ““) qqline(lv, col =”#888”, lwd = 1.4, lty = 2) if (k == 2) { # Annotate the hook n_hook <- sum(d_zero == 0) y_hook <- sort(log_zero)[n_hook] x_hook <- qq$x[order(lv)[n_hook]] arrows(x_hook - 0.5, y_hook - 0.3, x_hook - 0.1, y_hook - 0.05, length = 0.08, col = “#c62828”, lwd = 1.5) text(x_hook - 0.55, y_hook - 0.38, “downward”, col = “#c62828”, cex = 0.70, adj = 1) } if (k == 1) mtext(“Q–Q plot (log scale)”, side = 2, line = 2.6, cex = 0.78, col = “#555”, las = 0) }
mtext(“How data-handling strategy distorts distributional diagnostics”, outer = TRUE, cex = 0.92, font = 2, col = “#1e3a2f”)
par(mfrow = c(1, 1))
Read each column top-to-bottom. The **true values** (left column) show a right-skewed
histogram, a symmetric bell on the log scale, and a straight Q–Q line — exactly what a
log-normal distribution looks like. The **naïve zero-substitution** column (centre)
tells a completely different story: a spike of zeros at the left of the histogram, a
distorted bimodal kernel density with an artefactual left lobe, and a severe downward
hook in the Q–Q plot. That hook is entirely a product of the data-handling choice — it
does not reflect any real feature of the environment. A statistician who saw only this
column might wrongly conclude the data are bimodal, or that the log-normal model is
inappropriate, or that some kind of mixture distribution is needed. The **ROS-imputed**
column (right) closely recovers the true picture.
<div class="callout callout-warn">
<strong>⚠️ The Q–Q hook is a diagnostic red flag</strong>
If you see a sharp downward hook at the left tail of a Q–Q plot on log-transformed
environmental data, check your non-detects first. The hook almost always means zeros
(or a clump of identical substituted values like MDL/2) have been included in the
log-transformed dataset. It is not a distribution problem — it is a censoring problem
masquerading as one.
</div>
### 1.4 How the bias compounds: a numerical experiment
The figure below makes the bias concrete across a range of censoring fractions, using
1,000 simulated datasets at each level.
```{r fig-bias-vs-censoring, echo=TRUE, fig.height=4.0,
fig.cap="Bias in the estimated mean as a function of censoring fraction, for three data-handling strategies. Each point is the average of 1,000 simulated log-normal datasets (n = 50, same parameters throughout); the horizontal dashed line is the true mean. Substituting zero produces severe negative bias that grows linearly with censoring fraction. Substituting RL/2 is better but still biased above 20–30%. ROS (using cenros()) tracks the true mean closely across the full range shown."}
set.seed(999)
n_sim <- 500 # simulations per censoring level (kept low for knit speed)
n_obs <- 50
mu_b <- 2.0; sd_b <- 0.85
true_mn <- exp(mu_b + sd_b^2 / 2) # log-normal mean
cen_targets <- seq(0.05, 0.70, by = 0.05)
results <- lapply(cen_targets, function(target_cen) {
# Find RL that produces approximately target censoring fraction
rl_b <- qlnorm(target_cen, mu_b, sd_b)
means_zero <- means_half <- means_ros <- numeric(n_sim)
for (i in seq_len(n_sim)) {
tv <- rlnorm(n_obs, mu_b, sd_b)
cen <- tv < rl_b
obs <- ifelse(cen, 0, tv)
means_zero[i] <- mean(obs)
means_half[i] <- mean(ifelse(cen, rl_b / 2, obs))
r <- tryCatch(cenros(obs, cen), error = function(e) NULL)
means_ros[i] <- if (!is.null(r)) mean(r$modeled) else NA
}
data.frame(
cen_frac = target_cen,
bias_zero = mean(means_zero, na.rm = TRUE) - true_mn,
bias_half = mean(means_half, na.rm = TRUE) - true_mn,
bias_ros = mean(means_ros, na.rm = TRUE) - true_mn
)
}) |> do.call(rbind, args = _)
par(mar = c(4.5, 5, 3.5, 1.5), las = 1)
plot(results$cen_frac * 100, results$bias_zero,
type = "b", pch = 19, lwd = 2,
col = "#c62828",
xlab = "Censoring fraction (%)",
ylab = expression(widehat(mu) - mu[true]~~(bias~~in~~mean)),
main = "Bias in estimated mean vs. censoring fraction",
ylim = range(c(results$bias_zero, results$bias_half,
results$bias_ros), na.rm = TRUE) * c(1.1, 1.3),
cex.main = 0.93)
lines(results$cen_frac * 100, results$bias_half,
type = "b", pch = 17, lwd = 2, col = "#e65100")
lines(results$cen_frac * 100, results$bias_ros,
type = "b", pch = 15, lwd = 2, col = "#2d5a3d")
abline(h = 0, lty = 2, col = "#555", lwd = 1.5)
legend("bottomleft",
legend = c("Substitute 0", "Substitute RL/2", "ROS (principled)"),
col = c("#c62828", "#e65100", "#2d5a3d"),
pch = c(19, 17, 15), lwd = 2, bty = "n", cex = 0.88)Two take-aways from this simulation. First, substituting zero is never unbiased: even at 5% censoring, the mean is already measurably underestimated. Second, RL/2 substitution is approximately unbiased up to about 15–20% censoring but degrades rapidly beyond that. ROS stays close to the dashed zero-bias line across the full range shown. This is why the EPA, USGS, and Helsel (2012) all recommend moving beyond substitution methods for datasets with more than 15% non-detects.
1.5 The regulatory and scientific context
Environmental agencies routinely encounter datasets with 30–80% non-detects:
- EPA Method 200.8 for metals in water has MDLs of 0.02–2 µg/L for trace metals; clean background groundwater frequently falls below these limits.
- E. coli monitoring in recreational waters: when conditions are good, many samples return “below 1 CFU/100 mL.”
- PM2.5 daily averages near rural background sites can fall below 2–3 µg/m³, the typical FRM instrument noise floor.
- Total phosphorus in oligotrophic reference lakes: 5–20 µg/L MDLs mean half the spring samples may be non-detects.
Getting the statistics right matters for regulatory compliance (is this site exceeding the 303(d) listing criterion?), trend detection (is lead declining post-ban?), and risk assessment (what is the 95th percentile exposure?).
2. Core Concepts
Detection limit vocabulary
| Term | Abbreviation | Definition |
|---|---|---|
| Method Detection Limit | MDL | Lowest concentration distinguishable from blank at 99% confidence (EPA) |
| Reporting Limit | RL | Practical limit used in lab reports; often 3–5 × MDL |
| Limit of Detection | LOD | Generic term; use carefully — meaning varies by lab |
| Limit of Quantitation | LOQ | Level above which measurement is acceptably precise; typically 3–10 × MDL |
| Minimum Reporting Level | MRL | USGS term equivalent to RL |
In this tutorial we use reporting limit (RL) as the
generic term for the censoring threshold. A measurement reported as
< RL is a non-detect.
Left censoring versus right censoring
In survival analysis (right censoring), we observe time up to some cutoff c and know only that T > c. Here, with left censoring, we observe that a concentration C is known only to satisfy C < RL. The Kaplan–Meier estimator handles both through a symmetric reversal: we flip the data, apply the standard KM algorithm, then flip the estimated CDF back.
The fraction censored matters
The appropriate method depends heavily on what fraction of your data is below the reporting limit:
| Censoring fraction | Recommended approaches |
|---|---|
| < 15% | Any substitution method is likely acceptable; ROS preferred |
| 15–50% | ROS or MLE; avoid substitution |
| 50–80% | MLE or Kaplan–Meier; be cautious about distributional assumptions |
| > 80% | Kaplan–Meier only; summary statistics become unreliable |
⚠️ Multiple detection limits
A single dataset often contains measurements from different laboratories, different years, or different analytical methods — each with its own RL. This is called multiply-censored data. The Kaplan–Meier and ROS methods handle multiple RLs naturally. Simple substitution methods (½ MDL) cannot, because there is no single substitution value to use.
3. Packages and Data Setup
Required packages
# Core statistical methods for non-detects
library(NADA) # cenros(), cenfit(), cencdf(), cenmle(), and hypothesis tests
library(NADA2) # updated companion: cenperm2(), cen2means(), cenQQ()
library(survival) # underlying Surv() infrastructure
# Data wrangling and visualization
library(dplyr)
library(tidyr)
library(ggplot2)
library(patchwork) # combine ggplot panels
library(knitr)
library(kableExtra)💡 Installing NADA and NADA2
NADA is on CRAN: install.packages("NADA").
NADA2 may require the development version:
install.packages("NADA2") or
remotes::install_github("SwampThingPaul/NADA2"). Both
implement Dennis Helsel’s methods from Statistics for Censored
Environmental Data (Helsel 2012), the authoritative reference for
this topic.
Dataset 1: E. coli in recreational water (water quality)
We simulate monthly E. coli monitoring at three beach sites over two years. Units are CFU/100 mL. The state recreational water standard is 235 CFU/100 mL (single-sample). The laboratory reporting limit is 1 CFU/100 mL; samples below this are recorded as 0 in the raw data (a common field convention — zero means non-detect, not absence).
set.seed(2024)
# Helper: simulate log-normal samples with left-censoring at RL
sim_censored <- function(n, meanlog, sdlog, rl) {
true_vals <- rlnorm(n, meanlog = meanlog, sdlog = sdlog)
detected <- true_vals >= rl
# Observed value: true value if detected, 0 if non-detect (field convention)
obs <- ifelse(detected, true_vals, 0)
data.frame(
true = true_vals,
observed = obs,
censored = !detected, # TRUE = non-detect
rl = rl
)
}
# Three beach sites with different E. coli levels
set.seed(101)
ecoli <- bind_rows(
sim_censored(24, meanlog = 3.2, sdlog = 1.1, rl = 1) |> mutate(site = "Beach A"),
sim_censored(24, meanlog = 2.4, sdlog = 1.3, rl = 1) |> mutate(site = "Beach B"),
sim_censored(24, meanlog = 4.1, sdlog = 0.9, rl = 1) |> mutate(site = "Beach C")
)
# Summarise censoring fractions
ecoli |>
group_by(site) |>
summarise(
n = n(),
n_nd = sum(censored),
pct_nd = round(100 * mean(censored), 1),
min_det = round(min(observed[!censored]), 1),
max_det = round(max(observed[!censored]), 1)
) |>
kable(col.names = c("Site", "n", "Non-detects", "% ND", "Min detected", "Max detected"),
caption = "E. coli dataset summary (CFU/100 mL; RL = 1 CFU/100 mL)")| Site | n | Non-detects | % ND | Min detected | Max detected |
|---|---|---|---|---|---|
| Beach A | 24 | 0 | 0.0 | 2.6 | 118.0 |
| Beach B | 24 | 2 | 8.3 | 1.3 | 51.8 |
| Beach C | 24 | 0 | 0.0 | 15.9 | 259.3 |
Dataset 2: Total phosphorus in lake monitoring (water quality)
Total phosphorus (TP) in a temperate lake monitored quarterly at surface and deep stations. Units are µg/L. The RL is 5 µg/L; values below this are coded as 0. Elevated TP indicates eutrophication risk.
set.seed(202)
tp <- bind_rows(
sim_censored(30, meanlog = 2.1, sdlog = 0.8, rl = 5) |> mutate(depth = "Surface"),
sim_censored(30, meanlog = 2.8, sdlog = 0.7, rl = 5) |> mutate(depth = "Deep")
)
tp |>
group_by(depth) |>
summarise(n = n(), n_nd = sum(censored),
pct_nd = round(100 * mean(censored), 1)) |>
kable(col.names = c("Depth", "n", "Non-detects", "% ND"),
caption = "Total phosphorus dataset summary (µg/L; RL = 5 µg/L)")| Depth | n | Non-detects | % ND |
|---|---|---|---|
| Deep | 30 | 1 | 3.3 |
| Surface | 30 | 7 | 23.3 |
Dataset 3: Lead (Pb) in drinking water (water quality)
Lead sampling at 50 household taps under a Tier 2 monitoring plan. Units are µg/L. The EPA action level is 15 µg/L; the RL is 0.5 µg/L. Many modern homes with no lead service lines return non-detects.
set.seed(303)
pb <- sim_censored(50, meanlog = 0.8, sdlog = 1.4, rl = 0.5) |>
mutate(tap_id = paste0("Tap_", sprintf("%02d", 1:50)))
cat("Lead dataset: n =", nrow(pb),
" | Non-detects:", sum(pb$censored),
" (", round(100 * mean(pb$censored), 1), "%)\n")## Lead dataset: n = 50 | Non-detects: 6 ( 12 %)cat("Detection range:", round(range(pb$observed[!pb$censored]), 2), "µg/L\n")## Detection range: 0.56 56.5 µg/LDataset 4: PM2.5 daily concentrations (air quality)
Daily PM2.5 at a rural background monitoring site over one year. Units are µg/m³. The RL is 2 µg/m³; days with very clean air register as 0. The 24-hour NAAQS is 35 µg/m³.
set.seed(404)
# PM2.5: slight seasonal signal (higher in winter), mostly below RL in summer
pm25 <- data.frame(
doy = 1:365
) |>
mutate(
# true seasonal log-mean: peaks in winter (day ~15, ~350), dips in summer
meanlog = 1.6 + 0.5 * cos(2 * pi * (doy - 15) / 365),
true_conc = rlnorm(365, meanlog = meanlog, sdlog = 0.9),
observed = ifelse(true_conc >= 2, true_conc, 0),
censored = true_conc < 2,
rl = 2,
season = cut(doy, breaks = c(0, 79, 171, 263, 355, 366),
labels = c("Winter","Spring","Summer","Fall","Winter2"),
include.lowest = TRUE)
) |>
mutate(season = recode(season, "Winter2" = "Winter"))
cat("PM2.5 dataset: n =", nrow(pm25),
" | Non-detects:", sum(pm25$censored),
" (", round(100 * mean(pm25$censored), 1), "%)\n")## PM2.5 dataset: n = 365 | Non-detects: 59 ( 16.2 %)4. Why Naïve Approaches Fail
Before introducing principled methods, let us see exactly how badly the common shortcuts perform using the lead dataset, where we have the true values available for comparison.
# Compute summary stats under each substitution strategy
sub_zero <- pb$observed # non-detects treated as 0
sub_half <- ifelse(pb$censored, pb$rl / 2, pb$observed) # RL/2
sub_rl <- ifelse(pb$censored, pb$rl, pb$observed) # substitute RL
true_vals <- pb$true
# Use ROS (Section 6) for comparison — computed here, explained later
ros_fit <- cenros(pb$observed, pb$censored)
sub_ros <- ros_fit$modeled
methods_df <- bind_rows(
data.frame(value = true_vals, method = "True values"),
data.frame(value = sub_zero, method = "Substitute 0"),
data.frame(value = sub_half, method = "Substitute RL/2"),
data.frame(value = sub_rl, method = "Substitute RL"),
data.frame(value = sub_ros, method = "ROS (principled)")
) |>
mutate(method = factor(method, levels = c("True values","Substitute 0",
"Substitute RL/2","Substitute RL",
"ROS (principled)")))
ggplot(methods_df, aes(x = log(value + 0.01), fill = method, colour = method)) +
geom_density(alpha = 0.30, linewidth = 0.7) +
scale_fill_manual(values = c("#1e3a2f","#c62828","#e65100","#c8832a","#2d5a3d")) +
scale_colour_manual(values = c("#1e3a2f","#c62828","#e65100","#c8832a","#2d5a3d")) +
labs(x = "log(Pb concentration, µg/L)", y = "Density",
title = "Lead (Pb): substitution methods vs. true distribution",
fill = NULL, colour = NULL) +
theme_minimal(base_size = 13) +
theme(legend.position = "top",
plot.title = element_text(family = "serif", face = "bold", colour = "#1e3a2f"))
Distribution of lead concentrations by substitution method vs. the true distribution. Substituting zero collapses the lower tail entirely; substituting RL/2 is closer but systematically underestimates the lower tail. The ROS method (introduced in Section 6) closely recovers the true distribution.
bias_table <- data.frame(
Method = c("True values", "Substitute 0", "Substitute RL/2",
"Substitute RL", "ROS (principled)"),
Mean = round(c(mean(true_vals), mean(sub_zero), mean(sub_half),
mean(sub_rl), mean(sub_ros)), 3),
Median = round(c(median(true_vals), median(sub_zero), median(sub_half),
median(sub_rl), median(sub_ros)), 3),
SD = round(c(sd(true_vals), sd(sub_zero), sd(sub_half),
sd(sub_rl), sd(sub_ros)), 3),
P95 = round(c(quantile(true_vals, 0.95), quantile(sub_zero, 0.95),
quantile(sub_half, 0.95), quantile(sub_rl, 0.95),
quantile(sub_ros, 0.95)), 3)
)
kable(bias_table,
col.names = c("Method", "Mean (µg/L)", "Median (µg/L)", "SD", "95th pct"),
caption = "Summary statistics under different non-detect treatments. True values serve as the benchmark. ROS best recovers the true distribution statistics.") |>
kable_styling(bootstrap_options = "striped") |>
row_spec(1, bold = TRUE, color = "#1e3a2f") |>
row_spec(5, bold = TRUE, color = "#2d5a3d")| Method | Mean (µg/L) | Median (µg/L) | SD | 95th pct |
|---|---|---|---|---|
| True values | 5.206 | 2.413 | 8.828 | 17.02 |
| Substitute 0 | 5.163 | 2.413 | 8.853 | 17.02 |
| Substitute RL/2 | 5.193 | 2.413 | 8.836 | 17.02 |
| Substitute RL | 5.223 | 2.413 | 8.819 | 17.02 |
| ROS (principled) | 5.194 | 2.413 | 8.835 | 17.02 |
⚠️ The substitution bias is systematic, not random
Substituting zero always underestimates the mean and standard deviation. This is not sampling error that averages out — it is structural bias that worsens as the censoring fraction increases. With 40% non-detects, substituting zero can underestimate the true mean by 20–50%.
5. Kaplan–Meier for Left-Censored Data
The Kaplan–Meier (KM) estimator is the non-parametric workhorse for censored data. It makes no distributional assumption — it simply estimates the empirical CDF from the data as they are, treating non-detects as left-censored observations.
The flip trick
Standard KM is designed for right-censoring (T > c). To handle
left-censoring (C < RL), we flip the data: let x’ =
max(x) − x and apply the right-censored KM estimator, then flip the
estimated CDF back. The NADA package’s
cenfit() does this automatically.
# Kaplan-Meier for E. coli at Beach A
beach_a <- ecoli |> filter(site == "Beach A")
# cenfit(obs, censored_indicator) — censored = TRUE means non-detect
km_a <- cenfit(beach_a$observed, beach_a$censored)
# Print summary statistics estimated from KM
print(km_a)## n n.cen median mean sd
## 24.00000 0.00000 19.83262 31.09971 28.21992# Plot the KM CDF
par(mar = c(5, 4.5, 3.5, 2), las = 1,
col.axis = "#3a3a38", col.lab = "#1a1a18", cex.axis = 0.9)
plot(km_a,
xlab = "E. coli (CFU/100 mL)",
ylab = "Cumulative probability",
main = "Kaplan–Meier CDF — E. coli Beach A",
col = "#4a7c5e", lwd = 2)
# Mark the standard
abline(v = 235, lty = 2, col = "#c8832a", lwd = 1.8)
text(250, 0.15, "235 CFU/100 mL\n(state standard)", col = "#c8832a",
adj = 0, cex = 0.82)
# Mark the median
med_a <- median(km_a)
abline(h = 0.5, lty = 3, col = "#c62828", lwd = 1.4)
abline(v = med_a, lty = 3, col = "#c62828", lwd = 1.4)
text(med_a + 5, 0.05, paste0("Median ≈ ", round(med_a, 1)), col = "#c62828",
adj = 0, cex = 0.80)
Kaplan–Meier empirical CDF for E. coli concentrations at Beach A. The shaded step function drops at each detected value; the flat segments between detection limits are where non-detects contribute. The red dashed line marks the 50th percentile (median). The dashed vertical line marks the state standard of 235 CFU/100 mL.
Comparing three beach sites with KM
One of KM’s great strengths is comparing groups without any
distributional assumption. We fit separate KM curves and then apply the
non-parametric log-rank-type test
(cendiff()) for equality of the distributions.
km_all <- cenfit(ecoli$observed, ecoli$censored, groups = factor(ecoli$site))
par(mar = c(5, 4.5, 3.5, 2), las = 1)
plot(km_all,
xlab = "E. coli (CFU/100 mL)",
ylab = "Cumulative probability",
main = "Kaplan–Meier CDF — E. coli by Site",
col = c("#1e3a2f", "#4a7c5e", "#c8832a"), lwd = 2)
abline(v = 235, lty = 2, col = "#c62828", lwd = 1.5)
abline(h = 0.50, lty = 3, col = "#999", lwd = 1)
legend("bottomright",
legend = c("Beach A", "Beach B", "Beach C", "235 CFU standard"),
col = c("#1e3a2f", "#4a7c5e", "#c8832a", "#c62828"),
lty = c(1, 1, 1, 2), lwd = 2, bty = "n", cex = 0.85)
Kaplan–Meier CDFs for E. coli at three beach sites. The curves are well separated: Beach C (high E. coli) stacks far right, Beach B is intermediate, and Beach A has the lowest concentrations. The horizontal dashed line at p = 0.50 reads off each median; the vertical line marks the 235 CFU/100 mL standard.
# Non-parametric test for equality of distributions across sites
cendiff(ecoli$observed, ecoli$censored, groups = factor(ecoli$site))## N Observed Expected (O-E)^2/E (O-E)^2/V
## factor(ecoli$site)=Beach A 24 11.0 13.29 0.381 0.895
## factor(ecoli$site)=Beach B 24 7.9 16.35 4.362 12.157
## factor(ecoli$site)=Beach C 24 17.5 6.82 16.771 28.499
##
## Chisq= 30.2 on 2 degrees of freedom, p= 3e-07📌 Reading the output of cendiff()
cendiff() performs a Peto–Peto modified log-rank
test adapted for left-censored data. The test statistic and
p-value test H₀: all groups have the same distribution. A small p-value
(< 0.05) tells you the distributions differ significantly — it does
not tell you which pairs differ. For pairwise
comparisons, use cenperm2() from NADA2 (see
Section 8).
6. Regression on Order Statistics (ROS)
Regression on Order Statistics (ROS) is a semi-parametric method that estimates summary statistics by:
- Ranking all observations together (detected and non-detected).
- Computing expected normal (or log-normal) scores for the detected values.
- Fitting a regression of detected values on their normal scores.
- Extrapolating that regression line to predict what the non-detected values would have been, then substituting those modeled values.
The result is a complete dataset of modeled values that you can summarise or plot like ordinary data, as long as you remember that the non-detect slots have been imputed, not measured.
Fitting ROS for phosphorus
# ROS fit for total phosphorus — surface station
surf_tp <- tp |> filter(depth == "Surface")
ros_tp <- cenros(surf_tp$observed, surf_tp$censored)
# The returned object has $modeled — the full imputed dataset
cat("ROS-estimated mean TP (surface):", round(mean(ros_tp$modeled), 2), "µg/L\n")## ROS-estimated mean TP (surface): 13.37 µg/Lcat("ROS-estimated SD TP (surface):", round(sd(ros_tp$modeled), 2), "µg/L\n")## ROS-estimated SD TP (surface): 14.41 µg/Lcat("95th percentile (surface): ", round(quantile(ros_tp$modeled, 0.95), 2), "µg/L\n")## 95th percentile (surface): 35.77 µg/Lpar(mar = c(5, 4.5, 3.5, 2), las = 1)
# cenros() $modeled contains the full imputed dataset (NDs replaced by fitted values).
# We compute z-scores from Blom plotting positions on the complete modeled vector
# rather than using $zscores, which can contain NAs when many values are censored.
n_all <- length(ros_tp$modeled)
pp <- (rank(ros_tp$modeled) - 0.375) / (n_all + 0.25)
z_all <- qnorm(pp)
is_det <- !ros_tp$censored
x_detected <- ros_tp$modeled[is_det]
y_detected <- z_all[is_det]
x_modeled <- ros_tp$modeled[!is_det]
y_modeled <- z_all[!is_det]
# Compute safe axis limits (finite values only)
x_rng <- range(c(x_detected, x_modeled), finite = TRUE)
y_rng <- range(c(y_detected, y_modeled), finite = TRUE)
x_rng <- x_rng + c(-0.05, 0.05) * diff(x_rng)
y_rng <- y_rng + c(-0.05, 0.05) * diff(y_rng)
plot(x_detected, y_detected,
pch = 19, col = "#4a7c5e", cex = 1.1,
xlab = "Total phosphorus (µg/L)",
ylab = "Normal score (z)",
main = "ROS Probability Plot — Total Phosphorus (Surface)",
xlim = x_rng, ylim = y_rng)
if (length(x_modeled) > 0)
points(x_modeled, y_modeled, pch = 1, col = "#c62828", cex = 1.2, lwd = 1.5)
abline(v = surf_tp$rl[1], lty = 2, col = "#999", lwd = 1.5)
text(surf_tp$rl[1] + diff(x_rng) * 0.02, y_rng[1] + diff(y_rng) * 0.05,
"RL = 5 µg/L", col = "#777", adj = 0, cex = 0.80)
legend("topleft",
legend = c("Detected values", "ROS-modeled (imputed)"),
pch = c(19, 1),
col = c("#4a7c5e", "#c62828"),
bty = "n", cex = 0.88)
Probability plot for total phosphorus (surface samples) from the ROS fit. Blue circles are detected values; open red circles are the modeled (imputed) values extrapolated below the reporting limit (grey dashed vertical line). The modeled values extend the linear trend established by the detected values into the censored region.
Comparing surface and deep TP with ROS
deep_tp <- tp |> filter(depth == "Deep")
ros_deep <- cenros(deep_tp$observed, deep_tp$censored)
ros_summary <- data.frame(
Depth = c("Surface", "Deep"),
n = c(nrow(surf_tp), nrow(deep_tp)),
pct_nd = round(c(mean(surf_tp$censored), mean(deep_tp$censored)) * 100, 1),
mean_ros = round(c(mean(ros_tp$modeled), mean(ros_deep$modeled)), 2),
sd_ros = round(c(sd(ros_tp$modeled), sd(ros_deep$modeled)), 2),
p95_ros = round(c(quantile(ros_tp$modeled, 0.95),
quantile(ros_deep$modeled, 0.95)), 2)
)
kable(ros_summary,
col.names = c("Depth", "n", "% ND", "Mean (µg/L)", "SD", "95th pct"),
caption = "ROS-estimated total phosphorus statistics by depth stratum.")| Depth | n | % ND | Mean (µg/L) | SD | 95th pct |
|---|---|---|---|---|---|
| Surface | 30 | 23.3 | 13.37 | 14.41 | 35.77 |
| Deep | 30 | 3.3 | 24.32 | 18.08 | 59.84 |
7. Maximum Likelihood Estimation (MLE)
Maximum Likelihood Estimation fits a parametric distribution (usually log-normal) to the observed data, explicitly accounting for censoring in the likelihood function. Each detected value contributes the density f(x); each non-detect at RL contributes the probability F(RL) — the probability that the true value fell below the detection limit.
The log-likelihood for a log-normal distribution with n detected values and m non-detects is:
\[\ell(\mu, \sigma) = \sum_{i=1}^{n} \log f(x_i;\, \mu, \sigma) \;+\; \sum_{j=1}^{m} \log F(RL_j;\, \mu, \sigma)\]
Maximising this over (µ, σ) gives the MLE of the log-normal parameters, from which means, percentiles, and confidence intervals follow analytically.
MLE for lead (Pb) in drinking water
# cenmle() / survreg() require strictly positive values even for non-detects.
# Non-detects are coded as 0 in pb$observed; replace with the RL so that
# survreg sees a valid left-censored value (the true value is known to be < RL).
pb_obs_mle <- ifelse(pb$censored, pb$rl, pb$observed)
mle_pb <- cenmle(pb_obs_mle, pb$censored, dist = "lognormal")
print(mle_pb)## n n.cen median mean sd
## 50.000000 6.000000 2.290607 5.489231 11.954398# Extract mean and confidence interval
mle_stats <- summary(mle_pb)
cat("MLE log-normal fit for Pb:\n")## MLE log-normal fit for Pb:cat(" Mean: ", round(mean(mle_pb), 3), "µg/L\n")## Mean: 5.489 1.691 3.332 9.044 µg/Lcat(" Median: ", round(median(mle_pb),3), "µg/L\n")## Median: 2.291 µg/Lcat(" 95th pct: ", round(quantile(mle_pb, p = 0.95), 3), "µg/L\n")## 95th pct: 20.155 µg/LChecking the distributional assumption with a Q–Q plot
MLE depends on the assumed distribution being correct. Always check with a Q–Q plot.
par(mar = c(5, 4.5, 3.5, 2), las = 1)
cenQQ(pb_obs_mle, pb$censored)
title(main = "Log-normal Q–Q Plot — Lead (Pb)", cex.main = 0.95)
Censored Q–Q plot for lead (Pb) concentrations. Detected values (filled circles) are plotted against log-normal quantiles; a good fit follows the diagonal reference line. Non-detects are shown at the reporting limit boundary (0.5 µg/L).
When MLE is preferred over ROS
| Criterion | Kaplan–Meier | ROS | MLE |
|---|---|---|---|
| Distributional assumption | None | Log-normal (soft) | Log-normal (strict) |
| Handles multiple RLs | Yes | Yes | Yes |
| Censoring fraction | Any | < 80% | Any |
| Provides CIs | Bootstrap only | Bootstrap | Analytical |
| Easy to visualise | Yes (CDF) | Yes (prob plot) | Moderate |
| Sensitive to outliers | Low | Moderate | Moderate |
8. Putting It Together: PM2.5 Analysis
The PM2.5 dataset combines seasonal structure, a meaningful detection limit, and a regulatory threshold. We walk through the complete analysis workflow.
Step 1 — Visualise the raw data
pm25_plot <- pm25 |>
mutate(season = factor(season, levels = c("Winter","Spring","Summer","Fall")))
ggplot(pm25_plot, aes(x = doy, y = observed)) +
geom_point(data = filter(pm25_plot, !censored),
aes(colour = season), size = 1.5, alpha = 0.7) +
geom_point(data = filter(pm25_plot, censored),
shape = 25, colour = "#c62828", fill = "#c62828",
size = 1.8, alpha = 0.6) +
geom_hline(yintercept = 35, linetype = "dashed", colour = "#c8832a", linewidth = 0.8) +
annotate("text", x = 10, y = 37, label = "35 µg/m³ NAAQS", colour = "#c8832a",
hjust = 0, size = 3.5) +
scale_colour_manual(values = c("#1e3a2f","#4a7c5e","#c8832a","#8aad96"),
name = "Season") +
labs(x = "Day of year", y = expression(PM[2.5]~(µg/m^3)),
title = expression(Daily~PM[2.5]~concentrations~(rural~background~site))) +
theme_minimal(base_size = 13) +
theme(legend.position = "top",
plot.title = element_text(family = "serif", face = "bold", colour = "#1e3a2f"))
Daily PM2.5 concentrations over one year. Grey points are detected values; red triangles pointing down indicate non-detect days (coded as 0 in the raw data). The 35 µg/m³ NAAQS is shown as a dashed horizontal line. Non-detects cluster in summer, when clean-air conditions are most frequent.
Step 2 — Fit KM and compute seasonal summary statistics
# Compute KM-estimated mean and median by season
seasons <- c("Winter","Spring","Summer","Fall")
pm25_stats <- lapply(seasons, function(s) {
d <- filter(pm25, season == s)
# Replace 0 (non-detect) with RL for cenmle; KM/cenfit handle 0s natively
obs_mle <- ifelse(d$censored, d$rl, d$observed)
km <- cenfit(d$observed, d$censored)
mle <- tryCatch(cenmle(obs_mle, d$censored, dist = "lognormal"),
error = function(e) NULL)
data.frame(
Season = s,
n = nrow(d),
pct_nd = round(100 * mean(d$censored), 1),
km_mean = round(mean(km), 2),
km_med = round(median(km), 2),
mle_mean = if (!is.null(mle)) round(mean(mle), 2) else NA
)
}) |> bind_rows()
kable(pm25_stats,
col.names = c("Season", "n", "% ND", "KM mean", "KM median", "MLE mean"),
caption = "PM2.5 seasonal summary statistics (µg/m³) from Kaplan–Meier and MLE fits.")| Season | n | % ND | KM mean | KM median | MLE mean | |
|---|---|---|---|---|---|---|
| mean…1 | Winter | 89 | 5.6 | 11.25 | 8.51 | 11.64 |
| se…2 | Winter | 89 | 5.6 | 0.96 | 8.51 | 1.28 |
| 0.95LCL…3 | Winter | 89 | 5.6 | 9.36 | 8.51 | 9.47 |
| 0.95UCL…4 | Winter | 89 | 5.6 | 13.14 | 8.51 | 14.32 |
| mean…5 | Spring | 92 | 17.4 | 6.31 | 4.41 | 6.24 |
| se…6 | Spring | 92 | 17.4 | 0.59 | 4.41 | 0.65 |
| 0.95LCL…7 | Spring | 92 | 17.4 | 5.15 | 4.41 | 5.12 |
| 0.95UCL…8 | Spring | 92 | 17.4 | 7.47 | 4.41 | 7.61 |
| mean…9 | Summer | 92 | 34.8 | 4.06 | 2.47 | 3.80 |
| se…10 | Summer | 92 | 34.8 | 0.34 | 2.47 | 0.40 |
| 0.95LCL…11 | Summer | 92 | 34.8 | 3.40 | 2.47 | 3.11 |
| 0.95UCL…12 | Summer | 92 | 34.8 | 4.72 | 2.47 | 4.63 |
| mean…13 | Fall | 92 | 6.5 | 10.09 | 6.41 | 9.80 |
| se…14 | Fall | 92 | 6.5 | 1.15 | 6.41 | 1.10 |
| 0.95LCL…15 | Fall | 92 | 6.5 | 7.84 | 6.41 | 7.94 |
| 0.95UCL…16 | Fall | 92 | 6.5 | 12.35 | 6.41 | 12.09 |
Step 3 — Compare seasons with a non-parametric test
# Peto-Peto test across seasons
cendiff(pm25$observed, pm25$censored, groups = factor(pm25$season))## N Observed Expected (O-E)^2/E (O-E)^2/V
## factor(pm25$season)=Winter 89 57.5 31.0 22.69 38.41
## factor(pm25$season)=Spring 92 41.0 48.6 1.18 2.37
## factor(pm25$season)=Summer 92 26.5 60.5 19.05 44.49
## factor(pm25$season)=Fall 92 53.1 38.1 5.89 10.63
##
## Chisq= 72.3 on 3 degrees of freedom, p= 1e-15Step 4 — Exceedance probability
What fraction of days exceed the 35 µg/m³ NAAQS? Because detected values > 35 are always detects, we can compute this directly from the detected data.
exceed_stats <- pm25 |>
group_by(season) |>
summarise(
n_days = n(),
n_exceed = sum(observed > 35, na.rm = TRUE),
pct_exceed = round(100 * mean(observed > 35, na.rm = TRUE), 1),
annual_frac = round(n_exceed / 365 * 100, 1)
)
kable(exceed_stats,
col.names = c("Season", "Days", "Exceedances", "% of season", "% of year"),
caption = "PM2.5 exceedances of the 35 µg/m³ 24-hour NAAQS by season.")| Season | Days | Exceedances | % of season | % of year |
|---|---|---|---|---|
| Winter | 89 | 3 | 3.4 | 0.8 |
| Spring | 92 | 0 | 0.0 | 0.0 |
| Summer | 92 | 0 | 0.0 | 0.0 |
| Fall | 92 | 5 | 5.4 | 1.4 |
pm25 |>
filter(!censored) |>
mutate(season = factor(season, levels = c("Winter","Spring","Summer","Fall"))) |>
ggplot(aes(x = season, y = observed, fill = season)) +
geom_violin(trim = FALSE, alpha = 0.6, linewidth = 0.5) +
geom_boxplot(width = 0.12, outlier.size = 1, fill = "white", alpha = 0.8) +
geom_hline(yintercept = 35, linetype = "dashed", colour = "#c62828", linewidth = 0.9) +
annotate("text", x = 4.4, y = 37, label = "NAAQS 35 µg/m³",
colour = "#c62828", hjust = 1, size = 3.4) +
scale_fill_manual(values = c("#1e3a2f","#4a7c5e","#c8832a","#8aad96"), guide = "none") +
labs(x = NULL, y = expression(PM[2.5]~(µg/m^3)),
title = expression(PM[2.5]~by~season~(detected~values~only))) +
theme_minimal(base_size = 13) +
theme(plot.title = element_text(family = "serif", face = "bold", colour = "#1e3a2f"))
Distribution of PM2.5 concentrations by season. Detected values are shown as violin plots; the number of non-detect days appears at the bottom of each panel. The 35 µg/m³ NAAQS dashed line illustrates that exceedances occur primarily in winter. Summer shows the most non-detects (clean-air days).
9. Hypothesis Tests for Censored Data
Standard t-tests and ANOVA are invalid when a substantial fraction of
observations are non-detects. The NADA2 package provides
censored analogues.
Two-sample test: Lead levels in older vs. newer homes
# Simulate: older homes (pre-1978) have higher Pb than newer homes
set.seed(505)
pb2 <- bind_rows(
sim_censored(30, meanlog = 1.5, sdlog = 1.2, rl = 0.5) |>
mutate(era = "Pre-1978"),
sim_censored(30, meanlog = 0.3, sdlog = 1.1, rl = 0.5) |>
mutate(era = "Post-1978")
)
cat("Pre-1978 non-detects: ", sum(pb2$censored[pb2$era == "Pre-1978"]),
"/30 (", round(100*mean(pb2$censored[pb2$era == "Pre-1978"]),1),"%)\n")## Pre-1978 non-detects: 3 /30 ( 10 %)cat("Post-1978 non-detects:", sum(pb2$censored[pb2$era == "Post-1978"]),
"/30 (", round(100*mean(pb2$censored[pb2$era == "Post-1978"]),1),"%)\n\n")## Post-1978 non-detects: 7 /30 ( 23.3 %)# Censored Wilcoxon / Score test via cendiff()
cendiff(pb2$observed, pb2$censored, groups = factor(pb2$era))## N Observed Expected (O-E)^2/E (O-E)^2/V
## factor(pb2$era)=Post-1978 30 10.3 18.9 3.87 15.5
## factor(pb2$era)=Pre-1978 30 19.2 10.7 6.83 15.5
##
## Chisq= 15.5 on 1 degrees of freedom, p= 8e-05Pairwise comparisons for E. coli across sites
When comparing more than two groups, use cenperm2() from
NADA2 for permutation-based pairwise tests with multiple
comparison adjustment.
# Permutation-based pairwise test (BH adjustment by default)
cenperm2(ecoli$observed, ecoli$censored, factor(ecoli$site))## Permutation test of mean CensData: ecoli$observed by Factor: factor(ecoli$site)
## 9999 Permutations alternative = two.sided
## mean of Beach A = 31.1 to 31.1 mean of Beach B = 17.75 to 17.75
## Mean (Beach A - Beach B) = 13.35 to 13.35 p = 0.3687 to 0.3687
## Correlation with censored data
Does phosphorus increase with turbidity? Standard Pearson r is biased
when either variable contains non-detects. Use Kendall’s τ via
cenken().
set.seed(606)
# Simulate turbidity (NTU) — positively correlated with TP, no non-detects
n_corr <- 40
turb_ntu <- rexp(n_corr, rate = 0.08) + 2
# TP correlated with turbidity, with non-detects
tp_corr <- 3 + 1.2 * turb_ntu + rnorm(n_corr, 0, 4)
obs_tp <- pmax(tp_corr, 0)
nd_tp <- tp_corr < 5
obs_tp[nd_tp] <- 0
# Kendall's tau for censored data
cenken(obs_tp, nd_tp, turb_ntu)## Censored Kendall's Tau Regression:
## Slope : 1.20433
## Intercept : 3.907174
## Tau : 0.7717949
## p-value : 2.420508e-12📌 Censored Kendall’s τ
cenken() computes Akritas–Murphy–LaValley (AML)
Kendall’s τ, a rank-based correlation coefficient that discards
or downweights pairs where at least one member is a non-detect. The test
is valid as long as the censoring fraction is below ~70–80%. The sign
and magnitude of τ are interpreted exactly like ordinary Kendall’s τ:
positive means the two variables tend to increase together.
10. Reporting Non-Detect Statistics
What to always report
Reporting non-detects transparently is both a scientific and a regulatory obligation. Every summary table of censored data should include:
- Sample size (total n)
- Number and fraction of non-detects (e.g., “12/40, 30%”)
- Reporting limit(s) (and whether multiple RLs were present)
- The method used (KM, ROS, MLE, substitution) with a citation
- Estimates with uncertainty (mean ± SE or 95% CI; percentiles)
Complete example: Water quality summary table
# Compute KM statistics for all four contaminants in a unified format
datasets <- list(
"E. coli (Beach A, CFU/100 mL)" = list(obs = beach_a$observed,
cen = beach_a$censored, rl = 1),
"E. coli (Beach C, CFU/100 mL)" = list(obs = filter(ecoli, site=="Beach C")$observed,
cen = filter(ecoli, site=="Beach C")$censored, rl = 1),
"Total P – Surface (µg/L)" = list(obs = surf_tp$observed,
cen = surf_tp$censored, rl = 5),
"Lead Pb (µg/L)" = list(obs = pb$observed,
cen = pb$censored, rl = 0.5),
"PM2.5 – Winter (µg/m³)" = list(obs = filter(pm25, season=="Winter")$observed,
cen = filter(pm25, season=="Winter")$censored, rl = 2)
)
final_tbl <- lapply(names(datasets), function(nm) {
d <- datasets[[nm]]
km <- cenfit(d$obs, d$cen)
data.frame(
Contaminant = nm,
n = length(d$obs),
ND = paste0(sum(d$cen), " (", round(100*mean(d$cen),0), "%)"),
RL = d$rl,
KM_mean = round(mean(km), 2),
KM_median = round(median(km), 2),
Method = "Kaplan–Meier"
)
}) |> bind_rows()
kable(final_tbl,
col.names = c("Contaminant (units)", "n", "Non-detects",
"RL", "Mean", "Median", "Method"),
caption = "Summary statistics for environmental monitoring data with non-detects. All estimates are from Kaplan–Meier using the NADA package (Helsel et al. 2020).") |>
kable_styling(bootstrap_options = "striped")| Contaminant (units) | n | Non-detects | RL | Mean | Median | Method | |
|---|---|---|---|---|---|---|---|
| mean…1 | E. coli (Beach A, CFU/100 mL) | 24 | 0 (0%) | 1.0 | 31.10 | 19.83 | Kaplan–Meier |
| se…2 | E. coli (Beach A, CFU/100 mL) | 24 | 0 (0%) | 1.0 | 5.76 | 19.83 | Kaplan–Meier |
| 0.95LCL…3 | E. coli (Beach A, CFU/100 mL) | 24 | 0 (0%) | 1.0 | 19.81 | 19.83 | Kaplan–Meier |
| 0.95UCL…4 | E. coli (Beach A, CFU/100 mL) | 24 | 0 (0%) | 1.0 | 42.39 | 19.83 | Kaplan–Meier |
| mean…5 | E. coli (Beach C, CFU/100 mL) | 24 | 0 (0%) | 1.0 | 78.80 | 56.41 | Kaplan–Meier |
| se…6 | E. coli (Beach C, CFU/100 mL) | 24 | 0 (0%) | 1.0 | 13.71 | 56.41 | Kaplan–Meier |
| 0.95LCL…7 | E. coli (Beach C, CFU/100 mL) | 24 | 0 (0%) | 1.0 | 51.94 | 56.41 | Kaplan–Meier |
| 0.95UCL…8 | E. coli (Beach C, CFU/100 mL) | 24 | 0 (0%) | 1.0 | 105.66 | 56.41 | Kaplan–Meier |
| mean…9 | Total P – Surface (µg/L) | 30 | 7 (23%) | 5.0 | 13.92 | 9.99 | Kaplan–Meier |
| se…10 | Total P – Surface (µg/L) | 30 | 7 (23%) | 5.0 | 2.57 | 9.99 | Kaplan–Meier |
| 0.95LCL…11 | Total P – Surface (µg/L) | 30 | 7 (23%) | 5.0 | 8.88 | 9.99 | Kaplan–Meier |
| 0.95UCL…12 | Total P – Surface (µg/L) | 30 | 7 (23%) | 5.0 | 18.96 | 9.99 | Kaplan–Meier |
| mean…13 | Lead Pb (µg/L) | 50 | 6 (12%) | 0.5 | 5.23 | 2.38 | Kaplan–Meier |
| se…14 | Lead Pb (µg/L) | 50 | 6 (12%) | 0.5 | 1.25 | 2.38 | Kaplan–Meier |
| 0.95LCL…15 | Lead Pb (µg/L) | 50 | 6 (12%) | 0.5 | 2.78 | 2.38 | Kaplan–Meier |
| 0.95UCL…16 | Lead Pb (µg/L) | 50 | 6 (12%) | 0.5 | 7.68 | 2.38 | Kaplan–Meier |
| mean…17 | PM2.5 – Winter (µg/m³) | 89 | 5 (6%) | 2.0 | 11.25 | 8.51 | Kaplan–Meier |
| se…18 | PM2.5 – Winter (µg/m³) | 89 | 5 (6%) | 2.0 | 0.96 | 8.51 | Kaplan–Meier |
| 0.95LCL…19 | PM2.5 – Winter (µg/m³) | 89 | 5 (6%) | 2.0 | 9.36 | 8.51 | Kaplan–Meier |
| 0.95UCL…20 | PM2.5 – Winter (µg/m³) | 89 | 5 (6%) | 2.0 | 13.14 | 8.51 | Kaplan–Meier |
Common mistakes checklist
cenfit() or
cenros(). Even 15% non-detects can meaningfully bias a mean
when substituting zero.
cendiff() for two-group comparisons or
cenperm2() for multiple groups.
🔑 Key Takeaways
- Non-detects are left-censored observations, not zeros and not missing data. A result of “< RL” carries real information: the true value lies in [0, RL].
- Substitution methods (0, RL/2, RL) are biased whenever the censoring fraction exceeds ~15%. The bias is systematic and worsens with higher censoring.
-
Kaplan–Meier (
cenfit()) is non-parametric, handles multiple reporting limits, and is the safest choice when you cannot verify distributional assumptions. It is the baseline method. -
ROS (
cenros()) produces a complete imputed dataset useful for visualization and regression; valid up to ~80% censoring when data are approximately log-normal. -
MLE (
cenmle()) gives analytical confidence intervals and is most efficient when the distributional assumption is correct. Check with a censored Q–Q plot (cenQQ()). -
Use
cendiff()for group comparisons andcenken()for censored correlations. Never apply standard t-tests or ANOVA directly to substituted data. - Always report: n, number and fraction of non-detects, the reporting limit, and the method used — with a citation to Helsel (2012).
What’s next?
- Survival Analysis — the right-censored twin of this problem; time-to-event methods with the same KM and Cox machinery
- Time Series — detecting trends in long-term monitoring data that also contain non-detects (Mann–Kendall with censoring)
- Multiple Regression & GLMs — building predictive models when the response variable contains non-detects
Further reading
- Helsel, D.R. (2012). Statistics for Censored Environmental Data Using Minitab and R, 2nd ed. Wiley. — The authoritative reference; all methods here follow Helsel’s framework.
- Helsel, D.R., Hirsch, R.M., et al. (2020). Statistical Methods in Water Resources, USGS Techniques and Methods 4–A3. Freely available at pubs.usgs.gov.
- USEPA (2009). Statistical Analysis of Groundwater Monitoring Data at RCRA Facilities: Unified Guidance. EPA 530/R-09-007. — regulatory context for non-detect handling.