What is already known?

• • Meta-analyzing correlation matrices allows meta-analysts to test theories and hypotheses on the relations between multiple variables.

• • Pooling correlation matrices meta-analytically across studies often relies on univariate random-effects models or multivariate random-effects models (MV-REMs).

• • Hierarchical effect size multiplicity can arise when meta-analysts extract multiple correlation matrices from a study, often due to the inclusion of multiple samples, study sites, or countries.

• • Existing pooling approaches either largely ignore hierarchical effect size multiplicity or only partially account for it.

What is new?

• • Multilevel extensions of MV-REMs provide a framework for addressing hierarchical effect size multiplicity in the meta-analytic pooling of correlation matrices.

• • Dependencies among correlations within correlation matrices can be specified in the asymptotic sampling covariance matrices.

• • Multilevel, multivariate, and random-effects models (MLMV-REMs) enable meta-analysts to model the heterogeneity in correlation matrices within and between studies.

• • These models can be estimated using the multilevel and multivariate features in the R package “metafor.”

Potential impact for RSM readers

• • Meta-analysts who wish to synthesize correlation matrices in the presence of hierarchical effect size multiplicity may take an MLMV-REM approach.

• • We aim to encourage the development and evaluation of models addressing effect size multiplicity in the meta-analytic pooling of correlation matrices.

• • We offer an illustrative example along with the corresponding analytical R code that meta-analysts can adapt to their research.

2.1 Dependencies among effect sizes

While pooling correlation matrices across studies has gained popularity in meta-analysis, it still faces significant challenges. One of these challenges refers to effect size multiplicity that occurs due to the dependencies among effect sizes both within and across studies.Reference Stolwijk, Jak, Eichelsheim and Hoeve ^{10 ,} Reference Wilson, Polanin and Lipsey ^{12 ,} Reference Sheng, Kong, Cortina and Hou ¹³ CheungReference Cheung ³ summarized effect size dependencies as (a) dependencies among the sampling errors that occur due to the administration of multiple measures to the same sample in each study, (b) dependencies among true effect sizes at the population level across studies, and (c) dependencies due to hierarchies in the data, such as multiple correlation coefficients that are nested in studies.

In the context of pooling correlation matrices meta-analytically, such dependencies can occur not only because multiple correlation coefficients are stored in correlation matrices, but also because multiple correlation matrices are extracted from one study. For example, in a meta-analysis of psychological traits—such as anxiety, depression, and stress—studies might report multiple correlation matrices for samples defined by different clinical statuses (e.g., clinical vs. non-clinical samples). In another meta-analysis, authors may extract correlation matrices describing the relations among student achievement, motivation, and engagement in reading, mathematics, and science, using data from international large-scale assessments in education, such as the Programme for International Student Achievement or the Trends in Mathematics and Science Study.Reference Campos, Cheung and Scherer ^{7 ,} Reference Scherer, Siddiq and Nilsen ¹⁴ Given the inclusion of several countries in these assessments, multiple correlation matrices are available per assessment cycle.

These examples illustrate that dependencies among multiple correlation coefficients may not only exist within a correlation matrix, because multiple correlations were derived from the same sample, but also across multiple constructs or measures (i.e., correlational effect size multiplicity).Reference Pustejovsky and Tipton ⁵ The inclusion of multiple correlation matrices per study that are derived from independent samples also creates hierarchical effect size multiplicity.Reference López-López, Page, Lipsey and Higgins ¹⁵

Failing to account for such dependencies can bias meta-analytic results, lead to inaccurate pooled estimates, and compromise statistical inferences.Reference Pustejovsky and Tipton ^{5 ,} Reference López-López, Page, Lipsey and Higgins ^{15 ,} Reference Cheung ¹⁶ Several statistical approaches have been proposed to address effect size multiplicity when pooling correlation matrices meta-analytically, which we describe briefly in the next sections. Notably, we exclude approaches that ignore effect size multiplicity—be it the fact that multiple correlation coefficients are stored in correlation matrices or multiple correlation matrices are available per study. For an overview of these approaches, we kindly refer readers to Stolwijk et al.Reference Stolwijk, Jak, Eichelsheim and Hoeve ¹⁰

2.2 Univariate approaches

Univariate approaches pool each correlation coefficient within a correlation matrix separately across studies via a series of standard, univariate, fixed- or random-effects models.Reference Becker, Aloe, Cheung, Schmid, Stijnen and White ^{1 ,} Reference Cheung ^{3 ,} Reference Borenstein, Hedges, Higgins and Rothstein ¹⁷ For instance, taking a univariate approach, the correlation between anxiety and depression would be meta-analyzed independently of the correlations between anxiety and stress and between stress and depression. In situations where multiple correlations between the same constructs or measures are available due to the inclusion of multiple samples, univariate approaches can account for this hierarchical effect size multiplicity via multilevel random-effects models. Such models estimate the variation of correlation coefficients within and between studies and allow meta-analysts to explore possible causes of this heterogeneity.Reference Van den Noortgate, López-López, Marín-Martínez and Sánchez-Meca ^{18 ,} Reference Harrer, Cuijpers, Furukawa and Ebert ¹⁹ However, this approach neglects the inherent multivariate structure of correlation matrices, where multiple correlations share common constructs and are derived from the same sample (i.e., correlational effect size multiplicity). Ignoring this dependence among correlations within matrices leads to an underestimation of standard errors and, consequently, inflated Type I error rates.Reference Moeyaert, Ugille, Beretvas, Ferron, Bunuan and Van Den Noortgate ^{20 ,} Reference Riley ²¹ Hence, univariate approaches are not recommended for synthesizing correlation matrices with or without hierarchical effect size multiplicity.Reference Jak ⁴

