Equivalencies Between Ad Hoc Strategies and Multivariate Models for Meta-Analysis of Dependent Effect Sizes

Abstract

Meta-analyses of educational research findings frequently involve statistically dependent effect size estimates. Meta-analysts have often addressed dependence issues using ad hoc approaches that involve modifying the data to conform to the assumptions of models for independent effect size estimates, such as by aggregating estimates to obtain one summary estimate per study, conducting separate analyses of distinct subgroups of estimates, or combinations thereof. We show that these ad hoc approaches correspond exactly to certain multivariate models for dependent effect sizes. Specifically, we describe classes of multivariate random effects models that have likelihoods equivalent to those of models for effect sizes that have been averaged by study, classified into subgroups, or both. The equivalencies also apply to robust variance estimation methods.

Keywords

meta-analysis dependent effect sizes aggregated effects shifting unit of analysis

1. Dependent Effect Sizes

Across fields ranging from education and psychology to clinical medicine and health sciences, it is very common that some of the primary studies included in a meta-analysis report more than one effect size estimate based on the same sample (Ahn et al., 2012; Page et al., 2015; Tipton et al., 2019). For example, in a meta-analysis of effects of an educational intervention on reading outcomes, some studies might have collected multiple measures of reading performance (each of which meets the inclusion criteria of the synthesis) or may have measured outcomes at multiple follow-up times. Some studies might have tested more than one version of an intervention, such that it is of interest to include effect sizes comparing each version to a no-intervention control condition. Studies might even have several of these features, potentially contributing a large number of relevant effect size estimates.

When a meta-analysis includes multiple effect size estimates from the same sample, the estimates will typically be dependent, or correlated, because they are based on data from the same set of individuals. The meta-analysis will therefore need to use analytic methods that account for the dependence structure. In a review of methods for handling dependent effect size estimates, López-López and colleagues (2018) delineated two types of approaches: “reductionist” strategies, which modify the data structure to conform to the assumptions of a model appropriate for independent effect size estimates, and “integrative” strategies, which involve developing models for the full data, including dependent effect sizes. In this article, we shall examine the connections between these two classes of approaches.

1.1. Reductionist Strategies

Meta-analysts have often handled dependent effect size estimates by making ad hoc modifications to the data in order to meet the assumptions of a simple model for independent effect size estimates. One way to do so is to average the dependent effect size estimates from each study, yielding a dataset with a single summary effect size estimate from each independent sample (Borenstein et al., 2009; Rosenthal & Rubin, 1986). Another pragmatic approach is to classify effect sizes from each sample into subgroups (e.g., based on the type of outcome measure), where each subgroup includes at most one effect size from each sample (Becker, 2000). For situations where it is not feasible to distinctly classify effect sizes from every sample, Cooper (1998, 2016) described a “shifting unit of analysis” strategy. This entails creating distinct subgroups to the extent possible, then averaging together any dependent effect size estimates within each subgroup. The analyst then either conducts separate analyses of each subgroup or estimates a single model to compare results across subgroups.

These reductionist strategies originated out of practical necessity, before the development of software tools that could be used to model dependent effects (Becker, 2000). Even after such tools became widely available, aggregating, subgrouping, and shifting the unit of analysis have remained in regular use. In a review of meta-analytic practices based on syntheses published in 2016, Tipton and colleagues (2019) found that averaging effect sizes and shifting the unit of analysis were used frequently across several fields. For example, Furenes et al. (2021) aggregated effect sizes to the study level in a meta-analysis examining difference in children’s learning outcomes when using digital versus paper books. Other recent meta-analytic reviews on topics within education and psychology have applied the shifting unit of analysis approach (Graham et al., 2020; Jung et al., 2020; Koenka et al., 2021). Thus, meta-analysts have continued to use reductionist strategies in reviews published within the last 5 years.

1.2. Integrative Strategies

Integrative strategies entail developing analytic models for the full, unmodified data, which necessitates using methods that accommodate dependent observations. Multivariate meta-analysis models for dependent effect size estimates were described early in the history of meta-analysis methodology (Hedges & Olkin, 1985, Chapter 10; Raudenbush et al., 1988) and further developed during the 1990s (Becker, 1992; Berkey et al., 1998; Kalaian & Raudenbush, 1996; van Houwelingen et al., 2002). These models required specifying the degree of correlation between effect size estimates drawn from the same study. Although formulas and methods for estimating such correlations are available (Gleser & Olkin, 2009; Wei & Higgins, 2013), their use requires knowing the correlation between outcomes in each sample—information that is infrequently reported in primary studies. This difficulty, combined with the lack of software for estimating multivariate models, long limited the use of multivariate meta-analytic models (Becker, 2000).

More recent methodological developments in meta-analysis have brought renewed interest in multivariate models that can accommodate dependent effect sizes (Jackson et al., 2011). Innovations such as robust variance estimation (RVE; Hedges et al., 2010), pseudo-likelihood methods (Chen et al., 2014; Hong et al., 2018), multilevel meta-analytic models (Konstantopoulos, 2011; Van den Noortgate et al., 2013, 2015), and combinations thereof (Fernández-Castilla et al., 2020; Pustejovsky & Tipton, 2022) provide tools for directly modeling the effect size estimates, without requiring ad hoc modifications to the data structure or statistical information that is seldom reported in primary studies. As a result, these more recently developed approaches can help researchers to develop richer and more nuanced insights from meta-analytic reviews (Cheung, 2019). Although these more recently developed integrative methods have potential advantages over reductionist approaches, they are also more complex and less accessible to researchers without specialized training.

1.3. Aims

In light of the continued use of ad hoc approaches and the greater complexity of more recent innovations, meta-analysts might well have cause to wonder: Are newer methods like RVE or multilevel meta-analysis always necessary? Do they make obsolete the older approaches for handling-dependent effects? Or are simpler approaches such as averaging effect sizes, subgrouping, or shifting the unit of analysis still sometimes useful and reasonable to apply?

In this article, we aim to shed some partial light on these questions by demonstrating equivalencies between older, reductionist approaches that involve aggregating or subgrouping of effect size estimates and newer, multivariate modeling approaches for dependent effect sizes. In Section 2, we describe conditions under which estimates from a random effects meta-regression of averaged effect size estimates correspond precisely to those from a certain multivariate model called the correlated effects model. In Section 3, we describe conditions under which the shifting unit of analysis approach is equivalent to the subgroup correlated effects model for dependent effect sizes, as described by Pustejovsky and Tipton (2022). These equivalence relations concern the likelihoods of the models, implying that full and restricted maximum likelihood estimates from the reductionist models correspond exactly to estimates from the integrative models. The equivalencies also apply to RVE methods for each model, as we detail in Section 4.

The equivalence relations provide a basis for guidance about when the simpler, reductionist methods may be validly applied in practice. The representations of the simpler approaches as multivariate models also provide a way to directly compare them to more complex models, which is useful both for further theoretical development and for practical application, such as in likelihood ratio tests for comparison of nested models. We discuss further implications in Section 5.

1.4. General Notation

Certain matrix operations will be essential in the subsequent development. We will use $[\cdot]$ to denote vector or matrix stacking, such that $x = {[x_{i}]}_{i = 1}^{n}$ is a column vector with scalar entries $x_{1}, \dots, x_{n}$ . Similarly for a set of $n_{j} \times m$ matrices $X_{1}, \dots, X_{J}$ , $X = {[X_{j}]}_{j = 1}^{J}$ is the $(\sum_{j = 1}^{J} n_{j}) \times m$ matrix that stacks the component matrices vertically. We will use $\oplus$ to denote the Kronecker sum of scalars or matrices, such that $\oplus_{i = 1}^{n} x_{i} = diag (x_{1}, \dots, x_{n})$ and $\oplus_{j = 1}^{J} X_{j} = diag (X_{1}, \dots, X_{J})$ , and $\otimes$ to denote the Kronecker product.

We shall consider a collection of J studies, where study j contributes $k_{j} \geq 1$ effect size estimates. Let $T_{i j}$ denote effect size estimate i in study j, with sampling variance $s_{i j}^{2}$ . Effect size estimates from study j have correlated sampling errors, with correlation $ρ_{h i j}$ between effect size estimates h and i from study j. Throughout, we assume that the meta-analyst knows the true correlations, although in practice they might instead make a rough estimate of the degree of correlation, such as by assuming $ρ_{h i j} = 0.7$ for all pairs of estimates from each included study. Let $T_{j} = {[T_{i j}]}_{i = 1}^{k_{j}}$ be the $k_{j} \times 1$ vector of effect size estimates from study j. Let $S_{j}$ be the $k_{j} \times k_{j}$ sampling covariance matrix for $T_{j}$ , given by

S_{j} = [\begin{matrix} s_{1 j}^{2} & ρ_{12 j} s_{1 j} s_{2 j} & \dots & ρ_{1 k_{j} j} s_{1 j} s_{k_{j} j} \\ ρ_{21 j} s_{2 j} s_{1 j} & s_{2 j}^{2} & \dots & ρ_{2 k_{j} j} s_{2 j} s_{k_{j} j} \\ ⋮ & ⋮ & ⋱ & ⋮ \\ ρ_{k_{j} 1 j} s_{k_{j} j} s_{1 j} & ρ_{k_{j} 2 j} s_{k_{j} j} s_{2 j} & \dots & s_{k_{j} j}^{2} \end{matrix}] .

In the following sections, we will make reference to the likelihood and restricted likelihood based on a model for normally distributed effect size estimates. With a multivariate normal random vector $T$ of length K, let $X$ be a $K \times p$ matrix of predictors, let $β$ be a vector of p regression coefficients; and let $Ψ (θ)$ be the marginal variance matrix of $T$ under a given model, which is a function of a parameter vector $θ$ . Consider the regression model

T \sim N (X β, Ψ (θ)) .

For this model, the log of the likelihood is

ln L (β, θ) = - \frac{1}{2} ln | Ψ (θ) | - \frac{1}{2} (T - X β) Ψ^{- 1} (θ) (T - X β) - K c,

where $c = \frac{1}{2} ln (2 π)$ (Searle et al., 2009).

Because maximum likelihood estimators of variance components are known to be biased, researchers often prefer to use variance estimators that maximize the restricted likelihood. The restricted likelihood can be understood as the likelihood of a linear combination of $T$ that does not depend on $β$ , or as the marginal likelihood for $θ$ after integrating over the regression coefficients $β$ (Searle et al., 2009). The log of the restricted likelihood is

