Package 'metafor'

metafor: A Meta-Analysis Package for R

Description

The metafor package provides a comprehensive collection of functions for conducting meta-analyses in R. The package can be used to calculate various effect sizes or outcome measures and then allows the user to fit equal-, fixed-, and random-effects models to these data. By including study-level variables (‘moderators’) as predictors in these models, (mixed-effects) meta-regression models can also be fitted. For meta-analyses of \(2 \times 2\) tables, proportions, incidence rates, and incidence rate ratios, the package also provides functions that implement specialized methods, including the Mantel-Haenszel method, Peto's method, and a variety of suitable generalized linear mixed-effects models (i.e., mixed-effects logistic and Poisson regression models). For non-independent effects/outcomes (e.g., due to correlated sampling errors, correlated true effects or outcomes, or other forms of clustering), the package also provides a function for fitting multilevel and multivariate models to meta-analytic data.

Various methods are available to assess model fit, to identify outliers and/or influential studies, and for conducting sensitivity analyses (e.g., standardized residuals, Cook's distances, leave-one-out analyses). Advanced techniques for hypothesis testing and obtaining confidence intervals (e.g., for the average effect or outcome or for the model coefficients in a meta-regression model) have also been implemented (e.g., the Knapp and Hartung method, permutation tests, cluster-robust inference methods / robust variance estimation).

The package also provides functions for creating forest, funnel, radial (Galbraith), normal quantile-quantile, L'Abbé, Baujat, bubble, and GOSH plots. The presence of publication bias (or more precisely, funnel plot asymmetry or ‘small-study effects’) and its potential impact on the results can be examined via the rank correlation and Egger's regression test, the trim and fill method, the test of excess significance, and by applying a variety of selection models.

The escalc Function

[escalc] Before a meta-analysis can be conducted, the relevant results from each study must be quantified in such a way that the resulting values can be further aggregated and compared. The escalc function can be used to compute a wide variety of effect sizes or ‘outcome measures’ (and the corresponding sampling variances) that are often used in meta-analyses (e.g., risk ratios, odds ratios, risk differences, mean differences, standardized mean differences, response ratios / ratios of means, raw or r-to-z transformed correlation coefficients). Measures for quantifying some characteristic of individual groups (e.g., in terms of means, proportions, or incidence rates and transformations thereof), measures of change (e.g., raw and standardized mean changes), and measures of variability (e.g., variability ratios and coefficient of variation ratios) are also available.

The rma.uni Function

[rma.uni] The various meta-analytic models that are typically used in practice are special cases of the general linear (mixed-effects) model. The rma.uni function (with alias rma) provides a general framework for fitting such models. The function can be used in combination with any of the effect sizes or outcome measures computed with the escalc function or, more generally, any set of estimates (with corresponding sampling variances or standard errors) one would like to analyze. The notation and models underlying the rma.uni function are explained below.

For a set of \(i = 1, \ldots, k\) independent studies, let \(y_i\) denote the observed value of the effect size or outcome measure in the \(i\textrm{th}\) study. Let \(\theta_i\) denote the corresponding (unknown) true effect/outcome, such that \[y_i \mid \theta_i \sim N(\theta_i, v_i).\] In other words, the observed effect sizes or outcomes are assumed to be unbiased and normally distributed estimates of the corresponding true effects/outcomes with sampling variances equal to \(v_i\) (where \(v_i\) is just the square of the standard errors of the estimates). The \(v_i\) values are assumed to be known. Depending on the outcome measure used, a bias correction, normalizing, and/or variance stabilizing transformation may be necessary to ensure that these assumptions are (at least approximately) true (e.g., the log transformation for odds/risk ratios, the bias correction for standardized mean differences, Fisher's r-to-z transformation for correlations; see escalc for more details).

According to the random-effects model, we further assume that \(\theta_i \sim N(\mu, \tau^2)\), that is, the true effects/outcomes are normally distributed with \(\mu\) denoting the average true effect/outcome and \(\tau^2\) the variance in the true effects/outcomes (\(\tau^2\) is therefore often referred to as the amount of ‘heterogeneity’ in the true effects/outcomes or the ‘between-study variance’). The random-effects model can also be written as \[y_i = \mu + u_i + \varepsilon_i,\] where \(u_i \sim N(0, \tau^2)\) and \(\varepsilon_i \sim N(0, v_i)\). The fitted model provides estimates of \(\mu\) and \(\tau^2\), that is, \[\hat{\mu} = \frac{\sum_{i=1}^k w_i y_i}{\sum_{i=1}^k w_i},\] where \(w_i = 1/(\hat{\tau}^2 + v_i)\) and \(\hat{\tau}^2\) denotes an estimate of \(\tau^2\) obtained with one of the many estimators that have described in the literature for this purpose (this is the standard ‘inverse-variance’ method for random-effects models).

A special case of the model above is the equal-effects model (also sometimes called the common-effect(s) model) which arises when \(\tau^2 = 0\). In this case, the true effects/outcomes are homogeneous (i.e., \(\theta_1 = \theta_2 = \ldots = \theta_k \equiv \theta\)) and hence we can write the model as \[y_i = \theta + \varepsilon_i,\] where \(\theta\) denotes the true effect/outcome in the studies, which is estimated with \[\hat{\theta} = \frac{\sum_{i=1}^k w_i y_i}{\sum_{i=1}^k w_i},\] where \(w_i = 1/v_i\) (again, this is the standard ‘inverse-variance’ method as described in the meta-analytic literature). Note that the commonly-used term ‘fixed-effects model’ is not used here - for an explanation, see here.

Study-level variables (often referred to as ‘moderators’) can also be included as predictors in meta-analytic models, leading to so-called ‘meta-regression’ analyses (to examine whether the effects/outcomes tend to be larger/smaller under certain conditions or circumstances). When including moderator variables in a random-effects model, we obtain a mixed-effects meta-regression model. This model can be written as \[y_i = \beta_0 + \beta_1 x_{i1} + \beta_2 x_{i2} + \ldots + \beta_{p'} x_{ip'} + u_i + \varepsilon_i,\] where \(u_i \sim N(0, \tau^2)\) and \(\varepsilon_i \sim N(0, v_i)\) as before and \(x_{ij}\) denotes the value of the \(j\textrm{th}\) moderator variable for the \(i\textrm{th}\) study (letting \(p = p' + 1\) denote the total number of coefficients in the model including the model intercept). Therefore, \(\beta_j\) denotes how much the average true effect/outcome differs for studies that differ by one unit on \(x_{ij}\) and the model intercept \(\beta_0\) denotes the average true effect/outcome when the values of all moderator variables are equal to zero. The value of \(\tau^2\) in the mixed-effects model denotes the amount of ‘residual heterogeneity’ in the true effects/outcomes (i.e., the amount of variability in the true effects/outcomes that is not accounted for by the moderators included in the model). In matrix notation, the model can also be written as \[y = X\beta + u + \varepsilon,\] where \(y\) is a column vector with the observed effect sizes or outcomes, \(X\) is the \(k \times p\) model matrix (with the first column equal to 1s for the intercept term), \(\beta\) is a column vector with the model coefficients, and \(u\) and \(\varepsilon\) are column vectors for the random effects and sampling errors.

The rma.mh Function

[rma.mh] The Mantel-Haenszel method provides an alternative approach for fitting equal-effects models when dealing with studies providing data in the form of \(2 \times 2\) tables or in the form of event counts (i.e., person-time data) for two groups (Mantel & Haenszel, 1959). The method is particularly advantageous when aggregating a large number of studies with small sample sizes (the so-called sparse data or increasing strata case). The Mantel-Haenszel method is implemented in the rma.mh function. It can be used in combination with risk ratios, odds ratios, risk differences, incidence rate ratios, and incidence rate differences.

The rma.peto Function

[rma.peto] Yet another method that can be used in the context of a meta-analysis of \(2 \times 2\) table data is Peto's method (see Yusuf et al., 1985), implemented in the rma.peto function. The method provides an estimate of the (log) odds ratio under an equal-effects model. The method is particularly advantageous when the event of interest is rare, but see the documentation of the function for some caveats.

The rma.glmm Function

[rma.glmm] Dichotomous response variables and event counts (based on which one can calculate outcome measures such as odds ratios, incidence rate ratios, proportions, and incidence rates) are often assumed to arise from binomial and Poisson distributed data. Meta-analytic models that are directly based on such distributions are implemented in the rma.glmm function. These models are essentially special cases of generalized linear mixed-effects models (i.e., mixed-effects logistic and Poisson regression models). For \(2 \times 2\) table data, a mixed-effects conditional logistic model (based on the non-central hypergeometric distribution) is also available. Random/mixed-effects models with dichotomous data are often referred to as ‘binomial-normal’ models in the meta-analytic literature. Analogously, for event count data, such models could be referred to as ‘Poisson-normal’ models.

The rma.mv Function

[rma.mv] Standard meta-analytic models assume independence between the observed effect sizes or outcomes obtained from a set of studies. This assumption is often violated in practice. Dependencies can arise for a variety of reasons. For example, the sampling errors and/or true effects/outcomes may be correlated in multiple treatment studies (e.g., when multiple treatment groups are compared with a common control/reference group, such that the data from the control/reference group is used multiple times to compute the observed effect sizes or outcomes) or in multiple endpoint studies (e.g., when more than one effect size estimate or outcome is calculated based on the same sample of subjects due to the use of multiple endpoints or response variables). Correlations in the true effects/outcomes can also arise due to other forms of clustering (e.g., when multiple effects/outcomes derived from the same author, lab, or research group may be more similar to each other than effects/outcomes derived from different authors, labs, or research groups). In ecology and related fields, the shared phylogenetic history among the organisms studied (e.g., plants, fungi, animals) can also induce correlations among the effects/outcomes. The rma.mv function can be used to fit suitable meta-analytic multivariate/multilevel models to such data, so that the non-independence in the effects/outcomes is accounted for. Network meta-analyses (also called multiple/mixed treatment comparisons) can also be carried out with this function.

Future Plans and Updates

The metafor package is a work in progress and is updated on a regular basis with new functions and options. With metafor.news(), you can read the ‘NEWS’ file of the package after installation. Comments, feedback, and suggestions for improvements are always welcome.

Citing the Package

To cite the package, please use the following reference:

Viechtbauer, W. (2010). Conducting meta-analyses in R with the metafor package. Journal of Statistical Software, 36(3), 1-48. doi:10.18637/jss.v036.i03

Getting Started with the Package

The paper mentioned above is a good starting place for those interested in using the package. The purpose of the article is to provide a general overview of the package and its capabilities (as of version 1.4-0). Not all of the functions and options are described in the paper, but it should provide a useful introduction to the package. The paper can be freely downloaded from the URL given above or can be directly loaded with the command vignette("metafor").

In addition to reading the paper, carefully read this page and then the help pages for the escalc and the rma.uni functions (or the rma.mh, rma.peto, rma.glmm, and/or rma.mv functions if you intend to use these methods). The help pages for these functions provide links to many additional functions, which can be used after fitting a model. You can also read the entire documentation online at https://wviechtb.github.io/metafor/ (where it is nicely formatted and the output from all examples is provided).

A (pdf) diagram showing the various functions in the metafor package (and how they are related to each other) can be opened with the command vignette("diagram", package="metafor").

Finally, additional information about the package, several detailed analysis examples, examples of plots and figures provided by the package (with the corresponding code), some additional tips and notes, and a FAQ can be found on the package website at https://www.metafor-project.org.

Author(s)

Wolfgang Viechtbauer [email protected]
package website: https://www.metafor-project.org
author homepage: ⁠https://www.wvbauer.com⁠

Suggestions on how to obtain help with using the package can found on the package website at: https://www.metafor-project.org/doku.php/help

References

Cooper, H., Hedges, L. V., & Valentine, J. C. (Eds.) (2009). The handbook of research synthesis and meta-analysis (2nd ed.). New York: Russell Sage Foundation.

Hedges, L. V., & Olkin, I. (1985). Statistical methods for meta-analysis. San Diego, CA: Academic Press.

Mantel, N., & Haenszel, W. (1959). Statistical aspects of the analysis of data from retrospective studies of disease. Journal of the National Cancer Institute, 22(4), 719–748. ⁠https://doi.org/10.1093/jnci/22.4.719⁠

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⁠

Yusuf, S., Peto, R., Lewis, J., Collins, R., & Sleight, P. (1985). Beta blockade during and after myocardial infarction: An overview of the randomized trials. Progress in Cardiovascular Disease, 27(5), 335–371. ⁠https://doi.org/10.1016/s0033-0620(85)80003-7⁠

---

Add Polygons to Forest Plots

Description

Function to add polygons (sometimes called ‘diamonds’) to a forest plot, for example to show summary estimates for subgroups of studies or to show fitted/predicted values based on models involving moderators.

Usage

addpoly(x, ...)

Arguments

x	either an object of class "rma", an object of class "predict.rma", or the values at which polygons should be drawn. See ‘Details’.
...	other arguments.

Details

Currently, methods exist for three types of situations.

In the first case, object x is a fitted model coming from the rma.uni, rma.mh, rma.peto, rma.glmm, or rma.mv functions. The model must either be an equal- or a random-effects model, that is, the model should not contain any moderators. The corresponding method is addpoly.rma. It can be used to add a polygon to an existing forest plot (usually at the bottom), showing the summary estimate (with its confidence interval) based on the fitted model.

Alternatively, x can be an object of class "predict.rma" obtained with the predict function. In this case, polygons based on the predicted values are drawn. The corresponding method is addpoly.predict.rma.

Alternatively, object x can be a vector with values at which one or more polygons should be drawn. The corresponding method is addpoly.default.

Author(s)

Wolfgang Viechtbauer [email protected] https://www.metafor-project.org

References

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⁠

See Also

addpoly.rma, addpoly.predict.rma, and addpoly.default for the specific method functions.

forest for functions to draw forest plots to which polygons can be added.

---

Add Polygons to Forest Plots (Default Method)

Description

Function to add one or more polygons to a forest plot.

Usage

## Default S3 method:
addpoly(x, vi, sei, ci.lb, ci.ub, pi.lb, pi.ub,
        rows=-1, level, annotate, digits, width, mlab,
        transf, atransf, targs, efac, col, border, lty, fonts, cex,
        constarea=FALSE, ...)

Arguments

x	vector with the values at which the polygons should be drawn.
vi	vector with the corresponding variances.
sei	vector with the corresponding standard errors (note: only one of the two, vi or sei, needs to be specified).
ci.lb	vector with the corresponding lower confidence interval bounds. Not needed if vi or sei is specified. See ‘Details’.
ci.ub	vector with the corresponding upper confidence interval bounds. Not needed if vi or sei is specified. See ‘Details’.
pi.lb	optional vector with the corresponding lower prediction interval bounds.
pi.ub	optional vector with the corresponding upper prediction interval bounds.
rows	vector to specify the rows (or more generally, the horizontal positions) for plotting the polygons (defaults is -1). Can also be a single value to specify the row (horizontal position) of the first polygon (the remaining polygons are then plotted below this starting row).
level	optional numeric value between 0 and 100 to specify the confidence interval level (see here for details).
annotate	optional logical to specify whether annotations should be added to the plot for the polygons that are drawn.
digits	optional integer to specify the number of decimal places to which the annotations should be rounded.
width	optional integer to manually adjust the width of the columns for the annotations.
mlab	optional character vector with the same length as x giving labels for the polygons that are drawn.
transf	optional argument to specify a function to transform the x values and confidence interval bounds (e.g., transf=exp; see also transf).
atransf	optional argument to specify a function to transform the annotations (e.g., atransf=exp; see also transf).
targs	optional arguments needed by the function specified via transf or atransf.
efac	optional vertical expansion factor for the polygons.
col	optional character string to specify the color of the polygons.
border	optional character string to specify the border color of the polygons.
lty	optional character string to specify the line type for the prediction interval.
fonts	optional character string to specify the font for the labels and annotations.
cex	optional symbol expansion factor.
constarea	logical to specify whether the height of the polygons (when adding multiple) should be adjusted so that the area of the polygons is constant (the default is FALSE).
...	other arguments.

Details

The function can be used to add one or more polygons to an existing forest plot created with the forest function. For example, summary estimates based on a model involving moderators can be added to the plot this way (see ‘Examples’).

To use the function, one should specify the values at which the polygons should be drawn (via the x argument) together with the corresponding variances (via the vi argument) or with the corresponding standard errors (via the sei argument). Alternatively, one can specify the values at which the polygons should be drawn together with the corresponding confidence interval bounds (via the ci.lb and ci.ub arguments). Optionally, one can also specify the bounds of the corresponding prediction interval bounds via the pi.lb and pi.ub arguments.

If unspecified, arguments level, annotate, digits, width, transf, atransf, targs, efac (only if the forest plot was created with forest.rma), fonts, cex, annosym, and textpos are automatically set equal to the same values that were used when creating the forest plot.

Author(s)

Wolfgang Viechtbauer [email protected] https://www.metafor-project.org

References

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⁠

See Also

forest for functions to draw forest plots to which polygons can be added.

Examples

### calculate log risk ratios and corresponding sampling variances
dat <- escalc(measure="RR", ai=tpos, bi=tneg, ci=cpos, di=cneg, data=dat.bcg)

### fit mixed-effects model with absolute latitude as a moderator
res <- rma(yi, vi, mods = ~ ablat, slab=paste(author, year, sep=", "), data=dat)

### forest plot of the observed risk ratios
forest(res, addfit=FALSE, atransf=exp, xlim=c(-9,5), ylim=c(-4.5,15), cex=0.9,
       order=ablat, ilab=ablat, ilab.lab="Lattitude", ilab.xpos=-4.5,
       header="Author(s) and Year", top=2)

### predicted average log risk ratios for 10, 30, and 50 degrees absolute latitude
x <- predict(res, newmods=c(10, 30, 50))

### add predicted average risk ratios to the forest plot
addpoly(x$pred, sei=x$se, rows=-2, mlab=c("- at 10 Degrees", "- at 30 Degrees", "- at 50 Degrees"))
abline(h=0)
text(-9, -1, "Model-Based Estimates:", pos=4, cex=0.9)

---

Add Polygons to Forest Plots (Method for 'predict.rma' Objects)

Description

Function to add one or more polygons to a forest plot based on an object of class "predict.rma".

Usage

## S3 method for class 'predict.rma'
addpoly(x, rows=-2, annotate,
        addpred=FALSE, digits, width, mlab, transf, atransf, targs,
        efac, col, border, lty, fonts, cex, constarea=FALSE, ...)

Arguments

x	an object of class "predict.rma".
rows	vector to specify the rows (or more generally, the horizontal positions) for plotting the polygons (defaults is -2). Can also be a single value to specify the row (horizontal position) of the first polygon (the remaining polygons are then plotted below this starting row).
annotate	optional logical to specify whether annotations should be added to the plot for the polygons that are drawn.
addpred	logical to specify whether the bounds of the prediction interval should be added to the plot (the default is FALSE).
digits	optional integer to specify the number of decimal places to which the annotations should be rounded.
width	optional integer to manually adjust the width of the columns for the annotations.
mlab	optional character vector with the same length as x giving labels for the polygons that are drawn.
transf	optional argument to specify a function to transform the x values and confidence interval bounds (e.g., transf=exp; see also transf).
atransf	optional argument to specify a function to transform the annotations (e.g., atransf=exp; see also transf).
targs	optional arguments needed by the function specified via transf or atransf.
efac	optional vertical expansion factor for the polygons.
col	optional character string to specify the color of the polygons.
border	optional character string to specify the border color of the polygons.
lty	optional character string to specify the line type for the prediction interval.
fonts	optional character string to specify the font for the labels and annotations.
cex	optional symbol expansion factor.
constarea	logical to specify whether the height of the polygons (when adding multiple) should be adjusted so that the area of the polygons is constant (the default is FALSE).
...	other arguments.

Details

The function can be used to add one or more polygons to an existing forest plot created with the forest function. For example, summary estimates based on a model involving moderators can be added to the plot this way (see ‘Examples’).

To use the function, one should specify the values at which the polygons should be drawn (via the x argument) together with the corresponding variances (via the vi argument) or with the corresponding standard errors (via the sei argument). Alternatively, one can specify the values at which the polygons should be drawn together with the corresponding confidence interval bounds (via the ci.lb and ci.ub arguments). Optionally, one can also specify the bounds of the corresponding prediction interval bounds via the pi.lb and pi.ub arguments.

If unspecified, arguments annotate, digits, width, transf, atransf, targs, efac (only if the forest plot was created with forest.rma), fonts, cex, annosym, and textpos are automatically set equal to the same values that were used when creating the forest plot.

Author(s)

Wolfgang Viechtbauer [email protected] https://www.metafor-project.org

References

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⁠

See Also

forest for functions to draw forest plots to which polygons can be added.

Examples

### calculate log risk ratios and corresponding sampling variances
dat <- escalc(measure="RR", ai=tpos, bi=tneg, ci=cpos, di=cneg,
              data=dat.bcg, slab=paste(author, year, sep=", "))

### forest plot of the observed risk ratios
with(dat, forest(yi, vi, atransf=exp, xlim=c(-9,5), ylim=c(-4.5,15),
                 at=log(c(0.05, 0.25, 1, 4)), cex=0.9, order=alloc,
                 ilab=alloc, ilab.lab="Allocation", ilab.xpos=-4.5,
                 header="Author(s) and Year", top=2))

### fit mixed-effects model with allocation method as a moderator
res <- rma(yi, vi, mods = ~ 0 + alloc, data=dat)

### predicted log risk ratios for the different allocation methods
x <- predict(res, newmods=diag(3))

### add predicted average risk ratios to the forest plot
addpoly(x, efac=1.3, col="gray", addpred=TRUE,
        mlab=c("Alternate Allocation", "Random Allocation", "Systematic Allocation"))
abline(h=0)
text(-9, -1, "Model-Based Estimates:", pos=4, cex=0.9, font=2)

---

Add Polygons to Forest Plots (Method for 'rma' Objects)

Description

Function to add a polygon to a forest plot showing the summary estimate with corresponding confidence interval based on an object of class "rma".

Usage

## S3 method for class 'rma'
addpoly(x, row=-2, level=x$level, annotate,
        addpred=FALSE, digits, width, mlab, transf, atransf, targs,
        efac, col, border, lty, fonts, cex, ...)

Arguments

x	an object of class "rma".
row	numeric value to specify the row (or more generally, the horizontal position) for plotting the polygon (the default is -2).
level	numeric value between 0 and 100 to specify the confidence interval level (see here for details). The default is to take the value from the object.
annotate	optional logical to specify whether annotations for the summary estimate should be added to the plot.
addpred	logical to specify whether the bounds of the prediction interval should be added to the plot (the default is FALSE).
digits	optional integer to specify the number of decimal places to which the annotations should be rounded.
width	optional integer to manually adjust the width of the columns for the annotations.
mlab	optional character string giving a label for the summary estimate polygon. If unspecified, the function sets a default label.
transf	optional argument to specify a function to transform the summary estimate and confidence interval bound (e.g., transf=exp; see also transf).
atransf	optional argument to specify a function to transform the annotations (e.g., atransf=exp; see also transf).
targs	optional arguments needed by the function specified via transf or atransf.
efac	optional vertical expansion factor for the polygon.
col	optional character string to specify the color of the polygon.
border	optional character string to specify the border color of the polygon.
lty	optional character string to specify the line type for the prediction interval.
fonts	optional character string to specify the font for the label and annotations.
cex	optional symbol expansion factor.
...	other arguments.

Details

The function can be used to add a four-sided polygon, sometimes called a summary ‘diamond’, to an existing forest plot created with the forest function. The polygon shows the summary estimate (with its confidence interval bounds) based on an equal- or a random-effects model. Using this function, summary estimates based on different types of models can be shown in the same plot. Also, summary estimates based on a subgrouping of the studies can be added to the plot this way. See ‘Examples’.

If unspecified, arguments annotate, digits, width, transf, atransf, targs, efac (only if the forest plot was created with forest.rma), fonts, cex, annosym, and textpos are automatically set equal to the same values that were used when creating the forest plot.

Author(s)

Wolfgang Viechtbauer [email protected] https://www.metafor-project.org

References

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⁠

See Also

forest for functions to draw forest plots to which polygons can be added.

Examples

### meta-analysis of the log risk ratios using the Mantel-Haenszel method
res <- rma.mh(measure="RR", ai=tpos, bi=tneg, ci=cpos, di=cneg, data=dat.bcg,
              slab=paste(author, year, sep=", "))

### forest plot of the observed risk ratios with summary estimate
forest(res, atransf=exp, xlim=c(-8,6), ylim=c(-2.5,15), header=TRUE, top=2)

### meta-analysis of the log risk ratios using a random-effects model
res <- rma(measure="RR", ai=tpos, bi=tneg, ci=cpos, di=cneg, data=dat.bcg)

### add the summary estimate from the random-effects model to the forest plot
addpoly(res)

### forest plot with subgrouping of studies and summaries per subgroup
dat <- escalc(measure="RR", ai=tpos, bi=tneg, ci=cpos, di=cneg, data=dat.bcg,
              slab=paste(author, year, sep=", "))
res <- rma(yi, vi, data=dat)
tmp <- forest(res, xlim=c(-16, 4.6), at=log(c(0.05, 0.25, 1, 4)), atransf=exp,
              ilab=cbind(tpos, tneg, cpos, cneg), ilab.lab=c("TB+","TB-","TB+","TB-"),
              ilab.xpos=c(-9.5,-8,-6,-4.5), cex=0.75, ylim=c(-1, 27), order=alloc,
              rows=c(3:4,9:15,20:23), mlab="RE Model for All Studies",
              header="Author(s) and Year")
op <- par(cex=tmp$cex)
text(c(-8.75,-5.25), tmp$ylim[2], c("Vaccinated", "Control"), font=2)
text(-16, c(24,16,5), c("Systematic Allocation", "Random Allocation",
                        "Alternate Allocation"), font=4, pos=4)
par(op)
res <- rma(yi, vi, data=dat, subset=(alloc=="systematic"))
addpoly(res, row=18.5, mlab="RE Model for Subgroup")
res <- rma(yi, vi, data=dat, subset=(alloc=="random"))
addpoly(res, row=7.5, mlab="RE Model for Subgroup")
res <- rma(yi, vi, data=dat, subset=(alloc=="alternate"))
addpoly(res, row=1.5, mlab="RE Model for Subgroup")

---

Aggregate Multiple Effect Sizes or Outcomes Within Studies

Description

Function to aggregate multiple effect sizes or outcomes belonging to the same study (or to the same level of some other clustering variable) into a single combined effect size or outcome.

Usage

## S3 method for class 'escalc'
aggregate(x, cluster, time, obs, V, struct="CS", rho, phi,
          weighted=TRUE, checkpd=TRUE, fun, na.rm=TRUE,
          addk=FALSE, subset, select, digits, var.names, ...)

Arguments

x	an object of class "escalc".
cluster	vector to specify the clustering variable (e.g., study).
time	optional vector to specify the time points (only relevant when struct="CAR", "CS+CAR", or "CS*CAR").
obs	optional vector to distinguish different observed effect sizes or outcomes measured at the same time point (only relevant when struct="CS*CAR").
V	optional argument to specify the variance-covariance matrix of the sampling errors. If unspecified, argument struct is used to specify the variance-covariance structure.
struct	character string to specify the variance-covariance structure of the sampling errors within the same cluster (either "ID", "CS", "CAR", "CS+CAR", or "CS*CAR"). See ‘Details’.
rho	value of the correlation of the sampling errors within clusters (when struct="CS", "CS+CAR", or "CS*CAR"). Can also be a vector with the value of the correlation for each cluster.
phi	value of the autocorrelation of the sampling errors within clusters (when struct="CAR", "CS+CAR", or "CS*CAR"). Can also be a vector with the value of the autocorrelation for each cluster.
weighted	logical to specify whether estimates within clusters should be aggregated using inverse-variance weighting (the default is TRUE). If set to FALSE, unweighted averages are computed.
checkpd	logical to specify whether to check that the variance-covariance matrices of the sampling errors within clusters are positive definite (the default is TRUE).
fun	optional list with three functions for aggregating other variables besides the effect sizes or outcomes within clusters (for numeric/integer variables, for logicals, and for all other types, respectively).
na.rm	logical to specify whether NA values should be removed before aggregating values within clusters (the default is TRUE). Can also be a vector with two logicals (the first pertaining to the effect sizes or outcomes, the second to all other variables).
addk	logical to specify whether to add the cluster size as a new variable (called ki) to the dataset (the default is FALSE).
subset	optional (logical or numeric) vector to specify the subset of rows to include when aggregating the effect sizes or outcomes.
select	optional vector to specify the names of the variables to include in the aggregated dataset.
digits	optional integer to specify the number of decimal places to which the printed results should be rounded (the default is to take the value from the object).
var.names	optional character vector with two elements to specify the name of the variable that contains the observed effect sizes or outcomes and the name of the variable with the corresponding sampling variances (when unspecified, the function attempts to set these automatically based on the object).
...	other arguments.

Details

In many meta-analyses, multiple effect sizes or outcomes can be extracted from the same study. Ideally, such structures should be analyzed using an appropriate multilevel/multivariate model as can be fitted with the rma.mv function. However, there may occasionally be reasons for aggregating multiple effect sizes or outcomes belonging to the same study (or to the same level of some other clustering variable) into a single combined effect size or outcome. The present function can be used for this purpose.

The input must be an object of class "escalc". The error ‘Error in match.fun(FUN): argument "FUN" is missing, with no default’ indicates that a regular data frame was passed to the function, but this does not work. One can turn a regular data frame (containing the effect sizes or outcomes and the corresponding sampling variances) into an "escalc" object with the escalc function. See the ‘Examples’ below for an illustration of this.

The cluster variable is used to specify which estimates/outcomes belong to the same study/cluster.

In the simplest case, the estimates/outcomes within clusters (or, to be precise, their sampling errors) are assumed to be independent. This is usually a safe assumption as long as each study participant (or whatever the study units are) only contributes data to a single estimate/outcome. For example, if a study provides effect size estimates for male and female subjects separately, then the sampling errors can usually be assumed to be independent. In this case, one can set struct="ID" and multiple estimates/outcomes within the same cluster are combined using standard inverse-variance weighting (i.e., using weighted least squares) under the assumption of independence.

In other cases, the estimates/outcomes within clusters cannot be assumed to be independent. For example, if multiple effect size estimates are computed for the same group of subjects (e.g., based on different scales to measure some construct of interest), then the estimates are likely to be correlated. If the actual correlation between the estimates is unknown, one can often still make an educated guess and set argument rho to this value, which is then assumed to be the same for all pairs of estimates within clusters when struct="CS" (for a compound symmetric structure). Multiple estimates/outcomes within the same cluster are then combined using inverse-variance weighting taking their correlation into consideration (i.e., using generalized least squares). One can also specify a different value of rho for each cluster by passing a vector (of the same length as the number of clusters) to this argument.

If multiple effect size estimates are computed for the same group of subjects at different time points, then it may be more sensible to assume that the correlation between estimates decreases as a function of the distance between the time points. If so, one can specify struct="CAR" (for a continuous-time autoregressive structure), set phi to the autocorrelation (for two estimates one time-unit apart), and use argument time to specify the actual time points corresponding to the estimates. The correlation between two estimates, \(y_{it}\) and \(y_{it'}\), in the \(i\textrm{th}\) cluster, with time points \(\textrm{time}_{it}\) and \(\textrm{time}_{it'}\), is then given by \(\phi^{|\textrm{time}_{it} - \textrm{time}_{it'}|}\). One can also specify a different value of phi for each cluster by passing a vector (of the same length as the number of clusters) to this argument.

One can also combine the compound symmetric and autoregressive structures if there are multiple time points and multiple observed effect sizes or outcomes at these time points. One option is struct="CS+CAR". In this case, one must specify the time argument and both rho and phi. The correlation between two estimates, \(y_{it}\) and \(y_{it'}\), in the \(i\textrm{th}\) cluster, with time points \(\textrm{time}_{it}\) and \(\textrm{time}_{it'}\), is then given by \(\rho + (1 - \rho) \phi^{|\textrm{time}_{it} - \textrm{time}_{it'}|}\).

Alternatively, one can specify struct="CS*CAR". In this case, one must specify both the time and obs arguments and both rho and phi. The correlation between two estimates, \(y_{ijt}\) and \(y_{ijt'}\), with the same value for obs but different values for time, is then given by \(\phi^{|\textrm{time}_{ijt} - \textrm{time}_{ijt'}|}\), the correlation between two estimates, \(y_{ijt}\) and \(y_{ij't}\), with different values for obs but the same value for time, is then given by \(\rho\), and the correlation between two estimates, \(y_{ijt}\) and \(y_{ij't'}\), with different values for obs and different values for time, is then given by \(\rho \times \phi^{|\textrm{time}_{ijt} - \textrm{time}_{ijt'}|}\).

Finally, if one actually knows the correlation (and hence the covariance) between each pair of estimates (or has an approximation thereof), one can also specify the entire variance-covariance matrix of the estimates (or more precisely, their sampling errors) via the V argument (in this case, arguments struct, time, obs, rho, and phi are ignored). Note that the vcalc function can be used to construct such a V matrix and provides even more flexibility for specifying various types of dependencies. See the ‘Examples’ below for an illustration of this.

Instead of using inverse-variance weighting (i.e., weighted/generalized least squares) to combine the estimates within clusters, one can set weighted=FALSE in which case the estimates are averaged within clusters without any weighting (although the correlations between estimates as specified are still taken into consideration).

Other variables (besides the estimates) will also be aggregated to the cluster level. By default, numeric/integer type variables are averaged, logicals are also averaged (yielding the proportion of TRUE values), and for all other types of variables (e.g., character variables or factors) the most frequent category/level is returned. One can also specify a list of three functions via the fun argument for aggregating variables belonging to these three types.

Argument na.rm controls how missing values should be handled. By default, any missing estimates are first removed before aggregating the non-missing values within each cluster. The same applies when aggregating the other variables. One can also specify a vector with two logicals for the na.rm argument to control how missing values should be handled when aggregating the estimates and when aggregating all other variables.

Value

An object of class c("escalc","data.frame") that contains the (selected) variables aggregated to the cluster level.

The object is formatted and printed with the print function.

Author(s)

Wolfgang Viechtbauer [email protected] https://www.metafor-project.org

References

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⁠

See Also

escalc for a function to create escalc objects.

Examples

### copy data into 'dat' and examine data
dat <- dat.konstantopoulos2011
head(dat, 11)

### aggregate estimates to the district level, assuming independent sampling
### errors for multiples studies/schools within the same district
agg <- aggregate(dat, cluster=district, struct="ID", addk=TRUE)
agg

### copy data into 'dat' and examine data
dat <- dat.assink2016
head(dat, 19)

### note: 'dat' is an 'escalc' object
class(dat)

### turn 'dat' into a regular data frame
dat <- as.data.frame(dat)
class(dat)

### turn data frame into an 'escalc' object
dat <- escalc(measure="SMD", yi=yi, vi=vi, data=dat)
class(dat)

### aggregate the estimates to the study level, assuming a CS structure for
### the sampling errors within studies with a correlation of 0.6
agg <- aggregate(dat, cluster=study, rho=0.6)
agg

### use vcalc() and then the V argument
V <- vcalc(vi, cluster=study, obs=esid, data=dat, rho=0.6)
agg <- aggregate(dat, cluster=study, V=V)
agg

### use a correlation of 0.7 for effect sizes corresponding to the same type of
### delinquent behavior and a correlation of 0.5 for effect sizes corresponding
### to different types of delinquent behavior
V <- vcalc(vi, cluster=study, type=deltype, obs=esid, data=dat, rho=c(0.7, 0.5))
agg <- aggregate(dat, cluster=study, V=V)
agg

### reshape 'dat.ishak2007' into long format
dat <- dat.ishak2007
dat <- reshape(dat.ishak2007, direction="long", idvar="study", v.names=c("yi","vi"),
               varying=list(c(2,4,6,8), c(3,5,7,9)))
dat <- dat[order(study, time),]
dat <- dat[!is.na(yi),]
rownames(dat) <- NULL
head(dat, 8)

### aggregate the estimates to the study level, assuming a CAR structure for
### the sampling errors within studies with an autocorrelation of 0.9
agg <- aggregate(dat, cluster=study, struct="CAR", time=time, phi=0.9)
head(agg, 5)

---

Likelihood Ratio and Wald-Type Tests for 'rma' Objects

Description

For two (nested) models of class "rma.uni" or "rma.mv", the function provides a full versus reduced model comparison in terms of model fit statistics and a likelihood ratio test. When a single model is specified, a Wald-type test of one or more model coefficients or linear combinations thereof is carried out.

Usage

## S3 method for class 'rma'
anova(object, object2, btt, X, att, Z, rhs, adjust, digits, refit=FALSE, ...)

Arguments

object	an object of class "rma.uni" or "rma.mv".
object2	an (optional) object of class "rma.uni" or "rma.mv". Only relevant when conducting a model comparison and likelihood ratio test. See ‘Details’.
btt	optional vector of indices (or list thereof) to specify which coefficients should be included in the Wald-type test. Can also be a string to grep for. See ‘Details’.
X	optional numeric vector or matrix to specify one or more linear combinations of the coefficients in the model that should be tested. See ‘Details’.
att	optional vector of indices (or list thereof) to specify which scale coefficients should be included in the Wald-type test. Can also be a string to grep for. See ‘Details’. Only relevant for location-scale models (see rma.uni).
Z	optional numeric vector or matrix to specify one or more linear combinations of the scale coefficients in the model that should be tested. See ‘Details’. Only relevant for location-scale models (see rma.uni).
rhs	optional scalar or vector of values for the right-hand side of the null hypothesis when testing a set of coefficients (via btt or att) or linear combinations thereof (via X or Z). If unspecified, this defaults to a vector of zeros of the appropriate length. See ‘Details’.
adjust	optional argument to specify (as a character string) a method for adjusting the p-values of Wald-type tests for multiple testing. See p.adjust for possible options. Can be abbreviated. Can also be a logical and if TRUE, then a Bonferroni correction is used.
digits	optional integer to specify the number of decimal places to which the printed results should be rounded. If unspecified, the default is to take the value from the object.
refit	logical to specify whether models fitted with REML estimation and differing in their fixed effects should be refitted with ML estimation when conducting a likelihood ratio test (the default is FALSE).
...	other arguments.

Details

The function can be used in three different ways:

1. When a single model is specified (via argument object), the function provides a Wald-type test of one or more model coefficients, that is, \[\mbox{H}_0{:}\; \beta_{j \in \texttt{btt}} = 0,\] where \(\beta_{j \in \texttt{btt}}\) is the set of coefficients to be tested (by default whether the set of coefficients is significantly different from zero, but one can specify a different value under the null hypothesis via argument rhs).

   In particular, for equal- or random-effects models (i.e., models without moderators), this is just the test of the single coefficient of the model (i.e., \(\mbox{H}_0{:}\; \theta = 0\) or \(\mbox{H}_0{:}\; \mu = 0\)). For models including moderators, an omnibus test of all model coefficients is conducted that excludes the intercept (the first coefficient) if it is included in the model. If no intercept is included in the model, then the omnibus test includes all coefficients in the model including the first.

   Alternatively, one can manually specify the indices of the coefficients to test via the btt (‘betas to test’) argument. For example, with btt=c(3,4), only the third and fourth coefficients from the model are included in the test (if an intercept is included in the model, then it corresponds to the first coefficient in the model). Instead of specifying the coefficient numbers, one can specify a string for btt. In that case, grep will be used to search for all coefficient names that match the string (and hence, one can use regular expressions to fine-tune the search for matching strings). Using the btt argument, one can for example select all coefficients corresponding to a particular factor to test if the factor as a whole is significant. One can also specify a list of indices/strings, in which case tests of all list elements will be conducted. See ‘Examples’.

   For location-scale models fitted with the rma.uni function, one can use the att argument in an analogous manner to specify the indices of the scale coefficients to test (i.e., \(\mbox{H}_0{:}\; \alpha_{j \in \texttt{att}} = 0\), where \(\alpha_{j \in \texttt{att}}\) is the set of coefficients to be tested).

2. When a single model is specified (via argument object), one can use the X argument\(^1\) to specify a linear combination of the coefficients in the model that should be tested using a Wald-type test, that is, \[\mbox{H}_0{:}\; \tilde{x} \beta = 0,\] where \(\tilde{x}\) is a (row) vector of the same length as there are coefficients in the model (by default whether the linear combination is significantly different from zero, but one can specify a different value under the null hypothesis via argument rhs). One can also specify a matrix of linear combinations via the X argument to test \[\mbox{H}_0{:}\; \tilde{X} \beta = 0,\] where each row of \(\tilde{X}\) defines a particular linear combination to be tested (if rhs is used, then it should either be a scalar or of the same length as the number of combinations to be tested). If the matrix is of full rank, an omnibus Wald-type test of all linear combinations is also provided. Linear combinations can also be obtained with the predict function, which provides corresponding confidence intervals. See also the pairwise function for constructing a matrix of pairwise contrasts for testing the levels of a categorical moderator against each other.

   For location-scale models fitted with the rma.uni function, one can use the Z argument in an analogous manner to specify one or multiple linear combinations of the scale coefficients in the model that should be tested (i.e., \(\mbox{H}_0{:}\; \tilde{Z} \alpha = 0\)).

3. When specifying two models for comparison (via arguments object and object2), the function provides a likelihood ratio test (LRT) comparing the two models. The two models must be based on the same set of data, must be of the same class, and should be nested for the LRT to make sense. Also, LRTs are not meaningful when using REML estimation and the two models differ in terms of their fixed effects (setting refit=TRUE automatically refits the two models using ML estimation). Also, the theory underlying LRTs is only really applicable when comparing models that were fitted with ML/REML estimation, so if some other estimation method was used to fit the two models, the results should be treated with caution.

———

\(^1\) This argument used to be called L, but was renamed to X (but using L in place of X still works).

Value

An object of class "anova.rma". When a single model is specified (without any further arguments or together with the btt or att argument), the object is a list containing the following components:

QM	test statistic of the Wald-type test of the model coefficients.
QMdf	corresponding degrees of freedom.
QMp	corresponding p-value.
btt	indices of the coefficients tested by the Wald-type test.
k	number of outcomes included in the analysis.
p	number of coefficients in the model (including the intercept).
m	number of coefficients included in the Wald-type test.
...	some additional elements/values.

When btt or att was a list, then the object is a list of class "list.anova.rma", where each element is an "anova.rma" object as described above.

When argument X is used, the object is a list containing the following components:

QM	test statistic of the omnibus Wald-type test of all linear combinations.
QMdf	corresponding degrees of freedom.
QMp	corresponding p-value.
hyp	description of the linear combinations tested.
Xb	values of the linear combinations.
se	standard errors of the linear combinations.
zval	test statistics of the linear combinations.
pval	corresponding p-values.

When two models are specified, the object is a list containing the following components:

fit.stats.f	log-likelihood, deviance, AIC, BIC, and AICc for the full model.
fit.stats.r	log-likelihood, deviance, AIC, BIC, and AICc for the reduced model.
parms.f	number of parameters in the full model.
parms.r	number of parameters in the reduced model.
LRT	likelihood ratio test statistic.
pval	corresponding p-value.
QE.f	test statistic of the test for (residual) heterogeneity from the full model.
QE.r	test statistic of the test for (residual) heterogeneity from the reduced model.
tau2.f	estimated \(\tau^2\) value from the full model. NA for "rma.mv" objects.
tau2.r	estimated \(\tau^2\) value from the reduced model. NA for "rma.mv" objects.
R2	amount (in percent) of the heterogeneity in the reduced model that is accounted for in the full model (NA for "rma.mv" objects). This can be regarded as a pseudo \(R^2\) statistic (Raudenbush, 2009). Note that the value may not be very accurate unless \(k\) is large (Lopez-Lopez et al., 2014).
...	some additional elements/values.

The results are formatted and printed with the print function. To format the results as a data frame, one can use the as.data.frame function.

Note

The function can also be used to conduct a likelihood ratio test (LRT) for the amount of (residual) heterogeneity in random- and mixed-effects models. The full model should then be fitted with either method="ML" or method="REML" and the reduced model with method="EE" (or with tau2=0). The p-value for the test is based on a chi-square distribution with 1 degree of freedom, but actually needs to be adjusted for the fact that the parameter (i.e., \(\tau^2\)) falls on the boundary of the parameter space under the null hypothesis (see Viechtbauer, 2007, for more details).

LRTs for variance components in more complex models (as fitted with the rma.mv function) can also be conducted in this manner (see ‘Examples’).

Author(s)

Wolfgang Viechtbauer [email protected] https://www.metafor-project.org

References

Hardy, R. J., & Thompson, S. G. (1996). A likelihood approach to meta-analysis with random effects. Statistics in Medicine, 15(6), 619–629. ⁠https://doi.org/10.1002/(sici)1097-0258(19960330)15:6%3C619::aid-sim188%3E3.0.co;2-a⁠

Huizenga, H. M., Visser, I., & Dolan, C. V. (2011). Testing overall and moderator effects in random effects meta-regression. British Journal of Mathematical and Statistical Psychology, 64(1), 1–19. ⁠https://doi.org/10.1348/000711010X522687⁠

López-López, J. A., Marín-Martínez, F., Sánchez-Meca, J., Van den Noortgate, W., & Viechtbauer, W. (2014). Estimation of the predictive power of the model in mixed-effects meta-regression: A simulation study. British Journal of Mathematical and Statistical Psychology, 67(1), 30–48. ⁠https://doi.org/10.1111/bmsp.12002⁠

Raudenbush, S. W. (2009). Analyzing effect sizes: Random effects models. In H. Cooper, L. V. Hedges, & J. C. Valentine (Eds.), The handbook of research synthesis and meta-analysis (2nd ed., pp. 295–315). New York: Russell Sage Foundation.

Viechtbauer, W. (2007). Hypothesis tests for population heterogeneity in meta-analysis. British Journal of Mathematical and Statistical Psychology, 60(1), 29–60. ⁠https://doi.org/10.1348/000711005X64042⁠

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., & López-López, J. A. (2022). Location-scale models for meta-analysis. Research Synthesis Methods. 13(6), 697–715. ⁠https://doi.org/10.1002/jrsm.1562⁠

See Also

rma.uni and rma.mv for functions to fit models for which likelihood ratio and Wald-type tests can be conducted.

print for the print method and as.data.frame for the method to format the results as a data frame.

Examples

### calculate log risk ratios and corresponding sampling variances
dat <- escalc(measure="RR", ai=tpos, bi=tneg, ci=cpos, di=cneg, data=dat.bcg)

### fit random-effects model
res1 <- rma(yi, vi, data=dat, method="ML")
res1

### fit mixed-effects model with two moderators (absolute latitude and publication year)
res2 <- rma(yi, vi, mods = ~ ablat + year, data=dat, method="ML")
res2

### Wald-type test of the two moderators
anova(res2)

### alternative way of specifying the same test
anova(res2, X=rbind(c(0,1,0), c(0,0,1)))

### corresponding likelihood ratio test
anova(res1, res2)

### Wald-type test of a linear combination
anova(res2, X=c(1,35,1970))

### use predict() to obtain the same linear combination (with its CI)
predict(res2, newmods=c(35,1970))

### Wald-type tests of several linear combinations
anova(res2, X=cbind(1,seq(0,60,by=10),1970))

### adjust for multiple testing with the Bonferroni method
anova(res2, X=cbind(1,seq(0,60,by=10),1970), adjust="bonf")

### mixed-effects model with three moderators
res3 <- rma(yi, vi, mods = ~ ablat + year + alloc, data=dat, method="ML")
res3

### Wald-type test of the 'alloc' factor
anova(res3, btt=4:5)

### instead of specifying the coefficient numbers, grep for "alloc"
anova(res3, btt="alloc")

### specify a list for the 'btt' argument
anova(res3, btt=list(2,3,4:5))

### adjust for multiple testing with the Bonferroni method
anova(res3, btt=list(2,3,4:5), adjust="bonf")

############################################################################

### an example of doing LRTs of variance components in more complex models
dat <- dat.konstantopoulos2011
res <- rma.mv(yi, vi, random = ~ 1 | district/school, data=dat)

### likelihood ratio test of the district-level variance component
res0 <- rma.mv(yi, vi, random = ~ 1 | district/school, data=dat, sigma2=c(0,NA))
anova(res, res0)

### likelihood ratio test of the school-level variance component
res0 <- rma.mv(yi, vi, random = ~ 1 | district/school, data=dat, sigma2=c(NA,0))
anova(res, res0)

### likelihood ratio test of both variance components simultaneously
res0 <- rma.mv(yi, vi, data=dat)
anova(res, res0)

############################################################################

### an example illustrating a workflow involving cluster-robust inference
dat <- dat.assink2016

### assume that the effect sizes within studies are correlated with rho=0.6
V <- vcalc(vi, cluster=study, obs=esid, data=dat, rho=0.6)

### fit multilevel model using this approximate V matrix
res <- rma.mv(yi, V, random = ~ 1 | study/esid, data=dat)
res

### likelihood ratio tests of the two variance components
res0 <- rma.mv(yi, V, random = ~ 1 | study/esid, data=dat, sigma2=c(0,NA))
anova(res, res0)
res0 <- rma.mv(yi, V, random = ~ 1 | study/esid, data=dat, sigma2=c(NA,0))
anova(res, res0)

### use cluster-robust methods for inferences about the fixed effects
sav <- robust(res, cluster=study, clubSandwich=TRUE)
sav

### examine if 'deltype' is a potential moderator
res <- rma.mv(yi, V, mods = ~ deltype, random = ~ 1 | study/esid, data=dat)
sav <- robust(res, cluster=study, clubSandwich=TRUE)
sav

### note: the (denominator) dfs for the omnibus F-test are very low, so the results
### of this test may not be trustworthy; consider using cluster wild bootstrapping
## Not run: 
library(wildmeta)
Wald_test_cwb(res, constraints=constrain_zero(2:3), R=1000, seed=1234)

## End(Not run)

---

Baujat Plots for 'rma' Objects

Description

Function to create Baujat plots for objects of class "rma".

Usage

baujat(x, ...)

## S3 method for class 'rma'
baujat(x, xlim, ylim, xlab, ylab, cex, symbol="ids", grid=TRUE, progbar=FALSE, ...)

Arguments

x	an object of class "rma".
xlim	x-axis limits. If unspecified, the function sets the x-axis limits to some sensible values.
ylim	y-axis limits. If unspecified, the function sets the y-axis limits to some sensible values.
xlab	title for the x-axis. If unspecified, the function sets an appropriate axis title.
ylab	title for the y-axis. If unspecified, the function sets an appropriate axis title.
cex	symbol/character expansion factor.
symbol	either an integer to specify the pch value (i.e., plotting symbol), or "slab" to plot the study labels, or "ids" (the default) to plot the study id numbers.
grid	logical to specify whether a grid should be added to the plot. Can also be a color name.
progbar	logical to specify whether a progress bar should be shown (the default is FALSE).
...	other arguments.

Details

The model specified via x must be a model fitted with either the rma.uni, rma.mh, or rma.peto functions.

Baujat et al. (2002) proposed a diagnostic plot to detect sources of heterogeneity in meta-analytic data. The plot shows the contribution of each study to the overall \(Q\)-test statistic for heterogeneity on the x-axis versus the influence of each study (defined as the standardized squared difference between the overall estimate based on an equal-effects model with and without the study included in the model fitting) on the y-axis. The same type of plot can be produced by first fitting an equal-effects model with either the rma.uni (using method="EE"), rma.mh, or rma.peto functions and then passing the fitted model object to the baujat function.

For models fitted with the rma.uni function (which may be random-effects or mixed-effects meta-regressions models), the idea underlying this type of plot can be generalized as described by Viechtbauer (2021): The x-axis then corresponds to the squared Pearson residual of a study, while the y-axis corresponds to the standardized squared difference between the predicted/fitted value for the study with and without the study included in the model fitting.

By default, the points plotted are the study id numbers, but one can also plot the study labels by setting symbol="slab" (if study labels are available within the model object) or one can specify a plotting symbol via the symbol argument that gets passed to pch (see points for possible options).

Value

A data frame with components:

x	the x-axis coordinates of the points that were plotted.
y	the y-axis coordinates of the points that were plotted.
ids	the study id numbers.
slab	the study labels.

Note that the data frame is returned invisibly.

Author(s)

Wolfgang Viechtbauer [email protected] https://www.metafor-project.org

References

Baujat, B., Mahe, C., Pignon, J.-P., & Hill, C. (2002). A graphical method for exploring heterogeneity in meta-analyses: Application to a meta-analysis of 65 trials. Statistics in Medicine, 21(18), 2641–2652. ⁠https://doi.org/10.1002/sim.1221⁠

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. (2021). Model checking in meta-analysis. In C. H. Schmid, T. Stijnen, & I. R. White (Eds.), Handbook of meta-analysis (pp. 219–254). Boca Raton, FL: CRC Press. ⁠https://doi.org/10.1201/9781315119403⁠

See Also

rma.uni, rma.mh and rma.peto for functions to fit models for which Baujat plots can be created.

influence for other model diagnostics.

Examples

### copy data from Pignon et al. (2000) into 'dat'
dat <- dat.pignon2000

### calculate estimated log hazard ratios and sampling variances
dat$yi <- with(dat, OmE/V)
dat$vi <- with(dat, 1/V)

### meta-analysis based on all 65 trials
res <- rma(yi, vi, data=dat, method="EE", slab=trial)

### create Baujat plot
baujat(res)

### some variations of the plotting symbol
baujat(res, symbol=19)
baujat(res, symbol="slab")

### label only a selection of the more 'extreme' points
sav <- baujat(res, symbol=19, xlim=c(0,20))
sav <- sav[sav$x >= 10 | sav$y >= 0.10,]
text(sav$x, sav$y, sav$slab, pos=1, cex=0.8)

---

Construct Block Diagonal Matrix

Description

Function to construct a block diagonal matrix from (a list of) matrices.

Usage

bldiag(..., order)

Arguments

...	individual matrices or a list of matrices.
order	optional argument to specify a variable based on which a square block diagonal matrix should be ordered.

Author(s)

Posted to R-help by Berton Gunter (2 Sep 2005) with some further adjustments by Wolfgang Viechtbauer

See Also

rma.mv for the model fitting function that can take such a block diagonal matrix as input (for the V argument).

blsplit for a function that can split a block diagonal matrix into a list of sub-matrices.

Examples

### copy data into 'dat'
dat <- dat.berkey1998
dat

### construct list with the variance-covariance matrices of the observed outcomes for the studies
V <- lapply(split(dat[c("v1i","v2i")], dat$trial), as.matrix)
V

### construct block diagonal matrix
V <- bldiag(V)
V

### if we split based on 'author', the list elements in V are in a different order than tha data
V <- lapply(split(dat[c("v1i","v2i")], dat$author), as.matrix)
V

### can use 'order' argument to reorder the block-diagonal matrix into the correct order
V <- bldiag(V, order=dat$author)
V

---

Split Block Diagonal Matrix

Description

Function to split a block diagonal matrix into a list of sub-matrices.

Usage

blsplit(x, cluster, fun, args, sort=FALSE)

Arguments

x	a block diagonal matrix.
cluster	vector to specify the clustering variable to use for splitting.
fun	optional argument to specify a function to apply to each sub-matrix.
args	optional argument to specify any additional argument(s) for the function specified via fun.
sort	logical to specify whether to sort the list by the unique cluster values (the default is FALSE).

Value

A list of one or more sub-matrices.

Author(s)

Wolfgang Viechtbauer [email protected] https://www.metafor-project.org

See Also

bldiag for a function to create a block diagonal matrix based on sub-matrices.

vcalc for a function to construct a variance-covariance matrix of dependent effect sizes or outcomes, which often has a block diagonal structure.

Examples

### copy data into 'dat'
dat <- dat.assink2016

### assume that the effect sizes within studies are correlated with rho=0.6
V <- vcalc(vi, cluster=study, obs=esid, data=dat, rho=0.6)

### split V matrix into list of sub-matrices
Vs <- blsplit(V, cluster=dat$study)
Vs[1:2]
lapply(Vs[1:2], cov2cor)

### illustrate the use of the fun and args arguments
blsplit(V, cluster=dat$study, cov2cor)[1:2]
blsplit(V, cluster=dat$study, round, 3)[1:2]

---

Best Linear Unbiased Predictions for 'rma.uni' Objects

Description

Function to compute best linear unbiased predictions (BLUPs) of the study-specific true effect sizes or outcomes (by combining the fitted values based on the fixed effects and the estimated contributions of the random effects) for objects of class "rma.uni". Corresponding standard errors and prediction interval bounds are also provided.

Usage

blup(x, ...)

## S3 method for class 'rma.uni'
blup(x, level, digits, transf, targs, ...)

Arguments

x	an object of class "rma.uni".
level	numeric value between 0 and 100 to specify the prediction interval level (see here for details). If unspecified, the default is to take the value from the object.
digits	optional integer to specify the number of decimal places to which the printed results should be rounded. If unspecified, the default is to take the value from the object.
transf	optional argument to specify a function to transform the predicted values and interval bounds (e.g., transf=exp; see also transf). If unspecified, no transformation is used.
targs	optional arguments needed by the function specified under transf.
...	other arguments.

Value

An object of class "list.rma". The object is a list containing the following components:

pred	predicted values.
se	corresponding standard errors.
pi.lb	lower bound of the prediction intervals.
pi.ub	upper bound of the prediction intervals.
...	some additional elements/values.

The object is formatted and printed with the print function. To format the results as a data frame, one can use the as.data.frame function.

Note

For best linear unbiased predictions of only the random effects, see ranef.

For predicted/fitted values that are based only on the fixed effects of the model, see fitted and predict.

For conditional residuals (the deviations of the observed effect sizes or outcomes from the BLUPs), see rstandard.rma.uni with type="conditional".

Equal-effects models do not contain random study effects. The BLUPs for these models will therefore be equal to the fitted values, that is, those obtained with fitted and predict.

When using the transf argument, the transformation is applied to the predicted values and the corresponding interval bounds. The standard errors are then set equal to NA and are omitted from the printed output.

By default, a standard normal distribution is used to construct the prediction intervals. When the model was fitted with test="t", test="knha", test="hksj", or test="adhoc", then a t-distribution with \(k-p\) degrees of freedom is used.

To be precise, it should be noted that the function actually computes empirical BLUPs (eBLUPs), since the predicted values are a function of the estimated value of \(\tau^2\).

Author(s)

Wolfgang Viechtbauer [email protected] https://www.metafor-project.org

References

Kackar, R. N., & Harville, D. A. (1981). Unbiasedness of two-stage estimation and prediction procedures for mixed linear models. Communications in Statistics, Theory and Methods, 10(13), 1249–1261. ⁠https://doi.org/10.1080/03610928108828108⁠

Raudenbush, S. W., & Bryk, A. S. (1985). Empirical Bayes meta-analysis. Journal of Educational Statistics, 10(2), 75–98. ⁠https://doi.org/10.3102/10769986010002075⁠

Robinson, G. K. (1991). That BLUP is a good thing: The estimation of random effects. Statistical Science, 6(1), 15–32. ⁠https://doi.org/10.1214/ss/1177011926⁠

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⁠

See Also

rma.uni for the function to fit models for which BLUPs can be extracted.

predict and fitted for functions to compute the predicted/fitted values based only on the fixed effects and ranef for a function to compute the BLUPs based only on the random effects.

Examples

### calculate log risk ratios and corresponding sampling variances
dat <- escalc(measure="RR", ai=tpos, bi=tneg, ci=cpos, di=cneg, data=dat.bcg)

### meta-analysis of the log risk ratios using a random-effects model
res <- rma(yi, vi, data=dat)

### BLUPs of the true risk ratios for each study
blup(res, transf=exp)

### illustrate shrinkage of BLUPs towards the (estimated) population average
res <- rma(yi, vi, data=dat)
blups <- blup(res)$pred
plot(NA, NA, xlim=c(.8,2.4), ylim=c(-2,0.5), pch=19,
     xaxt="n", bty="n", xlab="", ylab="Log Risk Ratio")
segments(rep(1,13), dat$yi, rep(2,13), blups, col="darkgray")
points(rep(1,13), dat$yi, pch=19)
points(rep(2,13), blups, pch=19)
axis(side=1, at=c(1,2), labels=c("Observed\nValues", "BLUPs"), lwd=0)
segments(0, res$beta, 2.15, res$beta, lty="dotted")
text(2.3, res$beta, substitute(hat(mu)==muhat, list(muhat=round(res$beta[[1]], 2))), cex=1)

---

Extract the Model Coefficients and Variance-Covariance Matrix from 'matreg' Objects

Description

Methods for objects of class "matreg".

Usage

## S3 method for class 'matreg'
coef(object, ...)
## S3 method for class 'matreg'
vcov(object, ...)

Arguments

object	an object of class "matreg".
...	other arguments.

Details

The coef function extracts the estimated model coefficients from objects of class "matreg". The vcov function extracts the corresponding variance-covariance matrix.

Value

Either a vector with the estimated model coefficients or a variance-covariance matrix.

Author(s)

Wolfgang Viechtbauer [email protected] https://www.metafor-project.org

References

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⁠

See Also

matreg for the function to create matreg objects.

Examples

### fit a regression model with lm() to the 'mtcars' dataset
res <- lm(mpg ~ hp + wt + am, data=mtcars)
coef(res)
vcov(res)

### covariance matrix of the dataset
S <- cov(mtcars)

### fit the same regression model using matreg()
res <- matreg(y="mpg", x=c("hp","wt","am"), R=S, cov=TRUE,
              means=colMeans(mtcars), n=nrow(mtcars))
coef(res)
vcov(res)

---

Extract the Model Coefficient Table from 'permutest.rma.uni' Objects

Description

Function to extract the estimated model coefficients, corresponding standard errors, test statistics, p-values (based on the permutation tests), and confidence interval bounds from objects of class "permutest.rma.uni".

Usage

## S3 method for class 'permutest.rma.uni'
coef(object, ...)

Arguments

object	an object of class "permutest.rma.uni".
...	other arguments.

Value

A data frame with the following elements:

estimate	estimated model coefficient(s).
se	corresponding standard error(s).
zval	corresponding test statistic(s).
pval	p-value(s) based on the permutation test(s).
ci.lb	lower bound of the (permutation-based) confidence interval(s).
ci.ub	upper bound of the (permutation-based) confidence interval(s).

When the model was fitted with test="t", test="knha", test="hksj", or test="adhoc", then zval is called tval in the data frame that is returned by the function.

Author(s)

Wolfgang Viechtbauer [email protected] https://www.metafor-project.org

References

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⁠

See Also

permutest for the function to conduct permutation tests and rma.uni for the function to fit models for which permutation tests can be conducted.

Examples

### calculate log risk ratios and corresponding sampling variances
dat <- escalc(measure="RR", ai=tpos, bi=tneg, ci=cpos, di=cneg, data=dat.bcg)

### fit mixed-effects model with absolute latitude and publication year as moderators
res <- rma(yi, vi, mods = ~ ablat + year, data=dat)

### carry out permutation test
## Not run: 
set.seed(1234) # for reproducibility
sav <- permutest(res)
coef(sav)

## End(Not run)

---

Extract the Model Coefficients and Coefficient Table from 'rma' and 'summary.rma' Objects

Description

Function to extract the estimated model coefficients from objects of class "rma". For objects of class "summary.rma", the model coefficients, corresponding standard errors, test statistics, p-values, and confidence interval bounds are extracted.

Usage

## S3 method for class 'rma'
coef(object, ...)
## S3 method for class 'summary.rma'
coef(object, ...)

Arguments

object	an object of class "rma" or "summary.rma".
...	other arguments.

Value

Either a vector with the estimated model coefficient(s) or a data frame with the following elements:

estimate	estimated model coefficient(s).
se	corresponding standard error(s).
zval	corresponding test statistic(s).
pval	corresponding p-value(s).
ci.lb	corresponding lower bound of the confidence interval(s).
ci.ub	corresponding upper bound of the confidence interval(s).

When the model was fitted with test="t", test="knha", test="hksj", or test="adhoc", then zval is called tval in the data frame that is returned by the function.

Author(s)

Wolfgang Viechtbauer [email protected] https://www.metafor-project.org

References

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⁠

See Also

rma.uni, rma.mh, rma.peto, rma.glmm, and rma.mv for functions to fit models for which model coefficients/tables can be extracted.

Examples

### calculate log risk ratios and corresponding sampling variances
dat <- escalc(measure="RR", ai=tpos, bi=tneg, ci=cpos, di=cneg, data=dat.bcg)

### fit mixed-effects model with absolute latitude and publication year as moderators
res <- rma(yi, vi, mods = ~ ablat + year, data=dat)

### extract model coefficients
coef(res)

### extract model coefficient table
coef(summary(res))

---

Confidence Intervals for 'rma' Objects

Description

Functions to compute confidence intervals for the model coefficients, variance components, and other parameters in meta-analytic models.

Usage

## S3 method for class 'rma.uni'
confint(object, parm, level, fixed=FALSE, random=TRUE, type,
        digits, transf, targs, verbose=FALSE, control, ...)

## S3 method for class 'rma.mh'
confint(object, parm, level, digits, transf, targs, ...)

## S3 method for class 'rma.peto'
confint(object, parm, level, digits, transf, targs, ...)

## S3 method for class 'rma.glmm'
confint(object, parm, level, digits, transf, targs, ...)

## S3 method for class 'rma.mv'
confint(object, parm, level, fixed=FALSE, sigma2, tau2, rho, gamma2, phi,
        digits, transf, targs, verbose=FALSE, control, ...)

## S3 method for class 'rma.uni.selmodel'
confint(object, parm, level, fixed=FALSE, tau2, delta,
        digits, transf, targs, verbose=FALSE, control, ...)

## S3 method for class 'rma.ls'
confint(object, parm, level, fixed=FALSE, alpha,
        digits, transf, targs, verbose=FALSE, control, ...)

Arguments

object	an object of class "rma.uni", "rma.mh", "rma.peto", "rma.mv", "rma.uni.selmodel", or "rma.ls". The method is not yet implemented for objects of class "rma.glmm".
parm	this argument is here for compatibility with the generic function confint, but is (currently) ignored.
fixed	logical to specify whether confidence intervals for the model coefficients should be returned.
random	logical to specify whether a confidence interval for the amount of (residual) heterogeneity should be returned.
type	optional character string to specify the method for computing the confidence interval for the amount of (residual) heterogeneity (either "QP", "GENQ", "PL", or "HT").
sigma2	integer to specify for which \(\sigma^2\) parameter a confidence interval should be obtained.
tau2	integer to specify for which \(\tau^2\) parameter a confidence interval should be obtained.
rho	integer to specify for which \(\rho\) parameter the confidence interval should be obtained.
gamma2	integer to specify for which \(\gamma^2\) parameter a confidence interval should be obtained.
phi	integer to specify for which \(\phi\) parameter a confidence interval should be obtained.
delta	integer to specify for which \(\delta\) parameter a confidence interval should be obtained.
alpha	integer to specify for which \(\alpha\) parameter a confidence interval should be obtained.
level	numeric value between 0 and 100 to specify the confidence interval level (see here for details). If unspecified, the default is to take the value from the object.
digits	optional integer to specify the number of decimal places to which the results should be rounded. If unspecified, the default is to take the value from the object.
transf	optional argument to specify a function to transform the model coefficients and interval bounds (e.g., transf=exp; see also transf). If unspecified, no transformation is used.
targs	optional arguments needed by the function specified under transf.
verbose	logical to specify whether output should be generated on the progress of the iterative algorithms used to obtain the confidence intervals (the default is FALSE). See ‘Details’.
control	list of control values for the iterative algorithms. If unspecified, default values are used. See ‘Note’.
...	other arguments.

Details

Confidence intervals for the model coefficients can be obtained by setting fixed=TRUE and are simply the usual Wald-type intervals (which are also shown when printing the fitted object).

Other parameter(s) for which confidence intervals can be obtained depend on the model object:

• For objects of class "rma.uni" obtained with the rma.uni function, a confidence interval for the amount of (residual) heterogeneity (i.e., \(\tau^2\)) can be obtained by setting random=TRUE (which is the default). The interval is obtained iteratively either via the Q-profile method or via the generalized Q-statistic method (Hartung and Knapp, 2005; Viechtbauer, 2007; Jackson, 2013; Jackson et al., 2014). The latter is automatically used when the model was fitted with method="GENQ" or method="GENQM", the former is used in all other cases. Either method provides an exact confidence interval for \(\tau^2\) in random- and mixed-effects models. The square root of the interval bounds is also returned for easier interpretation. Confidence intervals for \(I^2\) and \(H^2\) are also provided (Higgins & Thompson, 2002). Since \(I^2\) and \(H^2\) are just monotonic transformations of \(\tau^2\) (for details, see print), the confidence intervals for \(I^2\) and \(H^2\) are also exact. One can also set type="PL" to obtain a profile likelihood confidence interval for \(\tau^2\) (and corresponding CIs for \(I^2\) and \(H^2\)), which would be more consistent with the use of ML/REML estimation, but is not exact (see ‘Note’). For models without moderators (i.e., random-effects models), one can also set type="HT", in which case the ‘test-based method’ (method III in Higgins & Thompson, 2002) is used to construct confidence intervals for \(\tau^2\), \(I^2\), and \(H^2\) (see also Borenstein et al., 2009, chapter 16). However, note that this method tends to yield confidence intervals that are too narrow when the amount of heterogeneity is large.

• For objects of class "rma.mv" obtained with the rma.mv function, confidence intervals are obtained by default for all variance and correlation components of the model. Alternatively, one can use the sigma2, tau2, rho, gamma2, or phi arguments to specify for which variance/correlation parameter a confidence interval should be obtained. Only one of these arguments can be used at a time. A single integer is used to specify the number of the parameter. The function provides profile likelihood confidence intervals for these parameters. It is a good idea to examine the corresponding profile likelihood plots (via the profile function) to make sure that the bounds obtained are sensible.

• For selection model objects of class "rma.uni.selmodel" obtained with the selmodel function, confidence intervals are obtained by default for \(\tau^2\) (for models where this is an estimated parameter) and all selection model parameters. Alternatively, one can choose to obtain a confidence interval only for \(\tau^2\) by setting tau2=TRUE or for one of the selection model parameters by specifying its number via the delta argument. The function provides profile likelihood confidence intervals for these parameters. It is a good idea to examine the corresponding profile likelihood plots (via the profile function) to make sure that the bounds obtained are sensible.

• For location-scale model objects of class "rma.ls" obtained with the rma.uni function, confidence intervals are obtained by default for all scale parameters. Alternatively, one can choose to obtain a confidence interval for one of the scale parameters by specifying its number via the alpha argument. The function provides profile likelihood confidence intervals for these parameters. It is a good idea to examine the corresponding profile likelihood plots (via the profile function) to make sure that the bounds obtained are sensible.

The methods used to find confidence intervals for these parameters are iterative and require the use of the uniroot function. By default, the desired accuracy (tol) is set equal to .Machine$double.eps^0.25 and the maximum number of iterations (maxiter) to 1000. These values can be adjusted with control=list(tol=value, maxiter=value), but the defaults should be adequate for most purposes. If verbose=TRUE, output is generated on the progress of the iterative algorithms. This is especially useful when model fitting is slow, in which case finding the confidence interval bounds can also take considerable amounts of time.

When using the uniroot function, one must also set appropriate end points of the interval to be searched for the confidence interval bounds. The function sets some sensible defaults for the end points, but it may happen that the function is only able to determine that a bound is below/above a certain limit (this is indicated in the output accordingly with < or > signs). It can also happen that the model cannot be fitted or does not converge especially at the extremes of the interval to be searched. This will result in missing (NA) bounds and corresponding warnings. It may then be necessary to adjust the end points manually (see ‘Note’).

Finally, it is also possible that the lower and upper confidence interval bounds for a variance component both fall below zero. Since both bounds then fall outside of the parameter space, the confidence interval then consists of the null/empty set. Alternatively, one could interpret this as a confidence interval with bounds \([0,0]\) or as indicating ‘highly/overly homogeneous’ data.

Value

An object of class "confint.rma". The object is a list with either one or two elements (named fixed and random) with the following elements:

estimate	estimate of the model coefficient, variance/correlation component, or selection model parameter.
ci.lb	lower bound of the confidence interval.
ci.ub	upper bound of the confidence interval.

When obtaining confidence intervals for multiple components, the object is a list of class "list.confint.rma", where each element is a "confint.rma" object as described above.

The results are formatted and printed with the print function. To format the results as a data frame, one can use the as.data.frame function.

Note

When computing a CI for \(\tau^2\) for objects of class "rma.uni", the estimate of \(\tau^2\) will usually fall within the CI bounds provided by the Q-profile method. However, this is not guaranteed. Depending on the method used to estimate \(\tau^2\) and the width of the CI, it can happen that the CI does not actually contain the estimate. Using the empirical Bayes or Paule-Mandel estimator of \(\tau^2\) when fitting the model (i.e., using method="EB" or method="PM") usually ensures that the estimate of \(\tau^2\) falls within the CI (for method="PMM", this is guaranteed). When method="GENQ" was used to fit the model, the corresponding CI obtained via the generalized Q-statistic method also usually contains the estimate \(\tau^2\) (for method="GENQM", this is guaranteed). When using ML/REML estimation, the profile likelihood CI (obtained when setting type="PL") is guaranteed to contain the estimate of \(\tau^2\).

When computing a CI for \(\tau^2\) for objects of class "rma.uni", the end points of the interval to be searched for the CI bounds are \([0,100]\) (or, for the upper bound, ten times the estimate of \(\tau^2\), whichever is greater). The upper bound should be large enough for most cases, but can be adjusted with control=list(tau2.max=value). One can also adjust the lower end point with control=list(tau2.min=value). You should only play around with this value if you know what you are doing.

For objects of class "rma.mv", the function provides profile likelihood CIs for the variance/correlation parameters in the model. For variance components, the lower end point of the interval to be searched is set to 0 and the upper end point to the larger of 10 and 100 times the value of the component. For correlations, the function sets the lower end point to a sensible default depending on the type of variance structure chosen, while the upper end point is set to 1. One can adjust the lower and/or upper end points with control=list(vc.min=value, vc.max=value). Also, the function adjusts the lower/upper end points when the model does not converge at these extremes (the end points are then moved closer to the estimated value of the component). The total number of tries for setting/adjusting the end points in this manner is determined via control=list(eptries=value), with the default being 10 tries.

For objects of class "rma.uni.selmodel" or "rma.ls", the function also sets some sensible defaults for the end points of the interval to be searched for the CI bounds (of the \(\tau^2\), \(\delta\), and \(\alpha\) parameter(s)). One can again adjust the end points and the number of retries (as described above) with control=list(vc.min=value, vc.max=value, eptries=value).

The Q-profile and generalized Q-statistic methods are both exact under the assumptions of the random- and mixed-effects models (i.e., normally distributed observed and true effect sizes or outcomes and known sampling variances). In practice, these assumptions are usually only approximately true, turning CIs for \(\tau^2\) also into approximations. Profile likelihood CIs are not exact by construction and rely on the asymptotic behavior of the likelihood ratio statistic, so they may be inaccurate in small samples, but they are inherently consistent with the use of ML/REML estimation.

Author(s)

Wolfgang Viechtbauer [email protected] https://www.metafor-project.org

References

Borenstein, M., Hedges, L. V., Higgins, J. P. T., & Rothstein, H. (2009). Introduction to meta-analysis. Chichester, UK: Wiley.

Hardy, R. J., & Thompson, S. G. (1996). A likelihood approach to meta-analysis with random effects. Statistics in Medicine, 15(6), 619–629. ⁠https://doi.org/10.1002/(sici)1097-0258(19960330)15:6%3C619::aid-sim188%3E3.0.co;2-a⁠

Hartung, J., & Knapp, G. (2005). On confidence intervals for the among-group variance in the one-way random effects model with unequal error variances. Journal of Statistical Planning and Inference, 127(1-2), 157–177. ⁠https://doi.org/10.1016/j.jspi.2003.09.032⁠

Higgins, J. P. T., & Thompson, S. G. (2002). Quantifying heterogeneity in a meta-analysis. Statistics in Medicine, 21(11), 1539–1558. ⁠https://doi.org/10.1002/sim.1186⁠

Jackson, D. (2013). Confidence intervals for the between-study variance in random effects meta-analysis using generalised Cochran heterogeneity statistics. Research Synthesis Methods, 4(3), 220–229. ⁠https://doi.org/10.1186/s12874-016-0219-y⁠

Jackson, D., Turner, R., Rhodes, K., & Viechtbauer, W. (2014). Methods for calculating confidence and credible intervals for the residual between-study variance in random effects meta-regression models. BMC Medical Research Methodology, 14, 103. ⁠https://doi.org/10.1186/1471-2288-14-103⁠

Viechtbauer, W. (2007). Confidence intervals for the amount of heterogeneity in meta-analysis. Statistics in Medicine, 26(1), 37–52. ⁠https://doi.org/10.1002/sim.2514⁠

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., & López-López, J. A. (2022). Location-scale models for meta-analysis. Research Synthesis Methods. 13(6), 697–715. ⁠https://doi.org/10.1002/jrsm.1562⁠

See Also

rma.uni, rma.mh, rma.peto, rma.glmm, rma.mv, and selmodel for functions to fit models for which confidence intervals can be computed.

profile for functions to create profile likelihood plots corresponding to profile likelihood confidence intervals.

Examples

### calculate log risk ratios and corresponding sampling variances
dat <- escalc(measure="RR", ai=tpos, bi=tneg, ci=cpos, di=cneg, data=dat.bcg)

### meta-analysis of the log risk ratios using a random-effects model
res <- rma(yi, vi, data=dat, method="REML")

### confidence interval for the total amount of heterogeneity
confint(res)

### mixed-effects model with absolute latitude in the model
res <- rma(yi, vi, mods = ~ ablat, data=dat)

### confidence interval for the residual amount of heterogeneity
confint(res)

### multilevel random-effects model
res <- rma.mv(yi, vi, random = ~ 1 | district/school, data=dat.konstantopoulos2011)

### profile plots and confidence intervals for the variance components
## Not run: 
par(mfrow=c(2,1))
profile(res, sigma2=1, steps=40, cline=TRUE)
sav <- confint(res, sigma2=1)
sav
abline(v=sav$random[1,2:3], lty="dotted")
profile(res, sigma2=2, steps=40, cline=TRUE)
sav <- confint(res, sigma2=2)
sav
abline(v=sav$random[1,2:3], lty="dotted")

## End(Not run)

### multivariate parameterization of the model
res <- rma.mv(yi, vi, random = ~ school | district, data=dat.konstantopoulos2011)

### profile plots and confidence intervals for the variance component and correlation
## Not run: 
par(mfrow=c(2,1))
profile(res, tau2=1, steps=40, cline=TRUE)
sav <- confint(res, tau2=1)
sav
abline(v=sav$random[1,2:3], lty="dotted")
profile(res, rho=1, steps=40, cline=TRUE)
sav <- confint(res, rho=1)
sav
abline(v=sav$random[1,2:3], lty="dotted")

## End(Not run)

---

Construct Contrast Matrix for Two-Group Comparisons

Description

Function to construct a matrix that indicates which two groups have been contrasted against each other in each row of a dataset.

Usage

contrmat(data, grp1, grp2, last, shorten=FALSE, minlen=2, check=TRUE, append=TRUE)

Arguments

data	a data frame in wide format.
grp1	either the name (given as a character string) or the position (given as a single number) of the first group variable in the data frame.
grp2	either the name (given as a character string) or the position (given as a single number) of the second group variable in the data frame.
last	optional character string to specify which group will be placed in the last column of the matrix (must be one of the groups in the group variables). If not given, the most frequently occurring second group is placed last.
shorten	logical to specify whether the variable names corresponding to the group names should be shortened (the default is FALSE).
minlen	integer to specify the minimum length of the shortened variable names (the default is 2).
check	logical to specify whether the variables names should be checked to ensure that they are syntactically valid variable names and if not, they are adjusted (by make.names) so that they are (the default is TRUE).
append	logical to specify whether the contrast matrix should be appended to the data frame specified via the data argument (the default is TRUE). If append=FALSE, only the contrast matrix is returned.

Details

The function can be used to construct a matrix that indicates which two groups have been contrasted against each other in each row of a data frame (with 1 for the first group, -1 for the second group, and 0 otherwise).

The grp1 and grp2 arguments are used to specify the group variables in the dataset (either as character strings or as numbers indicating the column positions of these variables in the dataset). Optional argument last is used to specify which group will be placed in the last column of the matrix.

If shorten=TRUE, the variable names corresponding to the group names are shortened (to at least minlen; the actual length might be longer to ensure uniqueness of the variable names).

The examples below illustrate the use of this function.

Value

A matrix with as many variables as there are groups.

Author(s)

Wolfgang Viechtbauer [email protected] https://www.metafor-project.org

References

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⁠

See Also

to.wide for a function to create ‘wide’ format datasets.

dat.senn2013, dat.hasselblad1998, dat.lopez2019 for illustrative examples.

Examples

### restructure to wide format
dat <- dat.senn2013
dat <- dat[c(1,4,3,2,5,6)]
dat <- to.wide(dat, study="study", grp="treatment", ref="placebo", grpvars=4:6)
dat

### add contrast matrix
dat <- contrmat(dat, grp1="treatment.1", grp2="treatment.2")
dat

### data in long format
dat <- dat.hasselblad1998
dat

### restructure to wide format
dat <- to.wide(dat, study="study", grp="trt", ref="no_contact", grpvars=6:7)
dat

### add contrast matrix
dat <- contrmat(dat, grp1="trt.1", grp2="trt.2", shorten=TRUE, minlen=3)
dat

---

Reconstruct Cell Frequencies of \(2 \times 2\) Tables

Description

Function to reconstruct the cell frequencies of \(2 \times 2\) tables based on other summary statistics.

Usage

conv.2x2(ori, ri, x2i, ni, n1i, n2i, correct=TRUE, data, include,
         var.names=c("ai","bi","ci","di"), append=TRUE, replace="ifna")

Arguments

ori	optional vector with the odds ratios corresponding to the tables.
ri	optional vector with the phi coefficients corresponding to the tables.
x2i	optional vector with the (signed) chi-square statistics corresponding to the tables.
ni	vector with the total sample sizes.
n1i	vector with the marginal counts for the outcome of interest on the first variable.
n2i	vector with the marginal counts for the outcome of interest on the second variable.
correct	optional logical (or vector thereof) to specify whether chi-square statistics were computed using Yates's correction for continuity (the default is TRUE).
data	optional data frame containing the variables given to the arguments above.
include	optional (logical or numeric) vector to specify the subset of studies for which the cell frequencies should be reconstructed.
var.names	character vector with four elements to specify the names of the variables for the reconstructed cell frequencies (the default is c("ai","bi","ci","di")).
append	logical to specify whether the data frame provided via the data argument should be returned together with the reconstructed values (the default is TRUE).
replace	character string or logical to specify how values in var.names should be replaced (only relevant when using the data argument and if variables in var.names already exist in the data frame). See the ‘Value’ section for more details.

Details

For meta-analyses based on \(2 \times 2\) table data, the problem often arises that some studies do not directly report the cell frequencies. The present function allows the reconstruction of such tables based on other summary statistics.

In particular, assume that the data of interest for a particular study are of the form:

		variable 2, outcome +		variable 2, outcome -		total
variable 1, outcome +		ai		bi		n1i
variable 1, outcome -		ci		di
total		n2i				ni

where ai, bi, ci, and di denote the cell frequencies (i.e., the number of individuals falling into a particular category), n1i (i.e., ai+bi) and n2i (i.e., ai+ci) are the marginal totals for the outcome of interest on the first and second variable, respectively, and ni is the total sample size (i.e., ai+bi+ci+di) of the study.

For example, if variable 1 denotes two different groups (e.g., treated versus control) and variable 2 indicates whether a particular outcome of interest has occurred or not (e.g., death, complications, failure to improve under the treatment), then n1i denotes the number of individuals in the treatment group, but n2i is not the number of individuals in the control group, but the total number of individuals who experienced the outcome of interest on variable 2. Note that the meaning of n2i is therefore different here compared to the escalc function (where n2i denotes ci+di).

If a study does not report the cell frequencies, but it reports the total sample size (which can be specified via the ni argument), the two marginal counts (which can be specified via the n1i and n2i arguments), and some other statistic corresponding to the table, then it may be possible to reconstruct the cell frequencies. The present function currently allows this for three different cases:

1. If the odds ratio \[OR = \frac{a_i d_i}{b_i c_i}\] is known, then the cell frequencies can be reconstructed (Bonett, 2007). Odds ratios can be specified via the ori argument.

2. If the phi coefficient \[\phi = \frac{a_i d_i - b_i c_i}{\sqrt{n_{1i}(n_i-n_{1i})n_{2i}(n_i-n_{2i})}}\] is known, then the cell frequencies can again be reconstructed (own derivation). Phi coefficients can be specified via the ri argument.

3. If the chi-square statistic from Pearson's chi-square test of independence is known (which can be specified via the x2i argument), then it can be used to recalculate the phi coefficient and hence again the cell frequencies can be reconstructed. However, the chi-square statistic does not carry information about the sign of the phi coefficient. Therefore, values specified via the x2i argument can be positive or negative, which allows the specification of the correct sign. Also, when using a chi-square statistic as input, it is assumed that it was computed using Yates's correction for continuity (unless correct=FALSE). If the chi-square statistic is not known, but its p-value, one can first back-calculate the chi-square statistic using qchisq(<p-value>, df=1, lower.tail=FALSE).

Typically, the odds ratio, phi coefficient, or chi-square statistic (or its p-value) that can be extracted from a study will be rounded to a certain degree. The calculations underlying the function are exact only for unrounded values. Rounding can therefore introduce some discrepancies.

If a marginal total is unknown, then external information needs to be used to ‘guestimate’ the number of individuals that experienced the outcome of interest on this variable. Depending on the accuracy of such an estimate, the reconstructed cell frequencies will be more or less accurate and need to be treated with due caution.

The true marginal counts also put constraints on the possible values for the odds ratio, phi coefficient, and chi-square statistic. If a marginal count is replaced by a guestimate which is not compatible with the given statistic, one or more reconstructed cell frequencies may be negative. The function issues a warning if this happens and sets the cell frequencies to NA for such a study.

If only one of the two marginal counts is unknown but a 95% CI for the odds ratio is also available, then the estimraw package can also be used to reconstruct the corresponding cell frequencies (Di Pietrantonj, 2006; but see Veroniki et al., 2013, for some cautions).

Value

If the data argument was not specified or append=FALSE, a data frame with four variables called var.names with the reconstructed cell frequencies.

If data was specified and append=TRUE, then the original data frame is returned. If var.names[j] (for \(\textrm{j} \in \{1, \ldots, 4\}\)) is a variable in data and replace="ifna" (or replace=FALSE), then only missing values in this variable are replaced with the estimated frequencies (where possible) and otherwise a new variable called var.names[j] is added to the data frame.

If replace="all" (or replace=TRUE), then all values in var.names[j] where a reconstructed cell frequency can be computed are replaced, even for cases where the value in var.names[j] is not missing.

Author(s)

Wolfgang Viechtbauer [email protected] https://www.metafor-project.org

References

Bonett, D. G. (2007). Transforming odds ratios into correlations for meta-analytic research. American Psychologist, 62(3), 254–255. ⁠https://doi.org/10.1037/0003-066x.62.3.254⁠

Di Pietrantonj, C. (2006). Four-fold table cell frequencies imputation in meta analysis. Statistics in Medicine, 25(13), 2299–2322. ⁠https://doi.org/10.1002/sim.2287⁠

Veroniki, A. A., Pavlides, M., Patsopoulos, N. A., & Salanti, G. (2013). Reconstructing 2 x 2 contingency tables from odds ratios using the Di Pietrantonj method: Difficulties, constraints and impact in meta-analysis results. Research Synthesis Methods, 4(1), 78–94. ⁠https://doi.org/10.1002/jrsm.1061⁠

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⁠

See Also

escalc for a function to compute various effect size measures based on \(2 \times 2\) table data.

Examples

### demonstration that the reconstruction of the 2x2 table works
### (note: the values in rows 2, 3, and 4 correspond to those in row 1)
dat <- data.frame(ai=c(36,NA,NA,NA), bi=c(86,NA,NA,NA), ci=c(20,NA,NA,NA), di=c(98,NA,NA,NA),
                  oddsratio=NA, phi=NA, chisq=NA, ni=NA, n1i=NA, n2i=NA)
dat$oddsratio[2] <- round(exp(escalc(measure="OR", ai=ai, bi=bi, ci=ci, di=di, data=dat)$yi[1]), 2)
dat$phi[3] <- round(escalc(measure="PHI", ai=ai, bi=bi, ci=ci, di=di, data=dat)$yi[1], 2)
dat$chisq[4] <- round(chisq.test(matrix(c(t(dat[1,1:4])), nrow=2, byrow=TRUE))$statistic, 2)
dat$ni[2:4]  <- with(dat, ai[1] + bi[1] + ci[1] + di[1])
dat$n1i[2:4] <- with(dat, ai[1] + bi[1])
dat$n2i[2:4] <- with(dat, ai[1] + ci[1])
dat

### reconstruct cell frequencies for rows 2, 3, and 4
dat <- conv.2x2(ri=phi, ori=oddsratio, x2i=chisq, ni=ni, n1i=n1i, n2i=n2i, data=dat)
dat

### same example but with cell frequencies that are 10 times as large
dat <- data.frame(ai=c(360,NA,NA,NA), bi=c(860,NA,NA,NA), ci=c(200,NA,NA,NA), di=c(980,NA,NA,NA),
                  oddsratio=NA, phi=NA, chisq=NA, ni=NA, n1i=NA, n2i=NA)
dat$oddsratio[2] <- round(exp(escalc(measure="OR", ai=ai, bi=bi, ci=ci, di=di, data=dat)$yi[1]), 2)
dat$phi[3] <- round(escalc(measure="PHI", ai=ai, bi=bi, ci=ci, di=di, data=dat)$yi[1], 2)
dat$chisq[4] <- round(chisq.test(matrix(c(t(dat[1,1:4])), nrow=2, byrow=TRUE))$statistic, 2)
dat$ni[2:4]  <- with(dat, ai[1] + bi[1] + ci[1] + di[1])
dat$n1i[2:4] <- with(dat, ai[1] + bi[1])
dat$n2i[2:4] <- with(dat, ai[1] + ci[1])
dat <- conv.2x2(ri=phi, ori=oddsratio, x2i=chisq, ni=ni, n1i=n1i, n2i=n2i, data=dat)
dat # slight inaccuracy in row 3 due to rounding

### demonstrate what happens when a true marginal count is guestimated
escalc(measure="PHI", ai=176, bi=24, ci=72, di=128)
conv.2x2(ri=0.54, ni=400, n1i=200, n2i=248) # using the true marginal counts
conv.2x2(ri=0.54, ni=400, n1i=200, n2i=200) # marginal count for variable 2 is guestimated
conv.2x2(ri=0.54, ni=400, n1i=200, n2i=50)  # marginal count for variable 2 is incompatible

### demonstrate that using the correct sign for the chi-square statistic is important
chisq <- round(chisq.test(matrix(c(40,60,60,40), nrow=2, byrow=TRUE))$statistic, 2)
conv.2x2(x2i=-chisq, ni=200, n1i=100, n2i=100) # correct reconstruction
conv.2x2(x2i=chisq, ni=200, n1i=100, n2i=100) # incorrect reconstruction

### demonstrate use of the 'correct' argument
tab <- matrix(c(28,14,12,18), nrow=2, byrow=TRUE)
chisq <- round(chisq.test(tab)$statistic, 2) # chi-square test with Yates' continuity correction
conv.2x2(x2i=chisq, ni=72, n1i=42, n2i=40) # correct reconstruction
chisq <- round(chisq.test(tab, correct=FALSE)$statistic, 2) # without Yates' continuity correction
conv.2x2(x2i=chisq, ni=72, n1i=42, n2i=40) # incorrect reconstruction
conv.2x2(x2i=chisq, ni=72, n1i=42, n2i=40, correct=FALSE) # correct reconstruction

### recalculate chi-square statistic based on p-value
pval <- round(chisq.test(tab)$p.value, 2)
chisq <- qchisq(pval, df=1, lower.tail=FALSE)
conv.2x2(x2i=chisq, ni=72, n1i=42, n2i=40)

---

Transform Observed Effect Sizes or Outcomes and their Sampling Variances using the Delta Method

Description

Function to transform observed effect sizes or outcomes and their sampling variances using the delta method.

Usage

conv.delta(yi, vi, ni, data, include, transf, var.names, append=TRUE, replace="ifna", ...)

Arguments

yi	vector with the observed effect sizes or outcomes.
vi	vector with the corresponding sampling variances.
ni	vector with the total sample sizes of the studies.
data	optional data frame containing the variables given to the arguments above.
include	optional (logical or numeric) vector to specify the subset of studies for which the transformation should be carried out.
transf	a function which should be used for the transformation.
var.names	character vector with two elements to specify the name of the variable for the transformed effect sizes or outcomes and the name of the variable for the corresponding sampling variances (if data is an object of class "escalc", the var.names are taken from the object; otherwise the defaults are "yi" and "vi").
append	logical to specify whether the data frame provided via the data argument should be returned together with the estimated values (the default is TRUE).
replace	character string or logical to specify how values in var.names should be replaced (only relevant when using the data argument and if variables in var.names already exist in the data frame). See the ‘Value’ section for more details.
...	other arguments for the transformation function.

Details

The escalc function can be used to compute a wide variety of effect sizes or ‘outcome measures’. In some cases, it may be necessary to transform one type of measure to another. The present function provides a general method for doing so via the delta method, which briefly works as follows.

Let \(y_i\) denote the observed effect size or outcome for a particular study and \(v_i\) the corresponding sampling variance. Then \(f(y_i)\) will be the transformed effect size or outcome, where \(f(\cdot)\) is the function specified via the transf argument. The sampling variance of the transformed effect size or outcome is then computed with \(v_i \times f'(y_i)^2\), where \(f'(y_i)\) denotes the derivative of \(f(\cdot)\) evaluated at \(y_i\). The present function computes the derivative numerically using the grad function from the numDeriv package.

The value of the observed effect size or outcome should be the first argument of the function specified via transf. The function can have additional arguments, which can be specified via the ... argument. However, due to the manner in which these additional arguments are evaluated, they cannot have names that match one of the arguments of the grad function (an error will be issued if such a naming clash is detected).

Optionally, one can use the ni argument to supply the total sample sizes of the studies. This has no relevance for the calculations done by the present function, but some other functions may use this information (e.g., when drawing a funnel plot with the funnel function and one adjusts the yaxis argument to one of the options that puts the sample sizes or some transformation thereof on the y-axis).

Value

If the data argument was not specified or append=FALSE, a data frame of class c("escalc","data.frame") with two variables called var.names[1] (by default "yi") and var.names[2] (by default "vi") with the transformed observed effect sizes or outcomes and the corresponding sampling variances (computed as described above).

If data was specified and append=TRUE, then the original data frame is returned. If var.names[1] is a variable in data and replace="ifna" (or replace=FALSE), then only missing values in this variable are replaced with the transformed observed effect sizes or outcomes (where possible) and otherwise a new variable called var.names[1] is added to the data frame. Similarly, if var.names[2] is a variable in data and replace="ifna" (or replace=FALSE), then only missing values in this variable are replaced with the sampling variances calculated as described above (where possible) and otherwise a new variable called var.names[2] is added to the data frame.

If replace="all" (or replace=TRUE), then all values in var.names[1] and var.names[2] are replaced, even for cases where the value in var.names[1] and var.names[2] is not missing.

Author(s)

Wolfgang Viechtbauer [email protected] https://www.metafor-project.org

References

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⁠

See Also

escalc for a function to compute various effect size measures.

Examples

############################################################################

### the following examples illustrate that the use of the delta method (with numeric derivatives)
### yields essentially identical results as the analytic calculations that are done by escalc()

### compute logit transformed proportions and corresponding sampling variances for two studies
escalc(measure="PLO", xi=c(5,12), ni=c(40,80))

### compute raw proportions and corresponding sampling variances for the two studies
dat <- escalc(measure="PR", xi=c(5,12), ni=c(40,80))
dat

### apply the logit transformation (note: this yields the same values as above with measure="PLO")
conv.delta(dat$yi, dat$vi, transf=transf.logit)

### using the 'data' argument
conv.delta(yi, vi, data=dat, transf=transf.logit, var.names=c("yi.t","vi.t"))

### or replace the existing 'yi' and 'vi' values
conv.delta(yi, vi, data=dat, transf=transf.logit, replace="all")

######################################

### use escalc() with measure D2ORN which transforms standardized mean differences (computed
### from means and standard deviations) into the corresponding log odds ratios
escalc(measure="D2ORN", m1i=m1i, sd1i=sd1i, n1i=n1i,
                        m2i=m2i, sd2i=sd2i, n2i=n2i, data=dat.normand1999)

### use escalc() to compute standardized mean differences (without the usual bias correction) and
### then apply the same transformation to the standardized mean differences
dat <- escalc(measure="SMD", m1i=m1i, sd1i=sd1i, n1i=n1i,
                             m2i=m2i, sd2i=sd2i, n2i=n2i, data=dat.normand1999, correct=FALSE)
conv.delta(yi, vi, data=dat, transf=transf.dtolnor.norm, replace="all")

######################################

### an example where the transformation function takes additional arguments

### use escalc() with measure RPB which transforms standardized mean differences (computed
### from means and standard deviations) into the corresponding point-biserial correlations
escalc(measure="RPB", m1i=m1i, sd1i=sd1i, n1i=n1i,
                      m2i=m2i, sd2i=sd2i, n2i=n2i, data=dat.normand1999)

### use escalc() to compute standardized mean differences (without the usual bias correction) and
### then apply the same transformation to the standardized mean differences
dat <- escalc(measure="SMD", m1i=m1i, sd1i=sd1i, n1i=n1i,
                             m2i=m2i, sd2i=sd2i, n2i=n2i, data=dat.normand1999, correct=FALSE)
conv.delta(yi, vi, data=dat, transf=transf.dtorpb, n1i=n1i, n2i=n2i, replace="all")

############################################################################

### a more elaborate example showing how this function could be used in the data
### preparation steps for a meta-analysis of standardized mean differences (SMDs)

dat <- data.frame(study=1:6,
         m1i=c(2.03,NA,NA,NA,NA,NA), sd1i=c(0.95,NA,NA,NA,NA,NA), n1i=c(32,95,145,NA,NA,NA),
         m2i=c(1.25,NA,NA,NA,NA,NA), sd2i=c(1.04,NA,NA,NA,NA,NA), n2i=c(30,99,155,NA,NA,NA),
         tval=c(NA,2.12,NA,NA,NA,NA), dval=c(NA,NA,0.37,NA,NA,NA),
         ai=c(NA,NA,NA,26,NA,NA), bi=c(NA,NA,NA,58,NA,NA),
         ci=c(NA,NA,NA,11,NA,NA), di=c(NA,NA,NA,74,NA,NA),
         or=c(NA,NA,NA,NA,2.56,NA), lower=c(NA,NA,NA,NA,1.23,NA), upper=c(NA,NA,NA,NA,5.30,NA),
         corr=c(NA,NA,NA,NA,NA,.32), ntot=c(NA,NA,NA,NA,NA,86))
dat

### study types:
### 1) reports means and SDs so that the SMD can be directly calculated
### 2) reports the t-statistic from an independent samples t-test (and group sizes)
### 3) reports the standardized mean difference directly (and group sizes)
### 4) dichotomized the continuous dependent variable and reports the resulting 2x2 table
### 5) dichotomized the continuous dependent variable and reports an odds ratio with 95% CI
### 6) treated the group variable continuously and reports a Pearson product-moment correlation

### use escalc() to directly compute the SMD and its variance for studies 1, 2, and 3
dat <- escalc(measure="SMD", m1i=m1i, sd1i=sd1i, n1i=n1i,
                             m2i=m2i, sd2i=sd2i, n2i=n2i, ti=tval, di=dval, data=dat)
dat

### use escalc() with measure OR2DN to compute the SMD value for study 4
dat <- escalc(measure="OR2DN", ai=ai, bi=bi, ci=ci, di=di, data=dat, replace=FALSE)
dat

### use conv.wald() to convert the OR and CI into the log odds ratio and its variance for study 5
dat <- conv.wald(out=or, ci.lb=lower, ci.ub=upper, data=dat,
                 transf=log, var.names=c("lnor","vlnor"))
dat

### use conv.delta() to transform the log odds ratio into the SMD value for study 5
dat <- conv.delta(lnor, vlnor, data=dat,
                  transf=transf.lnortod.norm, var.names=c("yi","vi"))
dat

### remove the lnor and vlnor variables (no longer needed)
dat$lnor  <- NULL
dat$vlnor <- NULL

### use escalc() with measure COR to compute the sampling variance of ri for study 6
dat <- escalc(measure="COR", ri=corr, ni=ntot, data=dat, var.names=c("ri","vri"))
dat

### use conv.delta() to transform the correlation into the SMD value for study 6
dat <- conv.delta(ri, vri, data=dat, transf=transf.rtod, var.names=c("yi","vi"))
dat

### remove the ri and vri variables (no longer needed)
dat$ri  <- NULL
dat$vri <- NULL

### now variable 'yi' is complete with the SMD values for all studies
dat

### fit an equal-effects model to the SMD values
rma(yi, vi, data=dat, method="EE")

############################################################################

### a more elaborate example showing how this function could be used in the data
### preparation steps for a meta-analysis of correlation coefficients

dat <- data.frame(study=1:6,
         ri=c(.42,NA,NA,NA,NA,NA),
         tval=c(NA,2.85,NA,NA,NA,NA),
         phi=c(NA,NA,NA,0.27,NA,NA),
         ni=c(93,182,NA,112,NA,NA),
         ai=c(NA,NA,NA,NA,61,NA), bi=c(NA,NA,NA,NA,36,NA),
         ci=c(NA,NA,NA,NA,39,NA), di=c(NA,NA,NA,NA,57,NA),
         or=c(NA,NA,NA,NA,NA,1.86), lower=c(NA,NA,NA,NA,NA,1.12), upper=c(NA,NA,NA,NA,NA,3.10),
         m1i=c(NA,NA,54.1,NA,NA,NA), sd1i=c(NA,NA,5.79,NA,NA,NA), n1i=c(NA,NA,66,75,NA,NA),
         m2i=c(NA,NA,51.7,NA,NA,NA), sd2i=c(NA,NA,6.23,NA,NA,NA), n2i=c(NA,NA,65,88,NA,NA))
dat

### study types:
### 1) reports the correlation coefficient directly
### 2) reports the t-statistic from a t-test of H0: rho = 0
### 3) dichotomized one variable and reports means and SDs for the two corresponding groups
### 4) reports the phi coefficient, marginal counts, and total sample size
### 5) dichotomized both variables and reports the resulting 2x2 table
### 6) dichotomized both variables and reports an odds ratio with 95% CI

### use escalc() to directly compute the correlation and its variance for studies 1 and 2
dat <- escalc(measure="COR", ri=ri, ni=ni, ti=tval, data=dat)
dat

### use escalc() with measure RBIS to compute the biserial correlation for study 3
dat <- escalc(measure="RBIS", m1i=m1i, sd1i=sd1i, n1i=n1i,
                              m2i=m2i, sd2i=sd2i, n2i=n2i, data=dat, replace=FALSE)
dat

### use conv.2x2() to reconstruct the 2x2 table for study 4
dat <- conv.2x2(ri=phi, ni=ni, n1i=n1i, n2i=n2i, data=dat)
dat

### use escalc() with measure RTET to compute the tetrachoric correlation for studies 4 and 5
dat <- escalc(measure="RTET", ai=ai, bi=bi, ci=ci, di=di, data=dat, replace=FALSE)
dat

### use conv.wald() to convert the OR and CI into the log odds ratio and its variance for study 6
dat <- conv.wald(out=or, ci.lb=lower, ci.ub=upper, data=dat,
                 transf=log, var.names=c("lnor","vlnor"))
dat

### use conv.delta() to estimate the tetrachoric correlation from the log odds ratio for study 6
dat <- conv.delta(lnor, vlnor, data=dat,
                  transf=transf.lnortortet.pearson, var.names=c("yi","vi"))
dat

### remove the lnor and vlnor variables (no longer needed)
dat$lnor  <- NULL
dat$vlnor <- NULL

### now variable 'yi' is complete with the correlations for all studies
dat

### fit an equal-effects model to the correlations
rma(yi, vi, data=dat, method="EE")

############################################################################

---

Estimate Means and Standard Deviations from Five-Number Summary Values

Description

Function to estimate means and standard deviations from five-number summary values.

Usage

conv.fivenum(min, q1, median, q3, max, n, data, include,
             method="default", dist="norm", transf=TRUE, test=TRUE,
             var.names=c("mean","sd"), append=TRUE, replace="ifna", ...)

Arguments

min	vector with the minimum values.
q1	vector with the lower/first quartile values.
median	vector with the median values.
q3	vector with the upper/third quartile values.
max	vector with the maximum values.
n	vector with the sample sizes.
data	optional data frame containing the variables given to the arguments above.
include	optional (logical or numeric) vector to specify the subset of studies for which means and standard deviations should be estimated.
method	character string to specify the method to use. Either "default" (same as "luo/wan/shi" which is the current default), "qe", "bc", "mln", or "blue". Can be abbreviated. See ‘Details’.
dist	character string to specify the assumed distribution for the underlying data (either "norm" for a normal distribution or "lnorm" for a log-normal distribution). Can also be a string vector if different distributions are assumed for different studies. Only relevant when method="default".
transf	logical to specify whether the estimated means and standard deviations of the log-transformed data should be back-transformed as described by Shi et al. (2020b) (the default is TRUE). Only relevant when dist="lnorm" and when method="default".
test	logical to specify whether a study should be excluded from the estimation if the test for skewness is significant (the default is TRUE, but whether this is applicable depends on the method; see ‘Details’).
var.names	character vector with two elements to specify the name of the variable for the estimated means and the name of the variable for the estimated standard deviations (the defaults are "mean" and "sd").
append	logical to specify whether the data frame provided via the data argument should be returned together with the estimated values (the default is TRUE).
replace	character string or logical to specify how values in var.names should be replaced (only relevant when using the data argument and if variables in var.names already exist in the data frame). See the ‘Value’ section for more details.
...	other arguments.

Details

Various effect size measures require means and standard deviations (SDs) as input (e.g., raw or standardized mean differences, ratios of means / response ratios; see escalc for further details). For some studies, authors may not report means and SDs, but other statistics, such as the so-called ‘five-number summary’, consisting of the minimum, lower/first quartile, median, upper/third quartile, and the maximum of the sample values (plus the sample sizes). Occasionally, only a subset of these values are reported.

The present function can be used to estimate means and standard deviations from five-number summary values based on various methods described in the literature (Bland, 2015; Cai et al. 2021; Hozo et al., 2005; Luo et al., 2016; McGrath et al., 2020; Shi et al., 2020a; Walter & Yao, 2007; Wan et al., 2014; Yang et al., 2022).

When method="default" (which is the same as "luo/wan/shi"), the following methods are used:

Case 1: Min, Median, Max

In case only the minimum, median, and maximum is available for a study (plus the sample size), then the function uses the method by Luo et al. (2016), equation (7), to estimate the mean and the method by Wan et al. (2014), equation (9), to estimate the SD.

Case 2: Q1, Median, Q3

In case only the lower/first quartile, median, and upper/third quartile is available for a study (plus the sample size), then the function uses the method by Luo et al. (2016), equation (11), to estimate the mean and the method by Wan et al. (2014), equation (16), to estimate the SD.

Case 3: Min, Q1, Median, Q3, Max

In case the full five-number summary is available for a study (plus the sample size), then the function uses the method by Luo et al. (2016), equation (15), to estimate the mean and the method by Shi et al. (2020a), equation (10), to estimate the SD.

———

The median is not actually needed in the methods by Wan et al. (2014) and Shi et al. (2020a) and hence it is possible to estimate the SD even if the median is unavailable (this can be useful if a study reports the mean directly, but instead of the SD, it reports the minimum/maximum and/or first/third quartile values).

Note that the sample size must be at least 5 to apply these methods. Studies where the sample size is smaller are not included in the estimation. The function also checks that min <= q1 <= median <= q3 <= max and throws an error if any studies are found where this is not the case.

Test for Skewness

The methods described above were derived under the assumption that the data are normally distributed. Testing this assumption would require access to the raw data, but based on the three cases above, Shi et al. (2023) derived tests for skewness that only require the reported quantile values and the sample sizes. These tests are automatically carried out. When test=TRUE (which is the default), a study is automatically excluded from the estimation if the test is significant. If all studies should be included, set test=FALSE, but note that the accuracy of the methods will tend to be poorer when the data come from an apparently skewed (and hence non-normal) distribution.

Log-Normal Distribution

When setting dist="lnorm", the raw data are assumed to follow a log-normal distribution. In this case, the methods as described by Shi et al. (2020b) are used to estimate the mean and SD of the log transformed data for the three cases above. When transf=TRUE (the default), the estimated mean and SD of the log transformed data are back-transformed to the estimated mean and SD of the raw data (using the bias-corrected back-transformation as described by Shi et al., 2020b). Note that the test for skewness is also carried out when dist="lnorm", but now testing if the log transformed data exhibit skewness.

Alternative Methods

As an alternative to the methods above, one can make use of the methods implemented in the estmeansd package to estimate means and SDs based on the three cases above. Available are the quantile estimation method (method="qe"; using the qe.mean.sd function; McGrath et al., 2020), the Box-Cox method (method="bc"; using the bc.mean.sd function; McGrath et al., 2020), and the method for unknown non-normal distributions (method="mln"; using the mln.mean.sd function; Cai et al. 2021). The advantage of these methods is that they do not assume that the data underlying the reported values are normally distributed (and hence the test argument is ignored), but they can only be used when the values are positive (except for the quantile estimation method, which can also be used when one or more of the values are negative, but in this case the method does assume that the data are normally distributed and hence the test for skewness is applied when test=TRUE). Note that all of these methods may struggle to provide sensible estimates when some of the values are equal to each other (which can happen when the data include a lot of ties and/or the reported values are rounded). Also, the Box-Cox method and the method for unknown non-normal distributions involve simulated data and hence results will slightly change on repeated runs. Setting the seed of the random number generator (with set.seed) ensures reproducibility.

Finally, by setting method="blue", one can make use of the BLUE_s function from the metaBLUE package to estimate means and SDs based on the three cases above (Yang et al., 2022). The method assumes that the underlying data are normally distributed (and hence the test for skewness is applied when test=TRUE).

Value

If the data argument was not specified or append=FALSE, a data frame with two variables called var.names[1] (by default "mean") and var.names[2] (by default "sd") with the estimated means and SDs.

If data was specified and append=TRUE, then the original data frame is returned. If var.names[1] is a variable in data and replace="ifna" (or replace=FALSE), then only missing values in this variable are replaced with the estimated means (where possible) and otherwise a new variable called var.names[1] is added to the data frame. Similarly, if var.names[2] is a variable in data and replace="ifna" (or replace=FALSE), then only missing values in this variable are replaced with the estimated SDs (where possible) and otherwise a new variable called var.names[2] is added to the data frame.

If replace="all" (or replace=TRUE), then all values in var.names[1] and var.names[2] where an estimated mean and SD can be computed are replaced, even for cases where the value in var.names[1] and var.names[2] is not missing.

When missing values in var.names[1] are replaced, an attribute called "est" is added to the variable, which is a logical vector that is TRUE for values that were estimated. The same is done when missing values in var.names[2] are replaced.

Attributes called "tval", "crit", "sig", and "dist" are also added to var.names[1] corresponding to the test statistic and critical value for the test for skewness, whether the test was significant, and the assumed distribution (for the quantile estimation method, this is the distribution that provides the best fit to the given values).

Note

A word of caution: Under the given distributional assumptions, the estimated means and SDs are approximately unbiased and hence so are any effect size measures computed based on them (assuming a measure is unbiased to begin with when computed with directly reported means and SDs). However, the estimated means and SDs are less precise (i.e., are more variable) than directly reported means and SDs (especially under case 1) and hence computing the sampling variance of a measure with equations that assume that directly reported means and SDs are available will tend to underestimate the actual sampling variance of the measure, giving too much weight to estimates computed based on estimated means and SDs (see also McGrath et al., 2023). It would therefore be prudent to treat effect size estimates computed from estimated means and SDs with caution (e.g., by examining in a moderator analysis whether there are systematic differences between studies directly reporting means and SDs and those where the means and SDs needed to be estimated and/or as part of a sensitivity analysis). McGrath et al. (2023) also suggest to use bootstrapping to estimate the sampling variance of effect size measures computed based on estimated means and SDs. See also the metamedian package for this purpose.

Also note that the development of methods for estimating means and SDs based on five-number summary values is an active area of research. Currently, when method="default", then this is identical to method="luo/wan/shi", but this might change in the future. For reproducibility, it is therefore recommended to explicitly set method="luo/wan/shi" (or one of the other methods) when running this function.

Author(s)

Wolfgang Viechtbauer [email protected] https://www.metafor-project.org

References

Bland, M. (2015). Estimating mean and standard deviation from the sample size, three quartiles, minimum, and maximum. International Journal of Statistics in Medical Research, 4(1), 57–64. ⁠https://doi.org/10.6000/1929-6029.2015.04.01.6⁠

Cai, S., Zhou, J., & Pan, J. (2021). Estimating the sample mean and standard deviation from order statistics and sample size in meta-analysis. Statistical Methods in Medical Research, 30(12), 2701–2719. ⁠https://doi.org/10.1177/09622802211047348⁠

Hozo, S. P., Djulbegovic, B. & Hozo, I. (2005). Estimating the mean and variance from the median, range, and the size of a sample. BMC Medical Research Methodology, 5, 13. ⁠https://doi.org/10.1186/1471-2288-5-13⁠

Luo, D., Wan, X., Liu, J. & Tong, T. (2016). Optimally estimating the sample mean from the sample size, median, mid-range, and/or mid-quartile range. Statistical Methods in Medical Research, 27(6), 1785–1805. ⁠https://doi.org/10.1177/0962280216669183⁠

McGrath, S., Zhao, X., Steele, R., Thombs, B. D., Benedetti, A., & the DEPRESsion Screening Data (DEPRESSD) Collaboration (2020). Estimating the sample mean and standard deviation from commonly reported quantiles in meta-analysis. Statistical Methods in Medical Research, 29(9), 2520–2537. ⁠https://doi.org/10.1177/0962280219889080⁠

McGrath, S., Katzenschlager, S., Zimmer, A. J., Seitel, A., Steele, R., & Benedetti, A. (2023). Standard error estimation in meta-analysis of studies reporting medians. Statistical Methods in Medical Research, 32(2), 373–388. ⁠https://doi.org/10.1177/09622802221139233⁠

Shi, J., Luo, D., Weng, H., Zeng, X.-T., Lin, L., Chu, H. & Tong, T. (2020a). Optimally estimating the sample standard deviation from the five-number summary. Research Synthesis Methods, 11(5), 641–654. ⁠https://doi.org/https://doi.org/10.1002/jrsm.1429⁠

Shi, J., Tong, T., Wang, Y. & Genton, M. G. (2020b). Estimating the mean and variance from the five-number summary of a log-normal distribution. Statistics and Its Interface, 13(4), 519–531. https://doi.org/10.4310/sii.2020.v13.n4.a9

Shi, J., Luo, D., Wan, X., Liu, Y., Liu, J., Bian, Z. & Tong, T. (2023). Detecting the skewness of data from the five-number summary and its application in meta-analysis. Statistical Methods in Medical Research, 32(7), 1338–1360. ⁠https://doi.org/10.1177/09622802231172043⁠

Walter, S. D. & Yao, X. (2007). Effect sizes can be calculated for studies reporting ranges for outcome variables in systematic reviews. Journal of Clinical Epidemiology, 60(8), 849-852. ⁠https://doi.org/10.1016/j.jclinepi.2006.11.003⁠

Wan, X., Wang, W., Liu, J. & Tong, T. (2014). Estimating the sample mean and standard deviation from the sample size, median, range and/or interquartile range. BMC Medical Research Methodology, 14, 135. ⁠https://doi.org/10.1186/1471-2288-14-135⁠

Yang, X., Hutson, A. D., & Wang, D. (2022). A generalized BLUE approach for combining location and scale information in a meta-analysis. Journal of Applied Statistics, 49(15), 3846–3867. ⁠https://doi.org/10.1080/02664763.2021.1967890⁠

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⁠

See Also

escalc for a function to compute various effect size measures based on means and standard deviations.

Examples

# example data frame
dat <- data.frame(case=c(1:3,NA), min=c(2,NA,2,NA), q1=c(NA,4,4,NA),
                  median=c(6,6,6,NA), q3=c(NA,10,10,NA), max=c(14,NA,14,NA),
                  mean=c(NA,NA,NA,7.0), sd=c(NA,NA,NA,4.2), n=c(20,20,20,20))
dat

# note that study 4 provides the mean and SD directly, while studies 1-3 provide five-number
# summary values or a subset thereof (corresponding to cases 1-3 above)

# estimate means/SDs (note: existing values in 'mean' and 'sd' are not touched)
dat <- conv.fivenum(min=min, q1=q1, median=median, q3=q3, max=max, n=n, data=dat)
dat

# check attributes (none of the tests are significant, so means/SDs are estimated for studies 1-3)
dfround(data.frame(attributes(dat$mean)), digits=3)

# calculate the log transformed coefficient of variation and corresponding sampling variance
dat <- escalc(measure="CVLN", mi=mean, sdi=sd, ni=n, data=dat)
dat

# fit equal-effects model to the estimates
res <- rma(yi, vi, data=dat, method="EE")
res

# estimated coefficient of variation (with 95% CI)
predict(res, transf=exp, digits=2)

############################################################################

# example data frame
dat <- data.frame(case=c(1:3,NA), min=c(2,NA,2,NA), q1=c(NA,4,4,NA),
                  median=c(6,6,6,NA), q3=c(NA,10,10,NA), max=c(14,NA,14,NA),
                  mean=c(NA,NA,NA,7.0), sd=c(NA,NA,NA,4.2), n=c(20,20,20,20))
dat

# try out different methods
conv.fivenum(min=min, q1=q1, median=median, q3=q3, max=max, n=n, data=dat)
set.seed(1234)
conv.fivenum(min=min, q1=q1, median=median, q3=q3, max=max, n=n, data=dat, method="qe")
conv.fivenum(min=min, q1=q1, median=median, q3=q3, max=max, n=n, data=dat, method="bc")
conv.fivenum(min=min, q1=q1, median=median, q3=q3, max=max, n=n, data=dat, method="mln")
conv.fivenum(min=min, q1=q1, median=median, q3=q3, max=max, n=n, data=dat, method="blue")

############################################################################

# example data frame
dat <- data.frame(case=c(1:3,NA), min=c(2,NA,2,NA), q1=c(NA,4,4,NA),
                  median=c(6,6,6,NA), q3=c(NA,10,14,NA), max=c(14,NA,20,NA),
                  mean=c(NA,NA,NA,7.0), sd=c(NA,NA,NA,4.2), n=c(20,20,20,20))
dat

# for study 3, the third quartile and maximum value suggest that the data have
# a right skewed distribution (they are much further away from the median than
# the minimum and first quartile)

# estimate means/SDs
dat <- conv.fivenum(min=min, q1=q1, median=median, q3=q3, max=max, n=n, data=dat)
dat

# note that the mean and SD are not estimated for study 3; this is because the
# test for skewness is significant for this study
dfround(data.frame(attributes(dat$mean)), digits=3)

# estimate means/SDs, but assume that the data for study 3 come from a log-normal distribution
# and back-transform the estimated mean/SD of the log-transformed data back to the raw data
dat <- conv.fivenum(min=min, q1=q1, median=median, q3=q3, max=max, n=n, data=dat,
                    dist=c("norm","norm","lnorm","norm"), replace="all")
dat

# this works now because the test for skewness of the log-transformed data is not significant
dfround(data.frame(attributes(dat$mean)), digits=3)

---

Convert Wald-Type Confidence Intervals and Tests to Sampling Variances

Description

Function to convert Wald-type confidence intervals (CIs) and test statistics (or the corresponding p-values) to sampling variances.

Usage

conv.wald(out, ci.lb, ci.ub, zval, pval, n, data, include,
          level=95, transf, check=TRUE, var.names, append=TRUE, replace="ifna", ...)

Arguments

out	vector with the observed effect sizes or outcomes.
ci.lb	vector with the lower bounds of the corresponding Wald-type CIs.
ci.ub	vector with the upper bounds of the corresponding Wald-type CIs.
zval	vector with the Wald-type test statistics.
pval	vector with the p-values of the Wald-type tests.
n	vector with the total sample sizes of the studies.
data	optional data frame containing the variables given to the arguments above.
include	optional (logical or numeric) vector to specify the subset of studies for which the conversion should be carried out.
level	numeric value (or vector) to specify the confidence interval level(s) (the default is 95; see here for details).
transf	optional argument to specify a function to transform out, ci.lb, and ci.ub (e.g., transf=log). If unspecified, no transformation is used.
check	logical to specify whether the function should carry out a check to examine if the point estimates fall (approximately) halfway between the CI bounds (the default is TRUE).
var.names	character vector with two elements to specify the name of the variable for the observed effect sizes or outcomes and the name of the variable for the corresponding sampling variances (if data is an object of class "escalc", the var.names are taken from the object; otherwise the defaults are "yi" and "vi").
append	logical to specify whether the data frame provided via the data argument should be returned together with the estimated values (the default is TRUE).
replace	character string or logical to specify how values in var.names should be replaced (only relevant when using the data argument and if variables in var.names already exist in the data frame). See the ‘Value’ section for more details.
...	other arguments.

Details

The escalc function can be used to compute a wide variety of effect sizes or ‘outcome measures’. However, the inputs required to compute certain measures with this function may not be reported for all of the studies. Under certain circumstances, other information (such as point estimates and corresponding confidence intervals and/or test statistics) may be available that can be converted into the appropriate format needed for a meta-analysis. The purpose of the present function is to facilitate this process.

The function typically takes a data frame created with the escalc function as input via the data argument. This object should contain variables yi and vi (unless argument var.names was used to adjust these variable names when the "escalc" object was created) for the observed effect sizes or outcomes and the corresponding sampling variances, respectively. For some studies, the values for these variables may be missing.

Converting Point Estimates and Confidence Intervals

In some studies, the effect size estimate or observed outcome may already be reported. If so, such values can be supplied via the out argument and are then substituted for missing yi values. At times, it may be necessary to transform the reported values (e.g., reported odds ratios to log odds ratios). Via argument transf, an appropriate transformation function can be specified (e.g., transf=log), in which case \(y_i = f(\textrm{out})\) where \(f(\cdot)\) is the function specified via transf.

Moreover, a confidence interval (CI) may have been reported together with the estimate. The bounds of the CI can be supplied via arguments ci.lb and ci.ub, which are also transformed if a function is specified via transf. Assume that the bounds were obtained from a Wald-type CI of the form \(y_i \pm z_{crit} \sqrt{v_i}\) (on the transformed scale if transf is specified), where \(v_i\) is the sampling variance corresponding to the effect size estimate or observed outcome (so that \(\sqrt{v_i}\) is the corresponding standard error) and \(z_{crit}\) is the appropriate critical value from a standard normal distribution (e.g., \(1.96\) for a 95% CI). Then \[v_i = \left(\frac{\textrm{ci.ub} - \textrm{ci.lb}}{2 \times z_{crit}}\right)^2\] is used to back-calculate the sampling variances of the (transformed) effect size estimates or observed outcomes and these values are then substituted for missing vi values in the dataset.

For example, consider the following dataset of three RCTs used as input for a meta-analysis of log odds ratios:

dat <- data.frame(study = 1:3,
                  cases.trt = c(23, NA, 4), n.trt = c(194, 183, 46),
                  cases.plc = c(38, NA, 7), n.plc = c(201, 188, 44),
                  oddsratio = c(NA, 0.64, NA), lower = c(NA, 0.33, NA), upper = c(NA, 1.22, NA))
dat <- escalc(measure="OR", ai=cases.trt, n1i=n.trt, ci=cases.plc, n2i=n.plc, data=dat)
dat

#   study cases.trt n.trt cases.plc n.plc oddsratio lower upper      yi     vi
# 1     1        23   194        38   201        NA    NA    NA -0.5500 0.0818
# 2     2        NA   183        NA   188      0.64  0.33  1.22      NA     NA
# 3     3         4    46         7    44        NA    NA    NA -0.6864 0.4437

where variable yi contains the log odds ratios and vi the corresponding sampling variances as computed from the counts and group sizes by escalc().

Study 2 does not report the counts (or sufficient information to reconstruct them), but the odds ratio and a corresponding 95% confidence interval (CI) directly, as given by variables oddsratio, lower, and upper. The CI is a standard Wald-type CI that was computed on the log scale (and whose bounds were then exponentiated). Then the present function can be used as follows:

dat <- conv.wald(out=oddsratio, ci.lb=lower, ci.ub=upper, data=dat, transf=log)
dat

#   study cases.trt n.trt cases.plc n.plc oddsratio lower upper      yi     vi
# 1     1        23   194        38   201        NA    NA    NA -0.5500 0.0818
# 2     2        NA   183        NA   188      0.64  0.33  1.22 -0.4463 0.1113
# 3     3         4    46         7    44        NA    NA    NA -0.6864 0.4437

Now variables yi and vi in the dataset are complete.

If the CI was not a 95% CI, then one can specify the appropriate level via the level argument. This can also be an entire vector in case different studies used different levels.

By default (i.e., when check=TRUE), the function carries out a rough check to examine if the point estimate falls (approximately) halfway between the CI bounds (on the transformed scale) for each study for which the conversion was carried out. A warning is issued if there are studies where this is not the case. This may indicate that a particular CI was not a Wald-type CI or was computed on a different scale (in which case the back-calculation above would be inappropriate), but can also arise due to rounding of the reported values (in which case the back-calculation would still be appropriate, albeit possibly a bit inaccurate). Care should be taken when using such back-calculated values in a meta-analysis.

Converting Test Statistics and P-Values

Similarly, study authors may report the test statistic and/or p-value from a Wald-type test of the form \(\textrm{zval} = y_i / \sqrt{v_i}\) (on the transformed scale if transf is specified), with the corresponding two-sided p-value given by \(\textrm{pval} = 2(1 - \Phi(\textrm{|zval|}))\), where \(\Phi(\cdot)\) denotes the cumulative distribution function of a standard normal distribution (i.e., pnorm). Test statistics and/or corresponding p-values of this form can be supplied via arguments zval and pval.

A given p-value can be back-transformed into the corresponding test statistic (if it is not already available) with \(\textrm{zval} = \Phi^{-1}(1 - \textrm{pval}/2)\), where \(\Phi^{-1}(\cdot)\) denotes the quantile function (i.e., the inverse of the cumulative distribution function) of a standard normal distribution (i.e., qnorm). Then \[v_i = \left(\frac{y_i}{\textrm{zval}}\right)^2\] is used to back-calculate a missing vi value in the dataset.

Note that the conversion of a p-value to the corresponding test statistic (which is then converted into sampling variance) as shown above assumes that the exact p-value is reported. If authors only report that the p-value fell below a certain threshold (e.g., \(p < .01\) or if authors only state that the test was significant – which typically implies \(p < .05\)), then a common approach is to use the value of the cutoff reported (e.g., if \(p < .01\) is reported, then assume \(p = .01\)), which is conservative (since the actual p-value was below that assumed value by some unknown amount). The conversion will therefore tend to be much less accurate.

Using the earlier example, suppose that only the odds ratio and the corresponding two-sided p-value from a Wald-type test (whether the log odds ratio differs significantly from zero) is reported for study 2.

dat <- data.frame(study = 1:3,
                  cases.trt = c(23, NA, 4), n.trt = c(194, 183, 46),
                  cases.plc = c(38, NA, 7), n.plc = c(201, 188, 44),
                  oddsratio = c(NA, 0.64, NA), pval = c(NA, 0.17, NA))
dat <- escalc(measure="OR", ai=cases.trt, n1i=n.trt, ci=cases.plc, n2i=n.plc, data=dat)
dat

  study cases.trt n.trt cases.plc n.plc oddsratio pval      yi     vi
1     1        23   194        38   201        NA   NA -0.5500 0.0818
2     2        NA   183        NA   188      0.64 0.17      NA     NA
3     3         4    46         7    44        NA   NA -0.6864 0.4437

Then the function can be used as follows:

dat <- conv.wald(out=oddsratio, pval=pval, data=dat, transf=log)
dat

#   study cases.trt n.trt cases.plc n.plc oddsratio pval      yi     vi
# 1     1        23   194        38   201        NA   NA -0.5500 0.0818
# 2     2        NA   183        NA   188      0.64 0.17 -0.4463 0.1058
# 3     3         4    46         7    44        NA   NA -0.6864 0.4437

Note that the back-calculated sampling variance for study 2 is not identical in these two examples, because the CI bounds and p-value are rounded to two decimal places, which introduces some inaccuracies. Also, if both (ci.lb, ci.ub) and either zval or pval is available for a study, then the back-calculation of \(v_i\) via the confidence interval is preferred.

Optionally, one can use the n argument to supply the total sample sizes of the studies. This has no relevance for the calculations done by the present function, but some other functions may use this information (e.g., when drawing a funnel plot with the funnel function and one adjusts the yaxis argument to one of the options that puts the sample sizes or some transformation thereof on the y-axis).

Value

If the data argument was not specified or append=FALSE, a data frame of class c("escalc","data.frame") with two variables called var.names[1] (by default "yi") and var.names[2] (by default "vi") with the (transformed) observed effect sizes or outcomes and the corresponding sampling variances (computed as described above).

If data was specified and append=TRUE, then the original data frame is returned. If var.names[1] is a variable in data and replace="ifna" (or replace=FALSE), then only missing values in this variable are replaced with the (possibly transformed) observed effect sizes or outcomes from out (where possible) and otherwise a new variable called var.names[1] is added to the data frame. Similarly, if var.names[2] is a variable in data and replace="ifna" (or replace=FALSE), then only missing values in this variable are replaced with the sampling variances back-calculated as described above (where possible) and otherwise a new variable called var.names[2] is added to the data frame.

If replace="all" (or replace=TRUE), then all values in var.names[1] and var.names[2] are replaced, even for cases where the value in var.names[1] and var.names[2] is not missing.

Note

A word of caution: Except for the check on the CI bounds, there is no possibility to determine if the back-calculations done by the function are appropriate in a given context. They are only appropriate when the CI bounds and tests statistics (or p-values) arose from Wald-type CIs / tests as described above. Using the same back-calculations for other purposes is likely to yield nonsensical values.

Author(s)

Wolfgang Viechtbauer [email protected] https://www.metafor-project.org

References

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⁠

See Also

escalc for a function to compute various effect size measures.

Examples

### a very simple example
dat <- data.frame(or=c(1.37,1.89), or.lb=c(1.03,1.60), or.ub=c(1.82,2.23))
dat

### convert the odds ratios and CIs into log odds ratios with corresponding sampling variances
dat <- conv.wald(out=or, ci.lb=or.lb, ci.ub=or.ub, data=dat, transf=log)
dat

############################################################################

### a more elaborate example based on the BCG vaccine dataset
dat <- dat.bcg[,c(2:7)]
dat

### with complete data, we can use escalc() in the usual way
dat1 <- escalc(measure="OR", ai=tpos, bi=tneg, ci=cpos, di=cneg, data=dat)
dat1

### random-effects model fitted to these data
res1 <- rma(yi, vi, data=dat1)
res1

### now suppose that the 2x2 table data are not reported in all studies, but that the
### following dataset could be assembled based on information reported in the studies
dat2 <- data.frame(summary(dat1))
dat2[c("yi", "ci.lb", "ci.ub")] <- data.frame(summary(dat1, transf=exp))[c("yi", "ci.lb", "ci.ub")]
names(dat2)[which(names(dat2) == "yi")] <- "or"
dat2[,c("or","ci.lb","ci.ub","pval")] <- round(dat2[,c("or","ci.lb","ci.ub","pval")], digits=2)
dat2$vi <- dat2$sei <- dat2$zi <- NULL
dat2$ntot <- with(dat2, tpos + tneg + cpos + cneg)
dat2[c(1,12),c(3:6,9:10)] <- NA
dat2[c(4,9), c(3:6,8)] <- NA
dat2[c(2:3,5:8,10:11,13),c(7:10)] <- NA
dat2$ntot[!is.na(dat2$tpos)] <- NA
dat2

### in studies 1 and 12, authors reported only the odds ratio and the corresponding p-value
### in studies 4 and 9, authors reported only the odds ratio and the corresponding 95% CI

### use escalc() first
dat2 <- escalc(measure="OR", ai=tpos, bi=tneg, ci=cpos, di=cneg, data=dat2)
dat2

### fill in the missing log odds ratios and sampling variances
dat2 <- conv.wald(out=or, ci.lb=ci.lb, ci.ub=ci.ub, pval=pval, n=ntot, data=dat2, transf=log)
dat2

### random-effects model fitted to these data
res2 <- rma(yi, vi, data=dat2)
res2

### any differences between res1 and res2 are a result of or, ci.lb, ci.ub, and pval being
### rounded in dat2 to two decimal places; without rounding, the results would be identical

---