Meta-Analysis in R with metafor
Sunghyun Chang
07 October 2026
Introduction
This code through explores meta-analysis in R using the metafor package. Meta-analysis is a statistical method that combines results from multiple independent studies on the same question to produce a single, more precise estimate of an effect.
Content Overview
Specifically, we’ll explain and demonstrate how to run a basic meta-analysis using dat.bcg, a built-in dataset of 13 clinical trials on the effectiveness of the BCG vaccine against tuberculosis. We will use four core functions: escalc(), rma(), forest(), and funnel().
We will then move on to advanced techniques: testing for publication bias with regtest() and trimfill(), explaining heterogeneity through meta-regression with rma(mods = ~ ablat) and regplot(), and checking the robustness of results with leave1out(), influence(), and baujat().
Why You Should Care
This topic is valuable because individual studies often disagree or are too small to reach firm conclusions. Meta-analysis helps decision-makers see what the overall body of research says, rather than relying on a single study.
Learning Objectives
Specifically, you’ll learn how to:
Basics:
Calculate effect sizes from raw study data with escalc()
Fit a random-effects model and interpret its key results with rma()
Visualize study results with a forest plot using forest()
Check for potential publication bias with a funnel plot using funnel()
Advanced:
Test and adjust for publication bias with regtest() and trimfill()
Explain heterogeneity with meta-regression using rma(mods = ) and regplot()
Assess the robustness of results with leave1out(), influence(), and baujat()
Running a Meta-Analysis with metafor
Here, we’ll show the standard workflow of a meta-analysis in six steps:
- calculate an effect size for each study
- combine them with a statistical model
- visualize the results
- check for publication bias
- Explain heterogeneity across studies
- Check the robustness of the results
Brief Theoretical Background
This is based on the work of Viechtbauer (2010), who developed metafor, and on the random-effects model of DerSimonian and Laird (1986).
Effect size: Each study reports results differently, so we first convert them into a common metric. Here we use the log risk ratio, which compares the risk of tuberculosis in vaccinated versus unvaccinated groups.
Random-effects model: This model assumes that the true effect may differ across studies (e.g., due to different populations or settings), and estimates both the average effect and how much the effects vary.
Publication Bias: Studies with significant results are more likely to be published, which can inflate the pooled effect (Egger et al., 1997).
Heterogeneity Analysis: When true effects vary across studies, meta-regression tests whether study characteristics can explain these differences (Thompson & Higgins, 2002).
Sensitivity Analysis: This checks whether the overall result depends on a few influential or outlying studies (Viechtbauer & Cheung, 2010).
Basic Example (1): escalc()
escalc() shows how to calculate effect sizes. Each row of dat.bcg contains the number of TB cases (tpos, cpos) and non-cases (tneg, cneg) in the treated and control groups.
dat <- escalc(measure = "RR",
ai = tpos, bi = tneg,
ci = cpos, di = cneg,
data = dat.bcg)
head(dat[, c("author", "year", "tpos", "tneg", "cpos", "cneg", "yi", "vi")])yi is the log risk ratio for each study, and vi is its sampling variance. For example, Aronson (1948) has yi = −0.89, which corresponds to a risk ratio of exp(−0.89) ≈ 0.41: vaccinated participants had about 59% lower risk of TB.
Negative yi values mean fewer TB cases in the vaccinated group.
Basic Example (2): rma()
rma() can be used for estimating the overall effect across all studies.
##
## Random-Effects Model (k = 13; tau^2 estimator: REML)
##
## tau^2 (estimated amount of total heterogeneity): 0.3132 (SE = 0.1664)
## tau (square root of estimated tau^2 value): 0.5597
## I^2 (total heterogeneity / total variability): 92.22%
## H^2 (total variability / sampling variability): 12.86
##
## Test for Heterogeneity:
## Q(df = 12) = 152.2330, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## -0.7145 0.1798 -3.9744 <.0001 -1.0669 -0.3622 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1##
## pred ci.lb ci.ub pi.lb pi.ub
## 0.4894 0.3441 0.6962 0.1546 1.5490The pooled log risk ratio is −0.71 (95% CI: −1.07 to −0.36, p < .0001). Converted with predict(), the pooled risk ratio is 0.49 (95% CI: 0.34 to 0.70), meaning vaccinated individuals had roughly half the risk of TB.
However, heterogeneity is very high: I² = 92.2% and Q(12) = 152.23 (p < .0001). The 95% prediction interval (0.15 to 1.55) includes 1, which suggests the vaccine may have little or no effect in some settings.
Basic Example (3): forest()
forest() can also be used for visualizing every study at once.
Each square is one study’s estimate (larger squares mean more weight), lines are 95% confidence intervals, and the diamond at the bottom is the pooled estimate. atransf = exp displays risk ratios instead of log values.
Most studies fall to the left of the vertical line (RR = 1), favoring the vaccine. Notably, the TPT Madras (1980) trial, one of the most heavily weighted studies (largest square), shows no effect (RR = 1.01), while Hart & Sutherland (1977) shows a strong effect (RR = 0.24). This wide spread is what drives the high I².
Basic Example (4): funnel()
funnel() is valuable for checking publication bias.
If studies are spread symmetrically around the pooled estimate, publication bias is less likely. An asymmetric funnel may suggest that small studies with null results were never published.
Since visual judgment can be subjective, we test this formally in the next section.
Advanced Functions: Publication Bias
Advanced Example (1): regtest()
While a funnel plot lets us see asymmetry, regtest() tests it statistically using Egger’s regression test. It checks whether smaller studies (with larger standard errors) tend to report systematically different effects.
##
## Regression Test for Funnel Plot Asymmetry
##
## Model: mixed-effects meta-regression model
## Predictor: standard error
##
## Test for Funnel Plot Asymmetry: z = -0.8033, p = 0.4218
## Limit Estimate (as sei -> 0): b = -0.5104 (CI: -1.1182, 0.0974)The test is not significant (z = −0.80, p = 0.42). A non-significant result (p > .05) means there is no statistical evidence of funnel plot asymmetry.
The limit estimate (b = −0.51), the expected effect in a hypothetical infinitely large study, is still negative, which supports a protective effect.
Advanced Example (2): trimfill()
The trim-and-fill method estimates how many studies might be “missing” from the funnel plot because of publication bias, imputes them, and recalculates the pooled effect.
##
## Estimated number of missing studies on the right side: 1 (SE = 2.4528)
##
## Random-Effects Model (k = 14; tau^2 estimator: REML)
##
## tau^2 (estimated amount of total heterogeneity): 0.3313 (SE = 0.1701)
## tau (square root of estimated tau^2 value): 0.5756
## I^2 (total heterogeneity / total variability): 92.14%
## H^2 (total variability / sampling variability): 12.72
##
## Test for Heterogeneity:
## Q(df = 13) = 154.6750, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## -0.6571 0.1785 -3.6805 0.0002 -1.0070 -0.3072 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
In the funnel plot, imputed studies appear as open (white) circles. If the adjusted estimate is close to the original, our conclusion is robust to publication bias.
Trim-and-fill estimates one missing study on the right side of the funnel. After adding it, the pooled log risk ratio changes from −0.71 to −0.66, which is RR = 0.52 (95% CI: 0.37 to 0.74), compared with 0.49 originally. Because the adjusted estimate remains significant (p = .0002), our conclusion is robust to potential publication bias.
Advanced Functions: Explaining Heterogeneity
Advanced Example (3): rma(…, mods = ~ ablat)
Our model showed very high heterogeneity (I² ≈ 92%). Meta-regression asks why: it adds a study-level variable (a moderator) to see whether it explains differences between studies. Here we use ablat, the absolute latitude of each study location, since climate and environment may affect how well the vaccine works.
##
## Mixed-Effects Model (k = 13; tau^2 estimator: REML)
##
## tau^2 (estimated amount of residual heterogeneity): 0.0764 (SE = 0.0591)
## tau (square root of estimated tau^2 value): 0.2763
## I^2 (residual heterogeneity / unaccounted variability): 68.39%
## H^2 (unaccounted variability / sampling variability): 3.16
## R^2 (amount of heterogeneity accounted for): 75.62%
##
## Test for Residual Heterogeneity:
## QE(df = 11) = 30.7331, p-val = 0.0012
##
## Test of Moderators (coefficient 2):
## QM(df = 1) = 16.3571, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## intrcpt 0.2515 0.2491 1.0095 0.3127 -0.2368 0.7397
## ablat -0.0291 0.0072 -4.0444 <.0001 -0.0432 -0.0150 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1The negative coefficient for ablat means the vaccine is more effective farther from the equator. The R² value shows how much of the heterogeneity is explained by latitude.
Latitude is a significant moderator (QM(1) = 16.36, p < .0001). The coefficient of −0.029 means that each additional degree of latitude lowers the log risk ratio by 0.029, so the vaccine is more effective farther from the equator.
Latitude explains 75.6% of the heterogeneity (R²): τ² drops from 0.31 to 0.08, and I² from 92.2% to 68.4%. Some heterogeneity remains (QE(11) = 30.73, p = .001), however, so other factors also matter.
Advanced Example (4): regplot()
regplot() draws a bubble plot of the meta-regression. Each bubble is a study, with larger bubbles for more precise studies, and the line shows the predicted effect across latitudes.
As latitude increases, the risk ratio drops further below 1, which means stronger protection.
The model predicts a risk ratio of about 0.88 at 13° latitude (near the equator) but about 0.26 at 55°, a dramatic difference in vaccine effectiveness.
Advanced Functions: Sensitivity Analysis
Advanced Example (5): leave1out()
A leave-one-out analysis re-runs the meta-analysis 13 times, each time removing one study. It shows whether any single study is driving the overall result.
##
## estimate zval pval ci.lb ci.ub Q
## -Aronson, 1948 0.4931 -3.7223 0.0002 0.3398 0.7155 151.5826
## -Ferguson & Simes, 1949 0.5199 -3.6195 0.0003 0.3649 0.7409 145.3176
## -Rosenthal et al, 1960 0.5038 -3.6916 0.0002 0.3501 0.7250 150.1970
## -Hart & Sutherland, 1977 0.5334 -3.5580 0.0004 0.3774 0.7541 96.5626
## -Frimodt-Moller et al, 1973 0.4657 -3.9845 0.0001 0.3198 0.6782 151.3200
## -Stein & Aronson, 1953 0.4912 -3.5499 0.0004 0.3318 0.7273 128.1867
## -Vandiviere et al, 1973 0.5193 -3.6307 0.0003 0.3646 0.7397 145.8296
## -TPT Madras, 1980 0.4517 -4.4184 0.0000 0.3175 0.6426 67.9858
## -Coetzee & Berjak, 1968 0.4765 -3.7686 0.0002 0.3241 0.7007 152.2051
## -Rosenthal et al, 1961 0.5205 -3.5439 0.0004 0.3627 0.7469 139.8271
## -Comstock et al, 1974 0.4687 -3.8708 0.0001 0.3193 0.6879 151.4655
## -Comstock & Webster, 1969 0.4678 -4.1727 0.0000 0.3274 0.6684 150.7868
## -Comstock et al, 1976 0.4595 -4.1908 0.0000 0.3194 0.6611 149.7884
## Qp tau2 I2 H2
## -Aronson, 1948 0.0000 0.3362 93.2259 14.7622
## -Ferguson & Simes, 1949 0.0000 0.2926 92.2540 12.9098
## -Rosenthal et al, 1960 0.0000 0.3207 92.9354 14.1551
## -Hart & Sutherland, 1977 0.0000 0.2628 90.4125 10.4302
## -Frimodt-Moller et al, 1973 0.0000 0.3278 92.7634 13.8187
## -Stein & Aronson, 1953 0.0000 0.3596 90.9118 11.0033
## -Vandiviere et al, 1973 0.0000 0.2930 92.2777 12.9495
## -TPT Madras, 1980 0.0000 0.2732 87.0314 7.7109
## -Coetzee & Berjak, 1968 0.0000 0.3495 93.2133 14.7346
## -Rosenthal et al, 1961 0.0000 0.2987 92.2322 12.8737
## -Comstock et al, 1974 0.0000 0.3405 91.8110 12.2114
## -Comstock & Webster, 1969 0.0000 0.3082 92.6782 13.6579
## -Comstock et al, 1976 0.0000 0.3037 92.3444 13.0623If the pooled risk ratio stays similar (and below 1) no matter which study is removed, the result is robust.
Removing any single study yields a pooled risk ratio between 0.45 (without TPT Madras, 1980) and 0.53 (without Hart & Sutherland, 1977). Every estimate remains significant with a confidence interval below 1, so no single study drives the overall conclusion. Removing TPT Madras also lowers I² the most (from 92.2% to 87.0%), which indicates it is a major source of heterogeneity.
Advanced Example (6): influence()
influence() calculates diagnostics such as Cook’s distance and standardized residuals to flag studies that are outliers or have unusually large influence on the results.
Studies exceeding the dotted threshold lines would be flagged as influential (highlighted in red).
No study exceeds the thresholds (dotted lines) for studentized residuals, DFFITS, or hat values, so none is formally flagged as influential. Study 4 (Hart & Sutherland, 1977) has the largest Cook’s distance, and removing study 8 (TPT Madras, 1980) produces the largest drop in Q (QE.del).
Advanced Example (7): baujat()
A Baujat plot shows each study’s contribution to overall heterogeneity (x-axis) against its influence on the pooled result (y-axis). Studies in the upper-right corner are both heterogeneous and influential.
This helps identify which specific studies are responsible for the high I² we observed earlier.
Studies 4 (Hart & Sutherland, 1977) and 8 (TPT Madras, 1980) appear in the upper-right corner, so they contribute the most to heterogeneity and also have the greatest influence on the pooled result. This matches what we found in leave1out() and influence().
Further Resources
Learn more about meta-analysis and the metafor package
with the following:
Package Documentation
The metafor Package Website – Official site by the package author, with tutorials, worked examples, and an FAQ.
metafor on CRAN – Reference manual listing every function and its arguments.
Tutorials
- Doing Meta-Analysis with R: A Hands-On Guide – A free online textbook that walks through each step of a meta-analysis in R.
Methodological Guidance
- Cochrane Handbook for Systematic Reviews of Interventions – The standard reference for conducting systematic reviews and meta-analyses, including how to assess heterogeneity and publication bias.
Works Cited
This code through references and cites the following sources:
Anthropic. (2026). Claude [Large language model]. https://claude.ai – Used to assist with structuring this tutorial and writing example code; all code was run and verified by the author.
Colditz, G. A., Brewer, T. F., Berkey, C. S., Wilson, M. E., Burdick, E., Fineberg, H. V., & Mosteller, F. (1994). Efficacy of BCG vaccine in the prevention of tuberculosis: Meta-analysis of the published literature. JAMA, 271(9), 698–702. https://doi.org/10.1001/jama.1994.03510330076038
DerSimonian, R., & Laird, N. (1986). Meta-analysis in clinical trials. Controlled Clinical Trials, 7(3), 177–188. https://doi.org/10.1016/0197-2456(86)90046-2
Egger, M., Davey Smith, G., Schneider, M., & Minder, C. (1997). Bias in meta-analysis detected by a simple, graphical test. BMJ, 315(7109), 629–634. https://doi.org/10.1136/bmj.315.7109.629
Thompson, S. G., & Higgins, J. P. T. (2002). How should meta-regression analyses be undertaken and interpreted? Statistics in Medicine, 21(11), 1559–1573. https://doi.org/10.1002/sim.1187
Viechtbauer, W. (2010). Conducting meta-analyses in R with the metafor package. Journal of Statistical Software, 36(3), 1–48. https://doi.org/10.18637/jss.v036.i03
Viechtbauer, W., & Cheung, M. W.-L. (2010). Outlier and influence diagnostics for meta-analysis. Research Synthesis Methods, 1(2), 112–125. https://doi.org/10.1002/jrsm.11