2.3 Multivariate approaches

Multivariate approaches can account for dependencies among correlation coefficients within studies, which arise from measuring multiple outcomes on the same sample.Reference Becker, Aloe, Cooper, Hedges and Valentine ^{22 ,} Reference Becker ²³ These models incorporate not only the sampling variances but also some information about the covariances of correlation coefficients, allowing for the estimation of a pooled correlation matrix that reflects the interrelationships between multiple outcomes.Reference Jak ⁴ By simultaneously leveraging all available correlations, multivariate approaches “borrow strength,” for instance, from the within-study correlations, leading to more precise estimates of mean effect sizes and more reliable assessments of between-study heterogeneity.Reference Riley ^{21 ,} Reference Jackson, White, Price, Copas and Riley ^{24 ,} Reference Jackson, White, Riley, Schmid, Stijnen and White ²⁵ Moreover, multivariate approaches often result in joint confidence intervals and distributions and allow meta-analysts to predict an effect controlling for another, correlated effect.Reference Jackson, White, Riley, Schmid, Stijnen and White ²⁵ However, the application of multivariate approaches in practice is often constrained by the limited availability of information about sampling covariances or correlations in primary studies, and estimation problems may occur due to the increased number of parameters and their joint distributions as compared to separate univariate meta-analyses.Reference Jackson, White, Riley, Schmid, Stijnen and White ^{25 ,} Reference Ahn, Ames and Myers ²⁶ In the following section, we briefly describe some multivariate approaches to pooling correlation matrices meta-analytically.

3.1 Specification of the MLMV-REM

Suppose we have a meta-analytic dataset that comprises correlation coefficients as effect sizes and their respective sampling variances.Footnote ⁱ Let ${\widehat{r}}_{ijk}$ denote the estimate of the ith correlation coefficient in the jth correlation matrix (or, respectively, the jth sample) extracted from the kth primary study. Given this terminology, the jth correlation matrix within the kth primary study that contains I correlation coefficients can then be written as ${\widehat{\mathbf{r}}}_{jk}={\left[{\widehat{r}}_{ijk}\right]}_{i=1}^I$ .

A generic version of an MLMV-REM can be specified by distinguishing between three levels of analysis. At Level 1, the effect size estimates ${\widehat{\mathbf{r}}}_{jk}$ are decomposed into true effect size ${\boldsymbol{\unicode{x3bb}}}_{jk}$ and sampling errorReference Harrer, Cuijpers, Furukawa and Ebert ¹⁹ :

$$\begin{align*}\mathrm{Level}\;1\left(\mathrm{sampling}\kern0.17em \mathrm{error}\right)\kern-1.5pt{:}\ {\widehat{\mathbf{R}}}_{jk}={\boldsymbol{\unicode{x3bb}}}_{jk}+{\mathbf{e}}_{jk}\;\mathrm{with}\;{\mathbf{e}}_{jk}\sim MVN\left(\mathbf{0},{\mathbf{V}}_{jk}\right)\!.\end{align*}$$

Given the multivariate nature of the data (i.e., multiple correlations within a correlation matrix due to multiple outcomes or constructs of interest), the sampling variances ${v}_{ijk}$ are stored in the diagonal of the sampling covariance matrices ${\mathbf{V}}_{jk}$ , and possible covariances among the sampling correlation coefficients are stored off the diagonal. ${\mathbf{V}}_{jk}$ take the following symmetric form:

$$\begin{align*}\mathrm{Sampling}\kern0.17em \mathrm{covariance}\kern0.17em \mathrm{matrix}{:}\ {\mathbf{V}}_{jk}=\left[\begin{array}{cccc}{v}_{1 jk}& & & \\ {}{\mathit{\operatorname{cov}}}_{\left(2,1\right) jk}& {v}_{2 jk}& & \\ {}\vdots & \vdots & \ddots & \\ {}{\mathit{\operatorname{cov}}}_{\left(I,1\right) jk}& {\mathit{\operatorname{cov}}}_{\left(I,2\right) jk}& \cdots & {v}_{Ijk}\end{array}\right]\end{align*}$$

To access the sampling covariances ${\mathit{\operatorname{cov}}}_{\left(i,{i}^{\prime}\right) jk}$ for the $i$ th and ${i}^{\prime }$ th variables, several solutions have been proposed.Footnote ⁱⁱ Wei and HigginsReference Wei and Higgins ³³ described a strategy in which the available information about the correlations among effect sizes is meta-analyzed, and the resultant, pooled correlation coefficient is then used to impute the correlations for studies or matrices that do not provide this information. Another solution is to specify values of the correlations among effect sizes, either informed by existing research or estimated correlations among constructs from the primary study data and then derive the resultant sampling covariances. However, oftentimes, this approach requires restrictive assumptions on correlations between effect sizes, such as the assumption of a constant sampling correlation.Reference Pustejovsky and Tipton ⁵ CheungReference Cheung ³ suggested to translate the meta-analytic models that pool the correlation matrices into the framework of SEM, which makes available an asymptotic sampling covariance matrix. Olkin and FinnReference Olkin and Finn ³⁴ reviewed the explicit expressions for constructing the covariances between two correlation coefficients that were proposed by Olkin and SiotaniReference Olkin and Siotani ³⁵ , who derived them based on large-sample theory. Specifically, skipping the abovementioned indices $j$ and $k$ for simplicity, for variables $A$ , $B$ , $C$ , $D$ , and $N$ representing the sample size, the asymptotic covariance between the correlation coefficients ${r}_{AB}$ and ${r}_{CD}$ isReference Olkin and Finn ³⁴

