[R-meta] interpreting and using the rma.mv output for different questions

Fri Apr 13 16:00:52 CEST 2018

Thank you Wolfgang for your time and effort,

My rma.mv model syntax is indeed: mod = rma.mv(yield, SEfromVARCOM, mods = ~ 0+treatExt, random = list(~ classOR|ref, ~ classOR|site/time), data = dtProClass, method="ML")

- Concerning the unbalancedness of time for fixed effects, please explain me if and why I’m wrong…

Time has 2 features, which I would think refer to random and fixed effects. Time can be a random variable because the different growing seasons will never be identical, i.e. at time 2 it can rain more than at time 3. The 4 different treatments at time 2 could be correlated. Moreover, the weather conditions at time 4 at Mozambique will also differ from time 4 at Senegal. Time is thus taking into account the weather conditions within a site (site/time).

Time can be regarded as a fixed effect, in the sense that for the organic input treatment, the build op of organic matter in the soil is hypothesised to increase the crop’s performances on the long run. So we expect higher yields at time 10 than at time 1 for the organic treatments (not for Control or mineral treatment). However, I’ve been told that although it is an interesting question, the data wasn’t designed for it and putting time as a fixed effect would be problematic. Because time has over 25 levels that are highly unbalanced, the estimates of the few long-term studies will suffer from high variance.

- Concerning my research questions… 

1) If I understand correctly, for hyp1 (differences between treatments) and hyp2 (differences between organic quality classes classOR) we should use the fixed effect estimates from the model’s output. But in the output the treatments (control, OR, MR, ORMR are mixed with the classOR levels (one, two, … Manure). So how do I disaggregate them or do a pairwise comparison on sets of predictor levels to test both hypotheses?

For pairwise comparison of each of the levels I know can use: summary(glht(mod, linfct=cbind(contrMat(rep(1,12), type = “Tukey”))), test = adjusted("none"))

2) Then for hyp3 (do treatments and classOR influence the yield variability over time?), where in the output above can I find the variance over site/time for each treatment separate? I expect to see a smaller variance for the ORMR treatments, but to different extents depending on organic matter quality classOR…

3) Finally, if I want to test the influences of for example soil texture classes 1 and 2, which are not included in the model, should I then compare the BLUP’s of each of the texture classes’ estimates? i.e. for every data observation in the dataset I compute its fixed effect estimate + ranef$classOR + ranef$classOR/site + ranef$classOR/site/time + ranef$classOR/ref.
Is there a function that does this automatically or does it need to be done manually?
Or how would you test if texture has a significant effect on the treatment effects on yield?

Thank you so much…

Best regards,
Gil

