rma.uni(yi, vi, sei, ai, bi, ci, di, n1i, n2i, m1i, m2i, sd1i, sd2i,
xi, mi, ri, ni, mods=NULL, data=NULL, intercept=TRUE, slab=NULL,
subset=NULL, measure="GEN", add=1/2, to="only0", vtype="LS",
method="REML", weighted=TRUE, level=95, digits=4, btt=NULL,
tau2=NULL, knha=FALSE, control=list())
rma(yi, vi, sei, ai, bi, ci, di, n1i, n2i, m1i, m2i, sd1i, sd2i,
xi, mi, ri, ni, mods=NULL, data=NULL, intercept=TRUE, slab=NULL,
subset=NULL, measure="GEN", add=1/2, to="only0", vtype="LS",
method="REML", weighted=TRUE, level=95, digits=4, btt=NULL,
tau2=NULL, knha=FALSE, control=list())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.TRUE).measure="GEN" (default), the observed effect sizes or outcomes and corresponding sampling variances (or standard errors) should be supplied to the function via the escalc function.escalc function.escalc function.method="FE". Random/mixed-effects models are fitted by setting methFALSE). See c("rma.uni", "rma"). The object is a list containing the following components:0 when method="FE".NA otherwise).length(yi) unless subset was used or if there are missing data).NA otherwise).print.rma.uni function. If you also want the fit statistics, use summary.rma (or use the fitstats.rma function to extract them).
Predicted/fitted values can be obtained with predict.rma.uni 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.
A confidence interval for the amount of (residual) heterogeneity in the random/mixed-effects model can be obtained with cint.rma.uni.
Forest, funnel, and radial plots (the latter only for models without moderators) are drawn with forest.rma, funnel.rma, and radial.rma. The qqnorm.rma.uni function provides a normal QQ plot of the standardized residuals. One can also just call plot.rma.uni on the fitted model object to obtain various plots at once.
Other assessor functions include coef.rma, vcov.rma, logLik.rma, and hatvalues.rma.uni.yi argument and the corresponding sampling variances via the vi argument (or supply the standard errors, the square root of sampling variances, via the sei argument).
Alternatively, the function can automatically calculate many effect size or outcome measures when supplied with the needed data. The escalc function describes which measures are currently implemented and what data/arguments should then be specified. The measure argument should then be set to the desired measure.
Specifying the Model
Assuming the observed outcomes and corresponding sampling variances are supplied via 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:
"HS"= Hunter-Schmidt estimator"HE"= Hedges estimator"DL"= DerSimonian-Laird estimator"SJ"= Sidik-Jonkman estimator"ML"= maximum-likelihood estimator"REML"= restricted maximum-likelihood estimator"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 a matrix with $k$ rows and $p'$ columns (e.g., using mods = cbind(mod1, mod2, mod3), where mod1, mod2, mod3 correspond to the names of the variables for the three moderator variables). 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. The intercept is automatically included in the model, unless intercept=FALSE is used.
Omnibus Test of Parameters
In models with more than one independent variable, an omnibus test of all the regression coefficients is conducted that excludes the intercept (the first coefficient) if the option intercept=TRUE is used (which is the default). If intercept=FALSE, then the omnibus test includes all of the coefficients in the model. Alternatively, one can 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.
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. You have to do the dummy coding yourself or use the model.matrix function to do the coding for you. An example is shown 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. Individual coefficients and confidence intervals are then based on the t-distribution with $k-p$ degrees of freedom ($p$ being the total number of coefficients in the model), while the omnibus test statistic then uses an F-distribution with $m$ and $k-p$ degrees of freedom. The Knapp and Hartung (2003) method is only meant to be used in the context of random- or mixed-effects models.### load BCG vaccine data
data(dat.bcg)
### calculate log risk rates and corresponding sampling variances
dat <- escalc(measure="RR", ai=tpos, bi=tneg, ci=cpos, di=cneg, data=dat.bcg)
dat.bcg$yi <- dat$yi
dat.bcg$vi <- dat$vi
### mixed-effects model with two moderators (absolute latitude and publication year)
rma(yi, vi, mods=cbind(ablat, year), data=dat.bcg, method="REML")
### supplying the raw data directly to the function
rma(ai=tpos, bi=tneg, ci=cpos, di=cneg, mods=cbind(ablat, year),
data=dat.bcg, measure="RR", method="REML")
### dummy coding of the allocation factor
alloc.random <- ifelse(dat.bcg$alloc == "random", 1, 0)
alloc.alternate <- ifelse(dat.bcg$alloc == "alternate", 1, 0)
alloc.systematic <- ifelse(dat.bcg$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=cbind(alloc.random, alloc.systematic, year, ablat),
data=dat.bcg, method="REML", btt=c(2,3))
### use model.matrix() to code the factor and set up the design matrix
### careful: X already includes the intercept, so need to use intercept=FALSE
X <- model.matrix(~ factor(alloc) + year + ablat, data=dat.bcg)
rma(yi, vi, mods=X, intercept=FALSE, data=dat.bcg, method="REML", btt=c(2,3))Run the code above in your browser using DataLab