Testing for an ecological threshold

hypothesis testing

ecology tutorial

ggplot2

Deciding whether a threshold is real in base R: the naive break test over-rejects, so bootstrap the null, and do not mistake a smooth curve for a kink.

Author

Tidy Ecology

Published

2026-05-31

Modified

2026-09-27

Corrected 27 September 2026: the closing section said the data support the choice between a threshold and a smooth bend only weakly; in the worked example AIC separates them clearly once the smooth curve is fitted, and the section now says so.

The previous post fitted a breakpoint by profiling the residual sum of squares. That procedure always returns a minimum, so a good broken-stick fit is not by itself evidence that a threshold exists. This post asks the harder question: is the break real, or would a straight line have done?

The tempting shortcut is to fit the two-piece model, read the p-value of the slope-change term, and call anything below 0.05 a threshold. That test is not valid. Under the hypothesis of no break the breakpoint is not a parameter of the model at all, so it is estimated from the noise, and the usual F reference does not hold. The fix is to build the null distribution by simulation, and to remember that rejecting a straight line is not the same as finding a kink.

Why the naive p-value is wrong

The statistic itself is fine: at each candidate breakpoint, compare the two-piece fit with the straight line by an F ratio, and take the largest F over all candidates. This maximised value is the natural test statistic.

seg_supF <- function(x, y, grid) {
  fit0 <- lm(y ~ x); rss0 <- sum(fit0$residuals^2); n <- length(y)
  best <- Inf; psi <- NA
  for (p in grid) {
    r <- sum(lm(y ~ x + pmax(x - p, 0))$residuals^2)
    if (r < best) { best <- r; psi <- p }
  }
  list(F = ((rss0 - best) / 1) / (best / (n - 3)), psi = psi)
}
grid <- seq(4, 26, by = 0.2)

The trouble is the reference. An ordinary F statistic for one added term follows an F distribution with one and $n - 3$ degrees of freedom. But this F was maximised over every candidate breakpoint, so it is systematically larger than a single F even when nothing is going on. Simulate straight-line data with no break, run the procedure, and read the naive p-value each time.

set.seed(4199)
Rrep <- 2000; n <- 80
naive_p <- numeric(Rrep)
for (r in seq_len(Rrep)) {
  x <- runif(n, 0, 30)
  y <- 20 - 0.15 * x + rnorm(n, 0, 3)         # a plain line, no break
  naive_p[r] <- pf(seg_supF(x, y, grid)$F, 1, n - 3, lower.tail = FALSE)
}
rej_naive <- mean(naive_p < 0.05)

Histogram of p-values under no break, strongly skewed towards zero rather than uniform, with the region below 0.05 shaded and a dashed line marking the flat height a valid test would show. — Figure 1: Naive p-values from 2000 data sets that contain no break. A valid test would give a flat histogram; instead the mass piles up near zero, and far more than 5% fall below 0.05.

The p-values in Figure 1 are not flat. They crowd towards zero, and about 18% fall below 0.05 rather than the 5% a valid test would give. Read literally, that is how often the naive break test would announce a threshold in data sets that have none.

A null distribution by bootstrap

The repair is to compare the observed statistic against its own null distribution rather than an off-the-shelf F. Fit the straight line, treat that as the truth, simulate new responses from it, and recompute the maximised F each time. Here is a data set with a genuine break to test.

set.seed(781)
n <- 100
x <- runif(n, 0, 30); psi0 <- 14
y <- 18 + 0.10 * x - 0.80 * pmax(x - psi0, 0) + rnorm(n, 0, 2.5)
obs <- seg_supF(x, y, grid)
naive_p_obs <- pf(obs$F, 1, n - 3, lower.tail = FALSE)