$$\begin{align*}\mathit{\operatorname{cov}}\left({r}_{AB},{r}_{CD}\right)&=\frac{1}{N-1}\bigg[\frac{1}{2}{\rho}_{AB}{\rho}_{CD}\left({\rho}_{AC}^2+{\rho}_{AD}^2+{\rho}_{BC}^2+{\rho}_{BD}^2\right)+\left({\rho}_{AC}{\rho}_{BD}+{\rho}_{AD}{\rho}_{BC}\right)\\&\quad-\left({\rho}_{AB}{\rho}_{AC}{\rho}_{AD}+{\rho}_{AB}{\rho}_{BC}{\rho}_{BD}+{\rho}_{AC}{\rho}_{BC}{\rho}_{CD}+{\rho}_{AD}{\rho}_{BD}{\rho}_{CD}\right)\bigg].\end{align*}$$

Conversely, the asymptotic covariance between ${r}_{AB}$ and ${r}_{AC}$ —two correlation coefficients sharing variable $A$ —isReference Olkin and Finn ³⁴

$$\begin{align*}\mathit{\operatorname{cov}}\left({r}_{AB},{r}_{AC}\right)=\frac{1}{N-1}\left[\left({\rho}_{BC}-\frac{1}{2}{\rho}_{AB}{\rho}_{AC}\right)\left(1-{\rho}_{AB}^2-{\rho}_{AC}^2\right)+\frac{1}{2}{\rho}_{AB}{\rho}_{AC}{\rho}_{BC}^2\right].\end{align*}$$

Finally, the asymptotic sampling variance of a correlation coefficient ${r}_{AB}$ reduces toReference Borenstein, Hedges, Higgins and Rothstein ¹⁷

$$\begin{align*}{v}_{AB}={\left(1-{\rho}_{AB}^2\right)}^2/\left(N-1\right).\end{align*}$$

Similar equations have also been derived for Fisher- ${Z}_r$ transformed correlations coefficients.Reference Card ^{27 ,} Reference Pustejovsky ³⁶ Irrespective of the choice for an approach to approximate or construct the sampling covariance matrix, the Level 1 specification of the MLMV-REM is the first instance in which meta-analysts must decide on how to represent the effect size multiplicity within correlation matrices.

To further address the nesting of multiple correlation matrices (or samples) in primary studies, two additional levels of analysis are specified. Level 2 contains information about the within-study heterogeneity of correlation coefficients, and Level 3 contains information about the between-study heterogeneity of correlation coefficients.Reference Fernández-Castilla, Jamshidi, Declercq, Beretvas, Onghena and Van den Noortgate ^{8 ,} Reference Van den Noortgate, López-López, Marín-Martínez and Sánchez-Meca ^{18 ,} Reference Harrer, Cuijpers, Furukawa and Ebert ¹⁹ The decomposition of effect sizes and the respective assumptions can be specified as follows:

$$\begin{align*}\mathrm{Level}\;2\;\left(\mathrm{within}\kern0.17em \mathrm{studies}\right)\kern-1.5pt{:}\ {\boldsymbol{\unicode{x3bb}}}_{jk}={\boldsymbol{\unicode{x3b7}}}_k+{\mathbf{u}}_{jk}\mathrm{with}\;{\mathbf{u}}_{jk}\sim MVN\left(\boldsymbol{0},\boldsymbol{\Sigma} \right)\!,\end{align*}$$

$$\begin{align*}\mathrm{Level}\;3\;\left(\mathrm{between}\kern0.17em \mathrm{studies}\right)\kern-1.5pt{:}\ {\boldsymbol{\unicode{x3b7}}}_k=\boldsymbol{\unicode{x3bc}} +{\mathbf{g}}_k\;\mathrm{with}\;{\mathbf{g}}_k\sim MVN\left(\boldsymbol{0},\boldsymbol{\Omega} \right)\!.\end{align*}$$

In this specification, the vector of correlation coefficients ${\boldsymbol{\unicode{x3bb}}}_{jk}$ is decomposed into an average study-level vector of correlations ${\boldsymbol{\unicode{x3b7}}}_k$ and its matrix- or sample-specific deviations from it, ${\mathbf{u}}_{jk}$ . At Level 2, these deviations may exhibit within-study heterogeneity, which is stored in the respective heterogeneity matrix $\boldsymbol{\Sigma}$ . Accordingly, ${\boldsymbol{\unicode{x3b7}}}_k$ is then further decomposed into an overall effect size estimate $\boldsymbol{\unicode{x3bc}}$ and the study-specific deviations from it, ${\mathbf{g}}_k$ . At Level 3, these deviations may exhibit between-study heterogeneity, which is stored in the respective heterogeneity matrix $\boldsymbol{\Omega}$ . Combining all levels of analysis results in the following overall model specification:

$$\begin{align*}{\widehat{\mathbf{R}}}_{jk}=\boldsymbol{\unicode{x3bc}} +{\mathbf{g}}_k+{\mathbf{u}}_{jk}+{\mathbf{e}}_{jk}\;\mathrm{with}\kern0.24em {\mathbf{e}}_{jk}\sim MVN\left(\mathbf{0},{\mathbf{V}}_{jk}\right),{\mathbf{u}}_{jk}\sim MVN\left(\mathbf{0},\boldsymbol{\Sigma} \right),{\mathbf{g}}_k\sim MVN\left(\mathbf{0},\boldsymbol{\Omega} \right)\!.\end{align*}$$

