From GLM coefficients to predicted values in R

GLM

prediction

count data

logistic regression

ecology tutorial

Turn count and logistic GLM coefficients into predictions in R: build the confidence band on the link scale, back-transform it, and read the effect correctly.

Author

Tidy Ecology

Published

2026-05-08

Modified

2026-10-15

Step 26 of the Tidy Ecology course, Part five: model a response.

Corrected 15 October 2026: the post said that a symmetric count band rarely dips below zero because the variance is low where the mean is low; on the response scale the standard error is the predicted mean times the log-scale one, so the band dips below zero only where the log-scale standard error is large, as it can be for sparse counts, and the text now says so.

Corrected 7 October 2026: the caption of the conditioning figure and the text before it said that holding elevation at another value shifts the canopy curve up or down; on the count scale it multiplies the curve by one factor, so the curves fan out, and both now say so.

Corrected 30 September 2026: the post attributed to Fox the practice of holding a factor at a chosen level; his effect displays hold a factor at its proportional distribution in the data, and the text now says so.

Updated 27 September 2026: a new section, Predicting for new sites with a standardised covariate, measures what standardising new plots by their own mean does to a predicted count, and the text no longer says that elev_z has mean zero by construction or is held at its mean: zero is the centre of its scale, not the mean of the sample.

A negative binomial model of abundance hands you a table of coefficients on the log scale. The figure a reviewer actually wants is predicted abundance against the gradient, with a confidence band, on the count scale. Getting from one to the other has a single rule, and skipping it produces a band that is the wrong shape and can stray outside the values a count or a probability is allowed to take.

The coefficient lines we wrote earlier (count GLMs, presence/absence) stopped at reading the numbers. This one turns those numbers into the predicted curve you publish.

A model to predict from

We simulate overdispersed counts that rise with canopy openness and fall with elevation, then fit a negative binomial GLM. The covariate elev_z is standardised elevation: it is drawn on a scale with mean zero and standard deviation one.

library(MASS)
library(ggplot2)

set.seed(53)
n      <- 200
canopy <- runif(n, 0.05, 0.95)        # canopy openness, the focal predictor
elev_z <- rnorm(n)                    # standardised elevation, a covariate
mu     <- exp(1.0 + 2.2 * canopy - 0.5 * elev_z)
count  <- rnbinom(n, size = 2, mu = mu)
dat    <- data.frame(count, canopy, elev_z)

m <- glm.nb(count ~ canopy + elev_z, data = dat)
round(coef(m), 3)

(Intercept)      canopy      elev_z 
      1.024       2.058      -0.599

Colours and plot theme used below

ink<-"#16241d"; body<-"#2c3a31"; forest<-"#275139"; gold<-"#c9b458"
paper<-"#f5f4ee"; line<-"#dad9ca"; faint<-"#5d6b61"
theme_te <- function(){
  theme_minimal(base_size = 12) +
  theme(plot.background=element_rect(fill=paper,colour=NA),
        panel.background=element_rect(fill=paper,colour=NA),
        panel.grid.minor=element_blank(),
        panel.grid.major=element_line(colour=line,linewidth=.3),
        axis.title=element_text(colour=ink), axis.text=element_text(colour=faint),
        plot.title=element_text(colour=ink,size=12,face="bold"),
        strip.text=element_text(colour=ink,face="bold"),
        legend.position="bottom", legend.title=element_text(colour=ink),
        legend.text=element_text(colour=body))
}
crit <- qnorm(0.975)

The canopy coefficient is about 2.06. That is not a count and it is not a slope in individuals per unit; it lives on the log scale, where the model is linear. Exponentiating it gives 7.83 (95% CI 5.15 to 11.90), a multiplicative factor we will come back to. To get an actual predicted count we have to leave the log scale, and that step is where the band is usually mishandled.

Build the interval on the link scale, then transform

predict() will return the fitted values and their standard errors on either scale. The rule is to ask for the link scale, form the interval there with a symmetric fit +/- 1.96 * se, and only then apply the inverse link (exp for a log link). The interval inherits the link’s range, so a back-transformed count band can never go below zero.

grid <- data.frame(canopy = seq(0.05, 0.95, length.out = 100), elev_z = 0)
pl   <- predict(m, newdata = grid, type = "link", se.fit = TRUE)

grid$mu <- exp(pl$fit)                          # predicted mean count
grid$lo <- exp(pl$fit - crit * pl$se.fit)       # lower 95% limit
grid$hi <- exp(pl$fit + crit * pl$se.fit)       # upper 95% limit
head(round(grid, 2), 3)

  canopy elev_z   mu   lo   hi
1   0.05      0 3.09 2.47 3.86
2   0.06      0 3.15 2.52 3.92
3   0.07      0 3.21 2.58 3.98

Plotting the same prediction on both scales shows what the inverse link does to the interval. On the log scale the line is straight and the band is symmetric about it, narrowest in the middle of the canopy range. On the count scale the line bends upward and the band becomes lopsided: tighter below the curve, wider above it.

two <- rbind(
  data.frame(canopy=grid$canopy, fit=pl$fit, lo=pl$fit-crit*pl$se.fit,
             hi=pl$fit+crit*pl$se.fit, scale="link scale: log(mean count)"),
  data.frame(canopy=grid$canopy, fit=grid$mu, lo=grid$lo, hi=grid$hi,
             scale="response scale: mean count"))
