emmprep.RdFunction to create a reference grid for use with the emmeans function from the package of the same name.
emmprep(x, verbose=FALSE, ...)an object of class "rma".
logical to specify whether information on some (extracted) settings should be printed when creating the reference grid (the default is FALSE).
other arguments that will be passed on to the qdrg function.
The emmeans package is a popular package that facilitates the computation of 'estimated marginal means'. The function is a wrapper around the qdrg function from the emmeans package to make "rma" objects compatible with the latter. Unless one needs to pass additional arguments to the qdrg function, one simply applies this function to the "rma" object and then the emmeans function (or one of the other functions that can be applied to "emmGrid" objects) to the resulting object to obtain the desired estimated marginal means.
An "emmGrid" object as created by the qdrg function from the emmeans package.
The resulting object will typically be used in combination with the emmeans function.
When creating the reference grid, the function extracts the degrees of freedom for tests/confidence intervals from the model object (if the model was fitted with test="t", test="knha", test="hksj", or test="adhoc"; otherwise the degrees of freedom are infinity). In some cases, there is not just a single value for the degrees of freedom, but an entire vector (e.g., for models fitted with rma.mv). In this case, the smallest value will be used (as a conservative option). One can set a different/custom value for the degrees of freedom with emmprep(..., df=value).
When the model object contains information about the outcome measure used in the analysis (which should be the case if the observed outcomes were computed with escalc or if the measure argument was set when fitting the model), then information about the appropriate back-transformation (if available) is stored as part of the returned object. If so, the back-transformation is automatically applied when calling emmeans with type="response".
The function also tries to extract the estimated value of \(\tau^2\) (or more precisely, its square root) from the model object (when the model is a random/mixed-effects model). This value is only needed when computing prediction intervals (i.e., when interval="predict" in predict.emmGrid) or when applying the bias adjustment in the back-transformation (i.e., when bias.adjust=TRUE in summary.emmGrid). For some models (e.g., those fitted with rma.mv), it is not possible to automatically extract the estimate. In this case, one can manually set the value with emmprep(..., sigma=value) (note: the argument is called sigma, following the conventions of summary.emmGrid and one must supply the square root of the \(\tau^2\) estimate).
By default, the reference grid is created based on the data used for fitting the original model (which is typically the sensible thing to do). One can specify a different dataset with emmprep(..., data=obj), where obj must be a data frame that contains the same variables as used in the original model fitted.
If the original model fitted involved redundant predictors that were dropped from the model (due to ‘rank deficiencies’), then the function cannot be used. In this case, one should remove any redundancies in the original model fitted before using this function.
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
### calculate log risk ratios and corresponding sampling variances
dat <- escalc(measure="RR", ai=tpos, bi=tneg, ci=cpos, di=cneg, data=dat.bcg)
### fit meta-regression model with absolute latitude as predictor
res <- rma(yi, vi, mods = ~ ablat, data=dat)
res
#>
#> Mixed-Effects Model (k = 13; tau^2 estimator: REML)
#>
#> tau^2 (estimated amount of residual heterogeneity): 0.0764 (SE = 0.0591)
#> tau (square root of estimated tau^2 value): 0.2763
#> I^2 (residual heterogeneity / unaccounted variability): 68.39%
#> H^2 (unaccounted variability / sampling variability): 3.16
#> R^2 (amount of heterogeneity accounted for): 75.62%
#>
#> Test for Residual Heterogeneity:
#> QE(df = 11) = 30.7331, p-val = 0.0012
#>
#> Test of Moderators (coefficient 2):
#> QM(df = 1) = 16.3571, p-val < .0001
#>
#> Model Results:
#>
#> estimate se zval pval ci.lb ci.ub
#> intrcpt 0.2515 0.2491 1.0095 0.3127 -0.2368 0.7397
#> ablat -0.0291 0.0072 -4.0444 <.0001 -0.0432 -0.0150 ***
#>
#> ---
#> Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#>
### create reference grid
sav <- emmprep(res, verbose=TRUE)
#>
#> Extracted formula: ~ ablat
#> Value of tau^2: 0.0764
#> Transformation: log
#>
### estimated marginal mean (back-transformed to the risk ratio scale)
if (require(emmeans))
emmeans(sav, specs="1", type="response")
#> Loading required package: emmeans
#> Welcome to emmeans.
#> Caution: You lose important information if you filter this package's results.
#> See '? untidy'
#> 1 response SE df asymp.LCL asymp.UCL
#> overall 0.486 0.0523 Inf 0.393 0.6
#>
#> Confidence level used: 0.95
#> Intervals are back-transformed from the log scale
### same as the predicted effect at the mean absolute latitude
predict(res, newmods = mean(model.matrix(res, asdf=TRUE)$ablat), transf=exp, digits=3)
#>
#> pred ci.lb ci.ub pi.lb pi.ub
#> 0.486 0.393 0.600 0.272 0.868
#>
### fit meta-regression model with allocation factor
res <- rma(yi, vi, mods = ~ alloc, data=dat)
res
#>
#> Mixed-Effects Model (k = 13; tau^2 estimator: REML)
#>
#> tau^2 (estimated amount of residual heterogeneity): 0.3615 (SE = 0.2111)
#> tau (square root of estimated tau^2 value): 0.6013
#> I^2 (residual heterogeneity / unaccounted variability): 88.77%
#> H^2 (unaccounted variability / sampling variability): 8.91
#> R^2 (amount of heterogeneity accounted for): 0.00%
#>
#> Test for Residual Heterogeneity:
#> QE(df = 10) = 132.3676, p-val < .0001
#>
#> Test of Moderators (coefficients 2:3):
#> QM(df = 2) = 1.7675, p-val = 0.4132
#>
#> Model Results:
#>
#> estimate se zval pval ci.lb ci.ub
#> intrcpt -0.5180 0.4412 -1.1740 0.2404 -1.3827 0.3468
#> allocrandom -0.4478 0.5158 -0.8682 0.3853 -1.4588 0.5632
#> allocsystematic 0.0890 0.5600 0.1590 0.8737 -1.0086 1.1867
#>
#> ---
#> Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#>
### create reference grid
sav <- emmprep(res)
### estimated marginal mean using proportional cell weighting
if (require(emmeans))
emmeans(sav, specs="1", type="response", weights="proportional")
#> 1 response SE df asymp.LCL asymp.UCL
#> overall 0.481 0.092 Inf 0.331 0.7
#>
#> Results are averaged over the levels of: alloc
#> Confidence level used: 0.95
#> Intervals are back-transformed from the log scale
### estimated marginal mean using equal cell weighting (this is actually the default)
if (require(emmeans))
emmeans(sav, specs="1", type="response", weights="equal")
#> 1 response SE df asymp.LCL asymp.UCL
#> overall 0.529 0.109 Inf 0.352 0.793
#>
#> Results are averaged over the levels of: alloc
#> Confidence level used: 0.95
#> Intervals are back-transformed from the log scale
### same as the predicted effect using cell proportions as observed in the data
### or using equal proportions for the three groups
predict(res, newmods = colMeans(model.matrix(res))[-1], transf=exp, digits=3)
#>
#> pred ci.lb ci.ub pi.lb pi.ub
#> 0.481 0.331 0.700 0.140 1.657
#>
predict(res, newmods = c(1/3,1/3), transf=exp, digits=3)
#>
#> pred ci.lb ci.ub pi.lb pi.ub
#> 0.529 0.352 0.793 0.152 1.838
#>
### fit meta-regression model with absolute latitude and allocation as predictors
res <- rma(yi, vi, mods = ~ ablat + alloc, data=dat)
res
#>
#> Mixed-Effects Model (k = 13; tau^2 estimator: REML)
#>
#> tau^2 (estimated amount of residual heterogeneity): 0.1446 (SE = 0.1124)
#> tau (square root of estimated tau^2 value): 0.3803
#> I^2 (residual heterogeneity / unaccounted variability): 70.11%
#> H^2 (unaccounted variability / sampling variability): 3.35
#> R^2 (amount of heterogeneity accounted for): 53.84%
#>
#> Test for Residual Heterogeneity:
#> QE(df = 9) = 26.2034, p-val = 0.0019
#>
#> Test of Moderators (coefficients 2:4):
#> QM(df = 3) = 11.0605, p-val = 0.0114
#>
#> Model Results:
#>
#> estimate se zval pval ci.lb ci.ub
#> intrcpt 0.2932 0.4050 0.7239 0.4691 -0.5006 1.0870
#> ablat -0.0273 0.0092 -2.9650 0.0030 -0.0453 -0.0092 **
#> allocrandom -0.2675 0.3504 -0.7633 0.4453 -0.9543 0.4193
#> allocsystematic 0.0585 0.3795 0.1540 0.8776 -0.6854 0.8023
#>
#> ---
#> Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#>
### create reference grid
sav <- emmprep(res)
### estimated marginal mean using equal cell weighting
if (require(emmeans))
emmeans(sav, specs="1", type="response")
#> 1 response SE df asymp.LCL asymp.UCL
#> overall 0.502 0.0718 Inf 0.379 0.664
#>
#> Results are averaged over the levels of: alloc
#> Confidence level used: 0.95
#> Intervals are back-transformed from the log scale
### same as the predicted effect at the mean absolute latitude and using equal proportions
### for the allocation factor
predict(res, newmods = c(mean(model.matrix(res, asdf=TRUE)$ablat),1/3,1/3), transf=exp, digits=3)
#>
#> pred ci.lb ci.ub pi.lb pi.ub
#> 0.502 0.379 0.664 0.226 1.113
#>
### create reference grid with ablat set equal to 10, 30, and 50 degrees
sav <- emmprep(res, at=list(ablat=c(10,30,50)))
### estimated marginal means at the three ablat values
if (require(emmeans))
emmeans(sav, specs="1", by="ablat", type="response")
#> ablat = 10:
#> 1 response SE df asymp.LCL asymp.UCL
#> overall 0.952 0.2300 Inf 0.593 1.529
#>
#> ablat = 30:
#> 1 response SE df asymp.LCL asymp.UCL
#> overall 0.552 0.0784 Inf 0.418 0.729
#>
#> ablat = 50:
#> 1 response SE df asymp.LCL asymp.UCL
#> overall 0.320 0.0712 Inf 0.207 0.495
#>
#> Results are averaged over the levels of: alloc
#> Confidence level used: 0.95
#> Intervals are back-transformed from the log scale
### same as the predicted effect at the chosen absolute latitude values and using equal
### proportions for the allocation factor
predict(res, newmods = cbind(c(10,30,50),1/3,1/3), transf=exp, digits=3)
#>
#> pred ci.lb ci.ub pi.lb pi.ub
#> 1 0.952 0.593 1.529 0.394 2.303
#> 2 0.552 0.418 0.729 0.249 1.223
#> 3 0.320 0.207 0.495 0.135 0.758
#>
############################################################################
### copy data into 'dat' and examine data
dat <- dat.mcdaniel1994
head(dat)
#> study ni ri type struct
#> 1 1 123 0.00 j s
#> 2 2 95 0.06 p u
#> 3 3 69 0.36 j s
#> 4 4 1832 0.15 j s
#> 5 5 78 0.14 j s
#> 6 6 329 0.06 j s
### calculate r-to-z transformed correlations and corresponding sampling variances
dat <- escalc(measure="ZCOR", ri=ri, ni=ni, data=dat)
### mixed-effects model with interview type as factor
res <- rma(yi, vi, mods = ~ factor(type), data=dat, test="knha")
#> Warning: 3 studies with NAs omitted from model fitting.
res
#>
#> Mixed-Effects Model (k = 157; tau^2 estimator: REML)
#>
#> tau^2 (estimated amount of residual heterogeneity): 0.0282 (SE = 0.0049)
#> tau (square root of estimated tau^2 value): 0.1681
#> I^2 (residual heterogeneity / unaccounted variability): 79.62%
#> H^2 (unaccounted variability / sampling variability): 4.91
#> R^2 (amount of heterogeneity accounted for): 1.92%
#>
#> Test for Residual Heterogeneity:
#> QE(df = 154) = 738.4411, p-val < .0001
#>
#> Test of Moderators (coefficients 2:3):
#> F(df1 = 2, df2 = 154) = 2.6775, p-val = 0.0719
#>
#> Model Results:
#>
#> estimate se tval df pval ci.lb ci.ub
#> intrcpt 0.2474 0.0196 12.6425 154 <.0001 0.2088 0.2861 ***
#> factor(type)p -0.1228 0.0608 -2.0210 154 0.0450 -0.2429 -0.0028 *
#> factor(type)s 0.0573 0.0625 0.9176 154 0.3603 -0.0661 0.1807
#>
#> ---
#> Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#>
### create reference grid
sav <- emmprep(res, verbose=TRUE)
#>
#> Extracted formula: ~ factor(type)
#> Degrees of freedom: 154
#> Value of tau^2: 0.0282
#> Transformation: r-to-z
#>
### estimated marginal mean (back-transformed to the correlation scale)
if (require(emmeans))
emmeans(sav, specs="1", type="response")
#> 1 response SE df lower.CL upper.CL
#> overall 0.222 0.0269 154 0.168 0.274
#>
#> Results are averaged over the levels of: type
#> Confidence level used: 0.95
#> Intervals are back-transformed from the r-to-z scale
### same as the predicted correlation using equal cell proportions
predict(res, newmods = c(1/3,1/3), transf=transf.ztor, digits=3)
#>
#> pred ci.lb ci.ub pi.lb pi.ub
#> 0.222 0.168 0.274 -0.111 0.510
#>
### estimated marginal means for the three interview types
if (require(emmeans))
emmeans(sav, specs="type", type="response")
#> type response SE df lower.CL upper.CL
#> j 0.243 0.0184 154 0.2058 0.279
#> p 0.124 0.0567 154 0.0109 0.234
#> s 0.296 0.0541 154 0.1854 0.399
#>
#> Confidence level used: 0.95
#> Intervals are back-transformed from the r-to-z scale
### same as the predicted correlations
predict(res, newmods = rbind(c(0,0), c(1,0), c(0,1)), transf=transf.ztor, digits=3)
#>
#> pred ci.lb ci.ub pi.lb pi.ub
#> 1 0.243 0.206 0.279 -0.087 0.524
#> 2 0.124 0.011 0.234 -0.223 0.443
#> 3 0.296 0.185 0.399 -0.047 0.576
#>
### illustrate use of the 'df' and 'sigma' arguments
res <- rma.mv(yi, vi, mods = ~ factor(type), random = ~ 1 | study,
data=dat, test="t", dfs="contain")
#> Warning: 3 rows with NAs omitted from model fitting.
res
#>
#> Multivariate Meta-Analysis Model (k = 157; method: REML)
#>
#> Variance Components:
#>
#> estim sqrt nlvls fixed factor
#> sigma^2 0.0282 0.1681 157 no study
#>
#> Test for Residual Heterogeneity:
#> QE(df = 154) = 738.4411, p-val < .0001
#>
#> Test of Moderators (coefficients 2:3):
#> F(df1 = 2, df2 = 154) = 2.9228, p-val = 0.0568
#>
#> Model Results:
#>
#> estimate se tval df pval ci.lb ci.ub
#> intrcpt 0.2474 0.0187 13.2089 154 <.0001 0.2104 0.2844 ***
#> factor(type)p -0.1228 0.0582 -2.1115 154 0.0363 -0.2377 -0.0079 *
#> factor(type)s 0.0573 0.0598 0.9587 154 0.3392 -0.0608 0.1754
#>
#> ---
#> Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#>
### create reference grid
sav <- emmprep(res, verbose=TRUE, df=154, sigma=0.1681)
#>
#> Extracted formula: ~ factor(type)
#> Degrees of freedom: 154
#> Value of tau^2: 0.0283
#> Transformation: r-to-z
#>
### estimated marginal mean (back-transformed to the correlation scale)
if (require(emmeans))
emmeans(sav, specs="1", type="response")
#> 1 response SE df lower.CL upper.CL
#> overall 0.222 0.0258 154 0.17 0.272
#>
#> Results are averaged over the levels of: type
#> Confidence level used: 0.95
#> Intervals are back-transformed from the r-to-z scale