Random-effects meta-analysis in R

meta-analysis
ecology tutorial
R
effect size
Pool effect sizes across ecological studies with a random-effects model in base R: Hedges’ g, DerSimonian-Laird and REML tau-squared, forest plots and coverage.
Author

Tidy Ecology

Published

2026-07-24

A meta-analysis combines effect sizes from many studies into a single summary, with a measure of how much they disagree. In ecology the studies usually differ in site, species, and design, so the pooled estimate has to carry that spread rather than paper over it. This post builds a random-effects meta-analysis from scratch in base R and shows why the simpler fixed-effect model gives intervals that are too narrow once real between-study heterogeneity is present.

We use the standardised mean difference (Hedges’ g) as the effect metric: the difference between a treatment group and a control group in pooled standard-deviation units, with a small-sample correction. Each study reports its own g and a within-study sampling variance v that we treat as known. The two models differ in one assumption. The fixed-effect (common-effect) model says every study estimates the same true effect and the only scatter is sampling noise. The random-effects model says each study has its own true effect drawn from a distribution with mean mu and variance tau-squared, so the summary target is the mean of that distribution.

Simulating a set of studies

We simulate 20 studies. Each has a true effect drawn from a normal distribution with mean 0.4 and between-study standard deviation 0.25 (so tau-squared is 0.0625). Sample sizes vary from study to study, which is what gives a meta-analysis its mix of precise and imprecise estimates. For each study we generate the raw group data, compute Cohen’s d from the pooled standard deviation, apply the Hedges small-sample correction J, and record the large-sample variance of g.

set.seed(151)
mu_true <- 0.40; tau_true <- 0.25; k <- 20
n <- sample(15:120, k, replace = TRUE)          # per-study sample size (per arm)

sim_study <- function(theta, ni) {
  x1 <- rnorm(ni, 0, 1)                          # control arm
  x2 <- rnorm(ni, theta, 1)                      # treatment arm
  sp <- sqrt(((ni - 1) * sd(x1)^2 + (ni - 1) * sd(x2)^2) / (2 * ni - 2))
  d  <- (mean(x2) - mean(x1)) / sp               # Cohen's d
  J  <- 1 - 3 / (4 * (2 * ni - 2) - 1)           # small-sample correction
  g  <- J * d                                    # Hedges' g
  vd <- (2 * ni) / (ni * ni) + d^2 / (2 * (2 * ni))
  c(g = g, v = J^2 * vd)                         # g and its sampling variance
}

theta_i <- rnorm(k, mu_true, tau_true)           # each study's own true effect
gv <- t(mapply(sim_study, theta_i, n))
g  <- gv[, "g"]; v <- gv[, "v"]
se <- sqrt(v)
round(cbind(n, g = g, se = se)[1:6, ], 3)
       n     g    se
[1,]  76 0.964 0.171
[2,]  27 0.295 0.270
[3,]  63 0.391 0.179
[4,] 114 0.642 0.135
[5,]  32 0.677 0.254
[6,]  48 1.133 0.218

Fixed-effect model

The fixed-effect estimate is a precision-weighted average, with weights equal to the inverse sampling variance. Precise studies (large n, small v) count more.

wf    <- 1 / v
mu_fe <- sum(wf * g) / sum(wf)
se_fe <- sqrt(1 / sum(wf))
ci_fe <- mu_fe + c(-1, 1) * 1.96 * se_fe
c(estimate = mu_fe, se = se_fe, lower = ci_fe[1], upper = ci_fe[2])
  estimate         se      lower      upper 
0.42620400 0.03763708 0.35243533 0.49997268

The fixed-effect summary is 0.426 with a 95% interval from 0.352 to 0.5, a width of 0.148. That interval is only honest if the studies really do share one true effect. Cochran’s Q tests that assumption: it is the weighted sum of squared deviations of each study from the fixed-effect summary, compared with a chi-squared distribution on k minus one degrees of freedom.

Q   <- sum(wf * (g - mu_fe)^2)
dfQ <- k - 1
pQ  <- pchisq(Q, dfQ, lower.tail = FALSE)
c(Q = Q, df = dfQ, p = pQ)
           Q           df            p 
6.785947e+01 1.900000e+01 2.079198e-07

Here Q is 67.9 on 19 degrees of freedom (p = 2.1^{-7}), so the studies disagree by more than sampling noise alone. The fixed-effect interval ignores that disagreement and will be too narrow.

Random-effects: DerSimonian-Laird

The random-effects model adds a between-study variance tau-squared to every study’s weight, so the weights become one over v plus tau-squared. The classic moment estimator of DerSimonian and Laird (1986) reads tau-squared straight off the Q statistic.

