Model-averaged coefficients: a caution

multimodel inference

model averaging

collinearity

ggplot2

ecology tutorial

Model-averaged regression coefficients mislead under collinearity in ecology (Cade 2015). Average predictions, not coefficients: a worked R demonstration.

Author

Tidy Ecology

Published

2026-07-23

Averaging predictions across a candidate set is safe. Averaging coefficients across the same set is where multi-model inference most often goes wrong in ecology, and it goes wrong precisely when it looks most useful: when predictors are correlated and you want a single number for each one’s effect.

The problem is not a coding slip. It is that a regression coefficient means different things in different models. In a model that also contains a correlated predictor, a coefficient is a partial effect, holding that other predictor fixed. In a model that omits it, the same coefficient absorbs the shared variation and becomes something larger. Averaging those two numbers, weighted by Akaike weights, blends quantities that are not on the same footing. Cade set this out in detail in 2015; this post reproduces the core of the argument in R and shows why the prediction stays trustworthy even as the averaged coefficient does not.

The Akaike weights post is the prerequisite, and the collinearity and variance inflation post covers the mechanism this one builds on.

Two correlated gradients

Eighty forest plots, log above-ground biomass as the response. Two environmental predictors carry real effects: soil moisture and soil organic carbon, which in real landscapes tend to rise together. A third predictor, slope aspect, has no effect. We control the correlation between moisture and organic carbon directly so we can watch what it does.

gen_data <- function(rho, N = 80, seed = 505) {
  if (!is.na(seed)) set.seed(seed)          # seed = NA gives a fresh draw each call
  z1 <- rnorm(N); z2 <- rnorm(N)
  moisture <- z1
  carbon   <- rho * z1 + sqrt(1 - rho^2) * z2   # correlation with moisture = rho
  aspect   <- rnorm(N)
  biomass  <- 2 + 1.0 * moisture + 1.4 * carbon + 0 * aspect + rnorm(N, 0, 2.0)
  data.frame(biomass, moisture, carbon, aspect)
}

The true coefficient of moisture is 1.0 throughout, regardless of how strongly it correlates with carbon. Any change we see in the estimate is an artefact of the modelling, not the ecology.

One coefficient, three different answers

Take a single dataset with strong correlation, and estimate the moisture coefficient two ways: from the model that also includes carbon (a partial effect), and from the model with moisture alone (a marginal effect).

d09 <- gen_data(rho = 0.9)
b_partial  <- coef(lm(biomass ~ moisture + carbon, data = d09))["moisture"]
b_marginal <- coef(lm(biomass ~ moisture, data = d09))["moisture"]
c(partial = as.numeric(b_partial), marginal = as.numeric(b_marginal))

  partial  marginal 
0.6885105 2.5707228

The two estimates disagree sharply. The partial effect is 0.69; the marginal is 2.57, well over three times larger, because with carbon left out the moisture term stands in for both gradients at once. On this one dataset the partial estimate sits below the true value of 1.0, a reminder that any single fit is a noisy draw, especially with strong collinearity inflating its standard error; the systematic behaviour comes from averaging over many datasets, which we do shortly. What holds on every draw is the gap: neither estimate is wrong for what it is, they simply answer different questions, and a model-averaged coefficient mixes them.

To see the averaged coefficient itself, we need the full candidate set. With three predictors there are eight models; the averaged moisture coefficient adds up each model’s moisture estimate weighted by its Akaike weight, counting models without moisture as contributing zero.

aicc <- function(m) { k <- length(coef(m)) + 1; N <- nobs(m); AIC(m) + 2 * k * (k + 1) / (N - k - 1) }
vars <- c("moisture", "carbon", "aspect")
subsets <- unlist(lapply(0:3, function(k) combn(vars, k, simplify = FALSE)), recursive = FALSE)
forms   <- vapply(subsets, function(s)
  if (length(s) == 0) "biomass ~ 1" else paste("biomass ~", paste(s, collapse = " + ")), character(1))