> On 11 Apr 2018, at 23:56, Viechtbauer, Wolfgang (SP) <wolfgang.viechtbauer at maastrichtuniversity.nl> wrote:
> 
> Dear Gil,
> 
> See some responses below.
> 
> Best,
> Wolfgang
> 
>> -----Original Message-----
>> From: R-sig-meta-analysis [mailto:r-sig-meta-analysis-bounces at r-
>> project.org <http://project.org/>] On Behalf Of Gil Gram
>> Sent: Thursday, 29 March, 2018 9:21
>> To: R-sig-meta-analysis at r-project.org <mailto:R-sig-meta-analysis at r-project.org>
>> Subject: [R-meta] interpreting and using the rma.mv output for different
>> questions
>> 
>> Hi all,
>> 
>> I’m seeking feedback on the way I’m trying to perform my meta-analysis
>> with metafor. I’ve been trying to find out by myself, but guess now I
>> need ‘personal’ advice to get me going again.
>> 
>> Data
>> 
>> The studies are about agronomic experiments with 4 treatments (control,
>> organic input OR, mineral input MR and the combined input ORMR) on maize
>> yield.
>> 
>> Primary parameters:
>> ref (48 references)
>> site (76)
>> yield (5020 observations) + SE
>> time (the growing season in the experiment —> unbalanced, with many
>> short-term trials, and few long-term)
>> classOR (6 quality classes of OR: zero (control and MR treatment), I, II,
>> III, IV, V —> slightly unbalanced)
>> treatment (4 levels —> 100% balanced)
>> —> treatExt (treatment extended = treatment and classOR combined, 8
>> levels: Control, MR, ORclass1, ORclass2, … ORMRclass1, ORMRclass2, …)
>> 
>> Secundary parameters (categorical):
>> soil texture, aridity, application rates, …
>> 
>> 
>> Hypotheses & Model
>> 
>> The three hypotheses to be tested are:
>> Does ORMR have significant higher yields than the other treatments?
>> Are these yield effects dependent on the quality of the organic N
>> resource (classOR)?
>> Are these yield increases more reliable (ie less variable in space and
>> time)?
>> Do the secondary parameters influence the effect
>> 
>> The multivariate model I came up with, to run with rma.mv:
>> 
>> yield ~ treatExt + (classOR | ref) + (classOR | site/time)
>> 
>> with treatExt as only fixed factor, and random intercepts for reference
>> (study variance), site (spatial variance), and time within site (temporal
>> variance).
> 
> This isn't 'metafor' syntax, but maybe that wasn't your intention.
> 
>> Questions/remarks
>> 
>> Does ‘treatment' need to be 100% balanced? Can we include some studies
>> with no control treatment for instance?
> 
> That shouldn't be a problem.
> 
>> I do expect differences over time for some treatments, but since ‘time’
>> is so unbalanced I cannot use it as a fixed effect.
> 
> Not sure why you think that.
> 
>> Why does Metafor only allow for ≤ 2 random elements?
> 
> rma.mv() allows for an unlimited number of random effects. But it only allows for two terms of the form '~ inner | outer'.
> 
>> What different model output information can I use to test the different
>> hypotheses?
>> hyp1: I look at the fixed effects only, but have to somehow disaggregate
>> the classOR from treatExt, so by doing an omnibus moderator test for OR
>> treatments separate etc.?
> 
> Sorry, I don't understand.
> 
>> hyp2: I look at the fixed effects only, and the (significant) differences
>> between the different levels of treatExt?
> 
> Sorry, I don't understand.
> 
>> hyp3: where in the output can I find the variance according to the
>> different levels of treatExt? (see below an example of my output)
> 
> Sorry, I don't understand. Also, I see no example output below.
> 
>> hyp4: I look at the predicted yield values (fixed+ran) and plot them in
>> boxplots, according to each level of the secondary parameters separately?
> 
> Sorry, I don't understand.
> 
>> Trying to figure question 4 out, made me confused about the following:
>> 
>> For lmer; predicted yield (predict())or model coefficient (coef()) =
>> fixed effect estimate (fitted() or summary() output, BLUE) + random
>> effect predictions (ranef(), BLUP)
>> 
>> For rma.mv;, however, I read this:
>> 
>> ranef <https://rdrr.io/cran/metafor/man/ranef.html <https://rdrr.io/cran/metafor/man/ranef.html>>() documentation:
>> predict(), coef(), fitted(), and the estimates of the summary() output,
>> produce the same estimates which are based only on the fixed effects of
>> the model
>> ranef() calculates best linear unbiased predictions (BLUPs) of only the
>> random effects
>> 
>> But from the blup <https://rdrr.io/cran/metafor/man/blup.html <https://rdrr.io/cran/metafor/man/blup.html>>()
>> documentation:
>> blup() for BLUPs that combine the fitted values based on the fixed
>> effects and the estimated contributions of the random effects (not for
>> rma.mv)
>> 
>> In other words,
>> in rma.mv, predict() and coef() are the same as fitted(): 
> 
> That's not correct. coef() gives the model coefficients. predict() gives predicted values. Those are the same as fitted(), but with predict() one gets SE and CIs and one can specify for which combination of moderator values predicted values should be computed.
> 
>> does it mean
>> that the model results estimates we see in summary() are the fixed effect
>> + random effects? 
> 
> No. Those are the fixed effects.
> 
>> And so if I want to know the fixed effects only, I
>> should subtract them with the random effects (fitted - ranef)?
>> you have two definitions for BLUM? (BLUP’s of random effects, and BLUP’s
>> of fixed and random effects combined
> 
> ranef() gives the BLUPs of the random effects only. blup() gives fitted values + random effects.
> 
>> Sorry if this seems much to ask, I’ve tried to compile all my questions
>> compactly…
>> 
>> Your advice would be genuinely appreciated!
>> 
>> Gil (PhD student at KUL/IITA)

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <https://stat.ethz.ch/pipermail/r-sig-meta-analysis/attachments/20180413/d61ce920/attachment-0001.html>

-------------- next part --------------
A non-text attachment was scrubbed...
Name: PastedGraphic-7.png
Type: image/png
Size: 593338 bytes
Desc: not available
URL: <https://stat.ethz.ch/pipermail/r-sig-meta-analysis/attachments/20180413/d61ce920/attachment-0002.png>

-------------- next part --------------
A non-text attachment was scrubbed...
Name: PastedGraphic-8.png
Type: image/png
Size: 224556 bytes
Desc: not available
URL: <https://stat.ethz.ch/pipermail/r-sig-meta-analysis/attachments/20180413/d61ce920/attachment-0003.png>