[R-meta] Non-positive definite covariance matrix for rma.mv
Pia-Magdalena Schmidt
p|@-m@gd@|en@@@chm|dt @end|ng |rom un|-bonn@de
Wed Sep 4 15:05:29 CEST 2024
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 using 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 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
> _______________________________________________
> 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
More information about the R-sig-meta-analysis
mailing list