Meta-Analysis in R with metafor: rma(), forest(), and funnel()

The metafor package, maintained by Wolfgang Viechtbauer and distributed on CRAN, is the R ecosystem’s most widely used tool for running a meta-analysis from raw effect-size data through to a publication-ready forest plot. Most narrative treatments of meta-analysis explain the statistics — heterogeneity, inverse-variance weighting, fixed- versus random-effects models — without showing what actually happens on screen when you run the analysis. This guide walks the three functions a researcher touches most: rma() to fit the model, forest() to plot it, and funnel() to check for small-study effects.

Before rma(): getting effect sizes into yi and vi

rma() does not take raw study data (means, counts, sample sizes) directly — it expects a vector of observed effect sizes (yi) and their corresponding sampling variances (vi, or standard errors via sei). metafor’s companion function escalc() computes these from raw study-level data for a wide range of effect-size types, set via its measure argument: "SMD" (standardized mean difference / Hedges’ g, for continuous outcomes), "MD" (raw mean difference), "OR", "RR", and "RD" (odds ratio, risk ratio, and risk difference, for 2×2 binary-outcome tables), and "ZCOR" (Fisher’s z-transformed correlation), among others. A typical call looks like:

dat <- escalc(measure = "SMD", m1i = mean_treat, sd1i = sd_treat, n1i = n_treat,
              m2i = mean_ctrl, sd2i = sd_ctrl, n2i = n_ctrl, data = raw_data)

This appends yi and vi columns to your data frame — the two arguments rma() actually needs. If your source studies already report an effect size and its variance or CI directly (common when pooling from a systematic review’s extraction table), you can skip escalc() and hand rma() the columns directly, as long as the effect-size definition is consistent across every included study. See CASRAI’s guide on standardized mean difference and Hedges’ g for the underlying formula.

Fitting the model with rma()

rma() (formally rma.uni()) fits equal-, fixed-, and random-effects meta-analytic models, and mixed-effects meta-regression models when moderators are added via the mods argument. The method argument controls which model is fit:

• method = "FE" fits a fixed-effect model, weighting each study by the inverse of its variance.
• Any of "DL" (DerSimonian-Laird), "HE" (Hedges), "SJ" (Sidik-Jonkman), "ML" (maximum likelihood), "PM" (Paule-Mandel), or "REML" (restricted maximum likelihood) fits a random-effects model using that method to estimate the between-study variance, τ². REML is the package default when method is left unspecified.

A minimal random-effects call:

res <- rma(yi, vi, data = dat)          # REML random-effects, the default
res <- rma(yi, vi, data = dat, method = "FE")   # fixed-effect instead

You can also pass a vector of methods (e.g. method = c("REML", "DL")) and metafor will try each in turn until one converges — useful when a small or unusual dataset makes REML fail to converge on its own. summary(res) or simply printing res returns the pooled estimate and its 95% CI, the test of the pooled effect (z-test), and the heterogeneity statistics: Cochran’s Q and its p-value, τ², and I² (the percentage of total variability attributable to between-study heterogeneity rather than sampling error). Read those heterogeneity numbers alongside CASRAI’s dedicated heterogeneity in meta-analysis guide before deciding whether pooling the studies at all was appropriate — a high I² is a reason to investigate subgroup or sensitivity analysis, not necessarily a reason to abandon the pooled estimate.

Choosing between the two model families is a substantive decision, not a technical default: a fixed-effect model assumes every included study is estimating the exact same true effect and differences are pure sampling error; a random-effects model assumes the true effect varies across studies (different populations, protocols, dosages) and estimates the pooled mean of that distribution of effects. See CASRAI’s fixed-effect vs. random-effects meta-analysis guide for how to make that call before you touch method.

Visualizing the pooled result with forest()

Once you have a fitted rma object, forest(res) produces a forest plot directly from it — one row per included study, a summary diamond for the pooled estimate at the bottom, and (by default) study labels and weights alongside the point estimates. Because it plots straight from the model object, the weights and pooled estimate shown are guaranteed to match what rma() actually fit — there is no separate step where the plot and the model can drift out of sync. Common customizations: slab to supply your own study labels (author/year) instead of row numbers, xlim/alim to control the plotted and displayed axis range separately, and addpoly() to layer a second summary polygon onto an existing forest plot — useful for showing a subgroup result alongside the overall pooled estimate on one plot. If you’re not yet comfortable reading the plot itself (what the square sizes, the horizontal lines, and the diamond represent), start with CASRAI’s how to read a forest plot guide and the forest plot dictionary entry — this section assumes that vocabulary and focuses only on producing the plot in R.

