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

[Gil Gram]

Dear Wolfgang and all,

I realized the inserted output was not shown in my mail below. Please
consider my adapted mail, and if not clear, please let me know:

On Fri, Apr 13, 2018 at 5:00 PM, Gil Gram <gilgram.fm at gmail.com> wrote:

> 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"))
>
> (first screenshot)
https://we.tl/uj8qClWmrm


>
> 2) Then for hyp3 (do treatments and classOR influence the yield
> variability over time?), where in the output link 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?
>
>
> (second screenshot)
https://we.tl/uj8qClWmrm


>
>
>
> 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] On Behalf Of Gil Gram
> Sent: Thursday, 29 March, 2018 9:21
> To: 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>() 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>()
> 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)
>
>
>

	[[alternative HTML version deleted]]