Using Your Own Data

If you want, you can conduct the examples below with your own data. Note that the interactive questions won't work in this case (the correct answers are based on the curry data).

If you want to use your own data, you can use the function check_data() to help make sure that your data are suitable for the tutorial examples. Install the function by running:

source("https://raw.githubusercontent.com/cjvanlissa/meta_workshop/refs/heads/master/check_function.R")

Then, apply the function to your data object by running check_data(your_data_object). It will give helpful suggestions to make it easier to follow the tutorial with your own data:

#> ✔ Looks like you're all set to do the workshop!

df <- curry
check_data(df)

#> There is no column named 'yi' in your data. If you have an effect size column, it will be easier to do the tutorial if you rename it, using syntax like this:
#>   
#>   names(df)[which(names(df) == 'youreffectsize')] <- 'yi'
#>   
#>   If you do not yet have an effect size column, you may need to compute it first. Run ?metafor::escalc to see the help for this function.

Note that the message suggests that we rename the effect size column, which is called d in the curry dataset, to yi. Let's do that:

names(df)[which(names(df) == 'd')] <- 'yi'

Using Different Demo Data

If you do not have your own data but are up for a challenge, you can use a different demo dataset to complete the tutorial. Two suggestions are:

metaforest::fukkink_lont
pema::bonapersona

Fixed-effects Model with rma

The rma() function has many arguments, all of which you can accessed by typing ?rma in your console once the metafor package is loaded, or selecting the function and pressing F1.

Here is a table with the most important arguments for our code:

Let's conduct our first fixed-effects-model Meta-Analysis. We we will give the results of this analysis the simple name m.

m <- rma(yi = df$yi,     # The yi-column of the df, which contains Cohen's d
         vi = df$vi,    # The vi-column of the df, which contains the variances
         method = "FE") # Run a fixed-effect model
m

#> Fixed-Effects Model (k = 56)
#> 
#> I^2 (total heterogeneity / total variability):   64.95%
#> H^2 (total variability / sampling variability):  2.85
#> 
#> Test for Heterogeneity:
#> Q(df = 55) = 156.9109, p-val < .0001
#> 
#> Model Results:
#> 
#> estimate      se    zval    pval   ci.lb   ci.ub      
#>   0.2059  0.0219  9.4135  <.0001  0.1630  0.2487  *** 
#> 
#> ---
#> Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

We now see a summary of the results of our Meta-Analysis, including

• The total number of included studies (k)
• The overall effect, its confidence interval and p-value
• The Q-test of heterogeneity

Testing Moderators' Significance

To evaluate the statistical significance of a predictor, we conduct a t-test (or Z-test) of its \(\beta\)-weight.

\[ t=\frac{\beta}{SE_{\beta}}\]

Which provides a \(p\)-value telling us if a variable significantly predicts effect size differences in our regression model.

Assessing Meta-Regression Model Fit

In conventional regression, \(R^2\) is commonly used to quantify the percentage of variance in the data explained by the model, as a percentage of total variance (0-100%). As this measure is commonly used, and many researchers know how to to interpret it, we can also calculate a \(R^2\) analog for meta-regression using this formula:

\[R_2=\frac{\hat\tau^2_{REM}-\hat\tau^2_{MEM}}{\hat\tau^2_{REM}}\]

Where \(\hat\tau^2_{REM}\) is the estimated total heterogenetiy based on the random-effects-model and \(\hat\tau^2_{REM}\) the total heterogeneity of our mixed-effects regression model.

NOTE however, that \(R^2\) refers to variance explained in the observed data. In linear regression, the more predictors you add, the better your model will explain the observed data. But this can decrease the generalizability of you model, through a process called overfitting: You're capturing noise in your dataset, not true effects that exist in the real world. In meta-regression, because the estimate of heterogeneity may differ, we may not actually see the \(R^2\) increasing with more predictors - but because the samples are often small, be very mindful of overfitting nonetheless.

Meta-regression in R

Meta-regressions can be conducted in R using the rma function in metafor. First, let's conduct a regression with a single categorical predictor, donorcode, which has two levels: Anxious and Typical:

m_cat <- rma(yi ~ donorcode, vi = vi, data = df)
m_cat

