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

[Gil Gram]

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).


Questions/remarks

Does ‘treatment' need to be 100% balanced? Can we include some studies with no control treatment for instance?
I do expect differences over time for some treatments, but since ‘time’ is so unbalanced I cannot use it as a fixed effect. 
Why does Metafor only allow for ≤ 2 random elements?
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.?
hyp2: I look at the fixed effects only, and the (significant) differences between the different levels of treatExt?
hyp3: where in the output can I find the variance according to the different levels of treatExt? (see below an example of my output)
hyp4: I look at the predicted yield values (fixed+ran) and plot them in boxplots, according to each level of the secondary parameters separately?

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(): does it mean that the model results estimates we see in summary() are the fixed effect + random 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


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]]