3.1.1 Fixed-effect model

The fixed-effect or common-effect model usually assumes all studies share the same true or population effect size $\mu $ (cf. the studiesReference Bonett ^{24 ,} Reference Bonett and Price ²⁵ for an alternative fixed-effects model that does not assume equal true effect sizes). The fixed-effect model is

(3.1) $$ \begin{align} \begin{aligned} y_i & = \mu + e_i, \\ \mathsf{Var} \left( e_i \right) &= v_i. \end{aligned} \end{align} $$

Figure 1(a) displays a graphical representation of the fixed-effect model. To analyze a meta-analysis as a structural equation model, we can consider the k studies in a meta-analysis as the k subjects in SEM. Since there are 61 studies in Jaramillo et al.,Reference Jaramillo, Mulki and Marshall ²² we can consider this as having 61 subjects from the SEM perspective. In this perspective, $y_i$ is an observed variable, while a latent variable represents the sampling error $e_i$ . Because the sampling variance $v_i$ is treated as known and varies across studies, we use a definition variable to assign each observed $v_i$ to the variance of the ith study, which is represented by data.vi in OpenMx. The model-implied variance and mean of the ith study are

(3.2) $$ \begin{align} \begin{aligned} \mu_i(\theta) & = \mu, \\ \Sigma_i(\theta) &= v_i. \end{aligned} \end{align} $$

This model has only one parameter ( $\theta = \mu $ ). After fitting the model with FIML, we may obtain the parameter estimate, its SE, and construct a 95% Wald CI by $\hat {\theta } \pm 1.96SE$ .

The typical analysis workflow involves specifying a model using lavaan syntax. Next, the proposed model is converted to the RAM specification using the lavaan2RAM function, and the data are fitted using the sem function. Optionally, the graphical model can be displayed using the plot function, and the model-implied matrices can be obtained using the impliedS function. The following R code specifies a fixed-effect model. The first model, m1a (lines 15–19), explicitly specifies $e_i$ as a latent variable. The mean of $y_i$ is labeled mu in the model (line 16), whereas the known variance $v_i$ is assigned as the sampling variance of $y_i$ using data.vi (line 18). Since the variances are assumed to be present on the observed variables in SEM, the second model, m1b (lines 26–28), fits the same model without specifying $e_i$ . Because these two models are equivalent, it is up to researchers to decide how to present the model.

After specifying the model, we may convert the lavaan equations to RAM specification with the lavaan2RAM function (line 32) and compute the model-implied variance and mean (Equation (3.1)) symbolically with the impliedS function (line 38). We specify std.lv=FALSE (line 32) when converting the model to RAM specification with the lavaan2RAM function because the default is std.lv=TRUE, which indicates that the variances of the latent variables are standardized at 1. However, we want to assign the known sampling variance data.vi in our model (lines 18 and 27) rather than fixing it at 1. Then, we may fit the model using the sem function (lines 41 and 42). Moreover, we may also plot the conceptual model (line 23) and the fitted model (line 46) with the parameter estimates using the plot function. After running the analysis, the parameter estimate and its 95% Wald CIs are $\hat {\mu }=0.1938, [0.1782, 0.2093]$ . We will discuss and compare the results after introducing all models.

## Load the required libraries
library(metaSEM)
library(symSEM)

## Prepare a sample data
## xi: a predictor
## zi: an outcome variable
dat1 <- data.frame(yi=Jaramillo05$r,
                         vi=Jaramillo05$r_v,
                         xi=scale(Jaramillo05$IDV, scale=FALSE),
                         zi=scale(Jaramillo05$Citations, scale=TRUE))

## Fixed-effect model with a latent variable ei
## data.vi means that vi in the data are fixed at each case in the model.
m1a <- "ei =~ 1*yi                 ## e is a latent variable and yi is its indicator
                      yi ~ mu*1                  ## Mean(yi) = mu
                      yi ~~ 0*yi          ## Set the measurement error variance of yi at 0
                      ei ~~ data.vi*ei    ## Var(yi) = vi, which is assumed known in meta-analysis
"

## Display the conceptual model using the semPlot package
## sizeInt=7: Set 7 as the size of the intercept, see help(semPaths) in R
plot(m1a, color = "yellow", sizeInt = 7)

## Fixed-effect model without specifying ei
m1b <- "yi ~ mu*1                   ## Mean(yi) = mu
                      yi ~~ data.vi*yi     ## Var(yi) = vi, the measurement error variance is assigned from 
                                          the data
"

## Convert the lavaan model to RAM specification
## obs.variables = "yi": "yi" is the observed variable
ram1 <- lavaan2RAM(m1a, obs.variables = "yi", std.lv = FALSE)

## Display the RAM model
ram1

## Get the model-implied mean and variance structures
impliedS(ram1)

## Fit the model to the data
fit1 <- sem("Fixed", RAM = ram1, data = dat1)
summary(fit1)

## Display the model with results using the semPlot package
## nDigits=4: Display 4 digits
plot(fit1, color = "green", sizeInt = 7, nDigits = 4)

3.1.2 Random-effects model

Studies can differ in various characteristics, such as samples, research designs, and measurement techniques. These differences can lead to variability in the effect sizes observed across different studies beyond sampling error. The random-effects model addresses this issue by allowing the possibility that each study estimates its own unique population effect size. Specifically, the random-effects model assumes that the observed effect size in each study ( $y_i$ ) is influenced by two sources of variability: true heterogeneity ( $\tau ^2$ ), which represents the actual differences in effect sizes across studies, and sampling error ( $v_i$ ). This approach allows researchers to generalize the findings beyond the individual studies in the meta-analysis, making the results more applicable to a broader context.Reference Hedges and Vevea ²⁶ The random-effects model of the ith study is written as

(3.3) $$ \begin{align} \begin{aligned} y_i & = \mu + u_i + e_i, \\ \mathsf{Var} \left(e_i \right) &= v_i, \\ \mathsf{Var} \left(u_i \right) &= \tau^2, \end{aligned} \end{align} $$

where $\mu $ is the average effect, and $u_i$ is the study-specific true effect deviated from the average effect. The model-implied variance and mean of the ith study are

(3.4) $$ \begin{align} \begin{aligned} \mu_i(\boldsymbol{\theta}) & = \mu, \\ \Sigma_i(\boldsymbol{\theta}) &= \tau^2 + v_i. \end{aligned} \end{align} $$

Now, there are two parameters in the model ( $\boldsymbol {\theta }=[\mu , \tau ^2]^{\mathrm {T}{}}$ ). The random-effects model becomes the fixed-effect model when $\tau ^2=0$ . Several estimators can be used to fit the model and make statistical inferences, such as the unweighted and weighted methods of moments, ML, and restricted maximum likelihood (REML) estimation. Readers can refer to, for example, ViechtbauerReference Viechtbauer ²⁷ and Veroniki et al.,Reference Veroniki, Jackson and Bender ²⁸ for comparing these estimators. The current tutorial uses FIML estimation, though REML estimation can also be implemented in SEM.Reference Cheung ²⁹

Figure 1(b) displays the random-effects model. The study-specific true effect deviating from the average effect is conceptualized as a latent variable $u_i$ . Its variance $\tau ^2$ represents an additive heterogeneity variance (lines 3 and 4 in the following R code). Therefore, there are two latent variables $u_i$ and $e_i$ in model m2a. Figure 1(c) shows the same model by simplifying it with a single latent variable $u_i$ , whereas $e_i$ is conceptualized as the measurement error. The model is shown in m2b in R. To save space, we would not repeat the steps similar to those in the fixed-effect model. Readers can refer to the Supplementary Material for the complete R code and output. After fitting the model, the estimated parameters and their 95% Wald CIs are $\hat {\mu }=0.1866\ [0.1487, 0.2245]$ and $\hat {\tau }^2=0.0170 [0.0089, 0.0251]$ .

