Run different types of meta-analyses

This wrapper function allows to simultaneously pool effect sizes using different meta-analytic approaches.

Usage

runMetaAnalysis(data,

                # Models to run
                which.run = c("overall", "combined",
                              "lowest.highest", "outliers",
                              "influence", "rob", "threelevel",
                              "threelevel.che"),
                              
                # Effect size measure
                es.measure = c("g", "RR", "EER", "CER", "ROM"),
                es.type = c("precalculated", "raw"),
                es.var = ifelse(identical(es.measure[1], "RR"), 
                                          ".log_rr", ".g"),
                se.var = ifelse(identical(es.measure[1], "RR"), 
                                          ".log_rr_se", ".g_se"),
                es.binary.raw.vars = 
                  c(".event_arm1", ".event_arm2",
                    ".totaln_arm1", ".totaln_arm2"),
                    
                # Estimator of the heterogeneity variance
                method.tau = "REML",
                method.tau.ci = "Q-Profile",
                i2.ci.boot = FALSE,
                nsim.boot = 5e3,
                hakn = TRUE,
                
                # Data specifications
                study.var = "study",
                arm.var.1 = "condition_arm1",
                arm.var.2 = "condition_arm2",
                measure.var = "instrument",
                low.rob.filter = "rob > 2",
                round.digits = 2,
                rob.data = NULL,
                
                # Model specifications
                which.combine = c("arms", "studies"),
                which.combine.var = "multi_arm1",
                which.outliers = c("overall", "combined"),
                which.influence = c("overall", "combined"),
                which.rob = c("overall", "combined"),
                which.waap.wls = c("overall", "combined"),
                nntCer = 0.2,
                rho.within.study = 0.6,
                phi.within.study = 0.9,
                power.within.study = 0.8,
                w1.var = ifelse(identical(es.measure[1], "g"), 
                                "n_arm1", "totaln_arm1"),
                w2.var = ifelse(identical(es.measure[1], "g"), 
                                "n_arm2", "totaln_arm2"),
                time.var = "time_weeks",
                vcov = c("simple", "complex"),
                near.pd = FALSE,
                use.rve = TRUE,
                
                # Output
                html = TRUE,
                
                # Additional arguments
                ...)

Arguments

