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

[Viechtbauer, Wolfgang (SP)]

Hi Gil,

How about using struct="HCS"? That will give you different tau^2 values for each treatment, but assumes that there is a single correlation parameter. See help(rma.mv) and search for "heteroscedastic compound symmetric structure". That might be more feasible with your data, even when distinguishing the sub-treatments (since the number of correlation parameters doesn't 'explode' then).

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 14:04
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,

Indeed, let me try to explain briefly but better;

We have ±2500 maize yield data under 4 different treatments, Control (no input), OR (organic input), MR (mineral input), and ORMR (combined organic and mineral input). However, we also have sub-treatments as for all the organic treatments (OR and ORMR) we have different 4 different kinds of organic input, i.e. class 1, 2, 3, and Manure.

We are interested in (1) the effect of treatment on yield, but also in (2) the effect of treatment on yield variability (variance across each treatment). We were hoping to get the latter from the random variance components, but as you have seen, quite impossible to get the variances at sub-treatment level. So I’m trying to find alternative ways, and one I considered was having the full model for question (1) but separate models for question (2) where I run the model on OR type subsets of the data.

I am aware that ideally we want to stick to one model, but maybe for what I need this is OK? Of course I’m interested in any other suggestions!

thanks,
Gil

> On 5 Sep 2019, at 14:40, Viechtbauer, Wolfgang (SP) <wolfgang.viechtbauer using maastrichtuniversity.nl> wrote:
> 
> Hi Gil,
> 
> I cannot answer the first question ("Does that seem OK to you?") without more details.
> 
> As for the complexity of the model, it has (as far as I can tell), 6 (rho) + 6 (phi) + 4 (tau2) + 4 (gamma2) + 2 (sigma2) = 22 parameters for the random effects structure. Still a lot, but might be feasible with 2500 data points. Of course, that alone isn't sufficient to determine whether this is the case; it also matter how those data points are distributed over the various levels and combinations.
> 
> As for struct="UN" -- that is fine. It automatically gets expanded to struct=c("UN","UN") inside of the function anyway.
> 
> 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 11:00
> 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 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