Understanding Meta-Analysis Through Data Simulation With Applications to Power Analysis

Abstract

Meta-analysis is a powerful tool to combine evidence from existing literature. Despite several introductory and advanced materials about organizing, conducting, and reporting a meta-analysis, to our knowledge, there are no introductive materials about simulating the most common meta-analysis models. Data simulation is essential for developing and validating new statistical models and procedures. Furthermore, data simulation is a powerful educational tool for understanding a statistical method. In this tutorial, we show how to simulate equal-effects, random-effects, and metaregression models and illustrate how to estimate statistical power. Simulations for multilevel and multivariate models are available in the Supplemental Material available online. All materials associated with this article can be accessed on OSF (https://osf.io/54djn/).

If you do not simulate it, you have not understood it.

A meta-analysis is an essential tool for combining knowledge from multiple studies quantitatively. Meta-analysis is commonly used together with a systematic review of the literature. The meta-analysis has several advantages. First, it allows combining evidence from multiple studies, assigning more weight to studies with lower estimation variability. Then using metaregression, it is possible to include variables (i.e., moderators) to explain the observed heterogeneity (Borenstein et al., 2009, pp. 187–203). More recently, location-scale models have been developed to include predictors also on the residual heterogeneity (i.e., Viechtbauer & López-López, 2022). Finally, given the replication crisis, there are statistical methods to determine the presence and extent of the publication bias. Despite the advantages, meta-analysis implementation is not always straightforward, especially for complex data structures. In addition, even though there are several introductory and advanced resources to understand meta-analysis (Borenstein et al., 2009; Harrer et al., 2021; Schmid et al., 2022), to our knowledge, there are no introductory resources about how to simulate realistic meta-analytic data.

Simulating data has several advantages because it requires understanding the statistical method and the data-generation process. Furthermore, data simulation is the primary tool when it comes to evaluating a new analysis method, estimating the statistical power, or understanding the long-run behavior of one’s data-generation process (Gelman et al., 2020, pp. 69–76; Gelman & Hill, 2006, pp. 155–176; Ingalls, 2011). A recent article by DeBruine and Barr (2021), which deeply inspired the current work, proposed a stimulating way to understand linear mixed-effects models via data simulation. Simulating data is also a powerful educational tool within this framework.

For these reasons, in this work, we aim to introduce the basic concepts of meta-analysis and Monte Carlo simulations for equal-effects, random-effects, and metaregression models with applications also to statistical power calculation. In the first section, we introduce basic concepts of the meta-analysis that are useful for setting up the simulation. We evaluate the effect size, variance calculation, and the equal- versus random-effects-model distinction. Then we describe how to simulate data for these models and simulate a metaregression with categorical and numerical predictors. Finally, we introduce the power analysis extending the previous examples to estimate the statistical power. We used the R statistical programming language (Version 4.3.1; R Core Team, 2023).

The aim of the tutorial is not to provide a complete theoretical introduction to meta-analysis but, rather, to present core topics using a simulation-based approach. Readers experienced in conducting meta-analyses can benefit from the proposed approach in terms of simulation setup and coding strategies. Readers without prior experience in meta-analysis can benefit from both the theoretical introduction and the simulation approach. However, for a more comprehensive overview of meta-analysis topics, the reader may refer to meta-analysis textbooks (e.g., Borenstein et al., 2009; Harrer et al., 2019).

We assume the reader is familiar with the basic concepts of R, but core functions will be explained. Code and materials are available on the OSF repository (https://osf.io/54djn/). For a theoretical introduction, simulation examples for multivariate and multilevel models, and more details about the coding approach, see the Supplemental Material available online.

Meta-Analysis Introduction

The meta-analysis is a statistical procedure to combine multiple studies (i.e., “primary studies”) into a single statistical analysis (Borenstein et al., 2009). The idea is that combining numerous preliminary studies improves the estimation of a particular phenomenon more efficiently compared with conducting a single study. In statistical terms, the concept of the meta-analysis is to switch the statistical unit from the single participant or observation (i.e., “Level 1”) to the study (i.e., “Level 2”). Given that some studies give more information because their estimation variability is smaller (e.g., higher sample size), the meta-analysis combines the studies, assigning more weight as a function of the precision (i.e., the inverse of the variance).

As an example that will be used throughout the article, we consider the efficacy of memory training in improving memory performance during a cognitive task. The typical primary study will collect data from a group of participants receiving the memory training (“experimental group”) and another group receiving a control treatment (“control group”). The focus of the meta-analysis is collecting multiple studies with similar aims and methods and estimating the average effect of memory training. Despite differences in the type of cognitive task or experimental setup, each primary study collects an experimental group ( $n_T$ ) and control group ( $n_C$ ) and computes the average performance (experimental group $T¯$ and control group $\overline{C}$ ) and standard deviations (experimental group $s_T$ and control group $s_C$ ).

Effect Size and Variance

The first step of a meta-analysis is to extract information from included studies. This common measure should give an immediate idea of the direction (i.e., the treatment improves or reduces performance) and the size of the effect. Standardized effect-size measures (for an overview, see Lakens, 2013), such as the standardized mean difference, which could be estimated by Cohen’s $d$ (Cohen, 1988), or the Pearson correlation coefficient $\rho$ , which could be calculated using the associated sample estimator r, are commonly used to compare heterogeneous outcome variables. ¹ If all studies used the same raw measure, such as reaction times, it is possible to directly meta-analyze the studies without standardizing, for example, using the unstandardized mean difference (UMD; Borenstein et al., 2009, pp. 21–24).

As reported in the previous section, beyond the effect size of each study, one needs to assign a weight according to the precision. For this reason, one needs to calculate the sampling variability of the effect size or the raw measure that will represent the estimation precision. Each raw or standardized effect-size measure (e.g., raw mean difference, Cohen’s $d$ , and Pearson’s correlation) has a different formula to calculate the sampling variability. The idea is to choose the appropriate measure considering the study design (e.g., between- vs. within-subjects) and available information and find the proper formula to compute the sampling variability. In addition, there are formulas and approaches to convert from one effect-size measure to another (Borenstein et al., 2009; Lakens, 2013; Lipsey & Wilson, 2001). ² Usually, the effect-size sampling variability depends mainly on the sample size that determines the weight assigned during the meta-analytic estimation.

Equal- Versus Fixed- Versus Random-Effects Model

The core of a meta-analysis is combining the results of multiple studies, giving more weight to studies that provide a more precise effect estimation. Essentially, there are three meta-analysis models: the equal-effects, the fixed-effects, and the random-effects model (Hedges & Vevea, 1998; Laird & Mosteller, 1990). The equal-effects model assumes that each study included in the meta-analysis is a more or less precise estimation of the true underlying effect (θ). In other terms, one is assuming that there is no variability (i.e., heterogeneity) among true effect sizes. On the other side, suppose some study-level characteristics (e.g., participants’ age, sex, or socioeconomic status) or the experimental paradigm (e.g., type of memory task or difficulty) could affect the treatment effect. In this case, there is true variability (i.e., heterogeneity) among studies. The random-effects model assumes a distribution of real effects with mean μ_θ and variance ${\tau }^2$ , thus estimating the heterogeneity. Finally, the fixed-effects model estimates the average effect of the included pool of studies ignoring the presence of heterogeneity. Although the distinction between the equal- and fixed-effects model is theoretical, the estimated model is the exactly the same. ³ The random-effects model assumes and estimates heterogeneity among effect sizes, leading to different model parameters. From an inferential point of view, the random-effects model provides unconditional inference to the population of effect sizes, whereas the fixed-effects model estimates the average effect of the selected studies and not population-level parameters. The equal-effects model assumes a single true population effect to be estimated. Figure 1 depicts the theoretical distinction between the equal- and random-effects models. In the presence of heterogeneity, moderators (e.g., the type of experimental setup) can be included into a metaregression model to explain the effect-size heterogeneity. An interesting proposal to combine evidence from these different models into a single analysis is called Bayesian model-averaged meta-analysis (Berkhout et al., 2023; Gronau et al., 2021). However, this model is beyond the scope of the tutorial.

Fig. 1.

The difference between the assumptions of the equal- and random-effects models. Each distribution depicts the sampling distribution of $k=4$ hypothetical studies ( $i=1$ , 2, . . ., 4) with a certain observed effect size $y_i$ (pink squares), sampling variability ${\sigma }_{{ϵ}_i}^2$ , and the 95% confidence interval (black segment). The equal-effects plot on the left suggests that each observed effect size has the same underlying true effect θ (each distribution has the same mean) with a different degree of precision (e.g., Study 1 is more precise than Study 3). In practice, each study has a different observed effect size such that studies with high precision (i.e., narrow sampling distributions) will be close to the real effect $(\theta )$ . The random-effects model on the right suggests that beyond the error term $({ϵ}_i)$ , each real effect size is composed of a fixed part (now μ_θ) and a random part $({\delta }_i)$ sampled from a normal distribution with mean zero and variance ${\tau }^2$ . When ${\tau }^2$ is zero, the random-effects model reduces to a equal-effects model.

Simulation

Monte Carlo simulations

The Monte Carlo methods are controlled experiments (Gentle, 2009). Given a set of fixed parameters, probability distributions, and the possibility of generating random numbers, it is possible to simulate the behavior of an empirical system. Monte Carlo simulations are used for statistical and mathematical problems that cannot be solved analytically.

A straightforward example regards estimating the sampling variability of the mean difference. When calculating the mean difference between two samples, one estimates the true mean difference at the population level with a certain degree of error (i.e., the standard error of the mean difference). The central-limit theorem states that the difference between the means of two random samples ( $\overline{X}$ and $\overline{Y}$ ) is approximately normally distributed with mean ${\mu }_x-{\mu }_y$ and standard error $\sqrt{\frac{{\sigma }_x^2}{n_x}+\frac{{\sigma }_y^2}{n_y}}$ . The same results can be obtained using Monte Carlo simulations using the following procedure:

1. Generate two random samples from two normal distributions with a fixed mean difference.

2. Calculate the mean difference.

3. Repeat the same process many times.

4. Calculate the standard deviation of the simulated values.

Using the method, one is estimating via simulation the standard deviation of the sampling distribution of the mean difference (i.e., the standard error of the mean difference). Increasing the number of simulations will produce more stable results:

The standard error estimated solving analytically is 0.26, and using the Monte Carlo simulation, one arrives at the same result (i.e., 0.26).

Simulation setup

This tutorial uses several R packages for the simulations, meta-analysis fitting, and figures/tables. For the data manipulation, we used the tidyverse (Wickham, 2023) package. For the models fitting, we used the metafor (Viechtbauer, 2010) package. For figures and tables, we used the ggplot2 (Wickham et al., 2023) or the metafor::forest() function, kableExtra (Zhu, 2021), and papaja (Aust & Barth, 2022) packages. We set the seed for reproducibility of the simulation environment.

Before diving into the specific simulations, in this section, we define the common aspects of all simulations in the following sections. In the current article, we focus on the two-level, equal-effects and random-effects models. All the examples refer to primary studies that assess the efficacy of a treatment by comparing a control and an experimental group. In simulation studies in which the purpose is not evaluating effect-size estimators, it is convenient to simulate unstandardized effect-size measures (see Viechtbauer, 2005, 2007). The estimator is unbiased, thus not requiring small-sample correction (Hedges, 1981, 1989). Furthermore, the effect size and the sampling variance are independent. Similar to the simulation approach by Viechtbauer (2005, 2007), the experimental group ( $T$ ) and control groups ( $C$ ) are sampled from normal distributions, respectively, $T_i~\mathcal{N}(\Delta ,1)$ and $C_i~\mathcal{N}(0,1)$ , where $\Delta$ is the UMD. The UMD and the sampling variance are calculated using Equations 1 and 2, where $y$ is the estimated value of $\Delta$ :

$$D=\overline{T}-\overline{C}.$$

(1)

$$\begin{array}{l}{\sigma }_{ϵD}^2=\frac{s_T^2}{n_T}+\frac{s_C^2}{n_C}\hfill \end{array}.$$

(2)

We can use the following algorithm implemented in the sim_study() function to simulate a single study. Before using the sim_study() function, we can create a data frame for the simulation using the make_data() function. ⁴ Table 1 depicts an example of the make_data() output:

Table 1.

Example of Data Generated With the make_data() function

id	nt	nc	es
1	20	20	0.3
2	20	20	0.3
. . .	. . .	. . .	. . .
29	20	20	0.3
30	20	20	0.3

Note: The id column is the identifier for each study. nt = number of participants in the experimental and control group; nc = number of participants in the experimental and control group; es = effect size.

1. Choose a $\Delta$ , $n_T$ , and $n_C$ value.

2. Simulate $n_T$ observations from a Gaussian distribution with $\mu =\Delta$ and ${\sigma }^2=1$ and $n_T$ observations from a Gaussian distribution with $\mu =0$ and ${\sigma }^2=1$ . In this way, the expected difference between groups will be $\Delta$ , and the expected variance for each group is 1.

3. Calculate the observed effect size $y$ and the sampling variance ${\sigma }_{ϵ}^2$ .

The simulation approach can be easily extended to calculating a standardized effect-size measure (e.g., Cohens’ $d$ ; Cohen, 1988) and the corresponding sampling variance. For example, after generating data for the two groups, the mean difference can be standardized using the pooled standard deviation and applying the appropriate correction (Hedges, 1981; e.g., Hedges, 1989):

The es is the true effect size (μ $({\mu }_{\theta }$ for the random-effects model and $\theta$ for the equal-effects model), nc and nt are the sample size for the control and experimental groups, aggregate controls if returning the effect size and the corresponding sampling variance or the participant-level data. We can generate a single study with the desired parameters with this function. A suggestion for each simulation step is to generate a large $n$ to reduce the sampling error and check the recovery of simulated parameters. This is a general strategy that can be applied to every simulation. The sim_studies() function will iterate through variables in the . . . argument, creating the meta-analysis data frame. The mapply() function is clearly explained in the Supplemental Material.

The following code simulates a single study with $n=10,000$ and checks the estimated mean and standard deviation:

The control group has a mean of −0.008 (SD = 0.996), and the experimental group has a mean of 0.306 (SD = 1.013), which are remarkably close to the simulated values.

Using the sim_study() function multiple times (with the appropriate adjustments), we can generate a series of studies simulating a data set for a meta-analysis (using the sim_studies() function). After each example, we will compute the appropriate model (e.g., equal or random effects) using the metafor package (Viechtbauer, 2010) to check the recovery of simulated parameters. ⁵ Table 2 summarizes the notation used in simulations and code in equations and code:

Table 2.

Variables in the Data-Generating Model and Associated R Code

Equation	Code	Description
$\mathrm{\theta }$	theta	Equal-effects model true effect size
${\mathrm{\mu }}_{\mathrm{\theta }}$ μ	mu_theta	Random-effects model true effect size
${\mathrm{\tau }}^2$	tau2	Effect sizes’ heterogeneity
${\mathrm{\tau }}_r^2$	tau2r	Residual effect sizes’ heterogeneity
$\overline{T},\overline{C}$		Mean of the experimental/control group
$s_{T,C}$		Standard deviation of the experimental/control group
$n_{T,C}$	nt, nc	Sample size of the experimental/control group
$y_i$	yi	Observed effect size
${\sigma }_{{ϵ}_i}^2$	vi	Observed sampling variance
${\mathrm{\delta }}_i$	deltai	Random effect for the study i
${\mathrm{\beta }}_0$	b0	Metaregression intercept
${\mathrm{\beta }}_1$	b1	Metaregression slope
$k$	k	Number of studies
${ϵ}_i$		Sampling error for the study i

Note: The yi and vi notation for the observed effect size $(y_i)$ and sampling variance $({\sigma }_{{ϵ}_i}^2)$ has been used to be consistent with the metafor notation.

Equal-effects model

The most basic model to simulate is the equal-effects model. As reported in the previous sections, the equal-effects model assumes the presence of a single true effect ( $\theta$ ), and the observed variability is caused by each effect size being a more or less imprecise estimation of the true effect. In other words, the only source of variability is the sampling variability that depends on the variance of the primary studies (i.e., the sample size). Equations 3 and 4 formalize the equal-effects model:

$$\begin{array}{l}{y}_i=\mathrm{\theta }+{ϵ}_i\hfill \end{array}.$$

(3)

$$\begin{array}{l}{ϵ}_i~\mathcal{N}(0,{\sigma }_{{ϵ}_i}^2)\hfill \end{array}.$$

(4)

Each observed effect size $(y_i)$ is composed of the real effect size $\theta$ plus an error term $({\epsilon }_i)$ that is sampled from a normal distribution with $\mu =0$ and ${\sigma }^2={\sigma }_{{ϵ}_i}^2$ (i.e., the known sampling variance of the study $i$ ). As demonstrated in Equation 1, the increase in sample size will decrease the sampling variability. A study with an extremely large sample size will essentially have $\theta$ as the observed effect size. Because we are sampling participants’ level data, the error component is already included in the sim_study() function. To simulate this model in R, we can just call the sim_study() multiple times according to the number of desired studies ( $k$ ). In addition, we simulate that each primary study will have a sample size of $n=20$ for both groups. We discuss later the appropriateness of this assumption:

We are using the $\theta$ parameter to generate random data using the sim_study() function that will introduce the random error component ${ϵ}_{\mathrm{i}}$ . Then we can fit the equal-effects meta-analysis model with the rma function and method = "EE" of the metafor package. The parameter we are estimating is $\theta$ , which is close to our simulation value (see Table 3). Increasing the number of studies $(k)$ and/or the number of participants in each study ( $n$ ) will improve estimation, reducing the standard error.

Table 3.

Summary of the Simulated Equal-Effects Model

	β	95% CI	z	p
Overall	0.339 (SE = 0.057)	$[0.228,0.451]$	5.959	< .001
$k=30$

Note: The only estimated parameter is the average effect $\theta$ $(\beta )$ with the standard error, 95% confidence interval, and the Wald z test. The z test evaluates the null hypothesis that the real effect equals zero. CI = confidence interval.

Increasing the sample size for each study will increase the estimation precision; thus, the variability among studies will be reduced. This can be easily demonstrated by simulating studies with a high sample size, as reported in Figure 2. As the sample size increases, the only source of variability (i.e., the error component ${ϵ}_i$ ) is close to zero, so each study is closer to the true simulated value.

Fig. 2.

Forest plots of two simulated equal-effects models. On the left, the simulated model has $nt,nc=$ $30$ for each included study, and on the right, the sample size for each study is $nt,nc=500$ . Given that the data were generated under an equal-effects model, when the sample size is high (on the right), each study is aligned on the real effect size because the error component $({ϵ}_i)$ is close to zero. The average effect size is similar (depending on the random-numbers generation) between the two scenarios, whereas the estimation precision (the width of the black diamond) is narrower on the right.

Random-effects model

The random-effects model can be considered an extension of the equal-effects model. The equal-effects model assumes that the real effect is a single value. The random-effects model relaxes this assumption, allowing the true effect size to vary across studies. For example, the difference between groups we are simulating could be influenced by the type of experiment or the participants’ age. Now, $\theta$ is no longer a single value but a distribution of values. Because of the effect size being a distribution, we need to estimate both the mean and the variance. The parameter ${\mu }_{\theta }$ is mean of the distribution, interpreted as the average effect size across different true effect sizes, and ${\tau }^2$ is the variance of the distribution, interpreted as variability or heterogeneity of effect sizes. In practical terms, we now have two sources of variability: ${\tau }^2$ , which expresses the real difference among effect sizes, and ${\sigma }_{{ϵ}_i}^2$ , which is the known sampling variance of each study, as in the equal-effects model. We can easily extend Equation 3 with Equations 5 through 7:

$$\begin{array}{l}{y}_i={\mathrm{\mu }}_{\theta }+{\mathrm{\delta }}_i+{\mathrm{ϵ}}_i\hfill \end{array}.$$

(5)

$$\begin{array}{l}{\delta }_i~\mathcal{N}(0,{\tau }^2)\hfill \end{array}.$$

(6)

$$\begin{array}{l}{ϵ}_i~\mathcal{N}\left(0,{\sigma }_{{ϵ}_i}^2\right)\hfill \end{array}.$$

(7)

Compared with the equal-effects model, we need to generate another adjustment to the overall effect ${\mu }_{\theta }$ from a normal distribution with mean 0 and variance ${\tau }^2$ . The real effect size for the study $i$ will be ${\mu }_{\theta }+{\delta }_i$ , where ${\delta }_i$ are the random-effects regulated by the ${\tau }^2$ parameter:

Then we can fit the random-effects meta-analysis model with the rma function and method = "REML". Table 4 depicts the model results ⁶ of the metafor package. We are now estimating two parameters: ${\mu }_{\theta }$ and ${\tau }^2$ . As for the equal-effects model, increasing the number of studies and/or the number of participants in each study will improve estimation, reducing the standard error.

Table 4.

Summary of the Simulated Random-Effects Model

	$\beta$	95% CI	z	p
Overall	0.387 (SE = 0.107)	$[0.177,0.597]$	3.616	< .001
$k=30$
${\tau }^2=0.246$ ( $SE=0.090$ )
$I^2=71.925\%$

Note: Compared with the equal-effects model, there are more parameters. The $\beta$ is the average effect ( $({\mu }_{\theta })$ ) with the standard error, 95% confidence interval, and the Wald z test. The ${\tau }^2$ is the estimated heterogeneity, and the $I^2$ (explained in the Random-Effects Model section) represents the percentage of total variability due to between-studies heterogeneity.

An important aspect of the random-effects model is the interplay between the heterogeneity $({\tau }^2)$ and the sampling variability $({\sigma }_i^2)$ . As the number of studies $k$ increases, the estimation of ${\mu }_{\theta }$ and ${\tau }^2$ becomes more precise (Blázquez-Rincón et al., 2023; Rubio-Aparicio et al., 2018). In addition, as the sample size of each study decreases, each ${\delta }_i$ will be estimated with higher precision, but as long as ${\tau }^2\ne 0$ , there will be variability among effect sizes. In other words, increasing the sample size of each study or the number of studies will not affect the value of ${\tau }^2$ but only the estimation precision (for a clear explanation, see Borenstein et al., 2009, Chapter 16, Figure 16.6). This can be easily demonstrated using the previous simulation, increasing each study sample size. Figure 3 depicts the same meta-analysis but with different precision in estimating the study $y_i$ . Compared with Figure 2, increasing the sample size of primary studies improves the estimation of each study without reducing the between-studies heterogeneity.

Fig. 3.

Forest plots of two simulated random-effects models. On the left, the simulated model has $nt,nc=30$ for each included study, and on the right, the sample size for each study is $nt,nc=$ $500$ . Compared with Figure 2, data were generated under a random-effects model. The estimated average effect is similar between the two scenarios regarding average effect and precision. However, compared with the equal-effects simulation, increasing the sample size of primary studies affects only the precision without reducing the true heterogeneity (i.e., ${\tau }^2\ne 0$ ).

The relationship between the sampling error and the heterogeneity can be expressed using the $I^2$ statistics (e.g., Higgins & Thompson, 2002), which is the percentage of the total variability ${\tau }^2+\tilde{v}$ that is attributable to real heterogeneity between studies $({\tau }^2)$ . In Equation 8, $\tilde{v}$ is the typical within-studies sampling variability (i.e., a summary statistics of sampling variances of included studies) as proposed by Higgins and Thompson (2002) ⁷ and implemented in Equation 9, with $w_i=1/{\sigma }_{{ϵ}_i}^2$ , and $k$ is the number of studies:

$$\begin{array}{l}{I}^2=\frac{{\tau }^2}{{\tau }^2+\tilde{v}}\hfill \end{array}.$$

(8)

$$\begin{array}{l}\tilde{v}=\frac{(k-1){{\displaystyle \sum }}^{\text{}}{w}_i}{{({{\displaystyle \sum }}^{\text{}}{w}_i)}^2-{{\displaystyle \sum }}^{\text{}}{w}_i^2}\hfill \end{array}.$$

(9)

From Equation 8 and Figure 3, it is clear that if each included study has a considerable sample size, the sampling variability $(\tilde{v})$ will be reduced, and the total variability will be mainly driven by real heterogeneity $({\tau }^2)$ . This is the crucial difference between the equal-effects and the random-effects models (see also Borenstein et al., 2009, pp. 117–122).

Given the interpretation of $I^2$ , it is possible to simulate a meta-analysis fixing a certain $I^2$ value. The only caveat is fixing the $\tilde{v}$ . When the sample size of each study is the same, ⁸ $\tilde{v}={\sigma }_{{ϵ}_i}^2$ . In the other case, $\tilde{v}$ needs to be calculated from sampling variances as reported in Higgins and Thompson (2002) and cannot be easily fixed a priori. With the assumption of homogeneous sample size across studies, ⁹ we can solve Equation 8 for ${\tau }^2$ , obtaining the heterogeneity value associated with a certain $I^2$ as reported in Equation 10.

Table 5 depicts the results of the random-effects model fixing the $I^2$ value. Choosing a meaningful ${\tau }^2$ for the simulation can be sometimes difficult. Beyond fixing the $I^2$ value, a possibility is choosing plausible ${\tau }^2$ values from the literature. For example, van Erp and colleagues (2017) estimated an empirical ${\tau }^2$ distribution across several published meta-analyses that can be used to simulate a plausible scenario:

Table 5.

Summary of the Random-Effects Model Fixing the $I^2$ Value

	$\beta$	95% CI	z	p
Overall	0.305 (SE = 0.072)	$[0.164,0.447]$	4.239	< .001
$k=30$
${\tau }^2=0.090$ ( $SE=0.041$ )
$I^2=58.183\%$

Note: The β is the average effect (m_θ)with the standard error, 95% confidence interval, and the Wald z test. The t2 is the estimated heterogeneity, and the I ² (explained in the Random-Effects Model section) represents the percentage of total variability due to between-studies heterogeneity. CI = confidence interval.

$$\begin{array}{l}{\tau }^2=-\frac{I^2\tilde{v}}{I^2-1}\hfill \end{array}.$$

(10)

Metaregression

From a linear-regression perspective, both the equal-effects and the random-effects models can be seen as intercept-only models in which only the mean (i.e., the linear regression intercept) is estimated. As reported in the meta-analysis introduction, the between-studies heterogeneity usually represents the true variability of the effect due to differences among primary studies. A natural extension of the intercept-only analysis is a model that includes variables (i.e., moderators) that could explain the observed heterogeneity among effect sizes. For example, a group of studies could use a particular memory task in which the expected effect is higher than another. In this way, considering the type of task will explain part of the observed heterogeneity, as in standard regression models. Figures 4 and 5 depict a random-effects metaregression model for a categorical and numerical predictor.

Fig. 4.

Graphical representation of a random-effects metaregression model with a categorical predictor (Condition A and Condition B). Each gray distribution represents the sampling distribution of included studies. The dotted line is the average effect (i.e., random-effects model without moderators). The effect size differs between Conditions A and B, including the condition moderator, explaining part of the total heterogeneity (pink plus green segments). The green segments depict the explained heterogeneity, and the pink segments depict the residual (unexplained) heterogeneity. The pink squares are simulated observed effect sizes from the sampling distributions.

Fig. 5.

Graphical representation of a random-effects metaregression model with a numerical predictor ( $x$ ). Each gray distribution represents the sampling distribution of included studies. The dotted line is the average effect (i.e., random-effects model without moderators). The effect size increases as a function of the $x$ variable. Therefore, including $x$ as a predictor explains the heterogeneity. The green segments depict the explained heterogeneity, and the pink segments depict the residual (unexplained) heterogeneity. The pink squares are simulated observed effect sizes from the sampling distributions.

Metaregression with a categorical moderator

A common example of metaregression is by including a categorical predictor, including information about study-level features. In our example, a group of studies uses an online memory task, whereas others use a standard lab-based task. Equation 5 can be easily extended for a metaregression model by including a variable encoding the type of task (online vs. lab-based) and the expected difference between the two levels of the moderator (i.e., the lab vs. online effect). In regression terms (see Equation 11), we could use a dummy variable ( $X_1$ ) that takes the value of 0 for the lab-based task (L) and a value of 1 for an online task (O). ¹⁰ Now, we fix ${\beta }_1$ to be the lab versus online effect (i.e., the expected mean difference between the two groups of studies) and the product between ${\beta }_1{X}_{1_i}$ will consider the lab versus online effect. Crucially, even though ${\tau }^2$ is still the heterogeneity between effect sizes, now we need to fix ${\tau }^2$ , considering that we included a moderator. In other terms, ${\delta }_i~\mathcal{N}(0,{\tau }_r^2)$ , where ${\tau }_r^2$ is the residual heterogeneity after including the moderator. We can describe our model using Equation 12 according to the value of $X_{1_i}$ :

$$\begin{array}{l}{y}_i={\beta }_0+{\delta }_i+{\beta }_1{X}_{1_i}+{ϵ}_i\hfill \end{array}.$$

(11)

$$\begin{array}{l}{y}_{L_i}={\beta }_0+{\delta }_i+{\beta }_1\times 0+{ϵ}_i\hfill \\ y_{O_i}={\beta }_0+{\delta }_i+{\beta }_1\times 1+{ϵ}_i\hfill \end{array}.$$

(12)

We can simulate the same scenario of the random-effects model with $k_O=15$ (online tasks) and $k_L=15$ (lab-based tasks). Then we fix the ${\beta }_1=0.2$ and ${\tau }_r^2=0.1$ ( $r$ for residual):

Now we can fit the metaregression model with the rma function, as for the random-effects model, with the addition of mods = ~ exp, which indicates which variables to consider as moderators. The results are presented in Table 6. Now the model will estimate an intercept parameter $(\mathrm{i}.\mathrm{e}.,{\beta }_0)$ that is the value of $y$ when $X_1$ is zero (i.e., for lab-based studies), or in other words, the expected value for lab-based studies. Then the ${\beta }_1$ parameter represents the estimated difference in $y$ between the values of $X_1$ (i.e., lab-based vs. online experiments). As said before, ${\tau }_r^2$ is now the residual heterogeneity that is interpreted as the variability between effect sizes after controlling for the moderator $X_1$ .

Table 6.

Summary of the Random-Effects Model With a Categorical Predictor (Lab vs. Online Experiments)

	$\beta$	95% CI	z	p
Intercept	0.043 (SE = 0.103)	$[-0.158,0.244]$	0.42	.674
Exponline	0.323 (SE = 0.144)	$[0.040,0.606]$	2.24	.025
$k=30$
${\tau }_r^2=0.091$ ( $SE=0.042$ )
$R^2=19.824\%$
$I^2=58.323\%$

Note: The intercept is the average effect for lab-based experiments, and exponline is the difference between lab-based and online experiments. The $R^2$ is the percentage of explained heterogeneity, and ${\tau }_r^2$ is the Estimated residual heterogeneity. The b is the average effect (m_θ)with the standard error, 95% confidence interval, and the Wald z test. The I ² (explained in the Random-Effects Model section) represents the percentage of total variability due to between-studies heterogeneity. CI = confidence interval.

Metaregression with a numerical moderator

The same approach can be used for a continuous predictor. For example, we can simulate that the average participant’s age within each study could explain part of the observed heterogeneity. Now, $X_1$ is a continuous predictor representing the average age for each study, and ${\beta }_1$ is the effect-size increase for a unit increase (i.e., 1 year) of average age. Sometimes guessing a plausible ${\beta }_1$ value with a continuous predictor is not straightforward. A first strategy could be to use values estimated from the literature. Another approach consists of setting up the model and simulating several expected $y_i$ and calculating the range of simulated values. A third possibility is fixing the proportion of explained heterogeneity and calculating the ${\beta }_1$ value accordingly. As in standard regression analysis, we can use the $R^2$ statistic to describe the amount of heterogeneity explained by the included moderators. Equation 13 reports how to calculate the $R^2$ for a metaregression model. The ${\tau }_r^2$ is the residual heterogeneity after considering the moderators, and ${\tau }_f^2$ is the heterogeneity estimated without considering the moderators. In the next sections, we present an example of the simulation-based and the $R^2$ -based approaches:

$$\begin{array}{l}{R}^2=1-\frac{{\tau }_r^2}{{\tau }_f^2}\hfill \end{array}.$$

(13)

In terms of regression parameters, now the intercept $({\beta }_0)$ is no longer the overall effect or the average of one category, as in the previous example, but the estimated value for a specific $X_1$ , thus for a specific age. If $X_1$ is a variable representing the average age for each study, then the ${\beta }_0$ is the average effect size when the age is zero. Depending on the moderator, the intercept is interpreted in different ways. For example, with the age, the intercept has no empirical meaning given that no studies could have a participant average age of zero. A strategy could be to mean-center the age (i.e., subtracting from each study age the average age across studies). Now the intercept is still the average effect size when the age is zero, but now zero is the average age. Note that the contrast coding for categorical predictors or centering numerical variables does not affect the overall model but only parameters’ values and interpretation.

Assessing the impact of ${\beta }_1$

As reported in the previous section, a strategy to guess plausible values for ${\beta }_1$ is by simulating several expected $y_i$ given the metaregression equation and summarizing or plotting the effect-size range. The range of simulated $y_i$ values are also affected by the simulated age values across studies. However, it is probably more intuitive to guess a plausible range of moderator values compared with the ${\beta }_1$ value. In this specific example, if all studies target a specific population (e.g., adults below 50 years), the expected average age range can be easily simulated. In our case, we simulated $k$ average age values from a uniform distribution ${\mathrm{age}}_i~U\left(20,40\right)$ . Then we can plot the distribution of $y_i$ values to check the plausibility of simulated values. As shown in Figure 6, with the same range for the moderator, a ${\beta }_1=0.1$ gives a plausible effect sizes range, whereas a ${\beta }_1=0.7$ predicts very extreme values:

Fig. 6.

Scatter plots with marginal histograms for the range of simulated $y_i$ with two ${\beta }_1$ values. The x-axis depicts the average age of simulated studies, and the y-axis depicts the simulated effect size. On the left, the majority of simulated values range between −1.5 and 1.5; thus, ${\beta }_1=0.1$ can be considered a plausible value. On the right, $(\beta =0.7)$ values range between −10 and 10, which, despite being theoretically possible, can be considered highly implausible values in real meta-analyses for psychological data.

Simulating using ${\mathrm{R}}^2.$

A more intuitive way to simulate a continuous predictor is fixing the desired $R^2$ value and finding the coefficient that produces the desired value. This approach has been implemented by López-López and colleagues (2014). We can use Equations 14 and 15 to find the ${\beta }_1$ value that is associated with a certain $R^2$ :

$$\begin{array}{l}{\mathrm{\beta }}_1^2={\mathrm{\tau }}^2{R}^2\hfill \end{array}.$$

(14)

$$\begin{array}{l}{\mathrm{\tau }}_r^2={\mathrm{\tau }}^2-{\mathrm{\beta }}_1^2\hfill \end{array}.$$

(15)

Now, we can simulate the regression model using $\sqrt{{\beta }_1^2}$ as coefficient and ${\tau }_r^2$ as residual heterogeneity. Results from the fitted model fixing the $R^2$ values are presented in Table 7 and Figure 7. As López-López and colleagues (2014) demonstrated, to reliably estimate $R^2$ , the number of studies needs to be large. ¹¹ López-López and colleagues generated the moderator $(X_1)$ values from a standard normal distribution. In the following example, we standardized the moderator (scale()) after simulating values on the age scale (e.g., runif(k, 20, 40)):

Table 7.

Summary of the Random-effects Model Fixing the $R^2$ Value

	$\beta$	95% CI	z	p
Intercept	0.340 (SE = 0.051)	$[0.240,0.439]$	6.680	< .001
Age0	0.168 (SE = 0.051)	$[0.068,0.269]$	3.296	< .001
$k=100$
${\tau }_r^2=0.194$ ( $SE=0.037$ )
$R^2=11.422\%$
$I^2=75.417\%$

Note: The intercept is the effect size for the average age (given that age is mean-centered). The age0 parameter is the slope between age and the effect size, interpreted as an increase in effect size for a unit increase in the average age. Both parameters report the standard error, 95% confidence interval, and the Wald z test. The R² is the percentage of explained heterogeneity and τ²_r is the estimated residual heterogeneity. I ² represents the percentage of total variability due to between-studies heterogeneity. CI = confidence interval.

Fig. 7.

Metaregression results for the random-effects model with a numerical moderator. Each effect size is represented with a black dot where the dimension represents the weight according to the inverse of the variance. The line represents the estimated metaregression slope with the 95% confidence interval (gray bands).

Power Analysis

The previous simulation examples can be easily implemented for multiple purposes. For example, we can use different effect sizes and variance estimators when using the sim_study() function to check the impact on the fitted meta-analysis model. However, one of the most critical applications is estimating the power of a specific statistical model. One of the purposes of power analysis by simulations is to estimate the required number of studies to detect a hypothetical effect size when planning a meta-analysis. At the same time, the meta-analysis that we presented (two-level equal- or random-effects model) can be considered as a multilab (e.g., Klein et al., 2018) study in which experiments are planned and not collected from the literature. We can use the same simulation approach to optimize the number of participants and studies when planning a multilab project.

As explained in the introduction, there are several approaches and tools to estimate the power of equal- and random-effects models (see Borenstein et al., 2009, Chapter 29; Harrer et al., 2019, Chapter 14). These methods are easy to implement but made strong assumptions, such as the homogeneity of sample size, and did not consider the uncertainty in estimating ${\tau }^2$ . Jackson and Turner (2017) partially solved the issue by developing an interesting method that takes into account the uncertainty in estimating ${\tau }^2$ without using simulations. However, complex simulation scenarios can be handled only using Monte Carlo methods.

A general Monte Carlo simulation for the power analysis can be implemented with the following steps:

1. Choose the model that generates the data (e.g., equal- or random-effects model).

2. Fix the relevant parameters (e.g., ${\tau }^2$ and $\theta$ ).

3. Simulate a data set.

4. Fit the appropriate model.

5. Store the p value associated with the parameter of interest.

6. Repeat Steps 3 to 5 a large number of times (e.g., 10,000).

7. Calculate the power as the proportion of p values below the $\alpha$ level.

For example, we can estimate the power of a random-effects model by repeating the simulation presented in the Random-Effects Model section many times. We simulated heterogeneity of sample sizes by sampling $n_T$ and $n_C$ values from a Poisson distribution with $\lambda =20$ . In this way, on average, the sample size is 20 for primary studies with a certain amount of heterogeneity. Usually, it is more informative to simulate different scenarios according to the relevant parameters, such as sample sizes, number of studies, or heterogeneity. For example, we can estimate the power with a different number of studies k. We define the do_sim() function that, according to the input parameter, repeats the simulation a certain number of times (i.e., nsim). ¹² Increasing the number of simulations will increase the power analysis estimation precision. Then the summary_sim() function analyzes each simulation, returning the relevant values. We repeat the simulation of the random-effects model several times with different parameters.

Simulation results are presented in Figure 8 and Table 8 showing that to reach 80% power (usually considered an appropriate level) with a = 0.05 we need ~35 studies. The same approach could be used to estimate the power of a meta-regression by simply modifying the do_sim() function simulating the effect of a moderator and extracting the relevant p-value.

Fig. 8.

Results from the random-effects model power analysis. The x-axis depicts the number of studies ( $k$ ), and the y-axis depicts the estimated power. The pink dotted line is the 80% power level, usually considered a good value for power analysis.

Table 8.

The Results From the Power Analysis Simulation

k	$\Delta$	${\tau }^2$	$n_{\mathrm{avg}}$	$n_{\min }$	n sim	Power
5	0.30	0.30	20	10	5,000	0.26
15	0.30	0.30	20	10	5,000	0.47
25	0.30	0.30	20	10	5,000	0.65
35	0.30	0.30	20	10	5,000	0.79
50	0.30	0.30	20	10	5,000	0.92

Note: The table depicts the simulation parameters, the estimated power using the summary_sim() function, and the average sample size ( $n$ ) across the simulations.

Conclusions

In the present work, we introduced the basic concepts of the meta-analysis regarding equal-effects, random-effects, and metaregression models with a simulation-based approach. We believe the presented examples are useful to implement alternative or more complex models. For example, the sim_study() function can be easily modified to simulate another effect-sizes index, such as correlations or odds ratios. In addition, more complex models, such as multivariate or multilevel models, can be simulated following a similar approach (see the Supplemental Material). The multilevel (e.g., three-level) model estimates another heterogeneity component representing the variability of multiple independent effect sizes within the same study. Likewise, the multivariate model includes the correlation between multiple outcomes and the correlation between sampling errors.

The present work did have a few limitations. First, we introduced only basic concepts about meta-analysis and Monte Carlo simulations, whereas setting up complex simulations requires more knowledge and complexity of the simulation setup. We decided to give the foundations to understand meta-analyses with a simulation approach because more complex models are still based on the same principles. Second, there are limitations concerning simulating participant-level data. We decided to simulate the meta-analysis data starting from the participant level to maximize the flexibility and clearness of each step. The downside concerns the efficiency and scalability of the simulation setup. For large-scale simulations (e.g., many conditions, iterations, or complex models), simulating from aggregated statistics is probably more efficient (for an example, see Heuvel et al., 2020) to improve the simulation efficency. ¹³

In conclusion, data simulation is a very powerful tool for each step of a data-analysis process, starting from the learning phase, in which simulating data can be used to understand the statistical model in terms of assumptions and the data-generation process, to the estimation of statistical power. Moreover, we believe that data simulation as part of a standard research workflow could improve overall research quality. Data simulation requires understanding the statistical model, setting appropriate and reasoned parameters, and realizing how the chosen analysis method behaves across different scenarios.

Supplemental Material

sj-pdf-1-amp-10.1177_25152459231209330 – Supplemental material for Understanding Meta-Analysis Through Data Simulation With Applications to Power Analysis

Supplemental material, sj-pdf-1-amp-10.1177_25152459231209330 for Understanding Meta-Analysis Through Data Simulation With Applications to Power Analysis by Filippo Gambarota and Gianmarco Altoè in Advances in Methods and Practices in Psychological Science