## ei is conceptualized as a latent variable.
## data.vi means that vi in the data are fixed at each case in the model.
m2a <- "ui =~ 1*yi                 ## ui is the random effect
                      ui ~~ tau2*ui       ## Var(ui) = tau2
                      ei =~ 1*yi                 ## ei, the error term, is conceptualized as a latent variable
                      ei ~~ data.vi*ei    ## Var(ei) = vi, which is assigned from data
                      yi ~ mu*1                  ## Mean(yi) = mu
                      yi ~~ 0*yi          ## Set the measurement error variance of yi at 0
"

## ei is conceptualized as the measurement error of yi
## data.vi means that vi in the data are fixed at each case in the model.
m2b <- "ui =~ 1*yi                 ## ui is the random effect
                      ui ~~ tau2*ui       ## Var(ui) = tau2
                      yi ~ mu*1                  ## Mean(yi) = mu
                      yi ~~ data.vi*yi    ## Var(yi) = vi, the measurement error variance is assigned from the data
"

3.1.3 Multiplicative-error model

As defined in Equation (3.3), the random-effects model conceptualizes heterogeneity additively through the $u_i$ . However, several alternative models have been proposed to account for heterogeneity in different ways. Notably, Thompson and SharpReference Thompson and Sharp ³⁰ introduced a multiplicative heterogeneity model to meta-analysis. The key idea behind this model is that the sampling error variance $v_i$ is not treated as a fixed and known quantity but rather as something that is estimated up to a proportion. This approach acknowledges the potential uncertainty in the precision of the effect size estimates across studies, offering a different perspective on handling heterogeneity in meta-analysis. Several recent studies, for example,Reference Stanley and Doucouliagos ^{31 –} Reference Stanley, Ioannidis, Maier, Doucouliagos, Otte and Bartoš ³³ suggest this model performs well, especially in handling potential publication bias. The model is

(3.5) $$ \begin{align} \begin{aligned} y_i & = \mu + \sqrt{\phi} e_i, \\ \mathsf{Var} \left( e_i \right) &= v_i, \end{aligned} \end{align} $$

where $\mu $ is the common effect and $\phi $ is a scaling or dispersion factor. The model-implied variance and mean are

(3.6) $$ \begin{align} \begin{aligned} \mu_i(\boldsymbol{\theta}) & = \mu,\\ \Sigma_i(\boldsymbol{\theta}) &= \phi v_i. \end{aligned} \end{align} $$

Figure 2(a) displays the graphical model. The factor loading $\sqrt {\phi }$ scales the impact of the latent factor $e_i$ on $y_i$ . The multiplicative-error model simplifies to the fixed-effect model when $\phi =1$ . However, when the studies exhibit heterogeneity, the estimated $\hat {\phi }$ is typically greater than 1, a situation known as the over-dispersion model.Reference Kulinskaya and Olkin ³⁴ In this context, the estimated SEs are larger than those in the fixed-effect model, reflecting the additional variability due to heterogeneity. It is also possible for the estimated $\hat {\phi }$ to be less than 1, which is referred to as the under-dispersion model. In this case, the estimated SEs are smaller than those in the fixed-effect model, indicating the estimated variability among the studies is even smaller than that assumed by the fixed-effect model. Some researchers opt to truncate $\hat {\phi }$ to 1 to ensure that the estimated SEs of the multiplicative-error model are not smaller than those in the fixed-effect model. This truncation helps maintain a conservative estimate of variability and prevents underestimating SEs. Readers may refer to van Aert and Jackson,Reference van Aert and Jackson ³⁵ which discusses the implications and methodologies regarding the adjustment of $\hat {\phi }$ in the multiplicative-error model.

Figure 2

Structural equation models of meta-analytic models with multiplicative and additive heterogeneity.

Note: Figure (a) illustrates a meta-analytic model characterized by multiplicative heterogeneity, while Figure (b) depicts a model that incorporates both multiplicative and additive heterogeneity. $y_i$ represents the effect size for the ith study, and $data.vi$ indicates the known sampling error variance associated with that effect size. The symbol $\mu $ denotes the average effect across studies. Additionally, $u_i$ represents the additive random effect, with $\tau ^2$ indicating the variance associated with heterogeneity. The parameter $\phi $ serves as the scaling factor for the multiplicative heterogeneity.

Figure 2 Long description

The image consists of two panels, A and B, representing S E M meta-analytic models.

Panel A: Multiplicative Heterogeneity

* At the top, a triangle containing the number 1 has a downward arrow labeled mu pointing to a gray rectangle labeled y sub i.

* To the left, a circle labeled e sub i has a self-loop arrow labeled data dot v i.

* A horizontal arrow labeled square root of phi points from the e sub i circle to the y sub i rectangle.

Panel B: Multiplicative and Additive Heterogeneity

* At the top, a triangle containing the number 1 has a downward arrow labeled mu pointing to a gray rectangle labeled y sub i.

* On the left side, there are two circles stacked vertically.

* The top circle is labeled u sub i and has a self-loop arrow labeled tau squared. An arrow labeled 1 points from u sub i down to the y sub i rectangle.

* The bottom circle is labeled e sub i and has a self-loop arrow labeled data dot v i. An arrow labeled square root of phi points from e sub i up to the y sub i rectangle.

The following R code specifies the multiplicative-error model. The factor loading from the latent factor $e_i$ to $y_i$ is represented by phi_sqrt ( $\sqrt {\phi }$ ) in line 2 of the code. To derive the estimate and CIs on $\phi $ , we introduce a new function of parameters phi, which is defined as the square of phi_sqrt using the := operator (line 6). This transformation allows us to request a Wald CI using the Delta method or an LBCI on phi. After fitting the model, the estimated parameters and their 95% Wald CIs are as follows: the point estimate for $\mu $ is $0.1938$ , with a 95% Wald CI of $[0.1571, 0.2304]$ , and the point estimate for $\sqrt {\phi }$ is $2.3588$ , with a 95% Wald CI of $[1.9402, 2.7773]$ . Consequently, the point estimate for $\phi $ is $5.5637$ , with a 95% Wald CI of $[3.5891, 7.5384]$ . Since the 95% Wald CI for $\phi $ does not include 1, this indicates that the data exhibit over-dispersion with excess heterogeneity, suggesting that the variability in the data is greater than the sampling error.

## Multiplicative error, same as the UWLS model
m3 <- "ei =~ phi_sqrt*yi            ## phi (phi_sqrt^2) is the heterogeneity scaling factor
                     ei ~~ data.vi*ei      ## Var(ei) = vi
                     yi ~ mu*1                    ## Mean(yi) = mu
                     yi ~~ 0*yi            ## Set the default error variance of yi at 0
                     phi := phi_sqrt^2                   ## Define phi as a function of phi_sqrt
"

This multiplicative-error model has been applied in large-scale empirical meta-analyses. For example, Stanley et al.Reference Stanley, Ioannidis, Maier, Doucouliagos, Otte and Bartoš ³³ compared the UWLS and random-effects models across 67,308 Cochrane meta-analyses, finding that the UWLS model often yielded superior representation of the heterogeneity structure in medical research. In addition, Stanley et al.Reference Stanley, Doucouliagos and Havranek ³⁶ presented evidence that UWLS can help reduce bias in meta-analyses of correlations, particularly in situations where the effect sizes and their sampling variances are correlated.

