Imputing Covariance for Meta-Analysis in the Presence of Interaction

1. Introduction

When conducting an individual patient data (IPD) meta-analysis, researchers may need to examine interactions between two predictors. [1,2] In a one-stage meta-analysis,[3] all data are combined into a single dataset before analysis, allowing for the calculation of the covariance between predictors. However, this approach may not be feasible in situations where data cannot be consolidated into a single dataset, such as when data-sharing restrictions prevent data transfer. In such cases, a two-stage IPD meta-analysis is necessary.[4]

When interactions are present, it is crucial to calculate confidence intervals for both the combined main effect and the interaction effect.[5] Additionally, determining the overall covariance between the two coefficients is vital. To the best of our knowledge, no existing studies have addressed the calculation of confidence intervals for combined effects in two-stage IPD meta-analyses involving interactions.

To fill this gap, we propose a method for calculating the overall covariance in the presence of interactions between two effects. This method is applicable to both fixed-effect and random-effect meta-analysis models, providing a practical solution for researchers analyzing interactions within a meta-analytic framework.

2. Regression with Interaction

Suppose we aim to conduct a random-effects model meta-analysis to explore the interaction between falls and age and investigate whether the effect of falls on hip fractures varies across different age levels. (Note: This example is for illustrative purposes only and is not intended to provide any evidence on the topic.) A Poisson regression model can be utilized to examine these effects.

The Poisson model is specified as follows:

$$Y_i\sim Poisson\left({\mu }_i\right)$$

(1)

with the following log-linear relationship:

$$\begin{array}{cc}\hfill log\left({\mu }_i\right)& ={\widehat{\beta }}_0+{\widehat{\beta }}_1*fall+{\widehat{\beta }}_2*(fall*age)+{\widehat{\beta }}_3*age\hfill \\ & ={\widehat{\beta }}_0+({\widehat{\beta }}_1+{\widehat{\beta }}_2*age)*fall+{\widehat{\beta }}_3*age\hfill \end{array}$$

The effect of falls can be calculated as:

$$\begin{array}{cc}\hfill \beta & =log\left({\mu }_{fall=1}\right)-log\left({\mu }_{fall=0}\right)\hfill \\ & ={\widehat{\beta }}_0+({\widehat{\beta }}_1+{\widehat{\beta }}_2*age)*1+{\widehat{\beta }}_3*age-[{\widehat{\beta }}_0+({\widehat{\beta }}_1+{\widehat{\beta }}_2*age)*0+{\widehat{\beta }}_3*age]\hfill \\ & =({\widehat{\beta }}_1+{\widehat{\beta }}_2*age)\hfill \end{array}$$

If exponentiatinge $\beta$, yields the rate ratio (RR) for the effect.

The effect of falls depends on $({\widehat{\beta }}_1+{\widehat{\beta }}_2*age)$, and its $(1-\alpha )100\%$ confidence interval can be calculated as:[5]

$$({\widehat{\beta }}_1+{\widehat{\beta }}_2*age)\pm t_{n-p-1,\alpha /2}\sqrt{var\left({\widehat{\beta }}_1\right)+ag{e}^2*var\left({\widehat{\beta }}_2\right)+2*age*cov({\widehat{\beta }}_1,{\widehat{\beta }}_2)}$$

(2)

In this model, n represents the sample size, p denotes the number of predictors, and for large samples, the t score can be approximated by the z score. The variable "age" serves as an effect modifier.

For a two-stage meta-analysis, the overall effects and variances of ${\widehat{\beta }}_1$ and ${\widehat{\beta }}_2$ can be calculated using the inverse-weighted average method.[6,7,8,9] Moreover, it is essential to combine these two regression coefficients to construct confidence intervals for the effect of falls at different age levels.

We present a method for combining coefficients and their covariance in a two-stage random-effects meta-analysis when interactions are present.

Suppose we have data from 10 cohort studies. For each cohort, we fit the specified Poisson regression model and obtain the following results:

