Dear all, although this issue has been raised recently (see below), I would highly appreciate your advice, particularly regarding whether I might use nearPD to approximate a PD V matrix. I am running a meta-analysis with 30 studies, 11 of which report more than one effect size (2 or 3). Therefore, I would like to fit a random-effects model using the rma.mv function to account for these dependencies. I used vcalc to approximate the var cov matrix and received the following warning for all clusters (11) with more than one effect size: ?The var-cov matrix appears to be not positive definite in cluster x?. The correlation (r= 0.6) used is an approximation based on available raw data. I inspected the eigenvalues of the according submatrices and found 3 negative eigenvalues close to 0, else values near 0 (see the values below). If I round the values accordingly to your previous emails, the matrix is semi-PD. I am wondering whether I might use the nearPD function to be able to fit the rma.mv model to my data? Best, Pia ES_all <- escalc(measure="SMCC", m1i=DatMA$'1_drug_mean', sd1i=DatMA$'1_drug_sd', ni = DatMA$'1_n', m2i=DatMA$'1_plc_mean', sd2i=DatMA$'1_plc_sd', pi=DatMA$'1_p_value', ri = DatMA$'1_ri') V <- vcalc(vi=ES_all $vi, cluster=DatMA$id_database) res <- rma.mv(yi=ES_all$yi, V, random = ~ 1 | id_database/effect_id, data = DatMA) eigenvalues: $`1` [1] 0.335876 $`2` [1] 6.788500e+00 3.552714e-15 0.000000e+00 $`3` [1] 0.2049939 $`4` [1] 0.1275193 $`5` [1] 2.792778e-01 1.387779e-17 $`6` [1] 0.1974481 $`7` [1] 0.2392147 $`8` [1] 0.09561453 $`9` [1] 0.1023009 $`10` [1] 7.202168e-02 6.938894e-18 $`11` [1] 0.08839118 $`12` [1] 0.05667013 $`13` [1] 2.160967e-01 5.551115e-17 -1.387779e-17 $`14` [1] 1.906503e-01 5.551115e-17 0.000000e+00 $`15` [1] 0.2379626 $`16` [1] 1.060134e+00 5.551115e-17 $`17` [1] 4.879074e-02 1.734723e-18 $`18` [1] 0.1956611 $`19` [1] 0.1319664 $`20` [1] 4.783856e-01 1.665335e-16 2.775558e-17 $`21` [1] 0.07701826 $`22` [1] 0.08400946 $`23` [1] 0.4219658 $`24` [1] 3.485124e+00 -5.551115e-17 $`25` [1] 4.095845e-01 -2.775558e-17 $`26` [1] 0.1124492 $`27` [1] 0.09964694 $`28` [1] 2.130419e-01 6.938894e-18 $`29` [1] 0.2391303 $`30` [1] 0.122081 round(eigen(V)$values, 8) [1] 6.78849973 3.48512378 1.06013363 0.47838557 0.42196575 0.40958453 0.33587598 0.27927780 0.23921468 0.23913035 0.23796263 0.21609670 0.21304191 [14] 0.20499394 0.19744806 0.19566106 0.19065030 0.13196637 0.12751926 0.12208097 0.11244916 0.10230087 0.09964694 0.09561453 0.08839118 0.08400946 [27] 0.07701826 0.07202168 0.05667013 0.04879074 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 [40] 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 On Mo, 1 Jul 2024 12:57:50 +0000 Viechtbauer, Wolfgang (NP) via R-sig-meta-analysis
<r-sig-meta-analysis at r-project.org> wrote:
Hi Gerta,
Actually, the var-cov matrix of all pairwise comparisons
is not PD. Consider these data:
dat <- data.frame(
study = c(1,1,1),
m1i = c(10,10,8),
m2i = c(8,7,7),
sd1i = c(1.2,1.2,1.4),
sd2i = c(1.4,1.3,1.3),
n1i = c(45,45,42),
n2i = c(42,40,40),
grp1 = c("A", "A", "B"),
grp2 = c("B", "C", "C"))
dat
We can compute all pairwise mean differences with:
dat <- escalc(measure="MD", m1i=m1i, sd1i=sd1i, n1i=n1i,
m2i=m2i, sd2i=sd2i, n2i=n2i,
data=dat)
dat
The variance-covariance matrix of the mean differences
can be constructed with:
V <- vcalc(vi, cluster=study, grp1=grp1, grp2=grp2,
w1=n1i/sd1i^2, w2=n2i/sd2i^2, data=dat)
V
We can double-check that these values are correct:
- yi[1] and yi[2] share group A, so the covariance
between them must be:
dat$sd1i[1]^2 / dat$n1i[1]
- yi[1] and yi[3] share group B (but in different
positions), so the covariance is:
-dat$sd2i[1]^2 / dat$n2i[1]
- yi[2] and yi[3] share group C, so the covariance is:
dat$sd2i[2]^2 / dat$n2i[2]
So these are all correct. But V is not PD. It is
semi-PD, with a third eigenvalue exactly equal to 0:
round(eigen(V)$values, 8)
Best,
Wolfgang
-----Original Message----- From: R-sig-meta-analysis <r-sig-meta-analysis-bounces at r-project.org> On Behalf Of Dr. Gerta R?cker via R-sig-meta-analysis Sent: Monday, July 1, 2024 14:18 To: R Special Interest Group for Meta-Analysis <r-sig-meta-analysis at r- project.org>; James Pustejovsky <jepusto at gmail.com> Cc: Dr. Gerta R?cker <gerta.ruecker at uniklinik-freiburg.de> Subject: Re: [R-meta] Non-positive definite covariance matrix for rma.mv Hi, I surely do not fully comprehend the problem. But what I do not understand is why fully consistent treatment effects should lead to a non-positive definite covariance matrix. In network meta-analysis, we always estimate fully consistent treatment effects for all pairwise comparisons, and the full covariance matrix will always be positive definite. This has nothing to do with the (by the way, unnecessary) choice of a baseline treatment. Or am I misunderstanding anything? Best, Gerta UNIVERSIT?TSKLINIKUM FREIBURG Institute for Medical Biometry and Statistics Dr. Gerta R?cker Guest Scientist Stefan-Meier-Stra?e 26 ? 79104 Freiburg gerta.ruecker at uniklinik-freiburg.de https://www.uniklinik-freiburg.de/imbi-en/employees.html?imbiuser=ruecker -----Urspr?ngliche Nachricht----- Von: R-sig-meta-analysis <r-sig-meta-analysis-bounces at r-project.org> Im Auftrag von Andreas Voldstad via R-sig-meta-analysis Gesendet: Montag, 1. Juli 2024 00:23 An: James Pustejovsky <jepusto at gmail.com>; R Special Interest Group for Meta- Analysis <r-sig-meta-analysis at r-project.org> Cc: Andreas Voldstad <andreas.voldstad at kellogg.ox.ac.uk> Betreff: Re: [R-meta] Non-positive definite covariance matrix for rma.mv Hi James, thank you for responding. Indeed, I am fitting a model across all the contrasts, and metafor did yell at me for the NPD V matrix. Is this a case where it might be appropriate to use the nearpd = T argument in vcalc? Best wishes, Andreas
________________________________
From: James Pustejovsky <jepusto at gmail.com>
Sent: Sunday, June 30, 2024 5:13:32 PM
To: R Special Interest Group for Meta-Analysis
<r-sig-meta-analysis at r-
project.org>
Cc: Andreas Voldstad <andreas.voldstad at kellogg.ox.ac.uk>
Subject: Re: [R-meta] Non-positive definite covariance
matrix for rma.mv
Hi Andreas,
I think the NPD issue arises from having treatment
contrasts that are perfectly
collinear. Consider that if you know the contrasts Tx1
vs Passive, Tx2 vs
Passive, and Tx1 vs Active, then you can infer the Tx2
vs Active contrast:
(Tx2 vs Passive) - (Tx1 vs Passive) + (Tx1 vs Active) =
(Tx2 - Tx1) + (Tx1 vs
Active) = (Tx2 vs Active)
As far as I can see, the only way to avoid this is to
calculate all effect size
contrasts relative to a common comparison condition.
Depending on the meta-analytic model in which you use
this covariance matrix,
the NPD issue might or might not be important?it really
depends on the structure
of the model. For example, if you are going to run
separate models with the
active control contrasts and passive control contrasts,
then the fact that the
whole covariance matrix is NPD is NBD?not a big deal.
But if you?re fitting some
sort of multilevel meta-analysis model across all the
contrasts, then metafor
will yell at you if you use an NBD matrix in the V
argument.
James
On Jun 30, 2024, at 8:37?AM, Andreas Voldstad via
R-sig-meta-analysis <r-sig-
meta-analysis at r-project.org> wrote:
Dear everyone,
I have one study in my meta-analysis with a covariance
matrix that is non-
positive definite.
I am wondering what consequences this might have for
my meta-analysis with
rma.mv, and whether you have any suggestions of how to
deal with the issue.
The study had four arms and effect sizes are
additionally correlated over time
and within dyads.
The correlations over time and dyad members used for
creating the covariance
matrix was based on the corresponding correlations
between participants' scores,
estimated from the study dataset.
Because of the four arms, some effect sizes have no
overlap, leading to blocks
of zeros in the covariance matrix.
The covariance matrix is reproduced below:
study_escalc <- data.frame( StudyID = c("StudyA",
"StudyA", "StudyA",
"StudyA", "StudyA", "StudyA", "StudyA", "StudyA",
"StudyA", "StudyA", "StudyA",
"StudyA", "StudyA", "StudyA", "StudyA", "StudyA"), Role
= c("Husband",
"Husband", "Wife", "Wife", "Husband", "Husband", "Wife",
"Wife", "Husband",
"Husband", "Wife", "Wife", "Husband", "Husband", "Wife",
"Wife"), Timenum =
c(2, 3, 2, 3, 2, 3, 2, 3, 2, 3, 2, 3, 2, 3, 2, 3), Tx =
c("Tx1", "Tx1", "Tx1",
"Tx1", "Tx2", "Tx2", "Tx2", "Tx2", "Tx1", "Tx1", "Tx1",
"Tx1", "Tx2", "Tx2",
"Tx2", "Tx2"), Control = c("Passive", "Passive",
"Passive", "Passive",
"Passive", "Passive", "Passive", "Passive", "Active",
"Active", "Active",
"Active", "Active", "Active", "Active", "Active"),
Contrast = c("Tx1 vs
Passive", "Tx1 vs Passive", "Tx1 vs Passive", "Tx1 vs
Passive", "Tx2 vs
Passive", "Tx2 vs Passive", "Tx2 vs Passive", "Tx2 vs
Passive", "Tx1 vs Active",
"Tx1 vs Active", "Tx1 vs Active", "Tx1 vs Active", "Tx2
vs Active", "Tx2 vs
Active", "Tx2 vs Active", "Tx2 vs Active"), delta =
c(0.467304, 0.429311,
0.145277, 0.248136, 0.684523, 0.537425, 0.380137,
0.335112, 0.432723, 0.480095,
0.145031, 0.252715, 0.628890, 0.572478, 0.373187,
0.338829), v = c(0.053375,
0.053154, 0.051417, 0.051710, 0.056543, 0.056865,
0.052208, 0.052756, 0.053863,
0.053454, 0.050774, 0.053052, 0.056766, 0.057132,
0.051522, 0.054108), N_tx =
c(38, 39, 39, 38, 36, 35, 39, 37, 38, 39, 39, 38, 36,
35, 39, 37), N_Control =
c(39, 38, 39, 40, 39, 38, 39, 40, 38, 38, 40, 38, 38,
38, 40, 38), Category =
c("Husband_2", "Husband_3", "Wife_2", "Wife_3",
"Husband_2", "Husband_3",
"Wife_2", "Wife_3", "Husband_2", "Husband_3", "Wife_2",
"Wife_3", "Husband_2",
"Husband_3", "Wife_2", "Wife_3"), Effect = c("Tx1 vs
Passive_Husband_2", "Tx1
vs Passive_Husband_3", "Tx1 vs Passive_Wife_2", "Tx1 vs
Passive_Wife_3", "Tx2 vs
Passive_Husband_2", "Tx2 vs Passive_Husband_3", "Tx2 vs
Passive_Wife_2", "Tx2 vs
Passive_Wife_3", "Tx1 vs Active_Husband_2", "Tx1 vs
Active_Husband_3", "Tx1 vs
Active_Wife_2", "Tx1 vs Active_Wife_3", "Tx2 vs
Active_Husband_2", "Tx2 vs
Active_Husband_3", "Tx2 vs Active_Wife_2", "Tx2 vs
Active_Wife_3"))
study_cormat <- matrix(c( 1.000000, 0.781447,
0.689538, 0.639455, 0.781447,
1.000000, 0.565192, 0.653545, 0.689538, 0.565192,
1.000000, 0.820650,
0.639455, 0.653545, 0.820650, 1.000000), nrow=4,
byrow=TRUE, dimnames = list(
c("Husband_2", "Husband_3", "Wife_2", "Wife_3"),
c("Husband_2", "Husband_3",
"Wife_2", "Wife_3")))
Study_Vmat <- vcalc(vi = v, cluster =
StudyID, obs = Category,
grp1 = Tx, grp2 = Control, w1 =
N_tx, w2 =
N_Control, rho = study_cormat,
data = study_escalc)
# Warning message: The var-cov matrix appears to be
not positive definite in
cluster StudyA. cov2cor(Study_Vmat)
Warm wishes,
Andreas Voldstad (he/him)
PhD student in Psychiatry
University of Oxford
_______________________________________________ R-sig-meta-analysis mailing list @ R-sig-meta-analysis at r-project.org To manage your subscription to this mailing list, go to: https://stat.ethz.ch/mailman/listinfo/r-sig-meta-analysis