This MLMV-REM allows meta-analysts to quantify correlation-specific heterogeneity within and between primary studies (see Model 1 in Table 1). In this context, “correlation-specific” means that the parameters are estimated for each of the correlations within a correlation matrix. For instance, if meta-analysts wish to synthesize the correlations among four variables, then each correlation matrix contains six correlations off the diagonal, which represent the associations among the six pairs of variables. In our specification of the MLMV-REM, each of these six types of correlations has their own estimates of the within- and between-study heterogeneity variances. In the next section, we describe a range of structures that the heterogeneity matrices $\boldsymbol{\Sigma}$ and $\boldsymbol{\Omega}$ may have.

Table 1 Specification of meta-analytic working models with the MLMV-REM approach

Note. The variable “Cell” is defined as a factor indicating the different types of correlations in the correlation matrix (e.g., “I2–I4” and “I11–I12”). “Correlation-specific” means that the parameters are estimated for each of the 10 correlations. “Constrained” means that the parameter is constrained to be equal across the types of correlations. “Correlation” contains the extracted Pearson correlation coefficients. “V” represents the sampling covariance matrix. “Cell” represents a nominal variable of the names of the different types of correlations. The random-effects correlations $\rho$ and $\phi$ refer to the correlations among random effects within and, respectively, between studies. These correlations were fixed to zero here. “HCS” specifies a heteroscedastic compound symmetric (CS) structure of the random effects, while “CS” specifies a CS structure.Reference Viechtbauer ²⁸

The MLMV-REM approach can be further extended to a mixed-effects meta-regression model by adding explanatory variables, such as sample, study, publication, or measurement characteristics. In such an extension, ${\mathbf{u}}_{jk}$ and ${\mathbf{g}}_k$ may become residuals in a regression sense, depending on the level at which the explanatory variables are defined.Reference Snijders and Bosker ³⁷ Moreover, heterogeneity could be explained specifically for each correlation coefficient at two levels of analysis (i.e., within and between studies). This entails that multiple measures of variance explanations are accessible for each correlation coefficient, once explanatory variables are added to the MLMV-REM (i.e., ${R}_{(2)i}^2$ and ${R}_{(3)i}^2$ for each of the $i=1,\dots, I$ ) correlations).Reference Cheung ^{3 ,} Reference Rights and Sterba ³⁸

Finally, the MLMV-REM approach has the ability to pool correlation matrices with partially missing correlations using likelihood estimation. In this approach, missing correlation coefficients are typically omitted from the likelihood calculation, while all available correlations for each study are retained.Reference McShane and Böckenholt ^{39 ,} Reference Lu ⁴⁰ Thus, rather than imputing unreported values or excluding entire studies, each unobserved correlation simply contributes no data to the joint likelihood, and only the observed correlations are modeled. This effectively results in pairwise deletion at the level of individual correlation coefficients. However, unlike naive pairwise methods that estimate each correlation separately, a joint multivariate model “borrows strength” across the various outcome measures and shared study-level effects, thereby leveraging the observed correlations from each study to improve the estimation of the weighted average correlations.Reference Jackson, White, Price, Copas and Riley ^{24 ,} Reference Riley, Thompson and Abrams ⁴¹

3.2 Specification of the Level 2 and Level 3 random effects

As noted earlier, heterogeneity can be modeled with the MLMV-REM both within and between studies for each type of correlation. Specifically, the symmetric heterogeneity matrices $\boldsymbol{\Sigma}$ and $\boldsymbol{\Omega}$ store the correlation-specific heterogeneity variances in the diagonal and the heterogeneity covariances in the off-diagonal. For I correlation coefficients, these matrices can take the following, often referred to as “unstructured” forms:

$$\begin{align*}\mathrm{Level}-2\;\mathrm{heterogeneity}\kern0.17em \mathrm{covariance}\kern0.17em \mathrm{matrix}{:}\ \boldsymbol{\Sigma} =\left[\begin{array}{cccc}{\tau}_1^2& & & \\ {}{\rho}_{\left(2,1\right)}{\tau}_2{\tau}_1& {\tau}_2^2& & \\ {}\vdots & \vdots & \ddots & \\ {}{\rho}_{\left(I,1\right)}{\tau}_I{\tau}_1& {\rho}_{\left(I,2\right)}{\tau}_I{\tau}_2& \cdots & {\tau}_I^2\end{array}\right],\end{align*}$$

$$\begin{align*}\mathrm{Level}-3\;\mathrm{heterogeneity}\kern0.17em \mathrm{covariance}\kern0.17em \mathrm{matrix}{:}\ \boldsymbol{\Omega} =\left[\begin{array}{cccc}{\gamma}_1^2& & & \\ {}{\phi}_{\left(2,1\right)}{\gamma}_2{\gamma}_1& {\gamma}_2^2& & \\ {}\vdots & \vdots & \ddots & \\ {}{\phi}_{\left(I,1\right)}{\gamma}_I{\gamma}_1& {\phi}_{\left(I,2\right)}{\gamma}_I{\gamma}_2& \cdots & {\gamma}_I^2\end{array}\right].\end{align*}$$

These random-effects covariance matrices contain the within-study variances ${\tau}_i^2$ , the between-study variances ${\gamma}_i^2$ for each correlation, and the covariances that are composed of the within- or between-study standard deviations, ${\tau}_i$ and ${\gamma}_i$ , and the correlations among the within- or between-level random effects of the correlation coefficients, ${\rho}_{\left(i,{i}^{\prime}\right)}$ and ${\phi}_{\left(i,{i}^{\prime}\right)}$ , for $i,{i}^{\prime}=1,\dots, I$ .