In Table 1, ${\widehat{\beta }}_{1i}$ represents the coefficient for falls (i.e., the main effect) in the $i^{\mathrm{th}}$ cohort, while ${\widehat{\beta }}_{2i}$ corresponds to the interaction effect of fall * age. $\mathrm{var}\left({\widehat{\beta }}_{1i}\right)$ and $\mathrm{var}\left({\widehat{\beta }}_{2i}\right)$ denote the variances of ${\widehat{\beta }}_{1i}$ and ${\widehat{\beta }}_{2i}$, respectively, and $\mathrm{cov}({\widehat{\beta }}_{1i},{\widehat{\beta }}_{2i})$ is the covariance between these coefficients. Finally, $n_i$ represents the sample size for each cohort.

3. Meta Analysis for the Main Effect and Interaction Effect

Typically, the inverse-variance weighted method is employed for both fixed-effect and random-effect model meta-analyses.

3.1. Fixed effect meta-analysis

For the fixed-effect meta-analysis, the overall main effect, ${\beta }_{1,\mathrm{overall}}$, can be calculated as:

$${\beta }_{1,\mathrm{overall}}={\displaystyle \frac{{\sum }_{i=1}^k{w}_{1i}{\widehat{\beta }}_{1i}}{{\sum }_{i=1}^k{w}_{1i}}}$$

(3)

where :

$$w_{1i}={\displaystyle \frac{1}{var\left({\widehat{\beta }}_{1i}\right)}}$$

and k is the number of studies.

The overall variance of the main effect, ${\beta }_{1,\mathrm{overall}}$, from the fixed-effect model can be calculated as:

$$V_{1,\mathrm{overall}}={\displaystyle \frac{1}{{\sum }_{i=1}^k{w}_{1i}}}$$

(4)

For the interaction effect, ${\beta }_{2i}$, the calculations of the overall effect, ${\beta }_{2,\mathrm{overall}}$, and the overall variance, $V_{2,\mathrm{overall}}$, are analogous to those of the main effect.

3.2. Random-effect meta-analysis

For the random-effect meta-analysis, the overall main effect, ${\beta }_{1,\mathrm{overall}}$, can be calculated as:

$${\beta }_{1,\mathrm{overall}}={\displaystyle \frac{{\sum }_{i=1}^k{w}_{1i}^{*}{\widehat{\beta }}_{1i}}{{\sum }_{i=1}^k{w}_{1i}^{*}}}$$

(5)

where:

$$w_{1i}^{*}={\displaystyle \frac{1}{var\left({\widehat{\beta }}_{1i}\right)+{\tau }^2}}$$

The overall variance of the main effect, ${\beta }_{1,\mathrm{overall}}$, from the random-effects model can be calculated as:

$$V_{1,\mathrm{overall}}={\displaystyle \frac{1}{{\sum }_{i=1}^k{w}_{1i}^{*}}}$$

(6)

The between-study variance, ${\tau }^2$, for the random-effects meta-analysis is calculated in the following steps:

Step 1: Calculate the Q-statistic

$$Q=\sum _{i=1}^k{w}_i{\beta }_i^2-{\displaystyle \frac{{\left({\sum }_{i=1}^k{w}_i{\beta }_i\right)}^2}{{\sum }_{i=1}^k{w}_i}}$$

(7)

Step 2: Calculate c

$$c=\sum _{i=1}^k{w}_i-{\displaystyle \frac{{\sum }_{i=1}^k{w}_i^2}{{\sum }_{i=1}^k{w}_i}}$$

(8)

Step 3: calculate ${\tau }^2$