two$scale <- factor(two$scale, levels=unique(two$scale))
ggplot(two, aes(canopy, fit)) +
  geom_ribbon(aes(ymin=lo, ymax=hi), fill=forest, alpha=.22) +
  geom_line(colour=forest, linewidth=.9) +
  facet_wrap(~scale, scales="free_y") +
  labs(x="canopy openness", y=NULL) +
  theme_te()

Two panels. Left, on the log scale, a straight line with a band symmetric about it, narrowest in the middle. Right, on the count scale, an upward-bending curve whose band is narrow below and wide above. — Figure 1: The fitted canopy effect on the link scale (left) and the response scale (right), with elevation held at zero, the centre of its scale. A symmetric interval on the left maps to a lopsided one on the right.

The asymmetry is not a quirk of the plot; it is the honest shape of the interval. A 95% interval for the log mean is symmetric, and the exponential of a symmetric interval is not symmetric. Reporting predicted +/- something on the count scale throws that shape away.

Why the symmetric shortcut misbehaves

The shortcut is to ask predict() for the response scale directly and draw fit +/- 1.96 * se. Because the standard error is now on the response scale and the band is forced symmetric, it is not held inside the range of the response. A logistic model of occurrence makes this visible: near the top of the gradient the symmetric upper limit pushes past a probability of 1, which is not a probability at all.

set.seed(206)
nb_   <- 180
moist <- runif(nb_, 0, 1)
occ   <- rbinom(nb_, 1, plogis(-2.5 + 6.0 * moist))
ml    <- glm(occ ~ moist, family = binomial)

gl  <- data.frame(moist = seq(0, 1, length.out = 100))
plk <- predict(ml, gl, type = "link",     se.fit = TRUE)
prr <- predict(ml, gl, type = "response", se.fit = TRUE)
gl$p   <- plogis(plk$fit)
gl$loc <- plogis(plk$fit - crit * plk$se.fit)   # correct: link CI, then plogis
gl$hic <- plogis(plk$fit + crit * plk$se.fit)
gl$lon <- prr$fit - crit * prr$se.fit           # naive: symmetric on response
gl$hin <- prr$fit + crit * prr$se.fit
c(naive_points_above_1 = sum(gl$hin > 1), max_naive_upper = round(max(gl$hin), 3))

naive_points_above_1      max_naive_upper 
               4.000                1.001

bandb <- rbind(
  data.frame(moist=gl$moist, lo=gl$loc, hi=gl$hic, kind="link-scale CI, then plogis() (correct)"),
  data.frame(moist=gl$moist, lo=gl$lon, hi=gl$hin, kind="response-scale SE, symmetric (naive)"))
ggplot() +
  geom_hline(yintercept=c(0,1), colour=faint, linewidth=.4, linetype="dotted") +
  geom_point(data=data.frame(moist,occ), aes(moist, occ), colour=faint, alpha=.3, size=1.1) +
  geom_ribbon(data=bandb, aes(moist, ymin=lo, ymax=hi, fill=kind), alpha=.20) +
  geom_line(data=bandb, aes(moist, hi, colour=kind), linewidth=.6) +
  geom_line(data=bandb, aes(moist, lo, colour=kind), linewidth=.6) +
  geom_line(data=gl, aes(moist, p), colour=ink, linewidth=.9) +
  scale_fill_manual(values=c(forest, gold)) +
  scale_colour_manual(values=c(forest, gold)) +
  coord_cartesian(ylim=c(-0.08, 1.08)) +
  labs(x="soil moisture", y="probability of occurrence", fill=NULL, colour=NULL) +
  theme_te()

A logistic curve from low to high probability. A gold band rises to meet the dotted line at probability one near the high end, while a forest band stays just below it. — Figure 2: Predicted probability of occurrence with two bands. The link-scale band (forest) stays inside zero and one by construction; the symmetric response-scale band (gold) passes one at the top of the gradient, reaching 1.001, which is too small a margin to see at this scale.

The overshoot here is small, only a little past 1, and where it appears is the informative part. Through the middle of the gradient the inverse logit is nearly straight across the width of the interval, so the symmetric band and the link-scale band all but coincide. At the two ends, where the curve bends hardest, they separate, and it is there that the symmetric upper limit runs past 1 with nothing to hold it back. A count band behaves the same way: on the response scale the standard error is the predicted mean times the log-scale standard error, so the symmetric band dips below zero only where that log-scale error is large, as it can be for sparse counts, but it is still the wrong shape everywhere. The link-scale route sidesteps all of this for free.

Marginal predictions: hold the other predictors somewhere

The grid above fixed elev_z at zero, the centre of its scale. That choice is part of the prediction, not a detail. To show the canopy effect you vary canopy and hold every other predictor at a representative value: the mean or median for a covariate, or, as here, the centre of its scale; for a factor, a chosen level, or, in the effect displays of Fox (2003), its proportional distribution in the data. Holding elevation at its tenth, fiftieth, and ninetieth percentile gives three versions of the same canopy curve, each a constant multiple of the others, so they are parallel on the log scale but fan out on the count scale.