This form of $\boldsymbol{\Sigma}$ and $\boldsymbol{\Omega}$ is “unstructured” and allows meta-analysts to include the covariances among random effects that are specific to each pair $\left(i,{i}^{\prime}\right)$ of correlation coefficients. However, although the specification of this form has a large degree of flexibility, meta-analysts may not be able to estimate all parameters within the heterogeneity covariance matrices. For instance, the likely small number of available correlations and primary studies could lead to unstable variances and small statistical power of the respective estimates. Moreover, information about the correlations among the random effects within and between studies may be lacking. Hence, meta-analysts often must simplify the structures of $\boldsymbol{\Sigma}$ and $\boldsymbol{\Omega}$ . We would like to highlight two approaches to simplifying the heterogeneity matrices. The first approach concerns the correlations among the random effects, ${\rho}_{\left(i,{i}^{\prime}\right)}$ and ${\phi}_{\left(i,{i}^{\prime}\right)}$ . Like the specification of the sampling covariance matrices ${\mathbf{V}}_{jk}$ , meta-analysts could fix ${\rho}_{\left(i,{i}^{\prime}\right)}$ and ${\phi}_{\left(i,{i}^{\prime}\right)}$ to some nonzero constant values, $\rho$ and $\phi$ , that are the same for each random effect or approximate them in some way. The structure with constant random-effects correlations is referred to as the “heteroscedastic compound symmetry (HCS)” structure. If $\rho$ or $\phi$ are fixed to zero, then the matrices $\boldsymbol{\Sigma}$ and $\boldsymbol{\Omega}$ have a diagonal structure—a special case of the HCS structure in which meta-analysts can model only the within- and between-study heterogeneity variances and assume that random effects at Levels 2 and 3 are independent. This diagonal structure is often used in the first stage of MASEM, in which multiple correlation matrices are pooled via an MV-REM.Reference Cheung ³ The heterogeneity matrices with an HCS or diagonal structure are as follows:

$$\begin{align*}\mathrm{HCS}\;\mathrm{structures}{:}\ \boldsymbol{\Sigma} =\left[\begin{array}{cccc}{\tau}_1^2& & & \\ {}{\rho \tau}_2{\tau}_1& {\tau}_2^2& & \\ {}\vdots & \vdots & \ddots & \\ {}{\rho \tau}_I{\tau}_1& {\rho \tau}_I{\tau}_2& \cdots & {\tau}_I^2\end{array}\right],\boldsymbol{\Omega} =\left[\begin{array}{cccc}{\gamma}_1^2& & & \\ {}\phi {\gamma}_2{\gamma}_1& {\gamma}_2^2& & \\ {}\vdots & \vdots & \ddots & \\ {}\phi {\gamma}_I{\gamma}_1& \phi {\gamma}_I{\gamma}_2& \cdots & {\gamma}_I^2\end{array}\right],\end{align*}$$

$$\begin{align*}\mathrm{Diagonal}\kern0.17em \mathrm{structures}{:}\ \boldsymbol{\Sigma} =\left[\begin{array}{cccc}{\tau}_1^2& & & \\ {}0& {\tau}_2^2& & \\ {}\vdots & \vdots & \ddots & \\ {}0& 0& \cdots & {\tau}_I^2\end{array}\right],\boldsymbol{\Omega} =\left[\begin{array}{cccc}{\gamma}_1^2& & & \\ {}0& {\gamma}_2^2& & \\ {}\vdots & \vdots & \ddots & \\ {}0& 0& \cdots & {\gamma}_I^2\end{array}\right]\;.\end{align*}$$

Notably, these structures contain heterogeneity variances that are specific to the random effects at Levels 2 and 3 and thus to the correlation coefficients. In other words, for each correlation coefficient (i.e., effect size), within- and between-study variance components are estimated.

A second approach to simplifying $\boldsymbol{\Sigma}$ and $\boldsymbol{\Omega}$ is to constrain these variance components to be equal across random effects. Imposing this constraint to an HCS structure results in a so-called “CS” structure. In the special case where $\rho$ or $\phi$ are fixed to zero, the structure is that of a “scaled identity” structure.Reference Viechtbauer ²⁸ Notably, the term “scaled identity structure” is based on the fact that the two heterogeneity matrices have the forms $\boldsymbol{\Sigma} ={\tau}^2{\boldsymbol{I}}_I$ and $\boldsymbol{\Omega} ={\gamma}^2{\boldsymbol{I}}_I$ , where ${\boldsymbol{I}}_I$ represents an $I\times I$ identity matrix:

$$\begin{align*}\mathrm{CS}\;\mathrm{structures}{:}\ \boldsymbol{\Sigma} =\left[\begin{array}{cccc}{\tau}^2& & & \\ {}{\rho \tau}^2& {\tau}^2& & \\ {}\vdots & \vdots & \ddots & \\ {}{\rho \tau}^2& {\rho \tau}^2& \cdots & {\tau}^2\end{array}\right],\boldsymbol{\Omega} =\left[\begin{array}{cccc}{\gamma}^2& & & \\ {}\phi {\gamma}^2& {\gamma}^2& & \\ {}\vdots & \vdots & \ddots & \\ {}\phi {\gamma}^2& \phi {\gamma}^2& \cdots & {\gamma}^2\end{array}\right],\end{align*}$$

$$\begin{align*}\mathrm{Scaled}\kern0.17em \mathrm{identity}\kern0.17em \mathrm{structures}{:}\ \boldsymbol{\Sigma} =\left[\begin{array}{cccc}{\tau}^2& & & \\ {}0& {\tau}^2& & \\ {}\vdots & \vdots & \ddots & \\ {}0& 0& \cdots & {\tau}^2\end{array}\right],\boldsymbol{\Omega} =\left[\begin{array}{cccc}{\gamma}^2& & & \\ {}0& {\gamma}^2& & \\ {}\vdots & \vdots & \ddots & \\ {}0& 0& \cdots & {\gamma}^2\end{array}\right].\end{align*}$$

