[R-meta] setting certain covariances to 0 for the random LHS term

[Viechtbauer, Wolfgang (SP)]

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 using r-project.org] On Behalf Of Gram, Gil (IITA)
Sent: Thursday, 05 September, 2019 8:20
To: r-sig-meta-analysis using 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 using maastrichtuniversity.nl<mailto:wolfgang.viechtbauer using 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 using r-project.org] On Behalf Of Gram, Gil (IITA)
Sent: Thursday, 29 August, 2019 9:57
To: r-sig-meta-analysis using 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