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

[Gram, Gil (IITA)]

Hi Wolfgang,

Thanks again. I understand the issue, and I was indeed afraid the model would be getting far too complex.

My alternative was to fit separate models for each class of organic resource (OR), only in order to get their individual variances. Does that seem OK to you?

And lastly, my full model would then look like this, with treatment having 4 levels and 2500 data points. As far as you can tell, not too complex?

res = rma.mv(sqrt(yi), vi, method = 'REML', struct=c(“UN”, “UN”), sparse=TRUE, data=dat,
                         mods = ~ rateORone + rateORtwo + rateORthree + rateORManure + kgMN
                         + I(rateORone^2) + I(rateORtwo^2) + I(rateORthree^2) + I(rateORManure^2) + I(kgMN^2)
                         + rateORone:kgMN + rateORtwo:kgMN + rateORthree:kgMN + rateORManure:kgMN
                         + I(rateORone^2):I(kgMN^2) + I(rateORtwo^2):I(kgMN^2) + I(rateORthree^2):I(kgMN^2) + I(rateORManure^2):I(kgMN^2)
                         + cropSys + idF,
                         random = list(~1|ref, ~1|idRow, ~ treatment | idSite, ~ treatment | idSite.time))

I was indeed only using struct=“UN” instead of c(“UN”, “UN”), will change that, thanks!

Appreciate your help, as always.
Gil


On 5 Sep 2019, at 11:36, Viechtbauer, Wolfgang (SP) <wolfgang.viechtbauer using maastrichtuniversity.nl<mailto:wolfgang.viechtbauer using maastrichtuniversity.nl>> wrote:

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<mailto: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><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<mailto: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


	[[alternative HTML version deleted]]