If meta-analysts choose one of the structures that assume fixed correlations among random effects, $\rho$ or $\phi$ , in our view, they should examine the extent to which the meta-analytic findings and inferences are sensitive to the choice of these fixed values.Reference Hedges, Tipton and Johnson ^{9 ,} Reference Harrer, Cuijpers, Furukawa and Ebert ¹⁹

3.3 Decision scheme and working models

3.3.1 Analytical decisions

To address hierarchical effect size multiplicity when pooling correlation matrices, meta-analysts must make several analytic decisions. Figure 1 provides an overview of these decisions and their possible outcomes. First, when meta-analysts have extracted the correlation matrices and the respective correlation coefficients from the primary study data, they must decide on the structure of the sampling covariance matrices ${\mathbf{V}}_{jk}$ . Specifically, assuming that correlations within correlation matrices/samples are independent results in diagonal sampling covariance matrices that essentially only contain the sampling variances of the correlation coefficients. Alternatively, meta-analysts could assume some dependence among the sampling errors and then specify ${\mathbf{V}}_{jk}$ with sampling covariances off the diagonals. In scenarios where multiple measures of the same sample have been taken over time or multiple, correlated constructs were assessed for each sample, a dependency structure with some correlation among sampling errors is conceptually reasonable and preferable over an independence structure.Reference Cheung ^{3 ,} Reference Pustejovsky and Tipton ^{5 ,} Reference Harrer, Cuijpers, Furukawa and Ebert ¹⁹ The extracted correlation coefficients ${\widehat{\mathbf{R}}}_{jk}$ and sampling covariance matrices can then be submitted to meta-analytic modeling.

Figure 1 Decision scheme for selecting a working model.

Second, if hierarchical effect size multiplicity is present in the meta-analytic data, meta-analysts must decide if they wish to account for it. In case they wish to ignore effect size multiplicity—for instance, due to close-to-zero within-study heterogeneity in the correlation coefficients—they may choose an MV-REM to pool the correlation matrices.Reference Becker, Aloe, Cheung, Schmid, Stijnen and White ^{1 ,} Reference Stolwijk, Jak, Eichelsheim and Hoeve ^{10 ,} Reference Sheng, Kong, Cortina and Hou ¹³ To model hierarchical effect size multiplicity, meta-analysts may choose a working model within the MLMV-REM approach (see the following section).

Third, meta-analysts must decide on the structures of the within- and between-study heterogeneity matrices. Specifically, do they allow for correlation-specific within- and between-study variances in correlation coefficients, and which assumption do they make on the random-effects correlations within and between studies (i.e., independence assumption with $\rho =\phi =0$ or dependence assumptions with some nonzero values of $\rho$ and $\phi$ )? While the former decision determines which of the four working models is selected, the latter decision must be made for any model that has been selected.

4.1 Description of the meta-analytic data

For illustrative purposes, we chose the data focusing on participants’ self-reported difficulties with describing feelings. This dimension was measured by the following five items of the TAS-20 (see Schroeders et al.Reference Schroeders, Kubera and Gnambs ⁴³ ):

• • Item 2 (I2): It is difficult for me to find the right words for my feelings.

• • Item 4 (I4): I am able to describe my feelings easily. (reverse-coded)

• • Item 11 (I11): I find it hard to describe how I feel about people.

• • Item 12 (I12): People tell me to describe my feelings more.

• • Item 17 (I17): It is difficult for me to reveal my innermost feeling, even to close friends.

Given the focus on these five items, the correlation matrices we intended to pool meta-analytically were 5  $\times$  5 matrices with 10 correlation coefficients each. The data were accessed from a csv file, stored in an R object named tas20, and then supplemented by an effect size identifier ESID.

Output:

To facilitate the meta-analytic modeling and to indicate the different analytical levels, we created identifiers of the effect sizes (ESID), the samples (SampleID), the primary studies (StudyID), and the type of correlation (Cell; e.g., “I2–I4” indicates the correlation between Items 2 and 4). The correlation coefficients were stored in the variable Correlation, and the names of the corresponding variables are indicated by Var1 and Var2. The dataset contains 880 correlation coefficients that were extracted from 88 samples in 62 primary studies, totaling 69,722 participants. On average, the samples comprised 792 participants (M = 792.3, SD = 1,508.8, Mdn = 327, Min = 99, Max = 1,2706), and each study included 1.4 correlation matrices (SD = 0.6, Mdn = 1, Min = 1, Max = 4). Overall, 22 primary studies contained more than two correlation matrices (i.e., samples).

4.2 Constructing sampling covariance matrices

If meta-analysts assume independent sampling errors, then they may directly use the sampling variances of the correlation coefficients as input for the meta-analytic working model. If meta-analysts decide to account for the dependencies among sampling errors within correlation matrices for a specific sample, they may construct sampling covariance matrices using the explicit expressions proposed by Olkin and Siotani.Reference Olkin and Siotani ³⁵ There are several ways to generate the sampling matrices. For instance, the R package “metafor” contains the rcalc() function, which creates a list of sampling covariance matrices (V) for each sample-specific correlation matrix.

