[R-sig-eco] partial RDA as MANOVA: appropriate permutations for a complex multi-level experimental design

[Sybille Haeussler]

Hello everyone:


I am using partial RDA to carry out a MANOVA of community data (plants, mycorrhizal fungi, beetles) in vegan for a 22-yr old forestry experiment (part of the North American Long Term Soil Productivity Study) with a multi-level randomized block design with 1 random factor (BLOCK = R) and 3 fixed factors ( A, B and C).  Combinations of factors A & B (with three levels each) are completely randomized on 9 plots within each block.  Factor C (with two levels) are randomized on one half of each of the 9 plots per block.  This is a split plot design.   I am analysing just 1 set of observations on each split plot.  The experiment is fully balanced and there are no missing data.  Thus 54 total observations.  Note that the split-plot observations could be analysed in the same way as a time series with 2 observations per plot.


Here is the M/ANOVA table:

Source of Variation     df or q SS=variance     MS =SS/df       F or F# df for F
Total   (N – 1) = 53    SSTot   MSTot
Block   (r – 1) = 2     SSR     MSR
Factor A        (a – 1) = 2     SSA     MSA     MSA/MSEAB       2,16
Factor B        (b – 1) = 2     SSB     MSB     MSB/MSEAB       2,16
A x B interaction       (a-1)(b-1) = 4  SSAB    MSAB    MSAB/MSEAB      4,16
Plot-level Error AB     (ab-1)(r-1) = 16        SSEAB   MSEAB

Factor C        (c-1) = 1       SSC     MSC     MSC/ MSEABC     1,18
OM x SP interaction     (a-1)(c-1) = 2  SSAC    MSAC    MSAC/ MSEABC    2, 18
C x SP interaction      (b-1)(c-1 )= 2  SSAB    MSAB    MSAB/ MSEABC    2, 18
OMxCxSP interaction     (a-1)(b-1)(c-1) = 4     SSABC   MSABC   MSABC/ MSEABC   4, 18
Split-Plot  Error ABC   ab(c-1)(r-1) = 18       SSEABC  MSEABC



Tests of split-plot level factors and their interactions (anything including Factor C) should have a Residual (denominator term) with 18 degrees of freedom.  Tests of plot level factors and their interactions (A, B, AxB) should have a Residual with 16 degrees of freedom.  I use this as a guide to see whether the permutation test is using the correct F statistic.


Following the methods for partial RDA as MANOVA developed by Pierre Legendre, Marti Anderson & colleagues and using Borcard et al. (2011) as a guide for the R code,  I set up a matrix of orthogonal dummy vectors that sum to zero (Helmert contrasts).  The number of vectors for each factor & interaction term matches the df listed above.


To test each interaction and main effect I will put the vectors for the factor to be tested in the explanatory (X) matrix, and put all of the other vectors into the covariance matrix (W).  The response matrix (Y) is the Hellinger-transformed species abundance values.  I first tested for the homogeneity of variance-covariance matrices using betadisper().


I have read Anderson (2001) and concluded that I need to permute the residuals rather than the exchangeable units for 3 reasons: (1) an exact test is not possible for interaction terms, (2) there are too few possible permutations at the split-plot level , (3) permuting residuals gives higher power than an exact test for factorial experiments.  On reviewing vegan supporting information I came to understand that vegan automatically permutes residuals with partial RDA.


I also have carefully read Simpson’s guide to the permutation tests and following his guidelines concluded that the proper design for my permutation tests should be as follows.  The data are entered in a single column as split-plots (1 or 2) within plots (1-9) within blocks (1-3).


For plot-level permutations:  I want to average the data from the two “samples” per plot and freely permute the 9 plots within each of the 3 blocks).  With 9 plots the set of all permutations should be very large:

R> anova (FACTORA.rda, permutations= how(nperm=999, within=Within(type="none"), plots = Plots(strata = env$PLOT), blocks = env$BLOCK))



I received the following error message:

'nperm' > set of all permutations; Resetting 'nperm'.

Set of permutations < 'minperm'. Generating entire set.

Error in doAllPerms(spl[[i]], strataP, typeW, typeP, mirrorW, mirrorP,  :

  object 'res' not found



For split-plot permutations:

R> anova (FACTORC.rda, permutations= how(within=Within(type="series"), plots = Plots(strata = env$PLOT), blocks = env$BLOCK))

I expected to get an error message saying that the number of permutation was too small, given only two samples within the “series” to permute.  Instead I got a result but the F# was incorrect because the denominator had df=34 rather than df=18.  Using type="free" generates the same F# but different p-value.


QUESTIONS:

1.  With partial RDA, assuming that vegan permutes residuals after all of the effects of all of the conditioning variables in W have been removed, can I just freely permute the residuals, permutations  =how(nperm = 999), and not worry about the complicated multi-level hierarchy of my experimental design?  Is it redundant to use the coding of Simpson for partial RDA?  If so, do the degrees of freedom matter? Borcard et al. (2011, p. 187) use restricted permutations for their simpler example of partial RDA as MANOVA for a 2-factor CRD experiment.


2.  If (1) is incorrect, how should I set up my how() statements to get the correct degrees of freedom and sufficient permutations ?  I have tried many possibilities other than those shown, but clearly don't understand it.


3. If there are not sufficient exchangeable units to permute at the split-plot level, what is the best work-around?  I am able to generate correct df at the plot level by using a 27-line dataset that pools the data from the two split-plots.


References:

Anderson, M.J. 2001. Permutation tests for univariate or multivariate analysis of variance and regression.  Can. J. Fish. Aquat. Sci. 58: 626-639.


Borcard, D., F. Gillet and P. Legendre. 2011. Numerical Ecology with R. Springer, NY.  (pages 185-188).


Simpson, G.L. (undated). Restricted permutations: using the permute package. https://cran.r-project.org/web/packages/permute/vignettes/permutations.pdf


Thanks for any help you can provide.

Sybille ?



Sybille Haeussler PhD RPF
2041 Monckton Rd.
Smithers, BC V0J 2N4
phone: 250 847-6082