Message-ID: <f9b17d6e3fec428197427917d5455436@UM-MAIL3214.unimaas.nl>
Date: 2019-09-05T08:36:12Z
From: Wolfgang Viechtbauer
Subject: [R-meta] setting certain covariances to 0 for the random LHS term
In-Reply-To: <92D045D1-CBE6-4F89-8E73-C07D84174E48@cgiar.org>
Hi Gil,
As far as I can tell, you have specified rho and phi correctly, at least dimension wise, as otherwise you should get a different error. Also, the positions of the 0s appear correct.
The error indicates that there are problems with the optimization. Apparently, fixing those correlations to zero introduces non-positive definiteness into the correlation matrix. Let's take a simple example illustrating how this can happen:
M <- matrix(.8, nrow=5, ncol=5)
diag(M) <- 1
M
eigen(M)
# all eigenvalues > 0, so everything is fine
# now fix some elements in M to 0
M[3,4:5] <- M[4:5,3] <- 0
M[4,5] <- M[5,4] <- 0
M
eigen(M)
# last eigenvalue < 0, so M is not PD (or even semi PD)
So, in essence, this suggests that the data and the resulting correlations are not compatible with all those correlations being equal to 0.
I should also point out that you are trying to fit a rather complex model, with 17 (rhos) + 17 (phis) + 10 (tau^2) + 10 (gamma^2) + 2 (sigma^2) = 56 parameters for the random effects structure. I hope you have tons of data, as otherwise this is a futile exercise.
Best,
Wolfgang
-----Original Message-----
From: R-sig-meta-analysis [mailto:r-sig-meta-analysis-bounces at r-project.org] On Behalf Of Gram, Gil (IITA)
Sent: Thursday, 05 September, 2019 8:20
To: r-sig-meta-analysis at r-project.org
Subject: Re: [R-meta] setting certain covariances to 0 for the random LHS term
Hi Wolfgang,
Thanks for your reply, and my apologies for not having found it in the documentation myself!
However, I must be doing something wrong because I get an "error in .ll.rma.mv(opt.res$par, reml = reml, Y = Y, M = V, A = A, X.fit = X, : Final variance-covariance matrix not positive definite? and warnings "In cov2cor(E) : diag(.) had 0 or NA entries; non-finite result is doubtful?.
With classMR.classOR being the ?inner? factor with 10 levels = Control, MR, ORone, ORtwo, ORthree, ORManure, ORMRone, ORMRtwo, ORMRthree, ORMRManure
and two random effect terms ~classMR.classOR|idSite, ~ classMR.classOR|idSite.time
My rho (and phi) vectors in column-wise order are as such:
?12, ?13, ?23, ?14, ?24, ?34, ?15, ?25, ?35, ?45, ?16, ?26, ?36, ?46, ?56, ?17, ?27, ?37, ?47, ?57, ?67, ?18, ?28, ?38, ?48, ?58, ?68, ?78, ?19, ?29, ?39, ?49, ?59, ?69, ?79, ?89, ?1-10, ?2-10, ?3-10, ?4-10, ?5-10, ?6-10, ?7-10, ?8-10, ?9-10
I?m only interested in correlations involving 1 (Control) and 2 (MR), so those I set to ?NA? and the rest to 0.
So please correct me where I?m wrong:
res = rma.mv(sqrt(yi), vi, method='REML', struct=c("UN","UN"), sparse=TRUE, data=dat
, rho=c(NA,NA,NA,NA,NA,0,NA,NA,0,0,NA,NA,0,0,0,NA,NA,0,0,0,0,NA,NA,0,0,0,0,0,NA,NA,0,0,0,0,0,0,NA,NA,0,0,0,0,0,0,0)
, phi=c(NA,NA,NA,NA,NA,0,NA,NA,0,0,NA,NA,0,0,0,NA,NA,0,0,0,0,NA,NA,0,0,0,0,0,NA,NA,0,0,0,0,0,0,NA,NA,0,0,0,0,0,0,0)
, mods = ~ rateORone + rateORtwo + rateORthree + rateORfour + rateORManure + kgMN
+ I(rateORone^2) + I(rateORtwo^2) + I(rateORthree^2) + I(rateORfour^2) + I(rateORManure^2) + I(kgMN^2)
+ rateORone:kgMN + rateORtwo:kgMN + rateORthree:kgMN + rateORfour:kgMN + rateORManure:kgMN
+ I(rateORone^2):I(kgMN^2) + I(rateORtwo^2):I(kgMN^2) + I(rateORthree^2):I(kgMN^2) + I(rateORfour^2):I(kgMN^2) + I(rateORManure^2):I(kgMN^2)
+ cropSys + idF,
random = list(~1|ref, ~1|idRow, ~ classMR.classOR|idSite, ~ classMR.classOR|idSite.time))
Thanks a lot!
Gil Gram
PhD researcher | Natural resource management/CCAFS
International Institute of Tropical Agriculture (IITA)
East African Hub/Kampala/Country office
Address: IITA-Uganda Plot 15B, Naguru East Road. P.O Box 7878, Kampala
Mobile: +256 755 315236 | Belgium Mobile: +32 484 981200
Skype: gil.gram
On 5 Sep 2019, at 00:14, Viechtbauer, Wolfgang (SP) <wolfgang.viechtbauer at maastrichtuniversity.nl<mailto:wolfgang.viechtbauer at maastrichtuniversity.nl>> wrote:
Hi Gil,
I don't understand the specifics of your particular case, but yes, you can set covariances (or rather, correlations) to zero. Argument 'rho' and 'phi' are used for that. I assume you are using struct="UN" -- or rather, struct=c("UN","UN") to be precise -- for those two random effects terms. Take a look at:
https://wviechtb.github.io/metafor/reference/rma.mv.html
and find the part that starts with: "For struct="UN", the values specified under 'rho' should be given in column-wise order." Setting an element of that vector to NA means that the corresponding correlation will be estimated. Setting it to 0 sets it to ... 0!
To illustrate, using a very simple case where there are only two levels for the 'inner' term (so there really is only one correlation):
dat <- dat.berkey1998
V <- lapply(split(dat[,c("v1i", "v2i")], dat$trial), as.matrix)
V <- bldiag(V)
res <- rma.mv(yi, V, mods = ~ outcome - 1, random = ~ outcome | trial, struct="UN", data=dat)
res
# So, if I set rho = NA, then the correlation is estimated (as above):
res <- rma.mv(yi, V, mods = ~ outcome - 1, random = ~ outcome | trial, struct="UN", data=dat, rho=c(NA))
res
# But I can set the correlation also to 0:
res <- rma.mv(yi, V, mods = ~ outcome - 1, random = ~ outcome | trial, struct="UN", data=dat, rho=c(0))
res
# which is actually identical here to:
res <- rma.mv(yi, V, mods = ~ outcome - 1, random = ~ outcome | trial, struct="DIAG", data=dat)
res
Now if you have more than one correlation, then rho will be a vector. And you can pick and choose which correlations to constrain to 0 and which ones will be estimated.
Argument 'phi' works the same way -- just for the second random effects term.
Best,
Wolfgang
-----Original Message-----
From: R-sig-meta-analysis [mailto:r-sig-meta-analysis-bounces at r-project.org] On Behalf Of Gram, Gil (IITA)
Sent: Thursday, 29 August, 2019 9:57
To: r-sig-meta-analysis at r-project.org
Subject: [R-meta] setting certain covariances to 0 for the random LHS term
Dear Wolfgang and others,
In my current model I have the following random structures: (treatment | site) and (treatment | site.time), i.e. random variance across site and across time-within-site for the different treatments.
Treatments here are defined as (1) control, (2) mineral input MR, (3) organic input OR, and (4) combined mineral and organic input ORMR.
Now I extend the number of treatments by the organic input sources, i.e. control, MR, ORone, ORtwo, ORthree, ORMRone, ORMRtwo, ORMRthree etc.
I wish to know the variances for each of these extended treatments, but I don?t want to see covariances between combinations of organic input sources, only with Control (and MR). Is it possible to let the model know that, by setting these to 0 somehow?
Thanks!
Gil Gram
PhD researcher | Natural resource management/CCAFS
International Institute of Tropical Agriculture (IITA)
East African Hub/Kampala/Country office
Address: IITA-Uganda Plot 15B, Naguru East Road. P.O Box 7878, Kampala
Mobile: +256 755 315236 | Belgium Mobile: +32 484 981200
Skype: gil.gram