# How can I perform bivariate random-effects meta-regression if the implied between-study covariance matrix is not positive definite?

I have described my problem here:

How can one get consistent (i.e. direct+indirect=total) effects in a Meta-Analytic SEM model with latent variables?

Unfortunately, the solution I've found to this problem led me to another problem. In fact, I wish to perform a bivariate random-effects (RE) meta-regression, but in this way I get a between-study (BS) covariance matrix not positive definite (due to a BS correlation below -1).

How can I consider the effect estimates jointly in a setting like that?

See below for details.

I have 6 studies, each of which led me to a variance/covariance matrix of 2 effects. All of such matrices imply indeed a negative correlation (between -0.04 and -0.27) between the two effects (given they are a direct and an indirect effect, it was perfectly expected, since they somehow compete with each other - if the total effect were fixed, their correlation would actually be -1). See below:

> VCA<-matrix(c( 0.003844, -0.0002555, -0.0002555, 0.001156), nrow=2, byrow=TRUE)
> VCB<-matrix(c(0.005184, -0.000274, -0.000274, 0.001764), nrow=2, byrow=TRUE)
> VCC<-matrix(c(0.000784, -5.2e-05, -5.2e-05, 0.000049), nrow=2, byrow=TRUE)
> VCD<-matrix(c( 0.002401, -0.000203, -0.000203, 0.001369), nrow=2, byrow=TRUE)
> VCE<-matrix(c(0.001156 , -0.000186, -0.000186, 0.000441), nrow=2, byrow=TRUE)
> VCF<-matrix(c(0.009801 ,  -4.05e-05,   -4.05e-05, 0.000081), nrow=2, byrow=TRUE)
>
> cov2cor(VCA)
[,1]       [,2]
[1,]  1.0000000 -0.1212049
[2,] -0.1212049  1.0000000
> cov2cor(VCB)
[,1]        [,2]
[1,]  1.00000000 -0.09060847
[2,] -0.09060847  1.00000000
> cov2cor(VCC)
[,1]       [,2]
[1,]  1.0000000 -0.2653061
[2,] -0.2653061  1.0000000
> cov2cor(VCD)
[,1]       [,2]
[1,]  1.0000000 -0.1119691
[2,] -0.1119691  1.0000000
> cov2cor(VCE)
[,1]       [,2]
[1,]  1.0000000 -0.2605042
[2,] -0.2605042  1.0000000
> cov2cor(VCF)
[,1]        [,2]
[1,]  1.00000000 -0.04545455
[2,] -0.04545455  1.00000000


If I perform a multivariate random-effects meta-regression on the two effects in Stata, I have the problem that the between-study covariance is estimated to be -1, i.e. a value on the boundary (technically, I think that if the BS correlation really were -1 we would have a semi-positive definite covariance matrix). This leads to estimation problems. For example, in R (by using the "metaSEM" library) the RE covariance matrix is simply not estimated, so we have an identification problem:

> random.ma1 <- meta(y=cbind(direct,ind), v=cbind(v_direct, c_dind, v_ind),  data=covmatrix, model.name="Random effects model")


summary(random.ma1)

Call:
meta(y = cbind(direct, ind), v = cbind(v_direct, c_dind, v_ind),
data = covmatrix, model.name = "Random effects model")

95% confidence intervals: z statistic approximation
Coefficients:
Estimate   Std.Error      lbound      ubound z value Pr(>|z|)
Intercept1 -6.2550e-02  5.6259e+01 -1.1033e+02  1.1020e+02 -0.0011   0.9991
Intercept2 -2.3011e-02  6.5564e+02 -1.2850e+03  1.2850e+03  0.0000   1.0000
Tau2_1_1    5.6628e-03          NA          NA          NA      NA       NA
Tau2_2_1   -2.5553e-03          NA          NA          NA      NA       NA
Tau2_2_2    1.0153e-03          NA          NA          NA      NA       NA

Q statistic on the homogeneity of effect sizes: 38.39708
Degrees of freedom of the Q statistic: 10
P value of the Q statistic: 3.236141e-05

Heterogeneity indices (based on the estimated Tau2):
Estimate
Intercept1: I2 (Q statistic)   0.7187
Intercept2: I2 (Q statistic)   0.8151

Number of studies (or clusters): 6
Number of observed statistics: 12
Number of estimated parameters: 5
Degrees of freedom: 7
-2 log likelihood: -37.73167
OpenMx status1: 5 ("0" or "1": The optimization is considered fine.
Other values may indicate problems.)


The fact that the BS-covariance matrix is not positive definite is confirmed to me by the following two analyses. If I estimate the between-study correlation, I get an impossible value (-1.07):

> meta1 <- meta(cbind(direct,ind), cbind(v_direct, c_dind, v_ind),data=covmatrix)
> coef1 <- coef(meta1, select="random")
> my.cov <- vec2symMat(coef1, byrow=TRUE)
> cov2cor(my.cov)
[,1]      [,2]
[1,]  1.000000 -1.065654
[2,] -1.065654  1.000000


The same problem occurs (in this case, with an estimated BS correlation of -1.10) if I estimate the variance/covariance matrix of random effects from three separated single-variable meta-analyses (one for each effect, and one for their sum, i.e. for the direct, indirect and total effect separately):

> A<-matrix(c(0.0079, -0.002592, -0.002592, 0.0007), ncol=2, nrow=2)
> cov2cor(A)
[,1]      [,2]
[1,]  1.000000 -1.102231
[2,] -1.102231  1.000000


EDIT

Here:

(pdf)

is an explanation of the reasons why the problem occurs with a small number of studies. With complete data and an estimated BS correlation of +1 or -1, bivariate RE meta-analysis becomes less efficient than the univariate one. However, no solution to this problem is given.

• I am not sure whether you have already tried it... but multivariate meta-regression may be implemented without prespecified variances, for instance defining a priory a given correlation between point effect estimates. I would suggest you to try this approach Nov 9, 2018 at 12:13
• For each study, there are w estimated effects and the 2x2 covariance matrix of the estimates. If I then perform MMR by leaving the BS covariance unstructured, the model either doesn't run or estimates a BS correlation of -1. If you mean to just use individual-patient data, I haven't found any command allowing me to do what you say, either in R or in Stata. However I've done something similar by hand, by estimating the two individual effects separately, and then assuming a null correlation between the two to estimate the total effect (given the correlation is negative, a conservative approach) Nov 9, 2018 at 14:28