3.1.4 Hybrid models of additive and multiplicative errors

Several researchersReference Nam, Mengersen and Garthwaite ^{37 ,} Reference DuMouchel and Berry ³⁸ have integrated additive and multiplicative error meta-analysis into hybrid models within a Bayesian framework. More recently, Baker and JacksonReference Baker and Jackson ³⁹ proposed a hybrid model by combining two traditions of meta-analysis—the multiplicative model in particle physics and the additive model in the meta-analysis literature. Given the complexity of their approach, it is instructive first to present a simplified version suggested by Schmid.Reference Schmid ⁴⁰

Hybrid model 1: Essentially, the simplest hybrid model combines two sources of variability: study-specific variance $\tau ^2$ and scaled sampling variance $\phi v_i$ Reference Schmid ⁴⁰ :

(3.7) $$ \begin{align} \begin{aligned} y_i & = \mu + u_i + \sqrt{\phi} e_i, \\ \mathsf{Var} \left( e_i \right) &= v_i, \\ \mathsf{Var} \left( u_i \right) &= \tau^2. \end{aligned} \end{align} $$

The model-implied variance and mean of the ith study are

(3.8) $$ \begin{align} \begin{aligned} \mu_i(\boldsymbol{\theta}) & = \mu, \\ \Sigma_i(\boldsymbol{\theta}) &= \tau^2 + \phi v_i. \end{aligned} \end{align} $$

Figure 2(b) displays the hybrid model. Conceptually, the two sources of variability are represented as two latent factors. The square root of the scaling factor $\sqrt {\phi }$ (phi_sqrt in the R code) is depicted by the factor loading of the latent factor $e_i$ to $y_i$ . The hybrid model becomes the random-effects model in Equation (3.3) when $\phi =1$ , whereas it is equivalent to the multiplicative model in Equation (3.5) when $\tau ^2=0$ . The following R code specifies the hybrid model. The latent factors ui (lines 5 and 6) and ei (lines 8 and 9) represent the study-specific effect and the sampling error, respectively. Similar to the R code for the multiplicative error model, we define a new function to compute the original $\phi $ (line 10).

After fitting the model, the parameter estimates and their 95% Wald CIs are as follows: $\hat {\mu } = 0.1874$ with a 95% Wald CI of $[0.1497, 0.2251]$ , $\sqrt {\hat {\phi }} = 1.2870$ with a 95% Wald CI of $[0.3825, 2.1916]$ , and $\hat {\tau }^2 = 0.0137$ with a 95% Wald CI of $[0.0007, 0.0267]$ . Additionally, the estimated value of $\phi $ and its 95% Wald CI is $1.6565$ with a 95% Wald CI of $[-0.6720, 3.9849]$ . Since the Wald CI for $\phi $ (or $\sqrt {\hat {\phi }}$ ) includes 1, it is not statistically significant at the 0.05 level for the hypothesis test $H_0: \phi = 1$ .

## Hybrid model 1: multiplicative and additive error
m4 <- "yi ~ mu*1                       ## Mean(yi) = mu
                     yi ~~ 0*yi               ## Set the measurement error variance of yi at 0
                     ## additive error
                     ui =~ 1*yi
                     ui ~~ tau2*ui            ## Var(ui) = tau2
                     ## multiplicative error
                     ei =~ phi_sqrt*yi               ## phi is a scaling factor
                     ei ~~ data.vi*ei         ## Var(ei) = vi
                     phi := phi_sqrt^2                      ## Define phi as a function of phi_sqrt
"

Hybrid model 2: Now, we discuss the original hybrid model proposed by Baker and Jackson.Reference Baker and Jackson ³⁹ The model is

(3.9) $$ \begin{align} \begin{aligned} y_i & = \mu + \tilde{e}_i, \\ \mathsf{Var} \left( \tilde{e}_i \right) &= v_i + \tau^2 v^{\omega}_i, \end{aligned} \end{align} $$

where $0 < \omega < 1$ . The model-implied variance and mean in the ith study are

(3.10) $$ \begin{align} \begin{aligned} \mu_i(\boldsymbol{\theta}) & = \mu, \\ \Sigma_i(\boldsymbol{\theta}) &= v_i + \tau^2 v^{\omega}_i. \end{aligned} \end{align} $$

This hybrid model simplifies to the random-effects model (Equation (3.3)) when $\omega = 0$ . When $\omega = 1$ , the model-implied variance is $\Sigma _i(\boldsymbol {\theta }) = (1 + \tau ^2)v_i$ , where $(1 + \tau ^2)$ acts as a non-negative scaling factor. This contrasts with the hybrid model in Equation (3.7), where the scaling factor $\phi $ can be smaller than 1. Please also note that $\tau ^2$ is no longer the additive heterogeneity variance because it is multiplied by $v_i$ . For values of $0 < \omega < 1$ , the model exhibits intermediate behavior between additive and multiplicative models.

It is not straightforward to represent this hybrid model graphically. Here, we introduce a method to specify a structural equation model with arbitrary functions of parameters. In lavaan syntax, the == operator ensures a linear constraint on the parameters. For example, a1==a2*a3 enforces that the parameter a1 will equal the product of the parameters a2 and a3 after the optimization. It should be noted that all three parameters (a1, a2, and a3) must have been defined before imposing the constraint. Moreover, there are only two free parameters here, as the third one is constrained by the equation.

The following R code specifies a generic model with mu as the mean and Vary as the variance of $y_i$ (lines 2 and 3). The line Vary==data.vi+tau2*data.vi^omega (line 5) implements $\Sigma _i(\boldsymbol {\theta }) = v_i + \tau ^2 v^{\omega }_i$ . However, the parameters (tau2 and omega) have not yet been defined in the model. To resolve this, we need to add the replace.constraints=TRUE argument (lines 11 and 14) to tell the functions to replace the left-hand side parameter (Vary) with the right-hand side function and automatically create the two parameters (tau2 and omega). With this trick, we may fit models with arbitrary functions on the mean and variance structure (see the next section). After fitting the model, the parameter estimates and their 95% Wald CIs are $\hat {\mu }=0.1869, [0.1491, 0.2248]$ , $\hat {\tau }^2 = 0.0331, [-0.0733, 0.1394]$ , and $\hat {\omega } = 0.1204, [-0.4562, 0.6971]$ . The $\hat {\omega }$ appears closer to 0 (indicating an additive error model) than 1 (indicating a multiplicative error model).

## Hybrid model 2: multiplicative and additive error
m5 <- "yi ~ mu*1                ## Mean(yi) = mu
                     yi ~~ Vary*yi     ## Var(yi) = vary
                     ## Define Vary as a function of vi, tau2, and omega
                     Vary == data.vi + tau2*data.vi^omega
"

ram5 <- lavaan2RAM(m5, obs.variables = "yi", std.lv = FALSE)

## replace.constraints = TRUE: Replace the constraint "Vary" with the new parameters
impliedS(ram5, replace.constraints = TRUE)

## replace.constraints = TRUE: Replace the constraint "Vary" with the new parameters
fit5 <- sem("Hybrid2", RAM = ram5, data = dat1, replace.constraints = TRUE)
summary(fit5)