fit0 <- lm(y ~ x); sig0 <- summary(fit0)$sigma
set.seed(55)
B <- 1200
null_F <- replicate(B, seg_supF(x, fitted(fit0) + rnorm(n, 0, sig0), grid)$F)
p_boot <- (1 + sum(null_F >= obs$F)) / (B + 1)
crit_boot  <- quantile(null_F, 0.95)
crit_naive <- qf(0.95, 1, n - 3)
size_naive <- mean(null_F > crit_naive)     # the naive test's true size

The maximised F here is 57.4, which the F table would call significant at about $2.1 \times 10^{-11}$. The bootstrap agrees that the break is real but is far more measured: its p-value is 0.001. The difference is entirely in the threshold. A valid 5% test needs an F above 6.6 (the 95th percentile of the null distribution), whereas the naive table uses 3.9. Using that lower bar, the actual rejection rate under no break is 17%, matching the 18% from the direct simulation above. Two routes, one conclusion: the F table is roughly three to four times too generous.

A significant test is not a threshold

There is a deeper limit, and it is the honest note of the cluster. The break test asks whether a straight line is adequate. It does not ask whether the true shape is a kink. A smooth curve that bends also fails the straight-line test, and it will fail it emphatically.

set.seed(4299)
n <- 90
x <- runif(n, 0, 30)
y <- 30 * x / (6 + x) + rnorm(n, 0, 2.2)        # smooth saturation, no kink
sm <- seg_supF(x, y, grid)
naive_p_smooth <- pf(sm$F, 1, n - 3, lower.tail = FALSE)

fit_seg <- lm(y ~ x + pmax(x - sm$psi, 0))
fit_mm  <- nls(y ~ SSmicmen(x, Vm, K))          # the smooth alternative
aic_seg <- AIC(fit_seg) + 2                      # +2 for the estimated breakpoint
aic_mm  <- AIC(fit_mm)

Scatter rising and levelling off; a solid smooth saturating curve and a dashed broken-stick line lie almost on top of each other through the points. — Figure 2: Smooth saturating data. The break test rejects a straight line overwhelmingly, yet a smooth curve (solid) fits as well as the broken stick (dashed): there is no kink to find.

The data in Figure 2 are a smooth saturating curve with no breakpoint anywhere. The maximised F is 122, and the naive p-value is about $3 \times 10^{-18}$: an overwhelming vote against the straight line. A practitioner who compares only a line with a broken stick will report a sharp threshold. But when the smooth saturating curve is fitted as well, AIC prefers it by about 12.4 points, revealing the apparent break as a gentle bend.

The honest limit

A break test controls a single decision: is a straight line enough? Do not run it with the F table, which is calibrated for a fixed breakpoint and rejects three to four times too often once the breakpoint is estimated; build the null distribution by simulation instead. And once a line is rejected, the work is not finished. Ask AIC, or a smooth alternative, whether the response truly turns at a point or simply curves. Distinguishing a genuine threshold from a smooth bend is a choice between shapes, and the break test alone cannot make it: it takes a smooth alternative fitted alongside, as here, where AIC separated the two by about 12 points. The next post carries the threshold idea into count and presence data, where the same distinction between a step and a slope has to be made explicitly.

References

Toms JD, Lesperance ML 2003. Ecology 84(8):2034-2041 (10.1890/02-0472)
Toms JD, Villard MA 2015. Avian Conservation and Ecology 10(1):2 (10.5751/ACE-00715-100102)
Davies RB 1987. Biometrika 74(1):33-43 (10.1093/biomet/74.1.33)
Andrews DWK 1993. Econometrica 61(4):821-856 (10.2307/2951764)
Muggeo VMR 2003. Statistics in Medicine 22(19):3055-3071 (10.1002/sim.1545)

