question regarding agricolae package for split-plot design with significant main-plot x sub-plot interaction

947 views
Skip to first unread message

Maegen Simmonds

unread,
Aug 20, 2015, 7:21:24 PM8/20/15
to davi...@googlegroups.com
Hello!

I've been using the agricolae package for ANOVA and LSD tests for data from split-plot designs with randomized complete block design of main-plots. I would greatly appreciate any advice on how to analyze simple effects when there's a significant main-plot x sub-plot interaction. In the R tutorial (https://cran.r-project.org/web/packages/agricolae/vignettes/tutorial.pdf) on p 36-37, a multiple mean comparisons with a factorial treatment structure is performed for a model with a significant interaction (clone x nitrogen):

outAOV <-aov (yield ~ block + clone * nitrogen, data = A)
# this does LSD on the "simple effects" as opposed to "main effects" of clone and nitrogen separately:
outFactorial <-LSD.test (outAOV, c("clone", "nitrogen"), + main = "Yield ~ block + nitrogen + clone + clone:nitrogen",console=TRUE) 

My questions: 1. Does this also apply to split-plot designs with RCBD mainplots? (see p 8 for a visual example of this general layout, https://stat.ethz.ch/education/semesters/as2011/anova/Slid10.pdf) 2. Do I have the correct error term specified for the mainplot in Error() below, or is it supposed to be Error(block:mainplot)?

# General model structure for split-plot with RCBD mainplots
outAOV <-aov (yield ~ block + mainplot * subplot + Error(mainplot:block), data = my.data)  
# When significant interaction, analyze simple effects 
# However, the following produces an error message, and I think it's because of the error term. I could take it out, but I'm not sure that's correct 
outSplitPlot <-LSD.test (outAOV, c("mainplot", "subplot"), console=TRUE)

I could also subset the data for each level of mainplot and analyze the subplot effects (and vice versa), but that reduces the power of the test. 

Thanks! 

-Maegen

--
Maegen Simmonds, Ph.D.

Guillaume Théroux Rancourt

unread,
Aug 20, 2015, 8:02:40 PM8/20/15
to davi...@googlegroups.com
Hi Maegen,

A split-plot design should be analyzed as a mixed model with your main plot and sub-plots in the random effects. There is an issue with getting the most reliable estimates when using only a aov() or lm(), especially when there is some special blocking like in a split-plot.

Here is an example from the nlme library I think (it's from class notes):

library(nlme)
data(Oats) #yield of oats as a function of 3 varieties and nitrogen concentration
Oats$nitro<-as.factor(Oats$nitro)
mod1<-lme(yield~nitro+Variety+nitro:Variety, random=~1 | Block/Variety, data=Oats)  #Variety is the sub-plot
summary(mod1)
anova(mod1)

There are several way to get both you means and SE, as well as your multiple comparisons. For the model above, you could rewrite the model to keep nitro because it is the only significant effect. Then, using the multcomp package, you could run something like that:
summary(glht(mod1, linfct = mcp(nitro = "Tukey")))
which also gives you Tukey corrected P value for the comparison, so you can then reconstruct the letters (a, b ...) that you would get in agricolae.

You can also get your means and SE with the AICcmodavg package, function predictSE (see the help for examples).

For more on mixed model and the lme package, see the book by Pinheiro and Bates (2000), which is probably available as a PDF from the library.

HTH,

Guillaume


--
Check out our R resources at http://www.noamross.net/davis-r-users-group.html
---
You received this message because you are subscribed to the Google Groups "Davis R Users' Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email to davis-rug+...@googlegroups.com.
Visit this group at http://groups.google.com/group/davis-rug.
For more options, visit https://groups.google.com/d/optout.

Maegen Simmonds

unread,
Aug 20, 2015, 8:22:00 PM8/20/15
to davi...@googlegroups.com
Thanks for your reply, Guillaume! 

Since the mainplot and subplots are fixed treatments, I'm not choosing to have them set as random. I'm interested in their effects. If the experiment was performed at multiple sites, then I would do a mixed effects model with site as random and block nested within site as random, as well. However, this is a simple scenario at one site, split-plot design. I think a fixed effect model would be valid. In the example you gave, I think block was the only random variable specified. That makes sense to me, but I think it's correct to have block as fixed as well.

My main question is how do I analyze the simple effects in a split-plot design with RCBD main-plot using aov() and lsd.test()? I think it's possible, but I could totally be wrong!

Thanks,
Maegen

Guillaume Théroux Rancourt

unread,
Aug 21, 2015, 2:43:47 PM8/21/15
to davi...@googlegroups.com
Hi Maegen,

Just to get over the mixed model question first, I have analyzed split-plot only using lme() because I had missing data. aov() and lm() do not produce reliable estimates when missing data is present (actually, the help file of aov() say that aov() and lm() are statistically inefficient when there is 2 or more error strata, i.e. a split-plot, and the data is unbalanced. Also, a previous mentor and professor responsible for the experimental design course told me that some agricultural science journals where accepting only mixed model for split-plot (although I do not know those journal and should check). So, this is why I'm kind of stuck to my lme(). ;) In the split-plot example I provided, Variety was both a fixed and a random effect, so you can have them both in your mixed model.

Even if I'm biased toward lme(), this page seems to give the info you need with the statistical explanation and R code: http://stats.stackexchange.com/questions/13788/split-plot-in-r.

Your Error term should have been
Error(mainplot/block)
if block is within the mainplot. This is not clear based on your formula in. I would guess that you actually mean:
Error(block/mainplot)
which would mean that your main plots are within your blocks, thus that "block > mainplot > subplot"

HTH,

Guillaume


Tim Bowles

unread,
Aug 21, 2015, 4:48:39 PM8/21/15
to davi...@googlegroups.com
Hi Maegen,

I would also recommend running your split plot as a mixed effects model. In my experience, using packages lme4 to fit the model, lmerTest for hypothesis tests of fixed factors with adjusted degrees of freedom (Satterthwaite or Kenward-Roger), and lsmeans for post-hoc hypothesis tests (like LSD), will work well.

For instance, for a split plot with the main plot as factor A and the subplot as factor B:
library(lme4)
model1 <- lmer(responsevariable ~ A*B + (1|block) + (1|block:A), data=yourdata)

library(lmerTest)
anova(model1, ddf='Kenward-Roger')

library(lsmeans)
cld(lsmeans(model1, tukey ~ A + B)) # for Tukey test


Best,
Tim



PhD Candidate
Jackson Lab
Graduate Group in Ecology, UC Davis
Reply all
Reply to author
Forward
0 new messages