Hybrid model 3: The third version of the hybrid model is the HKSJ model.Reference Hartung ^{41 –} Reference Sidik and Jonkman ⁴³ These authors recognized that the average effect size’s SE does not perform well under the random-effects model when the number of studies is small, primarily because the heterogeneity variance $\tau ^2$ and sampling variance $v_i$ are estimated rather than known. They proposed adjusting the calculated SE, which has been shown to perform better than the random-effects model when there are only a small number of studies.Reference Röver, Knapp and Friede ^{44 –} Reference IntHout, Ioannidis and Borm ⁴⁶ At first glance, it may not be apparent how their adjustment can be considered a hybrid model of multiplicative and additive error. To understand this, we can compare the multiplicative error model with the fixed-effect model. In the multiplicative error model, the term $\phi v_i$ in Equation (3.5) is proportional to $v_i$ in the fixed-effect model (Equation (3.1)). This proportionality means that the estimated common effects $\mu $ are identical in both models, but they differ in their SEs due to the scaling factor $\phi $ . Conventionally, the multiplicative error models are estimated with UWLS, which does not impose $\phi =1$ .Reference Stanley and Doucouliagos ^{31 ,} Reference Stanley and Doucouliagos ³²

van Aert and JacksonReference van Aert and Jackson ³⁵ noted the similar relationship between the HKSJ and the random-effects models and suggested estimating the HKSJ model with UWLS. Implementing the HKSJ method involves two steps. First, a random-effects model is fitted to the data. This step produces an estimate of the heterogeneity variance, $\tilde {\tau }^2$ , which, combined with the known sampling variance $v_i$ , results in the total error variance $\tilde {v}_i = \tilde {\tau }^2 + v_i$ in the ith study. The HKSJ model is

(3.11) $$ \begin{align} \begin{aligned} y_i & = \mu + \tilde{e}_i, \\ \mathsf{Var} \left( \tilde{e}_i \right) &= \phi(\tilde{\tau}^2 + v_i) = \phi \tilde{v}_i. \end{aligned} \end{align} $$

And the model-implied variance and mean in the ith study are

(3.12) $$ \begin{align} \begin{aligned} \mu_i(\boldsymbol{\theta}) & = \mu, \\ \Sigma_i(\boldsymbol{\theta}) &= \phi(\tilde{\tau}^2 + v_i). \end{aligned} \end{align} $$

The random-effects model assumes that the $y_i$ is distributed as $\tilde {v}_i$ . On the other hand, the HKSJ model considers that $y_i$ is only distributed proportional to $\tilde {v}_i$ , that is, $\phi \tilde {v}_i$ . Thus, the average effect is identical under the random-effects and HKSJ models, whereas the SEs differ unless $\phi =1$ . The HKSJ-like multiplicative adjustment can be implemented as a structural equation model with the following R code. This implementation represents the heterogeneity effect and sampling error as two latent variables, $u_i$ and $e_i$ . Their impact on the observed effect size $y_i$ is scaled by the factor phi_sqrt, as shown in lines 5 and 8 of the R code. The heterogeneity variance, obtained from the output of the random-effects model, is treated as a fixed value (0.01703407 in this example). After fitting the model, the estimated average effect $\hat {\mu }$ and its 95% Wald CI are 0.1866, with an interval of [0.1484, 0.2248]. The estimated scaling factor $\hat {\phi }$ and its 95% Wald CI are 1.0173, with an interval of [0.6563, 1.3784]. It is important to note that this model has no estimate for $\tau ^2$ , as it is treated as a fixed value. Researchers can use the value and confidence interval from the random-effects model output for completeness.

## Hybrid model 3: Hartung-Knapp-Sidik-Jonkman method
m6 <- "yi ~ mu*1                   ## Mean(yi) = mu
                     yi ~~ 0*yi           ## Set the measurement error variance of yi at 0
                     ## additive error
                     ui =~ phi_sqrt*yi           ## phi is a scaling factor
                     ui ~~ 0.01703407*ui  ## Estimated tau2 is obtained from the random-effects model
                     ## multiplicative error
                     ei =~ phi_sqrt*yi           ## phi is a scaling factor
                     ei ~~ data.vi*ei     ## Var(ei) = vi
                     phi := phi_sqrt^2                  ## Define phi as a function of phi_sqrt
"

Table 2 summarizes the results of all six models. The estimate of heterogeneity variance ( $\hat {\tau }^2$ ) and its SE for the hybrid 3 (HKSJ) model were copied from those of the random-effects model. Although the hybrid 3 model has the lowest AIC and BIC values, these values are likely incorrect because they do not account for the estimation of $\tau ^2$ . For instance, these values were calculated based on two parameters (mu and phi_sqrt), while $\tau ^2$ was not considered in the model complexity.

Table 2

Results of six meta-analytic models

Table 2 Long description

The table consists of 7 columns: Name, Fixed, Random, Multiplicative, Hybrid 1, Hybrid 2, and Hybrid 3.

Key parameter estimates and their standard errors (S E) include:

* mu: Values range from 0.1866 (Random and Hybrid 3) to 0.1938 (Fixed and Multiplicative). S E ranges from 0.0079 to 0.0195.

* tau 2: N A for Fixed and Multiplicative; 0.0170 for Random and Hybrid 3; 0.0137 for Hybrid 1; 0.0331 for Hybrid 2. S E ranges from 0.0041 to 0.0543.

* phi_sqrt: N A for Fixed, Random, and Hybrid 2; 2.3588 for Multiplicative; 1.2870 for Hybrid 1; 1.0086 for Hybrid 3. S E ranges from 0.0913 to 0.4615.

* omega: Only present for Hybrid 2 at 0.1204 (S E 0.2942).

Model fit statistics:

* Number of parameters: Fixed (1), Random (2), Multiplicative (2), Hybrid 1 (3), Hybrid 2 (3), Hybrid 3 (2).

* minus 2 L L: Fixed (129.0697), Random (minus 55.4423), Multiplicative (minus 44.6262), Hybrid 1 (minus 55.7740), Hybrid 2 (minus 55.6036), Hybrid 3 (minus 55.4513).

* A I C: Fixed (131.0697), Random (minus 51.4423), Multiplicative (minus 40.6262), Hybrid 1 (minus 49.7740), Hybrid 2 (minus 49.6036), Hybrid 3 (minus 51.4513).

* B I C: Fixed (133.1806), Random (minus 47.2205), Multiplicative (minus 36.4045), Hybrid 1 (minus 43.4414), Hybrid 2 (minus 43.2710), Hybrid 3 (minus 47.2296).

Note: Fixed, Random, Multiplicative, Hybrid1, Hybrid2, and Hybrid3 represent the models in Equations (3.1), (3.3), (3.5), (3.7), (3.9), and (3.11), respectively. Names with (SE) are the standard errors of the parameter estimates. $-$ 2LL = $-$ 2*loglikelihood. AIC = Akaike information criterion. BIC = Bayesian information criterion.

By inspecting the table, we can confirm that the estimated average effect ( $\hat {\mu }$ ) is identical in both the fixed-effect and multiplicative error models. However, the SE of the multiplicative error model is larger than that of the fixed-effect model because $\sqrt {\hat {\phi }}> 1$ . The same pattern is observed when comparing the random-effects and HKSJ models. The SE is larger in the HKSJ model due to the adjustment factor $\phi $ . When we compare the AIC values, it becomes evident that the fixed-effect model is the worst fit for the data, unlike the other models. The random-effects model and the three hybrid models all show similar good fits. However, it is crucial to remember that this is not a definitive assessment of their statistical performance, as there is only one sample. Future research may explore the advantages and limitations of these models.

3.2.1 A fixed-X meta-regression