This function takes the variable containing the correlation coefficients (Correlation), the names of the respective variables (Var1, Var2), the clustering variable (SampleID), the sample sizes (ni = N), and the name of the dataset (data = tas20) as arguments. The first part specifies that multiple correlations are nested in sample-specific correlation matrices, because each independent sample provided one correlation matrix, and the last part specifies that sampling variances and covariances are generated using the raw correlation coefficients. The resultant object V contains the full sampling covariance matrix with a block design to indicate the correlation matrices.

## Sampling covariance matrix of the first sample

blsplit(V, tas20$SampleID)$`1`

Output:

For the R code using Fisher’s- ${Z}_r$ transformation (via the additional option rtoz = TRUE), please see Supplementary Material S3.

Another option to construct sampling covariance matrices is to use the asyCov() function in the R package “metaSEM.”Reference Cheung ⁴⁴ This function has similar utilities to the above-described function. Nonetheless, it contains some additional features, such as the possibility to also construct sampling correlation matrices (via the cor.analysis option).

V <- asyCov(cmatnew, dat$n, cor.analysis = FALSE)

Several coding steps are leading up to using this function, such as generating a list of correlation matrices from the study samples (cmatnew), along with their sample sizes (n). Please find the details of these steps in Supplementary Material S2. Notice that, from version 1.2.6 of the “metaSEM” package, the asyCov() function no longer implements Cheung and Chan’sReference Cheung and Chan ⁴⁵ SEM approach that rests on a multivariate normality assumption but the equations by Olkin and Finn.Reference Olkin and Finn ³⁴ In the following analyses, we used the V matrix generated from the rcalc() function.

4.3 Approaches to pooling correlation matrices

Once the correlation matrices and sampling (co-)variance matrices have been extracted or generated, meta-analysts can employ the working models we proposed to estimate a pooled correlation matrix and heterogeneity (co-)variances. In this section, we describe in detail the implementation of Working Model 1 and explain briefly the implementation of Models 2–4 in R.

4.3.1 Working Model 1

Given the multivariate nature of the data (i.e., multiple correlations per sample that are based on multiple, measured variables in the primary studies) and the hierarchical effect size multiplicity (i.e., multiple correlation matrices nested in primary studies), we used the rma.mv() function in “metafor” to estimate Model 1. Figure 2 shows an example code with explanations of its elements. Specifically, this function contains the specification of the correlation coefficients (yi), the sampling covariance matrices (V), and the dataset in which these elements are stored (data). Next, the random effects within (|ESID) and between studies (|StudyID) are specified as a list. Given that Model 1 assumes correlation-specific estimates of these random effects, the option ~factor(Cell) is added to this list. Note that the sampling covariances were generated in the previous step using SampleID as a clustering identifier to account for the nesting of multiple correlation coefficients in the samples’ correlation matrices. In the pooling step, the sample specificity of the correlation coefficients is represented by ESID (i.e., each sample contributes only one correlation coefficient within a correlation matrix), and the hierarchical nesting of multiple correlation matrices in primary studies is represented by StudyID.

Figure 2 Elements of the analytic code to specify Working Model 1 in the R package “metafor”.