elev_levels <- unname(quantile(dat$elev_z, c(0.1, 0.5, 0.9)))
elev_lab    <- c("low elevation","median","high elevation")
cond <- do.call(rbind, lapply(seq_along(elev_levels), function(k){
  g <- data.frame(canopy = seq(0.05, 0.95, length.out = 100), elev_z = elev_levels[k])
  p <- predict(m, newdata = g, type = "link", se.fit = TRUE)
  data.frame(canopy=g$canopy, elev_set=elev_lab[k], mu=exp(p$fit),
             lo=exp(p$fit-crit*p$se.fit), hi=exp(p$fit+crit*p$se.fit))
}))
cond$elev_set <- factor(cond$elev_set, levels=elev_lab)
ggplot(cond, aes(canopy, mu, colour=elev_set, fill=elev_set)) +
  geom_ribbon(aes(ymin=lo, ymax=hi), alpha=.16, colour=NA) +
  geom_line(linewidth=.9) +
  scale_colour_manual(values=c(forest, gold, "#b5534e")) +
  scale_fill_manual(values=c(forest, gold, "#b5534e")) +
  labs(x="canopy openness", y="predicted mean count",
       colour="elevation held at", fill="elevation held at") +
  theme_te()

Three upward curves of predicted count against canopy openness, one each for low, median, and high elevation, with the low-elevation curve highest. — Figure 3: The fitted canopy effect with elevation held at three values. On the count scale the covariate multiplies the whole curve by one factor, so the ratio between curves stays constant while the gap between them widens as canopy opens.

At a canopy openness of 0.5 the predicted count runs from about 16 at low elevation to about 3.4 at high elevation. None of those numbers is “the” prediction; each is a prediction conditional on where you parked the covariate, which is why the held value belongs in the figure caption.

Predicting for new sites with a standardised covariate

So far elev_z has been a number without units. Give it some: suppose it was made from elevations in metres with a centre of 900 m and a scale of 250 m, so that elev_z is exactly (elev_m - 900) / 250. Adding that column to dat changes nothing about m. Now a resurvey adds 25 plots at the upper end of the gradient, between 1150 and 1400 m, and we want their predicted counts. The slip is to standardise the new elevations with scale(), which uses the new plots’ own mean and standard deviation. The two safe routes are to reuse the constants the model was fitted with, or to put scale() inside the model formula.

dat$elev_m <- 900 + 250 * dat$elev_z                # elevation in metres behind elev_z
stopifnot(isTRUE(all.equal(dat$elev_z, (dat$elev_m - 900) / 250)))

set.seed(530)
new_sites <- data.frame(canopy = runif(25, 0.05, 0.95), elev_m = runif(25, 1150, 1400))
truth <- mean(exp(1.0 + 2.2 * new_sites$canopy - 0.5 * (new_sites$elev_m - 900) / 250))

z_own <- as.numeric(scale(new_sites$elev_m))        # the new plots' own mean and SD
z_fit <- (new_sites$elev_m - 900) / 250             # the constants the model was fitted with
p_own <- predict(m, data.frame(canopy = new_sites$canopy, elev_z = z_own), type = "response")
p_fit <- predict(m, data.frame(canopy = new_sites$canopy, elev_z = z_fit), type = "response")

m_sc  <- glm.nb(count ~ canopy + scale(elev_m), data = dat)   # scale() in the formula
p_in  <- predict(m_sc, newdata = new_sites, type = "response")
sc_call <- attr(terms(m_sc), "predvars")[[4]]
sc_call

scale(elev_m, center = 935.729120326921, scale = 241.71651311522)

stopifnot(isTRUE(all.equal(unname(p_in), unname(p_fit))), mean(p_own) > mean(p_fit), coef(m)[["elev_z"]] < 0,
          isTRUE(all.equal(as.numeric(logLik(m_sc)), as.numeric(logLik(m)))),
          isTRUE(all.equal(coef(m_sc)[[3]], coef(m)[["elev_z"]] * sc_call$scale / 250)),
          isTRUE(all.equal(c(sc_call$center, sc_call$scale), c(mean(dat$elev_m), sd(dat$elev_m)))))
round(c(true_mean = truth, own_scale = mean(p_own), fitted_scale = mean(p_fit),
        scale_in_formula = mean(p_in)), 2)

       true_mean        own_scale     fitted_scale scale_in_formula 
            5.16            11.74             4.19             4.19

The true mean count at these 25 plots, from the coefficients used to simulate the data, is 5.16. Standardised on themselves, the plots get a mean prediction of 11.74, 2.3 times that. On the fitted scale their elev_z runs from 1.04 to 1.98; on their own scale it runs from -1.59 to 2.12, so the model reads them as plots spread around zero, the centre of its scale, instead of at its upper end, and with a negative elevation coefficient that means higher counts. With the fitted constants the prediction is 4.19, so self-standardising by itself inflates the mean prediction 2.8 times. The gap between 4.19 and 5.16 is the estimation error of m itself, which any correct scaling inherits; it is not part of the trap.

The refit with scale(elev_m) in the formula gives the same 4.19; its 25 predictions equal the fitted-scale ones to within floating-point rounding. predict() is safe here because the model’s terms store the training centre and scale (the predvars call printed above; in R 4.3.3 the makepredictcall() method in stats fills them in for scale()) and apply them to the new rows. That centre is the sample mean, 935.7 m, not 900 m, and the predictions still agree, because a linear recoding of a covariate is the same model reparameterised (the slope rescaled, the intercept shifted): the two fits have the same log-likelihood.