The most common model is to use a covariate to predict the effect size in cases of excess heterogeneity among studies. This approach is called a mixed-effects meta-analysis or meta-regression. Similar to regression analysis, predictors in this context can be seen as either fixed or random. In a fixed-X approach, predictors are considered fixed values and are part of the design matrix. This approach is widely used in meta-analysis, as seen in works in the literature.Reference Sera, Armstrong, Blangiardo and Gasparrini ^{47 –} Reference Raudenbush, Cooper, Hedges and Valentine ⁴⁹ On the other hand, the SEM literature typically assumes a random-X approach, treating predictors as random variables. Both approaches have their own strengths and limitations. Treating $x_i$ as fixed values has the advantage of not requiring any distributional assumptions about $x_i$ . This makes the method versatile, as it can handle binary, categorical, and continuous covariates equally well. However, the primary limitation of this approach is that it does not accommodate missing values in the covariates; any missing data must be removed before analysis.

To better understand these two approaches, let us consider the fixed-X meta-regression, where the individualism score of the countries where the studies were conducted, $x_i$ , serves as the covariate. The model is

(3.13) $$ \begin{align} \begin{aligned} y_i &= \beta_0 + \beta_1 x_i + u_i + e_i, \\ \mathsf{Var} \left( u_i \right) &= \tau^2, \\ \mathsf{Var} \left( e_i \right) &= v_i. \end{aligned} \end{align} $$

The model-implied variance and mean of the ith study are

(3.14) $$ \begin{align} \begin{aligned} \mu_i(\theta) & = \beta_0 + \beta_1 x_i, \\ \Sigma_i(\theta) &= \tau^2 + v_i. \end{aligned} \end{align} $$

Figure 3(a) presents the graphical model in SEM, which is implemented in the following R code. This model introduces a phantom variable $\eta _i$ , a latent variable with a zero varianceReference Rindskopf ⁵⁰ to represent the covariate. Specifically, the value of $x_i$ in the ith study is assigned to the specific path by the definition variable data.xi (lines 7–9). Consequently, the mean of $y_i$ in the ith study is modeled as $\beta _0 + \beta _1 x_i$ . The variance structure of the model closely resembles that of a random effects model, with the exception that $\tau ^2$ is now the residual heterogeneity variance after accounting for the covariate $x_i$ . Since $x_i$ is treated as a design matrix or specified fixed values, the analysis has only one observed variable $y_i$ , as indicated on line 12.

Figure 3

Structural equation models of univariate meta-analysis with one study characteristic.

Note: Figures (a) and (b) represent fixed-x and random-x mixed-effects meta-analyses, respectively. In these panels, $y_i$ and $data.vi$ represent the effect size and its known sampling error variance in the ith study, respectively, while $u_i$ and $e_i$ denote the random effect and sampling error in the ith study, respectively. Additionally, $\mu $ represents the average effect, and $\tau ^2$ represents the random-effects (residual) heterogeneity variance.

Figure 3 Long description

Panel A, fixed-x model. A circle labeled u sub i with a self-loop tau-squared points to a shaded rectangle y sub i with a weight of 1. Above y sub i, a circle eta sub i points down to y sub i with weight beta sub 1. To the right, a triangle containing 1 points to eta sub i with weight data dot x i and to y sub i with weight beta sub 0. A self-loop below y sub i is labeled data dot v i. Panel B, random-x model. Similar to A, but eta sub i is replaced by a shaded rectangle x i. A self-loop above x i is labeled sigma sub x squared. The triangle points to x i with weight mu sub x and to y sub i with weight beta sub 0. The path from x i to y sub i is weighted beta sub 1.

## Mixed-effects model with the covariate as a design matrix
m7 <- "yi ~ b0*1        ## Intercept of yi: b0
             yi ~~ data.vi*yi  ## Known sampling variance vi
             ui =~ 1*yi               ## ui: random effects
             ui ~~ tau2*ui     ## Residual heterogeneity variance tau2
             ## etai: a phantom latent variable to represent the covariate x
             etai =~ b1*yi            ## yi = b1*etai
             etai ~ data.xi*1         ## etai = xi
             etai ~~ 0*etai    ## Var(etai) = 0
"
## Only yi is an observed variable
ram7 <- lavaan2RAM(m7, obs.variables = "yi", std.lv = FALSE)

3.2.2 A random-X meta-regression

A standard practice in SEM is to treat $x_i$ as another variable in the model. The random-X model of the ith study is shown below:

(3.15) $$ \begin{align} \begin{aligned} y_i &= \beta_0 + \beta_1 x_i + u_i + e_i, \\ x_i &= \mu_x + e_{xi}, \\ \mathsf{Var} \left( u_i \right) &= \tau^2, \\ \mathsf{Var} \left( e_i \right) &= v_i, \\ \mathsf{Var} \left( e_{xi} \right) &= \sigma^2_x \\ \mathsf{E} \left( x_i \right) &= \mu_x. \end{aligned} \end{align} $$

The model-implied covariance matrix and the mean vector are

(3.16)

The main strength of treating $x_i$ as another variable is its ability to handle missing effect sizes and covariates using FIML estimation without excluding any data,Reference Enders ⁵¹ for a discussion. However, a key limitation is that the data are assumed to follow a multivariate normal distribution. Another key distinction between the fixed-X and random-X approaches is how the explained variance $R^2$ is calculated. The fixed-X approach involves fitting two models: an unconditional model that excludes covariates and a conditional model that incorporates the covariates. The $R^2$ is defined byReference Raudenbush, Cooper, Hedges and Valentine ⁴⁹

(3.17) $$ \begin{align} R^2_{\mathrm{fixed\text{-}x}} = \frac{\hat{\tau}^2_0 - \hat{\tau}^2_1}{\hat{\tau}^2_0}, \end{align} $$

where $\hat {\tau }^2_0$ and $\hat {\tau }^2_1$ are the estimated heterogeneity variances without and with the covariates. In contrast, the random-X approach, commonly used in SEM, defines $R^2$ using only the conditional model. The formula is

(3.18) $$ \begin{align} R^2_{\mathrm{random\text{-}x}} = \frac{\hat{\boldsymbol{\beta}}^{\mathrm{T}{}} \boldsymbol{S}_x \hat{\boldsymbol{\beta}} }{\hat{\boldsymbol{\beta}}^{\mathrm{T}{}} \boldsymbol{S}_x \hat{\boldsymbol{\beta}} + \hat{\tau}^2_1}, \end{align} $$

where $\boldsymbol {S}_x$ is the sample covariance matrix of the covariates and $\hat {\boldsymbol {\beta }}$ is the vector of the regression coefficients.

The $R^2_{\mathrm {random\text {-}x}}$ has several advantages over the $R^2_{\mathrm {fixed\text {-}x}}$ . First, the random-X approach retains all available data in the analysis in the presence of missing covariates, whereas the fixed-X approach must exclude studies with missing covariates. Lee and BeretvasReference Lee and Beretvas ⁵² examined several methods for handling missing covariates in the meta-analysis, concluding that the random-X approach, using the FIML estimation, outperformed the listwise deletion approach under their simulation conditions. However, further commentary by WhiteReference White ⁵³ and responses from Lee and BeretvasReference Lee and Beretvas ⁵⁴ highlight additional considerations and debates surrounding this issue. Second, the $R^2_{\mathrm {fixed\text {-}x}}$ becomes negative when $\hat {\tau }^2_1 > \hat {\tau }^2_0$ . Although $R^2_{\mathrm {random\text {-}x}}$ can also be negative, it only occurs when $\hat {\tau }^2_1 < 0$ , which is less likely than $\hat {\tau }^2_1 > \hat {\tau }^2_0$ . Third, CIs for $R^2_{\mathrm {random\text {-}x}}$ can be constructed using both Wald CI and LBCI methods, as it is defined as a function of other estimated parameters within the model. In contrast, constructing CIs for $R^2_{\mathrm {fixed\text {-}x}}$ is less straightforward. One option is to apply a bootstrap method to the data. However, it may violate the assumption that the covariates are fixed values.

