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"),
                es.type = c("precalculated", "raw"),
                es.var = ifelse(identical(es.measure[1], "g"), 
                                ".g", ".log_rr"),
                se.var = ifelse(identical(es.measure[1], "g"), 
                                ".g_se", ".log_rr_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"),
                nnt.cer = 0.2,
                rho.within.study = 0.6,
                phi.within.study = 0.9,
                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), or using risk ratios ("RR")? Meta-analyses will only be conducted using comparisons that contain non NA 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".
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". The default settings correspond with the standard output of calculateEffectSizes().
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 automatically set to .g_se; if es.measure = "RR", ".log_rr_se" is used. The default settings correspond with the standard output of calculateEffectSizes().
es.binary.raw.vars: character vector, defining the column names 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, 2017), "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.
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. Please note that the name of the variable must be included as a column in data.
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. 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.
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.
nnt.cer: numeric. Value between 0 and 1, indicating the assumed control group event rate to be used for calculating NNTs via the Furukawa-Leucht method.
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 conditional 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".
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.
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".

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". Pools all effect sizes within one study (defined by study.var) before pooling. This ensures that all effect sizes are independent (i.e., unit-of-analysis error & double-counting is avoided). To combine the effects, one has to assume a correlation of effect sizes within studies, empirical estimates of which are typically not available.
"lowest.highest". Runs a meta-analysis, but with only (i) the lowest and (ii) highest effect size within each study included.
"outlier". Runs a meta-analysis without statistical outliers (i.e. effect sizes for which the confidence interval does not overlap with the confidence intervall 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. 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.

Internally, the overall, combined, lowest.highest, outlier, 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 and threelevel.che models are implemented using metafor::rma.mv() in {metafor} (Viechtbauer, 2005).

Outlier selection is implemented using the dmetar::find.outliers() function, and influence analyses using the dmetar::InfluenceAnalysis() function. The latter function is a wrapper for metafor::influence.rma.uni().

\(~\)

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" and "threelevel.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 studies' \(\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) {
data("depressionPsyCtr")
library(meta)

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

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

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

# Show summary
res

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

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

# 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)

# 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")

# Run meta-analysis using raw response rate data
data %>% 
  runMetaAnalysis(es.measure = "RR",
                  es.type = "raw")
}

Usage

Arguments

Value

Details

Simple or complex variance-covariance approximation

Replacement functions

Risk of Bias data

See also

Author

Examples