up <- dat$elev_m + 200                                  # every original plot 200 m higher
pred_at <- function(z) predict(m, data.frame(canopy = dat$canopy, elev_z = z), type = "response")
own_now <- pred_at(as.numeric(scale(dat$elev_m)))       # standardised on themselves
own_up  <- pred_at(as.numeric(scale(up)))
fit_up  <- pred_at((up - 900) / 250)                    # the fitted constants
round(c(own_scale_change = max(abs(own_up - own_now)),
        fitted_scale_change = mean(fit_up) - mean(fitted(m))), 3)

   own_scale_change fitted_scale_change 
              0.000              -3.522

Self-standardising cannot see a change that every plot shares. Moving all 200 original plots 200 m uphill leaves their self-scaled predictions where they were, up to floating-point rounding: scale() subtracts the sample mean, so any constant added to every value drops out. That is a closed form, not a finding. On the fitted scale the same move lowers the mean predicted count by 3.52. Even without the move, the original plots standardised on themselves get a mean prediction of 10.19 against the fitted 9.25, because the sample of elev_z has mean 0.143 and standard deviation 0.967: standardised describes the scale, not this sample. Data leakage in ecological model validation shows the forgiving case, where training and test rows share one set of constants and whose constants they are makes no difference beyond floating-point noise; new plots given constants of their own are the case that does matter. Keep the centre and scale with the model, or let the formula keep them for you.

Reading the slope: multiplicative, not additive

On the response scale a log-link effect is multiplicative. Exponentiating the coefficients gives those factors, with Wald intervals:

round(exp(cbind(estimate = coef(m), confint.default(m))), 2)

            estimate 2.5 % 97.5 %
(Intercept)     2.79  2.19   3.55
canopy          7.83  5.15  11.90
elev_z          0.55  0.49   0.62

The exponentiated canopy coefficient, 7.83, is the factor by which the mean count multiplies per unit of canopy; per 0.1 of canopy the factor is exp(0.1 * 2.06), about 1.23, and it is the same factor everywhere along the gradient.

g2 <- data.frame(canopy = c(0.3, 0.4, 0.7, 0.8), elev_z = 0)
p2 <- predict(m, newdata = g2, type = "response")
c(ratio_0.3_to_0.4 = unname(round(p2[2]/p2[1], 3)),
  ratio_0.7_to_0.8 = unname(round(p2[4]/p2[3], 3)),
  jump_0.3_to_0.4  = unname(round(p2[2]-p2[1], 2)),
  jump_0.7_to_0.8  = unname(round(p2[4]-p2[3], 2)))

ratio_0.3_to_0.4 ratio_0.7_to_0.8  jump_0.3_to_0.4  jump_0.7_to_0.8 
           1.229            1.229            1.180            2.690

The ratio is constant at 1.23, but the absolute change is not: moving canopy from 0.3 to 0.4 adds about one individual, while the same step from 0.7 to 0.8 adds nearly three. That is why a single sentence like “canopy adds N individuals” has no fixed answer on the count scale, and why the predicted curve, not a lone coefficient, is the honest summary.

What the band does and does not cover

The interval here is for the mean count at a given canopy and elevation, not a prediction interval for a new plot. Individual plots scatter around the mean by the full negative binomial variance, which is much wider; if you need to bracket a future observation, simulate from the fitted model instead of widening this band. Two more limits worth stating: the held values are a choice, so report them, and the prediction has no anchor outside the range of the data, so do not read the curve past where canopy was actually sampled.

In practice you will reach for a helper rather than hand-building grids every time. The effects package (Fox 2003) and ggeffects (Ludecke 2018) compute adjusted predictions on the response scale, holding non-focal predictors at means or representative values, and return tidy frames ready for ggplot2; marginaleffects does the same and adds contrasts. The mechanics above are what those packages run underneath, which is the part worth seeing once before you let a function do it. The companion to this is checking that the model is worth predicting from at all: residual diagnostics for GLMs.

References

Fox J 2003 Journal of Statistical Software 8(15):1-27 (10.18637/jss.v008.i15)

Ludecke D 2018 Journal of Open Source Software 3(26):772 (10.21105/joss.00772)

McCullagh P, Nelder JA 1989 Generalized Linear Models, 2nd ed. Chapman and Hall. ISBN 0-412-31760-5

Nelder JA, Wedderburn RWM 1972 Journal of the Royal Statistical Society A 135(3):370-384 (10.2307/2344614)

Bolker BM 2008 Ecological Models and Data in R. Princeton University Press. ISBN 978-0-691-12522-0

