[R-meta] Intercept-slope model & network meta-analysis
Viechtbauer Wolfgang (SP)
wolfgang.viechtbauer at maastrichtuniversity.nl
Tue Aug 29 15:48:57 CEST 2017
Much more legible - thanks.
As for your first question: " Question: do I need to include the moderator variables in random effects? Is it enough to use the AIC to test the goodness of fit of the models and likelihood ratio of them to select the best model?"
The moderator variables are constant within 'study', so you cannot include them as random effects.
As for model selection, there are many different approaches and opinions one can take. You could use information criteria (like AIC) to select the model, but make sure you use REML=FALSE then (since you are comparing models with different fixed effects). Or you could fit the 'sev * mod_A' model, test the interaction (and report the results), if significant, report the results from that model, if not fit the 'sev + mod_A+' model and report that model.
You may also want to consider including block as a random effect. So: (sev|study/rep)
As for the second part: I would stick to just analyzing the raw data. I see no benefit here for aggregating and analyzing the means.
Wolfgang Viechtbauer, Ph.D., Statistician | Department of Psychiatry and
Neuropsychology | Maastricht University | P.O. Box 616 (VIJV1) | 6200 MD
Maastricht, The Netherlands | +31 (43) 388-4170 | http://www.wvbauer.com
From: R-sig-meta-analysis [mailto:r-sig-meta-analysis-bounces at r-project.org] On Behalf Of Juan Pablo Edwards Molina
Sent: Tuesday, August 29, 2017 15:25
To: r-sig-meta-analysis at r-project.org
Subject: [R-meta] Intercept-slope model & network meta-analysis
I have a datset containing 36 field plots experiments testing the effect
of several fungicides to control a soybean fungic disease.
This is how my raw data looks like (36 independent studies - CRBDs):
study fungic rep Mod_A Mod_B sev yield
1 Check 1 2High 1Low 55 2918
1 Check 2 2High 1Low 50 3468
1 Check 3 2High 1Low 45 1626
1 Check 4 2High 1Low 40 2921
1 Trt_A 1 2High 1Low 35 2414
1 Trt_A 2 2High 1Low 40 3104
1 Trt_A 3 2High 1Low 25 1878
1 Trt_A 4 2High 1Low 30 1952
1 Trt_B 1 2High 1Low 40 2708
1 Trt_B 2 2High 1Low 50 2475
At each study, a set of fungicides are the treatments including a
Check (different combinations across the studies, that´s why I adopted
network MA), "rep" are the blocks, "sev" is the disease (%) and
"yield" is the grain mass.
The moderator variables are study-specific characteristics, like disease
pressure (Mod_A) or Yield potential (Mod_B)
I have two objectives:
1° estimate the intercept and slope of the relationship yield ~ sev and
test the inclusion of moderator variables (I´m not testing the effect of
the treatments in this case, I´m interested on the trends of yield ~ sev).
I started using a multivariate Two-Stage Analysis approach then,
following the tutorial (
I moved into a multi-level Mixed-Effects Models with very similiar
results (but much more time-efficiency)
I am trying this:
# Overall random intercept and slopes
m1 <- lmer(yield ~ sev + (sev|study), data=df)
# Including effect of moderators on the intercept and slopes
m2 <- lmer(yield ~ sev * mod_A+ (sev|study), data=df)
# Including effect of moderator A on the intercept
m3 <- lmer(yield ~ sev + mod_A+ (sev|study), data=df)
# Including effect of moderator A on the slope
m4 <- lmer(yield ~ sev : mod_A+ (sev|study), data=df)
# Including effect of moderator A on the slope and moderator B on the
m5 <- lmer(yield ~ sev : mod_A + mod_B + (sev|study), data=longs)
Question: do I need to include the moderator variables in random effects?
Is it enough to use the AIC to test the goodness of fit of the models and
likelihood ratio of them to select the best model?
2° Then I do wanted to test the effect of treatments on yield, considering
mean differences to the untreated checks within each study.
So I performed a network meta-analysis, agreggating the data and estimating
the Mean Square Error from each study ANOVA:
study fungic yield_m Mod_A Mod_B MSE
1 Check 2640 2_High 1_Low 88931.95
1 Trt_A 2733 2_High 1_Low 88931.95
1 Trt_B 2858 2_High 1_Low 88931.95
where yield_m is the within-study treatment mean and MSE is the
within-study mean square error from ANOVA
The model I tried is:
net_D <- rma.mv(yield_m, vi2,
mods = ~ fungic * Mod_A,
random = ~ fungic | study,
struct= "UN", method="ML",
control = list(optimizer="nlm"))
anova(net_D, btt=9:14) # to test the effect of moderators
where vi2: vi = MSE / bk #Sampling variance for yi (bk = 4)
My concern is if Am I going well with this model? or should I try to use
the raw data as well, considering the block effect?
Thanks for your help!
(Phd candidate at Plant disease epidemiology lab in Univ. Sao Paulo -
More information about the R-sig-meta-analysis