#> 
#> Mixed-Effects Model (k = 56; tau^2 estimator: REML)
#> 
#> tau^2 (estimated amount of residual heterogeneity):     0.0567 (SE = 0.0177)
#> tau (square root of estimated tau^2 value):             0.2382
#> I^2 (residual heterogeneity / unaccounted variability): 67.81%
#> H^2 (unaccounted variability / sampling variability):   3.11
#> R^2 (amount of heterogeneity accounted for):            0.49%
#> 
#> Test for Residual Heterogeneity:
#> QE(df = 54) = 155.0054, p-val < .0001
#> 
#> Test of Moderators (coefficient 2):
#> QM(df = 1) = 1.5257, p-val = 0.2168
#> 
#> Model Results:
#> 
#>                   estimate      se    zval    pval    ci.lb   ci.ub    
#> intrcpt             0.0829  0.1332  0.6222  0.5338  -0.1781  0.3439    
#> donorcodeTypical    0.1730  0.1401  1.2352  0.2168  -0.1015  0.4476    
#> 
#> ---
#> Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

By default, metafor::rma() will recode factor variables to make one group the reference category and create dummies for all other groups. The syntax above created a binary indicator variable (a dummy variable) for people who fall in the Typical category. It then estimated the model intercept - which is the population effect size in the Anxious category, which serves as reference category here, and the effect of the donorcodeTypical dummy, which tells us the difference between the two groups' population effects.

m_reg <- rma(yi ~sex,
             vi = vi,
             data = df)
m_reg

#> 
#> Mixed-Effects Model (k = 56; tau^2 estimator: REML)
#> 
#> tau^2 (estimated amount of residual heterogeneity):     0.0544 (SE = 0.0173)
#> tau (square root of estimated tau^2 value):             0.2333
#> I^2 (residual heterogeneity / unaccounted variability): 66.50%
#> H^2 (unaccounted variability / sampling variability):   2.98
#> R^2 (amount of heterogeneity accounted for):            4.53%
#> 
#> Test for Residual Heterogeneity:
#> QE(df = 54) = 149.5878, p-val < .0001
#> 
#> Test of Moderators (coefficient 2):
#> QM(df = 1) = 2.1607, p-val = 0.1416
#> 
#> Model Results:
#> 
#>          estimate      se    zval    pval    ci.lb   ci.ub    
#> intrcpt    0.0648  0.1253  0.5168  0.6053  -0.1808  0.3104    
#> sex        0.0050  0.0034  1.4699  0.1416  -0.0017  0.0116    
#> 
#> ---
#> Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

m_multi <- rma(yi ~ sex + donorcode,
               vi = vi,
               data = df)
m_multi

#> 
#> Mixed-Effects Model (k = 56; tau^2 estimator: REML)
#> 
#> tau^2 (estimated amount of residual heterogeneity):     0.0551 (SE = 0.0175)
#> tau (square root of estimated tau^2 value):             0.2347
#> I^2 (residual heterogeneity / unaccounted variability): 66.89%
#> H^2 (unaccounted variability / sampling variability):   3.02
#> R^2 (amount of heterogeneity accounted for):            3.42%
#> 
#> Test for Residual Heterogeneity:
#> QE(df = 53) = 148.6166, p-val < .0001
#> 
#> Test of Moderators (coefficients 2:3):
#> QM(df = 2) = 3.0602, p-val = 0.2165
#> 
#> Model Results:
#> 
#>                   estimate      se     zval    pval    ci.lb   ci.ub    
#> intrcpt            -0.0337  0.1626  -0.2074  0.8357  -0.3523  0.2849    
#> sex                 0.0043  0.0035   1.2310  0.2183  -0.0025  0.0111    
#> donorcodeTypical    0.1361  0.1421   0.9579  0.3381  -0.1424  0.4146    
#> 
#> ---
#> Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

mods <- c("sex", "age", "location", "donorcode", "interventioncode", "controlcode", "recipients", "outcomecode")

res_multi2 <- rma(as.formula(paste("yi ~", paste(mods, collapse = "+"))), vi = vi, data = df)

res_multi2

