Last chance! 50% off unlimited learning
Sale ends in
rma.uni(yi, vi, sei, weights, ai, bi, ci, di, n1i, n2i, x1i, x2i,
t1i, t2i, m1i, m2i, sd1i, sd2i, xi, mi, ri, ti, sdi, ni, mods,
measure="GEN", intercept=TRUE, data, slab, subset,
add=1/2, to="only0", drop00=FALSE, vtype="LS",
method="REML", weighted=TRUE, knha=FALSE,
level=95, digits=4, btt, tau2, verbose=FALSE, control)
rma(yi, vi, sei, weights, ai, bi, ci, di, n1i, n2i, x1i, x2i,
t1i, t2i, m1i, m2i, sd1i, sd2i, xi, mi, ri, ti, sdi, ni, mods,
measure="GEN", intercept=TRUE, data, slab, subset,
add=1/2, to="only0", drop00=FALSE, vtype="LS",
method="REML", weighted=TRUE, knha=FALSE,
level=95, digits=4, btt, tau2, verbose=FALSE, control)escalc function for more details.escalc function for more details.escalc function for more details.escalc function for more details.escalc function for more details.escalc function for more details.escalc function for more details.escalc function for more details.escalc function for more details.escalc function for more details.escalc function for more details.escalc function for more details.escalc function for more details.escalc function for more details.escalc function for more details.escalc function for more details.escalc function for more details.escalc function for more details.escalc function for more details.escalc function for more details.measure="GEN" (default), the observed effect sizes or outcomes and corresponding sampling variances (or standard errors) should be supplied to the function via the TRUE).escalc function.escalc function.escalc function.escalc function.method="FE". Random/mixed-effects models are fitted by setting methodFALSE). See FALSE). See c("rma.uni","rma"). The object is a list containing the following components:0 when method="FE".print.rma.uni for more details.print.rma.uni for more details.print.rma.uni for more details.print.rma.uni function. If fit statistics should also be given, use summary.rma (or use the fitstats.rma function to extract them). Full versus reduced model comparisons in terms of fit statistics and likelihoods can be obtained with anova.rma.uni. Wald-type tests for sets of model coefficients can be obtained with the same function. Permutation tests for the model coefficient(s) can be obtained with permutest.rma.uni.
Predicted/fitted values can be obtained with predict.rma and fitted.rma. For best linear unbiased predictions, see blup.rma.uni.
The residuals.rma, rstandard.rma.uni, and rstudent.rma.uni functions extract raw and standardized residuals. Additional case diagnostics (e.g., to determine influential studies) can be obtained with the influence.rma.uni function. For models without moderators, leave-one-out diagnostics can also be obtained with leave1out.rma.uni.
A confidence interval for the amount of (residual) heterogeneity in the random/mixed-effects model can be obtained with confint.rma.uni.
Forest, funnel, radial, and L'abbe plots (the latter two only for models without moderators) can be obtained with forest.rma, funnel.rma, radial.rma, and labbe.rma. The qqnorm.rma.uni function provides normal QQ plots of the standardized residuals. One can also just call plot.rma.uni on the fitted model object to obtain various plots at once.
Tests for funnel plot asymmetry (which may be indicative of publication bias) can be obtained with ranktest.rma and regtest.rma. For models without moderators, the trimfill.rma.uni method can be used to carry out a trim and fill analysis and code{hc.rma.uni} provides the results from the method by Henmi and Copas (2010).
For models without moderators, a cumulative meta-analysis (i.e., adding one obervation at a time) can be obtained with cumul.rma.uni.
Other assessor functions include coef.rma, vcov.rma, logLik.rma, deviance.rma, AIC.rma, BIC.rma, hatvalues.rma.uni, and weights.rma.uni.yi argument and the corresponding sampling variances via the vi argument. Instead of specifying vi, one can specify the standard errors (the square root of the sampling variances) via the sei argument or the inverse of the sampling variances via the weights argument.
Alternatively, the function can automatically calculate the values of a chosen effect size or outcome measure (and the corresponding sampling variances) when supplied with the necessary data. The escalc function describes which effect size or outcome measures are currently implemented and what data/arguments should then be specified/used. The measure argument should then be set to the desired effect size or outcome measure.
Specifying the Model
The function can be used to fit fixed- and random/mixed-effects models, as well as meta-regression models including moderators (the difference between the various models is described in detail in the introductory yi and vi, the fixed-effects model is fitted with rma(yi, vi, method="FE"). The random-effects model is fitted with the same code but setting method to one of the various estimators for the amount of heterogeneity:
method="DL"= DerSimonian-Laird estimatormethod="HE"= Hedges estimatormethod="HS"= Hunter-Schmidt estimatormethod="SJ"= Sidik-Jonkman estimatormethod="ML"= maximum-likelihood estimatormethod="REML"= restricted maximum-likelihood estimatormethod="EB"= empirical Bayes estimator.mods argument. A single moderator can be given as a (row or column) vector of length $k$ specifying the values of the moderator. Multiple moderators are specified by giving an appropriate design matrix with $k$ rows and $p'$ columns (e.g., using mods = cbind(mod1, mod2, mod3), where mod1, mod2, and mod3 correspond to the names of the variables for the three moderator variables). The intercept is included in the model by default unless intercept=FALSE.
Alternatively, one can use standard formula syntax to specify the model. In this case, the mods argument should be set equal to a one-sided formula of the form mods = ~ model (e.g., mods = ~ mod1 + mod2 + mod3). Interactions, polynomial terms, and factors can be easily added to the model in this manner. When specifying a model formula via the mods argument, the intercept argument is ignored. Instead, the inclusion/exclusion of the intercept term is controlled by the specified formula (e.g., mods = ~ mod1 + mod2 + mod3 - 1 would lead to the removal of the intercept term).
A fixed-effects with moderators model is then fitted by setting method="FE", while a mixed-effects model is fitted by specifying one of the estimators for the amount of (residual) heterogeneity given earlier.
When the observed outcomes and corresponding sampling variances are supplied via the yi and vi arguments (or sei or weights), one can also directly specify moderators via the yi argument (e.g., rma(yi ~ mod1 + mod2 + mod3, vi)). In that case, the mods argument is ignored and the inclusion/exclusion of the intercept term again is controlled by the specified formula.
Omnibus Test of Parameters
For models including moderators, an omnibus test of all the 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 of the coefficients in the model including the first. Alternatively, one can manually specify the indices of the coefficients to test via the btt argument. For example, use btt=c(3,4) to only include the third and fourth coefficient from the model in the test (if an intercept is included in the model, then it corresponds to the first coefficient in the model).
Categorical Moderators
Categorical moderator variables can be included in the model in the same way that appropriately (dummy) coded categorical independent variables can be included in linear models. One can either do the dummy coding manually or use a model formula together with the factor function to let Rhandle the coding automatically. An example to illustrate these different approaches is provided below.
Knapp & Hartung Adjustment
By default, the test statistics of the individual coefficients in the model (and the corresponding confidence intervals) are based on the normal distribution, while the omnibus test is based on a chi-square distribution with $m$ degrees of freedom ($m$ being the number of coefficients tested). The Knapp and Hartung (2003) method (knha=TRUE) is an adjustment to the standard errors of the estimated coefficients, which helps to account for the uncertainty in the estimate of the amount of (residual) heterogeneity and leads to different reference distributions. Tests of individual coefficients and confidence intervals are then based on the t-distribution with $k-p$ degrees of freedom, while the omnibus test statistic then uses an F-distribution with $m$ and $k-p$ degrees of freedom ($p$ being the total number of model coefficients including the intercept if it is present). The Knapp and Hartung (2003) adjustment is only meant to be used in the context of random- or mixed-effects models.rma.mh, rma.peto, rma.glmm### load BCG vaccine data
data(dat.bcg)
### calculate log relative risks and corresponding sampling variances
dat <- escalc(measure="RR", ai=tpos, bi=tneg, ci=cpos, di=cneg, data=dat.bcg)
### random-effects model (method="REML" is default, so technically not needed)
rma(yi, vi, data=dat, method="REML")
rma(yi, sei=sqrt(vi), data=dat, method="REML")
rma(yi, weights=1/vi, data=dat, method="REML")
### supplying the cell frequencies directly to the function
rma(measure="RR", ai=tpos, bi=tneg, ci=cpos, di=cneg, data=dat, method="REML")
### mixed-effects model with two moderators (absolute latitude and publication year)
rma(yi, vi, mods=cbind(ablat, year), data=dat, method="REML")
### using a model formula to specify the same model
rma(yi, vi, mods = ~ ablat + year, data=dat, method="REML")
### using a model formula in the yi argument
rma(yi ~ ablat + year, vi, data=dat, method="REML")
### manual dummy coding of the allocation factor
alloc.random <- ifelse(dat$alloc == "random", 1, 0)
alloc.alternate <- ifelse(dat$alloc == "alternate", 1, 0)
alloc.systematic <- ifelse(dat$alloc == "systematic", 1, 0)
### test the allocation factor (in the presence of the other moderators)
### note: "alternate" is the reference level of the allocation factor
### note: the intercept is the first coefficient, so btt=c(2,3)
rma(yi, vi, mods = ~ alloc.random + alloc.systematic + year + ablat,
data=dat, method="REML", btt=c(2,3))
### using a model formula to specify the same model
rma(yi, vi, mods = ~ factor(alloc) + year + ablat, data=dat, method="REML", btt=c(2,3))
### all pairwise differences using the 'multcomp' package (with Holm's method)
res <- rma(yi, vi, mods = ~ factor(alloc) - 1, data=dat, method="REML")
res
summary(glht(res, linfct=contrMat(rep(1,3), type=("Tukey"))), test=adjusted("holm"))
### subgrouping versus using a single model with a factor (subgrouping provides
### an estimate of tau^2 within each subgroup, but the number of studies in each
### subgroup get quite small; the model with the allocation factor provides a
### single estimate of tau^2 based on a larger number of studies, but assumes
### that tau^2 is the same within each subgroup)
res.a <- rma(yi, vi, data=dat, subset=(alloc=="alternate"))
res.r <- rma(yi, vi, data=dat, subset=(alloc=="random"))
res.s <- rma(yi, vi, data=dat, subset=(alloc=="systematic"))
res.a
res.r
res.s
res <- rma(yi, vi, mods = ~ factor(alloc) - 1, data=dat)
res
### demonstrating that Q_E + Q_M = Q_Total for fixed-effects models
### note: this does not work for random/mixed-effects models, since Q_E and
### Q_Total are calculated under the assumption that tau^2 = 0, while the
### calculation of Q_M incorporates the estimate of tau^2
res <- rma(yi, vi, data=dat, method="FE")
res ### this gives Q_Total
res <- rma(yi, vi, mods = ~ ablat + year, data=dat, method="FE")
res ### this gives Q_E and Q_M
res$QE + res$QM
### decomposition of Q_E into subgroup Q-values
res <- rma(yi, vi, mods = ~ factor(alloc), data=dat)
res
res.a <- rma(yi, vi, data=dat, subset=(alloc=="alternate"))
res.r <- rma(yi, vi, data=dat, subset=(alloc=="random"))
res.s <- rma(yi, vi, data=dat, subset=(alloc=="systematic"))
res.a$QE ### Q-value within subgroup
res.r$QE ### Q-value within subgroup
res.s$QE ### Q-value within subgroup
res$QE
res.a$QE + res.r$QE + res.s$QERun the code above in your browser using DataLab