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. _______________________________________________ R-sig-ecology mailing list R-sig-ecology@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-ecology