metafor (version 1.9-9)

rma.uni: Meta-Analysis via Linear (Mixed-Effects) Models

Description

Function to fit the meta-analytic fixed- and random/mixed-effects models with or without moderators via linear (mixed-effects) models. See the documentation of the metafor-package for more details on these models.

Usage

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, test="z", 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, test="z", level=95, digits=4, btt, tau2, verbose=FALSE, control, ...)

Arguments

yi
vector of length $k$ with the observed effect sizes or outcomes. See ‘Details’.
vi
vector of length $k$ with the corresponding sampling variances. See ‘Details’.
sei
vector of length $k$ with the corresponding standard errors (only relevant when not using vi). See ‘Details’.
weights
optional argument to specify a vector of length $k$ with user-defined weights. See ‘Details’.
ai
see below and the documentation of the escalc function for more details.
bi
see below and the documentation of the escalc function for more details.
ci
see below and the documentation of the escalc function for more details.
di
see below and the documentation of the escalc function for more details.
n1i
see below and the documentation of the escalc function for more details.
n2i
see below and the documentation of the escalc function for more details.
x1i
see below and the documentation of the escalc function for more details.
x2i
see below and the documentation of the escalc function for more details.
t1i
see below and the documentation of the escalc function for more details.
t2i
see below and the documentation of the escalc function for more details.
m1i
see below and the documentation of the escalc function for more details.
m2i
see below and the documentation of the escalc function for more details.
sd1i
see below and the documentation of the escalc function for more details.
sd2i
see below and the documentation of the escalc function for more details.
xi
see below and the documentation of the escalc function for more details.
mi
see below and the documentation of the escalc function for more details.
ri
see below and the documentation of the escalc function for more details.
ti
see below and the documentation of the escalc function for more details.
sdi
see below and the documentation of the escalc function for more details.
ni
see below and the documentation of the escalc function for more details.
mods
optional argument to include one or more moderators in the model. A single moderator can be given as a vector of length $k$ specifying the values of the moderator. Multiple moderators are specified by giving a matrix with $k$ rows and as many columns as there are moderator variables. Alternatively, a model formula can be used to specify the model. See ‘Details’.
measure
character string indicating the type of data supplied to the function. When measure="GEN" (default), the observed effect sizes or outcomes and corresponding sampling variances (or standard errors) should be supplied to the function via the yi, vi, and sei arguments (only one of the two, vi or sei, needs to be specified). Alternatively, one can set measure to one of the effect size or outcome measures described under the documentation for the escalc function and specify the needed data via the appropriate arguments.
intercept
logical indicating whether an intercept term should be added to the model (the default is TRUE).
data
optional data frame containing the data supplied to the function.
slab
optional vector with labels for the $k$ studies.
subset
optional vector indicating the subset of studies that should be used for the analysis. This can be a logical vector of length $k$ or a numeric vector indicating the indices of the observations to include.
add
see the documentation of the escalc function.
to
see the documentation of the escalc function.
drop00
see the documentation of the escalc function.
vtype
see the documentation of the escalc function.
method
character string specifying whether a fixed- or a random/mixed-effects model should be fitted. A fixed-effects model (with or without moderators) is fitted when using method="FE". Random/mixed-effects models are fitted by setting method equal to one of the following: "DL", "HE", "SJ", "ML", "REML", "EB", "HS", or "GENQ". Default is "REML". See ‘Details’.
weighted
logical indicating whether weighted (default) or unweighted estimation should be used to fit the model.
test
character string specifying how test statistics and confidence intervals for the fixed effects should be computed. By default (test="z"), Wald-type tests and CIs are obtained, which are based on a standard normal distribution. When test="knha", the method by Knapp and Hartung (2003) is used for adjusting test statistics and confidence intervals. See ‘Details’.
level
numerical value between 0 and 100 specifying the confidence interval level (the default is 95).
digits
integer specifying the number of decimal places to which the printed results should be rounded (the default is 4).
btt
optional vector of indices specifying which coefficients to include in the omnibus test of moderators. See ‘Details’.
tau2
optional numerical value to specify the amount of (residual) heterogeneity in a random- or mixed-effects model (instead of estimating it). Useful for sensitivity analyses (e.g., for plotting results as a function of \tau²). When unspecified, the value of \tau² is estimated from the data.
verbose
logical indicating whether output should be generated on the progress of the model fitting (the default is FALSE). Can also be an integer. Values > 1 generate more verbose output. See ‘Note’.
control
optional list of control values for the iterative estimation algorithms. If unspecified, default values are defined inside the function. See ‘Note’.
...
additional arguments.

