---
title: "Meta-regression with moderators"
description: "Explain heterogeneity with a study-level moderator in base R: weighted meta-regression, residual tau-squared, an R-squared analogue, and error control."
date: "2026-05-24 11:00"
categories: [meta-analysis, ecology tutorial, R, meta-regression]
image: thumbnail.png
image-alt: "Bubble plot of study effect against a moderator, point area scaled by precision, with the fitted line and its band sloping clearly upwards."
---
A random-effects model measures how much studies disagree, but not why. Meta-regression asks whether a study-level covariate, a moderator, accounts for some of the spread: latitude, study duration, mean body size, sampling method. It is a weighted regression of the effect sizes on the moderator, with the between-study variance kept as a residual term. This post builds it in base R, reports how much heterogeneity the moderator absorbs, and shows the common mistake of running a meta-regression without that residual term, which turns real heterogeneity into false moderator effects.
```{r}
#| label: setup
#| code-fold: true
#| code-summary: "Setup: packages, colours and plot theme"
#| results: hide
#| message: false
#| warning: false
library(ggplot2)
knitr::opts_chunk$set(echo = TRUE, message = FALSE, warning = FALSE,
fig.width = 7.5, fig.height = 4.5, dpi = 100)
pal <- list(ink = "#16241d", body = "#2c3a31", forest = "#275139",
label = "#46604a", sage = "#93a87f", paper = "#f5f4ee",
line = "#dad9ca", faint = "#5d6b61", gold = "#cda23f",
mown = "#2f8f63", warn = "#b5534e", low = "#c9b458", high = "#1d5b4e")
theme_te <- function(base_size = 12) {
theme_minimal(base_size = base_size) +
theme(panel.grid.minor = element_blank(),
panel.grid.major = element_line(colour = pal$line, linewidth = 0.3),
plot.background = element_rect(fill = pal$paper, colour = NA),
panel.background = element_rect(fill = pal$paper, colour = NA),
axis.text = element_text(colour = pal$body),
axis.title = element_text(colour = pal$ink),
plot.title = element_text(colour = pal$ink, face = "bold"),
plot.subtitle = element_text(colour = pal$faint),
legend.position = "bottom",
legend.text = element_text(colour = pal$body),
legend.title = element_text(colour = pal$ink))
}
```
```{r}
#| label: sci-tex-helper
#| include: false
sci_tex <- function(x, digits = 2) {
s <- formatC(x, format = "e", digits = digits)
paste0("$", sub("e.*", "", s), " \\times 10^{", as.integer(sub(".*e", "", s)), "}$")
}
```
## A moderator that drives the effect
We simulate `r 25` studies whose true effect depends on a standardised moderator x, with intercept `r 0.10` and slope `r 0.35`, plus residual between-study noise with standard deviation `r 0.15`. Each study reports a Hedges' g and its known sampling variance.
```{r}
#| label: simulate
sim_study <- function(theta, ni) {
x1 <- rnorm(ni, 0, 1); x2 <- rnorm(ni, theta, 1)
sp <- sqrt(((ni - 1) * sd(x1)^2 + (ni - 1) * sd(x2)^2) / (2 * ni - 2))
d <- (mean(x2) - mean(x1)) / sp
J <- 1 - 3 / (4 * (2 * ni - 2) - 1)
g <- J * d
vd <- (2 * ni) / (ni * ni) + d^2 / (2 * (2 * ni))
c(g = g, v = J^2 * vd)
}
set.seed(153)
k <- 25; b0 <- 0.10; b1 <- 0.35; tau_r <- 0.15
x <- scale(runif(k, -2, 2))[, 1] # standardised moderator
n <- sample(15:120, k, replace = TRUE)
theta_i <- b0 + b1 * x + rnorm(k, 0, tau_r) # true effect depends on x
gv <- t(mapply(sim_study, theta_i, n))
g <- gv[, "g"]; v <- gv[, "v"]
```
## Fitting the meta-regression
Given a residual between-study variance tau-squared, the coefficients are a generalised least-squares fit with weights equal to one over v plus tau-squared. We estimate tau-squared by REML, maximising the restricted log-likelihood over a one-parameter search, then read the coefficients and their standard errors from the weighted normal equations.
```{r}
#| label: meta-regression
mreg <- function(y, v, X) {
reml <- function(t2) {
W <- 1 / (v + t2); XtW <- t(X * W)
b <- solve(XtW %*% X, XtW %*% y); r <- y - X %*% b
-0.5 * sum(log(v + t2)) - 0.5 * log(det(XtW %*% X)) - 0.5 * sum(W * r^2)
}
t2 <- optimize(reml, c(0, 10), maximum = TRUE)$maximum
W <- 1 / (v + t2); XtW <- t(X * W)
Vb <- solve(XtW %*% X); b <- Vb %*% (XtW %*% y)
list(b = as.vector(b), se = sqrt(diag(Vb)), tau2 = t2, Vb = Vb)
}
X <- cbind(1, x)
fit <- mreg(g, v, X)
QM <- (fit$b[2] / fit$se[2])^2 # Wald test for the slope
pQM <- pchisq(QM, 1, lower.tail = FALSE)
c(intercept = fit$b[1], slope = fit$b[2], se_slope = fit$se[2], QM = QM, p = pQM)
```
The estimated slope is `r round(fit$b[2], 3)` with a standard error of `r round(fit$se[2], 3)`, recovering the true `r 0.35`. The Wald test gives QM equal to `r round(QM, 1)` on one degree of freedom (p = `r sci_tex(pQM, 1)`), so the moderator is clearly related to the effect.
## How much heterogeneity does it explain?
Fit an intercept-only random-effects model to get the total between-study variance, then compare it with the residual left after the moderator. The proportion removed is an R-squared analogue for meta-analysis.
```{r}
#| label: r2-analog
wf <- 1 / v; mu_fe <- sum(wf * g) / sum(wf)
Q <- sum(wf * (g - mu_fe)^2); Cc <- sum(wf) - sum(wf^2) / sum(wf)
tau2_tot <- max(0, (Q - (k - 1)) / Cc) # total heterogeneity
tau2_res <- fit$tau2 # residual after moderator
R2 <- max(0, (tau2_tot - tau2_res) / tau2_tot)
c(tau2_total = tau2_tot, tau2_residual = tau2_res, R2_analog = R2)
```
Total between-study variance is `r round(tau2_tot, 3)`; after fitting the moderator the residual falls to `r round(tau2_res, 3)`, so the moderator accounts for about `r sprintf("%.0f%%", 100 * R2)` of the heterogeneity. The rest is genuine residual variation among studies that share the same moderator value.
```{r}
#| label: fig-mr-bubble
#| fig-cap: "Bubble plot of study effect against the moderator, with point area proportional to study precision and the fitted meta-regression line with its 95% band. Larger, more precise studies pull the line more strongly."
#| fig-alt: "Scatter plot of Hedges' g against a standardised moderator. Points of varying size, larger for more precise studies, rise from lower left to upper right. A forest-green fitted line with a shaded confidence band runs through them with a clear positive slope."
xg <- seq(min(x), max(x), length.out = 100)
Xg <- cbind(1, xg)
fitv <- as.vector(Xg %*% fit$b)
sef <- sqrt(rowSums((Xg %*% fit$Vb) * Xg))
band <- data.frame(x = xg, y = fitv, lo = fitv - 1.96 * sef, hi = fitv + 1.96 * sef)
bub <- data.frame(x = x, g = g, prec = 1 / (v + tau2_res))
ggplot(bub, aes(x, g)) +
geom_ribbon(data = band, aes(x, ymin = lo, ymax = hi), inherit.aes = FALSE,
fill = pal$sage, alpha = 0.30) +
geom_line(data = band, aes(x, y), inherit.aes = FALSE,
colour = pal$forest, linewidth = 0.9) +
geom_point(aes(size = prec), colour = pal$forest, alpha = 0.55) +
scale_size_area(max_size = 8, guide = "none") +
labs(x = "moderator (standardised)", y = "Hedges' g",
title = "Meta-regression: effect against a moderator") +
theme_te()
```
## The mistake: ignoring the residual variance
A tempting shortcut is to weight the regression only by within-study precision, one over v, and leave out the residual between-study variance. That is a fixed-effect meta-regression. When there is real residual heterogeneity, its standard errors are too small and the moderator test rejects far too often. We check this with a null moderator, one that has no true effect, and count how often each method declares it significant as the residual variance grows.
```{r}
#| label: type-one
freg <- function(y, v, X) { # fixed-effect meta-regression
W <- 1 / v; XtW <- t(X * W); Vb <- solve(XtW %*% X); b <- Vb %*% (XtW %*% y)
list(b = as.vector(b), se = sqrt(diag(Vb)))
}
type_one <- function(tt, nsim = 1000) {
fx <- 0; rx <- 0
for (s in seq_len(nsim)) {
xx <- scale(runif(k, -2, 2))[, 1]; nn <- sample(15:120, k, replace = TRUE)
thh <- b0 + 0 * xx + rnorm(k, 0, tt) # NULL slope
gg <- numeric(k); vv <- numeric(k)
for (i in seq_len(k)) { r <- sim_study(thh[i], nn[i]); gg[i] <- r[1]; vv[i] <- r[2] }
XX <- cbind(1, xx)
ff <- freg(gg, vv, XX); if (abs(ff$b[2] / ff$se[2]) > 1.96) fx <- fx + 1
rr <- tryCatch(mreg(gg, vv, XX), error = function(e) NULL)
if (!is.null(rr) && abs(rr$b[2] / rr$se[2]) > 1.96) rx <- rx + 1
}
c(fixed = fx / nsim, random = rx / nsim)
}
set.seed(4242)
grid <- c(0, 0.1, 0.2, 0.3)
t1tab <- as.data.frame(t(sapply(grid, type_one)))
t1tab$tau <- grid
round(t1tab[, c("tau", "fixed", "random")], 3)
```
With no residual heterogeneity both methods sit at or below the nominal five per cent, `r sprintf("%.3f", t1tab$fixed[t1tab$tau == 0])` for the fixed-effect fit and `r sprintf("%.3f", t1tab$random[t1tab$tau == 0])` for the random-effects one, the latter already a little conservative because tau-squared has to be estimated from `r k` studies. As residual variation grows, the fixed-effect meta-regression climbs to a rejection rate of `r sprintf("%.1f%%", 100 * t1tab$fixed[t1tab$tau == 0.3])`, so a moderator with no real effect looks significant in more than a third of analyses. The random-effects meta-regression stays near five per cent throughout (`r sprintf("%.1f%%", 100 * t1tab$random[t1tab$tau == 0.3])` at the largest residual variance). Leaving the residual term out does not make the analysis sharper; it manufactures moderator effects out of ordinary heterogeneity.
```{r}
#| label: fig-mr-typeone
#| fig-cap: "False-positive rate for a null moderator against residual between-study standard deviation. The fixed-effect meta-regression inflates badly as residual heterogeneity grows; the random-effects version holds near the nominal 0.05."
#| fig-alt: "Line chart with residual between-study standard deviation on the x axis and false-positive rate on the y axis. The fixed-effect line starts near 0.05 and rises steeply to about 0.37. The random-effects line stays close to a dashed horizontal reference at 0.05 across the range."
tdat <- rbind(
data.frame(tau = grid, rate = t1tab$fixed, model = "Fixed-effect meta-reg"),
data.frame(tau = grid, rate = t1tab$random, model = "Random-effects meta-reg"))
ggplot(tdat, aes(tau, rate, colour = model)) +
geom_hline(yintercept = 0.05, linetype = "dashed", colour = pal$faint) +
geom_line(linewidth = 0.8) + geom_point(size = 2.2) +
scale_colour_manual(values = c("Fixed-effect meta-reg" = pal$warn,
"Random-effects meta-reg" = pal$forest), name = NULL) +
labs(x = "residual between-study SD", y = "false-positive rate",
title = "Null moderator: false positives") +
theme_te()
```
## Honest limits
Even the random-effects meta-regression is slightly liberal at the largest residual variance with only `r k` studies, the same small-sample effect seen with the pooled mean; the Knapp-Hartung adjustment, which uses a t-distribution and a rescaled variance, tightens the test and is the safer default (its coverage for the pooled mean with four to ten studies is measured in [random-effects meta-analysis](../random-effects-meta-analysis/)). Two cautions matter more. Meta-regression with few studies over-fits easily, and a handful of moderators tested one after another will turn up a spurious winner; keep the number of moderators small relative to the number of studies. And an observed moderator is not a manipulated one: study-level covariates are often correlated with each other and with design features, so a significant moderator is an association, not a demonstrated cause.
## Where to go next
Meta-regression assumes the studies in hand are a fair sample of the evidence. If small studies with null or inconvenient results never reached publication, every quantity so far, the pooled mean, the heterogeneity, and the moderator slopes, can be biased. The next tutorial in this thread builds the funnel plot, Egger's test, and trim-and-fill to check for that, and the one after it fits the selection mechanism directly.
## References
- Thompson SG, Higgins JPT 2002. Statistics in Medicine 21(11):1559-1573 (10.1002/sim.1187)
- Knapp G, Hartung J 2003. Statistics in Medicine 22(17):2693-2710 (10.1002/sim.1482)
- Viechtbauer W 2005. Journal of Educational and Behavioral Statistics 30(3):261-293 (10.3102/10769986030003261)
- Nakagawa S, Santos ESA 2012. Evolutionary Ecology 26(5):1253-1274 (10.1007/s10682-012-9555-5)
- Borenstein M, Hedges LV, Higgins JPT, Rothstein HR 2009. Introduction to Meta-Analysis. ISBN 978-0-470-05724-7
## Related tutorials
- [Random-effects meta-analysis in R](../random-effects-meta-analysis/)
- [Heterogeneity in meta-analysis](../heterogeneity-in-meta-analysis/)
- [Correlations from gradients of different length](../correlations-from-gradients-of-different-length/)
- [Year trends in effect sizes and the SE covariate](../year-trends-in-effect-sizes/)
