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
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 at r-project.org] On Behalf Of Ana Ben?tez
Sent: Thursday, 19 September, 2019 12:14
To: r-sig-meta-analysis at 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
Dear Wolfgang,
Apologies for sending the text in rich-text format. Here I copy-paste the
original message and I additionally attach a piece of code to solve the
issue. I got the code from Alfredo S?nchez-Tojar but he claims he got it
from someone else, maybe even you Wolfgang. The code is this:
# this function is to make sure that a is matrix positive-definitive
# it uses the function nearPD in the Matrix package to compute the
# nearest positive definite matrix to an approximate one
# Function reused from: https://osf.io/ayxrt/
PDfunc <- function(matrix){
require(Matrix)
if(corpcor::is.positive.definite(matrix) == T){
x <- corDiag
}else{
x <- Matrix::nearPD(matrix)
}
mat <- x$mat
return(mat)
}
# is the matrix positive-definitive? No, run PDfunc to make it
# so. If a matrix is not positive-definite, it won't work
is.positive.definite(varcovar) # FALSE
varcovar <-PDfunc(varcovar)
is.positive.definite(varcovar) # TRUE
And this is my original email in case someone is interested:
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
El mar., 24 sept. 2019 a las 13:12, Viechtbauer, Wolfgang (SP) (<
wolfgang.viechtbauer at maastrichtuniversity.nl>) escribi?:
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 at r-project.org] On Behalf Of Ana Ben?tez
Sent: Thursday, 19 September, 2019 12:14
To: r-sig-meta-analysis at 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 <-
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
Hi Ana,
No, this code is not from me, but I know the nearPD() function. But that's just slapping a bandaid on, not really fixing the issue. In principle, for the scenario you are describing, the V matrix should automatically be positive definite. The fact that it isn't suggest that there might be a more fundamental problem. If so, simply adjust the V matrix to the 'nearest' PD matrix isn't a proper fix for that.
Best,
Wolfgang
-----Original Message-----
From: Ana Ben?tez [mailto:abenitez81 at gmail.com]
Sent: Thursday, 26 September, 2019 13:33
To: Viechtbauer, Wolfgang (SP)
Cc: r-sig-meta-analysis at r-project.org
Subject: Re: [R-meta] Variance-covariance is not positive definite.
Dear Wolfgang,
Apologies for sending the text in rich-text format. Here I copy-paste the original message and I additionally attach a piece of code to solve the issue. I got the code from Alfredo S?nchez-Tojar but he claims he got it from someone else, maybe even you Wolfgang. The code is this:
# this function is to make sure that a is matrix positive-definitive
# it uses the function nearPD in the Matrix package to compute the
# nearest positive definite matrix to an approximate one
# Function reused from: https://osf.io/ayxrt/
PDfunc <- function(matrix){
? require(Matrix)
? if(corpcor::is.positive.definite(matrix) == T){
? ? x <- corDiag
? }else{
? ? x <- Matrix::nearPD(matrix)
? }
? mat <- x$mat
? return(mat)
}
# is the matrix positive-definitive? No, run PDfunc to make it
# so. If a matrix is not positive-definite, it won't work
is.positive.definite(varcovar) # FALSE
?varcovar ?<-PDfunc(varcovar)
is.positive.definite(varcovar) # TRUE
And this is my original email in case someone is interested:
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
El mar., 24 sept. 2019 a las 13:12, Viechtbauer, Wolfgang (SP) (<wolfgang.viechtbauer at maastrichtuniversity.nl>) escribi?:
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