[R-sig-eco] split-split-plot design and adonis

[Steve Brewer]

Sharon,

I've had difficulty trying to figure out how to do a split-split-plot
anova on bray-curtis distances in R and thus may not be much help to you.
If you didn't have the habitat effect and associated subplots, you could
do a simple split-plot analysis using two separate analyses (two-way
nested.npmanova with plot nested within treatment for the whole-plot test
and adonis and strata for the split-plot test: species~treatment*shade +
plot; strata = plot for the split-plot effects). See my February 25 2013
message entitled "Permanova with anova." It was intended for
repeated-measures but applies to split-plots as well.

This approach can't be used for split-split analysis, however, because
nested.npmanova won't handle two levels of nesting.

What we need is a way for R to calculate distances among centroids given a
particular grouping variable and then make a distance matrix from those.

For example, in your case, you need to be able to calculate distances
among centroids among each treatment by shade level combination (perhaps
call it treatshade, with each combination having a unique label), and then
create a matrix from those centroids. You could then use nested.npmanova
to test the main effect of treatment invert.hel ~ treatment + plot; you
ignore the test of plot. You could then analyze the effects of shade and
treatment:shade, while ignoring the test for treatment using invert.hel ~
treatment*shade + plot, strata = plot. Next, for the split-split plot, you
would analyze the original distance matrix (I.e., distance among
subplots). In this case, however, you would need to create a new
categorical predictor variable in which each treatment by shade
combination gets a unique label (treatshade), but this time it will be
repeated for each habitat and subplot. Make the following model -
invert.hel ~ treatshade + habitat + habitat:shade + habitat:treatment +
habitat:treatment:shade + plot, strata = plot. In this case, you can
evaluate the habitat effect, the habitat*shade interaction, the
habitat*treatment interaction, and the three-way interaction, and ignore
the treatshade effect.

Sorry I can't be of more help. Maybe somebody on the listserv knows how to
calculate distances among centroids (or for that matter calculate
Bray-curtis centroids). That would help a lot.


J. Stephen Brewer 
Professor 
Department of Biology
PO Box 1848
 University of Mississippi
University, Mississippi 38677-1848
 Brewer web page - http://home.olemiss.edu/~jbrewer/
FAX - 662-915-5144
Phone - 662-915-1077




On 4/26/13 1:40 AM, "s.graham" <sharon.graham at pg.canterbury.ac.nz> wrote:

>Hello all, 
>
>I've been trying to adapt the code suggested by Kay to a similar
>split-split-plot experiment.  My experimental design was as follows:
>
>6 plots, 3 of which received "treatment" and 3 which were "control"
>Each plot was split into "shade" and "no shade"
>within each shade/no shade there were 6 subplots, 3 of which had "habitat"
>and 3 which did not
>
>So in total I had 72 sub-plots (3 reps of each possible combination),
>hierarchically arranged: Treatment/Shade/Habitat.  I have done mixed
>effects
>models on measures of community composition, such as species richness, etc
>and found that Treatment, Shade, Treatment:Shade (interaction) and
>Treatment:Leaves were all significant predictors of richness.  Now I would
>like to test the effects of my experimental combinations
>(Trt/Shade/Habitat)
>and interactions on community composition.
>
>Following Steve's comments and Kay's suggestions to Arnaud, I set up my
>code
>as the following: 
>
># To test for Treatment:Shade
>adonis(invert.hel ~ Treatment*Shade*Habitat+Plot, method="bray",
>strata=Plot, perm=999)
>
># To test for Shade:Habitat and Treatment:Habitat
>Plot_Shade <- interaction(Plot,Shade)
>adonis(invert.hel ~ Treatment*Shade*Habitat+Plot_Shade, method="bray",
>strata=Plot_Shade, perm=999)
>
># To test for Shade, Habitat main effects and Shade:Habitat
>set.seed(123)
>adonis(invert.hel ~ Shade * Habitat, perm = 999, strata = Plot)
>** NOTE: I am a little confused whether this accounts for Habitat being
>nested within Shade?
>
># To test Treatment main effect:
>nobs = length(Plot)
>ctrl <- permControl(strata = Plot, within = Within(type = "free"))
>(fit <- adonis(invert.hel ~ Treatment, permutations = 1))
>
>### number of perms
>B <- 999
>
>### setting up frame which will be populated by random r2 values:
>pop <- rep(NA, B + 1)
>
>### the first entry will be the true r2:
>pop[1] <- fit$aov.tab[1, 5]
>
>for(i in 2:(B+1)){
>     idx <- shuffle(nobs, control = ctrl)
>     fit.rand <- adonis(invert.hel ~ Treatment[idx], permutations = 1)
>     pop[i] <- fit.rand$aov.tab[1,5]
>}
>
>### get the p-value (H0: Treatment has no effect on location of
>communities):
>print(pval <- sum(pop >= pop[1]) / (B + 1))
>
>** Here, I get a p-value of 1, but Treatment has the highest R2 of all
>variables, and others with lower R2 (Shade, Habitat) came out as
>significant.  Treatment was also highly significant in my mixed effects
>models.  So I'm confused by that.
>
>When I look at the "pop" dataframe generated by this code, values are
>identical in every cell - I think this might be the problem?  Any ideas
>why
>that is happening?  Arnaud, did you ever use this code on your data and if
>so did you encounter the same problem?
>
>Thanks in advance for your consideration.
>
>Cheers,
>Sharon
>
>
>
>
>
>--
>View this message in context:
>http://r-sig-ecology.471788.n2.nabble.com/split-split-plot-design-and-adon
>is-tp7332029p7578101.html
>Sent from the r-sig-ecology mailing list archive at Nabble.com.
>
>_______________________________________________
>R-sig-ecology mailing list
>R-sig-ecology at r-project.org
>https://stat.ethz.ch/mailman/listinfo/r-sig-ecology