#> 
#> Mixed-Effects Model (k = 56; tau^2 estimator: REML)
#> 
#> tau^2 (estimated amount of residual heterogeneity):     0.0038 (SE = 0.0076)
#> tau (square root of estimated tau^2 value):             0.0618
#> I^2 (residual heterogeneity / unaccounted variability): 12.55%
#> H^2 (unaccounted variability / sampling variability):   1.14
#> R^2 (amount of heterogeneity accounted for):            93.30%
#> 
#> Test for Residual Heterogeneity:
#> QE(df = 26) = 34.1905, p-val = 0.1303
#> 
#> Test of Moderators (coefficients 2:30):
#> QM(df = 29) = 105.3506, p-val < .0001
#> 
#> Model Results:
#> 
#> 
#> ---
#> Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Optional: Pitfalls of Meta-Regression

Skip this section if you are pressed for time.

#> Warning: Redundant predictors dropped from the model.

df[c(5, 24, 28, 46), c("age", "location", "interventioncode", "controlcode", "outcomecode")]
#>      age location   interventioncode      controlcode outcomecode
#> 5  21.00   Canada Prosocial Spending        Self Help       Other
#> 24 20.00   Canada Prosocial Spending        Self Help   Happiness
#> 28 18.55      USA   Acts of Kindness Neutral Activity       Other
#> 46 29.95      USA   Acts of Kindness Neutral Activity       Other

Overfitting

To better understand the risks of (multiple) meta-regression models, we have to understand the concept of overfitting. Overfitting occurs when we build a statistical model which fits the data too closely. In essence, this means that we build a statistical model which can predict the data at hand very well, but performs bad at predicting future data it has never seen before. This happens if our model assumes that some variation in our data stems from a true "signal" in our data, when in fact we only model random noise (Hastie, Tibshirani, and Friedman 2009). As a result, our statistical model produces false positive results: it sees relationships where there are none.

Illustration of an overfitted model vs. model with a good fit

Regression methods, which usually utilize minimization or maximization procedures such as Ordinary Least Squares or Maximum Likelihood estimation, can be prone to overfitting (Hastie, Tibshirani, and Friedman 2009). Unfortunately, the risk of building a non-robust model, which produces false-positive results, is even higher once we go from conventional regression to meta-regression (J. P. T. Higgins and Thompson 2004). There are several reasons for this:

1. In Meta-Regression, our sample is often small, as we can only use the synthesized data of all analyzed studies \(k\).
2. As our meta-analysis aims to be a comprehensive overview of all evidence, we have no additional data on which we can "test" how well our regression model can predict high or low effect sizes.
3. In meta-regressions, we have to deal with the potential presence of effect size heterogeneity. Imagine a case in which we have two studies with different effect sizes and non-overlapping confidence intervals. Every variable which has different values for the different studies might be a potential explanation for effect size difference we find, while it seems straightforward that most of such explanations are spurious (J. P. T. Higgins and Thompson 2004).
4. Meta-regression as such, and multiple meta-regression in particular, make it very easy to "play around" with predictors. We can test numerous meta-regression models, include many more predictors or remove them in an attempt to explain the heterogeneity in our data. Such an approach is of course tempting, and often found in practice, because we, as meta-analysts, want to find a significant model which explains why effect sizes differ (J. Higgins et al. 2002). However, such behavior massively increases the risk of spurious findings (J. P. T. Higgins and Thompson 2004), because we can change parts of our model indefinitely until we find a significant model, which is then very likely to be overfitted (i.e., it mostly models statistical noise).

Some guidelines have been proposed to avoid an excessive false positive rate when building meta-regression models:

• Minimize the number of investigated predictors a priori. Predictor selection should be based on predefined scientific or theoretical questions we want to answer in our meta-regression.
• When evaluating the fit of a meta-regression model, we prefer models which achieve a good fit with less predictors. Use fit indices that balance fit and parsimony, such as the Akaike and Bayesian information criterion, to determine which model to retain if you compare several models.
• When the number of studies is low (which is very likely to be the case), and we want to compute the significance of a predictor, you can use the Knapp-Hartung adjustment to obtain more reliable estimates (J. Higgins et al. 2002), by specifying test = "knha when calling rma().

Two-level model

First, for simplicity, we run a two-level BRMA model, ignoring the fact that certain effect sizes come from the same study.

df_brma <- df[, c("yi", "vi", mods)]

set.seed(1)
fit_lasso <- brma(yi ~ .,
                  data = df_brma,
                  vi = "vi",
                  method = "lasso",
                  prior = c(df = 1, scale = 1),
                  mute_stan = FALSE)

