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

[Eduard Szöcs]

Also have a look at the new permute package:
http://cran.r-project.org/web/packages/permute/index.html

AFAIK permute is not fully hooked up in vegan at the moment (at least 
not in adonis),
but you could hack the adonis function to use shuffle() instead of 
permuted.index().

The permutation design, as proposed by Kay, would then look like this:

ctrl <- permControl(strata = Site,
                     within = Within(type = "none"),
                     blocks = Blocks(type = "free"))
table(Site, perm.Sediment = Sediment[shuffle(nobs, control = ctrl)])


HTH,

Edi

  




On 05/03/12 11:45, Kay Cichini wrote:
> Hi Arnaud and list-members,
>
> what about the following:
>
> library(vegan)
>
> Sediment<-as.factor(rep(c("Sed1","Sed2","Sed3"),each=36))
> Site<-as.factor(rep(c("Site1","Site2","Site3","Site4","Site5","Site6","Site7","Site8","Site9"),each=12))
> Hydrology<-as.factor(rep(rep(c("Dry","Wet"),each=6),9))
> Depth<-as.factor(rep(rep(c("D1","D2"),each=3),18))
>
> nobs = length(Site)
>
> # depth and hydrology are nested within sites:
> table(Site, Depth)
> table(Site, Hydrology)
>
> # so the following call is ok as it restricts permutation on sites -
> # that is, levels of depth and hydrology are shuffled only within and not
> # across sites:
> set.seed(123)
> community<- matrix(rnorm(nobs * 3, 10, 1), nobs, 3)
> adonis(community ~ Hydrology * Depth, perm = 999, strata = Site)
>
> # sediment is different:
> table(Site, Sediment)
>
> # to allow for this structure you would need to shuffle sites as
> # a whole. the below restriction will randomly assign levels of sediment
> # to the sites - like this:
> ctrl = permControl(strata = Site, permute.strata = TRUE)
> table(Site, perm.Sediment = Sediment[permuted.index2(nobs, control = ctrl)])
>
> ### to test sediment you'd need a customized code:
> (fit<- adonis(community ~ Sediment, 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<- permuted.index2(nobs, control = ctrl)
>       fit.rand<- adonis(community ~ Sediment[idx], permutations = 1)
>       pop[i]<- fit.rand$aov.tab[1,5]
> }
>
> ### get the p-value (H0: Sediment has no effect on location of communities):
> print(pval<- sum(pop>= pop[1]) / (B + 1))
>
> ### make a histogram to see random R2-values and the true one:
> hist(pop, xlab = "Population R2")
> abline(v = pop[1], col = 2, lty = 3)
>
> What do you think?
> Kay
>
>
> 2012/3/4 Arnaud Foulquier<arnaud.foulquier at gmail.com>
>
>> I have run the two analysis. The only differences are the p-values that
>> appear in the anova tables. Other values (Df, SumsOfSqs, MeanSqs, F.Model)
>> are not modified because they are calculated on the "original" (not
>> permuted) data.
>>
>> The p-value is calculated by comparing the value of F obtained with the
>> original data to the distribution obtained by permutations. If i understand
>> well, the strata argument is used to constrain the permutations of the
>> original data within the groups specified in the strata argument. The
>> distribution of the F-statistic that is generated by permutations is not
>> the
>> same according to whether the strata argument is specified or not. This is
>> why the p-values are modified when the strata argument is specified.
>>
>> However, I'm not sure that strata=Site_Hydro correctly specify the nested
>> structure of my data.
>>
>>
>>
>>
>>
>> --
>> View this message in context:
>> http://r-sig-ecology.471788.n2.nabble.com/split-split-plot-design-and-adonis-tp7332029p7342690.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
>>
> 	[[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-ecology mailing list
> R-sig-ecology at r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-ecology
>