avg_moisture <- function(d) {
  fits <- lapply(forms, function(f) lm(as.formula(f), data = d))
  w <- exp(-0.5 * (vapply(fits, aicc, numeric(1)) - min(vapply(fits, aicc, numeric(1)))))
  w <- w / sum(w)
  b <- vapply(fits, function(m) { cc <- coef(m); if ("moisture" %in% names(cc)) cc["moisture"] else 0 }, numeric(1))
  in_m <- vapply(subsets, function(s) "moisture" %in% s, logical(1))
  c(avg_coef = sum(w * b), wplus = sum(w[in_m]))
}
avg_moisture(d09)

 avg_coef     wplus 
0.3098513 0.3927900

Watching it drift with collinearity

The single dataset is suggestive, but the honest claim is about expected behaviour, so we repeat the whole exercise across a range of correlations, averaging over many datasets at each level to smooth out the noise. At each correlation we record the partial estimate, the marginal estimate, the model-averaged estimate, the summed weight for moisture, and, as a check, the model-averaged prediction at a fixed site.

one_rep <- function(rho) {
  d <- gen_data(rho, seed = NA)
  fits <- lapply(forms, function(f) lm(as.formula(f), data = d))
  w <- exp(-0.5 * (vapply(fits, aicc, numeric(1)) - min(vapply(fits, aicc, numeric(1)))))
  w <- w / sum(w)
  in_m <- vapply(subsets, function(s) "moisture" %in% s, logical(1))
  bm <- vapply(fits, function(m) { cc <- coef(m); if ("moisture" %in% names(cc)) cc["moisture"] else 0 }, numeric(1))
  x0 <- data.frame(moisture = 1, carbon = 1, aspect = 0)
  c(partial  = as.numeric(coef(fits[[which(forms == "biomass ~ moisture + carbon")]])["moisture"]),
    marginal = as.numeric(coef(fits[[which(forms == "biomass ~ moisture")]])["moisture"]),
    averaged = sum(w * bm),
    wplus    = sum(w[in_m]),
    pred     = sum(w * vapply(fits, function(m) as.numeric(predict(m, x0)), numeric(1))))
}

rhos <- c(0, 0.3, 0.6, 0.9)
set.seed(505)
sweep <- t(vapply(rhos, function(r) {
  colMeans(t(vapply(seq_len(300), function(i) one_rep(r), numeric(5))))
}, numeric(5)))
sweep <- data.frame(rho = rhos, sweep)
round(sweep, 3)

  rho partial marginal averaged wplus  pred
1 0.0   0.991    0.963    0.980 0.980 4.386
2 0.3   1.014    1.429    1.001 0.976 4.398
3 0.6   1.005    1.825    0.972 0.942 4.386
4 0.9   1.025    2.274    0.929 0.674 4.398

lo <- sweep[sweep$rho == 0, ]
hi <- sweep[sweep$rho == 0.9, ]

At near-zero correlation, all three coefficient estimates agree near the true 1.0 and the summed weight for moisture sits at 0.98. As the correlation climbs to 0.9, the partial estimate stays put at 1.02, but the marginal estimate inflates to 2.27, and the model-averaged estimate drifts down to 0.93 as the weight shifts among models that mean different things by “moisture”. The summed weight for moisture, meanwhile, collapses from 0.98 to 0.67, which would read as moisture becoming less important even though its true effect never changed.

library(ggplot2)
long <- data.frame(
  rho = rep(sweep$rho, 3),
  estimate = c(sweep$partial, sweep$marginal, sweep$averaged),
  kind = rep(c("partial (full model)", "marginal (moisture only)", "model-averaged"), each = nrow(sweep)))