\begin{array}{l} ln R (θ) = - \frac{1}{2} ln | X' Ψ^{- 1} (θ) X | - \frac{1}{2} ln | Ψ (θ) | \\ - \frac{1}{2} (T - X \hat{β} (θ))^' Ψ^{- 1} (θ) (T - X \hat{β} (θ)) - (K - p) c, \end{array}

where $\hat{β} (θ) = {(X' Ψ^{- 1} (θ) X)}^{- 1} X' Ψ^{- 1} (θ) T$ is the generalized least squares estimator for $β$ , written as a function of $θ$ (Searle et al., 2009).

2. Aggregating Effect Sizes

One reductionist method to handle the dependence between effect size estimates is to average—or aggregate—the effects within each study and then develop a study-level meta-regression model. Let ${\bar{T}}_{j}$ denote a weighted average of the effects from study j, and let V_j denote the sampling variance of this averaged effect size estimate. Given the sampling variances and correlations between effect size estimates from study j, the weights that minimize the sampling variance of the averaged effect size are based on the inverse of the sampling variance matrix $S_{j}$ . Using inverse-variance weights for aggregation,

{\bar{T}}_{j} = \frac{1_{j}^{'} S_{j}^{- 1} T_{j}}{1_{j}^{'} S_{j}^{- 1} 1_{j}} and V_{j} = \frac{1}{1_{j}^{'} S_{j}^{- 1} 1_{j}},

where $1_{j}$ is a $k_{j} \times 1$ vector of 1s.

In many applications, the sampling variances for individual effect sizes within a given study j will all be equal (or nearly equal), so that one can assume $s_{h j} = s_{i j} = s_{j}$ . One might also assume that there is a constant sampling correlation between every pair of effect sizes within a given study, so that $ρ_{h i j} = ρ_{j}$ . Under these assumptions, $S_{j}$ has a compound symmetric structure and the inverse-variance weighted average of the effect size estimates is identical to the simple average:

{\bar{T}}_{j} = \frac{1}{k_{j}} \sum_{i = 1}^{k_{j}} T_{i j},

with sampling variance

V_{j} = \frac{(k_{j} - 1) ρ_{j} + 1}{k_{j}} \times s_{j}^{2} .

The approach to aggregating effect sizes based on Equations 5 and 6 has been described in introductory meta-analysis texts (cf. Borenstein et al., 2009, Eq. 24.6) and implemented in software (e.g., the aggregate() function in the metafor package for R, Viechtbauer, 2010). Aggregating has been used in several recent reviews (e.g., De Wit et al., 2021; Emslander & Scherer, 2022; Gaskell et al., 2023). For instance, De Wit et al. (2021) conducted a meta-analysis to examine the effects of procedural learning on individuals with Alzheimer’s disease or amnestic mild cognitive impairment. In order to account for the dependence arising from multiple outcomes within some studies, the authors aggregated the effect size estimates within each study, assuming a common correlation between each pair of effects within each study.

2.1. Meta-Regression Model

Once the meta-analyst has calculated study-level aggregated effect size estimates, they would then proceed to estimate summary meta-analysis or meta-regression models. Let $x_{j}$ be a row vector of predictor variables for study j. These predictor variables might encode characteristics, such as publication year, study design type (e.g., randomized experiment versus quasi-experiment), use of specific intervention components, or any other relevant features of included studies. A conventional random effects meta-regression model for the averaged effect sizes is

{\bar{T}}_{j} = x_{j} β + u_{j} + e_{j},

where u_j is a random effect for study j, assumed to be normally distributed with mean zero and variance $τ^{2}$ ; and $Var (e_{j}) = V_{j}$ , the sampling variance of the aggregated effect size estimate from study j. This model can be estimated using standard software routines (e.g., rma.uni() from metafor). Note that the predictors in this model can only represent study-level characteristics because the model pertains to the study-level aggregated effects.

Rather than averaging effect size estimates, the meta-analyst might estimate a multivariate meta-regression model for the raw effect size estimates. Consider the model

T_{i j} = x_{j} β + u_{j} + e_{i j},

where u_j is a normally distributed, between-study random effect with mean zero and variance $τ^{2}$ and where $e_{i j}$ is a normally distributed sampling error for effect size i from study j. The sampling errors from study j are assumed to be multivariate-normally distributed with mean zero and known variance matrix $S_{j}$ , so that $Var (e_{i j}) = s_{i j}^{2}$ and $cor (e_{h j}, e_{i j}) = ρ_{h i j}$ . Following the terminology of Pustejovsky and Tipton (2022), this specification is a “correlated effects” model. It can be estimated using the rma.mv() function in the metafor package (Viechtbauer, 2010).

If the aggregated effect size estimates are calculated using inverse-variance weights as given in Equation 4, then the likelihood of the random effects meta-regression model in Equation 7 is equivalent to the likelihood of the correlated effects model given in Equation 8. By equivalent, we mean that the likelihoods are proportional. The equivalence of these models follows from Theorem 2 in Section 3.2. Equivalence implies that full maximum likelihood estimators and restricted maximum likelihood estimators based on Model (7) are, by definition, identical to the corresponding estimators based on Model (8). Model-based standard errors based on Model (7) are also equal to the standard errors based on (8). Furthermore, in a Bayesian inference framework, the posterior distributions of both models will be identical so long as the same priors for $β$ and for $τ$ are used in both models.

Because of its equivalence with the correlated effects model, the strategy of aggregating effect size estimates by study can be understood not only as a reductionist strategy but also as a multivariate model for the raw effect size estimates, which can be compared to other multivariate models. Conversely, the correlated effects model can be understood as a model based on averaging effect size estimates to the study level. This suggests that one could use averaged effect sizes for the purposes of creating model diagnostics or graphical representations (forest plots, funnel plots, etc.) of meta-regression results based on the correlated effects model.

For equivalence between these models to hold exactly, the averaged effect sizes must be calculated using inverse-variance weights. If the analyst uses other weights or other methods for approximating the sampling variance of ${\bar{T}}_{j}$ , then the equivalence will not necessarily be exact, although in practice, estimates based on the two models will often be quite similar. In order for the multivariate model for the raw effect sizes to correspond exactly to the simpler model for the averaged effect sizes, it is necessary that Model (8) is limited to study-level predictors and study-level random effects. If the specification includes predictors or random effects terms that vary within study, then the multivariate model will not necessarily correspond to any model for the study-level averaged effect size estimates.

The equivalence between Models (7) and (8) covers estimators based on the full or restricted likelihood but does not necessarily extend to other classes of estimators. Hedges et al. (2010) proposed a method of estimating the correlated effects model using a moment estimator of $τ$ and an “approximately efficient” weighted least squares estimator of $β$ . These methods are implemented in the robumeta package for R (Fisher et al., 2017) and Stata (Hedberg, 2014). However, this estimation strategy does not correspond exactly to any conventional approach for estimating Model (7) with the aggregated effect size estimates due to the form of the moment estimator and to the use of weights that are not exactly inverse variance under the correlated effects model. We give an example of the differences in Section 2.3.

2.2. More Complex Meta-Regression Models

When inverse-variance weights are used for averaging effect sizes, the equivalence between models extends to specifications that are more complex than the random effects meta-regression from Model (7). More generally, consider a model for the averaged effect sizes that has additional random effects terms:

{\bar{T}}_{j} = x_{j} β + z_{j} u + e_{j},

where $z_{j}$ is a $1 \times q$ row vector describing the random effects structure for study j, $u$ is a $q \times 1$ vector of random effects, and e_j is the normally distributed sampling error of the averaged effect size estimate from study j, with mean zero and $Var (e_{j}) = V_{j}$ . The random effects vector $u$ is normally distributed with mean zero and $q \times q$ variance matrix $Φ (θ)$ , which is a function of an $r \times 1$ parameter $θ$ (and possibly also a further set of covariates for each study). We provide a detailed example of the structure of the $z_{j}$ vector and random effects variance matrix $Φ (θ)$ in Online Appendix A.

This more general model encompasses a variety of different random effects structures that are of interest in different contexts (cf. Sera et al., 2019). For example, the J studies included in a synthesis might have been conducted by a smaller number of investigators or lab groups, each of which tends to use similar operational procedures across studies (e.g., similar experimental manipulations, preferred outcome measurement instruments, etc.). In such a situation, the meta-analyst might use a hierarchical random effects model (Konstantopoulos, 2011) to partition the between-study heterogeneity into between-lab variation and within-lab variation. In this specification, the random effects vector $u$ would include lab-level random effects and study-level random effects, with the vector $z_{j}$ consisting of sets of indicators identifying the unique labs and unique studies. Other examples of random effects structures that fall within Model (9) include the random effects location-scale model (Viechtbauer & López-López, 2022), phylogenetic meta-analytic models used in the synthesis of ecological data across multiple species (Hadfield & Nakagawa, 2010; Lajeunesse, 2009; Nakagawa & Santos, 2012), and cross-classified random effects meta-analysis models (Fernández-Castilla et al., 2019).

Rather than averaging effect sizes before fitting Model (9), the meta-analyst might instead fit a multivariate model for the raw effect size estimates. By Theorem 2 in Section 3.2, Model (9) is equivalent to the model

T_{i j} = x_{j} β + z_{j} u + e_{i j},

where the predictors $x_{j}$ and the random effects structure $z_{j}$ vary at the level of the study but are constant across effect size estimates within a given study, and where the sampling errors for study j, $e_{j} = {[e_{i j}]}_{i = 1}^{k_{j}}$ , are normally distributed with mean zero and known variance $S_{j}$ . Thus, meta-regression models based on aggregated effect sizes are equivalent to correlated effects models for the raw effect size estimates under quite general conditions, including models that have complex random effects structures at the study level.

2.3. An Illustration

We present an empirical example to illustrate the equivalence between random effects meta-regression models for aggregated effect size estimates and multivariate meta-regression models for raw effect size estimates, as well as to highlight some differences in auxiliary statistical results typically reported with univariate and multivariate models. Clinton-Lisell (2021) conducted a meta-analysis of studies investigating the effects of multitasking on reading performance and reading time. For the purposes of illustration, we reanalyzed the effects for reading performance outcomes, which were drawn from 20 primary studies. Effect size estimates were expressed as standardized mean differences with positive values representing improvement in reading performance. The dataset includes 31 effect size estimates because some studies report results for more than one measure of reading performance, resulting in dependent effect sizes. We estimated a model that includes two study-level predictors (out of five moderators examined in the original paper): pacing control (self-paced or experimenter-paced) and reading medium (screen or paper).

We used the metafor R package (Viechtbauer, 2010) for all analyses. We used the aggregate function to calculate the study-level average effect sizes and sampling variances. Following the original study, we assumed a common sampling correlation of $0.8$ between each pair of effect sizes within a study. We then fit a random effects meta-regression model for the aggregated effects using the rma.uni function. We fit the corresponding correlated effects model for the raw effect sizes using the rma.mv function, with the degrees of freedom calculated using the “contain” method. For both models, we report profile likelihood confidence intervals for the between-study heterogeneity $τ$ . Data and code for replicating all calculations are available at https://osf.io/9u52e/.

The top panel of Table 1 reports the estimates of the random effects meta-regression for aggregated effects (first column) and the multivariate correlated effects model (second column). It can be seen that the point estimates, standard errors, and confidence intervals for the fixed and random effects are identical across both models.

Table 1.

Comparison of Aggregated Meta-Regression Model and Multivariate Meta-Regression Model

Coefficient	Aggregated Effects Model			Multivariate Correlated Effects Model
Coefficient	Est	SE	95% CI	Est	SE	95% CI
Random effects meta-regression
Intercept	−0.550	.148	[−0.862, −0.238]	−0.550	.148	[−0.862, −0.238]
Self-paced	0.355	.121	[0.100, 0.610]	0.355	.121	[0.100, 0.610]
Screen reading	0.102	.136	[−0.185, 0.390]	0.102	.136	[−0.185, 0.390]
τ	0.177		[0.106, 0.311]	0.177		[0.106, 0.311]
Log likelihood	−1.801			−2,892.491
Q	392.644			6,207.101
Location-scale meta-regression
Intercept	−0.549	.135	[−0.834, −0.263]	−0.549	.135	[−0.834, −0.263]
Self-paced	0.358	.099	[0.150, 0.566]	0.358	.099	[0.150, 0.566]
Screen reading	0.090	.134	[−0.192, 0.372]	0.090	.134	[−0.192, 0.372]
τ (self-paced)	0.000		[0.000, 0.284]	0.000		[0.000, 0.284]
τ (exp-paced)	0.204		[0.115, 0.398]	0.204		[0.115, 0.398]
Log likelihood	−0.771			−2,891.461
Q	392.644			6,207.101

To illustrate that this equivalence extends to more complex models, we also estimated specifications where the between-study heterogeneity is allowed to vary as a function of pacing control (self-paced vs. experimenter-paced), as in a random effects location-scale model (Viechtbauer & López-López, 2022). The bottom panel of Table 1 reports the estimates of the location-scale model estimated using both the aggregated effect sizes and the multivariate equivalent. As expected, the regression coefficients and variance component estimates are all identical across both approaches, as are the associated standard errors and confidence intervals.

Although the parameter estimates are identical regardless of whether the models are estimated based on the aggregated or raw effect sizes, the log likelihoods and Q statistics differ depending on the specification. Theorem 2 shows that the difference in log likelihoods is a constant that does not depend on the model parameters. As a result, the difference in log likelihoods is the same in both the random effects specification (top panel) and the location-scale specification (bottom panel). A likelihood ratio test comparing the models with and without the scale predictor supports the simpler model in the top panel ( $χ_{1}^{2} = 2.10$ , $p = .151$ ).

The difference in the Q statistics arises because the Q statistic for the aggregated random effects model captures only between-study heterogeneity, whereas the Q statistic for the correlated effects model captures both between-study and within-study heterogeneity in effect size estimates. The difference in Q statistics can therefore be interpreted as a measure of within-study heterogeneity. In this example, the large difference in Q statistics indicates that there is substantial excess heterogeneity in effect sizes drawn from the same sample ( $p < .001$ compared to a $χ^{2}$ distribution with $K - J = 11$ degrees of freedom). This suggests that aggregating by study is concealing variation in the effect sizes.

Another method of evaluating within-study heterogeneity is to fit a multivariate model that includes additional random effects for each estimate, as in the correlated-and-hierarchical-effects (CHE) model described by Pustejovsky and Tipton (2022). A likelihood ratio test comparing the CHE model to the correlated effects model supports the more complex CHE model ( $χ_{1}^{2} = 5, 773.60$ , $p < .001$ ), consistent with conclusions based on the difference in Q statistics.

The original analysis by Clinton-Lisell (2021) used a correlated effects model estimated using the robu function from the robumeta package (Fisher et al., 2017). Using robu to fit the same specification as in the top panel of Table 1 gives coefficient estimates of −0.587 (SE = .118, 95% CI = [−0.863, −0.312]) for the intercept, 0.389 (SE = .118, 95% CI = [0.131, 0.646]) for self-paced, and .081 (SE = .110, 95% CI = [−0.167, 0.329]) for screen reading medium, with a between-study heterogeneity estimate of $τ = 0.295$ . These estimates differ somewhat from the results reported in Table 1 because robu uses different estimation methods that are not likelihood-based.

3. Subgrouping and Shifting the Unit of Analysis

Meta-regression models based on aggregated effect size estimates are constrained to using study-level predictors, yet meta-analysts may also want to consider predictors that vary within studies. The meta-analyst might use subgrouping or shifting the unit of analysis to handle dependence when analyzing effect size estimates classified by a categorical predictor. For example, Koenka et al. (2021) reported a meta-analytic review of the effects of different forms of feedback on academic motivation. The authors examined a variety of moderators, several of which varied across effect sizes within a study. For example, many primary studies reported effect sizes for multiple measures of motivation that differ in focus (i.e., internal vs. external focus). In order to investigate whether motivation focus explains variation in effects, the authors applied the shifting unit of analysis technique, aggregating across effect sizes from the same study that have the same focus, and then using meta-regression to examine differences in effects depending on whether the outcomes assessed internal or external focus.

Consider a categorical predictor with G categories and let $k_{g j}$ denote the number of effects in category g from study j. If $k_{g j}$ is never more than one, then the set of effect size estimates in category g might be treated as mutually independent. If multiple effects from a given study fall into the same category g (so $k_{g j} > 1$ ), the meta-analyst could use the shifting unit of analysis strategy, taking a weighted average of effects within each subgroup to remove dependence within the subgroup (although not between subgroups). In either scenario, the meta-analyst could then use a basic random effects meta-regression to analyze the effect size estimates from category g, applying the same approach for each category of interest. Alternately, the meta-analyst might choose to ignore the dependence between effect size estimates from the same study that fall in different categories, fitting a single random effects meta-regression model to make the comparisons of results across categories.

These strategies can be formally described using a set of indicator variables for each category of the predictor. Let $c_{i j}^{g}$ be an indicator equal to 1 if effect size estimate $T_{i j}$ falls into category g and otherwise equal to zero, for $i = 1, ..., k_{j}$ , $j = 1, ..., J$ , and $g = 1, ..., G$ . Define the $k_{j} \times 1$ vector $c_{j}^{g} = {[c_{i j}^{g}]}_{i = 1}^{k_{j}}$ . Let $K_{g j}$ denote the indices of effect sizes in study j that are included in group g; let $G_{j}$ denote the set of category indices for which study j includes one or more effect size estimates, for $j = 1, ..., J$ ; and let $J_{g}$ denote the set of size J_g consisting of the studies that include one or more effect size estimates from category g, for $g = 1, ..., G$ . Define the $k_{g j} \times k_{j}$ matrix $C_{j}^{g} = {[({[I (h = i) \times c_{i j}^{g}]}_{h = 1}^{k_{j}})^']}_{i \in K_{g j}}$ , which reduces the k_j rows in study j to those included in group g, for $g \in G_{j}$ and $j = 1, ..., J$ . Let $p_{g}$ be a $1 \times G$ row vector indicating category g, so that entry g is equal to 1 and all remaining entries are zero, for $g = 1, ..., G$ . Let $P_{j} = {[p_{g}]}_{g \in G_{j}}$ be the set of indicators for subgroups appearing in study j, for $j = 1, ..., J$ .

Following the same reasoning as in the previous section, the most efficient method of calculating weighted average effect size estimates is based on inverse-variance weighting of the effects within each category. With inverse-variance weighting, the weighted average effect size estimate for category g of study j is calculated as

{\tilde{T}}_{g j} = \frac{c_{j}^{g'} S_{j}^{- 1} C_{j}^{g'} C_{j}^{g} T_{j}}{c_{j}^{g'} S_{j}^{- 1} c_{j}^{g}},

with sampling variance given by

{\tilde{V}}_{g j} = \frac{1}{c_{j}^{g'} S_{j}^{- 1} c_{j}^{g}} .

In the common situation where the meta-analyst assumes that effect size estimates from study j all have the same sampling variance $s_{j}^{2}$ and same degree of correlation $ρ_{j}$ , then the inverse-variance weighted average is the same as the simple average of the estimates within a given category:

{\tilde{T}}_{g j} = \frac{1}{k_{g j}} \sum_{i = 1}^{k_{j}} c_{i j}^{g} T_{i j},

with sampling variance

{\tilde{V}}_{g j} = \frac{(k_{g j} - 1) ρ_{j} + 1}{k_{g j}} \times s_{j}^{2} .

If there is only a single effect size per category from any given study, then the category-specific estimates are merely a reindexing of the effect sizes in each study.

3.1. Subgroup Analysis

After subgrouping or shifting the unit of analysis, one analytic approach is to fit separate models based on the estimates from each category. Following this approach, the meta-analyst would fit the random effects meta-regression

{\tilde{T}}_{g j} = x_{j} β_{g} + u_{g j} + e_{g j},

using only the estimates from category g (i.e., $j \in J_{g}$ ). In Model (15), $β_{g}$ is a vector of meta-regression coefficients for subgroup g; $u_{g j}$ is a random effect for study j, assumed to be normally distributed with mean zero and variance $τ_{g}^{2}$ ; and $e_{g j}$ is a normally distributed sampling error for the averaged effect size estimate from category g of study j, with mean zero and $Var (e_{g j}) = {\tilde{V}}_{g j}$ . Let $L_{g}^{S} (β_{g}, τ_{g})$ and $R_{g}^{S} (τ_{g})$ denote the likelihood and restricted likelihood of Model (15) for aggregated effect size estimates from subgroup g, with generalized least squares estimators ${\hat{β}}_{S g} (τ_{g})$ , for $g = 1, ..., G$ .

Fitting separate models to each subgroup of effect size estimates has the disadvantage that it does not provide any easy way of making formal comparisons across categories. However, the approach described above can be represented as a single model, following the principles of composite marginal likelihood (Chen et al., 2014; Varin, 2008; Varin et al., 2011). Following the composite likelihood approach, one can estimate a single model using all of the data, but allowing for distinct parameters corresponding to each category of effects. The specification then becomes

{\tilde{T}}_{g j} = {\overset{°}{x}}_{g j} β + u_{g j} + e_{g j},

where ${\overset{°}{x}}_{g j} = P_{j} \otimes x_{j}$ is the set of predictors $x_{j}$ interacted with indicators for each of the G categories; $β = {[β_{g}]}_{g = 1}^{G}$ is the set of regression coefficients across all categories; and the random effects have variances that differ by category, with $Var (u_{g j}) = τ_{g}^{2}$ ; and $Var (e_{g j}) = {\tilde{V}}_{g j}$ . Model (16) is an instance of a random effects location-scale model (Viechtbauer & López-López, 2022), where the scale parameters are a function of indicator variables for each category.

The composite model representation is useful in that it allows for direct comparison of meta-regression coefficient estimates from different subgroups. This might seem counterintuitive because the dependence structure of the effect size estimates is misspecified—although aggregating removes the dependence of effect size estimates within each category, there may still be dependence between the effect size estimates across categories. However, RVE methods can be applied to account for this misspecification and provide valid assessments of uncertainty regarding the regression coefficients, even for comparisons between regression coefficients from different categories (Pustejovsky & Tipton, 2022). Chen et al. (2014) investigated this strategy in the context of multivariate meta-analysis where information about the correlations between effect size estimates is not available. Their simulation results showed that the approach provides well-calibrated standard errors and confidence intervals. We discuss RVE further in Section 4.

Let $L^{C G} (β, τ_{1}, \dots, τ_{G})$ denote the likelihood of Model (16) based on the subgroup-aggregated effect size estimates from across all categories $g = 1, .., G$ . Denote the corresponding restricted likelihood as $R^{C G} (τ_{1}, \dots, τ_{G})$ and the generalized least squares estimator as ${\hat{β}}_{C G} (τ_{1}, \dots, τ_{G})$ . The following Theorem states that the composite model likelihood is equivalent to the likelihoods corresponding to each subgroup model. We provide a proof in Online Appendix B.

Theorem 1: $L^{C G} (β, τ_{1}, \dots, τ_{G}) = \prod_{g = 1}^{G} L_{g}^{S} (β_{g}, τ_{g})$ and $R^{C G} (τ_{1}, \dots, τ_{G}) = \prod_{g = 1}^{G} R_{g}^{S} (τ_{g})$ . Further, ${\hat{β}}_{C G} (τ_{1}, \dots, τ_{G}) = {[{\hat{β}}_{S g} (τ_{g})]}_{g = 1}^{G}$ .

Theorem 1 implies that full or restricted maximum likelihood estimators based on the composite Model (16) are identical by definition to those based on Model (15) for each subgroup; model-based standard errors and confidence intervals are also identical across specifications. We show in Section 4 that robust variance estimators based on Model (15) are identical to those based on Model (16).

3.2. Ignoring Cross-Category Dependence

The strategy of shifting the unit of analysis and then estimating a meta-regression of aggregated effect size estimates from each category is not limited to models with distinct parameters for each category of effects (i.e., models with separable parameter sets). Rather, a meta-analyst might want to estimate models in which some parameters are category-specific while others are shared across categories. For example, they might be willing to assume that the category-specific random effects in Model (16) have the same variance across all categories, so that $Var (u_{g j}) = τ^{2}$ rather than $τ_{g}^{2}$ . As another example, they might wish to specify a meta-regression that includes category-specific intercepts along with additional control variables, the coefficients for which are constant across categories.

These specifications can be described by the model

{\tilde{T}}_{g j} = x_{g j} β + z_{g j} u + e_{g j},

where $x_{g j}$ is a generic set of predictors for the summary effect in category g of study j; $β$ is a set of regression coefficients; $z_{g j}$ is a $1 \times q$ row vector describing the random effects structure for the summary effect in category g of study j; and $u$ is a $q \times 1$ vector of random effects with mean zero and variance matrix $Φ (θ)$ . The sampling errors $e_{g j}$ are normally distributed with mean zero and variance ${\tilde{V}}_{g j}$ and are assumed to be to be mutually independent and independent of $u$ . Let $L^{I} (β, θ)$ and $R^{I} (θ)$ denote the full and restricted likelihoods of Model (17) based on the subgroup-aggregated effect size estimates, with generalized least squares estimator ${\hat{β}}_{I} (θ)$ .

Similar to Model (9), the random effects specification in Model (17) encompasses a broad set of possible structures, including hierarchical random effects models (Konstantopoulos, 2011), random effects location-scale model (Viechtbauer & López-López, 2022), phylogenetic meta-analytic models (Hadfield & Nakagawa, 2010; Lajeunesse, 2009), and cross-classified random effects models (Fernández-Castilla et al., 2019). Model (17) includes Model (16) as a special case, in which the parameters $β$ and $θ$ are separable. It also includes Models (7) and (9) as special cases, in which the data include only a single group ( $G = 1$ ).

If the category-specific average effect size estimates are calculated using the inverse-variance weights given in Equations 11 and 12, then Model (17) is equivalent to a model for the raw data that treats effect size estimates from different categories as independent. Let ${\overset{⌣}{x}}_{i j}$ and ${\overset{⌣}{z}}_{i j}$ denote the same sets of predictors as in (17), where the predictor for effect size i in study j is determined from the category into which i was classified. These can be calculated as

{\overset{⌣}{x}}_{i j} = \sum_{g = 1}^{G} c_{i j}^{g} x_{g j} and {\overset{⌣}{z}}_{i j} = \sum_{g = 1}^{G} c_{i j}^{g} z_{g j} .

Let ${\overset{⌣}{S}}_{j}$ denote the $k_{j} \times k_{j}$ matrix of sampling variances and covariances for the estimates from study j, which is a modification of $S_{j}$ , where estimates from different subgroups are assumed to be uncorrelated, so that $cor (e_{h j}, e_{i j}) = ρ_{h i j} \sum_{g = 1}^{G} c_{h j}^{g} c_{i j}^{g}$ . The modified sampling variance matrix can be calculated as

{\overset{⌣}{S}}_{j} = \sum_{g = 1}^{G} C_{j}^{g'} C_{j}^{g} S_{j} C_{j}^{g'} C_{j}^{g} .

Now, consider the model

T_{i j} = {\overset{⌣}{x}}_{i j} b + {\overset{⌣}{z}}_{i j} u + e_{i j},

where the sampling errors for study j, $e_{j}$ , are normally distributed with mean zero and variance ${\overset{⌣}{S}}_{j}$ , and where the random effects are normally distributed with mean zero and variance matrix $Φ (θ)$ , just as in Model (11). Following the terminology of Pustejovsky and Tipton (2022), this specification is a “subgroup-correlated effects” model. Let $L^{S C E} (β, θ)$ and $R^{S C E} (θ)$ denote the full and restricted likelihoods of Model (18) based on the raw effect size estimates, with generalized least squares estimator ${\hat{β}}_{S C E} (θ)$ .

The following theorem describes the equivalence between the likelihoods of Model (17) and Model (18). We provide a proof in Online Appendix C.

Theorem 2: If ${\tilde{T}}_{g j}$ and ${\tilde{V}}_{g j}$ are calculated using the inverse-variance weights from Equations 11 and 12, then $L^{S C E} (β, θ) = κ \times L^{I} (β, θ)$ and $R^{S C E} (θ) = κ \times R^{I} (θ)$ , where $κ$ is a constant that does not depend on the parameters $β$ or $θ$ , and ${\hat{β}}_{I} (θ) = {\hat{β}}_{S C E} (θ)$ .

The implication of Theorem 2 is that a meta-regression model estimated after shifting the unit of analysis with categories $g = 1, ..., G$ is equivalent to the subgroup correlated effects model for the raw effect size estimates. This equivalence holds even for models with complex random effects structures at the subgroup level or study level, so long as the predictors do not vary across estimates within the same subgroup of the same study. Because of its equivalence with the subgroup correlated effects model, the shifting unit of analysis technique can be understood not only as a reductive strategy but also as a multivariate model for the raw effect size estimates, which facilitates comparison to other multivariate models.

Theorem 2 also establishes the equivalence between Models (7) and (8) and between Models (9) and (10) because all of these models are special cases with $G = 1$ .

3.3. Another Illustration

Roberts et al. (2022) reported a synthesis of experimental studies of the enactment effect, the phenomenon that physically performed action enhances memory of a word or phase more strongly than simply reading it. We reanalyzed the behavioral studies synthesized by Roberts et al. (2022), comprised of 145 studies and 443 effects, in order to demonstrate the equivalencies described in Theorems 1 and 2. Data and code for replicating all calculations are available at https://osf.io/9u52e/.

The original authors investigated multiple moderator variables, including the comparison task, test format, and use of objects. For illustrative purposes, we examined difference in effect size based on the type of comparison task. The data include 33 studies with 98 effects involving experimenter-performed tasks (EPT), 17 studies with 29 effects involving imagery tasks (IT), and 123 studies with 316 effects involving verbal tasks (VT); 28 of the studies included effects from two distinct types of comparison tasks. The number of effects reported in each study ranges from 1 to 12 for EPT, from 1 to 4 for IT, and from 1 to 8 for VT, leading to the dependence within the comparison task subgroups within some studies. Considering this dependence structure, a meta-analyst might follow the shifting unit-of-analysis strategy and aggregate effect sizes for each comparison task within each study.

Table 2 reports results from several different approaches to the moderator analysis. In the first column, we calculated subgroup-aggregated effect size estimates for each comparison task type in each study, assuming a constant correlation of $ρ = 0.80$ between effect size estimates from the same subgroup and same study. We then conducted subgroup analysis, estimating separate random effects models for each comparison task type. This yielded the estimates of average effect size and between-study heterogeneity for each subgroup, as well as log-likelihood and Q statistics for each subgroup.

Table 2.

Comparison of Aggregated Effects Subgroup Analysis, Aggregated Effects Composite Model, and Subgroup Correlated Effects (SCE) Model

Coefficient	Aggregation + Subgroup Analysis			Aggregation + Composite Model			Multivariate SCE Model
Coefficient	Est	SE	95% CI	Est	SE	95% CI	Est	SE	95% CI
Average effect size
EPT	0.423	.100	[0.220, 0.626]	0.423	.100	[0.226, 0.619]	0.423	.100	[0.226, 0.619]
IT	0.636	.159	[0.298, 0.973]	0.636	.159	[0.321, 0.950]	0.636	.159	[0.321, 0.950]
VT	1.010	.066	[0.878, 1.142]	1.010	.066	[0.879, 1.141]	1.010	.066	[0.879, 1.141]
Between-study SD ( $τ$ )
EPT	0.483		[0.353, 0.727]	0.483		[0.332, 0.690]	0.483		[0.332, 0.690]
IT	0.540		[0.314, 0.881]	0.540		[0.321, 0.884]	0.540		[0.321, 0.884]
VT	0.656		[0.585, 0.817]	0.656		[0.555, 0.777]	0.656		[0.555, 0.777]
Log likelihood				−185.166			−1,316.356
EPT	−28.571
IT	−15.641
VT	−140.954
Q				794.294			3,419.316
EPT	114.805
IT	55.315
VT	624.174

Note. EPT = experimenter-performed tasks; IT = imagery task; VT = verbal task.

In the second column, we estimated average effect sizes for each subgroup using data from all three comparison tasks in a single composite model. The model included separate intercepts for each task type and allowed the between-study heterogeneity to differ by task type (i.e., a location-scale model). Because the estimates are based on a single model, we obtained only a single log-likelihood, equal to the sum of the log-likelihoods for the three task types reported in the first column. Similarly, we obtained a single Q statistic that is the sum of the Q statistics from the subgroup-specific models.

The point estimates and standard errors from the composite model are equal to those from the subgroup-specific models, as implied by Theorem 1. However, the confidence intervals for the regression coefficients differ slightly because they are calculated using different degrees of freedom. The degrees of freedom for the t tests in the subgroup analyses are equal to one less than the total number of studies reporting effect sizes in a given category, whereas those for the composite model are 170, equal to three less than the total number of studies across all three subgroup categories after shifting the unit of analysis. The profile-likelihood confidence intervals for between-study heterogeneity also differ slightly. The confidence intervals in the first column are calculated by profiling the subgroup-specific log likelihoods, whereas those in the second column are calculated by profiling the combined likelihood, a task that is more computationally demanding and that does not take advantage of the separable structure of the log-likelihood.

In the third column of Table 2, we estimated a multivariate subgroup correlated effects model using the raw effect size estimates. Consistent with Theorem 2, the point estimates, standard errors, and confidence intervals match those from the composite model estimated on the subgroup aggregated effects. (One can obtain confidence intervals exactly equivalent to those from the subgroup analysis in the first column by manually setting the dfs argument in the rma.mv function.) The log likelihood estimates and Q statistics from the SCE model differ from the corresponding statistics from the composite model for the aggregated effects, but the differences are constants that do not depend on the parameters of the model.

As noted in Section 2.3, the discrepancy in the Q statistics can be interpreted as a measure of heterogeneity in the effect size estimates over which we have aggregated. Here, the large difference in Q statistics indicates that there is substantial excess heterogeneity in effect sizes drawn from the same subgroup and same sample ( $p < 0.001$ compared to a $χ^{2}$ distribution with $K - J = 270$ degrees of freedom). This suggests that the meta-analyst should consider a model that allows for within-subgroup variation in effect size (or perhaps that includes effect-level predictors that could explain such variation) rather than aggregating over heterogeneous effects.

4. Robust Variance Estimation

Many of the models that we have described are subject to the criticism that they involve arbitrary or potentially misspecified assumptions about the sampling errors or other aspects of the dependence structure. For instance, if primary studies fail to report correlations among outcome measures, the correlation used to calculate the sampling variance of the aggregated effect sizes based on Equation 6 or 14 might be selected based only on convention. If the assumed correlation is incorrect, then V_j will also be incorrect, which will affect the estimation of the heterogeneity variance components, as well as the standard errors and confidence intervals for meta-regression coefficients. Further, if a meta-analyst uses the shifting unit of analysis approach and estimates a model that includes effect sizes from multiple categories (as in Equations 16 or 17), then the dependence structure of the model will be misspecified due to the correlation between effect size estimates from different categories in the same study.

To address the issues of potential mis-specification, Hedges et al. (2010) proposed RVE methods for meta-regression models. RVE methods provide means to estimate standard errors, conduct hypothesis tests, and construct confidence intervals for meta-regression coefficients that do not rely on the assumed dependence structure being correctly specified. Instead, RVE uses sandwich estimators that work under the weaker assumption that the effect size estimates can be grouped into clusters of dependent observations, where effect sizes from distinct clusters are independent. The original form of RVE described by Hedges et al. (2010) involves approximations that require a large number of independent clusters. However, subsequent work (Tipton, 2015; Tipton & Pustejovsky, 2015) has provided small-sample corrections that reduce the bias of the variance estimators and improve the calibration of hypothesis tests and confidence intervals when the data have a limited number of clusters. RVE methods are therefore an attractive tool for inference in meta-regression models with dependent effect sizes.

4.1. Subgroup Correlated Effects

The equivalence between models for aggregated effect size estimates and models for disaggregated, dependent effect sizes extends to RVE methods, including methods that incorporate small-sample corrections as described by Hedges et al. (2010) and Tipton (2015). We now demonstrate this equivalence by showing the equality of RVE applied to Model (17) for effect sizes aggregated within subgroups and RVE applied to the subgroup correlated effects Model (18) for the disaggregated effect sizes.

Suppose that studies can be grouped into M independent clusters, where $S_{m}$ denotes the set of studies in cluster m for $m = 1, ..., M$ . The clusters could be identical to studies, in which case $M = J$ and $S_{m} = {m}$ . Let ${\bar{T}}_{j} = {[{\tilde{T}}_{g j}]}_{g \in G_{j}}$ be the vector of subgroup aggregated effect size estimates from study j and ${\bar{T}}_{m} = {[{\bar{T}}_{j}]}_{j \in S_{m}}$ . Let ${\bar{X}}_{j} = {[x_{g j}]}_{g \in G_{j}}$ and ${\bar{X}}_{m} = {[{\bar{X}}_{j}]}_{j \in S_{m}}$ be the corresponding matrices of predictors, and let ${\bar{Z}}_{j} = {[z_{g j}]}_{g \in G_{j}}$ and ${\bar{Z}}_{m} = {[{\bar{Z}}_{j}]}_{j \in S_{m}}$ be the matrices of stacked row-vectors of random effects predictors. Let ${\bar{Ω}}_{m} = {\bar{Z}}_{m} Φ (\hat{θ}) {\bar{Z}}_{m}^{'} + {\bar{V}}_{m}$ , where $\hat{θ}$ is the maximum likelihood or restricted maximum likelihood estimate of the variance components and ${\bar{V}}_{m} = \oplus_{j \in S_{m}} \oplus_{g \in G_{j}} {\tilde{V}}_{g j}$ is a diagonal matrix of the sampling variances. Let ${\bar{W}}_{m} = {\bar{Ω}}_{m}^{- 1}$ , ${\bar{e}}_{m} = {\bar{T}}_{m} - {\bar{X}}_{m} {\hat{β}}_{I} (\hat{θ})$ , and ${\bar{D}}_{m}$ be any square root of ${\bar{Ω}}_{m}$ , so that ${\bar{Ω}}_{m} = {\bar{D}}_{m}^{'} {\bar{D}}_{m}$ . Let ${\bar{I}}_{m}$ be an identity matrix that has the same dimension as ${\bar{Ω}}_{m}$ .

A robust variance estimator for Model (17) is given by

V_{I}^{(a)} = \bar{B} (\sum_{m = 1}^{M} {\bar{X}}_{m}^{'} {\bar{W}}_{m} {\bar{A}}_{m}^{(a)} {\bar{e}}_{m} {\bar{e}}_{m}^{'} {\bar{A}}_{m}^{(a)} {\bar{W}}_{m} {\bar{X}}_{m}) \bar{B},

where $\bar{B} = {(\sum_{m = 1}^{M} {\bar{X}}_{m}^{'} {\bar{W}}_{m} {\bar{X}}_{m})}^{- 1}$ and ${\bar{A}}_{1}^{(a)}, \dots, {\bar{A}}_{M}^{(a)}$ are a set of small-sample adjustment matrices. Several different forms of adjustment matrices have been proposed in the literature. The forms described by Hedges et al. (2010) are ${\bar{A}}_{m}^{(0)} = {\bar{I}}_{m}$ , which makes no small-sample adjustment, and ${\bar{A}}_{m}^{(1)} = \sqrt{\frac{M}{M - 1}} {\bar{I}}_{m}$ , which makes a rough correction based on the number of unique clusters. Tipton (2015) considered adjustments of the form

\begin{matrix} {\bar{A}}_{m}^{(2)} = {\bar{D}}_{m}^{'} {[{\bar{D}}_{m} ({\bar{Ω}}_{m} - {\bar{X}}_{m} \bar{B} {\bar{X}}_{m}^{'}) {\bar{D}}_{m}^{'}]}^{- 1 / 2} {\bar{D}}_{m}, \\ {\bar{A}}_{m}^{(3)} = {({\bar{I}}_{m} - {\bar{X}}_{m}^{'} \bar{B} {\bar{X}}_{m}^{'} {\bar{W}}_{m})}^{- 1}, \end{matrix}

where ${\bar{A}}_{m}^{(2)}$ is derived from the bias-reduced linearization method first proposed by McCaffrey et al. (2001) and ${\bar{A}}_{m}^{(3)}$ is an approximation to a jackknife variance estimator.

Let ${\overset{⌣}{T}}_{m} = {[T_{j}]}_{j \in S_{m}}$ be the vector of disaggregated effect size estimates included in cluster m, with corresponding predictor matrices ${\overset{⌣}{X}}_{j} = {[{\overset{⌣}{x}}_{i j}]}_{i = 1}^{k_{j}}$ and ${\overset{⌣}{X}}_{m} = {[{\overset{⌣}{X}}_{j}]}_{j \in S_{m}}$ . Similarly, let ${\overset{⌣}{Z}}_{j} = {[{\overset{⌣}{z}}_{i j}]}_{i = 1}^{k_{j}}$ and ${\overset{⌣}{Z}}_{m} = {[{\overset{⌣}{Z}}_{j}]}_{j \in S_{m}}$ be the matrices of random effects predictors for $j = 1, ..., J$ and $m = 1, ..., M$ , respectively. Let

{\overset{⌣}{Ω}}_{m} = {\overset{⌣}{Z}}_{m} Φ (\hat{θ}) {\overset{⌣}{Z}}_{m}^{'} + \underset{j \in S_{m}}{\oplus} {\overset{⌣}{S}}_{j},

an estimate of the variance of the errors in the effect size estimates from cluster m based on Model (18). Let

{\overset{⌣}{W}}_{m} = {\overset{⌣}{Ω}}_{m}^{- 1}

{\overset{⌣}{e}}_{m} = {\overset{⌣}{T}}_{m} - {\overset{⌣}{X}}_{m} {\hat{β}}_{S C E} (\hat{θ})

, and

{\overset{⌣}{D}}_{m}

be any square root of

{\overset{⌣}{Ω}}_{m}

, so that

{\overset{⌣}{Ω}}_{m} = {\overset{⌣}{D}}_{m}^{'} {\overset{⌣}{D}}_{m}

. Let

{\overset{⌣}{I}}_{m}

be an identity matrix of the same dimension as

{\overset{⌣}{Ω}}_{m}

Let $V_{S C E}^{(a)}$ denote a robust variance estimator for Model (18). $V_{S C E}^{(a)}$ is calculated just as in Equation 19, but substituting ${\overset{⌣}{X}}_{m}$ for ${\bar{X}}_{m}$ , ${\overset{⌣}{W}}_{m}$ for ${\bar{W}}_{m}$ , ${\overset{⌣}{e}}_{m}$ for ${\bar{e}}_{m}$ , and ${\overset{⌣}{A}}_{m}^{(a)}$ for ${\bar{A}}_{m}^{(a)}$ , where ${\overset{⌣}{A}}_{m}^{(0)} = {\overset{⌣}{I}}_{m}$ , ${\overset{⌣}{A}}_{m}^{(1)} = \sqrt{\frac{M}{M - 1}} {\overset{⌣}{I}}_{m}$ ,

\begin{matrix} {\overset{⌣}{A}}_{m}^{(2)} = {\overset{⌣}{D}}_{m^{'}} {[{\overset{⌣}{D}}_{m} ({\overset{⌣}{Ω}}_{m} - {\overset{⌣}{X}}_{m} \overset{⌣}{B} {\overset{⌣}{X}}_{m}^{'}) {\overset{⌣}{D}}_{m}^{'}]}^{- 1 / 2} {\overset{⌣}{D}}_{m} \\ {\overset{⌣}{A}}_{m}^{(3)} = {({\overset{⌣}{I}}_{m} - {\overset{⌣}{X}}_{m} \overset{⌣}{B} {\overset{⌣}{X}}_{m}^{'} {\overset{⌣}{W}}_{m})}^{- 1} \end{matrix}

and $\overset{⌣}{B} = {(\sum_{m = 1}^{M} {\overset{⌣}{X}}_{m}^{'} {\overset{⌣}{W}}_{m} {\overset{⌣}{X}}_{m})}^{- 1}$ .

Under quite general conditions, $V_{S C E}^{(a)}$ and $V_{I}^{(a)}$ are identical for $a = 0, 1, and 3$ . For $a = 2$ , the two forms are identical under more restrictive conditions pertaining to the sampling variance matrices for each subgroup in each study and the number of effect sizes in each subgroup of each study in a given cluster. The following theorem gives a formal statement of the equivalence of $V_{S C E}^{(a)}$ and $V_{I}^{(a)}$ . We give a proof in Online Appendix D.

Theorem 3: If ${\tilde{T}}_{g j}$ and ${\tilde{V}}_{g j}$ are calculated using the inverse-variance weights from Equations 11 and 12 and if both Models (17) and (18) are estimated using maximum likelihood or both estimated using restricted maximum likelihood, then $V_{I}^{(a)} = V_{S C E}^{(a)}$ for $a = 0, 1, and 3$ . If, additionally, ${\tilde{S}}_{g j} = C_{j}^{g} S_{j} C_{j}^{g'}$ has a compound symmetric structure for $g \in G_{j}$ and $j = 1, ..., J$ and if $k_{g j} = k_{m}$ is constant for $g \in G_{j}$ and $j \in S_{m}$ , then $V_{I}^{(2)} = V_{S C E}^{(2)}$ .

Equivalencies of RVE applied to several of the other models we have described can be established as special cases of Theorem 3. Specifically, RVE applied to Model (7) is equal to RVE applied to Model (8). Similarly, RVE applied to Model (9) is equal to RVE applied to Model (10).

Somewhat restrictive conditions are required for the RVE matrices with small-sample correction $a = 2$ to be equal. The compound symmetry condition will be satisfied when the effect size estimates within a given subgroup of a given study have equal sampling variances ( $s_{h j}^{2} = s_{i j}^{2}$ when $c_{h j}^{g} \times c_{i j}^{g} = 1$ ) and equal sampling correlation ( $ρ_{h i j} = ρ_{g j}$ when $c_{h j}^{g} \times c_{i j}^{g} = 1$ ). In many applications, this condition will be met in practice because the meta-analyst will need to make simplifying assumptions about the dependence structure of effect sizes, such as by assuming a single, common correlation between pairs of effect size estimates within a given study (Pustejovsky & Tipton, 2022).

The second condition requiring that $k_{g j} = k_{m}$ is more restrictive and artificial. It will be more readily satisfied if the clustering variable is identical to the study identifier, so $S_{m} = m$ . The condition will then be met, for instance, if effect sizes are aggregated by study (without subgrouping) or if the subgroups include at most one effect size estimate per study. Otherwise, the small-sample corrected RVE estimator based on the raw effect size estimates will not necessarily be exactly equal to the estimator based on the aggregated effects or subgroup-aggregated effects.

Online Appendix F demonstrates Theorem 3 by applying the various forms of RVE to the empirical examples from Sections 2.3 and 3.3. In both examples, the sampling variances are not compound symmetric, leading to small discrepancies between the $a = 2$ forms based on the aggregated effect sizes versus those based on the multivariate models.

4.2. Subgroup-Specific and Composite Models

Theorem 3 does not cover the relationship between the subgroup-specific Model (15) and the composite Model (16). Model (16) can be understood as a special case of Model (17), in which the regression coefficients and variance components are separable. More precisely, the composite meta-regression has predictors of the form ${\bar{X}}_{m} = {[P_{j} \otimes x_{j}]}_{j \in S_{m}}$ , where $x_{j}$ is a row vector of study-level predictors, and thr estimated marginal variance matrix for cluster m equals to

{\bar{Ω}}_{m} = (\underset{j \in S_{m}}{\oplus} P_{j} Φ (\hat{θ}) P_{g}^{'}) + {\bar{V}}_{m},

where $Φ (\hat{θ}) = \oplus_{g = 1}^{G} {\hat{τ}}_{g}^{2}$ . Let $V_{C G}^{(a)}$ denote a robust variance estimator for Model (16), which is calculated just as in Equation 19.

Applying RVE to the subgroup-specific model entails using only the averaged effect sizes in that category. Let $V_{S g}^{(a)}$ denote a robust variance estimator for Model (15) using the averaged effect sizes in category g. Let M_g denote the number of clusters that include one or more averaged effect sizes in category g, for $g = 1, ..., G$ . Let ${⌈ V_{C G}^{(a)} ⌋}_{g g}$ denote the $p \times p$ submatrix of $V_{C G}^{(a)}$ corresponding to the coefficients for group g, that is, ${⌈ V_{C G}^{(a)} ⌋}_{g g} = (p_{g} \otimes I_{p}) V_{C G}^{(a)} (p_{g}^{'} \otimes I_{p})$ , where $I_{p}$ is a $p \times p$ identity matrix. The following theorem establishes the equivalence between RVE applied to the subgroup-specific model and RVE applied to the composite model. We give a proof in Online Appendix E and an illustrative example in Online Appendix F.

Theorem 4: If the predictors in Model (16) have the form ${\bar{X}}_{m} = {[P_{j} \otimes x_{j}]}_{j \in S_{m}}$ for $m = 1, ..., M$ and if the marginal variance matrices have the form given in (20), then ${⌈ V_{C G}^{(a)} ⌋}_{g g} = V_{S g}^{(a)}$ for $a = 0, 2, and 3$ and $g = 1, ..., G$ . If these conditions hold and $M_{g} = M$ for $g = 1, ..., G$ , then ${⌈ V_{C G}^{(1)} ⌋}_{g g} = V_{S g}^{(1)}$ for $g = 1, ..., G$ .

The implication of Theorem 4 is that applying RVE (with the $a = 0$ , $a = 2$ , or $a = 3$ small-sample correction) to the composite model recovers the same result as applying RVE to the component subgroup model for a single category. Standard errors for entries in ${\hat{β}}_{g}$ and covariances for pairs of entries in ${\hat{β}}_{g}$ will be the same, regardless of which model is used. Notably, however, using the composite model provides the further possibility of calculating robust estimates of the covariances between regression coefficients for different groups, which are not readily available when applying RVE to the component subgroup models. These covariances are needed for drawing inferences about contrasts between different subgroups (e.g., comparing the average effect size for group $g = 1$ vs. group $g = 2$ or conducting a joint test for equality of average effect sizes across groups $g = 1, ..., G$ ). The possibility of conducting inference on such contrasts is one of the primary advantages of the composite model representation.

5. Discussion

In prior discussions of methods for meta-analysis of dependent effect size estimates, reductionist strategies such as aggregating effect sizes by study or shifting the unit of analysis have been characterized as “ad hoc” techniques (Becker, 2000). This may leave researchers with the impression that reductionist strategies are simply outdated or inferior to more recently developed integrative approaches, such as multilevel meta-analysis or RVE—strategies which might be interpreted as more advanced or sophisticated. Such casual impressions are belied by the theorems that we have presented, which show that reductionist methods are, in fact, exactly equivalent to certain multivariate models for dependent effect sizes.

Equivalencies between different representations of a model have proven to be useful in several related areas of meta-analysis. For example, Cheung (2013, 2014) demonstrated how a range of meta-analytic models can be expressed as structural equation models. Recognizing such connections is useful for estimating complex random effects models and for handling missingness in meta-regression predictors. As a further example, White et al. (2012) showed how inconsistency models for network meta-analysis can be represented as multivariate meta-regression models, which can be fit using standard software tools. Just as in these examples, the connections we have described between reductionist strategies and multivariate models provide both conceptual and practical insights.

Conceptually, the correspondence between reductionist strategies and multivariate models provides a unifying framework for meta-analysis of dependent effect size estimates. Rather than treating reductive strategies as separate from and incommensurable with newer models and methods, Theorems 1 and 2 show that reductive strategies can be understood as specific forms of multivariate models. Namely, aggregating effect sizes by study is equivalent to using a correlated effects model and subgrouping or shifting the unit of analysis is equivalent to using a subgroup correlated effects model. Under the conditions we have described, the equivalencies apply to likelihood-based estimation methods and to RVE methods. Theorems 3 and Theorem 4 imply that reductionist strategies can be interpreted as working models, which can be combined with RVE methods to provide protection from model misspecification. Thus, for purposes of developing further theory and evaluating the properties of different approaches to meta-analysis of dependent effect sizes, it is sufficient to consider the framework of multivariate working models.

The multivariate representations of reductionist approaches are also useful for the purposes of model comparison. As demonstrated in the examples from Sections 2.3 and 3.3, the difference in Q statistics between models based on raw effect sizes and models based on aggregated effect sizes provides a measure for evaluating the appropriateness of strategies that involve aggregating effect size estimates. Alternately, likelihood ratio tests can be used to compare the correlated effects model to a model that also includes within-study heterogeneity. Because of their equivalence, a test of the correlated effects model reflects directly on the appropriateness of aggregating effects by study. Similarly, one could evaluate the appropriateness of shifting the unit of analysis by testing the subgroup correlated effects model against a model that allows for within-study, within-subgroup heterogeneity.

Comparisons among multivariate working models might be of particular interest in syntheses where different subgroups of effects can be conceptualized as a multivariate outcome. For instance, in a synthesis of intervention studies on second language learning, researchers might conceptualize outcomes as consisting of several dimensions, such as listening, speaking, and writing abilities. Currently, researchers often use the ad hoc methods of subgrouping or shifting the unit of analysis to analyze such data, conducting what amount to separate meta-analyses of available data from each dimension. Recognizing the connection between these approaches and the corresponding multivariate model, the subgroup correlated effects model, is useful because it can facilitate comparisons with other multivariate models. Rather than treating different subgroups of effects as independent, researchers might entertain a model in which the study-specific effect sizes from each subgroup follow a correlated, multivariate random effects distribution (Jackson et al., 2011; Sera et al., 2019). Further research is needed to understand the advantages and drawbacks of such multivariate random effects models relative to models that treat different dimensions as independent—especially for meta-analytic data where not every study includes findings for every dimension.

On a practical level, aggregating effect sizes, subgrouping, or shifting the unit of analysis might be appealing because they make an analysis easier to follow, because using them makes it easier to run diagnostics or create illustrations of the results, or because of software limitations. If the corresponding multivariate working model is appropriate, then the results that we have presented provide formal justification for using such reductionist strategies. For instance, a graphical representation of meta-analysis results such as a forest plot could be based on the aggregated effect size estimates rather than the raw estimates, reducing the complexity of the figure. Carrying out calculations on aggregated or subgroup-aggregated effect size estimates will also tend to be more computationally efficient than fitting the corresponding multivariate model. Although differences in computation time might be trivial for estimation of summary meta-analytic models with moderately sized datasets, computational efficiency might matter in practice for syntheses of large datasets and for more demanding calculations such as sensitivity analyses, profile likelihood confidence regions, or bootstrap hypothesis tests (e.g., Du et al., 2023; Joshi et al., 2022), or applications of multiple imputation for missing predictor values (Lee & Beretvas, 2023).

The equivalencies that we have presented are defined under the premise that there is no missingness in the predictor variables included in the meta-regression model. In practice, it is common that the predictors of interest may be available only for a subset of included studies. Recent methodological research has examined various techniques for handling missing predictors, including complete-case analysis, shifting-case analysis, multiple imputation, and full information maximum likelihood (Cheung, 2008, 2014; Lee & Beretvas, 2023; Schauer et al., 2022). If the complete-case analysis or shifting-case analysis method is used to remove the effect sizes with missing predictors, the equivalencies between ad hoc strategies and multivariate models will hold for the subset of data with completely observed predictors. However, further research is needed to examine different representations of meta-regression models with missing predictors if multiple imputation or full information maximum likelihood (Cheung, 2008, 2014) methods are used to account for missingness. The equivalencies that we have described in the present article also highlight the need for missing data methods that work under the multilevel structure inherent in meta-analysis of dependent effect sizes.

The correspondences between models that we have described are limited in several respects. First, exact equivalence between reductionist approaches and the corresponding multivariate models holds only for estimation methods that are based on the likelihood or restricted likelihood of the models. Some other commonly used methods are not based on the model likelihood, such as moment estimators for variance components, confidence intervals based on generalized Q-statistics (Jackson, 2013), or the Hartung and Knapp (2001) correction to standard errors. They will not necessarily be equivalent (or even defined) across both representations. Further research is needed to better understand the relative strengths of likelihood-based estimation versus moment-estimation methods for models with dependent effect sizes, such as those proposed by Hedges et al. (2010) and Jackson et al. (2013).

Second, exact equivalence with multivariate models is limited to reductionist approaches that involve aggregating using inverse-variance weighted averages, with sampling variances calculated based on an assumed correlation between effect size estimates. Some meta-analysts might use more casual approximations for calculating aggregated effect sizes and variances, such as taking simple averages of effect size estimates and sampling variances by study. This could lead to difference between results based on aggregated effect size estimates and results based on the multivariate model, although we expect that such differences will usually be minor in practice.

Finally, the equivalencies that we have described are limited to meta-regression models where the predictors vary only at the level of aggregation—that is, at the study level when aggregating effects by study or at the subgroup level when aggregating effects by subgroup within study. A key advantage of integrative models is that they permit inclusion of predictors that vary across effect sizes within the same subgroup or same study, allowing the meta-analyst to control for potential methodological confounding factors, differences between outcome measures, or other features that may explain within-study variation (López-López et al., 2018; Tipton et al., 2019). There remains a need to investigate the benefits and limitations of using the correlated effects or subgroup correlated effects working models for meta-regressions that include such effect-level predictors.

The array of tools and approaches for synthesizing dependent effect sizes has expanded substantially over the past 15 years. As meta-analysts continue to work with data that involve complex dependence structures, there is a need for further guidance about how to select an appropriate modeling strategy. Recognizing the equivalences between reductionist strategies and multivariate working models will, we hope, help in establishing a more principled basis for such decisions.

Supplemental Material

Supplemental Material, sj-docx-1-jeb-10.3102_10769986241232524 - Equivalencies Between Ad Hoc Strategies and Multivariate Models for Meta-Analysis of Dependent Effect Sizes

Supplemental Material, sj-docx-1-jeb-10.3102_10769986241232524 for Equivalencies Between Ad Hoc Strategies and Multivariate Models for Meta-Analysis of Dependent Effect Sizes by James E. Pustejovsky and Man Chen in Journal of Educational and Behavioral Statistics

Footnotes

Authors’ Note

Code and data for replicating the empirical examples presented in this article are available on the Open Science Framework at

Declaration of Conflicting Interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) received no financial support for the research and/or authorship of this article.

ORCID iD

James E. Pustejovsky

References

Ahn

Ames

A. J.

Myers

N. D.

(2012). A review of meta-analyses in education: Methodological strengths and weaknesses. Review of Educational Research, 82(4), 436–476. https://doi.org/10.3102/0034654312458162

Becker

B. J.

(1992). Using results from replicated studies to estimate linear models. Journal of Educational Statistics, 17(4), 341–362. https://doi.org/10.3102/10769986017004341

Becker

B. J.

(2000). Multivariate meta-analysis. In Brown

S. D.

Tinsley

H. E. A.

(Eds.), Handbook of applied multivariate statistics and mathematical modeling (pp. 499–525). Academic Press. https://doi.org/10.1016/B978-012691360-6/50018-5

Berkey

C. S.

Hoaglin

D. C.

Antczak-Bouckoms

Mosteller

Colditz

G. A.

(1998). Meta-analysis of multiple outcomes by regression with random effects. Statistics in Medicine, 17, 2537–2550. https://doi.org/10.1002/(SICI)1097-0258(19981130)17:22%3C2537::AID-SIM953%3E3.0.CO;2-C

Borenstein

Hedges

L. V.

Higgins

J. P. T.

Rothstein

H. R.

(2009). Introduction to meta-analysis. John Wiley. https://doi.org/10.1002/9780470743386

Chen

Hong

Riley

R. D.

(2014). An alternative pseudolikelihood method for multivariate random-effects meta-analysis. Statistics in Medicine, 34(3), 361–380. https://doi.org/10.1002/sim.6350

Cheung

M. W.-L.

(2008). A model for integrating fixed-, random-, and mixed-effects meta-analyses into structural equation modeling. Psychological Methods, 13(3), 182–202. https://doi.org/10.1037/a0013163

Cheung

M. W.-L.

(2013). Multivariate meta-analysis as structural equation models. Structural Equation Modeling: A Multidisciplinary Journal, 20(3), 429–454. https://doi.org/10.1080/10705511.2013.797827

Cheung

M. W.-L.

(2014). Modeling dependent effect sizes with three-level meta-analyses: A structural equation modeling approach. Psychological Methods, 19(2), 211–229. https://doi.org/10.1037/a0032968

10.

Cheung

M. W.-L.

(2019). A guide to conducting a meta-analysis with non-independent effect sizes. Neuropsychology Review, 29(4), 387–396. https://doi.org/10.1007/s11065-019-09415-6

11.

Clinton-Lisell

(2021). Stop multitasking and just read: Meta-analyses of multitasking’s effects on reading performance and reading time. Journal of Research in Reading, 44(4), 787–816. https://doi.org/10.1111/1467-9817.12372

12.

Cooper

H. M.

(1998). Synthesizing research: A guide for literature reviews (3rd ed.). Sage.

13.

Cooper

H. M.

(2016). Research synthesis and meta-analysis (5th ed.). Sage.

14.

De Wit

Marsiske

O’Shea

Kessels

R. P.

Kurasz

A. M.

DeFeis

Schaefer

Smith

G. E.

(2021). Procedural learning in individuals with amnestic mild cognitive impairment and Alzheimer’s dementia: A systematic review and meta-analysis. Neuropsychology Review, 31, 103–114. https://doi.org/10.1007/s11065-020-09449-1

15.

Jiang

(2023). Bootstrap-based between-study heterogeneity tests in meta-analysis. Multivariate Behavioral Research, 58(3), 484–503. https://doi.org/10.1080/00273171.2021.1997701

16.

Emslander

Scherer

(2022). The relation between executive functions and math intelligence in preschool children: A systematic review and meta-analysis. Psychological Bulletin, 148(5–6), 337. https://doi.org/10.1037/bul0000369

17.

Fernández-Castilla

Aloe

A. M.

Declercq

Jamshidi

Beretvas

S. N.

Onghena

Van den Noortgate

(2020). Estimating outcome-specific effects in meta-analyses of multiple outcomes: A simulation study. Behavior Research Methods, 53, 702–717. https://doi.org/10.3758/s13428-020-01459-4

18.

Fernández-Castilla

Maes

Declercq

Jamshidi

Beretvas

S. N.

Onghena

Van den Noortgate

(2019). A demonstration and evaluation of the use of cross-classified random-effects models for meta-analysis. Behavior Research Methods, 51(3), 1286–1304. https://doi.org/10.3758/s13428-018-1063-2

19.

Fisher

Tipton

Zhipeng

(2017). Robumeta: Robust variance meta-regression (Version 2.0). https://github.com/zackfisher/robumeta

20.

Furenes

M. I.

Kucirkova

Bus

A. G.

(2021). A comparison of children’s reading on paper versus screen: A meta-analysis. Review of Educational Research, 91(4), 483–517. https://doi.org/10.3102/0034654321998074

21.

Gaskell

Simmonds-Buckley

Kellett

Stockton

Somerville

Rogerson

Delgadillo

(2023). The effectiveness of psychological interventions delivered in routine practice: Systematic review and meta-analysis. Administration and Policy in Mental Health and Mental Health Services Research, 50(1), 43–57. https://doi.org/10.1007/s10488-022-01225-y

22.

Gleser

L. J.

Olkin

(2009). Stochastically dependent effect sizes. In Cooper

Hedges

L. V.

Valentine

J. C.

(Eds.), The handbook of research synthesis and meta-analysis (2nd ed., pp. 357–376). Russell Sage Foundation.

23.

Graham

Hebert

Fishman

Ray

A. B.

Rouse

A. G.

(2020). Do children classified with specific language impairment have a learning disability in writing? A meta-analysis. Journal of Learning Disabilities, 53(4), 292–310. https://doi.org/10.1177/0022219420917338

24.

Hadfield

J. D.

Nakagawa

(2010). General quantitative genetic methods for comparative biology: Phylogenies, taxonomies and multi-trait models for continuous and categorical characters. Journal of Evolutionary Biology, 23(3), 494–508. https://doi.org/10.1111/j.1420-9101.2009.01915.x

25.

Hartung

Knapp

(2001). On tests of the overall treatment effect in meta-analysis with normally distributed responses. Statistics in Medicine, 20(12), 1771–1782. https://doi.org/10.1002/sim.791

26.

Hedberg

E. C.

(2014). ROBUMETA: Stata module to perform robust variance estimation in meta-regression with dependent effect size estimates. Statistical Software Components S457219. Boston College Department of Economics.

27.

Hedges

L. V.

Olkin

(1985). Statistical methods for meta-analysis. Academic Press.

28.

Hedges

L. V.

Tipton

Johnson

M. C.

(2010). Robust variance estimation in meta-regression with dependent effect size estimates. Research Synthesis Methods, 1(1), 39–65. https://doi.org/10.1002/jrsm.5

29.

Hong

Riley

R. D.

Chen

(2018). An improved method for bivariate meta-analysis when within-study correlations are unknown. Research Synthesis Methods, 9(1), 73–88. https://doi.org/10.1002/jrsm.1274

30.

Jackson

(2013). Confidence intervals for the between-study variance in random effects meta-analysis using generalised Cochran heterogeneity statistics. Research Synthesis Methods, 4(3), 220–229. https://doi.org/10.1002/jrsm.1081

31.

Jackson

Riley

R. D.

White

I. R.

(2011). Multivariate meta-analysis: Potential and promise. Statistics in Medicine, 30(20), 2481–2498. https://doi.org/10.1002/sim.4172

32.

Jackson

White

I. R.

Riley

R. D.

(2013). A matrix-based method of moments for fitting the multivariate random effects model for meta-analysis and meta-regression. Biometrical Journal, 55(2), 231–245. https://doi.org/10.1002/bimj.201200152

33.

Joshi

Pustejovsky

J. E.

Beretvas

S. N.

(2022). Cluster wild bootstrapping to handle dependent effect sizes in meta-analysis with a small number of studies. Research Synthesis Methods, 13(4), 457–477. https://doi.org/10.1002/jrsm.1554

34.

Jung

Seo

Han

Henderson

M. D.

Patall

E. A.

(2020). Prosocial modeling: A meta-analytic review and synthesis. Psychological Bulletin, 146(8), 635–663. https://doi.org/10.1037/bul0000235

35.

Kalaian

H. A.

Raudenbush

S. W.

(1996). A multivariate mixed linear model for meta-analysis. Psychological Methods, 1(3), 227–235. https://doi.org/10.1037/1082-989X.1.3.227

36.

Koenka

A. C.

Linnenbrink-Garcia

Moshontz

Atkinson

K. M.

Sanchez

C. E.

Cooper

(2021). A meta-analysis on the impact of grades and comments on academic motivation and achievement: A case for written feedback. Educational Psychology, 41(7), 922–947. https://doi.org/10.1080/01443410.2019.1659939

37.

Konstantopoulos

(2011). Fixed effects and variance components estimation in three-level meta-analysis. Research Synthesis Methods, 2(1), 61–76. https://doi.org/10.1002/jrsm.35

38.

Lajeunesse

M. J.

(2009). Meta-analysis and the comparative phylogenetic method. The American Naturalist, 174(3), 369–381. https://doi.org/10.1086/603628

39.

Lee

Beretvas

S. N.

(2023). Comparing methods for handling missing covariates in meta-regression. Research Synthesis Methods, 14(1), 117–136. https://doi.org/10.1002/jrsm.1585

40.

López-López

J. A.

Page

M. J.

Lipsey

M. W.

Higgins

J. P. T.

(2018). Dealing with effect size multiplicity in systematic reviews and meta-analyses. Research Synthesis Methods, 9(3), 336–351. https://doi.org/10.1002/jrsm.1310

41.

McCaffrey

D. F.

Bell

R. M.

Botts

C. H.

(2001). Generalizations of biased reduced linearization. Proceedings of the Annual Meeting of the American Statistical Association, August 5–9, 2001.

42.

Nakagawa

Santos

E. S. A.

(2012). Methodological issues and advances in biological meta-analysis. Evolutionary Ecology, 26(5), 1253–1274. https://doi.org/10.1007/s10682-012-9555-5

43.

Page

M. J.

McKenzie

J. E.

Chau

Green

S. E.

Forbes

(2015). Methods to select results to include in meta-analyses deserve more consideration in systematic reviews. Journal of Clinical Epidemiology, 68(11), 1282–1291. https://doi.org/10.1016/j.jclinepi.2015.02.009

44.

Pustejovsky

J. E.

Tipton

(2022). Meta-analysis with robust variance estimation: Expanding the range of working models. Prevention Science, 23(3), 425–438. https://doi.org/10.1007/s11121-021-01246-3

45.

Raudenbush

S. W.

Becker

B. J.

Kalaian

H. A.

(1988). Modeling multivariate effect sizes. Psychological Bulletin, 103(1), 111–120. https://doi.org/10.1037/0033-2909.103.1.111

46.

Roberts

B. R.

MacLeod

C. M.

Fernandes

M. A.

(2022). The enactment effect: A systematic review and meta-analysis of behavioral, neuroimaging, and patient studies. Psychological Bulletin, 148(5–6), 397–434. https://doi.org/10.1037/bul0000360

47.

Rosenthal

Rubin

D. B.

(1986). Meta-analytic procedures for combining studies with multiple effect sizes. Psychological Bulletin, 99(3), 400–406. https://doi.org/10.1037/0033-2909.99.3.400

48.

Schauer

J. M.

Lee

Diaz

Pigott

T. D.

(2022). On the bias of complete- and shifting-case meta-regressions with missing covariates. Research Synthesis Methods, 13(4), 489–507. https://doi.org/10.1002/jrsm.1558

49.

Searle

S. R.

Casella

McCulloch

C. E.

(2009). Variance components. John Wiley.

50.

Sera

Armstrong

Blangiardo

Gasparrini

(2019). An extended mixed-effects framework for meta-analysis. Statistics in Medicine, 38(29), 5429–5444. https://doi.org/10.1002/sim.8362

51.

Tipton

(2015). Small sample adjustments for robust variance estimation with meta-regression. Psychological Methods, 20(3), 375–393. https://doi.org/10.1037/met0000011

52.

Tipton

Pustejovsky

J. E.

(2015). Small-sample adjustments for tests of moderators and model fit using robust variance estimation in meta-regression. Journal of Educational and Behavioral Statistics, 40(6), 604–634. https://doi.org/10.3102/1076998615606099

53.

Tipton

Pustejovsky

J. E.

Ahmadi

(2019). Current practices in meta-regression in psychology, education, and medicine. Research Synthesis Methods, 10(2), 180–194. https://doi.org/10.1002/jrsm.1339

54.

Van den Noortgate

López-López

J. A.

Marn-Martnez

Sánchez-Meca

(2013). Three-level meta-analysis of dependent effect sizes. Behavior Research Methods, 45(2), 576–594. https://doi.org/10.3758/s13428-012-0261-6

55.

Van den Noortgate

López-López

J. A.

Marn-Martnez

Sánchez-Meca

(2015). Meta-analysis of multiple outcomes: A multilevel approach. Behavior Research Methods, 47(4), 1274–1294. https://doi.org/10.3758/s13428-014-0527-2

56.

van Houwelingen

H. C.

Arends

L. R.

Stijnen

(2002). Advanced methods in meta-analysis: Multivariate approach and meta-regression. Statistics in Medicine, 21(4), 589–624. https://doi.org/10.1002/sim.1040

57.

Varin

(2008). On composite marginal likelihoods. AStA Advances in Statistical Analysis, 92(1), 1–28. https://doi.org/10.1007/s10182-008-0060-7

58.

Varin

Reid

Firth

(2011). An overview of composite likelihood methods. Statistica Sinica, 21(1), 5–42. https://www-jstor-org-443.web.bisu.edu.cn/stable/24309261

59.

Viechtbauer

(2010). Conducting meta-analyses in R with the metafor package. Journal of Statistical Software, 36(3), 1–48. https://doi.org/10.18637/jss.v036.i03

60.

Viechtbauer

López-López

J. A.

(2022). Location-scale models for meta-analysis. Research Synthesis Methods, 13(6), 697–715. https://doi.org/10.1002/jrsm.1562

61.

Wei

Higgins

J. P.

(2013). Estimating within-study covariances in multivariate meta-analysis with multiple outcomes. Statistics in Medicine, 32(7), 1191–1205. https://doi.org/10.1002/sim.5679

62.

White

I. R.

Barrett

J. K.

Jackson

Higgins

J. P. T.

(2012). Consistency and inconsistency in network meta-analysis: Model estimation using multivariate meta-regression. Research Synthesis Methods, 3(2), 111–125. https://doi.org/10.1002/jrsm.1045

Supplementary Material

Please find the following supplemental material available below.

For Open Access articles published under a Creative Commons License, all supplemental material carries the same license as the article it is associated with.

For non-Open Access articles published, all supplemental material carries a non-exclusive license, and permission requests for re-use of supplemental material or any part of supplemental material shall be sent directly to the copyright owner as specified in the copyright notice associated with the article.

0.00 MB

0.68 MB