Figure 3(b) illustrates the SEM utilizing the random-X approach. In this model, $x_i$ is treated as a variable characterized by its own mean and variance, as indicated in lines 6 and 7 of the following R code. The $R^2_{\mathrm {random\text {-}x}}$ is defined in line 9. It is crucial to note that when converting the equations to the RAM specification, $x_i$ must be included as a variable, as demonstrated in line 12. Upon fitting the models, it is observed that the results of both the fixed-X and random-X approaches are identical, as there is no missing data. The estimated coefficients, $\hat {\beta }_0$ and $\hat {\beta }_1$ , are 0.1860 and $-$ 0.0013, respectively. The 95% CIs for these coefficients are [0.1479, 0.2239] for $\hat {\beta }_0$ and [ $-$ 0.0033, 0.0006] for $\hat {\beta }_1$ . Additionally, the $R^2_{\mathrm {random\text {-}x}}$ , along with its 95% LBCI, is reported as 0.0371 with a confidence interval of [0.0000, 0.1926].

## Mixed-effects model with the covariate as a variable
m8 <- "yi ~ b0*1 + b1*xi ## Intercept (b0) and slope (b1)
                     yi ~~ data.vi*yi  ## Known sampling variance vi
                     ui =~ 1*yi               ## ui: random effects
                     ui ~~ tau2*ui     ## Residual heterogeneity variance tau2
                     xi ~ Meanx*1             ## Mean(xi) = Meanx
                     xi ~~ Varx*xi     ## Var(xi) = Varx
                     ## Define R2: explained variance/total variance
                     R2 := b1^2*Varx/(b1^2*Varx + tau2)
"
## Both yi and xi are observed variables
ram8 <- lavaan2RAM(m8, obs.variables = c("yi", "xi"), std.lv = FALSE)

3.2.3 Regressing a covariate on a true effect size

In traditional meta-analysis, the effect size is typically treated as the outcome variable, while the observed variables are predictors. However, there are instances where these roles may be reversed, with the effect size being used to predict the observed variables. Several researchers, for example,Reference Hanel and Haase ^{55 ,} Reference Lortie, Aarssen, Budden and Leimu ⁵⁶ , have been interested in exploring the relationship between effect size and the future impact of studies, measured in terms of citations. Since effect sizes are reported before the accumulation of citations, it is logical to use the effect size as a predictor in such analyses. The common practice involves conducting a regression analysis by regressing the number of citations on the effect sizes.Reference Hanel and Haase ⁵⁵ A significant issue arises in this context: each $y_i$ is distributed with a known sampling variance $v_i$ , which is usually ignored in the applications. It is well-documented in the literature that ignoring the precision in the predictors can lead to biased estimates of the regression coefficients and SEs.Reference Fuller ⁵⁷ Moreover, the calculated $R^2$ is also biased because it ignores the known sampling variance $v_i$ .

This problem can be effectively addressed using SEM. In SEM, the external variable $z_i$ can be regressed on the true effect size $f_{yi}$ rather than the observed effect size $y_i$ . The following equations depict the model where the number of citations $z_i$ is regressed on the true effect size $f_{yi}$ :

(3.19) $$ \begin{align} \begin{aligned} y_i &= f_{yi} + e_i, \\ z_i &= \beta_0 + \beta_1 f_{yi} + e_{zi}, \\ \mathsf{Var} \left( e_i \right) &= v_i, \\ \mathsf{E} \left( f_{yi} \right) &= \mu_y, \\ \mathsf{Var} \left( f_{yi} \right) &= \tau^2, \\ \mathsf{Var} \left( e_z \right) &= \sigma^2_z. \end{aligned} \end{align} $$

The model-implied variance and mean of the ith study are

(3.20)

Figure 4 illustrates the graphical model where the true effect size is conceptualized as a latent variable $f_{yi}$ . This latent variable is then used to predict the number of citations $z_i$ . The following R code is employed to implement the model. To reduce the need for extensive typing, only the true effect size $f_{yi}$ is explicitly defined as a latent variable in line 2. In contrast, the sampling error $e_i$ is not explicitly specified. Instead, the known sampling variance $v_i$ is assigned to the observed effect size $y_i$ using the syntax yi~~data.vi*yi in line 6. The $R^2$ , which does not involve $v_i$ , is defined on line 10. The results from the model indicate that the parameter regression coefficient $\hat {\beta }_1$ is estimated to be $-$ 1.9755, with the $SE=1.0733$ . The $R^2$ value is estimated to be 0.0672, with a 95% Wald CI of [0.0000, 0.2043] after truncating the negative estimated lower bound to zero. In contrast, the result of the wrong approach by regressing the outcome on the sample effect size is $\hat {\beta }_1=-1.4114$ with the $SE=0.8051$ and an $R^2$ value of 0.0495.

Figure 4

A structural equation model of regressing one study characteristic $z_i$ on the true effect size $f_{yi}$ .

Note: $y_i$ and $data.vi$ represent the effect size and its known sampling error variance in the ith study, respectively. $f_{yi}$ and $z_i$ represent the true effect size and the outcome variable in the ith study, respectively. Furthermore, $\mu _y$ signifies the average effect, while $\tau ^2_y$ indicates the variance of random-effects heterogeneity.

## Regressing a covariate on the true effect size
m9 <- "fyi =~ 1*yi                ## fyi: true effect size
                     fyi ~~ tau2*fyi     ## Var(fyi) = tau2
                     fyi ~ Meany*1              ## Mean(fyi) = Meany
                     yi ~ 0*1
                     yi ~~ data.vi*yi    ## Var(ei) = vi
                     zi ~ b0*1 + b1*fyi  ## Intercept (b0) and slope (b1)
                     zi ~~ ErrVarz*zi    ## Error variance of zi
                     ## Define an R2
                     R2 := b1^2*tau2/(b1^2*tau2 + ErrVarz)
"

3.2.4 Meta-regression with partially standardized covariates

When the predictors’ scales differ, interpreting regression coefficients using raw scores (unstandardized solutions) can be challenging.Reference Schmidt ⁵⁸ For example, our previous meta-regression analysis noted that when the individualism score ( $x_i$ ) of the countries in the studies increases by one unit, the expected effect size $y_i$ decreases by 0.0013. Although this change is not statistically significant, we will continue to use it in our discussion. The standard deviation of $x_i$ is 19.1356, making it difficult to contextualize the value.

Standardized coefficients are typically preferred, especially when predictors are on different scales. Standardization facilitates comparison and interpretation by ensuring that all variables have the same variance, thereby allowing for more accurate comparison and interpretation. There are two types of standardization: complete and partial. Complete standardization involves standardizing the predictors and the outcome variable, while partial standardization only standardizes the predictors. In meta-analysis, partial standardization is preferred because it maintains the scale of the effect size, which is essential for indicating residual heterogeneity.

A naive approach might involve manually standardizing the variables before conducting a meta-regression. However, this method is known to produce biased SEs because it does not consider the uncertainty of the standardization process.Reference Yuan and Chan ^{59 ,} Reference Jones and Waller ⁶⁰ A more accurate approach is to treat the standardized variables as latent variables within the SEM framework.Reference Cheung ⁶¹ This method properly incorporates the uncertainty and provides more accurate SEs.