data: data.frame. Effect size data in the wide format, as created by calculateEffectSizes.
which.run: character. Selection of models to be calculated. See 'Details'.
es.measure: character. Should meta-analyses be calculated using the bias-corrected standardized mean difference ("g"; default), risk ratios ("RR"), logit-transformed experimental/control group event rates ("EER" or "CER"), or ratio of means ("ROM")? Meta-analyses will only be conducted using comparisons that contain non-missing values in the es.var and se.var columns.
es.type: character. Should pre-calculated or raw event data (i.e. the Mantel-Haenszel method) be used for meta-analyses of risk ratios? Can be set to "precalculated" (default) or "raw". Ignored for all other values of es.measure.
es.var: character. Specifies the name of the variable containing the (pre- calculated) effect size data in data. When es.measure = "g", this is set to .g by default; ".log_rr" is used when es.measure = "RR". When es.measure = "EER" or "CER", logit-transformed event rates are calculated from es.binary.raw.vars and saved under this name for internal use (the returned data is not altered). For es.measure = "ROM", this defaults to .log_rom (log ratio of means) and can be set to any column name if the user stores the data elsewhere.
se.var: character. Specifies the name of the variable containing the (pre-calculated) standard errors (square root of the variance) of the effect size metric defined in es.var. If es.measure = "g", this is set to .g_se; if es.measure = "RR", ".log_rr_se" is used. For es.measure = "ROM", this defaults to .log_rom_se and can be set to any column name. The default settings correspond with the standard output of calculateEffectSizes.
es.binary.raw.vars: character. A vector defining the column names in data in which the (1) raw event counts in the experimental group, (2) raw event counts in the control/reference group, (3) sample size in the experimental group, and (4) sample size of the control/reference group are stored. Defaults correspond with the standard output of calculateEffectSizes.
method.tau: character. A character string indicating which method is used to estimate the between-study variance (tau-squared) and its square root (tau). Either "REML" (default), "DL", "PM", "ML", "HS", "SJ", "HE", or "EB", can be abbreviated (see metagen). Use "FE" to use a fixed-effect/"common effect" model.
method.tau.ci: character. A character string indicating which method is used to estimate the confidence interval of the between-study heterogeneity variance \(\tau^2\). Either "Q-Profile" (default and recommended; Viechtbauer, 2007), "BJ", "J", or "PL" can be abbreviated. See metagen and rma.uni for details.
i2.ci.boot: logical. Confidence intervals for \(\tau^2\) as calculated by the Q-Profile method are not directly applicable for three-level models, in which two heterogeneity variance components are estimated. By default, this argument is therefore set to FALSE, and no confidence intervals around \(\tau^2\) and \(I^2\) are provided for the "threelevel" and "threelevel.che" model. If this argument is set to TRUE, parametric bootstrapping is used to calculate confidence intervals around the between- and within-study heterogeneity estimates (\(\tau\) and \(I^2\)). Please note that this can take several minutes, depending on the number of effect sizes. If correctPublicationBias is used and i2.ci.boot is TRUE, bootstrapping will also be used to calculate confidence intervals around the \(G^2\) statistic (Rücker et al., 2011) used in the limit meta-analysis (note that \(G^2\) is printed as \(I^2\) in this package).
nsim.boot: numeric Number of bootstrap samples to be drawn when i2.ci.boot is TRUE. Defaults to 5000.
hakn: logical. Should the Knapp-Hartung adjustment for effect size significance tests be used? Default is TRUE. For multilevel models, TRUE means that tests are based on a t-distribution instead of the standard normal distribution; this has no effect on confidence intervals and tests when use.rve = TRUE, but is still used for the prediction interval.
study.var: character. The name of the variable in data in which the study IDs are stored.
arm.var.1: character. The name of the variable in data in which the condition (e.g. "guided iCBT") of the first arm within a comparison are stored.
arm.var.2: character. The name of the variable in data in which the condition (e.g. "wlc") of the second arm within a comparison are stored.
measure.var: character. The name of the variable in data in which the instrument used for the comparison is stored.
low.rob.filter: character. A filtering statement by which to include studies for the "low RoB only" analysis (e.g. "rob > 2", or "rob %in% c(3,4)"). The statement is evaluated within data, so all variable names must be included as columns in data. Rows for which the statement evaluates to NA are not included.
round.digits: numeric. Number of digits to round the (presented) results by. Default is 2.
rob.data: list. Optional list detailing how risk of bias data should be appended to the model (see Details).
which.combine: character. Should multiple effect sizes within one study be pooled on an "arms" (default) or "studies" level? When a study is a multi-arm trial, setting which.combine = "arms" will aggregate the effect sizes for each trial arm individually before pooling; the which.combine.var argument can be used to control which effects within a study should be aggregated. When which.combine = "studies", one overall aggregated effect is created for each study. This setting is preferable from a statistical perspective, since it ensures that all pooled effects can be assumed to be independent.
which.combine.var: character. Additional grouping variable within studies to be used for the "combined" analysis when which.combine = "arms". If the specified variable differs within one study (as defined by study.var), effects will be aggregated separately for each unique value in which.combine.var. If this variable is not included in data, effects are aggregated within studies. Defaults to "multi_arm1", the variable which encodes multi-arm intervention conditions in the Metapsy data standard.
which.outliers: character. Which model should be used to conduct outlier analyses? Must be "overall" or "combined", with "overall" being the default. If "overall" is not included in which.run, but "combined" is, the "combined" model is used. The same applies to which.influence, which.rob and which.waap.wls.
which.influence: character. Which model should be used to conduct influence analyses? Must be "overall" or "combined", with "overall" being the default.
which.rob: character. Which model should be used to conduct the "low risk of bias only" analyses? Must be "overall" or "combined", with "overall" being the default.
which.waap.wls: character. Which model should be used to run the (non-default) WAAP-WLS model? Must be "overall" or "combined", with "overall" being the default.
nntCer: numeric. Value between 0 and 1, indicating the assumed control group event rate to be used for calculating NNTs via the Furukawa-Leucht method. Pre-calculated risk ratios use the same event rate (\(NNT = 1/|CER(RR-1)|\)). When es.type = "raw", NNTs are based on the pooled risk difference, and this value is not used. NNTs are always shown as positive numbers, based on the absolute effect size, and thus do not show if an effect favors the first or second trial arm.
rho.within.study: numeric. Value between 0 and 1, indicating the assumed correlation of effect sizes within studies. This is relevant to combine effect sizes for the "combined" analysis type, and used to estimate the variance-covariance matrices needed for the correlated and hierarchical effects three-level model. Default is 0.6.
phi.within.study: numeric. Value between 0 and 1, indicating the assumed one-week autocorrelation of effect sizes. This is only used when vcov="complex" to approximate the variance-covariance matrices needed for the "combined" and "threelevel.che" model. Default is 0.9. See "Details".
power.within.study: numeric. Value between 0 and 1, indicating the minimum power (\(1-\beta\)) required for a study to be considered sufficiently powered. This value is used to identify adequately powered studies in the WAAP-WLS analysis. See "Details".
w1.var: character. Name of the variable in data in which the sample sizes of the first arm are stored. See "Details".
w2.var: character. Name of the variable in data in which the sample sizes of the second arm are stored. See "Details".
time.var: character. Name of the variable in data in which the assessment time point is stored. Should be expressed as weeks since randomization; other units (e.g. days, months) are also possible if phi.within.study is specified accordingly. See "Details".
vcov: character. For the "combined" and "threelevel.che" model, should variance-covariance matrices (representing the dependency structure of the data) be approximated using a heterogeneous compound symmetry ("simple"; default) or unstructured matrix structure ("complex")? See "Details".
near.pd: logical. If at least one of the study variance-covariance matrices constructed when vcov="complex" is not positive definite/invertible, should the nearPD function be used to compute the nearest positive definite matrix? Default is FALSE.
use.rve: logical. Should robust variance estimation be used to calculate confidence intervals and tests of three-level models? TRUE by default. The CR2 estimator with Satterthwaite degrees of freedom is used for the "threelevel" and "threelevel.che" models, and the CR1 estimator (clustered on study.var) for the "ccrem" and "ccrem.che" models.
html: logical. Should an HTML table be created for the results? Default is TRUE.
...: Additional arguments.

