Segmented regression for a breakpoint

regression

ecology tutorial

ggplot2

Fit a broken-stick model to an ecological threshold in base R: profile the residual sum of squares to place the breakpoint, and see why one straight line misleads.

Author

Tidy Ecology

Published

2026-08-05

Many ecological responses do not bend gently. A grassland loses forbs once nitrogen deposition passes a level the soil can buffer; a stream community collapses above a sediment load; edge effects fade at a fixed distance into a forest. In each case the relationship is roughly linear on either side of a breakpoint, and the breakpoint itself is the quantity of interest. A single straight line cannot represent that shape, and fitting one does not just lose detail: it produces a slope that describes neither side.

This post fits a two-piece (“broken-stick”) model in base R, locates the breakpoint by profiling the residual sum of squares, and reads off the slope on each side. The honest limit at the end is the theme of the whole cluster: a breakpoint is a feature of the model you chose, and its location carries real uncertainty.

A threshold response

Here is a simulated richness response to nitrogen deposition. Below a threshold the community is stable; above it, richness falls.

set.seed(4198)
n <- 120
nitrogen <- runif(n, 0, 30)             # kg N/ha/yr
psi_true <- 12                          # the real breakpoint
richness <- 25 + 0.10 * nitrogen - 0.90 * pmax(nitrogen - psi_true, 0) +
  rnorm(n, 0, 2.5)
dat <- data.frame(nitrogen, richness)

The generating slope is a near-flat 0.10 below the threshold and a declining -0.80 above it, with the join at 12.

Why one line misleads

The reflex is a single linear model.

fit_lin <- lm(richness ~ nitrogen, data = dat)
lin_slope <- coef(fit_lin)[["nitrogen"]]
lin_r2 <- summary(fit_lin)$r.squared

That line has a slope of -0.51 and an R-squared of 0.59. The number looks like a tidy answer, but it is an average of two regimes that do not exist together anywhere in the data. Below the threshold the response is flat; above it, it drops almost twice as steeply as the single slope suggests. Reporting -0.51 as “the effect of nitrogen” describes a mixture, not a mechanism.

Fitting a broken stick

A continuous two-piece model adds one hinge term. For a candidate breakpoint $\psi$, the term $(x - \psi)_+$ equals $x - \psi$ above $\psi$ and zero below it, so the model

\[y = \beta_0 + \beta_1 x + \beta_2 (x - \psi)_+ + \varepsilon\]

is ordinary least squares given $\psi$: the slope is $\beta_1$ below the breakpoint and $\beta_1 + \beta_2$ above it. The one awkward parameter is $\psi$, which enters non-linearly. The direct approach is to profile it: fit the linear model at every candidate breakpoint and keep the value that minimises the residual sum of squares.

rss_at <- function(p) {
  hinge <- pmax(dat$nitrogen - p, 0)
  sum(lm(richness ~ nitrogen + hinge, data = dat)$residuals^2)
}
grid <- seq(3, 27, by = 0.1)
rss_grid <- vapply(grid, rss_at, numeric(1))
psi_hat <- optimize(rss_at, lower = grid[which.min(rss_grid)] - 0.2,
                    upper = grid[which.min(rss_grid)] + 0.2)$minimum

hinge <- pmax(dat$nitrogen - psi_hat, 0)
fit_seg <- lm(richness ~ nitrogen + hinge, data = dat)
cf <- coef(fit_seg)
below_slope <- cf[["nitrogen"]]
above_slope <- cf[["nitrogen"]] + cf[["hinge"]]
seg_r2 <- summary(fit_seg)$r.squared

Figure 1: Residual sum of squares against the candidate breakpoint. The fit is tightest at the estimate (gold); the dashed line marks the true value.

The profile in Figure 1 has a single clear valley. Its minimum puts the breakpoint at 13 kg N/ha/yr, close to the true 12. The fitted slopes are 0.17 below the breakpoint and -0.93 above it, recovering the flat-then-declining shape. The two-piece model explains far more of the scatter than the line (R-squared 0.77 against 0.59), and once the extra breakpoint parameter is counted its AIC is lower by about 66.

How uncertain is the breakpoint?

