Checking a threshold model

model diagnostics

ecology tutorial

ggplot2

Checking a fitted threshold in base R: read the residuals, see why the Wald breakpoint interval undercovers, and test how far the threshold shifts under a different assumed shape.

Author

Tidy Ecology

Published

2026-08-05

A fitted threshold model is persuasive. It reports a breakpoint to two decimals, draws two clean lines, and lowers the AIC against a straight line. This closing post is about the checks that stand between that output and a claim. There are three: whether the two-piece shape actually fits, whether the interval around the breakpoint is honest, and how much the whole story depends on assumptions the data cannot settle. Each mirrors the diagnostic that closed the earlier method clusters, and the last is the familiar limit of this one.

Do the residuals still carry a pattern?

The first check is the cheapest. If a straight line is wrong and a broken stick is right, the line will leave structure in its residuals and the broken stick will not. Fit both and look.

set.seed(4231)
n <- 90
x <- runif(n, 0, 30); psi0 <- 14
y <- 18 + 0.10 * x - 0.55 * pmax(x - psi0, 0) + rnorm(n, 0, 2.5)
grid <- seq(4, 26, by = 0.2)
psi_hat <- grid[which.min(vapply(grid, function(p) sum(lm(y ~ x + pmax(x - p, 0))$residuals^2), numeric(1)))]
fit_line <- lm(y ~ x)
fit_thr  <- lm(y ~ x + pmax(x - psi_hat, 0))

line_curv <- summary(lm(residuals(fit_line) ~ poly(fitted(fit_line), 2)))$coefficients[3, 4]
thr_curv  <- summary(lm(residuals(fit_thr)  ~ poly(fitted(fit_thr),  2)))$coefficients[3, 4]

`geom_smooth()` using formula = 'y ~ x'

Two residual panels; the single-line residuals bow into a curve highlighted by a smoother, the broken-stick residuals scatter flat around zero. — Figure 1: Residuals against fitted values. The single line leaves a curve (a quadratic trend is clear); the broken stick leaves none.

A quadratic fitted to the line’s residuals has a highly significant curvature term (p 3.5^{-4}), the fingerprint of a bend the line cannot follow. The same term on the broken stick’s residuals is not significant (p 0.17): the two-piece shape has absorbed the structure (Figure 1). If the broken stick’s residuals still curved, that would be a warning that the true shape is smooth rather than kinked, the distinction from the previous post.

Is the breakpoint interval honest?

The breakpoint comes with a standard error, and the quick one is not trustworthy. The usual Wald interval treats the breakpoint like any coefficient, but its sampling distribution is skewed and its quick standard error is too small, so the interval is too narrow. A profile interval, the set of breakpoints whose fit is statistically indistinguishable from the best, is wider and closer to honest.

fit_bs <- function(x, y, grid) {
  n <- length(y)
  rss <- vapply(grid, function(p) sum(lm(y ~ x + pmax(x - p, 0))$residuals^2), numeric(1))
  psi <- grid[which.min(rss)]; rss_min <- min(rss); sig2 <- rss_min / (n - 3)
  aug <- lm(y ~ x + pmax(x - psi, 0) + as.numeric(x > psi))    # Muggeo's augmented fit
  se_psi <- summary(aug)$coefficients[4, 2] / abs(coef(aug)[[3]])
  wald <- psi + c(-1, 1) * 1.96 * se_psi
  prof <- range(grid[rss <= rss_min + sig2 * qchisq(0.95, 1)])
  list(wald = wald, prof = prof)
}

set.seed(4201)
Rrep <- 1500; n <- 80; psi0 <- 14
cov_wald <- 0; cov_prof <- 0
for (r in seq_len(Rrep)) {
  x <- runif(n, 0, 30)
  y <- 18 + 0.10 * x - 0.60 * pmax(x - psi0, 0) + rnorm(n, 0, 3)
  f <- fit_bs(x, y, grid)
  cov_wald <- cov_wald + (psi0 >= f$wald[1] && psi0 <= f$wald[2])
  cov_prof <- cov_prof + (psi0 >= f$prof[1] && psi0 <= f$prof[2])
}
cov_wald <- cov_wald / Rrep; cov_prof <- cov_prof / Rrep

Figure 2: Breakpoint intervals from 40 data sets, sorted. Red intervals miss the true breakpoint (gold line at 14). The Wald intervals miss more often than the profile intervals.

Across 1500 data sets the Wald interval covers the true breakpoint only 79% of the time, well short of the nominal 95%, while the profile interval reaches 93% (Figure 2). The Wald interval is not merely imperfect; it is confidently wrong about a value it barely constrains. Report the profile interval, or a bootstrap interval as in the first post, and never the quick standard error alone.

How much does the threshold depend on assumptions?

