rma.uni | R Documentation |
Function to fit meta-analytic equal-, fixed-, and random-effects models and (mixed-effects) meta-regression models using a linear (mixed-effects) model framework. See below and the introduction to the metafor-package for more details on these models. \loadmathjax
rma.uni(yi, vi, sei, weights, ai, bi, ci, di, n1i, n2i, x1i, x2i, t1i, t2i,
m1i, m2i, sd1i, sd2i, xi, mi, ri, ti, pi, sdi, r2i, ni, mods, scale,
measure="GEN", intercept=TRUE, data, slab, subset,
add=1/2, to="only0", drop00=FALSE, vtype="LS",
method="REML", weighted=TRUE, test="z",
level=95, btt, att, tau2, verbose=FALSE, digits, control, ...)
rma(yi, vi, sei, weights, ai, bi, ci, di, n1i, n2i, x1i, x2i, t1i, t2i,
m1i, m2i, sd1i, sd2i, xi, mi, ri, ti, pi, sdi, r2i, ni, mods, scale,
measure="GEN", intercept=TRUE, data, slab, subset,
add=1/2, to="only0", drop00=FALSE, vtype="LS",
method="REML", weighted=TRUE, test="z",
level=95, btt, att, tau2, verbose=FALSE, digits, control, ...)
yi |
vector of length \mjseqnk with the observed effect sizes or outcomes. See ‘Details’. |
vi |
vector of length \mjseqnk with the corresponding sampling variances. See ‘Details’. |
sei |
vector of length \mjseqnk with the corresponding standard errors (only relevant when not using |
weights |
optional argument to specify a vector of length \mjseqnk with user-defined weights. See ‘Details’. |
ai |
see below and the documentation of the |
bi |
see below and the documentation of the |
ci |
see below and the documentation of the |
di |
see below and the documentation of the |
n1i |
see below and the documentation of the |
n2i |
see below and the documentation of the |
x1i |
see below and the documentation of the |
x2i |
see below and the documentation of the |
t1i |
see below and the documentation of the |
t2i |
see below and the documentation of the |
m1i |
see below and the documentation of the |
m2i |
see below and the documentation of the |
sd1i |
see below and the documentation of the |
sd2i |
see below and the documentation of the |
xi |
see below and the documentation of the |
mi |
see below and the documentation of the |
ri |
see below and the documentation of the |
ti |
see below and the documentation of the |
pi |
see below and the documentation of the |
sdi |
see below and the documentation of the |
r2i |
see below and the documentation of the |
ni |
see below and the documentation of the |
mods |
optional argument to include one or more moderators in the model. A single moderator can be given as a vector of length \mjseqnk specifying the values of the moderator. Multiple moderators are specified by giving a matrix with \mjseqnk rows and as many columns as there are moderator variables. Alternatively, a model |
scale |
optional argument to include one or more predictors for the scale part in a location-scale model. See ‘Details’. |
measure |
character string to specify the type of data supplied to the function. When |
intercept |
logical to specify whether an intercept should be added to the model (the default is |
data |
optional data frame containing the data supplied to the function. |
slab |
optional vector with labels for the \mjseqnk studies. |
subset |
optional (logical or numeric) vector to specify the subset of studies that should be used for the analysis. |
add |
see the documentation of the |
to |
see the documentation of the |
drop00 |
see the documentation of the |
vtype |
see the documentation of the |
method |
character string to specify whether an equal- or a random-effects model should be fitted. An equal-effects model is fitted when using |
weighted |
logical to specify whether weighted (default) or unweighted estimation should be used to fit the model (the default is |
test |
character string to specify how test statistics and confidence intervals for the fixed effects should be computed. By default ( |
level |
numeric value between 0 and 100 to specify the confidence interval level (the default is 95). |
btt |
optional vector of indices to specify which coefficients to include in the omnibus test of moderators. Can also be a string to |
att |
optional vector of indices to specify which scale coefficients to include in the omnibus test. Only relevant for location-scale models. See ‘Details’. |
tau2 |
optional numeric 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 \mjseqn\tau^2). When unspecified, the value of \mjseqn\tau^2 is estimated from the data. |
verbose |
logical to specify whether output should be generated on the progress of the model fitting (the default is |
digits |
optional integer to specify the number of decimal places to which the printed results should be rounded. If unspecified, the default is 4. See also here for further details on how to control the number of digits in the output. |
control |
optional list of control values for the iterative estimation algorithms. If unspecified, default values are defined inside the function. See ‘Note’. |
... |
additional arguments. |
The function can be used in combination with any of the usual effect sizes or outcome measures used in meta-analyses (e.g., log risk ratios, log odds ratios, risk differences, mean differences, standardized mean differences, log transformed ratios of means, raw correlation coefficients, correlation coefficients transformed with Fisher's r-to-z transformation), or, more generally, any set of estimates (with corresponding sampling variances) one would like to analyze. Simply specify the observed effect sizes or 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 sizes or 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 sizes 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.
The function can be used to fit equal-, fixed-, and random-effects models, as well as (mixed-effects) meta-regression models including one or multiple moderators (the difference between the various models is described in detail on the introductory metafor-package help page).
Assuming the observed effect sizes or outcomes and corresponding sampling variances are supplied via the yi
and vi
arguments, an equal-effects model can be fitted with rma(yi, vi, method="EE")
. Setting method="FE"
fits a fixed-effects model (see here for a discussion of this model). 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
(which is the same as setting the weights equal to a constant).
A random-effects model can be fitted with the same code but setting the method
argument to one of the various estimators for the amount of heterogeneity:
method="DL"
= DerSimonian-Laird estimator (DerSimonian & Laird, 1986; Raudenbush, 2009),
method="HE"
= Hedges estimator (Hedges, 1983, 1992),
method="HS"
= Hunter-Schmidt estimator (Hunter & Schmidt, 1990; Viechtbauer et al., 2015),
method="HSk"
= Hunter-Schmidt estimator with a small sample-size correction (Brannick et al., 2019),
method="SJ"
= Sidik-Jonkman estimator (Sidik & Jonkman, 2005b, 2007),
method="ML"
= maximum likelihood estimator (Hardy & Thompson, 1996; Raudenbush, 2009),
method="REML"
= restricted maximum likelihood estimator (Viechtbauer, 2005; Raudenbush, 2009)
method="EB"
= empirical Bayes estimator (Morris, 1983; Berkey et al. 1995),
method="PM"
= Paule-Mandel estimator (Paule & Mandel, 1982; Viechtbauer et al., 2015),
method="GENQ"
= generalized Q-statistic estimator (DerSimonian & Kacker, 2007; Jackson et al., 2014),
method="PMM"
= median-unbiased Paule-Mandel estimator (Viechtbauer, 2021),
method="GENQM"
= median-unbiased generalized Q-statistic estimator (Viechtbauer, 2021).
For a description of the various estimators, see Brannick et al. (2019), DerSimonian and Kacker (2007), Raudenbush (2009), Veroniki et al. (2016), 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’, the empirical Bayes estimator is actually identical to the Paule-Mandel estimator (Viechtbauer et al., 2015), and 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 sampling variance weights, respectively). Finally, the two median-unbiased estimators are versions of the Paule-Mandel and generalized Q-statistic estimators that equate the respective estimating equations not to their expected values, but to the medians of their theoretical distributions (Viechtbauer, 2021).
One or more moderators can be included in a model via the mods
argument. A single moderator can be given as a (row or column) vector of length \mjseqnk specifying the values of the moderator. Multiple moderators are specified by giving an appropriate model matrix (i.e., \mjseqnX) with \mjseqnk 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 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 is controlled by the specified formula (e.g., mods = ~ mod1 + mod2 + mod3 - 1
would lead to the removal of the intercept).
When the observed effect sizes or outcomes and corresponding sampling variances are supplied via the yi
and vi
(or sei
) arguments, one can also 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 again is controlled by the specified formula.
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 (i.e., to test \mjseqn\mboxH_0:\; \beta_j \in \textttbtt = 0, where \mjseqn\beta_j \in \textttbtt is the set of coefficients to be tested). 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. The omnibus test is called the \mjseqnQ_M-test and follows asymptotically a chi-square distribution with \mjseqnm degrees of freedom (with \mjseqnm denoting the number of coefficients tested) under the null hypothesis (that the true value of all coefficients tested is equal to 0).
Categorical moderator variables can be included in the model via the mods
argument in the same way that appropriately (dummy) coded categorical 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 automate the coding (note that string/character variables in a model formula are automatically converted to factors). An example to illustrate these different approaches is provided below.
By default, tests of 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 (see above). As an alternative, one can set test="t"
, in which case tests of individual coefficients and confidence intervals are based on a t-distribution with \mjseqnk-p degrees of freedom, while the omnibus test then uses an F-distribution with \mjseqnm and \mjseqnk-p degrees of freedom (with \mjseqnk denoting the total number of estimates included in the analysis and \mjseqnp the total number of model coefficients including the intercept if it is present). Furthermore, when test="knha"
, the method by Hartung (1999), Sidik and Jonkman (2002), and Knapp and Hartung (2003) (the Knapp-Hartung method; also referred to as the Hartung-Knapp-Sidik-Jonkman method) is used, which applies an adjustment to the standard errors of the estimated coefficients (to account for the uncertainty in the estimate of the amount of (residual) heterogeneity) and uses t- and F-distributions as described above (see also here). Finally, one can set test="adhoc"
, in which case the Knapp-Hartung method is used, but with the restriction that the adjustment to the standard errors can never result in adjusted standard errors that are smaller than the unadjusted ones (see Jackson et al., 2017, section 4.3).
A test for (residual) heterogeneity is automatically carried out by the function. Without moderators in the model, this is simply Cochran's \mjseqnQ-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/outcomes are heterogeneous. When moderators are included in the model, this is the \mjseqnQ_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.
The function can also be used to fit so-called ‘location-scale models’ (Viechtbauer & López-López, 2022). In such models, one can specify not only predictors for the size of the average true outcome (i.e., for their ‘location’), but also predictors for the amount of heterogeneity in the outcomes (i.e., for their ‘scale’). The model is given by \mjdeqny_i = \beta_0 + \beta_1 x_i1 + \beta_2 x_i2 + ... + \beta_p' x_ip' + u_i + \varepsilon_i,y_i = \beta_0 + \beta_1 x_i1 + \beta_2 x_i2 + ... + \beta_p' x_ip' + u_i + \epsilon_i, \mjdeqnu_i \sim N(0, \tau_i^2), \; \varepsilon_i \sim N(0, v_i),u_i ~ N(0, tau_i^2), \epsilon_i \sim N(0, v_i), \mjdeqn\ln(\tau_i^2) = \alpha_0 + \alpha_1 z_i1 + \alpha_2 z_i2 + ... + \alpha_q' z_iq',log(tau^2) = \alpha_0 + \alpha z_i1 + \alpha z_i2 + ... + \alpha_q' z_iq', where \mjeqnx_i1, ..., x_ip'x_i1, ..., x_ip' are the values of the \mjseqnp' predictor variables that may be related to the size of the average true outcome (letting \mjseqnp = p' + 1 denote the total number of location coefficients in the model including the model intercept \mjseqn\beta_0) and \mjeqnz_i1, ..., z_iq'z_i1, ..., z_iq' are the values of the \mjseqnq' scale variables that may be related to the amount of heterogeneity in the outcomes (letting \mjseqnq = q' + 1 denote the total number of scale coefficients in the model including the model intercept \mjseqn\alpha_0). Location variables can be specified via the mods
argument as described above (e.g., mods = ~ mod1 + mod2 + mod3
). Scale variables can be specified via the scale
argument (e.g., scale = ~ var1 + var2 + var3
). A log link is used for specifying the relationship between the scale variables and the amount of heterogeneity so that \mjseqn\tau_i^2 is guaranteed to be non-negative (one can also set (the undocumented) argument link="identity"
to use an identity link, but this is more likely to lead to estimation problems). Estimates of the location and scale coefficients can be obtained either with maximum likelihood (method="ML"
) or restricted maximum likelihood (method="REML"
) estimation. An omnibus test of the scale coefficients is conducted as described above (where the att
argument can be used to specify which scale coefficients to include in the test).
An object of class c("rma.uni","rma")
. The object is a list containing the following components:
beta |
estimated coefficients of the model. |
se |
standard errors of the coefficients. |
zval |
test statistics of the coefficients. |
pval |
corresponding p-values. |
ci.lb |
lower bound of the confidence intervals for the coefficients. |
ci.ub |
upper bound of the confidence intervals for the coefficients. |
vb |
variance-covariance matrix of the estimated coefficients. |
tau2 |
estimated amount of (residual) heterogeneity. Always |
se.tau2 |
standard error of the estimated amount of (residual) heterogeneity. |
k |
number of studies included in the analysis. |
p |
number of coefficients in the model (including the intercept). |
m |
number of coefficients included in the omnibus test of moderators. |
QE |
test statistic of the test for (residual) heterogeneity. |
QEp |
corresponding p-value. |
QM |
test statistic of the omnibus test of moderators. |
QMp |
corresponding p-value. |
I2 |
value of \mjseqnI^2. See |
H2 |
value of \mjseqnH^2. See |
R2 |
value of \mjseqnR^2. See |
int.only |
logical that indicates whether the model is an intercept-only model. |
yi, vi, X |
the vector of outcomes, the corresponding sampling variances, and the model matrix. |
fit.stats |
a list with the log-likelihood, deviance, AIC, BIC, and AICc values under the unrestricted and restricted likelihood. |
... |
some additional elements/values. |
For location-scale models, the object is of class c("rma.ls","rma.uni","rma")
and includes the following components in addition to the ones listed above:
alpha |
estimated scale coefficients of the model. |
se.alpha |
standard errors of the coefficients. |
zval.alpha |
test statistics of the coefficients. |
pval.alpha |
corresponding p-values. |
ci.lb.alpha |
lower bound of the confidence intervals for the coefficients. |
ci.ub.alpha |
upper bound of the confidence intervals for the coefficients. |
va |
variance-covariance matrix of the estimated coefficients. |
tau2 |
as above, but now a vector of values. |
q |
number of scale coefficients in the model (including the intercept). |
QS |
test statistic of the omnibus test of the scale coefficients. |
QSp |
corresponding p-value. |
... |
some additional elements/values. |
The results of the fitted model are formatted and printed with the print
function. If fit statistics should also be given, use summary
(or use the fitstats
function to extract them). Full versus reduced model comparisons in terms of fit statistics and likelihood ratio tests can be obtained with anova
. 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
. Tests and confidence intervals based on (cluster) robust methods can be obtained with robust
.
Predicted/fitted values can be obtained with predict
and fitted
. For best linear unbiased predictions, see blup
and ranef
.
The residuals
, rstandard
, and rstudent
functions extract raw and standardized residuals. Additional model diagnostics (e.g., to determine influential studies) can be obtained with the influence
function. For models without moderators, leave-one-out diagnostics can also be obtained with leave1out
. For models with moderators, variance inflation factors can be obtained with vif
.
A confidence interval for the amount of (residual) heterogeneity in the random/mixed-effects model can be obtained with confint
. For location-scale models, confint
can provide confidence intervals for the scale coefficients.
Forest, funnel, radial, L'Abbé, and Baujat plots can be obtained with forest
, funnel
, radial
, labbe
, and baujat
(radial and L'Abbé plots only for models without moderators). The qqnorm
function provides normal QQ plots of the standardized residuals. One can also just call plot
on the fitted model object to obtain various plots at once. For random/mixed-effects models, the profile
function can be used to obtain a plot of the (restricted) log-likelihood as a function of \mjseqn\tau^2. For location-scale models, profile
draws analogous plots based on the scale coefficients. For models with moderators, regplot
draws scatter plots / bubble plots, showing the (marginal) relationship between the observed outcomes and a selected moderator from the model.
Tests for funnel plot asymmetry (which may be indicative of publication bias) can be obtained with ranktest
and regtest
. For models without moderators, the trimfill
method can be used to carry out a trim and fill analysis and hc
provides a random-effects model analysis that is more robust to publication bias (based on the method by Henmi & Copas, 2010). The test of ‘excess significance’ can be carried out with the tes
function. Selection models can be fitted with the selmodel
function.
For models without moderators, a cumulative meta-analysis (i.e., adding one observation at a time) can be obtained with cumul
.
Other extractor functions include coef
, vcov
, logLik
, deviance
, AIC
, BIC
, hatvalues
, and weights
.
While the HS, HSk, HE, DL, SJ, and GENQ estimators of \mjseqn\tau^2 are based on closed-form solutions, the ML, REML, and EB estimators must be obtained iteratively. For this, the function makes use of the Fisher scoring algorithm, which is robust to poor starting values and usually converges quickly (Harville, 1977; Jennrich & Sampson, 1976). By default, the starting value is set equal to the value of the Hedges (HE) estimator and the algorithm terminates when the change in the estimated value of \mjseqn\tau^2 is smaller than \mjeqn10^-510^(-5) from one iteration to the next. The maximum number of iterations is 100 by default (which should be sufficient in most cases). Information on the progress of the algorithm can be obtained by setting verbose=TRUE
. One can also set verbose
to an integer (verbose=2
yields even more information and verbose=3
also sets option(warn=1)
temporarily).
A different starting value, threshold, and maximum number of iterations can be specified via the control
argument by setting control=list(tau2.init=value, threshold=value, maxiter=value)
. The step length of the Fisher scoring algorithm can also be adjusted by a desired factor with control=list(stepadj=value)
(values below 1 will reduce the step length). If using verbose=TRUE
shows the estimate jumping around erratically (or cycling through a few values), decreasing the step length (and increasing the maximum number of iterations) can often help with convergence (e.g., control=list(stepadj=0.5, maxiter=1000)
).
The PM, PMM, and GENQM estimators also involve iterative algorithms, which make 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 100
(as above). The upper bound of the interval searched (tau2.max
) is set to the larger of 100 and 10*mad(yi)^2
(i.e., 10 times the squared median absolute deviation of the observed effect sizes or outcomes computed with the mad
function). These values can be adjusted with control=list(tol=value, maxiter=value, tau2.max=value)
.
All of the heterogeneity estimators except SJ can in principle yield negative estimates for the amount of (residual) heterogeneity. However, negative estimates of \mjseqn\tau^2 are outside of the parameter space. For the HS, HSk, HE, DL, and GENQ estimators, negative estimates are therefore truncated to zero. For the ML, REML, and EB estimators, the Fisher scoring algorithm makes use of step halving (Jennrich & Sampson, 1976) to guarantee a non-negative estimate. Finally, for the PM, PMM, and GENQM estimators, the lower bound of the interval searched is set to zero by default. For those brave enough to step into risky territory, there is the option to set the lower bound for all these estimators to some other value besides zero (even a negative one) with control=list(tau2.min=value)
, but the lowest value permitted is -min(vi)
(to ensure that the marginal variances are always non-negative).
The Hunter-Schmidt estimator for the amount of heterogeneity is defined in Hunter and Schmidt (1990) only in the context of the random-effects model when analyzing correlation coefficients. A general version of this estimator for random- and mixed-effects models not specific to any particular outcome measure is described in Viechtbauer (2005) and Viechtbauer et al. (2015) and is implemented here.
The Sidik-Jonkman estimator starts with a crude estimate of \mjseqn\tau^2, which is then updated as described in Sidik and Jonkman (2005b, 2007). If, instead of the crude estimate, one wants to use a better a priori estimate, one can do so by passing this value via control=list(tau2.init=value)
.
One can also specify a vector of estimators via the method
argument (e.g., rma(yi, vi, method=c("REML","DL"))
). The various estimators are then applied in turn until one converges. This is mostly useful for simulation studies where an estimator (like the REML estimator) is not guaranteed to converge and one can then substitute one (like the DL estimator) that does not involve iterative methods and is guaranteed to provide an estimate.
Outcomes with non-positive sampling variances are problematic. If a sampling variance is equal to zero, then its weight will be \mjseqn1/0 for equal-effects models when using weighted estimation. Switching to unweighted estimation is a possible solution then. For random/mixed-effects model, some estimators of \mjseqn\tau^2 are undefined when there is at least one sampling variance equal to zero. Other estimators may work, but it may still be necessary to switch to unweighted model fitting, especially when the estimate of \mjseqn\tau^2 converges to zero.
When including moderators in the model, it is possible that the model matrix is not of full rank (i.e., there is a linear relationship between the moderator variables included in the model). The function automatically tries to reduce the model matrix to full rank by removing redundant predictors, but if this fails the model cannot be fitted and an error will be issued. Deleting (redundant) moderator variables from the model as needed should solve this problem.
Some general words of caution about the assumptions underlying the models:
The sampling variances (i.e., the vi
values) are treated as if they are known constants, even though in practice they are usually estimates themselves. This implies that the distributions of the test statistics and corresponding confidence intervals are only exact and have nominal coverage when the within-study sample sizes are large (i.e., when the error in the sampling variance estimates is small). Certain outcome measures (e.g., the arcsine square root transformed risk difference and Fisher's r-to-z transformed correlation coefficient) are based on variance stabilizing transformations that also help to make the assumption of known sampling variances much more reasonable.
When fitting a mixed/random-effects model, \mjseqn\tau^2 is estimated and then treated as a known constant thereafter. This ignores the uncertainty in the estimate of \mjseqn\tau^2. As a consequence, the standard errors of the parameter estimates tend to be too small, yielding test statistics that are too large and confidence intervals that are not wide enough. The Knapp and Hartung (2003) adjustment (i.e., using test="knha"
) can be used to counter this problem, yielding test statistics and confidence intervals whose properties are closer to nominal.
Most effect sizes or outcome measures do not have exactly normal sampling distributions as assumed under the various models. However, the normal approximation usually becomes more accurate for most effect sizes or outcome measures as the within-study sample sizes increase. Therefore, sufficiently large within-study sample sizes are (usually) needed to be certain that the tests and confidence intervals have nominal levels/coverage. Again, certain outcome measures (e.g., Fisher's r-to-z transformed correlation coefficient) may be preferable from this perspective as well.
For location-scale models, model fitting is done via numerical optimization over the model parameters. By default, nlminb
is used for the optimization. One can also chose a different optimizer from optim
via the control
argument (e.g., control=list(optimizer="BFGS")
or control=list(optimizer="Nelder-Mead")
). Besides nlminb
and one of the methods from optim
, one can also choose one of the optimizers from the minqa
package (i.e., uobyqa
, newuoa
, or bobyqa
), one of the (derivative-free) algorithms from the nloptr
package, the Newton-type algorithm implemented in nlm
, the various algorithms implemented in the dfoptim
package (hjk
for the Hooke-Jeeves, nmk
for the Nelder-Mead, and mads
for the Mesh Adaptive Direct Searches algorithm), the quasi-Newton type optimizers ucminf
and lbfgsb3c
and the subspace-searching simplex algorithm subplex
from the packages of the same name, the Barzilai-Borwein gradient decent method implemented in BBoptim
, the Rcgmin
and Rvmmin
optimizers, or the parallelized version of the L-BFGS-B algorithm implemented in optimParallel
from the package of the same name. When using an identity link with link="identity"
, constrained optimization (to ensure non-negative \mjseqn\tau_i^2 values) as implemented in constrOptim
is used by default. Alternative optimizers in this case are the solnp
solver from the Rsolnp
package, nloptr
, or the augmented Lagrangian adaptive barrier minimization algorithm constrOptim.nl
from the alabama
package.
The optimizer name must be given as a character string (i.e., in quotes). Additional control parameters can be specified via the control
argument (e.g., control=list(iter.max=1000, rel.tol=1e-8)
). For nloptr
, the default is to use the BOBYQA implementation from that package with a relative convergence criterion of 1e-8
on the function value (i.e., log-likelihood), but this can be changed via the algorithm
and ftop_rel
arguments (e.g., control=list(optimizer="nloptr", algorithm="NLOPT_LN_SBPLX", ftol_rel=1e-6)
) (note: when using optimizer="nloptr"
in combination with an identity link, the "NLOPT_LN_COBYLA"
algorithm is automatically used, since this one allows for inequality constraints). For optimParallel
, the control argument ncpus
can be used to specify the number of cores to use for the parallelization (e.g., control=list(optimizer="optimParallel", ncpus=2)
). With parallel::detectCores()
, one can check on the number of available cores on the local machine.
Under certain circumstances (e.g., when the amount of heterogeneity is very small for certain combinations of values for the scale variables and scale coefficients), the values of the scale coefficients may try to drift towards minus or plus infinity, which can lead to problems with the optimization. One can impose constraints on the scale coefficients via control=list(alpha.min=minval, alpha.max=maxval)
where minval
and maxval
are either scalars or vectors of the appropriate length.
Finally, for location-scale models, the standard errors of the scale coefficients are obtained by inverting the Hessian, which is numerically approximated using the hessian
function from the numDeriv
package. This may fail (especially when using an identity link), leading to NA
values for the standard errors and hence test statistics, p-values, and confidence interval bounds. One can set control argument hessianCtrl
to a list of named arguments to be passed on to the method.args
argument of the hessian
function (the default is control=list(hessianCtrl=list(r=8))
). One can also set control=list(hesspack="pracma")
in which case the hessian
function from the pracma
package is used instead for approximating the Hessian.
Even if the Hessian can be approximated and inverted, the standard errors may be unreasonably large when the likelihood surface is very flat around the estimated scale coefficients. This is more likely to happen when \mjseqnk is small and when the amount of heterogeneity is very small under some conditions as defined by the scale coefficients/variables. Setting constraints on the scale coefficients as described above can also help to mitigate this issue.
Wolfgang Viechtbauer wvb@metafor-project.org https://www.metafor-project.org
Berkey, C. S., Hoaglin, D. C., Mosteller, F., & Colditz, G. A. (1995). A random-effects regression model for meta-analysis. Statistics in Medicine, 14(4), 395–411. https://doi.org/10.1002/sim.4780140406
Brannick, M. T., Potter, S. M., Benitez, B., & Morris, S. B. (2019). Bias and precision of alternate estimators in meta-analysis: Benefits of blending Schmidt–Hunter and Hedges approaches. Organizational Research Methods, 22(2), 490–514. https://doi.org/10.1177/1094428117741966
Cochran, W. G. (1954). The combination of estimates from different experiments. Biometrics, 10(1), 101–129. https://doi.org/10.2307/3001666
DerSimonian, R., & Laird, N. (1986). Meta-analysis in clinical trials. Controlled Clinical Trials, 7(3), 177–188. https://doi.org/10.1016/0197-2456(86)90046-2
DerSimonian, R., & Kacker, R. (2007). Random-effects model for meta-analysis of clinical trials: An update. Contemporary Clinical Trials, 28(2), 105–114. https://doi.org/10.1016/j.cct.2006.04.004
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<619::AID-SIM188>3.0.CO;2-A
Hartung, J. (1999). An alternative method for meta-analysis. Biometrical Journal, 41(8), 901–916. https://doi.org/10.1002/(SICI)1521-4036(199912)41:8<901::AID-BIMJ901>3.0.CO;2-W
Harville, D. A. (1977). Maximum likelihood approaches to variance component estimation and to related problems. Journal of the American Statistical Association, 72(358), 320–338. https://doi.org/10.2307/2286796
Hedges, L. V. (1983). A random effects model for effect sizes. Psychological Bulletin, 93(2), 388–395. https://doi.org/10.1037/0033-2909.93.2.388
Hedges, L. V. (1992). Meta-analysis. Journal of Educational Statistics, 17(4), 279–296. https://doi.org/10.3102/10769986017004279
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(29), 2969–2983. https://doi.org/10.1002/sim.4029
Hunter, J. E., & Schmidt, F. L. (1990). Methods of meta-analysis: Correcting error and bias in research findings. Thousand Oaks, CA: Sage.
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
Jackson, D., Law, M., Rücker, G., & Schwarzer, G. (2017). The Hartung-Knapp modification for random-effects meta-analysis: A useful refinement but are there any residual concerns? Statistics in Medicine, 36(25), 3923–3934. https://doi.org/10.1002/sim.7411
Jennrich, R. I., & Sampson, P. F. (1976). Newton-Raphson and related algorithms for maximum likelihood variance component estimation. Technometrics, 18(1), 11–17. https://doi.org/10.2307/1267911
Knapp, G., & Hartung, J. (2003). Improved tests for a random effects meta-regression with a single covariate. Statistics in Medicine, 22(17), 2693–2710. https://doi.org/10.1002/sim.1482
Morris, C. N. (1983). Parametric empirical Bayes inference: Theory and applications. Journal of the American Statistical Association, 78(381), 47–55. https://doi.org/10.2307/2287098
Paule, R. C., & Mandel, J. (1982). Consensus values and weighting factors. Journal of Research of the National Bureau of Standards, 87(5), 377–385. https://doi.org/10.6028/jres.087.022
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. (2002). A simple confidence interval for meta-analysis. Statistics in Medicine, 21(21), 3153–3159. https://doi.org/10.1002/sim.1262
Sidik, K., & Jonkman, J. N. (2005a). A note on variance estimation in random effects meta-regression. Journal of Biopharmaceutical Statistics, 15(5), 823–838. https://doi.org/10.1081/BIP-200067915
Sidik, K., & Jonkman, J. N. (2005b). Simple heterogeneity variance estimation for meta-analysis. Journal of the Royal Statistical Society, Series C, 54(2), 367–384. https://doi.org/10.1111/j.1467-9876.2005.00489.x
Sidik, K., & Jonkman, J. N. (2007). A comparison of heterogeneity variance estimators in combining results of studies. Statistics in Medicine, 26(9), 1964–1981. https://doi.org/10.1002/sim.2688
Veroniki, A. A., Jackson, D., Viechtbauer, W., Bender, R., Bowden, J., Knapp, G., Kuss, O., Higgins, J. P., Langan, D., & Salanti, G. (2016). Methods to estimate the between-study variance and its uncertainty in meta-analysis. Research Synthesis Methods, 7(1), 55–79. https://doi.org/10.1002/jrsm.1164
Viechtbauer, W. (2005). Bias and efficiency of meta-analytic variance estimators in the random-effects model. Journal of Educational and Behavioral Statistics, 30(3), 261–293. https://doi.org/10.3102/10769986030003261
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). Median-unbiased estimators for the amount of heterogeneity in meta-analysis. European Congress of Methodology, Valencia, Spain. https://www.wvbauer.com/lib/exe/fetch.php/talks:2021_viechtbauer_eam_median_tau2.pdf
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
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(3), 360–374. https://doi.org/10.1037/met0000023
rma.mh
, rma.peto
, rma.glmm
, and rma.mv
for other model fitting functions.
### calculate log risk ratios and corresponding sampling variances
dat <- escalc(measure="RR", ai=tpos, bi=tneg, ci=cpos, di=cneg, data=dat.bcg)
### fit a random-effects model using the log risk ratios and sampling variances as input
### note: method="REML" is the default, so one could leave this out
rma(yi, vi, data=dat, method="REML")
### fit a random-effects model using the log risk ratios and standard errors as input
### note: the second argument of rma() is for the *sampling variances*, so we use the
### named argument 'sei' to supply the standard errors to the function
dat$sei <- sqrt(dat$vi)
rma(yi, sei=sei, data=dat)
### fit a random-effects model supplying the 2x2 table cell frequencies to the function
rma(measure="RR", ai=tpos, bi=tneg, ci=cpos, di=cneg, data=dat)
### fit a mixed-effects model with two moderators (absolute latitude and publication year)
rma(yi, vi, mods=cbind(ablat, year), data=dat)
### using a model formula to specify the same model
rma(yi, vi, mods = ~ ablat + year, data=dat)
### using a model formula as part of the yi argument
rma(yi ~ ablat + year, vi, data=dat)
### 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,
### since this is the dummy/level we leave out of the model
### note: the intercept is the first coefficient, so with btt=2:3 we test
### coefficients 2 and 3, corresponding to the coefficients for the
### allocation factor
rma(yi, vi, mods = ~ alloc.random + alloc.systematic + year + ablat, data=dat, btt=2:3)
### using a model formula to specify the same model
rma(yi, vi, mods = ~ factor(alloc) + year + ablat, data=dat, btt=2:3)
### factor() is not needed as character variables are automatically converted to factors
rma(yi, vi, mods = ~ alloc + year + ablat, data=dat, btt=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)
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
############################################################################
### an example of a location-scale model
dat <- dat.bangertdrowns2004
### fit a standard random-effects model
res <- rma(yi, vi, data=dat)
res
### fit the same model as a location-scale model
res <- rma(yi, vi, scale = ~ 1, data=dat)
res
### check that we obtain the same estimate for tau^2
predict(res, newscale=1, transf=exp)
### add the total sample size (per 100) as a location and scale predictor
dat$ni100 <- dat$ni/100
res <- rma(yi, vi, mods = ~ ni100, scale = ~ ni100, data=dat)
res
### variables in the location and scale parts can differ
res <- rma(yi, vi, mods = ~ ni100 + meta, scale = ~ ni100 + imag, data=dat)
res
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.