Cc      <- sum(wf) - sum(wf^2) / sum(wf)
tau2_dl <- max(0, (Q - dfQ) / Cc)                # DerSimonian-Laird moment estimate
ws      <- 1 / (v + tau2_dl)                     # random-effects weights
mu_re   <- sum(ws * g) / sum(ws)
se_re   <- sqrt(1 / sum(ws))
ci_re   <- mu_re + c(-1, 1) * 1.96 * se_re
c(tau2 = tau2_dl, estimate = mu_re, se = se_re, lower = ci_re[1], upper = ci_re[2])
      tau2   estimate         se      lower      upper 
0.07335959 0.42836456 0.07235156 0.28655551 0.57017361

Random-effects: REML

The moment estimator is quick but a maximum-likelihood approach is the modern default. Restricted maximum likelihood (REML) estimates tau-squared by maximising the restricted log-likelihood, which accounts for the estimation of mu (Viechtbauer 2005). It is a one-parameter search over tau-squared, easy to code with optimize.

reml_ll <- function(t2) {
  wi <- 1 / (v + t2)
  mh <- sum(wi * g) / sum(wi)
  -0.5 * sum(log(v + t2)) - 0.5 * log(sum(wi)) - 0.5 * sum(wi * (g - mh)^2)
}
opt       <- optimize(reml_ll, c(0, 5), maximum = TRUE)
tau2_reml <- opt$maximum
wr        <- 1 / (v + tau2_reml)
mu_reml   <- sum(wr * g) / sum(wr)
se_reml   <- sqrt(1 / sum(wr))
ci_reml   <- mu_reml + c(-1, 1) * 1.96 * se_reml
c(tau2 = tau2_reml, estimate = mu_reml, se = se_reml)
      tau2   estimate         se 
0.07302399 0.42835821 0.07223386

The two estimators of tau-squared agree closely here: 0.073 from DerSimonian-Laird and 0.073 from REML, both near the true 0.0625. The random-effects summary is 0.428 with a 95% interval from 0.287 to 0.57. That interval is 1.92 times as wide as the fixed-effect one, because the between-study variance flows into the standard error. The point estimate barely moves; what changes is the honesty of the uncertainty.

ord <- order(g)
fdat <- data.frame(
  lab = factor(paste0("S", seq_len(k))[ord], levels = paste0("S", seq_len(k))[ord]),
  est = g[ord], lo = g[ord] - 1.96 * se[ord], hi = g[ord] + 1.96 * se[ord])
sdat <- data.frame(
  lab = factor(c("RE (random)", "FE (fixed)"),
               levels = c("RE (random)", "FE (fixed)")),
  est = c(mu_re, mu_fe), lo = c(ci_re[1], ci_fe[1]), hi = c(ci_re[2], ci_fe[2]))

ggplot(fdat, aes(est, lab)) +
  geom_vline(xintercept = 0, linetype = "dashed", colour = pal$faint) +
  geom_pointrange(aes(xmin = lo, xmax = hi), colour = pal$forest,
                  linewidth = 0.4, fatten = 1.6) +
  geom_pointrange(data = sdat, aes(xmin = lo, xmax = hi), colour = pal$warn,
                  shape = 18, linewidth = 0.9, fatten = 3.2) +
  labs(x = "Hedges' g (standardised mean difference)", y = NULL,
       title = "Study effects and two pooled summaries") +
  theme_te()
Forest plot: twenty horizontal point-and-whisker rows for the studies, sorted by effect size, above two summary rows. The fixed-effect summary near 0.43 has a narrow interval; the random-effects summary at the same point has a visibly wider interval. A dashed vertical line marks zero effect.
Figure 1: Forest plot of 20 simulated studies (Hedges’ g with 95% intervals), ordered by effect size, with the fixed-effect and random-effects pooled estimates. The random-effects interval is wider because it carries the between-study variance.

Why the intervals differ: a coverage check

The claim that the fixed-effect interval is too narrow is testable. We repeat the whole exercise many times at several true values of tau, and record how often each model’s 95% interval contains the true mu. A well-behaved interval should contain the truth about 95% of the time. For speed the coverage loop uses the DerSimonian-Laird estimator, which has a closed form.