ggplot(long, aes(rho, estimate, colour = kind)) +
  geom_hline(yintercept = 1.0, linetype = "dashed", colour = "grey55") +
  geom_line(linewidth = 0.8) + geom_point(size = 2) +
  scale_colour_manual(values = c("partial (full model)" = "#4c72b0",
                                 "marginal (moisture only)" = "#c44e52",
                                 "model-averaged" = "#55a868"), name = NULL) +
  labs(x = "correlation between moisture and carbon", y = "estimated moisture coefficient") +
  theme_minimal(base_size = 12) + theme(legend.position = "top")

Line chart of three coefficient estimates against correlation. The partial line stays near one, the marginal line rises steeply, and the averaged line drifts, while a dashed line marks the true value of one. — Figure 1: The moisture coefficient under increasing collinearity: the partial estimate tracks the truth, the marginal estimate inflates, and the model-averaged estimate drifts between them. The true value (dashed) is fixed at 1.0.

The prediction is fine the whole time

Here is the reassuring half. The same candidate set, averaged the same way, gives a model-averaged prediction that barely moves as collinearity increases, because prediction does not depend on attributing shared variation to one gradient or the other. The models disagree about how to split the effect between moisture and carbon, but they agree about the total at any given site.

ggplot(sweep, aes(rho, pred)) +
  geom_hline(yintercept = 4.4, linetype = "dashed", colour = "grey55") +
  geom_line(colour = "#4c72b0", linewidth = 0.8) + geom_point(colour = "#4c72b0", size = 2) +
  coord_cartesian(ylim = c(3.8, 5.0)) +
  labs(x = "correlation between moisture and carbon", y = "model-averaged prediction at site") +
  theme_minimal(base_size = 12)

Line chart of the model-averaged prediction against correlation, staying flat near four point four across the whole range, with a dashed line at the true value. — Figure 2: The model-averaged prediction at a fixed site (moisture and carbon both one standard deviation up) stays close to its true value of 4.4 across all correlations, even as the averaged coefficient drifts.

Why coefficients cannot be averaged the way predictions can

Cade’s point, in plain terms, is that collinearity changes the scale on which a coefficient lives. The regression coefficient for moisture is a change in biomass per unit of moisture holding the other terms fixed, but when carbon is dropped, the relevant “unit” of moisture now carries some of carbon’s variation with it. The denominators are no longer the same across models, so neither the parameters nor their estimates share a common scale, and averaging numbers on different scales produces something with no clean interpretation. Standardising by the partial standard deviation of each predictor, as Cade suggests, can make the scales commensurate, which is a necessary step before averaging coefficients is even sensible, and often it reveals that averaging still should not be done.

The summed weight has the same disease from a different angle. Galipaud and colleagues showed it takes a wide range of values for predictors that are collinear or even unrelated to the response, so a falling w+ under collinearity says more about the candidate set than about the ecology.

Predictions escape all of this because they are evaluated at a point, on the response scale, where the split of variation between correlated predictors is irrelevant to the total.

What to do instead

If the goal is prediction, average predictions and report the unconditional interval; the coefficient pathology never touches you. If the goal is the effect of a specific predictor, do not read it off a model-averaged coefficient or a summed weight under collinearity. Report the partial effect from a single, well-justified model and be explicit that it is conditional on the other terms, or move to a method built for correlated predictors, such as ridge or elastic-net regression, which shrinks in a controlled way rather than averaging across incompatible scales. Model averaging is a tool for prediction under model uncertainty, not a shortcut to clean effect sizes.

References

Cade BS (2015) Ecology 96(9):2370-2382 (10.1890/14-1639.1)

Galipaud M, Gillingham MAF, David M, Dechaume-Moncharmont FX (2014) Methods in Ecology and Evolution 5(10):983-991 (10.1111/2041-210X.12251)

Grueber CE, Nakagawa S, Laws RJ, Jamieson IG (2011) Journal of Evolutionary Biology 24(4):699-711 (10.1111/j.1420-9101.2010.02210.x)

