---
title: "Checking a thermal performance analysis"
description: "Three checks for a thermal performance curve in R: the extrapolated critical maximum, curve-family dependence, and measurement error in body temperature."
date: "2026-07-29 12:00"
categories: [ecology tutorial, R, thermal ecology, model diagnostics, physiology]
image: thumbnail.png
image-alt: "Thermal performance curves fitted to the same data diverging sharply on the descending limb where measurements are sparse."
---
```{r setup}
#| include: false
knitr::opts_chunk$set(echo = TRUE, message = FALSE, warning = FALSE,
dev = "png", dpi = 120, fig.path = "figures/",
fig.width = 7, fig.height = 4.4)
library(ggplot2)
theme_te <- theme_minimal(base_size = 12) +
theme(panel.background = element_rect(fill = "white", colour = NA),
plot.background = element_rect(fill = "white", colour = NA),
panel.grid.minor = element_blank(),
plot.title = element_text(face = "bold", size = 12))
theme_set(theme_te)
col_full <- "#3a6b52"; col_cut <- "#b5651d"; col_err <- "#b5651d"; col_true <- "#3a6b52"; col_ref <- "#555555"
briere <- function(T, a, T0, Tm) { p <- a*T*(T-T0)*sqrt(pmax(Tm-T, 0)); p[T < T0 | T > Tm] <- 0; p }
a0 <- 8e-4; T0 <- 8; Tm <- 42
Tg <- seq(0, 45, 0.005); Topt <- Tg[which.max(briere(Tg, a0, T0, Tm))]; Pmax0 <- max(briere(Tg, a0, T0, Tm))
set.seed(394)
Td <- rep(seq(10, 40, 2.5), each = 4)
y <- pmax(briere(Td, a0, T0, Tm) + rnorm(length(Td), 0, 0.06), 0)
fitTm <- function(T, yy, st = 44) {
f <- nls(yy ~ aa*T*(T-t0)*sqrt(pmax(tm-T, 0)),
start = list(aa = 6e-4, t0 = 6, tm = st), control = nls.control(warnOnly = TRUE, maxiter = 200))
s <- summary(f)$coefficients; c(est.tm = s["tm","Estimate"], se.tm = s["tm","Std. Error"], coef(f))
}
full <- fitTm(Td, y)
r38 <- fitTm(Td[Td <= 37.5], y[Td <= 37.5])
r35 <- fitTm(Td[Td <= 35], y[Td <= 35])
n_above <- sum(Td > Topt); hot <- max(Td)
## curve-family disagreement on the upper limit
cb <- coef(nls(y ~ aa*Td*(Td-t0)*sqrt(pmax(tm-Td, 0)), start = list(aa = 6e-4, t0 = 6, tm = 44)))
cg <- coef(nls(y ~ A*exp(-((Td-mo)^2)/(2*s^2)), start = list(A = 2, mo = 34, s = 8)))
ctmax_bri <- cb["tm"]; ctmax_gau <- cg["mo"] + cg["s"]*sqrt(-2*log(0.05))
## measurement error in body temperature
set.seed(394)
Terr <- Td + rnorm(length(Td), 0, 2.5)
ce <- coef(nls(y ~ aa*Terr*(Terr-t0)*sqrt(pmax(tm-Terr, 0)),
start = list(aa = 6e-4, t0 = 6, tm = 46), control = nls.control(warnOnly = TRUE, maxiter = 200)))
TgE <- seq(ce["t0"], ce["tm"], 0.01); PgE <- ce["aa"]*TgE*(TgE-ce["t0"])*sqrt(pmax(ce["tm"]-TgE, 0))
PmaxE <- max(PgE); pmax_drop <- 100 * (PmaxE/Pmax0 - 1)
```
A [thermal performance curve](../thermal-performance-curves/) delivers three headline numbers, the optimum, the breadth, and the critical thermal maximum, and they carry very different weights of evidence. The optimum and breadth sit inside the data. CTmax usually does not. This post is three checks that a performance-curve result should survive before those numbers are used to rank species or predict range limits.
## Check 1: how much data supports CTmax?
CTmax is where the fitted curve reaches zero, and that point is almost always past the hottest temperature anyone measured. Refit the curve after dropping the warmest data and watch the estimate move.
```{r ctmax}
fitTm <- function(T, yy) {
f <- nls(yy ~ aa*T*(T-t0)*sqrt(pmax(tm-T, 0)),
start = list(aa = 6e-4, t0 = 6, tm = 44), control = nls.control(warnOnly = TRUE))
summary(f)$coefficients["tm", c("Estimate","Std. Error")]
}
rbind(full = round(fitTm(Td, y), 2),
drop_37.5 = round(fitTm(Td[Td <= 37.5], y[Td <= 37.5]), 2),
drop_35 = round(fitTm(Td[Td <= 35], y[Td <= 35]), 2))
```
The hottest measurement is `r sprintf("%.0f", hot)` degrees C, yet CTmax comes out near `r sprintf("%.1f", full["est.tm"])`: it is already an extrapolation `r sprintf("%.0f", full["est.tm"]-hot)` degrees beyond the data, and only `r n_above` of the `r length(Td)` points sit above the optimum at all. Drop the descending-limb points and the estimate drifts up to `r sprintf("%.1f", r35["est.tm"])`, but the real tell is the standard error, which grows from `r sprintf("%.1f", full["se.tm"])` to `r sprintf("%.1f", r35["se.tm"])` degrees. If sampling stops near the optimum, CTmax is not estimated, it is guessed. The check: report CTmax with its interval, and refuse to compare CTmax across species unless each is anchored by measurements on the descending limb.
```{r fig-ctmax}
#| fig-cap: "The upper limit is pinned by a handful of hot points. The full-data fit (green) reaches zero near 42 degrees; refit without the two hottest temperatures (ochre) and both the estimate and its uncertainty grow."
#| fig-alt: "Performance data with two fitted curves; the green full-data curve falls to zero near 42 degrees, the ochre reduced-data curve extends further right, showing the upper limit is poorly constrained."
cf <- function(cc) { g <- seq(cc["t0"], cc["tm"], 0.05); data.frame(T = g, P = cc["aa"]*g*(g-cc["t0"])*sqrt(pmax(cc["tm"]-g, 0))) }
cA <- cf(full[3:5]); cB <- cf(r35[3:5])
ggplot(data.frame(Td, y), aes(Td, y)) +
geom_point(colour = col_ref, size = 2, alpha = 0.7) +
geom_line(data = cA, aes(T, P), colour = col_full, linewidth = 1) +
geom_line(data = cB, aes(T, P), colour = col_cut, linewidth = 1, linetype = "22") +
geom_vline(xintercept = 35, colour = col_ref, linetype = "dotted") +
annotate("text", x = 35.3, y = 0.15, label = "data dropped above here", hjust = 0, colour = col_ref, size = 3.2) +
labs(x = "Temperature (C)", y = "Performance", title = "CTmax is an extrapolation")
```
## Check 2: does the answer depend on the curve family?
Because CTmax lives in the extrapolated region, the model chosen to get there matters more than anywhere else on the curve. Fit the asymmetric Briere and a symmetric Gaussian to the same data and compare where each says performance vanishes.
```{r family}
cb <- coef(nls(y ~ aa*Td*(Td-t0)*sqrt(pmax(tm-Td, 0)), start = list(aa = 6e-4, t0 = 6, tm = 44)))
cg <- coef(nls(y ~ A*exp(-((Td-mo)^2)/(2*s^2)), start = list(A = 2, mo = 34, s = 8)))
c(briere_CTmax = unname(cb["tm"]), gaussian_upper = unname(cg["mo"] + cg["s"]*sqrt(-2*log(0.05))))
```
The Briere model puts the upper limit at `r sprintf("%.1f", ctmax_bri)` degrees; the Gaussian, reaching five per cent of its peak, at `r sprintf("%.1f", ctmax_gau)`. The two agree on the interpolated peak and disagree by `r sprintf("%.0f", ctmax_gau - ctmax_bri)` degrees on the extrapolated limit, from identical data. When a number depends this strongly on a modelling choice, fit more than one family and report the spread, rather than trusting the decimal a single model returns.
## Check 3: is the temperature measured with error?
Performance is regressed on body temperature, but body temperature is hard to measure and is often replaced by air temperature. Error in the predictor flattens a nonlinear curve, the same regression dilution that attenuates a straight-line slope, and it biases every summary.
```{r measerror}
set.seed(394)
Terr <- Td + rnorm(length(Td), 0, 2.5) # body temperature measured with error
ce <- coef(nls(y ~ aa*Terr*(Terr-t0)*sqrt(pmax(tm-Terr, 0)),
start = list(aa = 6e-4, t0 = 6, tm = 46), control = nls.control(warnOnly = TRUE)))
round(c(true_CTmax = 42, fitted_CTmax = unname(ce["tm"])), 1)
```
A temperature standard deviation of `r sprintf("%.1f", 2.5)` degrees drops the apparent peak performance by `r sprintf("%.1f", -pmax_drop)` per cent, widens the curve, and pushes the apparent CTmax up to `r sprintf("%.1f", ce["tm"])` degrees. A curve fitted to air temperature will look flatter and more heat-tolerant than the animal really is. The check: use body temperature, not a proxy, and where a proxy is unavoidable, treat the estimated breadth as an upper bound and the optimum as approximate ([measurement error and regression dilution](../measurement-error-regression-dilution/) has the general version).
```{r fig-error}
#| fig-cap: "Measurement error in temperature flattens the curve. The true curve (green) fitted to error-laden temperatures (ochre) loses height, widens, and reports a higher critical maximum."
#| fig-alt: "Two thermal performance curves; the true green curve is taller and narrower, the ochre curve fitted with temperature error is flatter, wider and extends to a higher upper limit."
ggplot() +
geom_line(data = data.frame(T = Tg, P = briere(Tg, a0, T0, Tm)), aes(T, P), colour = col_true, linewidth = 1) +
geom_line(data = data.frame(T = TgE, P = PgE), aes(T, P), colour = col_err, linewidth = 1, linetype = "22") +
annotate("text", x = 14, y = Pmax0*0.96, label = "true", colour = col_true, size = 3.6) +
annotate("text", x = 14, y = PmaxE*0.9, label = "fitted with temperature error", colour = col_err, hjust = 0, size = 3.4) +
labs(x = "Temperature (C)", y = "Performance", title = "Temperature error flattens the curve")
```
## The deeper caveat
These three checks bound the numbers a curve returns. They do not fix what the curve is. A performance curve is measured by holding an animal at a constant temperature and reading a rate, and the two earlier posts showed why that is not the field: performance under [fluctuating temperature](../jensens-inequality-and-thermal-variability/) is not performance at the mean, and the curve itself shifts with acclimation and can be edited by behaviour (Sinclair et al. 2016 Ecology Letters 19(11):1372-1385). Run the three checks so the numbers are honest, then remember that even an honest constant-temperature curve is a laboratory object, and the animal lives outside.
## Related tutorials
- [Thermal performance curves in R](../thermal-performance-curves/)
- [Jensen's inequality and thermal variability](../jensens-inequality-and-thermal-variability/)
- [Thermal safety margins and warming](../thermal-safety-margins-and-warming/)
- [Measurement error and regression dilution](../measurement-error-regression-dilution/)
## References
- Kingsolver JG 2009. The American Naturalist 174(6):755-768 (10.1086/648310).
- Schoolfield RM et al. 1981. Journal of Theoretical Biology 88(4):719-731 (10.1016/0022-5193(81)90246-0).
- Sinclair BJ et al. 2016. Ecology Letters 19(11):1372-1385 (10.1111/ele.12686).