Checking for small-study effects with funnel()

funnel(res) produces a funnel plot — each study’s effect estimate plotted against a measure of its precision (standard error by default, though sample size and other precision measures can be substituted) — from the same fitted rma object, no re-specification needed. Asymmetry in the funnel — smaller (less precise) studies scattered off to one side rather than spread symmetrically around the pooled estimate — is the classic visual signal of possible publication bias or other small-study effects, though metafor’s own documentation is explicit that asymmetry has multiple possible causes beyond publication bias (genuine heterogeneity, methodological differences correlated with study size), so don’t over-read the plot alone. Two functions extend the visual check into a formal test:

• regtest(res) runs Egger’s regression test for funnel-plot asymmetry — regressing the standardized effect against precision and testing whether the intercept differs from zero. Like Egger’s test generally, it has low power with few included studies.
• trimfill(res) runs the trim-and-fill procedure, a rank-based method that estimates how many studies may be “missing” from one side of the funnel and imputes them to produce an adjusted pooled estimate. Calling funnel() on the object trimfill() returns plots the observed and imputed points together on one funnel plot.

Treat all three (visual funnel plot, regtest(), and trimfill()) as complementary evidence, not a single pass/fail gate — the same caution CASRAI’s funnel plot and publication bias entries give for interpreting asymmetry by eye.

Common mistakes when running this in R for the first time

• Mixing effect-size definitions across studies. rma() will happily pool a column of yi values even if half were computed as "SMD" and half as raw "MD" by a co-author working from a different extraction sheet — it has no way to detect this. Confirm every row’s effect size was computed with the same measure before pooling.
• Confusing variance and standard error. vi is a variance, sei is a standard error (its square root) — passing a standard-error column into the vi slot silently produces wrong weights and a wrong CI, with no error thrown, since both are just numeric vectors to rma().
• Reporting I² as if it were the same thing as τ². I² is a percentage describing the proportion of total variance due to heterogeneity; τ² is the actual estimated between-study variance on the scale of the effect size. Both are printed by rma()‘s summary output but they answer different questions.
• Treating a non-significant Q-test as proof of no heterogeneity. Cochran’s Q has notoriously low power with a small number of studies — a non-significant Q-test with five included studies does not confirm homogeneity; look at I² and the confidence interval around τ² as well.

Where metafor fits among evidence-synthesis tools

metafor is a script-based, open-source option, which trades a GUI for full reproducibility — the entire analysis, from raw data to the final forest plot, is one auditable R script rather than a series of manual point-and-click steps in commercial software. Researchers running a full systematic review (not just the pooling step) typically pair it with review-management platforms that handle screening and data extraction upstream, before the cleaned effect sizes ever reach escalc() and rma().

Frequently asked questions

Is metafor free?

Yes — metafor is released under the GPL license and distributed free through CRAN, R’s standard package repository (install.packages("metafor")).

What does the “rma” in rma() stand for?

Random-effects (and, more generally, mixed-effects) meta-analysis — the function’s formal name is rma.uni(), with rma() as its shorthand alias, and it fits fixed-effect models too via method = "FE", not only random-effects ones.

Do I need escalc() every time?

No. escalc() is only needed when you’re computing effect sizes from raw study data (means/SDs/Ns or 2×2 counts). If your extraction table already has a consistent effect size and variance per study, hand those columns to rma() directly.

Can forest() and funnel() be customized for a journal figure?

Yes — both accept extensive graphical parameters (label text, axis limits, colors, point shapes) since they’re built on R’s base graphics system, and the resulting plot can be exported with any standard R device (png(), pdf(), svg()) for direct use in a manuscript.

For the underlying methodology this package implements — what a systematic review is, when meta-analysis is and isn’t appropriate within one, and the PRISMA reporting standard — see CASRAI’s research tools hub and its guides on choosing between fixed- and random-effects models and interpreting heterogeneity.