The last check is the honest limit of the cluster. A breakpoint is reported as a number, but that number is conditional on choices the data support only weakly: the assumed shape, the number of breakpoints, and the noise. Change any of them and the threshold moves.

set.seed(9014)
n <- 120
x <- runif(n, 0, 30); psi0 <- 14
y <- 18 + 0.10 * x - 0.70 * pmax(x - psi0, 0) + rnorm(n, 0, 2.5)
psi_cont <- grid[which.min(vapply(grid, function(p) sum(lm(y ~ x + pmax(x - p, 0))$residuals^2), numeric(1)))]
psi_step <- grid[which.min(vapply(grid, function(p) sum(lm(y ~ x + as.numeric(x > p))$residuals^2), numeric(1)))]

est_psi <- function(xx, yy) {
  g <- seq(4, 26, by = 0.25)
  g[which.min(vapply(g, function(p) sum(lm(yy ~ xx + pmax(xx - p, 0))$residuals^2), numeric(1)))]
}
set.seed(31)
boot_psi <- replicate(1500, { i <- sample.int(n, n, replace = TRUE); est_psi(x[i], y[i]) })
bci <- quantile(boot_psi, c(0.025, 0.975))

On one data set, a continuous slope-change model puts the threshold at 15.4, near the true 14. Assume instead that the response steps to a new level, and the same data place the threshold at 20, more than four units away. The noise adds its own width: the bootstrap interval for the continuous estimate runs from 12.5 to 17.8. So a single reported breakpoint stands for a range that depends on the assumed shape and on chance, and the number of breakpoints (from the previous post) is a further choice on top.

The honest limit

A threshold model does not read a breakpoint out of nature. It fits the breakpoint of a shape you chose, with an interval that only the slower methods report honestly, and a location that shifts when the shape or the count changes. The disciplined summary is a checked model: residuals without remaining structure, a profile or bootstrap interval rather than the quick one, and a sensitivity range across the shapes and counts the data cannot rule out. That last step is the same move made throughout these clusters, where an unverifiable assumption is turned into a stated magnitude rather than hidden inside a point estimate: the E-value for an unmeasured confounder, the delta shift for missing-data mechanisms, the reference frame for compositional change, the extrapolation for a nonlinear form. A threshold is not a fact the data hand over; it is a claim, conditional on a shape, and worth only as much as the checks behind it.

References

Muggeo VMR 2017. Australian and New Zealand Journal of Statistics 59(3):311-322 (10.1111/anzs.12200)
Muggeo VMR 2003. Statistics in Medicine 22(19):3055-3071 (10.1002/sim.1545)
Toms JD, Lesperance ML 2003. Ecology 84(8):2034-2041 (10.1890/02-0472)
Ficetola GF, Denoel M 2009. Ecography 32(6):1075-1084 (10.1111/j.1600-0587.2009.05571.x)
Bolker BM 2008. Ecological Models and Data in R. Princeton University Press (ISBN 978-0-691-12522-0)

