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

[Viechtbauer, Wolfgang (SP)]

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)