[R-meta] Non-positive definite covariance matrix for rma.mv

[Viechtbauer, Wolfgang (NP)]

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 using 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 using r-
> project.org>; James Pustejovsky <jepusto using gmail.com>
> Cc: Dr. Gerta Rücker <gerta.ruecker using 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 using 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 using r-project.org> Im Auftrag
> von Andreas Voldstad via R-sig-meta-analysis
> Gesendet: Montag, 1. Juli 2024 00:23
> An: James Pustejovsky <jepusto using gmail.com>; R Special Interest Group for Meta-
> Analysis <r-sig-meta-analysis using r-project.org>
> Cc: Andreas Voldstad <andreas.voldstad using 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 using gmail.com>
> Sent: Sunday, June 30, 2024 5:13:32 PM
> To: R Special Interest Group for Meta-Analysis <r-sig-meta-analysis using r-
> project.org>
> Cc: Andreas Voldstad <andreas.voldstad using 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 using 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