Value

An object of class c("rma.uni","rma"). The object is a list containing the following components:The results of the fitted model are formated and printed with the 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. Wald-type tests for sets of model coefficients or linear combinations thereof can be obtained with the same function. Permutation tests for the model coefficient(s) can be obtained with permutest.rma.uni. Tests and confidence intervals based on (cluster) robust methods can be obtained with robust.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, L'abbé, and Baujat plots can be obtained with forest.rma, funnel.rma, radial.rma, labbe.rma, and baujat (radial and L'abbé plots only for models without moderators). 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. For random/mixed-effects models, the profile.rma.uni function can be used to obtain a plot of the (restricted) log-likelihood as a function of \tau².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 hc.rma.uni provides a random-effects model analysis that is more robust to publication bias (based on the method by Henmi & 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 extractor functions include coef.rma, vcov.rma, logLik.rma, deviance.rma, AIC.rma, BIC.rma, hatvalues.rma.uni, and weights.rma.uni.

Details

Specifying the Data

The function can be used in conjunction with any of the usual effect size or outcome measures used in meta-analyses (e.g., log odds ratios, log relative risks, risk differences, mean differences, standardized mean differences, raw correlation coefficients, correlation coefficients transformed with Fisher's r-to-z transformation, and so on). Simply specify the observed outcomes via the 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. The escalc function can be used to compute a wide variety of effect size and outcome measures (and the corresponding sampling variances) based on summary statistics.

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 metafor-package help file).

Assuming the observed outcomes and corresponding sampling variances are supplied via yi and vi, a fixed-effects model can be fitted with rma(yi, vi, method="FE"). Weighted estimation (with inverse-variance weights) is used by default. User-defined weights can be supplied via the weights argument. Unweighted estimation can be used by setting weighted=FALSE (same as setting the weights equal to a constant).

A random-effects model can be fitted with the same code but setting method to one of the various estimators for the amount of heterogeneity:

  • method="DL" = DerSimonian-Laird estimator
  • method="HE" = Hedges estimator
  • method="HS" = Hunter-Schmidt estimator
  • method="SJ" = Sidik-Jonkman estimator
  • method="ML" = maximum-likelihood estimator
  • method="REML" = restricted maximum-likelihood estimator
  • method="EB" = empirical Bayes estimator
  • method="PM" = Paule-Mandel estimator
  • method="GENQ" = generalized Q-statistic estimator
For a description of the various estimators, see DerSimonian and Kacker (2007), Raudenbush (2009), Viechtbauer (2005), and Viechtbauer et al. (2015). Note that the Hedges estimator is also called the ‘variance component estimator’ or ‘Cochran estimator’, the Sidik-Jonkman estimator is also called the ‘model error variance estimator’, and the empirical Bayes estimator is actually identical to the Paule-Mandel estimator (Paule & Mandel, 1982). Finally, the generalized Q-statistic estimator is a general method-of-moments estimator (DerSimonian & Kacker, 2007) requiring the specification of weights (the HE and DL estimators are just special cases with equal and inverse variance weights, respectively).

One or more moderators can be included in these models via the 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 model matrix (i.e., $X$) with $k$ rows and as many columns as there are moderator variables (e.g., mods = cbind(mod1, mod2, mod3), where mod1, mod2, and mod3 correspond to the names of the variables for three moderator variables). The intercept is added to the model matrix by default unless intercept=FALSE.

Alternatively, one can use the 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 (or sei) arguments, 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, with btt=c(3,4), only the third and fourth coefficient from the model would be included 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 via the mods argument 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 R handle 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 a standard 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 (test="knha") 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.

Test for (Residual) Heterogeneity

A test for (residual) heterogeneity is automatically carried out by the function. Without moderators in the model, this is simply Cochran's $Q$-test (Cochran, 1954), which tests whether the variability in the observed effect sizes or outcomes is larger than would be expected based on sampling variability alone. A significant test suggests that the true effects or outcomes are heterogeneous. When moderators are included in the model, this is the $Q_E$-test for residual heterogeneity, which tests whether the variability in the observed effect sizes or outcomes not accounted for by the moderators included in the model is larger than would be expected based on sampling variability alone.

References

Berkey, C. S., Hoaglin, D. C., Mosteller, F., & Colditz, G. A. (1995). A random-effects regression model for meta-analysis. Statistics in Medicine, 14, 395--411.

Cochran, W. G. (1954). The combination of estimates from different experiments. Biometrics, 10, 101--129.

DerSimonian, R., & Laird, N. (1986). Meta-analysis in clinical trials. Controlled Clinical Trials, 7, 177--188.

DerSimonian, R., & Kacker, R. (2007). Random-effects model for meta-analysis of clinical trials: An update. Contemporary Clinical Trials, 28, 105--114.

Harville, D. A. (1977). Maximum likelihood approaches to variance component estimation and to related problems. Journal of the American Statistical Association, 72, 320--338.

Hedges, L. V. (1983). A random effects model for effect sizes. Psychological Bulletin, 93, 388--395.

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

Henmi, M., & Copas, J. B. (2010). Confidence intervals for random effects meta-analysis and robustness to publication bias. Statistics in Medicine, 29, 2969--2983.

Hunter, J. E., & Schmidt, F. L. (2004). Methods of meta-analysis: Correcting error and bias in research findings (2nd ed.). Thousand Oaks, CA: Sage.

Jennrich, R. I., & Sampson, P. F. (1976). Newton-Raphson and related algorithms for maximum likelihood variance component estimation. Technometrics, 18, 11--17.

Knapp, G., & Hartung, J. (2003). Improved tests for a random effects meta-regression with a single covariate. Statistics in Medicine, 22, 2693--2710.

Morris, C. N. (1983). Parametric empirical Bayes inference: Theory and applications (with discussion). Journal of the American Statistical Association, 78, 47--65.

Paule, R. C., & Mandel, J. (1982). Consensus values and weighting factors. Journal of Research of the National Bureau of Standards, 87, 377--385.

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.

Sidik, K., & Jonkman, J. N. (2005a). A note on variance estimation in random effects meta-regression. Journal of Biopharmaceutical Statistics, 15, 823--838.

Sidik, K., & Jonkman, J. N. (2005b). Simple heterogeneity variance estimation for meta-analysis. Journal of the Royal Statistical Society, Series C, 54, 367--384.

Sidik, K., & Jonkman, J. N. (2007). A comparison of heterogeneity variance estimators in combining results of studies. Statistics in Medicine, 26, 1964--1981.

Viechtbauer, W. (2005). Bias and efficiency of meta-analytic variance estimators in the random-effects model. Journal of Educational and Behavioral Statistics, 30, 261--293.

Viechtbauer, W. (2010). Conducting meta-analyses in R with the metafor package. Journal of Statistical Software, 36(3), 1--48. http://www.jstatsoft.org/v36/i03/.

Viechtbauer, W., López-López, J. A., Sánchez-Meca, J., & Marín-Martínez, F. (2015). A comparison of procedures to test for moderators in mixed-effects meta-regression models. Psychological Methods, 20, 360--374.

See Also

rma.mh, rma.peto, rma.glmm, and rma.mv for other model fitting functions.

Examples

Run this code
### 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")

### supplying the 2x2 table cell frequencies directly to the rma() 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 as part of 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))

### test all pairwise differences with Holm's method (using the 'multcomp' package if installed)
res <- rma(yi, vi, mods = ~ factor(alloc) - 1, data=dat, method="REML")
res
if (require(multcomp))
   summary(glht(res, linfct=contrMat(c("alternate"=1,"random"=1,"systematic"=1),
           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 is 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 "alternate"
res.r$QE ### Q-value within subgroup "random"
res.s$QE ### Q-value within subgroup "systematic"

res$QE
res.a$QE + res.r$QE + res.s$QE

Run the code above in your browser using DataCamp Workspace