--- title: "Testing for an ecological threshold" description: "Deciding whether a threshold is real in base R: the naive break test over-rejects, so bootstrap the null, and do not mistake a smooth curve for a kink." date: "2026-05-31 11:00" categories: [R, hypothesis testing, ecology tutorial, ggplot2] image: thumbnail.png date-modified: "2026-09-27" image-alt: "Saturating data with a smooth curve and a broken-stick fit lying almost on top of each other." --- *Corrected 27 September 2026: the closing section said the data support the choice between a threshold and a smooth bend only weakly; in the worked example AIC separates them clearly once the smooth curve is fitted, and the section now says so.* ```{r setup} #| include: false library(ggplot2) te_ink <- "#16241d"; te_body <- "#2c3a31"; te_forest <- "#275139" te_label <- "#46604a"; te_sage <- "#93a87f"; te_paper <- "#f5f4ee" te_line <- "#dad9ca"; te_faint <- "#5d6b61"; te_gold <- "#cda23f" te_green <- "#2f8f63"; te_terra <- "#b5534e" theme_te <- function(base_size = 12) { theme_minimal(base_size = base_size) + 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(), axis.text = element_text(colour = te_faint, size = base_size - 2), axis.title = element_text(colour = te_body, size = base_size - 1), plot.title = element_text(colour = te_ink, face = "bold", size = base_size + 1), plot.subtitle = element_text(colour = te_label, size = base_size - 1), plot.caption = element_text(colour = te_faint, size = base_size - 3), legend.position = "top", plot.margin = margin(10, 12, 8, 10)) } ``` ```{r sci-tex-helper} #| include: false sci_tex <- function(x, digits = 2) { s <- formatC(x, format = "e", digits = digits) paste0("$", sub("e.*", "", s), " \\times 10^{", as.integer(sub(".*e", "", s)), "}$") } ``` The [previous post](../segmented-regression-breakpoints/) fitted a breakpoint by profiling the residual sum of squares. That procedure always returns a minimum, so a good broken-stick fit is not by itself evidence that a threshold exists. This post asks the harder question: is the break real, or would a straight line have done? The tempting shortcut is to fit the two-piece model, read the p-value of the slope-change term, and call anything below 0.05 a threshold. That test is not valid. Under the hypothesis of no break the breakpoint is not a parameter of the model at all, so it is estimated from the noise, and the usual F reference does not hold. The fix is to build the null distribution by simulation, and to remember that rejecting a straight line is not the same as finding a kink. ## Why the naive p-value is wrong The statistic itself is fine: at each candidate breakpoint, compare the two-piece fit with the straight line by an F ratio, and take the largest F over all candidates. This maximised value is the natural test statistic. ```{r supf} seg_supF <- function(x, y, grid) { fit0 <- lm(y ~ x); rss0 <- sum(fit0$residuals^2); n <- length(y) best <- Inf; psi <- NA for (p in grid) { r <- sum(lm(y ~ x + pmax(x - p, 0))$residuals^2) if (r < best) { best <- r; psi <- p } } list(F = ((rss0 - best) / 1) / (best / (n - 3)), psi = psi) } grid <- seq(4, 26, by = 0.2) ``` The trouble is the reference. An ordinary F statistic for one added term follows an F distribution with one and $n - 3$ degrees of freedom. But this F was *maximised* over every candidate breakpoint, so it is systematically larger than a single F even when nothing is going on. Simulate straight-line data with no break, run the procedure, and read the naive p-value each time. ```{r typeI} set.seed(4199) Rrep <- 2000; n <- 80 naive_p <- numeric(Rrep) for (r in seq_len(Rrep)) { x <- runif(n, 0, 30) y <- 20 - 0.15 * x + rnorm(n, 0, 3) # a plain line, no break naive_p[r] <- pf(seg_supF(x, y, grid)$F, 1, n - 3, lower.tail = FALSE) } rej_naive <- mean(naive_p < 0.05) ``` ```{r fig-pvalues} #| echo: false #| fig-cap: "Naive p-values from 2000 data sets that contain no break. A valid test would give a flat histogram; instead the mass piles up near zero, and far more than 5% fall below 0.05." #| fig-alt: "Histogram of p-values under no break, strongly skewed towards zero rather than uniform, with the region below 0.05 shaded and a dashed line marking the flat height a valid test would show." pd <- data.frame(p = naive_p) ggplot(pd, aes(p)) + geom_histogram(breaks = seq(0, 1, 0.05), fill = te_sage, colour = te_paper, linewidth = 0.3) + geom_hline(yintercept = Rrep / 20, colour = te_forest, linewidth = 0.8, linetype = "22") + annotate("rect", xmin = 0, xmax = 0.05, ymin = 0, ymax = Inf, fill = te_terra, alpha = 0.18) + annotate("text", x = 0.5, y = Rrep / 20 + 30, label = "uniform (a valid test)", colour = te_forest, size = 3.4, hjust = 0.5) + annotate("text", x = 0.08, y = max(table(cut(naive_p, seq(0, 1, 0.05)))) * 0.7, hjust = 0, label = sprintf("%.0f%% reject at 0.05", 100 * rej_naive), colour = te_terra, size = 3.4) + labs(x = "Naive p-value for a slope change (data have no break)", y = "Count", title = "The naive break test rejects far too often", subtitle = "Under no break the p-values should be flat; instead they pile up near zero") + theme_te() ``` The p-values in @fig-pvalues are not flat. They crowd towards zero, and about `r sprintf("%.0f%%", 100 * rej_naive)` fall below 0.05 rather than the 5% a valid test would give. Read literally, that is how often the naive break test would announce a threshold in data sets that have none. ## A null distribution by bootstrap The repair is to compare the observed statistic against its own null distribution rather than an off-the-shelf F. Fit the straight line, treat that as the truth, simulate new responses from it, and recompute the maximised F each time. Here is a data set with a genuine break to test. ```{r break-example} set.seed(781) n <- 100 x <- runif(n, 0, 30); psi0 <- 14 y <- 18 + 0.10 * x - 0.80 * pmax(x - psi0, 0) + rnorm(n, 0, 2.5) obs <- seg_supF(x, y, grid) naive_p_obs <- pf(obs$F, 1, n - 3, lower.tail = FALSE) fit0 <- lm(y ~ x); sig0 <- summary(fit0)$sigma set.seed(55) B <- 1200 null_F <- replicate(B, seg_supF(x, fitted(fit0) + rnorm(n, 0, sig0), grid)$F) p_boot <- (1 + sum(null_F >= obs$F)) / (B + 1) crit_boot <- quantile(null_F, 0.95) crit_naive <- qf(0.95, 1, n - 3) size_naive <- mean(null_F > crit_naive) # the naive test's true size ``` The maximised F here is `r round(obs$F, 1)`, which the F table would call significant at about `r sci_tex(naive_p_obs, 1)`. The bootstrap agrees that the break is real but is far more measured: its p-value is `r round(p_boot, 3)`. The difference is entirely in the threshold. A valid 5% test needs an F above `r round(crit_boot, 1)` (the 95th percentile of the null distribution), whereas the naive table uses `r round(crit_naive, 1)`. Using that lower bar, the actual rejection rate under no break is `r sprintf("%.0f%%", 100 * size_naive)`, matching the `r sprintf("%.0f%%", 100 * rej_naive)` from the direct simulation above. Two routes, one conclusion: the F table is roughly three to four times too generous. ## A significant test is not a threshold There is a deeper limit, and it is the honest note of the cluster. The break test asks whether a straight line is adequate. It does not ask whether the true shape is a kink. A smooth curve that bends also fails the straight-line test, and it will fail it emphatically. ```{r smooth} set.seed(4299) n <- 90 x <- runif(n, 0, 30) y <- 30 * x / (6 + x) + rnorm(n, 0, 2.2) # smooth saturation, no kink sm <- seg_supF(x, y, grid) naive_p_smooth <- pf(sm$F, 1, n - 3, lower.tail = FALSE) fit_seg <- lm(y ~ x + pmax(x - sm$psi, 0)) fit_mm <- nls(y ~ SSmicmen(x, Vm, K)) # the smooth alternative aic_seg <- AIC(fit_seg) + 2 # +2 for the estimated breakpoint aic_mm <- AIC(fit_mm) ``` ```{r fig-smooth} #| echo: false #| fig-cap: "Smooth saturating data. The break test rejects a straight line overwhelmingly, yet a smooth curve (solid) fits as well as the broken stick (dashed): there is no kink to find." #| fig-alt: "Scatter rising and levelling off; a solid smooth saturating curve and a dashed broken-stick line lie almost on top of each other through the points." cf <- coef(fit_seg); gx <- seq(0, 30, length.out = 400) pred <- data.frame(gx, seg = cf[[1]] + cf[[2]] * gx + cf[[3]] * pmax(gx - sm$psi, 0), mm = predict(fit_mm, newdata = data.frame(x = gx))) ggplot(data.frame(x, y), aes(x, y)) + geom_point(colour = te_forest, alpha = 0.5, size = 1.7) + geom_line(data = pred, aes(gx, seg), colour = te_terra, linewidth = 1.1, linetype = "22") + geom_line(data = pred, aes(gx, mm), colour = te_forest, linewidth = 1.1) + geom_vline(xintercept = sm$psi, colour = te_gold, linewidth = 0.6) + annotate("text", x = 16, y = 8, hjust = 0, label = "broken stick (dashed)", colour = te_terra, size = 3.4) + annotate("text", x = 16, y = 5, hjust = 0, label = "smooth curve (solid)", colour = te_forest, size = 3.4) + labs(x = "Environmental gradient", y = "Response", title = "A significant break test is not a threshold", subtitle = "Smooth saturating data: the test rejects a line, but there is no kink") + theme_te() ``` The data in @fig-smooth are a smooth saturating curve with no breakpoint anywhere. The maximised F is `r round(sm$F, 0)`, and the naive p-value is about `r sci_tex(naive_p_smooth, 0)`: an overwhelming vote against the straight line. A practitioner who compares only a line with a broken stick will report a sharp threshold. But when the smooth saturating curve is fitted as well, AIC prefers it by about `r round(aic_seg - aic_mm, 1)` points, revealing the apparent break as a gentle bend. ## The honest limit A break test controls a single decision: is a straight line enough? Do not run it with the F table, which is calibrated for a fixed breakpoint and rejects three to four times too often once the breakpoint is estimated; build the null distribution by simulation instead. And once a line is rejected, the work is not finished. Ask AIC, or a smooth alternative, whether the response truly turns at a point or simply curves. Distinguishing a genuine threshold from a smooth bend is a choice between shapes, and the break test alone cannot make it: it takes a smooth alternative fitted alongside, as here, where AIC separated the two by about `r round(aic_seg - aic_mm, 0)` points. The next post carries the threshold idea into count and presence data, where the same distinction between a step and a slope has to be made explicitly. ## References - Toms JD, Lesperance ML 2003. Ecology 84(8):2034-2041 (10.1890/02-0472) - Toms JD, Villard MA 2015. Avian Conservation and Ecology 10(1):2 (10.5751/ACE-00715-100102) - Davies RB 1987. Biometrika 74(1):33-43 (10.1093/biomet/74.1.33) - Andrews DWK 1993. Econometrica 61(4):821-856 (10.2307/2951764) - Muggeo VMR 2003. Statistics in Medicine 22(19):3055-3071 (10.1002/sim.1545) ## Related tutorials - [Segmented regression for a breakpoint](../segmented-regression-breakpoints/) - [Threshold regression with a GLM](../threshold-regression-glm/) - [Checking a threshold model](../checking-a-threshold-model/) - [Derivatives of a GAM trend: when did it change?](../derivatives-of-a-smooth-trend/)