$${\tau }^2=\left\{\begin{array}{cc}& {\displaystyle \frac{Q-df}{c}}\mathrm{if}Q>df\hfill \\ & 0\mathrm{if}Q<df\hfill \end{array}\right.$$

(9)

Here, $df$ represents the degrees of freedom, which is equal to the number of studies minus one, i.e. $k-1$.

Similarly, for the interaction effect, ${\beta }_{2i}$, the overall effect, ${\beta }_{2,\mathrm{overall}}$, and its variance, $V_{2,\mathrm{overall}}$, can be calculated using the same approach as for the main effect.

After calculating the Q-statistic, the $I^2$-statistic, which quantifies the percentage of variation across studies due to heterogeneity, can be determined using:

$$I^2=\left\{\begin{array}{cc}& {\displaystyle \frac{Q-df}{Q}}*100\%\mathrm{if}Q>df\hfill \\ & 0\mathrm{if}Q<df\hfill \end{array}\right.$$

(10)

4. Covariance in Meta-Analysis

4.1. Relationship Between Covariance and Correlation Coefficient r

The overall main effect, ${\beta }_{1,\mathrm{overall}}$, and the overall interaction effect, ${\beta }_{2,\mathrm{overall}}$, can be calculated using popular R packages such as ’meta’ and ’metafor’, or SAS macros. Consequently, the overall effect at different ages can be determined. However, to the best of our knowledge, no existing papers provide guidance on how to combine covariance in the presence of interactions in meta-analysis to calculate an overall covariance, as required by Equation(2).

In this context, we propose a method for calculating the overall covariance in meta-analysis. This method leverages the relationship between the correlation coefficient, variance, and covariance, as expressed in the following equation:[10]

$$r={\displaystyle \frac{cov(X,Y)}{{\sigma }_1{\sigma }_2}}$$

(11)

Here, $cov(X,Y)$ represents the covariance between two random variables X and Y, while ${\sigma }_1$ and ${\sigma }_2$ are their standard deviations. Specifically, X and Y can be replaced with the estimated coefficients ${\widehat{\beta }}_1$ and ${\widehat{\beta }}_2$, and the standard deviations with their corresponding standard errors. Based on this, the overall covariance can be calculated as:

$$cov{({\widehat{\beta }}_1,{\widehat{\beta }}_2)}_{overall}=r_{overall}\times \sqrt{V_{{\widehat{\beta }}_{1,\mathrm{overall}}}{V}_{{\widehat{\beta }}_{2,\mathrm{overall}}}}$$

(12)

In Equation (12), $cov{({\widehat{\beta }}_1,\widehat{\beta }2)}_{\mathrm{overall}}$ represents the overall covariance between ${\widehat{\beta }}_1$ and ${\widehat{\beta }}_2$, and $r_{overall}$ is the overall correlation coefficient. This overall covariance is derived from either a fixed-effect or random-effects meta-analysis. However, a key step involves calculating an overall correlation coefficient $r_{\mathrm{overall}}$, which is not straightforward in meta-analysis.

To calculate the overall correlation coefficient, $r_{\mathrm{overall}}$, the following steps can be followed:[6]

1. For each cohort, extract the standard errors of the two coefficients and their covariance from the regression model.

2.Calculate the correlation coefficient $r_i$ for each cohort using the Equation (11)

$$r_i={\displaystyle \frac{cov({\widehat{\beta }}_{1i},{\widehat{\beta }}_{2i})}{{\sigma }_{1i}{\sigma }_{2i}}}$$

(13)

where $cov({\widehat{\beta }}_{1i},{\widehat{\beta }}_{2i})$ represents the covariance between the two coefficients, and ${\sigma }_{1i}$, ${\sigma }_{2i}$ are the standard errors of ${\widehat{\beta }}_{1i}$ and ${\widehat{\beta }}_{2i}$, respectively.

3. Perform a random-effect meta-analysis (or fixed-effect meta-analysis) on $r_i$ values to estimate $r_{\mathrm{overall}}$.

Once $r_{\mathrm{overall}}$ is determined, it can be substituted into Equation(12) to calculate the overall covariance $cov{({\widehat{\beta }}_1,{\widehat{\beta }}_2)}_{overall}$

4.2. Fixed-Effect Meta-Analysis for Correlation Coefficient r

To conduct a meta-analysis for correlation coefficients, Fisher’s z-transformation is first applied, expressed as:

$$z_i=0.5\times ln\left({\displaystyle \frac{1+r_i}{1-r_i}}\right)$$

(14)

The variance of $z_i$ is calculated as:

$$v_i={\displaystyle \frac{1}{n_i-3}}$$

(15)

where $n_i$ represents the valid sample size for each cohort.

The meta-analysis for $z_i$ follows the same procedure as for other statistics, such as, such as ${\beta }_i$. Once the overall z value is obtained,it is transformed back to r using the equation:

$$r={\displaystyle \frac{e^{2z}-1}{e^{2z}+1}}$$

(16)

For the fixed-effect meta-analysis of $z_i$, the weight is given by:

$$w_i={\displaystyle \frac{1}{v_i}}$$

(17)

where $v_i$ is variance of $z_i$ which is ${\displaystyle \frac{1}{n_i-3}}$

The overall Fisher’Z is calculated as:

$$Z_{overall}={\displaystyle \frac{{\sum }_{i=1}^k{w}_i{z}_i}{{\sum }_{i=1}^k{w}_i}}$$

(18)

where k is the total number of studies included in the meta-analysis.

All above calculations are for fixed-effect meta-analysis. Then we can transform $Z_{overall}$ back to $r_{overall}$ using Equation (16).

4.3. Random-Effect Meta-Analysis for Correlation Coefficient r

For a random-effect meta-analysis, the following additional steps are required:

1. Calculate the Q-statistic:

$$Q=\sum _{i=1}^k{w}_i{z}_i^2-{\displaystyle \frac{{\left({\sum }_{i=1}^k{w}_i{z}_i\right)}^2}{{\sum }_{i=1}^k{w}_i}}$$

(19)

2. Calculate the between-study variance ${\tau }^2$ :

$${\tau }^2=\left\{\begin{array}{cc}{\displaystyle \frac{Q-df}{c}}\hfill & \mathrm{if}Q>df\hfill \\ 0\hfill & \mathrm{if}Q<df\hfill \end{array}\right.$$

(20)

where:

$$c=\sum _{i=1}^k{w}_i-{\displaystyle \frac{{\sum }_{i=1}^k{w}_i^2}{{\sum }_{i=1}^k{w}_i}}$$

(21)

and $df$ is the degrees of freedom, equal to $k-1$

3. Recalculate variance and weight

The variance is updated as:

$v_i^{*}=v_i+{\tau }^2$ .

The new weight becomes:

$w^{*}={\displaystyle \frac{1}{v_i^{*}}}$

4. Calculate the overall Fisher’s $Z_{overall}^{*}$

$$Z_{overall}^{*}={\displaystyle \frac{{\sum }_{i=1}^k{w}_i^{*}{z}_i}{{\sum }_{i=1}^k{w}_i^{*}}}$$

(22)

5. Transform back to $r_{\mathrm{overall}}$

Once $Z_{overall}^{*}$ is obtained, transform it back to overall $r_{overall}$ using Equation (16)

6. Final steps

After determining $r_{overall}$ using either fixed-effect or random-effect meta-analysis, the overall covariance can be calculated by Equation (12).

5. Example: Analyzing the 10 Cohort Studies in Table 1

5.1. Meta-Analysis for ${\beta }_1$

We wrote custom R functions to perform meta-analysis on the data in Table 1. The following R functions conduct both fixed-effect and random-effect meta-analyses for ${\beta }_1$.

studies <- read.csv(’https://raw.githubusercontent.com/enwuliu/

meta-analysis/main/random_effect_meta_sim.csv’, header = TRUE)

b1<-studies$b1  # beta1 from the 10 cohorts

var_b1 <- studies$var_b1  # variances of the coefficients

fixed_effect_meta <- function(B, V) {

  W <- 1 / V

  Beta <- sum(W * B) / sum(W)

  Var <- 1 / sum(W)

  resultlist <- list(’Overall beta’ = Beta, ’Overall variance’ = Var)

  return(resultlist)

}

random_effect_meta <- function(B, V) {

  W <- 1 / V

  Q <- sum(W * B^2) - (sum(W * B))^2 / sum(W)

  df <- length(B) - 1

  c <- sum(W) - sum(W^2) / sum(W)

  tau_square <- ifelse(Q > df, (Q - df) / c, 0)

  V_star <- V + tau_square

  W_star <- 1 / V_star

  Beta_star <- sum(W_star * B) / sum(W_star)

  Var_star <- 1 / sum(W_star)

  resultlist <- list(

 ’Overall beta’ = Beta_star,

 ’Overall variance’ = Var_star)

  return(resultlist)

}

fixed_b1 <- fixed_effect_meta(b1, var_b1)

fixed_b1

# $`Overall beta`

# [1] 1.040411

#

# $`Overall variance`

# [1] 0.6841809

random_b1 <- random_effect_meta(b1, var_b1)

random_b1

# $`Overall beta`

# [1] 1.014141

#

# $`Overall variance`

# [1] 0.7308388

From the results, the overall fixed-effect estimate for ${\beta }_1$ is 1.040411, with an overall variance of 0.6841809. The overall random-effect estimate is 1.014141, with an overall variance of 0.7308388.

5.2. Meta-Analysis for ${\beta }_2$

We used the same R functions to perform meta-analysis for ${\beta }_2$:

b2 <- studies$b2  # beta2 from the 10 cohorts

var_b2 <- studies$var_b2  # variances of the coefficients

fixed_b2<-fixed_effect_meta(b2, var_b2)

fixed_b2

# $`Overall beta`

# [1] -0.01140614

#

# $`Overall variance`

# [1] 0.0001403961

random_b2<-random_effect_meta(b2, var_b2)

random_b2

# $`Overall beta`

# [1] -0.01140614

#

# $`Overall variance`

# [1] 0.0001403961

Both fixed-effect and random-effect meta-analyses yield an overall interaction effect of ${\beta }_2$=-0.01140614, with a variance of 0.0001403961.

5.3. Meta-Analysis for r

The following R functions were used to conduct a meta-analysis of the correlation coefficients

fixed_effect_meta_r<-function(v1,v2,cov12,sample_size){

  r<-cov12/sqrt(v1*v2)                   #correlation coefficient

  z<-0.5*log((1+r)/(1-r))                #Fisher z transformation

  v<-1/(sample_size-3)                   #variance for z

  W<-1/v                                 #weight

  z_overall<-sum(W*z)/sum(W)             #overall random effect of z

  r_overall<-(exp(2*z_overall)-1)/(exp(2*z_overall)+1) #transform back to r

  return(r_overall)

}

random_effect_meta_r<-function(v1,v2,cov12,sample_size){

  r<-cov12/sqrt(v1*v2)                   #correlation coefficient

  z<-0.5*log((1+r)/(1-r))                #Fisher z transformation

  v<-1/(sample_size-3)                   #variance for z

  W<-1/v                                 #weight

  Q<-sum(W*z^2)-(sum(W*z))^2/sum(W)      #total variance or Q statistics

  c<-sum(W)-sum(W^2)/sum(W)              #c is the scaling factor

  df=length(v1)-1                        #degree of freedom

  tau_square<-ifelse(Q>df,(Q-df)/c,0)    #between study variance

  v_star=v+tau_square

  W_star=1/v_star                        #new weight

  z_star<-sum(W_star*z)/sum(W_star)      #overall random effect of z

  r_overall<-(exp(2*z_star)-1)/(exp(2*z_star)+1) #transform back to r

  return(r_overall)

}

v1<-studies$var_b1

v2<-studies$var_b2

v12<-studies$cov_b1b2

sample_size<-studies$sample_size

r<-v12/sqrt(v1*v2)

fixed_meta_r<-fixed_effect_meta_r(v1,v2,v12,sample_size)

fixed_meta_r

#[1] -0.9600286

random_meta_r<-random_effect_meta_r(v1,v2,v12,sample_size)

random_meta_r

#[1] -0.9493409

The overall r for fixed-effect meta-analysis is -0.9600286, while the random-effect meta-analysis gives -0.9493409.

5.4. Overall Covaraince for ${\beta }_1$ and ${\beta }_2$

Using Equation (12), the overall covariance is calculated as follows:

• Fixed-effect model:

  $$\begin{array}{cc}\hfill cov{({\widehat{\beta }}_1,{\widehat{\beta }}_2)}_{fixed}& =r_{fixed}\times \sqrt{V_{fixed{\widehat{\beta }}_1}{V}_{fixed{\widehat{\beta }}_2}}\hfill \\ \hfill & =-0.9600286\times \sqrt{0.6841809\times 0.0001403961}\hfill \\ \hfill & =-0.009409082\hfill \end{array}$$
• Random-effect model:

  $$\begin{array}{cc}\hfill cov{({\widehat{\beta }}_1,{\widehat{\beta }}_2)}_{random}& =r_{random}\times \sqrt{V_{random{\widehat{\beta }}_1}{V}_{random{\widehat{\beta }}_2}}\hfill \\ \hfill & =-0.9493409\times \sqrt{0.7308388\times 0.0001403961}\hfill \\ \hfill & =-0.009616357\hfill \end{array}$$

5.5. Calculating the Dffect of Falls on Hip Fracture at Different Ages and Its 95% Confidence Intervals

Using the information derived above, we can calculate the effect of falls on hip fracture at different ages and its corresponding 95% confidence intervals using Equation (2). Additionally, we generate a plot of the rate ratio (RR) for hip fracture, comparing fallers versus non-fallers at various ages, by exponentiating the estimated $\beta$ coefficients.

To illustrate, we calculate the effect of falls on hip fracture for ages 50 to 90. First, we conduct a fixed-effect meta-analysis, followed by a random-effect meta-analysis.

The following R code calculates the interaction effect of falls with age and the corresponding 95% confidence intervals, using Equation (2), based on the results of the fixed- and random-effect meta-analyses.

library(ggplot2)

library(ggpubr)

plot_with_interaction <- function(age, b1, v1, b2, v2, r) {

cov_b1b2 <- r * sqrt(v1 * v2)

RR <- exp(b1 + b2 * age)

RR_lower <- exp((b1+b2*age)-1.96 * sqrt(v1+age^2 * v2+2 * age*cov_b1b2))

RR_upper <- exp((b1+b2*age)+1.96*sqrt(v1+age^2 * v2+2*age*cov_b1b2))

ndata <- as.data.frame(cbind(RR, RR_lower, RR_upper, age))

ggplot(data = ndata, aes(x = age)) +

geom_line(aes(y = RR)) +

geom_line(aes(y = RR_lower), color = "steelblue", linetype = "dashed") +

geom_line(aes(y = RR_upper), color = "steelblue", linetype = "dashed") +

scale_y_continuous(name = "Rate Ratio (RR): Fallers vs. Non-Fallers",

                       limits = c(0, 4), expand = c(0, 0)) +

theme_bw() +

theme(axis.line = element_line(color = ’black’),

          axis.title.y = element_text(size = 8),

          axis.title.x = element_text(size = 8),

          axis.text.y = element_text(size = 8),

          axis.text.x = element_text(size = 8),

          plot.background = element_blank(),

          panel.grid.major = element_blank(),

          panel.grid.minor = element_blank(),

          panel.border = element_blank())

}

# Fixed-effect model

age <- seq(50, 90, 1)

b1 <- fixed_b1$`Overall beta`

v1 <- fixed_b1$`Overall variance`

b2 <- fixed_b2$`Overall beta`

v2 <- fixed_b2$`Overall variance`

r <- fixed_meta_r

p1 <- plot_with_interaction(age, b1, v1, b2, v2, r)

p1

# Random-effect model

b1_rand <- random_b1$`Overall beta`

v1_rand <- random_b1$`Overall variance`

b2_rand <- random_b2$`Overall beta`

v2_rand <- random_b2$`Overall variance`

r_rand <- random_meta_r

p2 <- plot_with_interaction(age, b1_rand, v1_rand, b2_rand, v2_rand, r_rand)

p2

p3 <- ggarrange(p1, p2,

                labels = c("Fixed Effect", "Random Effect"),

                ncol = 2, nrow = 1)

p3

6. An R Package for Calculating Overall Covariance in Meta-Analysis

To facilitate these calculations, we developed an R package called ’covmeta’, available on GitHub. This package can be installed and used within the R environment. Below is the R code to install the package and calculate the overall covariance using the data in Table 1:

#install the package

library(devtools)

install_github("enwuliu/covmeta")

library(covmeta)

# Load the dataset

studies <- read.csv(’https://raw.githubusercontent.com/enwuliu/

meta-analysis/main/random_effect_meta_sim.csv’, header = TRUE)

# Function arguments

b1 <- studies$b1          # Beta1 coefficients from the 10 cohorts

v1 <- studies$var_b1      # Variances of Beta1

b2 <- studies$b2          # Beta2 coefficients from the 10 cohorts

v2 <- studies$var_b2      # Variances of Beta2

cov_b1b2 <- studies$cov_b1b2  # Covariance between Beta1 and Beta2

sample_size <- studies$sample_size  # Sample sizes

# Calculate the overall main effect, interaction effect, and covariance

# Fixed-effect meta-analysis

cov_meta(b1, v1, b2, v2, cov_b1b2, sample_size, ’fixed’)

# Random-effect meta-analysis

cov_meta(b1, v1, b2, v2, cov_b1b2, sample_size, ’random’)

7. Discussion

This paper proposes a method for calculating overall covariance for regression coefficients in the presence of interactions in a two-stage meta-analysis. The approach utilizes the statistical relationship between the correlation coefficient and variances to estimate the covariance. By conducting separate meta-analyses, the variances and correlation coefficient can be obtained. While statistical software packages provide built-in capabilities for analyzing the overall main effect, interaction effects, and their variances,[11,12] no studies have explicitly demonstrated how to merge covariances in meta-analysis in the presence of interactions. We present a transparent and straightforward method to synthesize covariance in meta-analysis.

Detecting interactions is crucial in medical research, as interaction analyses can help determine whether an intervention is more effective for certain individuals.[13,14,15] In a two-stage meta-analysis, interactions between an exposure and covariates can be explored using meta-regression. However, the meta-regression method cannot study patient-level factors[16] and often suffers from low statistical power.[17] An alternative is dividing participants into subgroups (e.g., by age) and performing separate meta-analyses for each subgroup. While this avoids synthesizing interactions, it also tends to have low power.[18]

Another method for investigating covariances in meta-analysis is multivariate meta-analysis, which synthesizes correlated effects. For instance, in hypertension trials, systolic and diastolic blood pressure outcomes can be pooled using this approach.[19,20,21] In multivariate meta-analysis, the first-stage analysis estimates the effects (e.g., $\beta$ coefficients) and their covariances, which serve as inputs for the second-stage analysis conducted using mixed-model regression.[22,23] However, in this framework, the two correlated variables are treated as outcome variables rather than predictors.

Our method has some limitations. It assumes a linear relationship between the two regression coefficients, which may not always hold.[24]Additionally, the estimation of the correlation coefficient may be less reliable when the number of included studies is small.[25]

In conclusion, we have introduced a simple and transparent method for calculating covariance and confidence intervals in meta-analysis when interactions are present.