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

Thu Sep 19 12:13:46 CEST 2019

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/

	[[alternative HTML version deleted]]