cover_at <- function(tt, nsim = 1200) {
  fe <- 0; re <- 0
  for (s in seq_len(nsim)) {
    th <- rnorm(k, mu_true, tt)
    gg <- numeric(k); vv <- numeric(k)
    for (i in seq_len(k)) { r <- sim_study(th[i], n[i]); gg[i] <- r[1]; vv[i] <- r[2] }
    wF <- 1 / vv; mF <- sum(wF * gg) / sum(wF); sF <- sqrt(1 / sum(wF))
    if (abs(mF - mu_true) <= 1.96 * sF) fe <- fe + 1
    QQ <- sum(wF * (gg - mF)^2); CC <- sum(wF) - sum(wF^2) / sum(wF)
    t2 <- max(0, (QQ - (k - 1)) / CC); wS <- 1 / (vv + t2)
    mR <- sum(wS * gg) / sum(wS); sR <- sqrt(1 / sum(wS))
    if (abs(mR - mu_true) <= 1.96 * sR) re <- re + 1
  }
  c(fe = fe / nsim, re = re / nsim)
}
set.seed(2718)
grid   <- c(0, 0.1, 0.2, 0.3, 0.4)
covtab <- as.data.frame(t(sapply(grid, cover_at)))
covtab$tau <- grid
round(covtab[, c("tau", "fe", "re")], 3)
  tau    fe    re
1 0.0 0.958 0.968
2 0.1 0.905 0.946
3 0.2 0.788 0.932
4 0.3 0.643 0.939
5 0.4 0.542 0.918

At tau equal to zero the two models agree and both sit near the nominal 0.95. As between-study variance grows the fixed-effect interval decays: its coverage falls from 0.958 to 0.542, so at the largest heterogeneity it misses the truth in nearly half of all meta-analyses. The random-effects interval holds much closer to nominal across the range (0.918 at the largest tau). The lesson is direct: with heterogeneous studies, a fixed-effect interval understates uncertainty, and the size of that error grows with the heterogeneity.

cdat <- rbind(
  data.frame(tau = grid, cover = covtab$fe, model = "Fixed-effect"),
  data.frame(tau = grid, cover = covtab$re, model = "Random-effects"))
ggplot(cdat, aes(tau, cover, colour = model)) +
  geom_hline(yintercept = 0.95, linetype = "dashed", colour = pal$faint) +
  geom_line(linewidth = 0.8) + geom_point(size = 2.2) +
  scale_colour_manual(values = c("Fixed-effect" = pal$warn,
                                 "Random-effects" = pal$forest), name = NULL) +
  labs(x = "true between-study SD (tau)", y = "95% interval coverage",
       title = "Coverage of the pooled mean") +
  coord_cartesian(ylim = c(0.4, 1)) +
  theme_te()
Line chart with true between-study standard deviation on the x axis and interval coverage on the y axis. The fixed-effect line starts near 0.95 at zero and drops steeply to about 0.54. The random-effects line stays close to a dashed horizontal reference at 0.95 across the whole range.
Figure 2: Interval coverage of the pooled mean against true between-study standard deviation, from simulation. Fixed-effect coverage collapses as heterogeneity grows; random-effects coverage stays near the nominal 0.95.

Honest limits

The random-effects interval here uses a normal (Wald) approximation with an estimated tau-squared. That approximation is good but not perfect: at the largest heterogeneity with only 20 studies the coverage sits a little below nominal (0.918), the known slight anticonservatism of Wald intervals in small meta-analyses. The Knapp-Hartung adjustment, which uses a t-distribution and a rescaled variance, restores coverage and is worth preferring in practice. Two further points matter more than the estimator. First, the effect metric is a modelling choice: Hedges’ g, log response ratios, and correlation coefficients each have their own variance formula and their own interpretation, and they are not interchangeable. Second, the model assumes the studies are independent; multiple effects from one study or one research group break that assumption and call for a multilevel meta-analysis.

Where to go next

Heterogeneity is not a nuisance to be summarised in one number and forgotten. The next tutorial takes tau-squared apart with the Q statistic, I-squared, and the prediction interval, and shows why I-squared on its own can mislead. From there, meta-regression asks whether a study-level moderator explains some of the spread, and a funnel-plot check asks whether the set of studies is a fair sample of the evidence.

References

  • DerSimonian R, Laird N 1986. Controlled Clinical Trials 7(3):177-188 (10.1016/0197-2456(86)90046-2)
  • Hedges LV 1981. Journal of Educational Statistics 6(2):107-128 (10.3102/10769986006002107)
  • 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)
  • Gurevitch J, Koricheva J, Nakagawa S, Stewart G 2018. Nature 555(7695):175-182 (10.1038/nature25753)
  • Borenstein M, Hedges LV, Higgins JPT, Rothstein HR 2009. Introduction to Meta-Analysis. ISBN 978-0-470-05724-7