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

[James Pustejovsky]

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
> Please don’t feel obliged to read or respond to my email outside your own working hours.
> 
>    [[alternative HTML version deleted]]
> 
> _______________________________________________
> R-sig-meta-analysis mailing list @ R-sig-meta-analysis using r-project.org
> To manage your subscription to this mailing list, go to:
> https://stat.ethz.ch/mailman/listinfo/r-sig-meta-analysis