[R-meta] Variance-covariance is not positive definite.

Tue Sep 24 13:11:42 CEST 2019

Dear Ana,

Please send messages to this list in plain text (no HTML / rich-text emails). As you can see below, your message includes lots of superfluous empty lines, which makes it hard to read.

Based on what you posted, it is not possible to determine why the construction of the V matrix leads to a non-positive definite matrix. We would have to see the actual dataset (mamdata). The only thing I can suggest to check is that the common control group is always given by variables Mean_m, sd_m, and N_m and not the other three variables that are used for the comparison groups. Also, how was variable 'var' computed? Was it computed using escalc() or was this computed based on other software? In the latter case, the computation used may not match up with what the calc.v() function below uses for computing the covariances.

Best,
Wolfgang

-----Original Message-----
From: R-sig-meta-analysis [mailto:r-sig-meta-analysis-bounces using r-project.org] On Behalf Of Ana Benítez
Sent: Thursday, 19 September, 2019 12:14
To: r-sig-meta-analysis using r-project.org
Subject: [R-meta] Variance-covariance is not positive definite.

Dear Wolfgang and meta-analysis list subscribers,

I am running a muti-level meta-regression in which I have several
comparisons using a common control for some of the studies. My effect size
is the log response ratio (RR) and I have as random effects Study, Effect
size ID and Species nested in Family and Order. I am also rerunning the
same analysis using a phylogenetic correlation matrix instead of nesting
species in higher taxonomic taxa (see end of the message).

To calculate the var-covar matrix I have used Wolfgang’s code:

calc.v <- function(x) {

  v <- matrix(x$sd_m[1]^2 / (x$N_m[1] * x$Mean_m[1]^2), nrow=nrow(x),
ncol=nrow(x))

  diag(v) <- x$var

  v

}

#make sure we order the database by the variable we are splitting on
(CommonControl)

mamdata <- mamdata[order(mamdata$CommonControl),]

Vmam<- bldiag(lapply(split(mamdata, mamdata$CommonControl), calc.v))

head(Vmam)

Yet, when I run the meta-analysis, and although I get reasonable results, I
get the warning: 'V' appears to be not positive definite. See results here:

Multivariate Meta-Analysis Model (k = 705; method: REML)

   logLik   Deviance        AIC        BIC       AICc

 116.6641  -233.3282  -219.3282  -187.4407  -219.1671

Variance Components:

            estim    sqrt  nlvls  fixed                        factor

sigma^2.1  0.0149  0.1223    131     no                     Reference

sigma^2.2  0.0275  0.1657    705     no                            ID

sigma^2.3  0.0013  0.0363     13     no                         Order

sigma^2.4  0.0069  0.0831     37     no                  Order/Family

sigma^2.5  0.0242  0.1554    167     no  Order/Family/Species.nominal

Test of Moderators (coefficient(s) 2):

QM(df = 1) = 18.0854, p-val < .0001

Model Results:

                        estimate      se     zval    pval    ci.lb
ci.ub

intrcpt              0.1768  0.0636   2.7780  0.0055   0.0521   0.3015   **

logmass           -0.0880  0.0207  -4.2527  <.0001  -0.1286  -0.0475  ***

I have also tried computing the var-covar matrix using the script provided
by Moatt et al. 2016, The effect of dietary restriction on reproduction: a
meta-analytic perspective. BMC evolutionary biology, 16(1), 199, available
at https://datadryad.org/stash/dataset/doi:10.5061/dryad.3fc02. I slightly
modified the script to calculate the var-covar matrix for RR as:

# create square matrix matching N of ES, filled with zeros

V <- matrix(0,nrow = dim(mamdata)[1],ncol = dim(mamdata)[1])

rownames(V) <- mamdata$ID

colnames(V) <- mamdata$ID

# find start and end coordinates for the subsets

shared_coord <-
which(mamdata$CommonControl%in%mamdata$CommonControl[duplicated(mamdata$CommonControl)]==TRUE)

shared_coord

# matrix of combinations of coordinates for each experiment with shared
control

combinations <- do.call("rbind", tapply(shared_coord,
mamdata[shared_coord,"CommonControl"], function(x) t(combn(x,2))))

combinations

# calculate covariance values between ES values at the positions in
shared_list and place them on the matrix

for (i in 1:dim(combinations)[1]){

  p1 <- combinations[i,1]

  p2 <- combinations[i,2]

  p1_p2_cov <- mamdata$sd_m[1]^2 / (mamdata$N_m[1] * mamdata$Mean_m[1]^2)

  V[p1,p2] <- p1_p2_cov

  V[p2,p1] <- p1_p2_cov

}

# add the diagonal - use df$var as matrix diagonal

diag(V) <-  mamdata$var

Using this approach the results are pretty similar but I still get the
warning: 'V' appears to be not positive definite.

Multivariate Meta-Analysis Model (k = 705; method: REML)

   logLik   Deviance        AIC        BIC       AICc

 115.2436  -230.4872  -216.4872  -184.5997  -216.3261

Variance Components:

            estim    sqrt  nlvls  fixed                        factor

sigma^2.1  0.0153  0.1238    131     no                     Reference

sigma^2.2  0.0274  0.1656    705     no                            ID

sigma^2.3  0.0014  0.0370     13     no                         Order

sigma^2.4  0.0070  0.0836     37     no                  Order/Family

sigma^2.5  0.0242  0.1554    167     no  Order/Family/Species.nominal

Test of Moderators (coefficient(s) 2):

QM(df = 1) = 18.3750, p-val < .0001

Model Results:

                        estimate      se     zval    pval    ci.lb
ci.ub

intrcpt              0.1799  0.0639   2.8159  0.0049   0.0547   0.3051   **

logmass           -0.0890  0.0208  -4.2866  <.0001  -0.1298  -0.0483  ***

Needless to say that when I run the phylogenetic meta-analysis I also get
the warning. Results below:

Multivariate Meta-Analysis Model (k = 702; method: REML)

   logLik   Deviance        AIC        BIC       AICc

 115.2222  -230.4445  -218.4445  -191.1380  -218.3233

Variance Components:

            estim    sqrt  nlvls  fixed           factor    R

sigma^2.1  0.0135  0.1161    131     no        Reference   no

sigma^2.2  0.0276  0.1662    702     no               ID   no

sigma^2.3  0.0274  0.1654    164     no             SPID   no

sigma^2.4  0.0090  0.0951    161     no  Species.nominal  yes

Test of Moderators (coefficient(s) 2):

QM(df = 1) = 17.6159, p-val < .0001

Model Results:

                        estimate      se     zval    pval    ci.lb
ci.ub

intrcpt               0.1486  0.0744   1.9979  0.0457   0.0028   0.2943    *

logmass           -0.0792  0.0189  -4.1971  <.0001  -0.1162  -0.0422  ***

My question is whether I should be worried by the warning or not and if so,
how I can fix it.

Thanks in advance.

Best,

Ana

-- 
Ana Benítez López

Juan de la Cierva-Incorporación Postdoctoral Fellow
The Integrative Ecology Group (http://ebd10.ebd.csic.es/)
Estación Biológica de Doñana (CSIC)
Avenida Américo Vespucio, 26, 41092, Sevilla, Spain
Phone: 954 46 67 00 ext. 1451

 https://www.researchgate.net/profile/Ana_Benitez-Lopez

orcid.org/0000-0002-6432-1837

http://www.environmentalevidencejournal.org/