sum <- summary(fit_lasso)
sum$coefficients[, c("mean", "sd", "2.5%", "97.5%", "n_eff", "Rhat")]

Optional: Three-level model with Horseshoe Prior

The default prior in brma() is the horseshoe prior, because it performed best in simulation studies. Moreover, we can take into account the fact that some effect sizes might come from the same study by fitting a three-level model as follows:

df_brma_threelevel <- df[, c("yi", "vi", "study_id", mods)]
set.seed(7)
fit_hs <- brma(yi ~ .,
               data = df_brma_threelevel,
               vi = "vi",
               study = "study_id")

Check Convergence

For any random forest model, it is important to determine how many trees are needed to get stable results. Convergence is indicated by stabilization of the cumulative mean squared out-of-bag prediction error (MSEoob), as a function of the number of trees in the model. We run the analysis once with a very high number of trees, and pick a smaller number of trees, at which the model is also seen to have converged, to speed up computationally heavy steps, such as replication and model tuning. We re-examine convergence for the final model. Note that we can account for multiple effect sizes per study in MetaForest as well:

# Run model with many trees to check convergence
set.seed(74)
check_conv <- MetaForest(yi~.,
data = df_brma_threelevel,
study = "study_id",
whichweights = "random",
num.trees = 10000)
# Plot convergence trajectory
plot(check_conv)

This model converges at around 2500 trees.

Model Tuning

Machine learning models often have hyperparameters that must be tuned to get good performance for a specific learning task. We can tune the model using the R package caret. As tuning parameters, we consider all three types of weights (uniform, fixed-, and random-effects), the number of candidate variables at each split from 2–6, and a minimum node size from 2–6. We select the model with smallest root mean squared prediction error (RMSE) as the final model, based on 5-fold clustered cross-validation. Clustered cross-validation means that effect sizes from the same study are always included in the same fold, to account for the dependency in the data. Note that the number of folds cannot exceed the number of clusters in the data. Moreover, if the number of clusters is very small, one might have to resort to specifying the same number of folds as clusters. Model tuning is computationally intensive and might take a long time.

# Load caret
library(caret)
# Set up 5-fold clustered CV
grouped_cv <- trainControl(method = "cv",
index = groupKFold(df_brma_threelevel$study_id, k = 5))

# Set up a tuning grid
tuning_grid <- expand.grid(whichweights = c("random", "fixed", "unif"),
mtry = 2:6,
min.node.size = 2:6)
# X should contain only retained moderators, clustering variable, and vi
X <- df_brma_threelevel[, c("study_id", "vi", mods)]
# Train the model
set.seed(3)
mf_cv <- train(y = df_brma_threelevel$yi,
x = X,
study = "study_id", # Name of the clustering variable
method = ModelInfo_mf(),
trControl = grouped_cv,
tuneGrid = tuning_grid,
num.trees = 2500)
# Best model tuning parameters
mf_cv$results[which.min(mf_cv$results$RMSE), ]

#>    whichweights mtry min.node.size      RMSE  Rsquared       MAE     RMSESD
#> 32        fixed    3             3 0.2760904 0.2936028 0.2240218 0.08613236
#>    RsquaredSD      MAESD
#> 32  0.2951682 0.06719638

Based on the root mean squared error, the best combination of tuning parameters were fixed-effect weights, with three candidate variables per split, and a minimum of three cases per terminal node.

Interpreting the Results

The object returned by train already contains the final model, estimated with the best combination of tuning parameters.

# Extract final model
final <- mf_cv$finalModel
# Extract R^2_oob from the final model
r2_oob <- final$forest$r.squared

This model has a positive estimate of explained variance in new data, \(R^2_{oob} = .19\).

Now, we proceed to interpreting the moderator effects, by examining variable importance and partial dependence plots. Variable importance is determined by randomly permuting each variable in turn (so it loses all meaningful association with the outcome), and then assessing how much the model's predictive performance diminishes after this random permutation. For an important variable, the model performance drops after permutation. Partial dependence plots show the model's predicted values at different levels of one variable, averaging over the values of all other variables. Thus, you can think of this as the bivariate marginal association.

# Plot variable importance
VarImpPlot(final)

# Sort the variable names by importance
ordered_vars <- names(final$forest$variable.importance)[
order(final$forest$variable.importance, decreasing = TRUE)]

