Dear Pia,
You use vcalc() as follows:
V <- vcalc(vi=ES_all$vi, cluster=DatMA$id_database)
This assumes a correlation of 1 within clusters. An
example:
dat <- data.frame(study=c(1,1,1,2,3,3), vi=runif(6, .01,
.10))
V <- vcalc(vi, cluster=study, data=dat)
blsplit(V, dat$study, fun=cov2cor)
I don't think this is what you intended. If you have
multiple observations within clusters, then you indicate
this to the function via 'obs' combined with 'rho' to
provide the correlation:
dat <- data.frame(study=c(1,1,1,2,3,3),
obs=c(1,2,3,1,1,2), vi=runif(6, .01, .10))
V <- vcalc(vi, cluster=study, obs=obs, data=dat,
rho=0.6)
blsplit(V, dat$study, fun=cov2cor)
Best,
Wolfgang
-----Original Message-----
From: Pia-Magdalena Schmidt
<pia-magdalena.schmidt at uni-bonn.de>
Sent: Wednesday, September 4, 2024 15:05
To: R Special Interest Group for Meta-Analysis
<r-sig-meta-analysis at r-
project.org>; James Pustejovsky <jepusto at gmail.com>
Cc: Viechtbauer, Wolfgang (NP)
<wolfgang.viechtbauer at maastrichtuniversity.nl>
Subject: Re: [R-meta] Non-positive definite covariance
matrix for rma.mv
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