How to Make a Forest Plot in R with metafor (Step by Step)

Knowing how to read and interpret a forest plot is one thing; building your own so it comes out at publication quality is another. In this tutorial you will build a complete forest plot in R with the metafor package (the reference for meta-analysis in R) from raw means and standard deviations to a figure ready to submit to a Q1 journal. All the code is reproducible: copy it, swap in your data and it works.

What you need before you start

A forest plot summarizes a meta-analysis, so you need data from several studies on the same question. There are two common scenarios:

• You have the raw data: for each study, the means, standard deviations and sample sizes of the two groups (treatment and control). From these you will compute the effect size.
• You have effects already computed: for each study, the effect size (for example a Cohen's d or an odds ratio) and its variance or standard error. This usually comes from tables in prior articles.

In both cases, the workflow in metafor is the same: (1) obtain the effect and its variance with escalc(), (2) fit the model with rma() and (3) draw it with forest(). Let's go step by step.

Step 1: Install and load metafor

# First time only
install.packages("metafor")

library(metafor)

The metafor package (Viechtbauer, 2010) is by far the most widely used for meta-analysis in R and the one reviewers expect to see. It is actively maintained and covers everything from basic meta-analysis to meta-regression, multilevel models and robust-variance estimation.

Step 2: Prepare the data and compute the effect size

Suppose you have six studies comparing a psychological intervention against a control group on an anxiety measure. We build the data frame and compute the effect size with escalc():

mydata <- data.frame(
  study   = c("Garcia 2019", "Lopez 2020", "Martinez 2021",
              "Chen 2021", "Fernandez 2022", "Sanchez 2023"),
  m_trt   = c(22.1, 19.8, 24.5, 21.0, 26.2, 20.4),
  sd_trt  = c(5.2, 4.8, 6.1, 5.5, 5.9, 4.6),
  n_trt   = c(40, 65, 32, 58, 45, 70),
  m_ctrl  = c(24.8, 21.5, 27.0, 22.3, 29.1, 22.0),
  sd_ctrl = c(5.4, 5.0, 6.3, 5.2, 6.0, 4.9),
  n_ctrl  = c(38, 63, 30, 55, 44, 68)
)

# Standardized mean difference (Hedges' g, corrected for small samples)
dat <- escalc(measure = "SMD",
              m1i = m_trt, sd1i = sd_trt, n1i = n_trt,
              m2i = m_ctrl, sd2i = sd_ctrl, n2i = n_ctrl,
              data = mydata, slab = study)

After this, dat has two new columns: yi (the effect size, here Hedges' g) and vi (its sampling variance). The argument measure = "SMD" returns Hedges' g, which applies the small-sample bias correction to Cohen's d, the standard in psychology. For dichotomous data, use measure = "OR" (odds ratio) or "RR" (risk ratio) with the frequency arguments (ai, bi, ci, di); and if you already have the computed effects, skip this step and pass your yi and vi columns straight to the model.

The slab (study label) argument defines the per-study labels that appear on the left of the plot. Set it here and you won't have to repeat it later.

Step 3: Fit the random-effects model

res <- rma(yi, vi, data = dat, method = "REML")
res

The rma() function fits a random-effects model by default (REML estimator), which is the standard in psychology and health sciences: it assumes each study estimates a different true effect drawn from a distribution and incorporates the between-study variance (τ²) into the calculation. If for some justified reason you need fixed effects, use method = "FE", but except in very specific cases, keep REML.

Printing res gives you the pooled effect, its confidence interval, and the heterogeneity indicators: I², τ² and the Q test. You should understand what they mean before drawing anything; if I² is high, the forest-plot diamond is a poorly representative average. The guide on heterogeneity in meta-analysis (I²) explains how to interpret them and when the pooled effect stops being informative.

Step 4: Draw the basic forest plot

With the model fitted, a single line produces the plot:

forest(res)

This already produces a working forest plot: one row per study with its square proportional to the weight and its confidence interval, the pooled-effect diamond at the bottom and the axis with the effect scale. But for publication it is worth refining.

Step 5: Customize for publication quality

These arguments cover 95% of what you need for the figure to look like one in a Q1 article:

forest(res,
       header = c("Study", "g [95% CI]"),    # column headers
       xlab   = "Hedges' g",                  # axis label
       mlab   = "Random-effects model (REML)",# diamond label
       slab   = dat$study,                    # study labels
       order  = "obs",                        # order by effect size
       cex    = 0.9)                           # font size

Some additional useful arguments: xlim and alim control the limits of the drawing area and of the axis; at sets the axis ticks (e.g. at = c(-1, -0.5, 0, 0.5)); refline moves the line of no effect (default 0 for differences, use refline = 1 for odds ratios); and ilab lets you add extra columns such as sample sizes. For odds ratios or risk ratios, add atransf = exp so the axis shows the original scale instead of the log scale.

Step 6: Add the heterogeneity statistics to the plot

A publishable forest plot shows I², Q and τ², usually below the diamond. You can include them automatically with bquote and mlab:

forest(res, header = TRUE, xlab = "Hedges' g",
       mlab = bquote(paste("RE Model (Q = ",
              .(formatC(res$QE, digits = 2, format = "f")),
              ", df = ", .(res$k - res$p),
              ", p = ", .(formatC(res$QEp, digits = 3, format = "f")),
              "; ", I^2, " = ", .(formatC(res$I2, digits = 1, format = "f")), "%)")))

This way the reader sees at a glance whether heterogeneity undermines the diamond, without having to hunt for it in the text.

Forest plot with subgroups

If a moderator variable (for example, intervention format: in-person vs. online) might explain the heterogeneity, a subgroup forest plot is worth it. The cleanest way is to fit the model with the moderator and draw the studies grouped, adding a partial diamond per subgroup with addpoly():

# Model with categorical moderator
res.mod <- rma(yi, vi, mods = ~ format, data = dat)

# Submodels per level
res.in  <- rma(yi, vi, data = dat, subset = (format == "in-person"))
res.onl <- rma(yi, vi, data = dat, subset = (format == "online"))

# Draw ordered by subgroup and add partial diamonds
forest(res, order = order(dat$format), addfit = FALSE,
       header = TRUE, xlab = "Hedges' g")
addpoly(res.in,  row = 1, mlab = "Subtotal: in-person")
addpoly(res.onl, row = 5, mlab = "Subtotal: online")

Adjust the row values to the number of studies in each subgroup. The between-subgroups test (Q_between) comes from res.mod: if it is significant, the moderator explains part of the variation between studies.

Export the plot at publication quality

Journals require images at 300 dpi or higher, and many demand vector format (PDF, EPS, TIFF). Open a graphics device, draw and close it:

# PNG at 300 dpi
png("forest_plot.png", width = 8, height = 6, units = "in", res = 300)
forest(res, header = TRUE, xlab = "Hedges' g")
dev.off()

# Vector PDF (recommended by many journals)
pdf("forest_plot.pdf", width = 8, height = 6)
forest(res, header = TRUE, xlab = "Hedges' g")
dev.off()

Vector format (PDF/EPS) does not lose quality when scaled and is usually the production editor's preferred option. Adjust width and height until the labels don't overlap.

Four common mistakes when building the forest plot

1. Passing the standard deviation instead of the variance to vi. The second argument of rma() is the sampling variance of the effect, not the standard error nor the standard deviation. If you have the standard error, use sei instead (rma(yi, sei = my_se, ...)). Confusing them artificially inflates or shrinks the precision and distorts the weights.

2. Skipping escalc() and feeding raw means. A forest plot needs effects on a common scale. If you compare studies with means on different scales without standardizing, the plot is meaningless. escalc() does that conversion and also corrects the small-sample bias in Hedges' g.

3. Forgetting the sign of the effect. Decide from the start which direction is "improvement" and check that all studies are coded the same way. A study with the subtraction reversed (control − treatment instead of treatment − control) will appear on the wrong side and bias the pooled effect.

4. Using fixed effects without justifying it. Although rma() uses random effects by default, sometimes borrowed code carries method = "FE". Unless you have a solid argument (very few functionally identical studies), keep random effects: it is what reviewers expect in psychology and health.

Alternative: the meta package

If you prefer slightly more direct syntax, the meta package also produces publishable forest plots. metacont() for continuous data or metabin() for binary data compute the effect and the model in a single function, and forest() draws the plot with a slightly different style. Both metafor and meta are valid options accepted in indexed journals; metafor offers more control and flexibility for advanced analyses (multilevel, meta-regression, robust models), so it is the one I recommend if you plan to grow in complexity.

And once the plot is ready

Building the forest plot is the second-to-last step. Before closing the analysis, check for publication bias with a funnel plot and the Egger's test, and make sure you can interpret every element of the forest plot before writing the discussion. For the final write-up, the guide on how to report results in APA 7 covers the format journals expect.