# Plot partial dependence
PartialDependence(final, vars = ordered_vars,
rawdata = TRUE, pi = .95)

More Information

Reference: Van Lissa (2020)

For a more elaborate tutorial, see this open access book chapter on MetaForest.

Repeating the RMA

Let's repeat our previous random effects meta-analysis, and examine a funnel plot of the distribution of effect sizes:

res_rma <- rma(yi = yi, vi = vi, data = df)
funnel(res_rma)

An effect size over 0 means that helping others increased happiness, even a little bit.

With a Bayes Factor, we could test the informative hypothesis (just a directional null hypothesis) that helping increases happiness in one particular study, for example the first row of the dataset:

set.seed(1)
bain(x = c("yi" = df$yi[1]), # Select the first effect size
     Sigma = matrix(df$vi[1], 1, 1), # Make covariance matrix
     hypothesis = "yi > 0", # Define informative hypothesis
     n = df$n1i[1]+df$n2c[1]) # Calculate total sample size
#> Bayesian informative hypothesis testing for an object of class numeric:
#> 
#>    Fit   Com   BF.u  BF.c     PMPa  PMPb  PMPc 
#> H1 0.999 0.500 1.999 1508.434 1.000 0.667 0.999
#> Hu                                  0.333      
#> Hc 0.001 0.500 0.001                      0.001
#> 
#> Hypotheses:
#>   H1: yi>0
#> 
#> Note: BF.u denotes the Bayes factor of the hypothesis at hand versus the unconstrained hypothesis Hu. BF.c denotes the Bayes factor of the hypothesis at hand versus its complement. PMPa contains the posterior model probabilities of the hypotheses specified. PMPb adds Hu, the unconstrained hypothesis. PMPc adds Hc, the complement of the union of the hypotheses specified.

The results indicate that the first study provides 1508 times more evidence in favor of our informative hypothesis (helping increases happiness) than against it, \(BF = 1508.434\).

Product Bayes Factor

We now perform the PBF analysis, using the pbf() method for numeric input (see ?pbf). This interface is very similar to rma(), and is specifically designed for applications where PBF is applied to meta-analytic datasets. The yi and vi arguments are the same as those of rma(). Additional argument ni is used to construct the prior for the approximate Bayes factors (Hoijtink et al. 2019). Importantly, the hypothesis argument determines which informative hypotheses are tested.

set.seed(1)
res_pbf <- pbf(yi = df$yi,
               vi = df$vi,
               ni = sum(df$n1i+df$n2c),
               hypothesis = "y > 0")

The resulting output shows the product Bayes factors for our informative hypothesis (PBF column), as well as study-specific Bayes factors (remaining columns, trimmed here to keep the document tidy). The PBF in favor of our hypothesis is very high (a 7.45 with 69 zeroes), so the body of evidence provides consistent support for the notion that helping increases happiness.

This last question relates to generalizability. When we perform a normal fixed- or random-effects meta-analysis, we rely on sampling theory and assume that each effect size represents a random sample from the population of possible effect sizes.

If this assumption is justified, then we can claim that the pooled effect size generalized to the population those effect sizes were randomly sampled from. Thus, our pooled effect size tells us something about future studies.

Of course, this assumption is very questionable, due to publication bias, and the fact that researchers build upon previous work! So in practice, meta-analysis might not be as generalizable as we'd like.

The PBF does not rely on sampling theory and does not make such assumptions about random sampling from the population. As a result, it can only be interpreted as a quantitative summary of the evidence within the present sample of studies.

It would be incorrect to interpret the PBF as saying anything about the population; we should interpret it as a summary statistic of the studies at hand. Of course, the same could be said about conventional meta-analysis, if the assumption of random sampling is violated.

More Information

References: Van Lissa, Kuiper, and Clapper (2023), Van Lissa et al. (2020)

For a more elaborate tutorial on the Product Bayes Factor, see this paper, in press at Research Synthesis Methods.

Credit

This tutorial is partly based on the book Doing Meta-Analysis with R: A Hands-On Guide, which is also recommended as further reading:

Harrer, M., Cuijpers, P., Furukawa, T.A., & Ebert, D.D. (2021). Doing Meta-Analysis with R: A Hands-On Guide. Boca Raton, FL and London: Chapmann & Hall/CRC Press. ISBN 978-0-367-61007-4.