--- title: "From GLM coefficients to predicted values in R" description: "Turn count and logistic GLM coefficients into predictions in R: build the confidence band on the link scale, back-transform it, and read the effect correctly." date: "2026-05-08 10:00" date-modified: "2026-10-15" categories: [GLM, prediction, count data, logistic regression, ecology tutorial] image: thumbnail.png image-alt: "A fitted curve shown twice, as a straight line with a symmetric band on the log scale and as a bending curve with a lopsided band on the count scale." --- *Corrected 15 October 2026: the post said that a symmetric count band rarely dips below zero because the variance is low where the mean is low; on the response scale the standard error is the predicted mean times the log-scale one, so the band dips below zero only where the log-scale standard error is large, as it can be for sparse counts, and the text now says so.* *Corrected 7 October 2026: the caption of the conditioning figure and the text before it said that holding elevation at another value shifts the canopy curve up or down; on the count scale it multiplies the curve by one factor, so the curves fan out, and both now say so.* *Corrected 30 September 2026: the post attributed to Fox the practice of holding a factor at a chosen level; his effect displays hold a factor at its proportional distribution in the data, and the text now says so.* *Updated 27 September 2026: a new section, Predicting for new sites with a standardised covariate, measures what standardising new plots by their own mean does to a predicted count, and the text no longer says that elev_z has mean zero by construction or is held at its mean: zero is the centre of its scale, not the mean of the sample.* A negative binomial model of abundance hands you a table of coefficients on the log scale. The figure a reviewer actually wants is predicted abundance against the gradient, with a confidence band, on the count scale. Getting from one to the other has a single rule, and skipping it produces a band that is the wrong shape and can stray outside the values a count or a probability is allowed to take. The coefficient lines we wrote earlier ([count GLMs](../glm-count-data-abundance/), [presence/absence](../logistic-regression-presence-absence/)) stopped at reading the numbers. This one turns those numbers into the predicted curve you publish. ## A model to predict from We simulate overdispersed counts that rise with canopy openness and fall with elevation, then fit a negative binomial GLM. The covariate `elev_z` is standardised elevation: it is drawn on a scale with mean zero and standard deviation one. ```{r} #| label: data #| warning: false #| message: false library(MASS) library(ggplot2) set.seed(53) n <- 200 canopy <- runif(n, 0.05, 0.95) # canopy openness, the focal predictor elev_z <- rnorm(n) # standardised elevation, a covariate mu <- exp(1.0 + 2.2 * canopy - 0.5 * elev_z) count <- rnbinom(n, size = 2, mu = mu) dat <- data.frame(count, canopy, elev_z) m <- glm.nb(count ~ canopy + elev_z, data = dat) round(coef(m), 3) ``` ```{r} #| label: setup #| code-fold: true #| code-summary: "Colours and plot theme used below" #| results: hide #| message: false #| warning: false ink<-"#16241d"; body<-"#2c3a31"; forest<-"#275139"; gold<-"#c9b458" paper<-"#f5f4ee"; line<-"#dad9ca"; faint<-"#5d6b61" theme_te <- function(){ theme_minimal(base_size = 12) + theme(plot.background=element_rect(fill=paper,colour=NA), panel.background=element_rect(fill=paper,colour=NA), panel.grid.minor=element_blank(), panel.grid.major=element_line(colour=line,linewidth=.3), axis.title=element_text(colour=ink), axis.text=element_text(colour=faint), plot.title=element_text(colour=ink,size=12,face="bold"), strip.text=element_text(colour=ink,face="bold"), legend.position="bottom", legend.title=element_text(colour=ink), legend.text=element_text(colour=body)) } crit <- qnorm(0.975) ``` The canopy coefficient is about `r sprintf("%.2f", coef(m)[["canopy"]])`. That is not a count and it is not a slope in individuals per unit; it lives on the log scale, where the model is linear. Exponentiating it gives `r sprintf("%.2f", exp(coef(m)[["canopy"]]))` (95% CI `r sprintf("%.2f", exp(confint.default(m)["canopy", 1]))` to `r sprintf("%.2f", exp(confint.default(m)["canopy", 2]))`), a multiplicative factor we will come back to. To get an actual predicted count we have to leave the log scale, and that step is where the band is usually mishandled. ## Build the interval on the link scale, then transform `predict()` will return the fitted values and their standard errors on either scale. The rule is to ask for the **link** scale, form the interval there with a symmetric `fit +/- 1.96 * se`, and only then apply the inverse link (`exp` for a log link). The interval inherits the link's range, so a back-transformed count band can never go below zero. ```{r} #| label: predict-grid grid <- data.frame(canopy = seq(0.05, 0.95, length.out = 100), elev_z = 0) pl <- predict(m, newdata = grid, type = "link", se.fit = TRUE) grid$mu <- exp(pl$fit) # predicted mean count grid$lo <- exp(pl$fit - crit * pl$se.fit) # lower 95% limit grid$hi <- exp(pl$fit + crit * pl$se.fit) # upper 95% limit head(round(grid, 2), 3) ``` Plotting the same prediction on both scales shows what the inverse link does to the interval. On the log scale the line is straight and the band is symmetric about it, narrowest in the middle of the canopy range. On the count scale the line bends upward and the band becomes lopsided: tighter below the curve, wider above it. ```{r} #| label: fig-link-response #| fig-cap: "The fitted canopy effect on the link scale (left) and the response scale (right), with elevation held at zero, the centre of its scale. A symmetric interval on the left maps to a lopsided one on the right." #| fig-alt: "Two panels. Left, on the log scale, a straight line with a band symmetric about it, narrowest in the middle. Right, on the count scale, an upward-bending curve whose band is narrow below and wide above." #| fig-width: 8.6 #| fig-height: 3.7 two <- rbind( data.frame(canopy=grid$canopy, fit=pl$fit, lo=pl$fit-crit*pl$se.fit, hi=pl$fit+crit*pl$se.fit, scale="link scale: log(mean count)"), data.frame(canopy=grid$canopy, fit=grid$mu, lo=grid$lo, hi=grid$hi, scale="response scale: mean count")) two$scale <- factor(two$scale, levels=unique(two$scale)) ggplot(two, aes(canopy, fit)) + geom_ribbon(aes(ymin=lo, ymax=hi), fill=forest, alpha=.22) + geom_line(colour=forest, linewidth=.9) + facet_wrap(~scale, scales="free_y") + labs(x="canopy openness", y=NULL) + theme_te() ``` The asymmetry is not a quirk of the plot; it is the honest shape of the interval. A 95% interval for the log mean is symmetric, and the exponential of a symmetric interval is not symmetric. Reporting `predicted +/- something` on the count scale throws that shape away. ## Why the symmetric shortcut misbehaves The shortcut is to ask `predict()` for the response scale directly and draw `fit +/- 1.96 * se`. Because the standard error is now on the response scale and the band is forced symmetric, it is not held inside the range of the response. A logistic model of occurrence makes this visible: near the top of the gradient the symmetric upper limit pushes past a probability of 1, which is not a probability at all. ```{r} #| label: logistic set.seed(206) nb_ <- 180 moist <- runif(nb_, 0, 1) occ <- rbinom(nb_, 1, plogis(-2.5 + 6.0 * moist)) ml <- glm(occ ~ moist, family = binomial) gl <- data.frame(moist = seq(0, 1, length.out = 100)) plk <- predict(ml, gl, type = "link", se.fit = TRUE) prr <- predict(ml, gl, type = "response", se.fit = TRUE) gl$p <- plogis(plk$fit) gl$loc <- plogis(plk$fit - crit * plk$se.fit) # correct: link CI, then plogis gl$hic <- plogis(plk$fit + crit * plk$se.fit) gl$lon <- prr$fit - crit * prr$se.fit # naive: symmetric on response gl$hin <- prr$fit + crit * prr$se.fit c(naive_points_above_1 = sum(gl$hin > 1), max_naive_upper = round(max(gl$hin), 3)) ``` ```{r} #| label: fig-boundary #| fig-cap: "Predicted probability of occurrence with two bands. The link-scale band (forest) stays inside zero and one by construction; the symmetric response-scale band (gold) passes one at the top of the gradient, reaching 1.001, which is too small a margin to see at this scale." #| fig-alt: "A logistic curve from low to high probability. A gold band rises to meet the dotted line at probability one near the high end, while a forest band stays just below it." #| fig-width: 8.0 #| fig-height: 4.4 bandb <- rbind( data.frame(moist=gl$moist, lo=gl$loc, hi=gl$hic, kind="link-scale CI, then plogis() (correct)"), data.frame(moist=gl$moist, lo=gl$lon, hi=gl$hin, kind="response-scale SE, symmetric (naive)")) ggplot() + geom_hline(yintercept=c(0,1), colour=faint, linewidth=.4, linetype="dotted") + geom_point(data=data.frame(moist,occ), aes(moist, occ), colour=faint, alpha=.3, size=1.1) + geom_ribbon(data=bandb, aes(moist, ymin=lo, ymax=hi, fill=kind), alpha=.20) + geom_line(data=bandb, aes(moist, hi, colour=kind), linewidth=.6) + geom_line(data=bandb, aes(moist, lo, colour=kind), linewidth=.6) + geom_line(data=gl, aes(moist, p), colour=ink, linewidth=.9) + scale_fill_manual(values=c(forest, gold)) + scale_colour_manual(values=c(forest, gold)) + coord_cartesian(ylim=c(-0.08, 1.08)) + labs(x="soil moisture", y="probability of occurrence", fill=NULL, colour=NULL) + theme_te() ``` The overshoot here is small, only a little past 1, and where it appears is the informative part. Through the middle of the gradient the inverse logit is nearly straight across the width of the interval, so the symmetric band and the link-scale band all but coincide. At the two ends, where the curve bends hardest, they separate, and it is there that the symmetric upper limit runs past 1 with nothing to hold it back. A count band behaves the same way: on the response scale the standard error is the predicted mean times the log-scale standard error, so the symmetric band dips below zero only where that log-scale error is large, as it can be for sparse counts, but it is still the wrong shape everywhere. The link-scale route sidesteps all of this for free. ## Marginal predictions: hold the other predictors somewhere The grid above fixed `elev_z` at zero, the centre of its scale. That choice is part of the prediction, not a detail. To show the canopy effect you vary canopy and hold every other predictor at a representative value: the mean or median for a covariate, or, as here, the centre of its scale; for a factor, a chosen level, or, in the effect displays of Fox (2003), its proportional distribution in the data. Holding elevation at its tenth, fiftieth, and ninetieth percentile gives three versions of the same canopy curve, each a constant multiple of the others, so they are parallel on the log scale but fan out on the count scale. ```{r} #| label: fig-conditioning #| fig-cap: "The fitted canopy effect with elevation held at three values. On the count scale the covariate multiplies the whole curve by one factor, so the ratio between curves stays constant while the gap between them widens as canopy opens." #| fig-alt: "Three upward curves of predicted count against canopy openness, one each for low, median, and high elevation, with the low-elevation curve highest." #| fig-width: 8.0 #| fig-height: 4.4 elev_levels <- unname(quantile(dat$elev_z, c(0.1, 0.5, 0.9))) elev_lab <- c("low elevation","median","high elevation") cond <- do.call(rbind, lapply(seq_along(elev_levels), function(k){ g <- data.frame(canopy = seq(0.05, 0.95, length.out = 100), elev_z = elev_levels[k]) p <- predict(m, newdata = g, type = "link", se.fit = TRUE) data.frame(canopy=g$canopy, elev_set=elev_lab[k], mu=exp(p$fit), lo=exp(p$fit-crit*p$se.fit), hi=exp(p$fit+crit*p$se.fit)) })) cond$elev_set <- factor(cond$elev_set, levels=elev_lab) ggplot(cond, aes(canopy, mu, colour=elev_set, fill=elev_set)) + geom_ribbon(aes(ymin=lo, ymax=hi), alpha=.16, colour=NA) + geom_line(linewidth=.9) + scale_colour_manual(values=c(forest, gold, "#b5534e")) + scale_fill_manual(values=c(forest, gold, "#b5534e")) + labs(x="canopy openness", y="predicted mean count", colour="elevation held at", fill="elevation held at") + theme_te() ``` At a canopy openness of 0.5 the predicted count runs from about `r sprintf("%.0f", predict(m, data.frame(canopy = 0.5, elev_z = elev_levels[1]), type = "response"))` at low elevation to about `r sprintf("%.1f", predict(m, data.frame(canopy = 0.5, elev_z = elev_levels[3]), type = "response"))` at high elevation. None of those numbers is "the" prediction; each is a prediction conditional on where you parked the covariate, which is why the held value belongs in the figure caption. ## Predicting for new sites with a standardised covariate So far `elev_z` has been a number without units. Give it some: suppose it was made from elevations in metres with a centre of 900 m and a scale of 250 m, so that `elev_z` is exactly `(elev_m - 900) / 250`. Adding that column to `dat` changes nothing about `m`. Now a resurvey adds 25 plots at the upper end of the gradient, between 1150 and 1400 m, and we want their predicted counts. The slip is to standardise the new elevations with `scale()`, which uses the new plots' own mean and standard deviation. The two safe routes are to reuse the constants the model was fitted with, or to put `scale()` inside the model formula. ```{r} #| label: new-sites dat$elev_m <- 900 + 250 * dat$elev_z # elevation in metres behind elev_z stopifnot(isTRUE(all.equal(dat$elev_z, (dat$elev_m - 900) / 250))) set.seed(530) new_sites <- data.frame(canopy = runif(25, 0.05, 0.95), elev_m = runif(25, 1150, 1400)) truth <- mean(exp(1.0 + 2.2 * new_sites$canopy - 0.5 * (new_sites$elev_m - 900) / 250)) z_own <- as.numeric(scale(new_sites$elev_m)) # the new plots' own mean and SD z_fit <- (new_sites$elev_m - 900) / 250 # the constants the model was fitted with p_own <- predict(m, data.frame(canopy = new_sites$canopy, elev_z = z_own), type = "response") p_fit <- predict(m, data.frame(canopy = new_sites$canopy, elev_z = z_fit), type = "response") m_sc <- glm.nb(count ~ canopy + scale(elev_m), data = dat) # scale() in the formula p_in <- predict(m_sc, newdata = new_sites, type = "response") sc_call <- attr(terms(m_sc), "predvars")[[4]] sc_call stopifnot(isTRUE(all.equal(unname(p_in), unname(p_fit))), mean(p_own) > mean(p_fit), coef(m)[["elev_z"]] < 0, isTRUE(all.equal(as.numeric(logLik(m_sc)), as.numeric(logLik(m)))), isTRUE(all.equal(coef(m_sc)[[3]], coef(m)[["elev_z"]] * sc_call$scale / 250)), isTRUE(all.equal(c(sc_call$center, sc_call$scale), c(mean(dat$elev_m), sd(dat$elev_m))))) round(c(true_mean = truth, own_scale = mean(p_own), fitted_scale = mean(p_fit), scale_in_formula = mean(p_in)), 2) ``` The true mean count at these 25 plots, from the coefficients used to simulate the data, is `r sprintf("%.2f", truth)`. Standardised on themselves, the plots get a mean prediction of `r sprintf("%.2f", mean(p_own))`, `r sprintf("%.1f", mean(p_own) / truth)` times that. On the fitted scale their `elev_z` runs from `r sprintf("%.2f", min(z_fit))` to `r sprintf("%.2f", max(z_fit))`; on their own scale it runs from `r sprintf("%.2f", min(z_own))` to `r sprintf("%.2f", max(z_own))`, so the model reads them as plots spread around zero, the centre of its scale, instead of at its upper end, and with a negative elevation coefficient that means higher counts. With the fitted constants the prediction is `r sprintf("%.2f", mean(p_fit))`, so self-standardising by itself inflates the mean prediction `r sprintf("%.1f", mean(p_own) / mean(p_fit))` times. The gap between `r sprintf("%.2f", mean(p_fit))` and `r sprintf("%.2f", truth)` is the estimation error of `m` itself, which any correct scaling inherits; it is not part of the trap. The refit with `scale(elev_m)` in the formula gives the same `r sprintf("%.2f", mean(p_in))`; its 25 predictions equal the fitted-scale ones to within floating-point rounding. `predict()` is safe here because the model's terms store the training centre and scale (the `predvars` call printed above; in R 4.3.3 the `makepredictcall()` method in stats fills them in for `scale()`) and apply them to the new rows. That centre is the sample mean, `r sprintf("%.1f", sc_call$center)` m, not 900 m, and the predictions still agree, because a linear recoding of a covariate is the same model reparameterised (the slope rescaled, the intercept shifted): the two fits have the same log-likelihood. ```{r} #| label: fp-guard-1 #| include: false stopifnot(max(abs(p_in - p_fit)) / max(abs(p_fit)) < 1e-10) ``` ```{r} #| label: uniform-shift up <- dat$elev_m + 200 # every original plot 200 m higher pred_at <- function(z) predict(m, data.frame(canopy = dat$canopy, elev_z = z), type = "response") own_now <- pred_at(as.numeric(scale(dat$elev_m))) # standardised on themselves own_up <- pred_at(as.numeric(scale(up))) fit_up <- pred_at((up - 900) / 250) # the fitted constants round(c(own_scale_change = max(abs(own_up - own_now)), fitted_scale_change = mean(fit_up) - mean(fitted(m))), 3) ``` Self-standardising cannot see a change that every plot shares. Moving all 200 original plots 200 m uphill leaves their self-scaled predictions where they were, up to floating-point rounding: `scale()` subtracts the sample mean, so any constant added to every value drops out. That is a closed form, not a finding. On the fitted scale the same move lowers the mean predicted count by `r sprintf("%.2f", mean(fitted(m)) - mean(fit_up))`. Even without the move, the original plots standardised on themselves get a mean prediction of `r sprintf("%.2f", mean(own_now))` against the fitted `r sprintf("%.2f", mean(fitted(m)))`, because the sample of `elev_z` has mean `r sprintf("%.3f", mean(dat$elev_z))` and standard deviation `r sprintf("%.3f", sd(dat$elev_z))`: standardised describes the scale, not this sample. [Data leakage in ecological model validation](../data-leakage-in-model-validation/) shows the forgiving case, where training and test rows share one set of constants and whose constants they are makes no difference beyond floating-point noise; new plots given constants of their own are the case that does matter. Keep the centre and scale with the model, or let the formula keep them for you. ```{r} #| label: fp-guard-2 #| include: false stopifnot(max(abs(own_up - own_now)) / max(abs(own_now)) < 1e-10) ``` ## Reading the slope: multiplicative, not additive On the response scale a log-link effect is multiplicative. Exponentiating the coefficients gives those factors, with Wald intervals: ```{r} #| label: exp-coef round(exp(cbind(estimate = coef(m), confint.default(m))), 2) ``` The exponentiated canopy coefficient, `r sprintf("%.2f", exp(coef(m)[["canopy"]]))`, is the factor by which the mean count multiplies per unit of canopy; per 0.1 of canopy the factor is `exp(0.1 * 2.06)`, about `r sprintf("%.2f", exp(0.1 * coef(m)[["canopy"]]))`, and it is the same factor everywhere along the gradient. ```{r} #| label: multiplicative g2 <- data.frame(canopy = c(0.3, 0.4, 0.7, 0.8), elev_z = 0) p2 <- predict(m, newdata = g2, type = "response") c(ratio_0.3_to_0.4 = unname(round(p2[2]/p2[1], 3)), ratio_0.7_to_0.8 = unname(round(p2[4]/p2[3], 3)), jump_0.3_to_0.4 = unname(round(p2[2]-p2[1], 2)), jump_0.7_to_0.8 = unname(round(p2[4]-p2[3], 2))) ``` The ratio is constant at `r sprintf("%.2f", p2[2] / p2[1])`, but the absolute change is not: moving canopy from 0.3 to 0.4 adds about one individual, while the same step from 0.7 to 0.8 adds nearly three. That is why a single sentence like "canopy adds N individuals" has no fixed answer on the count scale, and why the predicted curve, not a lone coefficient, is the honest summary. ## What the band does and does not cover The interval here is for the mean count at a given canopy and elevation, not a prediction interval for a new plot. Individual plots scatter around the mean by the full negative binomial variance, which is much wider; if you need to bracket a future observation, simulate from the fitted model instead of widening this band. Two more limits worth stating: the held values are a choice, so report them, and the prediction has no anchor outside the range of the data, so do not read the curve past where canopy was actually sampled. In practice you will reach for a helper rather than hand-building grids every time. The `effects` package (Fox 2003) and `ggeffects` (Ludecke 2018) compute adjusted predictions on the response scale, holding non-focal predictors at means or representative values, and return tidy frames ready for `ggplot2`; `marginaleffects` does the same and adds contrasts. The mechanics above are what those packages run underneath, which is the part worth seeing once before you let a function do it. The companion to this is checking that the model is worth predicting from at all: [residual diagnostics for GLMs](../glm-residual-diagnostics/). ## References Fox J 2003 Journal of Statistical Software 8(15):1-27 (10.18637/jss.v008.i15) Ludecke D 2018 Journal of Open Source Software 3(26):772 (10.21105/joss.00772) McCullagh P, Nelder JA 1989 Generalized Linear Models, 2nd ed. Chapman and Hall. ISBN 0-412-31760-5 Nelder JA, Wedderburn RWM 1972 Journal of the Royal Statistical Society A 135(3):370-384 (10.2307/2344614) Bolker BM 2008 Ecological Models and Data in R. Princeton University Press. ISBN 978-0-691-12522-0 ## Related tutorials - [GLM residual diagnostics](../glm-residual-diagnostics/) - [Interaction terms in GLMs](../interaction-terms-glm/) - [Contrasts after a GLM](../contrasts-after-glm/) - [Effect plots that show the data](../effect-plots-that-show-the-data/)