orchaRd (version 2.2.1)

mod_results: mod_results

Description

Using a metafor model object of class rma or rma.mv, this function creates a table of model results containing the mean effect size estimates for all levels of a given categorical moderator, and their corresponding confidence and prediction intervals. The function is capable of calculating marginal means from meta-regression models, including those with multiple moderator variables of mixed types (i.e. continuous and categorical variables).

Usage

mod_results(
  model,
  mod = "1",
  group,
  N = NULL,
  weights = "prop",
  by = NULL,
  at = NULL,
  subset = FALSE,
  upper = TRUE,
  ...
)

Value

A data frame containing all the model results including mean effect size estimate, confidence and prediction intervals

Arguments

model

rma.mv model object

mod

Moderator variable of interest that one wants marginal means for. Defaults to the intercept, i.e. "1".

group

The grouping variable that one wishes to plot beside total effect sizes, k. This could be study, species, or any grouping variable one wishes to present sample sizes for.

N

The name of the column in the data specifying the sample size so that each effect size estimate is scaled to the sample size, N. Defaults to NULL, so that precision is used for scaling each raw effect size estimate instead of sample size.

weights

How to marginalize categorical variables. The default is weights = "prop", which weights moderator level means based on their proportional representation in the data. For example, if "sex" is a moderator, and males have a larger sample size than females, then this will produce a weighted average, where males are weighted more towards the mean than females. This may not always be ideal. In the case of sex, for example, males and females are roughly equally prevalent in a population. As such, you can give the moderator levels equal weight using weights = "equal".

by

Character vector indicating the name that predictions should be conditioned on for the levels of the moderator.

at

List of levels one wishes to predict at for the corresponding variables in by. Used when one wants marginalised means. This argument can also be used to suppress levels of the moderator when argument subset = TRUE. Provide a list as follows: list(mod = c("level1", "level2")).

subset

Used when one wishes to only plot a subset of levels within the main moderator of interest defined by mod. Default is FALSE, but use TRUE if you wish to subset levels of a moderator plotted (defined by mod) for plotting. Levels one wishes to plot are specified as a list, with the level names as a character string in the at argument. For subsetting to work, the at argument also needs to be specified so that the mod_results function knows what levels one wishes to plot.

upper

Logical, defaults to TRUE, indicating that the first letter of the character string for the moderator variable should be capitalized.

...

Additional arguments passed to emmeans::emmeans().

Author

Shinichi Nakagawa - [e-mail]

Daniel Noble - [e-mail]

Details

Prediction intervals with random-slope models (struct = "GEN" or "HCS"):

When your rma.mv model includes a random-slope term (e.g., random = ~ moderator | study with struct = "GEN" or "HCS"), prediction intervals require the full variance-covariance matrix of the random effects — not just the sum of the variance components. In a random-slope model the random intercept and slope are typically correlated, and a correct prediction interval at a given moderator value x must incorporate:

$$Var(u_{0j} + u_{1j} x) = \tau^2_0 + 2 x \rho \tau_0 \tau_1 + x^2 \tau^2_1$$

where \(\tau^2_0\) and \(\tau^2_1\) are the random intercept and slope variances, and \(\rho\) is their correlation. Simply summing tau2, gamma2, and sigma2 ignores the slope variance contribution and the covariance, producing intervals that may be too narrow or too wide depending on the moderator value.

metafor::predict.rma() does not currently return prediction intervals for models with GEN or HCS structures for this reason. orchaRd will issue a warning when it detects such a structure. Users should verify prediction intervals manually using the approach described by Pustejovsky (2024, R-sig-meta-analysis mailing list): https://stat.ethz.ch/pipermail/r-sig-meta-analysis/2024-March/005152.html.

A practical workaround is to re-centre the moderator at the value of interest (so the random slope term evaluates at zero) and then compare the resulting prediction interval with a manual calculation using the full variance-covariance matrix from the model.

Examples

# \donttest{
# Simple eklof data
data(eklof)
eklof<-metafor::escalc(measure="ROM", n1i=N_control, sd1i=SD_control,
m1i=mean_control, n2i=N_treatment, sd2i=SD_treatment, m2i=mean_treatment, data = eklof)
# Add the unit level predictor
eklof$Datapoint<-as.factor(seq(1, dim(eklof)[1], 1))
# fit a MLMR - accouting for some non-independence
eklof_MR<-metafor::rma.mv(yi=yi, V=vi, mods=~ Grazer.type, random=list(~1|ExptID,
~1|Datapoint), data = eklof)
results <- mod_results(eklof_MR, mod = "Grazer.type", group = "ExptID")

# Fish example demonstrating marginalised means
data(fish)
model <- metafor::rma.mv(yi = lnrr, V = lnrr_vi,
  random = list(~1 | group_ID, ~1 | es_ID),
  mods = ~ trait.type + deg_dif,
  method = "REML", test = "t", data = fish)
overall <- mod_results(model, group = "group_ID")
across_trait <- mod_results(model, group = "group_ID", mod = "trait.type")
# Marginalised means, conditioning on levels of a continuous moderator
across_trait_by_deg <- mod_results(model, group = "group_ID",
  mod = "trait.type", at = list(deg_dif = c(5, 10, 15)), by = "deg_dif")
# }