--- title: "Checking a threshold model" description: "Checking a fitted threshold in base R: read the residuals, see why the Wald breakpoint interval undercovers, and test how far the threshold shifts under a different assumed shape." date: "2026-08-05 12:00" categories: [R, model diagnostics, ecology tutorial, ggplot2] image: thumbnail.png image-alt: "Histogram of bootstrap breakpoint estimates with two vertical lines marking the estimates from a slope-change and a step model." --- ```{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)) } ``` A fitted threshold model is persuasive. It reports a breakpoint to two decimals, draws two clean lines, and lowers the AIC against a straight line. This closing post is about the checks that stand between that output and a claim. There are three: whether the two-piece shape actually fits, whether the interval around the breakpoint is honest, and how much the whole story depends on assumptions the data cannot settle. Each mirrors the diagnostic that closed the earlier method clusters, and the last is the familiar limit of this one. ## Do the residuals still carry a pattern? The first check is the cheapest. If a straight line is wrong and a broken stick is right, the line will leave structure in its residuals and the broken stick will not. Fit both and look. ```{r residuals} set.seed(4231) n <- 90 x <- runif(n, 0, 30); psi0 <- 14 y <- 18 + 0.10 * x - 0.55 * pmax(x - psi0, 0) + rnorm(n, 0, 2.5) grid <- seq(4, 26, by = 0.2) psi_hat <- grid[which.min(vapply(grid, function(p) sum(lm(y ~ x + pmax(x - p, 0))$residuals^2), numeric(1)))] fit_line <- lm(y ~ x) fit_thr <- lm(y ~ x + pmax(x - psi_hat, 0)) line_curv <- summary(lm(residuals(fit_line) ~ poly(fitted(fit_line), 2)))$coefficients[3, 4] thr_curv <- summary(lm(residuals(fit_thr) ~ poly(fitted(fit_thr), 2)))$coefficients[3, 4] ``` ```{r fig-residuals} #| echo: false #| fig-cap: "Residuals against fitted values. The single line leaves a curve (a quadratic trend is clear); the broken stick leaves none." #| fig-alt: "Two residual panels; the single-line residuals bow into a curve highlighted by a smoother, the broken-stick residuals scatter flat around zero." rd <- rbind( data.frame(fit = fitted(fit_line), res = residuals(fit_line), model = "single line"), data.frame(fit = fitted(fit_thr), res = residuals(fit_thr), model = "broken stick")) rd$model <- factor(rd$model, levels = c("single line", "broken stick")) ggplot(rd, aes(fit, res)) + geom_hline(yintercept = 0, colour = te_faint, linewidth = 0.4) + geom_point(colour = te_forest, alpha = 0.5, size = 1.5) + geom_smooth(method = "loess", se = FALSE, colour = te_terra, linewidth = 0.9, span = 0.9) + facet_wrap(~ model) + labs(x = "Fitted value", y = "Residual", title = "Do the residuals still carry a pattern?", subtitle = "The single line leaves a curve in its residuals; the broken stick does not") + theme_te() + theme(strip.text = element_text(colour = te_label, face = "bold")) ``` A quadratic fitted to the line's residuals has a highly significant curvature term (p `r signif(line_curv, 2)`), the fingerprint of a bend the line cannot follow. The same term on the broken stick's residuals is not significant (p `r round(thr_curv, 2)`): the two-piece shape has absorbed the structure (@fig-residuals). If the broken stick's residuals still curved, that would be a warning that the true shape is smooth rather than kinked, the distinction from the previous post. ## Is the breakpoint interval honest? The breakpoint comes with a standard error, and the quick one is not trustworthy. The usual Wald interval treats the breakpoint like any coefficient, but its sampling distribution is skewed and its quick standard error is too small, so the interval is too narrow. A profile interval, the set of breakpoints whose fit is statistically indistinguishable from the best, is wider and closer to honest. ```{r fit-bs} fit_bs <- function(x, y, grid) { n <- length(y) rss <- vapply(grid, function(p) sum(lm(y ~ x + pmax(x - p, 0))$residuals^2), numeric(1)) psi <- grid[which.min(rss)]; rss_min <- min(rss); sig2 <- rss_min / (n - 3) aug <- lm(y ~ x + pmax(x - psi, 0) + as.numeric(x > psi)) # Muggeo's augmented fit se_psi <- summary(aug)$coefficients[4, 2] / abs(coef(aug)[[3]]) wald <- psi + c(-1, 1) * 1.96 * se_psi prof <- range(grid[rss <= rss_min + sig2 * qchisq(0.95, 1)]) list(wald = wald, prof = prof) } ``` ```{r coverage} set.seed(4201) Rrep <- 1500; n <- 80; psi0 <- 14 cov_wald <- 0; cov_prof <- 0 for (r in seq_len(Rrep)) { x <- runif(n, 0, 30) y <- 18 + 0.10 * x - 0.60 * pmax(x - psi0, 0) + rnorm(n, 0, 3) f <- fit_bs(x, y, grid) cov_wald <- cov_wald + (psi0 >= f$wald[1] && psi0 <= f$wald[2]) cov_prof <- cov_prof + (psi0 >= f$prof[1] && psi0 <= f$prof[2]) } cov_wald <- cov_wald / Rrep; cov_prof <- cov_prof / Rrep ``` ```{r fig-coverage} #| echo: false #| fig-cap: "Breakpoint intervals from 40 data sets, sorted. Red intervals miss the true breakpoint (gold line at 14). The Wald intervals miss more often than the profile intervals." #| fig-alt: "Two panels of horizontal interval bars against a gold vertical line at 14; more of the Wald bars are red (missing the line) than the profile bars." rows <- list() for (r in seq_len(40)) { x <- runif(n, 0, 30); y <- 18 + 0.10 * x - 0.60 * pmax(x - psi0, 0) + rnorm(n, 0, 3) fr <- fit_bs(x, y, grid); pj <- grid[which.min(vapply(grid, function(p) sum(lm(y ~ x + pmax(x - p, 0))$residuals^2), numeric(1)))] rows[[length(rows) + 1]] <- data.frame(id = r, psi = pj, method = "Wald", lo = fr$wald[1], hi = fr$wald[2], cover = psi0 >= fr$wald[1] && psi0 <= fr$wald[2]) rows[[length(rows) + 1]] <- data.frame(id = r, psi = pj, method = "profile", lo = fr$prof[1], hi = fr$prof[2], cover = psi0 >= fr$prof[1] && psi0 <= fr$prof[2]) } cat_df <- do.call(rbind, rows) ord <- cat_df$id[cat_df$method == "Wald"][order(cat_df$psi[cat_df$method == "Wald"])] cat_df$rank <- match(cat_df$id, ord) cat_df$method <- factor(cat_df$method, levels = c("Wald", "profile")) ggplot(cat_df, aes(y = rank)) + geom_vline(xintercept = psi0, colour = te_gold, linewidth = 0.7) + geom_segment(aes(x = lo, xend = hi, yend = rank, colour = cover), linewidth = 0.7) + facet_wrap(~ method) + scale_colour_manual(values = c(`TRUE` = te_green, `FALSE` = te_terra), labels = c(`TRUE` = "covers", `FALSE` = "misses"), name = NULL) + labs(x = "Breakpoint interval", y = "Data set (sorted)", title = "Wald intervals miss the breakpoint too often", subtitle = "40 data sets; gold line is the true breakpoint at 14") + theme_te() + theme(strip.text = element_text(colour = te_label, face = "bold"), axis.text.y = element_blank()) ``` Across 1500 data sets the Wald interval covers the true breakpoint only `r sprintf("%.0f%%", 100 * cov_wald)` of the time, well short of the nominal 95%, while the profile interval reaches `r sprintf("%.0f%%", 100 * cov_prof)` (@fig-coverage). The Wald interval is not merely imperfect; it is confidently wrong about a value it barely constrains. Report the profile interval, or a bootstrap interval as in the [first post](../segmented-regression-breakpoints/), and never the quick standard error alone. ## How much does the threshold depend on assumptions? The last check is the honest limit of the cluster. A breakpoint is reported as a number, but that number is conditional on choices the data support only weakly: the assumed shape, the number of breakpoints, and the noise. Change any of them and the threshold moves. ```{r sensitivity} set.seed(9014) n <- 120 x <- runif(n, 0, 30); psi0 <- 14 y <- 18 + 0.10 * x - 0.70 * pmax(x - psi0, 0) + rnorm(n, 0, 2.5) psi_cont <- grid[which.min(vapply(grid, function(p) sum(lm(y ~ x + pmax(x - p, 0))$residuals^2), numeric(1)))] psi_step <- grid[which.min(vapply(grid, function(p) sum(lm(y ~ x + as.numeric(x > p))$residuals^2), numeric(1)))] est_psi <- function(xx, yy) { g <- seq(4, 26, by = 0.25) g[which.min(vapply(g, function(p) sum(lm(yy ~ xx + pmax(xx - p, 0))$residuals^2), numeric(1)))] } set.seed(31) boot_psi <- replicate(1500, { i <- sample.int(n, n, replace = TRUE); est_psi(x[i], y[i]) }) bci <- quantile(boot_psi, c(0.025, 0.975)) ``` On one data set, a continuous slope-change model puts the threshold at `r round(psi_cont, 1)`, near the true `r psi0`. Assume instead that the response steps to a new level, and the same data place the threshold at `r round(psi_step, 1)`, more than four units away. The noise adds its own width: the bootstrap interval for the continuous estimate runs from `r round(bci[[1]], 1)` to `r round(bci[[2]], 1)`. So a single reported breakpoint stands for a range that depends on the assumed shape and on chance, and the number of breakpoints (from the [previous post](../threshold-regression-glm/)) is a further choice on top. ## The honest limit A threshold model does not read a breakpoint out of nature. It fits the breakpoint of a shape you chose, with an interval that only the slower methods report honestly, and a location that shifts when the shape or the count changes. The disciplined summary is a checked model: residuals without remaining structure, a profile or bootstrap interval rather than the quick one, and a sensitivity range across the shapes and counts the data cannot rule out. That last step is the same move made throughout these clusters, where an unverifiable assumption is turned into a stated magnitude rather than hidden inside a point estimate: the E-value for an unmeasured confounder, the delta shift for missing-data mechanisms, the reference frame for compositional change, the extrapolation for a nonlinear form. A threshold is not a fact the data hand over; it is a claim, conditional on a shape, and worth only as much as the checks behind it. ## References - Muggeo VMR 2017. Australian and New Zealand Journal of Statistics 59(3):311-322 (10.1111/anzs.12200) - Muggeo VMR 2003. Statistics in Medicine 22(19):3055-3071 (10.1002/sim.1545) - Toms JD, Lesperance ML 2003. Ecology 84(8):2034-2041 (10.1890/02-0472) - Ficetola GF, Denoel M 2009. Ecography 32(6):1075-1084 (10.1111/j.1600-0587.2009.05571.x) - Bolker BM 2008. Ecological Models and Data in R. Princeton University Press (ISBN 978-0-691-12522-0) ## Related tutorials - [Segmented regression for a breakpoint](../segmented-regression-breakpoints/) - [Testing for an ecological threshold](../detecting-ecological-thresholds/) - [Threshold regression with a GLM](../threshold-regression-glm/) - [GLM residual diagnostics](../glm-residual-diagnostics/)