--- title: "Model-averaged coefficients: a caution" description: "Model-averaged regression coefficients mislead under collinearity in ecology (Cade 2015). Average predictions, not coefficients: a worked R demonstration." date: "2026-07-23 13:00" categories: [R, multimodel inference, model averaging, collinearity, ggplot2, ecology tutorial] image: thumbnail.png image-alt: "Three estimates of one predictor's coefficient diverging as collinearity increases, while the true value stays fixed." --- Averaging predictions across a candidate set is safe. Averaging coefficients across the same set is where multi-model inference most often goes wrong in ecology, and it goes wrong precisely when it looks most useful: when predictors are correlated and you want a single number for each one's effect. The problem is not a coding slip. It is that a regression coefficient means different things in different models. In a model that also contains a correlated predictor, a coefficient is a partial effect, holding that other predictor fixed. In a model that omits it, the same coefficient absorbs the shared variation and becomes something larger. Averaging those two numbers, weighted by Akaike weights, blends quantities that are not on the same footing. Cade set this out in detail in 2015; this post reproduces the core of the argument in R and shows why the prediction stays trustworthy even as the averaged coefficient does not. The [Akaike weights](../akaike-weights/) post is the prerequisite, and the [collinearity and variance inflation](../collinearity-and-vif/) post covers the mechanism this one builds on. ## Two correlated gradients Eighty forest plots, log above-ground biomass as the response. Two environmental predictors carry real effects: soil moisture and soil organic carbon, which in real landscapes tend to rise together. A third predictor, slope aspect, has no effect. We control the correlation between moisture and organic carbon directly so we can watch what it does. ```{r gen} gen_data <- function(rho, N = 80, seed = 505) { if (!is.na(seed)) set.seed(seed) # seed = NA gives a fresh draw each call z1 <- rnorm(N); z2 <- rnorm(N) moisture <- z1 carbon <- rho * z1 + sqrt(1 - rho^2) * z2 # correlation with moisture = rho aspect <- rnorm(N) biomass <- 2 + 1.0 * moisture + 1.4 * carbon + 0 * aspect + rnorm(N, 0, 2.0) data.frame(biomass, moisture, carbon, aspect) } ``` The true coefficient of moisture is 1.0 throughout, regardless of how strongly it correlates with carbon. Any change we see in the estimate is an artefact of the modelling, not the ecology. ## One coefficient, three different answers Take a single dataset with strong correlation, and estimate the moisture coefficient two ways: from the model that also includes carbon (a partial effect), and from the model with moisture alone (a marginal effect). ```{r oneset} d09 <- gen_data(rho = 0.9) b_partial <- coef(lm(biomass ~ moisture + carbon, data = d09))["moisture"] b_marginal <- coef(lm(biomass ~ moisture, data = d09))["moisture"] c(partial = as.numeric(b_partial), marginal = as.numeric(b_marginal)) ``` The two estimates disagree sharply. The partial effect is `r sprintf("%.2f", b_partial)`; the marginal is `r sprintf("%.2f", b_marginal)`, well over three times larger, because with carbon left out the moisture term stands in for both gradients at once. On this one dataset the partial estimate sits below the true value of 1.0, a reminder that any single fit is a noisy draw, especially with strong collinearity inflating its standard error; the systematic behaviour comes from averaging over many datasets, which we do shortly. What holds on every draw is the gap: neither estimate is wrong for what it is, they simply answer different questions, and a model-averaged coefficient mixes them. To see the averaged coefficient itself, we need the full candidate set. With three predictors there are eight models; the averaged moisture coefficient adds up each model's moisture estimate weighted by its Akaike weight, counting models without moisture as contributing zero. ```{r avgcoef} aicc <- function(m) { k <- length(coef(m)) + 1; N <- nobs(m); AIC(m) + 2 * k * (k + 1) / (N - k - 1) } vars <- c("moisture", "carbon", "aspect") subsets <- unlist(lapply(0:3, function(k) combn(vars, k, simplify = FALSE)), recursive = FALSE) forms <- vapply(subsets, function(s) if (length(s) == 0) "biomass ~ 1" else paste("biomass ~", paste(s, collapse = " + ")), character(1)) avg_moisture <- function(d) { fits <- lapply(forms, function(f) lm(as.formula(f), data = d)) w <- exp(-0.5 * (vapply(fits, aicc, numeric(1)) - min(vapply(fits, aicc, numeric(1))))) w <- w / sum(w) b <- vapply(fits, function(m) { cc <- coef(m); if ("moisture" %in% names(cc)) cc["moisture"] else 0 }, numeric(1)) in_m <- vapply(subsets, function(s) "moisture" %in% s, logical(1)) c(avg_coef = sum(w * b), wplus = sum(w[in_m])) } avg_moisture(d09) ``` ## Watching it drift with collinearity The single dataset is suggestive, but the honest claim is about expected behaviour, so we repeat the whole exercise across a range of correlations, averaging over many datasets at each level to smooth out the noise. At each correlation we record the partial estimate, the marginal estimate, the model-averaged estimate, the summed weight for moisture, and, as a check, the model-averaged prediction at a fixed site. ```{r sweep} one_rep <- function(rho) { d <- gen_data(rho, seed = NA) fits <- lapply(forms, function(f) lm(as.formula(f), data = d)) w <- exp(-0.5 * (vapply(fits, aicc, numeric(1)) - min(vapply(fits, aicc, numeric(1))))) w <- w / sum(w) in_m <- vapply(subsets, function(s) "moisture" %in% s, logical(1)) bm <- vapply(fits, function(m) { cc <- coef(m); if ("moisture" %in% names(cc)) cc["moisture"] else 0 }, numeric(1)) x0 <- data.frame(moisture = 1, carbon = 1, aspect = 0) c(partial = as.numeric(coef(fits[[which(forms == "biomass ~ moisture + carbon")]])["moisture"]), marginal = as.numeric(coef(fits[[which(forms == "biomass ~ moisture")]])["moisture"]), averaged = sum(w * bm), wplus = sum(w[in_m]), pred = sum(w * vapply(fits, function(m) as.numeric(predict(m, x0)), numeric(1)))) } rhos <- c(0, 0.3, 0.6, 0.9) set.seed(505) sweep <- t(vapply(rhos, function(r) { colMeans(t(vapply(seq_len(300), function(i) one_rep(r), numeric(5)))) }, numeric(5))) sweep <- data.frame(rho = rhos, sweep) round(sweep, 3) ``` ```{r sweepnums} lo <- sweep[sweep$rho == 0, ] hi <- sweep[sweep$rho == 0.9, ] ``` At near-zero correlation, all three coefficient estimates agree near the true 1.0 and the summed weight for moisture sits at `r sprintf("%.2f", lo$wplus)`. As the correlation climbs to 0.9, the partial estimate stays put at `r sprintf("%.2f", hi$partial)`, but the marginal estimate inflates to `r sprintf("%.2f", hi$marginal)`, and the model-averaged estimate drifts down to `r sprintf("%.2f", hi$averaged)` as the weight shifts among models that mean different things by "moisture". The summed weight for moisture, meanwhile, collapses from `r sprintf("%.2f", lo$wplus)` to `r sprintf("%.2f", hi$wplus)`, which would read as moisture becoming less important even though its true effect never changed. ```{r fig-coef} #| fig-cap: "The moisture coefficient under increasing collinearity: the partial estimate tracks the truth, the marginal estimate inflates, and the model-averaged estimate drifts between them. The true value (dashed) is fixed at 1.0." #| fig-alt: "Line chart of three coefficient estimates against correlation. The partial line stays near one, the marginal line rises steeply, and the averaged line drifts, while a dashed line marks the true value of one." #| fig-width: 6.2 #| fig-height: 3.8 library(ggplot2) long <- data.frame( rho = rep(sweep$rho, 3), estimate = c(sweep$partial, sweep$marginal, sweep$averaged), kind = rep(c("partial (full model)", "marginal (moisture only)", "model-averaged"), each = nrow(sweep))) ggplot(long, aes(rho, estimate, colour = kind)) + geom_hline(yintercept = 1.0, linetype = "dashed", colour = "grey55") + geom_line(linewidth = 0.8) + geom_point(size = 2) + scale_colour_manual(values = c("partial (full model)" = "#4c72b0", "marginal (moisture only)" = "#c44e52", "model-averaged" = "#55a868"), name = NULL) + labs(x = "correlation between moisture and carbon", y = "estimated moisture coefficient") + theme_minimal(base_size = 12) + theme(legend.position = "top") ``` ## The prediction is fine the whole time Here is the reassuring half. The same candidate set, averaged the same way, gives a model-averaged prediction that barely moves as collinearity increases, because prediction does not depend on attributing shared variation to one gradient or the other. The models disagree about how to split the effect between moisture and carbon, but they agree about the total at any given site. ```{r fig-pred} #| fig-cap: "The model-averaged prediction at a fixed site (moisture and carbon both one standard deviation up) stays close to its true value of 4.4 across all correlations, even as the averaged coefficient drifts." #| fig-alt: "Line chart of the model-averaged prediction against correlation, staying flat near four point four across the whole range, with a dashed line at the true value." #| fig-width: 6.2 #| fig-height: 3.4 ggplot(sweep, aes(rho, pred)) + geom_hline(yintercept = 4.4, linetype = "dashed", colour = "grey55") + geom_line(colour = "#4c72b0", linewidth = 0.8) + geom_point(colour = "#4c72b0", size = 2) + coord_cartesian(ylim = c(3.8, 5.0)) + labs(x = "correlation between moisture and carbon", y = "model-averaged prediction at site") + theme_minimal(base_size = 12) ``` ## Why coefficients cannot be averaged the way predictions can Cade's point, in plain terms, is that collinearity changes the scale on which a coefficient lives. The regression coefficient for moisture is a change in biomass per unit of moisture holding the other terms fixed, but when carbon is dropped, the relevant "unit" of moisture now carries some of carbon's variation with it. The denominators are no longer the same across models, so neither the parameters nor their estimates share a common scale, and averaging numbers on different scales produces something with no clean interpretation. Standardising by the partial standard deviation of each predictor, as Cade suggests, can make the scales commensurate, which is a necessary step before averaging coefficients is even sensible, and often it reveals that averaging still should not be done. The summed weight has the same disease from a different angle. Galipaud and colleagues showed it takes a wide range of values for predictors that are collinear or even unrelated to the response, so a falling w+ under collinearity says more about the candidate set than about the ecology. Predictions escape all of this because they are evaluated at a point, on the response scale, where the split of variation between correlated predictors is irrelevant to the total. ## What to do instead If the goal is prediction, average predictions and report the unconditional interval; the coefficient pathology never touches you. If the goal is the effect of a specific predictor, do not read it off a model-averaged coefficient or a summed weight under collinearity. Report the partial effect from a single, well-justified model and be explicit that it is conditional on the other terms, or move to a method built for correlated predictors, such as ridge or elastic-net regression, which shrinks in a controlled way rather than averaging across incompatible scales. Model averaging is a tool for prediction under model uncertainty, not a shortcut to clean effect sizes. ## Related tutorials - [Akaike weights and evidence ratios](../akaike-weights/) - [Model-averaged predictions](../model-averaged-predictions/) - [Collinearity and variance inflation](../collinearity-and-vif/) - [Ridge regression for collinearity](../ridge-regression-for-collinearity/) ## References Cade BS (2015) Ecology 96(9):2370-2382 (10.1890/14-1639.1) Galipaud M, Gillingham MAF, David M, Dechaume-Moncharmont FX (2014) Methods in Ecology and Evolution 5(10):983-991 (10.1111/2041-210X.12251) Grueber CE, Nakagawa S, Laws RJ, Jamieson IG (2011) Journal of Evolutionary Biology 24(4):699-711 (10.1111/j.1420-9101.2010.02210.x)