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

[s.graham]

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-adonis-tp7332029p7578101.html
Sent from the r-sig-ecology mailing list archive at Nabble.com.