The breakpoint is an estimate, not a landmark, and its standard error is not one of the numbers lm prints. A case bootstrap gives an honest interval: resample rows, refit the whole procedure, and collect the breakpoint each time.

est_psi <- function(x, y) {
  f <- function(p) sum(lm(y ~ x + pmax(x - p, 0))$residuals^2)
  g <- seq(3, 27, by = 0.25)
  p0 <- g[which.min(vapply(g, f, numeric(1)))]
  optimize(f, lower = p0 - 0.5, upper = p0 + 0.5)$minimum
}
set.seed(101)
boot_psi <- replicate(2000, {
  i <- sample.int(n, n, replace = TRUE)
  est_psi(dat$nitrogen[i], dat$richness[i])
})
psi_ci <- quantile(boot_psi, c(0.025, 0.975))
psi_se <- sd(boot_psi)

Species richness against nitrogen: a dashed straight line slopes gently down, while a solid broken-stick line is flat then turns sharply down at a breakpoint marked in gold with a shaded uncertainty band. — Figure 2: One straight line (dashed) against the broken-stick fit (solid). The gold line is the estimated breakpoint; the shaded band is its bootstrap 95% interval.

The bootstrap puts the breakpoint at roughly 8.9 to 14.3 kg N/ha/yr (standard error about 1.2), shown as the band in Figure 2. The interval covers the true value but spans several units: a study reporting a threshold of “13” without that band overstates how precisely the join is known.

The honest limit

Two cautions carry into the rest of this cluster. First, the breakpoint is a property of a model that assumes a break; the profile will always return a minimum, even for data with no genuine threshold, so a low RSS at some $\psi$ is not evidence that a threshold exists. Testing whether a break is real is a separate problem, because under the no-break hypothesis the breakpoint is not identified and the usual test on $\beta_2$ is invalid; that is the subject of the next post. Second, a sharp join is one shape among several. A smooth curve can fit a threshold-looking scatter just as well, and telling a genuine break from a gradual bend is a question the data alone rarely settle.

References

Toms JD, Lesperance ML 2003. Ecology 84(8):2034-2041 (10.1890/02-0472)
Muggeo VMR 2003. Statistics in Medicine 22(19):3055-3071 (10.1002/sim.1545)
Muggeo VMR 2008. R News 8(1):20-25 (ISSN 1609-3631)
Hudson DJ 1966. Journal of the American Statistical Association 61(316):1097-1129 (10.1080/01621459.1966.10482198)
Bolker BM 2008. Ecological Models and Data in R. Princeton University Press (ISBN 978-0-691-12522-0)