Figure 5(a) and the subsequent equations depict a model for a partially standardized solution. This model has two latent variables, $f_{yi}$ and $f_{xi}$ , representing the unstandardized true effect size and partially standardized predictor, respectively. The standard deviation of the predictor is represented by its factor loading $\sigma _{xi}$ (sd_xi in line 8), whereas the variance of $f_{xi}$ is fixed at 1 (line 9). The residual variance of the true effect size is represented by $\tau ^2$ (tau2_res in line 3), which is used to define the $R^2$ in line 14. After fitting the partially standardized solution, the estimated $R^2$ is $0.0371$ . The lower bound of the 95% Wald CI is truncated to zero, whereas the upper bound is 0.1425:

(3.21) $$ \begin{align} \begin{aligned} y_i &= f_{yi} + e_i, \\ x_i &= \sigma_x f_{xi}, \\ f_{yi} &= \beta_0 + \beta_1 f_{xi} + u_i, \\ \mathsf{Var} \left( e_i \right) &= v_i, \\ \mathsf{Var} \left( u_i \right) &= \tau^2_{f_y}, \\ \mathsf{Var} \left( f_{xi} \right) &= 1, \\ \mathsf{E} ( f_{yi}) &= \beta_0 + \beta_1 \mu_x/\sigma_{x}, \\ \mathsf{E} \left( x_i \right) &= \mu_x/\sigma_x. \end{aligned} \end{align} $$

The model-implied variance and mean of the ith study are

(3.22)

Figure 5

Structural equation models of partially standardized meta-regression (a) and correlation between a true effect size and a covariate (b).

Note: Figures (a) and (b) illustrate models with partially standardized effects and correlation between an effect size and a covariate. In these panels, $y_i$ and $data.vi$ represent the effect size and its known sampling error variance in the ith study, respectively, while $u_i$ and $e_i$ denote the random effect and sampling error in the ith study, respectively. Additionally, $\mu $ represents the average effect, and $\tau ^2$ represents the random-effects (residual) heterogeneity variance.

Figure 5 Long description

Two structural equation models labeled A and B.

Panel A: A triangle containing the number 1 is on the far left. Two arrows originate from it. The top arrow, labeled mu sub x forward slash sigma sub x, points to a circle labeled f sub x i. The bottom arrow, labeled beta sub 0, points to a circle labeled f sub y i. The circle f sub x i has a self-loop labeled 1 and an arrow labeled beta sub 1 pointing down to f sub y i. From f sub x i, an arrow labeled sigma sub x points right to a gray rectangle x sub i. From f sub y i, an arrow labeled 1 points right to a gray rectangle y sub i. Circle f sub y i has a self-loop labeled tau squared. Rectangle y sub i has a self-loop labeled data dot v i.

Panel B: Similar layout to A. The triangle 1 has two arrows. The top arrow is mu sub x forward slash sigma sub x pointing to f sub x i. The bottom arrow is mu sub y forward slash tau sub y pointing to f sub y i. Instead of a downward arrow, a double-headed arrow labeled rho (f sub y, f sub x) connects f sub x i and f sub y i. Circle f sub x i has a self-loop of 1. Circle f sub y i has a self-loop of 1. From f sub x i, an arrow labeled sigma sub x points to rectangle x sub i. From f sub y i, an arrow labeled tau sub y points to rectangle y sub i. Rectangle y sub i has a self-loop labeled data dot v i.

## Partial standardized xi but not yi
m10 <- "fyi =~ 1*yi                 ## fyi: Latent effect size of yi
                      fyi ~~ tau2_res*fyi  ## tau2_res: residual heterogeneity variance
                      fyi ~ b0*1 + b1*fxi  ## b0: intercept, b1: slope
                      yi ~ 0*1
                      yi ~~ data.vi*yi     ## Var(ei) = vi
                      ## fxi: Standardized score of xi
                      fxi =~ sd_xi*xi             ## sd_xi: sd of xi
                      fxi ~~ 1*fxi         ## Var(fxi) =1 for standardization
                      fxi ~ mu_x_sd*1
                      xi ~ 0*1
                      xi ~~ 0*xi
                      ## Define R2
                      R2 := b1^2/(b1^2 + tau2_res)
"

3.2.5 Correlating a covariate with a true effect size

The correlation coefficient is crucial in statistics. It is the cornerstone of regression analysis, SEM, and many multivariate statistics. It provides a quick understanding of how two variables are linearly related to one another, and it is a simple tool for detecting potential multicollinearity among predictors, which can be problematic in regression models. As a result, it is a common practice to present a correlation table of the variables in academic papers, offering a snapshot of the linear relationships between them.

In the context of meta-analyses, correlations between effect sizes and covariates are sometimes reported. However, researchers often calculate these correlations by correlating the sample effect sizes $y_i$ and the covariates.Reference Hanel and Haase ^{55 ,} Reference Lortie, Aarssen, Budden and Leimu ^{56 ,} Reference Barto and Rillig ⁶² Similarly, researchers sometimes calculate the correlation between two effect sizes by ignoring their sampling error. ⁶³ As discussed in the previous section, this approach can introduce bias into the estimated correlations and their SEs.Reference Schauer, Fitzgerald, Peko-Spicer, Whalen, Zejnullahi and Hedges ⁶⁴ The presence of sampling variances $v_i$ can distort the true relationships between variables, leading to inaccurate and potentially misleading conclusions. It is, therefore, essential to account for these variances to obtain more reliable and valid correlation estimates.

Figure 5(b) and the following equations define the model to measure the correlation between the true effect size $f_{yi}$ and a covariate $x_i$ . The variances of the standardized $f_{yi}$ and latent covariate $f_{xi}$ are fixed at one, where $\tau _{f_y}$ and $\sigma _x$ represent the standard deviations of true effect size and $x_i$ . Then the covariance between $f_{yi}$ and $f_{xi}$ is the required correlation $\rho _{(x,f_y)}$ . It should be noted that the means of $f_{yi}$ and $f_{xi}$ are also scaled by $1/\tau _{f_y}$ and $1/\sigma _x$ . The model is

(3.23) $$ \begin{align} \begin{aligned} y_i &= \tau_{f_y} f_{yi} + e_i, \\ x_i &= \sigma_x f_{xi}, \\ \mathsf{Var} \left( e_i \right) &= v_i, \\ \mathsf{Var} \left( f_{yi} \right) &= \mathsf{Var} \left( f_{xi} \right) = 1, \\ \mathsf{Cor} \left( f_{yi}, f_{xi} \right) &= \rho_{(x,f_y)} \\ \mathsf{E} \left( f_{yi} \right) &= \mu_y/\tau_{f_y}, \\ \mathsf{E} \left( f_{x_i} \right) &= \mu_x/\sigma_x. \end{aligned} \end{align} $$

The model-implied covariance matrix and mean vector of the ith study are

(3.24)

The following R code implements this model. The key idea is to define two standardized latent variables (fyi and fxi) to represent the standardized effects of the true effect and covariate. The corresponding factor loadings (tauy and sd_xi) now represent the standard deviations of the variables (lines 3 and 9). After fitting the model, the estimated correlation and its 95% Wald CI are $\hat {\rho }= -0.1925, [-0.4667, 0.0816]$ , which is not statistically significant at 0.05. If we disregard the known sample variance of the effect size and correlate the sample effect size with the covariate, the estimated correlation is $-$ 0.1734, with a 95% Wald CI of [ $-$ 0.4074, 0.0820].

## Correlation between fyi and fxi
m11 <- "## fyi: Standardized true effect size of yi
                      fyi =~ tauy*yi            ## tauy: sd of fyi
                      fyi ~~ 1*fyi       ## Var(fyi) = 1
                      fyi ~ y_mean_div_tauy*1                       ## rescaled mean
                      yi ~ 0*1
                      yi ~~ data.vi*yi
                      ## fxi: Standardized score of xi
                      fxi =~ sd_xi*xi           ## sd_xi: sd of xi
                      fxi ~~ 1*fxi       ## Var(fxi) = 1 for standardization
                      fxi ~~ Cor*fyi     ## Correlation between fyi and fxi
                      fxi ~ x_mean_div_sigmax*1 ## rescaled mean
                      xi ~ 0*1
                      xi ~~ 0*xi
