[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

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