Suppose that meta-analysts do have some evidence that the Level 2 and Level 3 random effects can be considered independent. Specifying an HCS structure within and between studies (struc = c(“HCS,”"HCS”)), they may assume a correlation of zero within and zero between studies (rho = 0 and, respectively, phi = 0). Next, the estimation procedure is set; in our example, we chose restricted maximum likelihood estimation (method = “REML”). Finally, we specified the variable indicating the type of correlation as a moderator and turned off the intercept of this meta-regression part (mods = ~factor(Cell)-1). This way, we get estimates of each pooled correlation coefficient directly. Labelling the rma.mv() object “tas20.mlmvrem1” and using the summary() function, we obtained the following output of this model estimation:

First, the output contains information about the log-likelihood, deviance, and information criteria, which can be used for further model comparisons. Second, the output displays the within-study variance estimates (tau^2.1-tau^2.10) and the between-study variance estimates (gamma^2.1-gamma^2.10) for each of the 10 correlation coefficients. These estimates exhibited a considerable range, with some variances being close to zero and others as high as 0.0088 (tau^2.7). Third, the Q statistics are shown, indicating, for instance, that statistically significant (residual) heterogeneity existed in the data ( ${Q}_E$ [870] = 6,986.58, p < .01). Fourth, the estimates of the weighted average correlation coefficients, their standard errors, z-statistics, p-values, and the respective confidence intervals are shown (section Model results). Notably, all pooled coefficients were positive and were statistically significantly different from zero.

Table 2 Weighted average effect sizes and variance components in Working Models 1–4 with independent random effects (with $\rho = \phi=0$ )

Note: $\overline{r}$  = Weighted average Pearson correlation coefficient with cluster-robust standard error and 95% confidence interval. “I” refers to the item in the TAS-20. Correlations are labeled using the item numbers (e.g., “I11–I12” refers to the correlation between Items 11 and 12). AIC = Akaike’s Information Criterion; BIC = Bayesian Information Criterion; AICc = Corrected AIC.Reference Viechtbauer ²⁸ All models assume independent random effects (i.e., $\rho =\phi =0$ ) and block-diagonal sampling covariance matrices.

If meta-analysts wish to further obtain cluster-robust standard errors of the pooled correlation coefficient estimates (with primary studies as clusters, indicated by StudyID), they may use the convenience function robust() in the R package “metafor” that calls the “clubSandwich” package:

Selected output:

Table 2 summarizes the key estimates of Model 1. However, the previous output does not contain any information about the preference of Model 1 over alternative models. Hence, in the next section, we describe the estimation of Models 2–4 and perform model comparisons to select an appropriate working model for the meta-analytic pooling of the correlation matrices.

4.3.2 Working Models 2–4

The specification of Working Models 2–4 is like that of Model 1. While Model 1 allowed for correlation-specific estimates of the heterogeneity variances, Models 2–4 restricted these estimates to be the same across correlations. This is reflected in the specification of the structures of Level 2 and Level 3 random effects. Table 1 shows sample codes to estimate Models 2–4 in “metafor.” For instance, in Model 2, the random-effects structures are specified as struc = c(“HCS,”"CS”), and the respective output contains the following variance estimates:

Notably, the within-study variance estimates are correlation-specific (tau^2.1-tau^2.10), and only one heterogeneity variance is estimated at the between level (gamma^2). Similarly, Model 3 can be specified using the option struc = c(“CS,”"HCS”) to implement the assumption of correlation-specific heterogeneity variances between studies but one uniform variance estimate within studies. Finally, the structure option in Model 4 specifies two compound structures struc = c(“CS,”"CS”). Please find the output of the estimation of Models 3 and 4 in Supplementary Material S2. Table 2 shows the parameters of Models 2–4.

4.3.3 Model selection

To select a working model, we compared the information criteria and performed likelihood-ratio tests across Models 1–4. Table 2 contains the respective values of the AIC, BIC, and a corrected AIC. Consistently across these information criteria, Model 2 showed the lowest values and thus the best fit. Moreover, a likelihood-ratio test between Models 1 and 2 suggested that these models did not differ statistically significantly in their fit, ${\unicode{x3c7}}^2$ (9) = 7.60, p = .57. We performed this test using the anova() function:

Output:

Further likelihood-ratio tests suggested the preference of Model 2 over Models 3 and 4. This preference of Model 2 indicated that the between-level variance estimates did not differ significantly and can thus be considered statistically equal across correlation coefficients. Hence, a meta-analyst may select Model 2 as her working model to pool the correlation matrices meta-analytically and conduct subsequent moderator analyses. Notice that our specifications of Models 1–4 assumed that $\rho =\phi =0$ .

4.4 Moderator analyses

As mentioned earlier, Models 1–4 can be extended to mixed-effects meta-regression models to explore moderator effects and examine the extent to which heterogeneity can be explained by the characteristics of the primary studies, samples, or measures. In this section, we illustrate this possible extension by a categorical (i.e., type of the sample, coded as 0 = non-clinical and 1 = clinical; ClinicalBinary) and a continuous moderator (i.e., the proportion of female participants in the primary study; PropFemale).

4.4.1 Categorical moderator

There are several ways of exploring the possible moderating effect of a categorical variable using the MLMVREM approach. One way is to extend the model specification of Model 2 in “metafor” by specifying a crossed effect between the moderator and the variable that indicates the type of correlation (i.e., factor(Cell)). This specification still contains the direct effect of factor(Cell), which serves as the intercepts, that is, the expected weighted average correlations for each type of correlation when the sample is nonclinical. Hence, the model specification no longer contains an overall intercept (−1).

Selected output of the model summary:

The estimates of the model parameters indicate statistically significant and negative moderator effects of the type of the sample on three correlation coefficients: I2–I4 (B = −0.084, SE = 0.032, p = .010), I4–I11 (B = −0.075, SE = 0.024, p = .002), and I4–I17 (B = −0.064, SE = 0.027, p = .017). Hence, these correlations tended to be smaller for clinical samples than for nonclinical samples. The cluster-robust standard errors supported this result. Next, we computed the variance explanations for each level of analysis and some selected correlations.

About 6.4% of the between-study variation in effect sizes was explained by the type of sample, and about 5.6% (6.7%, 0.5%) of the within-study variation in the correlations I2–I4 (I4–11, I4–I17) was explained.

Another way of exploring the effect sizes across the moderator’s categories is to estimate the subgroup-specific effect sizes directly. To achieve this, the moderator option in the above-described specification is changed into

mods =~ factor(Cell):factor(ClinicalBinary)-1.

Consequently, meta-analysts obtain the effect size estimates of each correlation for studies with nonclinical (code 0) and clinical samples (code 1). However, this way, the similarities and differences in effects are not tested statistically.

Selected output:

4.5 Sensitivity analyses

As noted earlier, some decisions on model parameters in the MLMV-REM may be informed by prior research or meta-analysts’ substantive expertise. To examine the sensitivity of the meta-analytic findings, we first explored how different a-priori fixed values of $\rho$ (0.0, 0.5, 1.0) and $\phi$ (0.0, 0.5, 1.0) affected (a) model selection and (b) key parameter estimates of Model 1, such as the weighted average correlation coefficients or heterogeneity variances. Second, we conducted the analyses using Fisher’s- ${Z}_r$ transformed scores.

For the raw correlations, Models 1 and 2 were suitable meta-analytic working models to pool the correlation matrices across the conditions. When $\rho$ and $\phi$ were small, Model 2 was likely to be preferred over Model 1. The relation between $\rho$ (or, respectively, $\phi$ ) and the within-sample variance estimates was mostly V-shaped, and the largest values of the respective variance estimates occurred for $\rho =\phi =1$ . The between-sample variance estimates showed the same pattern. Finally, pooled correlation coefficients increased with higher values of $\rho$ and $\phi$ . For Fisher’s- ${Z}_r$ transformed scores and given values of $\rho$ and $\phi$ , pooled correlation coefficients, within- and between-sample variance estimates were slightly larger than those obtained from the raw correlations.

The detailed results of these sensitivity analyses are shown in Supplementary Materials S2 and S3. Overall, the pooled correlations and heterogeneity estimates were sensitive to the choices of $\rho$ and $\phi$ , and simulation studies are needed to examine this sensitivity systematically.