Value

Returns an object of class "runMetaAnalysis". This object includes, among other things, a data.frame with the name summary, in which all results are summarized - including the studies which were removed for some analysis steps. Other objects are the "raw" model objects returned by all selected analysis types. This allows to conduct further operations on some models specifically (e.g. run a meta-regression by plugging one of the model objects in meta:::update.meta.

Details

The runMetaAnalysis function is a wrapper for several types of meta-analytic models that are typically used. It allows to run all of these models in one step in order to generate results that are somewhat closer to being "publication-ready". See plot.runMetaAnalysis and createRobSummary for additional functionality.

By default, the following models are calculated:

"overall". Runs a generic inverse-variance (random-effects) model. All included effect sizes are treated as independent. When es.measure = "RR" and es.type = "raw", the Mantel-Haenszel method is used for pooling instead.
"combined". Aggregates the effect sizes within studies (defined by study.var) before pooling. To combine the effects, one has to assume a correlation of effect sizes within studies, empirical estimates of which are typically not available. By default (which.combine = "arms"), the different arms of a multi-arm trial are aggregated separately, as defined by which.combine.var. All effect sizes are only independent (i.e., unit-of-analysis error & double-counting is avoided) when which.combine = "studies".
"lowest.highest". Runs a meta-analysis, but with only (i) the lowest and (ii) highest effect size within each study included.
"outliers". Runs a meta-analysis without statistical outliers (i.e. effect sizes for which the confidence interval does not overlap with the confidence interval of the overall effect).
"influence". Runs a meta-analysis without influential cases (see influence.rma.uni for details).
"rob". Runs a meta-analysis with only low-RoB studies included.
"threelevel". Runs a multilevel (three-level) meta-analysis model, with effect sizes nested in studies.
"threelevel.che". Runs a multilevel (three-level) meta-analysis model, with effect sizes nested in studies. Variance-covariance matrices of each study with two or more effect sizes are estimated using rho.within.study as the assumed overall within-study correlation, and the average sampling variance within a study for all of its effect sizes. This imputation allows to run a "correlated and hierarchical effects" (CHE) model, which is typically a good approximation for data sets with unknown and/or complex dependence structures.

The following non-default models can also be specified:

"ccrem". Runs a multilevel (three-level) cross-classified random effects model (CCREM), with measure.var modeled as a crossed random effect. Such models can be more appropriate when, across studies, effect sizes are nested in various instruments, subscales, or outcomes; and if we want to estimate (and generalize over) the variation due different instruments (or subscales, outcomes; Fernandez-Castilla et al., 2018).
"ccrem.che". Runs a multilevel (three-level) cross-classified random effects model (CCREM), with measure.var modeled as a crossed random effect. Additionally, variance-covariance matrices of each study with two or more effect sizes are approximated using rho.within.study as the assumed overall within-study correlation, similar to the "correlated and hierarchical effects" (CHE, threelevel.che) model.
"waap.wls". Runs a weighted average of adequately powered studies (WAAP) weighted least squares model (WLS). If at least 3 effect sizes are available, effects are only computed among adequately powered studies, which are determined using the pooled common-effect estimate as basis. If less than 3 studies show adequate power, a simple WLS model among all included effect sizes is returned. This type of analysis employing multiplicative error models has been proposed to obtain more robust effect estimates that guard against common biases introduced by selective publication and QRPs (Stanley, Doucouliagos & Ioannidis, 2017; Stanley & Doucouliagos, 2017; Carter et al., 2019)

Internally, the overall, combined, lowest.highest, outliers, influence and rob models are fitted by calling the meta::metagen() or meta::metabin() function, respectively, in {meta} (Balduzzi, Rücker & Schwarzer, 2019). The threelevel, threelevel.che, ccrem, and ccrem.che models are implemented using metafor::rma.mv() in {metafor} (Viechtbauer, 2010).

Influence analyses are based on metafor::influence.rma.uni().

Note that the \(I^2\) values of models fitted with {meta} are based on Cochran's \(Q\), while the three-level models and selection models use the "generalized" definition based on the estimated heterogeneity variance and a "typical" sampling variance. The values can therefore differ somewhat, even when based on the same data.

For event rates (es.measure = "EER" or "CER"), a continuity correction of 0.5 is added to the number of events and non-events in studies with an event rate of 0 or 1.

For ratio of means (es.measure = "ROM"), data must contain the columns specified by es.var and se.var (default .log_rom and .log_rom_se, as created by calculateEffectSizes when calculate.rom = TRUE). Custom es.var and se.var can be used if log-ROM and SE are stored in other columns. Values are used internally and results are back-transformed (exponentiated) for display.

\(~\)

Simple or complex variance-covariance approximation

The vcov argument controls if the effect size dependencies within the data should be approximated using a "simple" (default) or more "complex" (but potentially more accurate) method. This argument is only relevant for the "combined", "threelevel.che", and "ccrem.che" models. The default "simple" method constructs variance-covariance matrices \(\Sigma_k\) for each study using a constant sampling correlation \(\rho\) (defined by rho.within.study), which is identical across all studies, outcomes, and time points. This simplifying assumption is part of the formulation of the CHE model originally provided by Pustejovsky and Tipton (2022).

Naturally, employing a common value of \(\rho\) across all studies may not be reasonable in some analyses, and other information may be available to better approximate the effect size dependencies in the collected data. Setting vcov to "complex" allows to assume that correlations between effect sizes may differ conditional on the type of dependency. This means that the variance-covariance matrix \(\Sigma_k\) of some study \(k\) is approximated by an unstructured matrix with varying \(\rho_{ij}\) (instead of a heterogeneous compound symmetry matrix with fixed \(\rho\), as is used when vcov="simple").

\[\begin{array}{ccc}\texttt{vcov="simple"} & \texttt{vcov="complex"} & \\ \Sigma_k = \begin{bmatrix} \sigma^2_1 \\ \rho \sigma_2 \sigma_1 & \sigma^2_2 & & \\ \rho \sigma_3 \sigma_1 & \rho \sigma_3 \sigma_2 & \sigma^2_3 & \\ \rho \sigma_4 \sigma_1 & \rho \sigma_4 \sigma_2 & \rho \sigma_4 \sigma_3 & \sigma^2_4 \end{bmatrix} & \Sigma_k = \begin{bmatrix} \sigma^2_1 & & & \\ \rho_{21} \sigma_2 \sigma_1 & \sigma^2_2 & & \\ \rho_{31} \sigma_3 \sigma_1 & \rho_{32} \sigma_3 \sigma_2 & \sigma^2_3 & \\ \rho_{41} \sigma_4 \sigma_1 & \rho_{42} \sigma_4 \sigma_2 & \rho_{43} \sigma_4 \sigma_3 & \sigma^2_4 \end{bmatrix} & \end{array}\]

For example, setting vcov = "complex" allows to additionally incorporate assumed correlations specific to multiple testing over time (e.g. correlations between effects at post-test and long-term follow-up). The value provided in phi.within.study represents the (auto-)correlation coefficient \(\phi\), which serves as a rough estimate of the re-test correlation after 1 week. When a vector of follow-up lengths is provided in time.var, this allows to model the gradual decrease in correlation between measurements over time. Furthermore, it is possible to calculate a correlation coefficient \(\rho_w\) for multi-arm trials, which is directly proportional to the size of each individual trial arm. When all trial arms have the same size, meaning that each arm's weight \(w\) is identical, \(\rho_w\) is known to be 0.5. Multiarm weights \(w\) (and thus \(\rho_w\)) can be derived if the w1.var and w2.var variables, containing the sample size of each study arm, are provided.

Using the complex approximation method increases the risk that at least one study's \(\Sigma_k\) matrix is not positive definite. In this case, the function automatically switches back to the constant sampling correlation approximation.

\(~\)

Replacement functions

Once a model has been fitted using runMetaAnalysis, replacement functions are defined for each function argument. This allows to quickly tweak one or more analysis settings, which are implemented once the rerun function is called. Say that we saved the results of runMetaAnalysis in an object m. If, for example, we want to check the results using a different estimator of \(\tau^2\), leaving all other settings the same, we could run e.g. method.tau(m) <- "PM", followed by rerun(m). This would provide results using the Paule-Mandel estimator. A list of all available setting replacement functions is provided here.

\(~\)

Risk of Bias data

Using the rob.data argument, it is possible to specify which variables in the data set provided in data contain Risk of Bias (ROB) assessment information. If specified, this allows to generate forest plots with added ROB information when using plot.runMetaAnalysis. ROB information added this way can also be used to create ROB summary plots using createRobSummary.

The object provided to rob.data must be a list element with these elements:

domains: A vector of characters, specifying variables that contain ratings for different ROB domains (e.g. allocation concealment, missing data handling, ...).
domain.names (Optional): A vector of characters with the same length as domains, which provides long-format labels for each included domain.
overall.rob (Optional): A single character specifying the variable in data that contains the overall ROB rating.
categories: A vector of characters, specifying how ROB ratings were coded in the selected variables (e.g., "low, "high", "unclear").
symbols (Optional): A vector of single-letter characters (or symbols) that should be used when plotting ROB rating in the forest plot.
colors (Optional): A vector of characters of the same length as categories, specifying colors for each rating code.

For concrete usage examples with this functionality, see "Examples".

For more details see the Get Started vignette.

Author

Mathias Harrer mathias.h.harrer@gmail.com, Paula Kuper paula.r.kuper@gmail.com, Pim Cuijpers p.cuijpers@vu.nl

Examples

if (FALSE) { # \dontrun{
data("depressionPsyCtr")
library(meta)

depressionPsyCtr %>%
  checkDataFormat() %>%
  checkConflicts() %>%
  calculateEffectSizes() %>% 
  filterPoolingData(condition_arm2 %in% 
                      c("wl", "other ctr")) -> data

# Run the meta-analyses
runMetaAnalysis(data) -> res

# Check if variance components of "threelevel.che" model
# are identifiable
profile(res, "threelevel.che")

# Run non-default models
# - Cross-classified random effects model
runMetaAnalysis(data, 
                which.run = c("ccrem", "ccrem.che"), 
                vcov = "complex")

# - Weighted average of adequately powered studies (WAAP-WLS)
runMetaAnalysis(data, which.run = "waap.wls")
runMetaAnalysis(data, which.run = "waap.wls", power.within.study = .9)

# - Meta-analysis using raw response rate data
runMetaAnalysis(data, 
                es.measure = "RR",
                es.type = "raw")

# - Estimate intervention response rates, then pool 
data %>% 
  calculateEffectSizes(impute.response = TRUE) %>% 
  runMetaAnalysis(which.run = "combined", 
                  es.measure = "EER")

# - Run a meta-analysis of ratio-of-means (RoMs)
data %>% 
  calculateEffectSizes(calculate.rom = TRUE) %>% 
  runMetaAnalysis(which.run = c("combined", "threelevel"),
                  es.measure = "ROM")

# Use replacement function to show results for
# differing settings
method.tau(res) <- "PM"
hakn(res) <- FALSE
rerun(res)

# Show forest plot (by default, "overall" is used)
plot(res)

# Compare effects across models
plot(res, "summary")

# Show forest plot of specific analysis
plot(res, "outliers")
plot(res, "threelevel")
plot(res, "baujat")
plot(res, "influence")

# Show forest plot with empirical Bayes estimates
plot(res, "combined", eb = TRUE)

# Extract specific model and do further calculations
# (e.g. meta-regression on 'year')
metaRegression(res$model.overall, ~ scale(year))

# Conduct a subgroup analysis
subgroupAnalysis(res, country)

# Correct for publication bias/small-study effects
(correctPublicationBias(res) -> res.pb)
plot(res.pb, "trimfill")
plot(res.pb, "limitmeta")
plot(res.pb, "selection")

# For the combined analysis, set which.combine to
# "studies" here, so that all effects in a study are aggregated
# first before pooling
data %>% 
  runMetaAnalysis(which.combine = "studies") %>% 
  plot("combined")

# Define ROB data to be added to the models
robData <- list(
  # Names of ROB variables included in 'data'
  domains = c("sg", "ac", "ba", "itt"),
  # Long-format labels for each ROB domain
  domain.names = c("Sequence Generation", 
                   "Allocation Concealment", 
                   "Blinding of Assessors", 
                   "ITT Analyses"),
  # Codes used to rate the risk of bias (sr=self-report)
  categories = c("0", "1", "sr"),
  # Symbols that should be used for these codes in forest plots
  symbols = c("-", "+", "s"),
  # Colors to be used in forest plots for each of these codes
  colors = c("red", "green", "yellow"))

# Re-run model with appended ROB data
res <- runMetaAnalysis(data, rob.data = robData) 

# Generate forest plot with ROB data
plot(res, "combined")

# Create a summary plot
createRobSummary(res, 
                 name.low = "1", 
                 name.high = "0", 
                 name.unclear = "sr",
                 which.run = "combined")
} # }

Usage

Arguments

Value

Details

Simple or complex variance-covariance approximation

Replacement functions

Risk of Bias data

See also

Author

Examples