--- title: "Segmented regression for a breakpoint" description: "Fit a broken-stick model to an ecological threshold in base R: profile the residual sum of squares to place the breakpoint, and see why one straight line misleads." date: "2026-08-05 09:00" categories: [R, regression, ecology tutorial, ggplot2] image: thumbnail.png image-alt: "Scatter of species richness against nitrogen with a broken-stick fit that turns downwards at a breakpoint." --- ```{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)) } ``` Many ecological responses do not bend gently. A grassland loses forbs once nitrogen deposition passes a level the soil can buffer; a stream community collapses above a sediment load; edge effects fade at a fixed distance into a forest. In each case the relationship is roughly linear on either side of a **breakpoint**, and the breakpoint itself is the quantity of interest. A single straight line cannot represent that shape, and fitting one does not just lose detail: it produces a slope that describes neither side. This post fits a two-piece ("broken-stick") model in base R, locates the breakpoint by profiling the residual sum of squares, and reads off the slope on each side. The honest limit at the end is the theme of the whole cluster: a breakpoint is a feature of the model you chose, and its location carries real uncertainty. ## A threshold response Here is a simulated richness response to nitrogen deposition. Below a threshold the community is stable; above it, richness falls. ```{r data} set.seed(4198) n <- 120 nitrogen <- runif(n, 0, 30) # kg N/ha/yr psi_true <- 12 # the real breakpoint richness <- 25 + 0.10 * nitrogen - 0.90 * pmax(nitrogen - psi_true, 0) + rnorm(n, 0, 2.5) dat <- data.frame(nitrogen, richness) ``` The generating slope is a near-flat `r sprintf("%.2f", 0.10)` below the threshold and a declining `r sprintf("%.2f", 0.10 - 0.90)` above it, with the join at `r psi_true`. ## Why one line misleads The reflex is a single linear model. ```{r naive} fit_lin <- lm(richness ~ nitrogen, data = dat) lin_slope <- coef(fit_lin)[["nitrogen"]] lin_r2 <- summary(fit_lin)$r.squared ``` That line has a slope of `r round(lin_slope, 2)` and an R-squared of `r round(lin_r2, 2)`. The number looks like a tidy answer, but it is an average of two regimes that do not exist together anywhere in the data. Below the threshold the response is flat; above it, it drops almost twice as steeply as the single slope suggests. Reporting `r round(lin_slope, 2)` as "the effect of nitrogen" describes a mixture, not a mechanism. ## Fitting a broken stick A continuous two-piece model adds one hinge term. For a candidate breakpoint $\psi$, the term $(x - \psi)_+$ equals $x - \psi$ above $\psi$ and zero below it, so the model $$y = \beta_0 + \beta_1 x + \beta_2 (x - \psi)_+ + \varepsilon$$ is ordinary least squares *given* $\psi$: the slope is $\beta_1$ below the breakpoint and $\beta_1 + \beta_2$ above it. The one awkward parameter is $\psi$, which enters non-linearly. The direct approach is to profile it: fit the linear model at every candidate breakpoint and keep the value that minimises the residual sum of squares. ```{r segmented} rss_at <- function(p) { hinge <- pmax(dat$nitrogen - p, 0) sum(lm(richness ~ nitrogen + hinge, data = dat)$residuals^2) } grid <- seq(3, 27, by = 0.1) rss_grid <- vapply(grid, rss_at, numeric(1)) psi_hat <- optimize(rss_at, lower = grid[which.min(rss_grid)] - 0.2, upper = grid[which.min(rss_grid)] + 0.2)$minimum hinge <- pmax(dat$nitrogen - psi_hat, 0) fit_seg <- lm(richness ~ nitrogen + hinge, data = dat) cf <- coef(fit_seg) below_slope <- cf[["nitrogen"]] above_slope <- cf[["nitrogen"]] + cf[["hinge"]] seg_r2 <- summary(fit_seg)$r.squared ``` ```{r fig-profile} #| echo: false #| fig-cap: "Residual sum of squares against the candidate breakpoint. The fit is tightest at the estimate (gold); the dashed line marks the true value." #| fig-alt: "A U-shaped curve of residual sum of squares against candidate breakpoint, reaching a clear minimum near nitrogen 13, close to the dashed true value at 12." prof <- data.frame(grid, rss_grid) ggplot(prof, aes(grid, rss_grid)) + geom_vline(xintercept = psi_true, colour = te_sage, linewidth = 0.7, linetype = "22") + geom_line(colour = te_forest, linewidth = 1.0) + geom_point(data = data.frame(x = psi_hat, y = min(rss_grid)), aes(x, y), colour = te_gold, size = 2.6) + annotate("text", x = psi_true - 0.5, y = max(rss_grid) * 0.98, hjust = 1, label = "true", colour = te_label, size = 3.4) + annotate("text", x = psi_hat + 0.6, y = min(rss_grid) + 40, hjust = 0, label = "minimum", colour = te_label, size = 3.4) + labs(x = "Candidate breakpoint (kg N/ha/yr)", y = "Residual sum of squares", title = "The breakpoint is where the fit is tightest", subtitle = "Profile the RSS over every candidate, then take the minimum") + theme_te() ``` The profile in @fig-profile has a single clear valley. Its minimum puts the breakpoint at `r round(psi_hat, 1)` kg N/ha/yr, close to the true `r psi_true`. The fitted slopes are `r round(below_slope, 2)` below the breakpoint and `r round(above_slope, 2)` above it, recovering the flat-then-declining shape. The two-piece model explains far more of the scatter than the line (R-squared `r round(seg_r2, 2)` against `r round(lin_r2, 2)`), and once the extra breakpoint parameter is counted its AIC is lower by about `r round(AIC(fit_lin) - (AIC(fit_seg) + 2))`. ## How uncertain is the breakpoint? The breakpoint is an estimate, not a landmark, and its standard error is not one of the numbers `lm` prints. A case bootstrap gives an honest interval: resample rows, refit the whole procedure, and collect the breakpoint each time. ```{r bootstrap} est_psi <- function(x, y) { f <- function(p) sum(lm(y ~ x + pmax(x - p, 0))$residuals^2) g <- seq(3, 27, by = 0.25) p0 <- g[which.min(vapply(g, f, numeric(1)))] optimize(f, lower = p0 - 0.5, upper = p0 + 0.5)$minimum } set.seed(101) boot_psi <- replicate(2000, { i <- sample.int(n, n, replace = TRUE) est_psi(dat$nitrogen[i], dat$richness[i]) }) psi_ci <- quantile(boot_psi, c(0.025, 0.975)) psi_se <- sd(boot_psi) ``` ```{r fig-fits} #| echo: false #| fig-cap: "One straight line (dashed) against the broken-stick fit (solid). The gold line is the estimated breakpoint; the shaded band is its bootstrap 95% interval." #| fig-alt: "Species richness against nitrogen: a dashed straight line slopes gently down, while a solid broken-stick line is flat then turns sharply down at a breakpoint marked in gold with a shaded uncertainty band." gx <- seq(0, 30, length.out = 400) pred <- data.frame(gx, seg = cf[[1]] + cf[[2]] * gx + cf[[3]] * pmax(gx - psi_hat, 0), lin = coef(fit_lin)[[1]] + coef(fit_lin)[[2]] * gx) ggplot(dat, aes(nitrogen, richness)) + annotate("rect", xmin = psi_ci[[1]], xmax = psi_ci[[2]], ymin = -Inf, ymax = Inf, fill = te_sage, alpha = 0.20) + geom_point(colour = te_forest, alpha = 0.55, size = 1.7) + geom_line(data = pred, aes(gx, lin), colour = te_terra, linewidth = 1.0, linetype = "22") + geom_line(data = pred, aes(gx, seg), colour = te_forest, linewidth = 1.1) + geom_vline(xintercept = psi_hat, colour = te_gold, linewidth = 0.7) + annotate("text", x = 2, y = 8.5, hjust = 0, label = "single line (naive)", colour = te_terra, size = 3.4) + labs(x = "Nitrogen deposition (kg N/ha/yr)", y = "Species richness", title = "A single line hides the threshold", subtitle = "Broken-stick fit (solid) vs one straight line (dashed)") + theme_te() ``` The bootstrap puts the breakpoint at roughly `r round(psi_ci[[1]], 1)` to `r round(psi_ci[[2]], 1)` kg N/ha/yr (standard error about `r round(psi_se, 1)`), shown as the band in @fig-fits. The interval covers the true value but spans several units: a study reporting a threshold of "`r round(psi_hat, 1)`" without that band overstates how precisely the join is known. ## The honest limit Two cautions carry into the rest of this cluster. First, the breakpoint is a property of a model that *assumes* a break; the profile will always return a minimum, even for data with no genuine threshold, so a low RSS at some $\psi$ is not evidence that a threshold exists. Testing whether a break is real is a separate problem, because under the no-break hypothesis the breakpoint is not identified and the usual test on $\beta_2$ is invalid; that is the subject of the next post. Second, a sharp join is one shape among several. A smooth curve can fit a threshold-looking scatter just as well, and telling a genuine break from a gradual bend is a question the data alone rarely settle. ## References - Toms JD, Lesperance ML 2003. Ecology 84(8):2034-2041 (10.1890/02-0472) - Muggeo VMR 2003. Statistics in Medicine 22(19):3055-3071 (10.1002/sim.1545) - Muggeo VMR 2008. R News 8(1):20-25 (ISSN 1609-3631) - Hudson DJ 1966. Journal of the American Statistical Association 61(316):1097-1129 (10.1080/01621459.1966.10482198) - Bolker BM 2008. Ecological Models and Data in R. Princeton University Press (ISBN 978-0-691-12522-0) ## Related tutorials - [Detecting ecological thresholds](../detecting-ecological-thresholds/) - [Threshold regression with a GLM](../threshold-regression-glm/) - [Checking a threshold model](../checking-a-threshold-model/) - [Species-area relationships](../species-area-relationships/)