"

3.2.6 Location-scale and nonlinear models

Most meta-regression models in the literature are linear models. Figure 6(a) illustrates a linear relationship between the predictor $x_i$ and the true effect size $f_{yi}$ . The shaded area represents the 95% predictive intervals (PIs), which approximately include 95% of the true effect sizes. While nonlinear relationships, such as quadratic or cubic terms, can still be accommodated within linear models by including polynomial terms, there are cases where genuinely nonlinear models may be required.Reference Boucher and Bennetts ^{65 –} Reference Breit, Scherrer, Tucker-Drob and Preckel ⁶⁹ This section addresses these advanced models.

Figure 6

Some meta-regression models.

Note: Figure (a) shows a standard mixed-effects model on the true effect size ( $f_{y_i}$ ) with $f_{yi} = \beta _0 + \beta _1 x_i + u_i$ and $\mathsf {Var} \left ( u_i \right )=\tau ^2$ . Figure (b) shows a location-scale model on the heterogeneity variance with $f_{yi} = \mu + u_i$ and $\mathsf {Var} \left ( u_i \right ) = \exp (\alpha _0 + \alpha _1 x_i)$ . Panel (c) shows a location-scale model on both the mean and the heterogeneity variance with $f_{yi} = \beta _0 + \beta _1 x_i + u_i$ and $\mathsf {Var} \left ( u_i \right ) = \exp (\alpha _0 + \alpha _1 x_i)$ . The shaded areas represent the 95% predictive intervals of the true effect size.

Figure 6 Long description

A three-panel figure with a shared grid layout. Each panel has a horizontal x-axis labeled Predictor (x sub i) and a vertical y-axis labeled True effec size (f sub y i).

* Panel A shows a linear increase in the true effect size as the predictor increases. The shaded predictive interval maintains a constant width throughout the trend, representing a standard mixed-effects model.

* Panel B shows a horizontal flat line for the true effect size. The shaded predictive interval is wide on the left and narrows significantly as the predictor increases, representing a location-scale model on the heterogeneity variance.

* Panel C shows a linear increase in the true effect size similar to Panel A. However, the shaded predictive interval is wide at low predictor values and narrows as the predictor increases, representing a location-scale model applied to both the mean and the heterogeneity variance.

The location-scale models have been utilized in multilevel models.Reference Hedeker, Mermelstein and Demirtas ⁷⁰ and SEMReference McNeish ⁷¹ Recently, Viechtbauer and López-LópezReference Viechtbauer and López-López ²³ extended these models for meta-analysis. The core idea is to model the mean and variance structures separately. For example, Figure 6(b) depicts a model with an average effect where the heterogeneity decreases as $x_i$ increases. Conversely, Figure 6(c) illustrates that $x_i$ has a positive relationship with the expected mean but a negative relationship with the residual heterogeneity variances. These models are similar to the regression models with heteroskedasticity.

In the location-scale model with additive error, the model-implied variance and mean for the ith study are as follows:

(3.25) $$ \begin{align} \begin{aligned} \mu_i(\theta) & = \beta_0 + \beta_1 x_i, \\ \Sigma_i(\theta) &= \exp{(\alpha_0+\alpha_1 x_i)} + v_i, \end{aligned} \end{align} $$

where $\tau ^2_0=\exp {(\alpha _0)}$ is the residual heterogeneity variance when $x_i=0$ . An exponential function is typically employed to ensure that the predicted residual heterogeneity variance is nonnegative. The mean structure is identical to that in a standard meta-regression. The intriguing aspect is the variance structure. Although we use $x_i$ in this example’s mean and variance structures, the predictors can be different for these two structures.

The following R code implements the location-scale model in Equation (3.25). Lines 2–6 define a generic random-effects model with the average effect mu and residual heterogeneity variance tau2. Lines 8 and 9 define the linear relationship for mu and the nonlinear relationship for tau2. b0 ( $\beta _0$ ) and b1 ( $\beta _1$ ) are the regression coefficients in predicting the mean structure mu, whereas a0 ( $\alpha _0$ ) and a1 ( $\alpha _1$ ) are the regression coefficients in predicting the variance structure tau2. Since we aim to replace the parameters on the left-hand side with the expressions on the right-hand side, we need to add the replace.constraints = TRUE argument in the analysis (see the discussion in Section 3.1.4). After fitting the model, neither regression coefficient is statistically significant in predicting the mean and variance structures. The parameter estimates and their 95% Wald CIs are ( $\hat {\beta }_1=-0.0013, [-0.0031, 0.0005]$ ) and ( $\hat {\alpha }_1=0.0054, [-0.0204, 0.0312]$ ). The results are identical to those in the metafor package.

## Location-scale model with an additive heterogeneity variance
m12 <- "yi ~ mu*1                 ## Mean(yi) = mu
                      yi ~~ data.vi*yi   ## Var(ei) = vi
                      ## ui is the latent variable of the random effect
                      ui =~ 1*yi                ## yi = 1*ui
                      ui ~~ tau2*ui      ## Var(ui) = tau2 of the residual heterogenenity variance
                      ## Specify a linear relationship on the mean structure
                      mu == b0 + b1*data.xi
                      ## Specify a nonlinear relationship on the variance structure
                      tau2 == exp(a0 + a1*data.xi)
"

Since the location-scale model has been implemented in the metafor package, some readers may wonder about the necessity of implementing their own model given this existing functionality. The present SEM approach offers the flexibility to specify a wide range of arbitrary models tailored to specific research needs. For example, the location-scale model with two or more effect sizes can be easily accommodated within the SEM-based meta-analysis framework. Additionally, we can also explore models that have not yet been implemented in current software. For instance, to my knowledge, the location-scale model has not been extended to incorporate multiplicative error or hybrid models (cf. the studyReference Kulinskaya and Olkin ³⁴ for other ideas). As an example, the following model-implied variance and mean in the ith study shows a location-scale model with a multiplicative error:

(3.26) $$ \begin{align} \begin{aligned} \mu_i(\theta) & = \beta_0 + \beta_1 x_i, \\ \Sigma_i(\theta) &= \exp{(\alpha_0+\alpha_1 x_i)}v_i, \end{aligned} \end{align} $$

where $\phi _0=\exp {(\alpha _0)}$ is the heterogeneity scaling factor when $x_i=0$ . The implementation of the location-scale model with a multiplicative error is shown below. After fitting the model, the predictor is statistically significant in predicting the mean structure ( $\hat {\beta }_1=-0.0026$ , 95% Wald CI = $[-0.0047, -0.0005]$ ) but not for the variance structure ( $\hat {\alpha }_1=-0.0110$ , 95% Wald CI = $[-0.0313, 0.0093]$ ). In theory, we can also create a location-scale model for the hybrid model that incorporates both multiplicative and additive heterogeneity. Still, it is yet to be determined what the interpretations of the parameters are due to the large number of parameters involved.

## Location-scale model with a multiplicative heterogeneity variance
m13 <- "yi ~ mu*1                 ## Mean(yi) = mu
                      yi ~~ sigma2*yi    ## Var(yi) = sigma2
                      ## Specify nonlinear on mean and variance structures
                      mu == b0 + b1*data.xi
                      sigma2 == exp(a0 + a1*data.xi)*data.vi
"