Chapter 1: Permutational ANOVA and MANOVA (PERMANOVA) Key references: Method: Anderson (2001a), McArdle & Anderson (2001) Permutation techniques: Anderson (2001b), Anderson & ter Braak (2003) 1.1 General description Key references Method: Anderson (2001a) , McArdle & Anderson (2001) Permutation techniques: Anderson (2001b) , Anderson & ter Braak (2003) PERMANOVA is a routine for testing the simultaneous response of one or more variables to one or more factors in an analysis of variance (ANOVA) experimental design on the basis of any resemblance measure, using permutation methods. It is assumed that the user has relevant knowledge of multi-factorial ANOVA, which has the same basic logic in multivariate as in univariate analysis (see Underwood (1981) , Underwood (1997) and Quinn & Keough (2002) ), and an understanding of what it means to test a multivariate hypothesis (see Clarke (1993) ). A more complete description of the method is given in Anderson (2001a) and McArdle & Anderson (2001) . These papers merely elaborate the essential idea of partitioning for dissimilarity matrices which was originally (to our knowledge) presented by Brian McArdle ( McArdle (1990) , McArdle (1994) ) and which has also been articulated by Pillar & Orloci (1996) , Legendre & Anderson (1999) and Gower & Krzanowski (1999) . In essence, the routine performs a partitioning of the total sum of squares4 according to the full experimental design specified by the user, including appropriate treatment of factors that are fixed or random, crossed or nested (hierarchical), and all interaction terms. The routine will correctly calculate an appropriate distance-based pseudo-F statistic for each term in the model, based on the expectations of mean squares (EMS), in a fashion that is directly analogous to the construction of the F statistic for multi-factorial univariate ANOVA models ( Cornfield & Tukey (1956) , Hartley (1967) , Rao (1968) ). P-values are obtained using an appropriate permutation procedure for each term, and the user can specify whether permutation of raw data or residuals under either a full or reduced model are to be used ( Anderson (2001b) , Anderson & ter Braak (2003) ). Correct P-values may also be obtained through Monte Carlo random draws from the asymptotic permutation distribution ( Anderson & Robinson (2003) ) in the event that too few permutations are available for a given test. In addition to the main overall PERMANOVA partitioning and tests, the routine will also perform a posteriori pair-wise comparisons among levels of factors, including within individual levels of other factors (or cells) in the case of significant interaction terms. Other important features of the PERMANOVA routine include: (i) catering for unbalanced designs, including a choice regarding the type of sum of squares to be used for the partitioning; (ii) pooling or the exclusion of individual terms from a model; (iii) choice to include one or more quantitative covariates in the model; (iv) contrasts; (v) analysis of designs lacking replication; and (vi) analysis of asymmetric designs. Before using PERMANOVA, the user should have a fairly solid grasp of the basic issues and logic in experimental design (some good reference textbooks on this topic include Mead (1988) , Snedecor & Cochran (1989) , Winer, Brown & Michels (1991) , Underwood (1997) and Quinn & Keough (2002) ). Know your hypotheses and know your design! It is also best to store and import the factors and their levels along with the data in order to avoid making mistakes. Otherwise, use the tools available in PRIMER to either create or import factors and their levels (see chapter 2 in Clarke & Gorley (2006) ). 4 The total sum of squares is understood here to be defined by reference to the distance/dissimilarity measure of choice. It will only correspond to the traditional univariate total sum of squares when one variable is being analysed and Euclidean distance has been chosen as the basis of the analysis. 1.2 Partitioning We shall begin by considering the balanced one-way (single factor) ANOVA experimental design. A factor is defined as a categorical variable that identifies several groups, treatments or levels which we wish to compare. Imagine that we have one factor with a groups (or levels) and n observations (samples5) per group for a total of N = a × n samples. For each sample, we have recorded the values for each of p different variables. Recall that in univariate analysis of variance, the total sum of squares ($SS _ T$, the sum of squared deviations of observations from the overall mean) is partitioned into two parts that are meaningful for testing hypotheses about group differences: the within-group (or residual) sum of squares ($SS _ {Res}$, the sum of squared deviations of observations from their own group mean) and the among-group sum of squares ($SS _ A$, the sum of squared deviations of group means from the overall mean). A directly analogous partitioning is done in multivariate space by PERMANOVA. PERMANOVA may be thought of as a method that takes a geometrical approach to MANOVA ( Edgington (1995) ). Let each of the p variables be a dimension, and each of the N samples be represented by points in the p-dimensional space according to the values they take for each variable along each dimension. Now, the simplest of all multivariate systems has only p = 2 variables. It is good to consider this situation, because it is easy to draw just 2 dimensions! Now, imagine that there were n = 10 replicate samples in each of a = 3 groups, as shown in Fig. 1.1. Fig. 1.1. Plot of a hypothetical data set with p = 2 variables (dimensions) and n = 10 replicate samples in each of a = 3 groups. The three groups are identified by different symbols in the plot. The whole set of N = 30 samples taken together create a data cloud which is centred on a point called the overall centroid. For a Euclidean system, this is obtained as the arithmetic average for each of the (in this case 2) variables. Similarly, each of the groups also has its own group centroid, located in the centre of each of the clouds of points identified for each group. Fig. 1.2. Plots of the hypothetical data set from Fig. 1.1 showing the geometric partitioning. Just as in univariate ANOVA, we can consider the distance of any given point (sample) from the overall centroid in this space as being made up of two parts: the distance from the point to its group centroid (Fig. 1.2a) plus the distance from the group centroid to the overall centroid (Fig. 1.2b). This is the essence of the geometric approach to MANOVA in Euclidean space. We can calculate the sums of squares as: $SS _ T$ = the sum of squared distances from the samples to the overall centroid, $SS _ {Res}$ = the sum of squared distances from the samples to their own group centroid, and $SS _ A$ = the sum of squared distances from the group centroids to the overall centroid. The well-known univariate ANOVA identity: $SS _T = SS _ {Res} + SS _ A$ also holds for this geometric conception of MANOVA in Euclidean space. Verdonschot & ter Braak (1994) and Legendre & Anderson (1999) also remark on how these sums of squares are equal to the sum of the individual univariate sums of squares for each of the separate variables if Euclidean distance is used. 5 The word sample will be used throughout this manual in the manner that ecologists, and not statisticians, have come to understand the word. A sample shall mean a single unit used for sampling, such as a core, a transect, or a quadrat. This is consistent with the use of this word in PRIMER. 1.3 Huygens’ theorem This partitioning is fine and perfectly valid for Euclidean distances. But what happens if we wish to base the analysis on some other dissimilarity (or similarity) measure? This is important because Euclidean distance is generally regarded as inappropriate for analysing ecological species abundance data6, especially because of its lack of emphasis on species composition (e.g., Faith, Minchin & Belbin (1987) , Legendre & Legendre (1998) , Clarke (1993) , Clarke, Somerfield & Chapman (2006) ). Other measures, such as Bray-Curtis ( Bray & Curtis (1957) ), Kulczynski (see Faith, Minchin & Belbin (1987) ), Jaccard (see Legendre & Legendre (1998) ), binomial deviance ( Anderson & Millar (2004) ), Hellinger ( Rao (1995) ) or a modified Gower measure ( Anderson, Ellingsen & McArdle (2006) ), may be more appropriate, in different situations, for analysing community data. Although analyses on the basis of some measures (such as chi-squared or Hellinger) can be achieved by applying Euclidean distances to data that have been transformed in a particular way (see Legendre & Gallagher (2001) for details), we wish to retain the full flexibility of methods like ANOSIM or MDS and allow the analysis to be based on any reasonable distance or dissimilarity measure of choice. Unfortunately, for many resemblance measures favoured by ecologists, we cannot easily calculate distances to centroids. The reason for this is that the sample centroid (the point lying in the centre of the data cloud) in the space defined by these measures will not be the same as the vector of arithmetic averages for each variable. It is actually quite unclear just how one could go about calculating centroids for the majority of these non-Euclidean measures on the basis of the raw data alone. Thus, for these situations, we shall instead rely on what is known as Huygens’ theorem7 (Fig. 1.3). This theorem states (for Euclidean space) that the sum of squared distances from individual points to their group centroid is equal to the sum of squared inter-point distances, divided by the number of points in the group (e.g., Legendre & Anderson (1999) , see Appendix B therein). Fig. 1.3. Schematic diagram of Huygens’ theorem. This theorem appears to be very simple, and it is. In fact, it is perfectly safe to calculate these two quantities using Euclidean distances and you are encouraged to try a simple example (say in two dimensions) to prove it to yourself. Importantly, what it means for us is that we can use the inter-point distances (or dissimilarities) alone in order to calculate the sums of squares, without ever having to calculate the position of the centroids at all! Therefore, for the cases where we cannot calculate centroids directly (i.e., for most dissimilarity measures), we can use the inter-point values instead to do the partitioning. The only properties required of our dissimilarity measure in order to use this approach are that it be non-negative, symmetric (the distance or dissimilarity from point A to point B must be the same as from point B to point A), and that all self dissimilarities (from point A to itself) be zero. Virtually all measures that are worth using possess these three properties8. 6 See Warton & Hudson (2004) for an alternative point of view and counter-argument to this assumption. 7 Christian Huygens (1629-1695) was a Dutch mathematician, astronomer and physicist, famous for, among other things, patenting the first pendulum clock. 8 The mathematical and geometric properties of dissimilarity measures are reviewed by Gower & Legendre (1986) . 1.4 Sums of squares from a distance matrix We can now consider the structure of a distance/dissimilarity matrix and how sums of squares for a one-way multivariate ANOVA partitioning would be calculated (Fig. 1.4). Fig. 1.4. Calculation of sums of squares directly from a distance/dissimilarity matrix. If we let $d _ {ij}$ be the dissimilarity (or distance) between sample i and sample j, then the total sum of squares is the sum of the inter-point dissimilarities among all samples, divided by N: $$ SS _ T = \frac{1}{N} \sum _ {i=1} ^ {N-1} \sum _ {j=i+1} ^ {N} d _ {ij} ^2 \tag{1.1} $$ and the residual (within-group) sum of squares (assuming, for now, a balanced design having an equal sample size of n per group) is: $$ SS _ {Res} = \frac{1}{n} \sum _ {i=1} ^ {N-1} \sum _ {j=i+1} ^ {N} d _ {ij} ^2 \omega _ {ij} \tag{1.2} $$ where $\omega _ {ij}$ takes the value of 1 if samples i and j are in the same group, otherwise it takes the value of zero. This amounts to adding up the squares of all the dissimilarities between samples that occur within the same group. These quantities are shown schematically in Fig. 1.4. The among-group sum of squares can also be calculated directly or, more simply, as the difference: $SS _A = SS _T - SS _ {Res}$. Partitioning of distance matrices having Euclidean metric properties according to ANOVA experimental designs has been discussed previously by Edgington (1995) , Pillar & Orloci (1996) and Excoffier, Smouse & Quattro (1992) , whereas Gower & Krzanowski (1999) extended this idea to semi-metric dissimilarities having only the properties of symmetry (i.e., $d _{ij} = d _{ji}$) and that $d _ {ij} \ge 0$ and $d _ {ii} = 0$ for all samples. A related point is that PERMANOVA (or any of the other methods that partition variability in the routines offered by the PERMANOVA+ add-on) does not suffer from the problems recently identified by Legendre, Borcard & Peres-Neto (2005) for the Mantel and partial Mantel test ( Mantel (1967) , Smouse, Long & Sokal (1986) )9. Confusion might arise because Legendre, Borcard & Peres-Neto (2005) used the words “variation partitioning on distance matrices” to describe the general Mantel approach. PERMANOVA, however, is not a Mantel test. The essential issue regarding the partial Mantel approach (known for some time to be potentially problematic, see Dutilleul, Stockwell, Frigon et al. (2000) , Legendre (2000) , Raufaste & Rousset (2001) and Rousset (2002) ) is that it works on “unwound” distance matrices, where inter-point distance values are treated independently as a single vector. In contrast, the methods in the PERMANOVA+ add-on all work on the distance/dissimilarity matrix directly but, importantly, they retain its inherent structure; the values in the matrix are not unwound nor are they treated or modelled as independent of one another. The PERMANOVA+ methods are therefore directly akin to the so-called “canonical partitioning” methods referred to by Legendre, Borcard & Peres-Neto (2005) , and are correct for partitioning and analysing the actual variability inherent in multivariate data clouds. 9 The original simple Mantel test ( Mantel (1967) ) to relate two distance matrices, from which several of the PRIMER routines (RELATE, BIOENV, BEST, BVSTEP and 2-stage MDS) all drew some inspiration, is valid and does have utility in appropriate applications, as pointed out by Legendre, Borcard & Peres-Neto (2005) . 1.5 The pseudo-F statistic Once the partitioning has been done we are ready to calculate a test statistic associated with the general multivariate null hypothesis of no differences among the groups. For this, following R. A. Fisher’s lead, a pseudo-F ratio is defined as: $$ F = \frac{ SS _ A / \left( a - 1 \right)} { SS _ {Res} / \left( N - a \right) } \tag{1.3} $$ where (a – 1) are the degrees of freedom associated with the factor and (N – a) are the residual degrees of freedom. It is clear here that, as the pseudo-F statistic in (1.3) gets larger, the likelihood of the null hypothesis being true diminishes. Interestingly, if there is only one variable in the analysis and one has chosen to use Euclidean distance, then the resulting PERMANOVA F ratio is exactly the same as the original F statistic in traditional ANOVA10 ( Fisher (1924) ). In general, however, the PERMANOVA F ratio should be thought of as a “pseudo” F statistic, because it does not have a known distribution under a true null hypothesis. There is only one situation for which this distribution is known and corresponds to Fisher’s traditional F distribution, namely: (i) if the analysis is being done on a single response variable and (ii) the distance measure used was Euclidean distance and (iii) the single response variable is normally distributed. In all other cases (multiple variables, non-normal variables and/or non-Euclidean dissimilarities), all bets are off! Therefore, in general, we cannot rely on traditional tables of the F distribution to obtain a P-value for a given multivariate data set. Some other test statistics based on resemblance measures (and using randomization or permutation methods to obtain P-values, see the next section) have been suggested for analysing one-way ANOVA designs (e.g., such as the average between-group similarity divided by the average within-group similarity as outlined by Good (1982) and Smith, Pontasch & Cairns (1990) , see also all of the good ideas in the book by Mielke & Berry (2001) and references therein). Unlike pseudo-F, however, these can be limited in that they may not necessarily yield straightforward extensions to multi-way designs. 10 In fact, a nice way to familiarise oneself with the routine is to do a traditional univariate ANOVA using some other package and compare this with the outcome from the analysis of that same variable based on Euclidean distances using PERMANOVA. 1.6 Test by permutation An appropriate distribution for the pseudo-F statistic under a true null hypothesis is obtained by using a permutation (or randomization) procedure (e.g., Edgington (1995) , Manly (2006) ). The idea of a permutation test is this: if there is no effect of the factor, then it is equally likely that any of the individual factor labels could have been associated with any one of the samples. The null hypothesis suggests that we could have obtained the samples in any order by reference to the groups, if groups do not differ. So, another possible value of the test statistic under a true null hypothesis can be obtained by randomly shuffling the group labels onto different sample units. The random shuffling of labels is repeated a large number of times, and each time, a new value of pseudo-F, which we will call pseudo-$F ^ \pi$, is calculated (Fig. 1.5). Note that the samples themselves have not changed their positions in the multivariate space at all as a consequence of this procedure, it is only that the particular labels associated with each sample have been randomly re-allocated to them (Fig. 1.5). If the null hypothesis were true, then the pseudo-F statistic actually obtained with the real ordering of the data relative to the groups will be similar to the values obtained under permutation. If, however, there is a group effect, then the value of pseudo-F obtained with the real ordering will appear large relative to the distribution of values obtained under permutation. Fig. 1.5. After calculating an observed value of the test statistic, a value that might have been obtained under a true null hypothesis is calculated by permuting labels. The frequency distribution of the values of pseudo-$F ^ \pi$ is discrete: that is, the number of possible ways that the data could have been re-ordered is finite. The probability associated with the test statistic under a true null hypothesis is calculated as the proportion of the pseudo-$F ^ \pi$ values that are greater than or equal to the value of pseudo-F observed for the real data. Hence, $$ P = \frac{ \left( \text{No. of } F ^ \pi \ge F \right) +1}{ \left( \text{Total no. of } F ^ \pi \right) +1} \tag{1.4} $$ In this calculation, we include the observed value as a member of the distribution, appearing simply as “+1” in both the numerator and denominator of (1.4). This is because one of the possible random orderings of the data is the ordering we actually got and its inclusion makes the P-value err slightly on the conservative side ( Hope (1968) ) as is desirable. For multivariate data, the samples (either as whole rows or whole columns) of the data matrix (raw data) are simply permuted randomly among the groups. Note that permutation of the raw data for multivariate analysis does not mean that values in the data matrix are shuffled just anywhere. A whole sample (e.g., an entire row or column) is permuted as a unit; the exchangeable units are the labels associated with the sample vectors of the data matrix (e.g., Anderson (2001b) )11. For the one-way test, enumeration of all possible permutations (re-orderings of the samples) gives a P-value that yields an exact test of the null hypothesis. An exact test is one where the probability of rejecting a true null hypothesis is exactly equal to the a priori chosen significance level ($\alpha$). Thus, if $\alpha$ were chosen to be 0.05, then the chance of a type I error for an exact test is indeed 5%. In practice, the possible number of permutations is very large in most cases. A practical strategy, therefore, is to perform the test using a large random subset of values of pseudo-$F ^ \pi$ , drawn randomly, independently and with equal probability from the distribution of pseudo-$F ^ \pi$ for all possible permutations. PERMANOVA does not systematically do all permutations, but rather draws a random subset of them, with the number to be done being chosen by the user. Such a test is still exact ( Dwass (1957) ). However, separate runs will therefore result in slightly different P-values for a given test, but these differences will be very small for large numbers of permutations (e.g., in the 3rd decimal place for 9999 permutations) and should not affect interpretation. Once again, as a rather nice intuitive bonus, PERMANOVA done on one response variable alone and using Euclidean distance yields Fisher’s traditional univariate F statistic. So, PERMANOVA can also be used to do univariate ANOVA but where P-values are obtained by permutation (e.g., Anderson & Millar (2004) ), thus avoiding the assumption of normality. Note also that if the univariate data do happen to conform to the traditional assumptions (normality, etc.), then the permutation P-value converges on the traditional normal-theory P-value in any event. 11 Equivalently, permutations can be achieved by the simultaneous re-ordering of the rows and columns of the resemblance matrix (e.g., see Fig. 2 in Anderson (2001b) ). 1.7 Assumptions Recall that for traditional one-way ANOVA, the assumptions are that the errors are independent, that they are normally distributed with a mean of zero and a common variance, and that they are added to the treatment effects. In the case of a one-way analysis, the PERMANOVA test using permutations assumes only that the samples are exchangeable under a true null hypothesis12. The assumption of exchangeability is tantamount to assuming that the multivariate observations (samples) are independent and identically distributed (i.i.d.) under a true null hypothesis. Thus, although there are no explicit assumptions regarding the distributions of the original variables (they are certainly not assumed to be normally distributed), independence and homogeneity of dispersions (in the space of the resemblance measure) are directly implied by the permutation procedure. Clearly, if samples have very different dispersions in different groups, then they are not really exchangeable. Also, if the samples are unequally correlated with one another (e.g., temporally or spatially), then randomly shuffling them will destroy this kind of inherent structure. In contrast, it is not expected that the individual variables which have been measured on the same samples (in the multivariate case) are independent of one another and this is not assumed. When permutations are done, the values for different variables within a sample are kept together as a unit, so whatever correlation structure there might be among the variables is not altered under permutation. Fig. 1.6. PERMANOVA will be sensitive to differences in dispersions (left) but not differences in correlation structure (right) among groups. PERMANOVA, like ANOSIM ( Clarke (1993) ), will be sensitive to differences in dispersion among groups (Fig. 1.6). Indeed, the construction of pseudo-F ratios in PERMANOVA uses pooled estimates of within-group variability, so homogeneity of multivariate dispersions is also implicit in the partitioning. A separate test for homogeneity of dispersions, using the PERMDISP routine (see chapter 2) can be done prior to performing PERMANOVA (or, indeed, to investigate the null hypothesis of homogeneity in its own right). However, we consider that a non-significant result from PERMDISP is not strictly necessary to achieve prior to using PERMANOVA. It is likely that PERMDISP will detect differences in dispersion that, in many cases, are not substantial enough to “de-rail” (i.e. to inflate the error rates of) the PERMANOVA test13. This is analogous to the situation in univariate analysis; traditional ANOVA is quite robust to many forms of heterogeneity, especially with large sample sizes ( Box (1953) ). Instead, we can consider the homogeneity of dispersions to be included as part of the general null hypothesis of “no differences” among groups being tested by PERMANOVA (even though the focus of the PERMANOVA test is to detect location effects). If significant heterogeneity were detected by PERMDISP and differences among groups were also detected using PERMANOVA, then the latter could have been caused by differences in location, differences in dispersion, or some combination of the two. Thus, performing a test using PERMDISP, as well as examining the average within and between-group dissimilarities and the position of samples from different groups in unconstrained ordination plots (MDS or PCO), will help to uncover the nature of any differences among groups detected by PERMANOVA. Unlike many of the traditional MANOVA test statistics (e.g., Mardia, Kent & Bibby (1979) , Seber (1984) ), PERMANOVA will not be sensitive, however, to differences in correlation structure among groups (Fig. 1.6). See Krzanowski (1993) for a permutation test designed to compare correlation structures among variables across different groups. Fig. 1.7. Producing an MDS plot of the Ekofisk macrofauna data. 12 Exchangeability of multivariate observations (samples) is assured if we have done a random allocation of sample units to groups or treatments a priori ( Fisher (1935) ). For observational studies, where we cannot do this (i.e., the groups already occur in nature and we draw a random sample from them), we must assume exchangeability under a true null hypothesis ( Kempthorne (1966) ). 13 Certainly a worthy topic for future study is to discover the conditions under which PERMANOVA will show inflated rates of either Type I or Type II error in the face of heterogeneity in the distributions of multivariate samples among groups. 1.8 One-way example (Ekofisk oil-field macrofauna) Our first real example comes from a study by Gray, Clarke, Warwick et al. (1990) , who studied changes in community structure of soft-sediment benthic macrofauna in relation to oil-drilling activity at the Ekofisk oil platform in the North Sea. These data consist of p = 174 species sampled by bottom grabs at each of N = 39 stations. The stations were placed roughly along five transects radiating out from the centre of the oil platform. Stations have been grouped with labels according to their distance from the oil platform as A (> 3.5 km), B (1 km – 3.5 km), C (250 m – 1 km) and D (< 250 m). The data are located in the file ekma.pri in the ‘Ekofisk’ folder of the ‘Examples v6’ directory. Open the file by selecting File>Open from within PRIMER and using the browser. Next, select Edit>Factors to view the factor labels associated with each sample. See chapter 2 of Clarke & Gorley (2006) for detailed information concerning creating, importing and editing factors to identify sample groups within PRIMER. Of interest here is to test the null hypothesis of no differences among the communities inhabiting the benthic habitats in these four different groups. First, we may wish to visualise the relationships among the samples in terms of a relevant resemblance measure, using ordination. A well-known robust procedure for doing this is non-metric multidimensional scaling14 (MDS, Shepard (1962) , Kruskal (1964) , Kruskal & Wish (1978) , Minchin (1987) ). Produce a resemblance matrix among the samples by selecting Analyse > Pre-treatment > Transform (overall) > Transformation: fourth-root, followed by Analyse > Resemblance > (Analyse between$\bullet$Samples) & (Measure$\bullet$Bray-Curtis similarity). Next, produce an MDS plot by selecting Analyse > MDS and click ‘OK’ with all of the default options. Once the 2-dimensional graph is in view, show the samples according to their factor labels by selecting Graph > Data labels & symbols > (Labels > $\checkmark$Plot > $\checkmark$By factor Dist) and by removing the $\checkmark$ from the (Symbols > Plot) box (Fig. 1.7). The resulting ordination plot suggests that the communities are fairly distinct in the different distance groups, and also that they tend to occur along a gradient, from those lying closest to the platform (group D) to those lying furthest away (group A). 14 Here and throughout, we shall use the acronym MDS to denote non-metric (as opposed to metric) multi-dimensional scaling. More details about the method of MDS and its implementation in PRIMER can be found in chapter 5 of Clarke & Warwick (2001) and chapter 7 of Clarke & Gorley (2006) . 1.9 Creating a design file We shall formally test the hypothesis of no differences in community structure among the four groups (where, in this case, “differences in community structure” is defined by the Bray-Curtis measure on fourth-root transformed data) using PERMANOVA. For any analysis using PERMANOVA, there are necessarily two steps involved. First, one must create a design file, which provides all of the necessary information for the partitioning to be done according to the correct factors and experimental design. The design file is unique to the PERMANOVA+ add-on package, with its own special icon:. It has a direct link to the factor information associated with the samples in the resemblance matrix. (The factor information associated with the raw data matrix is inherited by any resemblance matrix produced from the data.) Once an appropriate design file has been created, the second step is to run the PERMANOVA analysis itself on the resemblance matrix, providing the name of the design file. Although the above general description of the method has used a distance or dissimilarity matrix at its base, PERMANOVA will also perform a correct analysis when given a resemblance matrix of similarities. For example, in the present case, the Bray-Curtis dissimilarity is defined simply as 100 minus the Bray-Curtis similarity. This is done automatically and the user need not perform any extra steps. To create a design file for the Ekofisk macrofauna example (Fig. 1.8), use the explorer tree in the left-hand panel and click on the resemblance matrix, then select PERMANOVA+ > Create PERMANOVA design. In the initial dialog box, entitled ‘PERMANOVA design properties’, specify (Title: Ekofisk oilfield macrofauna) & (Number of factors: 1). This will bring up a new design worksheet file with one row for each factor (so only one row in the present example) and four columns. The first column is used to indicate the name of each factor, the second column is for specifying any factors within which it may be nested, the third column is for specifying whether the factor is fixed or random and the fourth column is for specifying specific contrasts among levels of the factor. You can either begin by typing the name of the factor inside the cell of the first column, or by double clicking inside this cell and a list of all the factors associated with the resemblance matrix will be brought up automatically for you to choose from. For this data set, select the factor Dist, as shown (Fig. 1.8). This factor is not nested within any other factor (leave the cell in the second column blank), it is fixed (select ‘Fixed’ in the third cell, which is the default), and we do not wish at present to specify any specific contrasts (leave the cell in the fourth column blank). Once the design file is created, you may change its name, edit its properties and so on, as for any other PRIMER file within the workspace. For example, select Edit from the main menu to see how rows (i.e., factors) may be inserted, moved or deleted. Rows of the design file can also be edited by selecting them directly, or by right-clicking on the design file itself, as you would for editing a worksheet (e.g., see chapter 1 of Clarke & Gorley (2006) ). Like other PRIMER files, the design file can also be saved separately in its own format (*.ppd). For the present example, change the name of the design file for the Ekofisk data from the default (Design1) to One-way by right-clicking on its name in the explorer tree and selecting Rename item. Save the entire workspace with the name ekofisk.pwk for continued analysis. Fig. 1.8. Creating a design file for the Ekofisk one-way design. 1.10 Running PERMANOVA To run PERMANOVA on the Ekofisk data, click on the resemblance matrix and select PERMANOVA+ > PERMANOVA. In the PERMANOVA dialog box (Fig. 1.9), leave the defaults for all options, except (Num. permutations: 9999) & (Permutation method: $\bullet$Unrestricted permutation of raw data). Although any of the available methods of permutation offered by PERMANOVA are sound (see the section Methods of permutation), permutation of raw data will provide an exact test for a simple one-way design. In addition, although the default for the number of permutations is 999, it is clearly desirable to perform as many permutations as reasonable time will allow. Power and precision increase with increases in the number of permutations ( Hope (1968) ). Manly (2006) , (pp. 94-98), suggested that, to draw inferences at a significance level of 0.05, P-values should be calculated using at least 999 permutations, whereas 4999 permutations should be done to draw inferences at a level of 0.01. The ever-increasing speed of personal computers generally allows a large number to be chosen here for most designs with moderate sample sizes. Recall also that PERMANOVA obtains a random subset of all possible permutations, so will not necessarily reproduce exactly the same P-value for a given test if the analysis is run again. However, any such difference from one run to the next will be small (in the 3rd decimal place for 9999 permutations). Fig. 1.9. The dialog box for running PERMANOVA. The results of PERMANOVA are shown in a new separate window with text-form information for this analysis (Fig. 1.10). The first part of the file provides information regarding the choices made, such as transformations, the resemblance measure and the method and number of permutations. There is also information regarding the experimental design and whether any terms were excluded (see the section Pooling or excluding terms). The default ‘Type’ of sums of squares is Type III, which will be discussed in detail in the section Unbalanced designs. These three types of sums of squares are equivalent for one-way models or any balanced (equal replication) ANOVA designs, so need not concern us here. Fig. 1.10. Window of results from PERMANOVA on the Ekofisk macrofauna. The essential information of interest is provided in the ‘PERMANOVA table of results’, which contains the sources of variation in the model (‘Source’), the degrees of freedom (‘df’), sums of squares (‘SS’), mean squares (‘MS’), pseudo-F ratio (‘Pseudo-F’) and permutation P-value (‘P(perm)’). This table can be read and interpreted in a way that is directly analogous to a traditional ANOVA, but bear in mind that what is being partitioned here is multivariate variability based on the chosen resemblance measure, with P-values obtained using permutations15. The last column of the PERMANOVA table (‘Unique perms’) indicates how many unique values of the test statistic were obtained under permutation. Recall that PERMANOVA does not systematically do all permutations, but rather draws a random subset of them. In the example (Fig. 1.10), this value is very large (9862) and close to the number of random permutations that were chosen to be done by the user (9999). This means that only a few repeated values of pseudo-F$\pi$ were encountered under permutation, and the number of unique values is plenty enough to make reasonable inferences using the resulting permutation P-value, as shown. This information is important because, in some cases, the number of possible permutations is not large, and very few unique values of the test statistic are obtained. In such cases, a more meaningful (but approximate) P-value can be obtained by random sampling from the asymptotic permutation distribution instead (see the section Monte Carlo P-values). Underneath the PERMANOVA table of results are given further details regarding the analysis, including the expected mean squares (EMS) for each term in the model, the construction of the F ratio and the estimated sizes of components of variation (see the sections Components of variation, Expected mean squares, Constructing F from EMS and Estimating components of variation). For the Ekofisk example, there is indeed a significant effect of the distance groupings (Fig. 1.10, pseudo-F = 4.93, P = 0.0001). A natural next question to ask is: wherein do these significant differences lie? That is, which groups differ significantly from which other groups? 15 For the one-way case, such as this, the PERMANOVA permutation P-values are exact; however, for multi-factor models, especially mixed models, analyses with covariates or unbalanced designs, the P-values provided are, necessarily, approximate, and rely on the exchangeability of residuals from the linear ANOVA model being fitted. 1.11 Pair-wise comparisons Pair-wise comparisons among all pairs of levels of a given factor of interest are obtained by doing an additional separate run of the PERMANOVA routine. This is appropriate because which particular comparisons should be done, in most cases, is not known a priori, but instead will follow logically from the specific terms in the model found to be statistically significant in the main PERMANOVA analysis. We could use the F ratio for the pair-wise tests as well, but a more natural statistic to use here is a direct multivariate analogue to the univariate t statistic. In traditional univariate analysis, an F test to examine the effects of a factor having only two groups is equivalent to performing a two-sample (two-tailed) t test. In fact, in this case the t statistic for comparing two groups is simply equal to the square root of the F ratio. Similarly, when PERMANOVA performs pair-wise tests, it takes each pair of levels to be compared, in turn and, treating it like a specific contrast, calculates pseudo-t as the square root of pseudo-F. If a single response variable is being analysed and the resemblance measure chosen is Euclidean distance, then the t statistics calculated by PERMANOVA for the pair-wise tests correspond exactly to Gosset’s original t statistic ( Student (1908) ). However, unlike traditional statistics packages, P-values for all pair-wise tests in PERMANOVA are obtained using permutations, not tables. This is necessary because, if the data are not normal, if more than one variable is being analysed and/or if the distance measure used is not Euclidean distance, then the distribution of this pseudo-t statistic under a true null hypothesis is unknown, so (as with pseudo-F), we generate it using permutations. To run pair-wise comparisons for the Ekofisk example, click on the window containing the resemblance matrix once again and select PERMANOVA+ > PERMANOVA. All of the choices in the dialog box that were made for the main test should still be there. (PRIMER v6 has a good memory!) Most of these will remain as before. The same design file is needed (One-way), but this time choose (Test > $\bullet$Pair-wise test > For term: Dist > For pairs of levels of factor: Dist). The file of results will contain the same preliminary information regarding the design and other choices as was seen for the main PERMANOVA test, but then the individual comparisons, using pseudo-t and with P-values by permutation, are given for each pair of levels of the factor. Generally speaking, as with the F ratio, the larger the t statistic, the greater the evidence against the null hypothesis of no difference in community structure between the two groups. Fig. 1.11. Pair-wise tests among groups for the Ekofisk macrofauna. For the Ekofisk example, there is fairly strong evidence to suggest that all of the groups differ from one another (P < 0.001 for most comparisons, Fig. 1.11). However, the evidence against the null hypothesis for the comparison of group B with group C is weaker than the rest (Fig. 1.11, t = 1.24, P = 0.021). A relevant point to note here is that no corrections for multiple comparisons have been made to any of these tests. It is a well-known phenomenon in statistics that the more tests you do, the greater your chance of rejecting one or more null hypotheses simply by chance. For example, if we were to perform 20 tests using an a priori significance level of 0.05, we would expect to reject a true null hypothesis (i.e., to get a P-value smaller than 0.05) in one of those 20 tests by chance alone. However, the permutation P-values do provide an exact test of each individual null hypothesis of interest. In contrast, most ad hoc experiment-wise corrections that could be used here (such as Bonferroni) are inexact and known to be overly conservative (e.g., Day & Quinn (1989) ). Thus, our philosophy (here and elsewhere) is to report the exact permutation P-values directly and to let the user decide whether or not to apply any additional corrections. We recommend that one should consider the set of tests as a whole, in probabilistic terms. If there are a groups, then there will be a(a – 1)/2 tests. For the Ekofisk example, this is 4(4 – 1)/2 = 6 tests. Clearly, in the present case, if all null hypotheses were true, it would be highly unlikely to get six out of six P-values less than 0.05 (as we have done) simply by chance! If, however, the number of tests were very large, with few small P-values encountered (e.g., if one obtained only one or two P-values less than 0.05 out of 20 or so tests), then one might choose to apply a formal correction or, at the very least, exercise caution in interpreting the results16. The routine for pair-wise tests also provides a triangular matrix containing the average resemblances between samples that are either in the same group (along the diagonal) or in different groups (sub-diagonal elements, Fig. 1.11). This identifies the relative sizes of average similarities (or dissimilarities) between each pair of groups in the multivariate space, and also helps (along with the formal assessment provided in the PERMDISP routine, see chapter 2) to identify potential differences among the groups in terms of their within-group variability. 16 See Wheldon, Anderson & Johnson (2007) , in which an exact correction for family-wise error across a large set of non-independent simultaneous multivariate permutation tests was achieved using the permutation distributions. 1.12 Monte Carlo P-values (Victorian avifauna) In some situations, there are not enough possible permutations to get a reasonable test. Consider the case of two groups, with two observations per group. There are a total of 4 observations, so the total number of possible re-orderings (permutations) of the 4 samples is 4! = 24. However, with a groups and n replicates per group, the number of distinct possible outcomes for the F statistic in the one-way test is $(an)! / [ a! ( n ! ) ^ a ]$ (e.g., Clarke (1993) ), which in this case is: $(4)! / [2! (2!)^2] = 3 $ unique outcomes. This means that even if the observed value of pseudo-F is quite large, the smallest possible P-value that can be obtained is P = 0.333. This is clearly insufficient to make statistical inferences at a significance level of 0.05. An alternative is to use the result given in Anderson & Robinson (2003) regarding the asymptotic permutation of the numerator (or denominator) of the test statistic under permutation (see equations (1) and (4) on p. 305 of Anderson & Robinson (2003) ). It is demonstrated (under certain mild assumptions17 that each of the sums of squares has, under permutation, an asymptotic distribution that is a linear form in chi-square variables, where the coefficients are actually the eigenvalues from a PCO of the resemblance matrix (see chapter 3). Thus, chi-square variables can be drawn randomly and independently, using Monte Carlo sampling, and these can be combined with the eigenvalues to construct the asymptotic permutation distribution for each of the numerator and denominator and, thus, for the entire pseudo-F statistic, in the event that too few actual unique permutations exist. A case in point where such an issue arises is given by a study of Victorian avifauna by Mac Nally & Timewell (2005) . The data consist of counts of p = 27 nectarivorous bird species at each of eight sites of contrasting flowering intensity within the Rushworth State Forest in Victoria, Australia. One pair of sites had heavy flowering (‘good’ sites), another pair had intermediate flowering (‘medium’ sites), and a third pair had relatively little flowering (‘poor’ sites). Two sites near the good sites (‘adjacent’ sites), were also selected to explore possible “spill-over” effects. Sites were sampled using a strip transect method and data from four surveys were summed for each site ( Mac Nally & Timewell (2005) ). The data are located in the file vic.pri in the ‘VictAvi’ folder of the ‘Examples add-on’ directory. An MDS ordination of the 8 sites on the basis of the binomial deviance dissimilarity measure ( Anderson & Millar (2004) ) indicates potential differences among the four groups of sites in terms of their avifaunal community structure and an overall gradient from good to poor sites (Fig. 1.12). Fig. 1.12. MDS plot of Victorian avifauna at sites with different flowering intensities. We have a = 4 and n = 2, so there will be a total of $(8)!/[4!(2!)^4] = 105$ unique values of the pseudo-F ratio under permutation for the main PERMANOVA test. This will provide some basis for making inferences from the P-value. For the pair-wise tests, however, there will be only 3 unique values under permutation, so obtaining Monte Carlo P-values would clearly be useful for these. To obtain Monte Carlo results, in addition to the permutation P-values, simply choose ($\checkmark$Do Monte Carlo tests) in the PERMANOVA dialog. For the Victorian avifauna, the PERMANOVA main test suggests that there are significant differences in bird communities among the four types of sites (Fig. 1.13). The permutation and Monte Carlo P-values are quite similar in value (‘P(perm)’ = 0.02 and ‘P(MC)’ = 0.03) and note how the routine also faithfully shows that 105 unique values of pseudo-F were obtained under permutation (‘Unique perms’). Fig. 1.13. PERMANOVA results for the main test and pair-wise tests for the Victorian avifauna, including Monte Carlo P-values. The P-values for the pair-wise tests are virtually meaningless in this example, as there are only 3 unique possible values under permutation in each case. However, the Monte Carlo P-values are clearly much more valuable here and can be interpreted. Sites with medium flowering intensity appear to differ significantly from those with good flowering intensity (‘P(MC)’ = 0.02), which is also reflected in the relative size of the pseudo-t statistic (= 4.16). None of the other P-values were less than 0.05, although the comparison of poor with good sites came close (‘P(MC)’ = 0.06). One possible reason this latter comparison did not achieve a higher t statistic is because of the large within-group variation of the poor sites (average within-group distance = 66.0, Fig. 1.13). If the number of unique permutations is large (say, 100 or more), then the permutation P-value should be preferred over the Monte Carlo P-value, as a general rule, because it will provide a more accurate test. The Monte Carlo P-values are based on asymptotic theory (albeit a very robust theory that generates distributions for test statistics which are specific to each dataset), so they rely on large-sample approximations. When there are a large number of possible permutations, then the Monte Carlo and permutation P-values should be very similar, essentially converging on the same answer. When, on the other hand, there are very few possible permutations, then the two P-values may be quite different (as seen in the Victorian avifauna example), in which case the Monte Carlo P-value should probably be used in preference. Although the Monte Carlo P-value is, at least, interpretable in cases such as these, it is only an approximation which relies on the central limit theorem, so clearly we shouldn’t get too carried away with our statistical inferences and interpretation of P-values from studies having such small sample sizes18. In multi-factor designs (discussed below), it is not always easy to calculate how many unique permutations there are for a given term in the model. It depends not only on the permutation method chosen, but also potentially depends on other factors in the model, how many levels they have and whether they may, under permutation, coincide with the term being tested in serendipitous ways. Therefore, for each test, PERMANOVA simply keeps track of how many unique values of the test statistic it encounters out of the total number of random permutations done. Armed with this information (‘Unique perms’), it is then possible to judge how useful the permutation P-value given actually is, under the circumstances. 17 The essential assumptions here are (a) that the observations of multivariate observations under permutation are independent and identically distributed, (b) that the distances are not just governed by a couple of very large ones (so that the central limit theorem can apply) and (c) that the distances are not too discontinuous – so a small change in the data would not produce a large change in the distance (sensible distance functions satisfy this). See Anderson & Robinson (2003) for details. 18 Unfortunately, it is precisely in these situations where sample sizes are small and we would like to use the Monte Carlo P-value that the asymptotic approximations assumed by this approach are actually the most precarious! A clear topic for further study is to discover under what conditions, more specifically, the Monte Carlo P-values may become unreliable. 1.13 PERMANOVA versus ANOSIM The analysis of similarities (ANOSIM), described by Clarke (1993) is also available within PRIMER and can be used to analyse multivariate resemblances according to one-way and some limited two-way experimental designs19. Not surprisingly, ANOSIM and PERMANOVA will tend to give very similar results for the one-way design on a given resemblance matrix. There are two essential differences, however, between ANOSIM and PERMANOVA. First, ANOSIM ranks the values in the resemblance matrix before proceeding with the analysis. The rationale behind the ranking procedure in ANOSIM is that the information of interest is the relationships among the dissimilarities (i.e., whether a given dissimilarity is larger or smaller than another) and not the values of the dissimilarities themselves. This is consistent with the philosophy of non-metric MDS ordination, which seeks to preserve only the rank order of the dissimilarities among samples. In contrast, PERMANOVA takes the point of view that the information of interest is in the dissimilarity values themselves, which describe a cloud of samples in multivariate space. This means that for PERMANOVA one must take special care to choose a measure of resemblance that is meaningful for the data and the goals of the analysis. For example, squared Euclidean distance may give different PERMANOVA results than Euclidean distance itself, whilst such a monotonic transform of the resemblances does not change the ranks and therefore cannot change ANOSIM. The second essential difference is in the construction of the test statistic. The ANOSIM R statistic ( Clarke (1993) ) is scaled to take a value between -1 and +1. This is a very useful feature, as it makes it possible to interpret the R statistic directly as an absolute measure of the strength of the difference between groups. R values are also directly comparable among different studies. In contrast, the value of pseudo-F (or pseudo-t) is, first of all, necessarily reliant on the degrees of freedom of the analysis, so cannot necessarily be compared in value across studies. A value of pseudo-F = 2.0 (like its univariate analogue) will generally provide much stronger evidence against the null hypothesis if the residual degrees of freedom are 98 than if they are 5. Although values of pseudo-F may be comparable across different tests where the degrees of freedom are equal (for a given dissimilarity measure and original number of variables, that is), it is also worth bearing in mind that the variability among groups (as measured by the numerator of the statistic) is always scaled against the variability within groups (as measured by the denominator). Thus, the within-group variability has an important role to play in the value of pseudo-F (or pseudo-t). An example is provided by the Victorian avifauna comparisons (Fig. 1.13), where, despite the pattern shown on the MDS plot (that samples from poor sites are farther away from the good sites than are the medium sites), pseudo-t is actually larger for the difference between good and medium sites than it is between good and poor sites, simply because the within-group variability between the poor sites is so high. ANOSIM, in contrast, yields an R statistic value of 1.0 (its maximum possible value) in both cases. In summary, while ANOSIM’s R can be interpreted directly as a measure of the size of the between-group differences, PERMANOVA’s pseudo-F (or pseudo-t) cannot necessarily be interpreted in this way. The sizes of effects in PERMANOVA are measured and compared in other ways: either by the average similarities (or dissimilarities) among pairs of groups (provided by the pair-wise routine) or by examining the estimated sizes of components of variation (see the section Estimating components of variation). In addition, in PERMANOVA it is the P-values (either ‘P(perm)’ or, when necessary, ‘P(MC)’) which should be used as a measure of strength of evidence with respect to any particular null hypothesis. The Monte Carlo option available here also means that the power of the test need not especially rely on the number of possible permutations, as is the case for ANOSIM. Power in PERMANOVA will rely, however, on the number of replicates (more particularly, on the denominator degrees of freedom) available for the test. Unlike ANOSIM, PERMANOVA achieves a partitioning of multivariate variability. As discussed in the introduction (page 0.3), this means PERMANOVA can be used to analyse much more complex experimental designs than ANOSIM. Although one could conceivably rank the dissimilarities before proceeding with a PERMANOVA analysis, this is not generally advisable when the goal is to achieve a partitioning of multivariate variability. The reason is that ranking dissimilarities loses information and therefore may result in less power20. Another reason is that ranking the dissimilarities will tend to make the multivariate system highly non-metric, which can result in negative sums of squares and thus negative values of pseudo-F! The concept of negative variance which arises in non-metric or semi-metric geometric systems is discussed in more detail by Legendre & Legendre (1998) and McArdle & Anderson (2001) and in chapter 3 below on principal coordinates analysis (PCO). Suffice it for now to state simply that such results are confusing and difficult to interpret. They usually result from a poor choice of resemblance measure, or from ranking resemblances unnecessarily, so should be avoided if possible. The partitioning of variability described by the resemblance matrix on the basis of most reasonable dissimilarity measures (that have not been ranked) will generally produce a result where all of the SS (and pseudo-F ratios) are positive. 19 See chapter 12 in Clarke & Gorley (2006) and chapter 6 in Clarke & Warwick (2001) . 20 This is analogous to the way that non-parametric univariate statistics are less powerful than their more traditional parametric counterparts when the assumptions of the latter are fulfilled. Interestingly enough, distance-based permutation tests (using Euclidean distance) can achieve even greater power than the traditional MANOVA test statistics in some situations, even when the assumptions of the traditional tests are true ( Smith (1998) , Mielke & Berry (2001) , see also chapter 5 on CAP). 1.14 Two-way crossed design (Subtidal epibiota) The primary advantage of PERMANOVA is its ability to analyse complex experimental designs. The partitioning inherent in the routine allows interaction terms in crossed designs to be estimated and tested explicitly. As an example, consider a manipulative experiment done to test the effects of shading and proximity to the seafloor on the development of subtidal epibiotic assemblages in Sydney Harbour, Australia ( Glasby (1999) ). The data are located in the file sub.pri in the ‘SubEpi’ folder of the ‘Examples add-on’ directory. It was observed that assemblages colonising pilings at marinas were different from those colonising nearby natural sandstone rocky reefs. Glasby (1999) proposed that this difference might be due to the fact that assemblages on pilings are far from the seafloor, or it might be caused by them being shaded by the pier structure. The experiment he designed to test these ideas (Fig. 1.14a) therefore consisted of two factors: Factor A: Position (fixed with a = 2 levels, either near (N) or far (F) from the sea floor). Factor B: Shade (fixed with b = 3 levels, shaded by an opaque Perspex roof (S), open (O) or a procedural control (C), consisting of a clear Perspex roof). The design was balanced, with n = 4 replicate sandstone plates (15 × 15 cm2) deployed in each of the a × b = 6 combinations of the two factors for a total of N = a × b × n = 24 samples in the experiment. The six combinations of the treatment levels (Fig. 1.14b). A balanced design is defined as a design that has a complete cell structure (i.e., where all combinations of the factors of interest are represented) and an equal number of replicate samples within every one of its cells. The percentage cover21 of each of 45 different taxa colonising the plates after a period of 33 weeks were measured. (There are p = 46 variables in the sub.pri worksheet; one of the variables, ‘bare’, is a measure of the unoccupied space on each plate). Fig. 1.14. Schematic diagram of a two-way crossed experimental design (a) as originally described; (b) in terms of the cell structure; and (c) as in (a), but with the order of the factors swapped. This design is crossed because, for every level of factor A, we find all levels of factor B, and vice versa. A crossed design is also identifiable as one for which we could just as easily have swapped the order of the factors, either in the description or as drawn in a schematic diagram (Fig. 1.14c). A common mistake is to assume that one factor is nested within another (see Nested designs), just because of the order in which the factors happen to be listed, or the perception that one factor physically occurs within another (e.g., treatments within sites or blocks), when, in fact, the order of the factors is readily swapped (e.g., because all treatments occur at all sites and the treatments have the same name and meaning at all sites). Perhaps the most important additional feature of a crossed design is that the two factors can interact. By an interaction, we mean that the size and/or the direction of the effects for one factor are not consistently additive, but instead they depend on which level of the other factor you happen to be talking about. An interaction between two factors contributes an additional term to the model, so an additional sum of squares appears as a ‘Source’ in the PERMANOVA table. For the present experimental design, partitioning of the total variation therefore yields the following terms: $SS _ P = \text{sum of squares due to the Position factor}$; $SS _ S = \text{sum of squares due to the Shade factor}$; $SS _ {P \times S} = \text{sum of squares due to the interaction between the two factors}$; and $SS _ {Res} = \text{residual sum of squares} $. Also, because the design is balanced, these individual terms are independent of one another and add up to the total sum of squares: $SS _ T = SS _ P + SS _ S + SS _ {P \times S} + SS _ {Res}$. For these data, a traditional MANOVA analysis is out of the question. First, the variables clearly do not fulfill the assumptions (multivariate normality, homogeneity of variance-covariance matrices among groups, etc.). Second, the number of variables (p = 46) exceeds the total number of samples (N = 24). Third, Euclidean (or Mahalanobis) distance is not an appropriate choice here, given the large number of zeros. Thus, we shall analyse the variability in this multivariate system on the basis of Bray-Curtis dissimilarities, with a fourth-root transformation of the densities in order to downweight the influence of the more abundant taxa. Open the file sub.pri and look at the factors by choosing Edit > Factors. Do an overall fourth-root transformation, then obtain a resemblance matrix based on Bray-Curtis. Fig. 1.15. PERMANOVA analysis of the two-way crossed design for subtidal epibiota. Analysis by PERMANOVA proceeds according to the two steps previously outlined for the one-way case: (i) create the design file and (ii) run the PERMANOVA routine. For the first step, choose PERMANOVA+ > Create PERMANOVA design > (Title: Subtidal epibiota) > (Number of factors: 2). The design file will begin by showing 2 blank rows. Double click in the first cell of each row and select the name of each factor in turn: Position in row 1 and Shade in row 2. Choose Fixed for both of the factors in column 3 and leave columns 2 and 4 of the design file blank (Fig. 1.15). Choose File > Rename Design and rename the design file Two-way crossed, for reference. Save the workspace as sub.pwk. Unlike the PERMANOVA routine available in DOS, it is not necessary for the order of the factors in the design file to match the order in which the factors happen to appear in the data file. PERMANOVA+ for PRIMER uses label matching to identify factors and their levels and so the order of the factors with respect to one another as they are designated within either the data worksheet or the resemblance matrix has no consequence for the analysis. However, the order of the terms in the design file does have a consequence: by default the factors are fit in the order in which they are listed in the design file. For balanced experimental designs (like this one), changing the order of terms in the model has no effect on the results anyway22. The order will matter, however, for unbalanced designs or designs with covariables if Type I sums of squares are used (see the section on Unbalanced designs). Next, with the resemblance matrix highlighted, choose PERMANOVA+ > PERMANOVA > (Design worksheet: Two-way crossed) & (Test $\bullet$Main test) & (Sums of Squares $\bullet$Type III (partial)) & (Permutation method $\bullet$Permutation of residuals under a reduced model) & (Num. permutations: 9999). This is a more complex experimental design than the one-way case (indeed, most experimental designs involve more than one factor), so here we shall use the default method of permutation of residuals under a reduced model (see Methods of permutation for more details). The results (Fig. 1.15) indicate that there is no statistically significant interaction in the effects of Position and Shade on variability in these assemblages, although the P-value is not large (P = 0.09). Each of the main factors do appear, however, to have strong effects (P < 0.001 for each test). These results are also reflected in the patterns seen in the MDS plot of the two factors, where labels are used to denote different levels of the Position factor and symbols are used to denote different levels of the Shade factor (Fig. 1.16). Fig. 1.16. MDS plot of subtidal epibiota with labels for the two-way crossed design. There is a clear separation in the MDS plot between assemblages near (N) versus those far (F) from the seafloor, suggesting a strong effect of Position. There is also a tendency for the samples in the shaded treatments to occur in the lower left of the diagram, compared to those in either the procedural control or in open treatments, which occur more in the centre or upper right. Although the size of the shade effect appeared to be a bit larger for assemblages near to the seafloor than for those far from the seafloor (Fig 1.16), the direction was very similar and the test revealed no significant interaction. 21 Taxa that were present but occupied less than 1% cover were given an arbitrary value of 0.5% 22 To see this, run the PERMANOVA analysis where the first factor listed in the design file (in the first row) is Shade and the second factor in the design file is Position. 1.15 Interpreting interactions What do we mean by an “interaction” between two factors in multivariate space? Recall that for a univariate analysis, a significant interaction means that the effects of one factor (if any) are not the same across levels of the other factor. An interaction can be caused by the size and/or the direction of the effect being different within different levels of the other factor. Similarly, in multivariate space, we can consider that a single factor has an effect on the combined set of response variables when it causes a shift in the location of the data cloud. If the magnitude and/or the direction of that shift for a given factor (e.g., factor A) is different when considered separately within each level of the other factor (say, factor B), then we may consider that those two factors interact23. Fig. 1.17. Two-dimensional examples of two-way crossed designs in which: (a) effects of A occur in a similar direction and are of a similar magnitude within each level of B, and vice versa; (b) A apparently has no effects, while B has similar effects within each level of A; (c) effects of A depend on which level of B you are in, having either no effect (within B1), a clear effect (within B2) or a clear and large effect (within B3); and (d) interaction due to not just the size, but also the direction of the effects of A within each level of B. To clarify ideas, consider a hypothetical two-way crossed experimental design with factor A having 2 levels (A1 and A2), factor B having 3 levels (B1, B2 and B3) and which can be drawn in two dimensions in Euclidean space (Fig. 1.17). We can conceive of a situation where the effect (shift in location) due to factor A (e.g., from triangles to circles) is of a similar size and in a similar direction, regardless of which level of factor B we happen to be in (Fig. 1.17a). In this case, there is no interaction. Similarly, we can conceive of a situation where there are apparently no effects of factor A at all, regardless of factor B, which would also result in there being no significant interaction (Fig. 1.17b). In contrast, we can conceive of situations where the sizes (or directions) of factor A’s effects are different within different levels of factor B (Fig. 1.17c, d). Note that, for univariate analysis, we consider changes in direction only along a single dimension for a given response variable, with an effect being either positive or negative. In contrast, there are many possible ways that changes in direction can happen in multivariate space. The logic of the interpretation of a test for interaction in PERMANOVA follows, by direct analogy, the logic employed in univariate ANOVA. In particular, it is important to examine the test of the interaction term first in order to know how to proceed. As for univariate analysis, the presence of a significant interaction generally indicates that the test(s) of main effects (i.e., the test of each factor alone, ignoring the other factor) may not be meaningful. Appropriate logic dictates that the next step after obtaining a significant interaction is to do pair-wise comparisons for the factor of interest separately within each level of the other factor (i.e., to do pair-wise comparisons among levels of factor A within each level of factor B) and vice versa (if necessary). On the other hand, if the interaction term is not statistically significant, then this indicates effects (if any) of each factor do not depend on the other factor, and one may examine the individual pseudo-F tests of each factor alone (and do subsequent pair-wise comparisons, if appropriate), ignoring the other factor. For the subtidal epibiotic assemblages, the lack of a significant interaction meant that we could logically proceed to investigate the main effects, which were each highly statistically significant. The next logical step is to consider pair-wise comparisons separately for each factor. For the factor of Position, there are only 2 levels, so the pseudo-F ratio demonstrates already that there is a significant difference between these two groups (near and far). (We might decide nevertheless to do pair-wise comparisons for this term, if we wish to know something more about the sizes of the average within-group or between-group dissimilarities, for example.) On the other hand, the factor of Shade has three levels and pair-wise comparisons are desirable to discern wherein the significant differences may lie (i.e., between which pairs of groups). Run the PERMANOVA routine again with all of the same choices in the dialog as for the main test, except this time choose (Test $\bullet$Pair-wise test > For term: Shade > For pairs of levels of factor: Shade) (Fig. 1.18). The results show that there is no significant difference between assemblages in the open and procedural control treatments (C, O), although both of these differed significantly from assemblages in the shade treatment (S) (Fig. 1.18). Fig. 1.18. PERMANOVA dialog and results for the comparison among the three shade treatments for subtidal epibiota. In the event that the interaction term had been statistically significant, then the logical way to proceed would have been to do pair-wise tests among levels of the Shade factor separately for each of the near and far situations (i.e., within each level of factor Position). To do pair-wise comparisons like this, following on from an interaction term, one would choose (Test $\bullet$Pair-wise test > For term: PositionxShade > For pairs of levels of factor: Shade) in the PERMANOVA dialog. Similarly, one could also in that case consider doing comparisons among levels of the factor Position (near versus far) separately within each of the shade levels (S, O and C). For the latter, one would choose (Test $\bullet$Pair-wise test > For term: PositionxShade > For pairs of levels of factor: Position). In either case, one first chooses the interaction term of interest, followed by the particular factor involved in that interaction which is the focus for that set of pair-wise comparisons. In complex designs, there may be multiple sets of pair-wise comparisons that the user may wish to do in order to follow up the full model partitioning and results shown by the main test. As a final note, keep in mind that the patterns in an MDS (or PCO) ordination plot may or may not show the reasons for an interaction, if it is present, because the full dimensionality of the multivariate cloud has been reduced (usually to 2 or 3 dimensions). Importantly, PERMANOVA works on the underlying dissimilarities themselves for the test, so its results should be trusted over and above any patterns (or lack of patterns) apparent in the ordination. Usually an ordination will help, however, to interpret the PERMANOVA results, provided the stress is not too high24. 23 It was noted earlier that PERMANOVA is not sensitive to differences in correlation structure among groups (Fig. 1.6). Similarly, the PERMANOVA test of interaction is focused more on the additivity of the sizes of the effects, rather than on their direction, per se. Thus, differences in effects which are directional but which do not affect the additivity of effect sizes may go undetected. Studies examining the power of PERMANOVA to detect different kinds of interactions are needed to clarify this issue further. 24 Stress is a measure of how well inter-point distances in an MDS ordination represent the rank-ordered inter-sample dissimilarities in the original resemblance matrix. For a more complete description of the notion of stress, see previous references to MDS along with chapter 5 of Clarke & Warwick (2001) and chapter 7 of Clarke & Gorley (2006) . 1.16 Additivity Central to an understanding of what an interaction means for linear models25 is the idea of additivity. Consider the example of a two-way crossed design for a univariate response variable, where the cell means and marginal means are as shown in Fig. 1.19a. Note that the marginal means are the means of the levels of each factor ignoring the other factor. In an additive model, the difference between two levels of a factor (say between B1 and B2) between individual cells (i.e., within each level of A, that is to say, within each column) are equal to the differences in the marginal means (i.e., the difference between the mean of B1 and B2 if factor A were to be ignored). This can be contrasted with the situation where the differences in cell means are quite different from the differences in marginal means (e.g., Fig. 1.19b), in which case, there is an interaction between the factors. So, this is another way to articulate what is meant by a significant interaction: effects of factors within levels of other factors are non-additive and thus do not match the corresponding shifts in marginal means. The interaction term, in fact, measures the deviation of the cell means we actually got from what we would expect them to be if they were to follow the marginal means, as would be the case if the effects of the two factors were purely additive. Fig. 1.19. Marginal and cell means for a univariate crossed design showing examples of (a) additive effects, (b) multiplicative effects and (c) additivity after log10-transformation of (b). Clearly, the additivity (or not) of the effects of factors is also going to depend on whether or not the data have been transformed (or standardised or ranked) prior to analysis. This is as true for multivariate data as it is for univariate data. For example, if a log (base 10) transformation is applied to the means shown in Fig. 1.19b, then we would have an additive model with no significant interaction (Fig. 1.19c). Such a situation typifies phenomena where the true effects are multiplicative, rather than being additive. In univariate analysis, transformations can often be used to remove significant interaction terms, yielding additivity ( Tukey (1949) , Box & Cox (1964) , Kruskal (1965) , Winsberg & Ramsay (1980) ). For multivariate analysis of ecological data, however, transformations are usually applied neither to fulfill assumptions, nor in order to remove significant interactions, but rather as a method of changing the relative emphasis of the analysis on rare versus more abundant species (e.g., Clarke & Green (1988) , Clarke & Warwick (2001) ). In PRIMER, a blanket transformation can be applied to all variables by choosing Analyse > Pre-treatment > Transform (overall) and then choosing from a range, in increasing severity, from no transformation, square root, fourth root or log(x+1) down to a reduction of the values to binary presence (1) or absence (0). An approach using an intermediate-level transformation (square root or fourth root) has been recommended as a way to reduce the contribution of highly abundant species in relation to less abundant ones in the calculation of the Bray-Curtis measure; rare species will contribute more, the more severe the transformation ( Clarke & Green (1988) , Clarke & Warwick (2001) ). In addition to the transformation, additivity of effects in multivariate analysis is also going to depend on whether or not the dissimilarities are ranked before analysis (yet another reason why patterns in a non-metric MDS, which preserves ranks only, may not necessarily clearly reflect what is given in the PERMANOVA output). The choice of dissimilarity measure itself is also very important here. By performing the partitioning, PERMANOVA is effectively applying a linear model to a multivariate data cloud, as defined by these choices. So the presence of a significant interaction (or not) by PERMANOVA will naturally depend on them. Nevertheless, the choice of an appropriate dissimilarity measure (and also the choice of transformation, if any) should genuinely be driven by the biology and ecology (or other nature) of the system being studied and what is appropriate regarding your hypotheses, and not by reference to these statistical issues (unlike typical traditional univariate ANOVA). 25 The ANOVA models analysed by PERMANOVA are linear only in the space of the multivariate cloud defined by the dissimilarity measure of choice; they are not linear in the space of the original variables (unless the resemblance measure chosen was Euclidean distance). 1.17 Methods of permutations As for the one-way case, the distribution of each of the pseudo-F ratios in a multi-way design is generally unknown. Thus, a permutation test (or some other approach using re-sampling methods) is desirable. When there is more than one factor, situations commonly arise which prevent the possibility of obtaining an exact test of individual terms in the model using permutations. For example, there is no exact permutation test for an interaction (but see Pesarin (2001) , who describes a synchronised permutation method for testing interactions). In addition, restricted permutation methods for testing main effects in ANOVA models generally have low power ( Anderson & ter Braak (2003) ). However, several good approximate permutation methods can be used instead to get accurate P-values (e.g., Anderson & Legendre (1999) ). PERMANOVA provides three general options regarding the method of permutation to be used: (i) unrestricted permutation of raw data, (ii) permutation of residuals under a reduced model, or (iii) permutation of residuals under the full model. These methods and their properties are described in detail elsewhere (e.g., Anderson & Legendre (1999) , Anderson (2001b) , Anderson & Robinson (2001) , Anderson & ter Braak (2003) , Manly (2006) ). Although these methods do not give exact P-values for complex designs in all cases, they are asymptotically exact26 and give very reliable results. In practice, these three approaches will give very similar results, so there is (thankfully) no need to agonise much about making a choice here. All three of the methods are implemented in PERMANOVA so as to ensure that the correct exchangeable units (identified by the denominator of the pseudo-F ratio) are used for each individual test (see Anderson & ter Braak (2003) for details). Some of the known properties of the three methods are outlined below. (i) Unrestricted permutation of raw data. This is a good approximate test proposed for complex ANOVA designs by Manly (1997) . It will generally have type I error close to $\alpha$, although with larger sample sizes it tends to be more conservative (less powerful) than the tests that permute residuals ( Anderson & ter Braak (2003) )27. However, this method does not need large sample sizes to work well ( Gonzalez & Manly (1998) ). It is also, computationally, the fastest option. The method does suffer from a few problems, however, if there happen to be outliers in covariables (if present, see Kennedy & Cade (1996) ), so should not be used for such cases. (ii) Permutation of residuals under a reduced model. This approach was first described for linear models by Freedman & Lane (1983) and is the default option in PERMANOVA because it has excellent empirical and theoretical properties. Empirically, it yields the best power and the most accurate type I error for multi-factorial designs in the widest set of circumstances ( Anderson & Legendre (1999) , Anderson & ter Braak (2003) ). Also, this method is theoretically the closest to the conceptually exact test ( Anderson & Robinson (2001) ). The idea is to isolate the term of interest in the model for each test by fitting the other terms (the reduced model), obtaining residuals from that reduced model, and permuting those. In other words, the entities that are exchangeable under the null hypothesis for a particular term are the errors (estimated by the residuals) obtained after removing the terms in the model that are not of interest for that test. The definition of the reduced model (and therefore the residuals arising from them) thus depends on which term is being tested28. (iii) Permutation of residuals under the full model. This method was described by ter Braak (1992) . The idea is to obtain residuals of the full model by subtracting from each replicate the mean corresponding to its particular cell (combination of factor levels). These residuals are estimating the errors associated with each replicate. These are then permuted and the statistic is re-calculated for all terms using these residuals (as if they were the data) under permutation. This method mostly gives results highly comparable to method (ii). It relies somewhat more than (ii), however, on large within-cell sample sizes for precision. It has the advantage, however, of being faster than method (ii) for the analysis of the entire design, as the same residuals are permuted for all terms in the model under test. In general, we recommend using method (ii), which is the default. Method (i), however, does provide an exact test for the one-way case, so should be used for one-way ANOVA models. Otherwise, note that methods (ii) and (iii) both require estimation of parameters (means) in order to calculate residuals as deviations from fitted values. When sample sizes are small, these estimates are not very precise (i.e., they may not be very close to their “true” values), so the residuals being permuted, in turn, may not be good representatives of the “true” errors ( Anderson & Robinson (2001) ). Thus, in the case of relatively small sample sizes (say, n < 4 replicates per cell), method (i) is also recommended (provided there are no outliers in covariables, as mentioned above). Method (iii) is probably only advisable if you wish to use method (ii), but the time required is getting overly burdensome. 26 An asymptotically exact test is a test for which the type I error (probability of rejecting the null hypothesis when it is true) asymptotically approaches (converges on) the a priori chosen significance level ($\alpha$) with increases in the sample size (N). 27 Note that “less powerful” does not necessarily mean that using (i) will give you a smaller P-value than (ii) or (iii) for any particular data set. It means that, in repeated simulations, the empirical power (estimated probability of rejecting the null hypothesis when it is false) was, on average, smaller for method (i) than for either of the other two methods in most situations. 28 For unbalanced designs, it also depends on which Type of SS is chosen for the test. For Type I SS, the order in which the terms are fitted will also matter here. 1.18 Additional assumptions Recall that we assume for the analysis of a one-way design by PERMANOVA that the multivariate observations are independent and identically distributed (i.i.d.) under a true null hypothesis. When more complex (multi-way) designs are analysed, a few more assumptions are added, due to the way partitioning is done and the method of permutation used (generally, the method used is permutation of residuals under a reduced model). More specifically, for multi-way designs, PERMANOVA fits an additive linear model to the multivariate samples in the space of the chosen resemblance measure. This assumes that effects of factors and their interactions can be modeled meaningfully in this additive fashion, as opposed to using, say, a non-linear, multiplicative or other approach (e.g., see Millar, Anderson & Zunun (2005) ). It also assumes that the errors (which may be estimated using the residuals after fitting a given PERMANOVA model) are i.i.d. across the full design in the space of the chosen resemblance measure. At present, it is unknown to what extent and in what circumstances departures from these assumptions might affect error rates or interpretations of results from PERMANOVA. Plots showing distributions of residuals vs fitted values (e.g., from either PERMANOVA or DISTLM models, see chapter 4) for each of a series of PCO dimensions (see chapter 3) in order of decreasing importance, for example, and also multivariate ordination plots of residuals could provide helpful diagnostic tools. Further study is warranted to develop these tools and to investigate the effects of violations of the assumptions, especially for mixed models, models with covariates or unbalanced designs (which are all discussed in more detail later in this chapter) and even for simple designs that include interactions. In addition, the good asymptotic properties of the methods of permutation of residuals used by PERMANOVA (and DISTLM) require reasonable sample sizes within the cells. Small numbers of observations within cells will result in inter-correlations among residuals. Although exchangeability of errors (implying independence and homogeneity) and additive effects in the space of the resemblance measure are fairly modest assumptions, we nevertheless look forward to future studies where the full implications of the use of these models in different situations and with different resemblance measures may become clearer. In the meantime, simulations done with univariate data having highly non-normal error structures indicate that the permutation methods implemented by PERMANOVA are quite stable for reasonable sample sizes (as indicated above) and when observations are i.i.d. 1.19 Contrasts In some cases, what is of interest in a particular experimental design is not necessarily the comparisons among all pairs of levels of some factor found to be significant, but rather to compare one or more groups (or levels) together versus one or more other groups. A comparison such as this is called a contrast and is usually logically formulated at the design stage (a priori). Pair-wise comparisons, on the other hand, are generally done after the analysis, so are also sometimes called unplanned or a posteriori comparisons. For example, in the study of the subtidal epibiotic assemblages, we may especially wish to compare a priori the assemblages in the shaded treatment with those occurring in either the open or the procedural control treatments. That is, we would like to contrast group (S) with groups (C, O), effectively treating the latter two treatments together as a single group. Another possible contrast of interest would be the comparison of the procedural control with the open group (i.e., C versus O), as this contrast would identify possible artefacts of the structure used to create shade. Fig. 1.20. Dialog to create contrasts for the subtidal epibiota. PERMANOVA allows the user to specify particular contrasts of interest. Each of these has 1 degree of freedom and is used to further partition the sums of squares (SS) attributable to that factor. In addition, any interaction terms involving that factor are also partitioned according to the contrast. To see how this works for the subtidal data, click on the design file named Two-way crossed in the sub.pwk workspace and choose Tools>Duplicate, then choose File>Rename Design and name this new design file With Contrasts for reference. Next, in row 2 for the Shade factor, double click on the cell in the final column, headed ‘Contrasts’, which will bring up a new ‘Contrasts’ window (Fig. 1.20). Click on ‘Add’ in order to add a new row to the list of contrasts shown in the window. Double click within the cell in the first column (‘Name’) and type the name S-vs-(C,O). Double click in the cell in the second column (‘Contrast’). This brings up a dialog in which all of the levels available in that factor are listed. Identify the contrast by clicking on S, then on to place it on the left and then click on each of C and O in turn, each followed by , to place the latter two levels on the right, followed by ‘OK’. Next, add a second row in the contrasts window which specifies a contrast of the C and O treatments, to be called C-vs-O (Fig. 1.20). When you are finished, the design file will show the names of the contrasts you have specified in the ‘Contrasts’ column for the Shade factor (Fig. 1.20). Re-run the PERMANOVA analysis by selecting the resemblance matrix and choosing PERMANOVA+ > PERMANOVA > (Design worksheet: With Contrasts) & (Test $\bullet$Main test) & (Sums of Squares $\bullet$Type III (partial)) & (Permutation method $\bullet$Permutation of residuals under a reduced model) & (Num. permutations: 9999). The contrasts associated with a particular factor are indented in the output file, to emphasise that these are a further partitioning of the SS associated with a given factor (or interactions involving that factor) (Fig. 1.21). These are therefore “extra” tests, in addition to the test of the factor as a whole. In general, one may construct up to $df _A$ orthogonal (independent) contrasts, where $df _A$ is the number of degrees of freedom for the factor. For these situations (i.e., orthogonal a.k.a. independent contrasts), then the SS for the contrasts chosen will add up to the SS for the factor (that is, if the number of orthogonal contrasts is equal to $df _A$) and it is generally considered unnecessary to worry about doing any corrections for multiple tests. This is true of the example we have done. However, if the contrasts are not orthogonal, and especially if one has chosen to perform a great many contrasts that exceed $df _A$ in number, then the potential inflation of overall type I error due to performing multiple tests may be an issue (see the section on Pair-wise comparisons). For our example (Fig. 1.21), the results reveal the strong effect of the shading treatment versus the other two treatments (note how S-vs-(C,O) has a pseudo-F ratio even larger than the pseudo-F for the Shade factor overall) and also show the lack of any apparent procedural artefact (C-vs-O, P = 0.28). In addition, the interaction term of Position × S-vs-(C,O) is approaching statistical significance (P = 0.05), providing increased evidence that the effect of shading (which is tested more directly by this specific contrast than by the test comparing all three treatments) is different near the seafloor compared to far away from the seafloor (e.g., Glasby (1999) ). Fig. 1.21. Results for subtidal epibiota, including specified contrasts. A special situation where contrasts might be useful is in the comparison of, say, an impact location versus several control locations in an environmental impact study design ( Underwood (1992) ). If there is no replication of the impact state, but there are several control locations, then this actually generates what is known as an asymmetrical design. For such designs, the correct SS can be obtained using contrasts, however, it is necessary to treat the model explicitly as an asymmetrical design in order to get correct pseudo-F ratios and P-values for all of the tests (see the specific section on Asymmetrical designs). 1.20 Fixed vs random factors (Tasmanian meiofauna) All of the factors considered so far have been fixed, but factors can be either fixed or random. In univariate ANOVA the choice of whether a particular factor is fixed or random has important consequences for the assumptions underlying the model, the expected values of mean squares (thus, the construction of a given F ratio), particularly in more complex designs, the hypothesis tested by the F ratio and, perhaps most importantly, the extent and nature of the inferences. This is also true for PERMANOVA, which follows the analogous univariate ANOVA models in terms of the construction of pseudo-F ratios from expectations of mean squares (EMS). For a fixed factor, there is a finite set of levels which have been explicitly chosen to represent particular states. These states (the levels) are generally explicit because they have been manipulated (e.g., shade, open and control), because they already exist in nature (e.g., male and female), or because they bear some meaningful relationship to other chosen levels (e.g., high, medium and low). All of the levels of a fixed factor (or at least all of the ones we are effectively interested in) occur in the experiment. So, for a fixed factor, individual levels have a meaning in themselves, and generally, if we were going to repeat the experiment, we would choose the same levels to investigate again (unless we were going to change the hypothesis). The effects of a fixed factor are deemed to be constant values for each nominated state (or level)29. Furthermore, the component of variation attributable to a fixed factor in a given model is considered in terms of the sum of squared fixed effects (divided by the appropriate degrees of freedom). Finally, when we perform the test of the fixed factor, the resulting P-value and any statistical inferences to be drawn from it apply only to those levels and to no others. For a random factor, however, the particular levels included in the experiment are a random subset from a population of possible levels we could have included (e.g., sites, locations, blocks). We do not consider the individual levels (site 1, site 2, etc.) as representing any particular chosen state. The levels do not have any particular meaning in themselves; we are not interested in comparing, say, site 1 vs site 2, but rather multiple levels of the factor are included to provide us with a measure of the kind of variability that we might expect across the population of possible levels (e.g., among sites). Random factors, unlike fixed factors, actually contribute another source of random variance into the model (in addition to the error variance) and so the effects, rather than being fixed, are instead used to estimate the size of the variance component for that factor. Repeating the experiment would also probably not result in the same levels being chosen again. Importantly, when the F ratio is constructed and the P-value is calculated for a random factor, the statistical inference applies to the variance component for the whole population of possible levels (or effects) that could have been chosen, and not just to the levels included in the experiment. The difference in the inference space for the fixed versus the random factor is important. For example, if one obtained a P-value less than $\alpha = 0.05$ (the usual convention) for a fixed factor, one might state: “There were significant differences among (say) these three treatments in the structure of their assemblages”. In contrast, for a random factor, one might state: “There was significant variability among sites in the structure of the assemblages.” The fixed factor is more specific and the random factor is more general. Note also that the statement regarding the fixed factor often logically calls for more information, such as pair-wise comparisons between individual levels – which treatments were significantly different from one another and how? Whereas, the second statement does not require anything like this, because the individual levels (sites) are generally not of any further interest in and of themselves (we don’t care whether site 1 differs from site 3, etc.); it is enough to know that significant variability among sites is present. Thus, it is generally not logical to do pair-wise comparisons among levels of a random factor (although PERMANOVA will not prevent you from doing such tests if you insist on doing them)! A logical question to ask, instead, following the discovery of a significant random factor would be: how much of the overall variability is explained by that factor (e.g., sites)? For a random factor, we are therefore more interested to estimate the size of its component of variation and to compare this with other sources of variation in the model, including the residual (see Estimating components of variation). To clarify these ideas, consider an example of a two-way crossed experimental design used to study the effects of disturbance by soldier crabs on meiofauna at a sandflat in Eaglehawk Neck, Tasmania ( Warwick, Clarke & Gee (1990) ). The N = 16 samples consist of two replicates within each combination of four blocks (areas across the sandflat) and two natural ‘treatments’ (either disturbed or undisturbed by soldier crab burrowing activity). There were p = 56 taxa recorded in the study, consisting of 39 nematode taxa and 17 copepod taxa. We therefore have Factor A: Treatment (fixed with a = 2 levels, disturbed (D) or undisturbed (U) by crabs). Factor B: Block (random with b = 4 levels, labeled simply B1-B4). The first factor is fixed, because these two treatment levels do have a particular meaning, representing particular states in nature (disturbed or undisturbed) that are of interest to us and that we explicitly wish to compare. In contrast, the second factor in the experiment, Blocks, is random, because we do not have any particular hypotheses concerning the states of, say B1 or B2, but rather, these are included in the experimental design in order to estimate variability across the sandflat at a relevant spatial scale (i.e., among blocks) and to avoid pseudo-replication (sensu Hurlbert (1984) ) in the analysis of treatment effects. A design which has both random and fixed factors is called a mixed model. This particular design (where a fixed factor is crossed with a random one), also allows us to test for generality or consistency in treatment effects across the sandflats. The data are located in the file tas.pri in the ‘TasMei’ folder of the ‘Examples add-on’ directory. View the factors by choosing Edit>Factors. In order to obtain different symbols for each of the a × b = 2 × 4 = 8 cells in the design, create a new factor whose 8 levels are all combinations of Treatment × Block by choosing ‘Combine’ in the Factors dialog box and then followed by ‘OK’. Create a resemblance matrix on the basis of Bray-Curtis similarities after square-root transformation. An MDS of these data (as shown in chapter 6 of Clarke & Warwick (2001) ) shows a clear effect of disturbance on these assemblages, with some variability among the blocks as well (Fig. 1.22). Fig. 1.22. MDS of Tasmanian meiofauna showing variation among the blocks on the x-axis and separation of disturbed from undisturbed communities on the y-axis. Set up the design file according to the correct experimental design (Fig. 1.23) and rename it as Mixed model, for reference. The next step is to run PERMANOVA on the resemblance matrix according to the mixed model design. Due to the small number of observations per cell, choose (Permutation method ·Unrestricted permutation of raw data) & (Num. permutations: 9999). Save the workspace with the PERMANOVA results as tas.pwk. The PERMANOVA results show that there is a statistically significant (though borderline) interaction term, indicating that the treatment effects vary from one block to the next (P = 0.044, Fig. 1.23). The MDS plot shows that the direction of the treatment effects appears nevertheless to be fairly consistent across the blocks (at least insofar as this can be discerned using an ordination to represent the higher-dimensional cloud of points), so in this case the significant interaction may be caused by there being slight differences in the sizes of the treatment effects for different blocks. 29​ In univariate analysis, the effect for a given group is defined as the deviation of the group mean from the overall mean. Similarly, when using PERMANOVA for multivariate analysis, the effect for a given level is the deviation (distance) of the group’s centroid from the overall centroid in the multivariate space, as defined by the dissimilarity measure chosen. 1.21 Components of variation For any given ANOVA design, PERMANOVA identifies a component of variation for each term in the model, denoted as ‘S(*)’ for the fixed terms and ‘V(*)’ for the random terms. This notation is used because, in the analogous univariate case, components of variation for a fixed factor are sums of squared fixed effects (divided by appropriate degrees of freedom), while components of variation due to random factors are actual measures of variability or variance components. This is appropriate, because the hypotheses are also different in these two cases. For fixed effects, the hypothesis only concerns the effects of those levels that were included in the experiment, whereas for the random factor, the hypothesis is about the variability among a population of levels, of which the levels in the experiment are a random representative sample. Note that any interaction involving a random factor will also be random. The residual component is also denoted by ‘V(Res)’ because it too is a measure of variability. Thus, for example, in the two-way mixed model design for the Tasmanian meiofauna, the total variation is partitioned according to four sources, as follows: S(Tr): sum of squared treatment effects (divided by degrees of freedom); V(Bl): variation due to blocks; V(TrxBl): variation in treatment effects among blocks; and V(Res): residual variation. Note that if one is using PERMANOVA to analyse a single variable with Euclidean distance, then the components of variation are indeed true variance components (for random factors) and true sums of squared fixed effects divided by degrees of freedom (for fixed factors). Here and in PERMANOVA, they are simply called components of variation, however, in order to cover the more general cases, including the analysis of many variables on the basis of non-Euclidean resemblance measures. Fig. 1.23. Design file and PERMANOVA results for the Tasmanian meiofauna, including details of the expected mean squares and construction of F ratios for each term in the mixed model 1.22 Expected mean squares (EMS) An important consequence of the choice made for each factor as to whether it be fixed or random is identified by examining the expected mean squares (EMS) for each term in the resulting model. This is vitally important because the EMS’s are used to identify an appropriate denominator mean square that one must use for each particular term in the model in order to construct a correct pseudo-F ratio that will isolate that term of interest for the test. For univariate ANOVA, the underlying theory for deriving expectations of sums of squares and mean squares is covered well elsewhere (e.g., Cornfield & Tukey (1956) , Hartley (1967) , Rao (1968) , Winer, Brown & Michels (1991) , Searle, Casella & McCulloch (1992) ). PERMANOVA actually uses the same “rules” for constructing these expectations, implementing these as a direct multivariate analogue to the univariate approach. The default rules used by PERMANOVA assume that fixed effects sum to zero, following Cornfield & Tukey (1956) and Winer, Brown & Michels (1991) . Some constraint on fixed effects is necessary, due to the intrinsic over-parameterisation of the ANOVA model (e.g., Scheffé (1959) ), and the sum-to-zero constraint is a highly convenient one. The constraint chosen does not affect the sums of squares. It does, however, affect the EMS’s and thus the F ratios and P-values that are obtained for ANOVA mixed models. Some well-known computer packages for univariate statistics relax the sum-to-zero constraint for fixed effects in mixed interactions, including SPSS and the ‘proc GLM’ routine in SAS. To obtain EMS’s in accordance with these packages, remove the $\checkmark$ in front of the ‘Fixed effects sum to zero’ box in the PERMANOVA dialog. See Hartley & Searle (1969) , Searle (1971) , Hocking (1973) , McLean, Sanders & Stroup (1991) and Searle, Casella & McCulloch (1992) for further discussion and debate regarding this issue. Once the individual components of variation have been identified, then the expectations of the mean squares for each term in the model are determined precisely in terms of these components, and are provided in the output under the heading ‘Details of the expected mean squares (EMS) for the model’. Thus, we can see, for example (Fig. 1.23), that the expectation for the mean square calculated for the term ‘Block’ (abbreviated as ‘Bl’ in the output) is one times the residual variation plus four times the variation due to blocks (denoted as ‘1*V(Res) + 4*V(Bl)’ in the output). Each term in the model will have an expected mean square that consists of some linear combination of the components of variation in the model. 1.23 Constructing $F$ from EMS The determination of the EMS’s gives a direct indication of how the pseudo-F ratio should be constructed in order to isolate the term of interest to test a particular hypothesis (e.g., Table 1.1). Consider the test of the term ‘Block’ (‘Bl’) in the above mixed-model experimental design. The null hypothesis is that there is no significant variability among blocks. Another way of writing this, in the notation used by PERMANOVA, is H$_0$: V(Bl) = 0. To begin, we will construct an F ratio where the numerator is the mean square for the term of interest (e.g., to test blocks, then the numerator will be the mean square for blocks). Now, given the EMS for this numerator term of interest, we essentially need to find a denominator whose expectation would correspond to the numerator if the null hypothesis were true. In the case of the ‘Block’ term, if V(Bl) were equal to zero, then the EMS for blocks would just be 1*V(Res). The term whose EMS corresponds to 1*V(Res) is the residual. Therefore, the pseudo-F ratio for the test of the ‘Block’ term is the mean square for blocks divided by the residual mean square, viz: $F _ {Bl} = MS _ {Bl} / MS _ {Res}$. By following a similar logic, we can see that the test of the interaction term (‘TrxBl’) is provided by $F _ {Tr \times Bl} = MS _ {Tr \times Bl} / MS _ {Res}$ (Table 1.1). For the ‘Treatment’ term (‘Tr’), we see that its EMS is: 1*V(Res) + 2*V(TrxBl) + 8*S(Tr). The null hypothesis here is that there are no consistent treatment effects, or, equivalently, $ \text{H} _ 0 $: S(Tr) = 0. If the null hypothesis were true, then the EMS for treatments would be: 1*V(Res) + 2*V(TrxBl). The term whose EMS corresponds to this is the interaction term ‘TrxBl’. Therefore, the pseudo-F ratio for the test of no treatment effects is the mean square for treatments divided by the interaction mean square, i.e., $F _ {Tr} = MS_ {Tr} / MS_ {Tr \times Bl}$. Clearly, it would be incorrect to construct $F _ {Tr} = MS_ {Tr} / MS_ {Res}$, because then, if V(TrxBl) were non-zero, we might easily reject $\text{H} _0$: S(Tr) = 0 even if it were true30. In the PERMANOVA output, the terms which provide the numerator and denominator mean squares for each pseudo-F ratio are provided in the output under the heading ‘Construction of Pseudo-F ratio(s) from mean squares’ (Fig. 1.23). Also given here are the degrees of freedom associated with the numerator (‘Num.df’) and denominator (‘Den.df’) terms. Most of the time, the multipliers on these mean squares will simply be 1, and a single term can be found to provide an appropriate denominator mean square for each of the relevant hypotheses in the model. In some cases, however, a linear combination of mean squares must be sought in order to construct correct pseudo-F ratios to test particular terms (see the section Linear combinations of mean squares). PERMANOVA can deal with these situations and, in such cases, details of the linear combinations used are also provided. Table 1.1. Null hypothesis, construction of pseudo-F, and ratio of expectations for each term in the mixed model for the study of Tasmanian meiofauna. Note how the construction of pseudo-F isolates the component of interest under the null hypothesis (circled) so that, in each case, the numerator and denominator expectations will match one another if the null hypothesis were true. 30 Some might argue that the existence of an interaction necessarily implies the existence of at least some non-zero treatment effect(s), but we consider the null hypothesis for the test of the main effect in a mixed model such as this to include the concept of consistency in treatment effects, rather than simply the notion of whether there are any treatment effects at all. See the section on Inference space and power. 1.24 Exchangeable units The denominator mean square of the pseudo-F ratio for any particular term in the analysis is important not just because it isolates the component of interest in the numerator for the test: it also identifies the exchangeable units needed to obtain a correct test by permutation ( Anderson & ter Braak (2003) ). In the one-way case, it is clear that the units that are exchangeable under a true null hypothesis are the individual samples. These can be shuffled randomly among the groups (or, alternatively, the group labels can be randomly shuffled across all samples) if the groups have no effect and the null hypothesis is true. In fact, whenever a term has a pseudo-F ratio with the residual mean square as its denominator, then the exchangeable units for the test are the individual samples themselves (regardless of which of the three methods of permutation offered by PERMANOVA is to be employed). For more complex designs, the correct exchangeable units for a given test are identified by the term used as the denominator mean square of that particular term. This is sensible from the perspective that the denominator identifies what constitute the “errors” for a given null hypothesis. Thus, in the Tasmanian meiofauna example, the pseudo-F ratio for the test of the main effect of treatments (‘Tr’) is $F _ {Tr} = MS _ {Tr} / MS _ {Tr \times Bl}$. As the denominator here is the interaction term ‘TrxBl’, the exchangeable units for this test are the 8 cells that correspond to the 2 × 4 combinations of treatments by blocks. The samples within each of those 8 cells will be kept together as a unit under permutation. This yields 8! / [(2!)$^4$ × 4!] = 105 unique values of the numerator and 8! / [4! × 2!] = 840 unique values for the denominator and, thus, 840 unique values of the whole pseudo-F statistic under permutation (as shown in the ‘Unique perms’ column for the term ‘Tr’ in Fig. 1.23 above). For more details concerning exchangeable units for permutation tests in ANOVA designs, see Anderson & ter Braak (2003) . 1.25 Inference space and power It is worthwhile pausing to consider how the above tests correspond to meaningful hypotheses for the mixed model. What is being examined by $F _ {Tr}$ is the extent to which the sum of squared fixed effects can be detected as being non-zero over and above the potential variability in these effects among blocks, i.e., over and above the interaction variability (if present). Thus, in a crossed mixed model like this, $F _ {Tr \times Bl}$ first provides a test of generality (i.e., do the effects of treatments vary significantly among blocks?), whereas the test of the main effect of treatments ($F _ {Tr}$) provides a test of the degree of consistency. In other words, even if there is variation in the effects of treatments (i.e. V(TrxBl) ≠ 0), are these consistent enough in their size and direction that an overall effect can be detected over and above this? Given the pattern shown in the MDS plot, it is not surprising to learn that, in this case, the main treatment effect is indeed discernible over and above the variation in its effects from block to block (P = 0.0037). We can contrast these results with what would have happened if we had done this analysis but treated the blocks as fixed instead of random (Fig. 1.24). The consequence of this choice is a change to the EMS for the ‘Treatment’ term, which is now 1*V(Res) + 8*S(Tr) and therefore contains no component of variation for the interaction. As a consequence, the pseudo-F ratio for treatment main effects in this fully fixed model is constructed as $F _ {Tr} = MS _{Tr} / MS _ {Res}$ and, correspondingly, the denominator degrees of freedom for this test has increased from 3 to 8, the value of the pseudo-F ratio has changed from 4.67 to 8.08 and the P-value has decreased (cf. Figs. 1.23 & 1.24). Fig. 1.24. Design file and analysis of Tasmanian meiofauna, treating the ‘Blocks’ as fixed. On the face of it, we appear to have achieved a gain in power by using the fully fixed model in this case, as opposed to the mixed model that treated ‘Blocks’ as random. Power is the probability of rejecting the null hypothesis when it is false. When using a statistic like pseudo-F (or pseudo-t), power is generally increased by increases in the denominator degrees of freedom. Basically, the more information we have about a system, the easier it is to detect small effects. This choice of whether to treat a given factor as either fixed or random, however, doesn’t just affect the potential power of the test, it also rather dramatically affects the nature of our hypotheses and our inference space. If we choose to treat ‘Blocks’ as random (Fig. 1.23), then: (i) the test of ‘TrxBl’ is a test of the generality of disturbance effects across blocks; (ii) the test of ‘Bl’ is a test of the spatial variability among blocks; and (iii) the test of ‘Tr’ is a test of the consistency in treatment effects, over and above the potential variability in its effects among blocks. Importantly, the inference space for each test refers to the population of possible blocks from which we could have sampled. In contrast, if we choose to treat ‘Blocks’ as fixed (Fig. 1.24), then: (i) the test of ‘TrxBl’ is a test of whether disturbance effects differ among those four particular blocks; (ii) the test of ‘Bl’ is a test of whether there are any differences among those four particular blocks; and (iii) the test of ‘Tr’ is a test of treatment effects, ignoring any potential variation in its effects among blocks. Also, the inference space for each of these tests in the fully fixed model refers only to those four blocks included in our experiment and no others. In the end, it is up to the user to decide which hypotheses are the most relevant in a particular situation. The choice of whether to treat a given factor as fixed or random will dictate the EMS, the pseudo-F ratio, the extent of the inferences and the power of the tests in the model. Models with mixed and random effects will tend to have less power than models with only fixed effects. (Consider: in the context of the example, it would generally be easier to reject the null hypothesis that treatments have no effect whatsoever than it would be to reject the null hypothesis that treatments have no effect given some measured spatial variability in those effects.) However, random and mixed models can provide a much broader and therefore usually a more meaningful inference space (e.g., extending to the wider population of possible blocks across the sampled study area, and not just to those that were included in the experiment). Such models will therefore tend to correspond to much more logical and ecologically relevant hypotheses in many situations. As we have seen, PERMANOVA employs direct multivariate analogues to the univariate results for the derivation of the EMS and construction of the pseudo-F ratio, so all of the well-known issues regarding logical inferences in experimental design that would occur for the univariate case (e.g., Cornfield & Tukey (1956) ; Underwood (1981) ; Hurlbert (1984) ; Underwood (1997) ) necessarily need to be considered for any multivariate analysis to be done by PERMANOVA as well. As a final note regarding inference, the traditional randomization test is well known to be conditional on the order statistics of the data (e.g. Fisher (1935) ). In other words, the results (P-values) depend on the realised data values. This fact led Pitman ( Pitman (1937a) , Pitman (1937b) , Pitman (1937c) ) and Edgington (1995) to argue that the inferences from any test done using a randomization procedure can only extend to the actual data themselves, and can never extend to a wider population31. In a similar vein, Manly (1997) , section 7.6, stated that randomization tests must, by their very nature, only allow factors to be treated as fixed, because “Testing is conditional on the factor combinations used, irrespective of how these were chosen” (p. 142). PERMANOVA, however, uses permutation methods in order to obtain P-values, but it also (clearly) allows factors to be treated as random. How can this be? Fisher (1935) considered the validity of a permutation test to be ensured by virtue of the a priori random allocation of treatments to individual units in an experiment. That is, random allocation before the experiment justifies randomization of labels to the data afterwards in order to create alternative possible outcomes we could have observed. However, we shall consider that a permutation test gains its validity more generally (such as, for example, in observational studies, where no a priori random allocation is possible), by virtue of (i) random sampling and (ii) the assumption of exchangeability under a true null hypothesis. As stated earlier (see the section Assumptions), PERMANOVA assumes only the exchangeability of appropriate units under a true null hypothesis. Random sampling of the levels of a random factor from a population of possible levels (like random sampling of individual samples), coupled with the assumption that “errors” (whether they be individual units or cells at a higher level in the design) are independent and identically distributed (“i.i.d.”) (e.g., Kempthorne (1952) ) ensures the validity of permutation tests for observational studies. For further discussion, see Kempthorne (1966) , Kempthorne & Doerfler (1969) and Draper, Hodges, Mallows et al. (1993) . Thus, provided (i) the population from which levels have been chosen can be conceived of and articulated clearly, (ii) the exchangeability of levels can be asserted under a true null hypothesis and (iii) random sampling has been used, then the permutation test is valid for random factors, with an inference space that logically extends to that population. In section 7.3 of the more recent edition of Manly’s book ( Manly (2006) ), the validity of extending randomization procedures to all types of analysis of variance designs, including fixed and random factors, is also now acknowledged on the basis of this important notion of exchangeability ( Anderson & ter Braak (2003) ). 31 To be fair, Edgington (1995) also suggested that making such wider inferences using the normal-theory based tests was almost always just as much a “leap of faith” as it would be for a randomization test, due to the unlikely nature of the assumptions required by the former. 1.26 Testing the design Given the fact that so many important aspects of the results (pseudo-F ratios, P-values, power, the inference space, etc.) depend so heavily on the experimental design (information given in the design file), one might wish to examine the qualities of various designs, even before embarking on the serious task of actually gathering the data. In PERMANOVA+, a special routine is provided that allows the user to explore different designs without actually analysing data. All that is required to perform a test of a given design is (i) a (dummy) data worksheet file possessing relevant factor information and (ii) a design file. The dummy data file with factor information is needed in order to identify the number of replicates per cell and how many levels there are for all of the factors. For example, highlight the design file Mixed model in the tas.pwk file. Choose PERMANOVA+ > Test design > (Dummy data worksheet: tas) & (Sums of Squares $\bullet$Type III (partial)). The results will include details of the expected mean squares (EMS) for each term in the model and the construction of the pseudo-F ratio for each term, including the degrees of freedom for each test (Fig. 1.25). This allows the user to trial different experimental designs and to consider the best options, for relevant statistical inferences and for assessing power (on the basis of denominator degrees of freedom) for given hypotheses of interest. For example, in the output for the Tasmanian meiofauna, one can see that the construction for the pseudo-F ratio for ‘Tr’ (the most important test of interest to the experimenter here) is $F _{Tr} = MS _ {Tr} / MS _ {Tr \times Bl}$. Therefore, a more powerful experimental design for the test of disturbance effects (the term ‘Tr’) would actually be obtained by increasing the number of blocks in the experiment (and thus, the degrees of freedom associated with the term ‘TrxBl’), rather than increasing the number of replicates per block (residual degrees of freedom), which would be unlikely to have any immediate effect on the test of $F _ {Tr}$. Fig. 1.25. Testing the two-way mixed model design for the Tasmanian meiofauna. 1.27 Nested design (Holdfast invertebrates) We have seen how a crossed design is identifiable by virtue of every level of one factor being present in every level of the other factor, and vice versa (e.g., Fig. 1.14). We can contrast this situation with a nested design. A factor is nested within another (upper-level) factor if its levels take on different identities within each level of that upper-level factor. For example, consider a study by Anderson, Diebel, Blom et al. (2005) of the spatial variability in assemblages of invertebrates colonising holdfasts32 of the kelp Ecklonia radiata in northeastern New Zealand. An hierarchical sampling design was used (Fig. 1.26). Divers collected $n = 5$ individual holdfasts (separated by metres) within each of 2 areas (separated by tens of metres) within each of 2 sites (separated by hundreds of metres to kilometers), within each of 4 locations (separated by hundreds of kilometers) along the coast. The design was balanced and fully nested with three factors: Factor A: Locations (random with a = 4 levels). Factor B: Sites (random with b = 2 levels, nested in Locations). Factor C: Areas (random with c = 2 levels, nested in Sites and Locations) A design with a series of nested terms, like this one, is sometimes also called a fully hierarchical design. Such designs are especially useful for describing patterns and estimating variability at different temporal or spatial scales ( Andrew & Mapstone (1987) , Underwood, Chapman & Connell (2000) ). Consideration of a schematic diagram for this design (Fig. 1.26) indicates directly how it differs from the crossed designs seen earlier (cf. Fig. 1.14). Unlike the crossed design, we cannot swap the order of the factors. The fact that sites are nested within locations means that we are obliged to consider them in this order: locations first (at the top of the diagram), then sites within locations, and so on. Furthermore, the sites at the first location have nothing to do with the sites at the second location. Individual site levels actually belong to particular location levels – an important hallmark of a nested factor. Finally, for a nested design, the fact that particular levels of factor B belong to particular levels of factor A indicates that the individual levels of factor B do not have a particular meaning in and of themselves. For example, the specific identity or meaning of site 1 at location 1 in the design is not the same as that of site 3 at location 2. Instead, it is clear that the levels of factor B are used to measure variability at the correct spatial (or temporal) scale in order to test the upper-level factor in the design (e.g., Hurlbert (1984) ). Thus, any nested factor is also, necessarily, random. An upper-level factor, however, may be either fixed or random (in the present example, the ‘Locations’ factor happens to be random). Fig. 1.26. Schematic diagram of the sampling design for holdfast invertebrates. The data from this example are located in the file hold.pri in the ‘HoldNZ’ folder of the ‘Examples add-on’ directory. There were p = 351 variables recorded from a total of N = a × b × c × n = 80 holdfasts. Most of the variables were counts of abundances, but some species or taxa (primarily encrusting forms, such as sponges, bryozoans and ascidians) were recorded using a subjective ordinal rating (0 = absent, 1 = rare, 2 = present, 3 = common). For this example, we shall begin by concentrating on data obtained for molluscs only. Open the file and choose Select > Variables > (Indicator levels > Indicator name: Phylum), then click on the ‘Levels…’ button. Next, in the ‘Selection’ dialog, click on to move all of the phyla into the ‘Available’ box, then click on the name Mollusca, followed by to place it into the ‘Include’ box and click ‘OK’ (Fig. 1.27). Now with the molluscan species selected, choose Tools > Duplicate to obtain these in their own separate worksheet and rename this molluscs for reference. For this new worksheet, select Edit > Properties and change the ‘Title’ to NZ kelp holdfast molluscs. Fig. 1.27. Dialog for selection of a subset of variables (molluscs only) for the holdfast data. The objective here is to partition the variability in the species composition of molluscs according to the three-factor hierarchical experimental design. With our focus, for the moment, on composition alone (presence/absence), we will base the analysis on the Jaccard measure, which (when expressed as a dissimilarity) is directly interpretable as the percentage of unshared species between two sample units. We wish to determine if there is significant variability among areas, among sites and among locations in the composition of the molluscan assemblage. Furthermore, if significant variability is detected at any of these levels, it would then be logical (and interesting) to estimate and compare the sizes of each component of variation, which correspond to these different spatial scales. Click on the mollusc worksheet, then select Analyse > Resemblance > (Analyse between •Samples) & (Measure •More (tab)), click on the tab at the top of the dialog window that is labeled ‘More’, then choose (•Similarity P/A > S7 Jaccard) and click ‘OK’. Next, we need to create an appropriate design file for the analysis (Fig. 1.28). Highlight the resulting Jaccard resemblance matrix and select PERMANOVA+ > Create PERMANOVA design > (Title: NZ kelp holdfast molluscs) & (Number of factors: 3). With the design file open, in the first column, place the name of each factor into its appropriate row by double-clicking inside the cell (Location in row 1, Site in row 2 and Area in row 3). Next, in column 3 of the design file, specify that each of these factors is ‘Random’. Now, you will need to specify explicitly the nesting relationships among the factors. To specify that sites are nested within locations, double click inside the cell in row 2 (for Sites) and column 2 (headed ‘Nested in’). This will bring up a new dialog that allows you to choose the factor within which ‘Sites’ are nested. Click on Location in the ‘Available’ box, followed by to move it over into the ‘Include’ box, then click ‘OK’. Specify also that Area is nested in Site, using the same approach (Fig. 1.28)33. When you are finished specifying the design, rename the design file Nested design and save the workspace created so far as hold.pwk. Fig. 1.28. Creating the PERMANOVA design for the New Zealand holdfast data, including nesting. To run the analysis, highlight the resemblance matrix, choose PERMANOVA+ > PERMANOVA > (Design worksheet: Nested design) & (Num. permutations: 9999), leaving all other options as the defaults. The results show significant variability at each level in the design (Fig. 1.29). The EMS’s reveal the rationale for constructing correct pseudo-F ratios for each term in the model (Fig. 1.29). Fig. 1.29. Design file and PERMANOVA analysis of variability in mollusc composition (based on the Jaccard measure) from kelp holdfast assemblages. 32 A holdfast is a root-like structure at the base of the kelp that holds it to the substratum, which is usually a rocky reef. A great diversity of invertebrates inhabit the interstices of these complex structures. 33 PERMANOVA will do the correct analysis here if we choose to nest Areas within Sites only (because Sites have already been specified as being nested in Locations), as shown, or if we choose to specify explicitly that Areas are nested in Sites and also in Locations. 1.28 Estimating components of variation The EMS’s also yield another important insight: they provide a direct method to get unbiased estimates of each of the components of variation in the model. PERMANOVA estimates these components using mean squares, in a directly analogous fashion to the unbiased univariate ANOVA estimators of variance components (e.g., Searle, Casella & McCulloch (1992) ). In essence, this is achieved by setting the mean squares equal to their expectations and solving for the component of interest. For example, by setting $MS _ {Ar}$ and $MS _ {Res}$ equal to their respective expectations (placing “hats” on the parameters to indicate that we are now talking about estimates of these things, rather than their true parameter values), we have:   $ MS _ {Ar} = 1 \times \hat{\text{V}} \left( \text{Res} \right) + 5 \times \hat{\text{V}} \left( \text{Ar(Si(Lo))} \right) $   $ MS _ {Res} = 1 \times \hat{\text{V}} \left( \text{Res} \right)$ Thus,   $ \hat{\text{V}} \left( \text{Res} \right) = MS _ {Res} /1 $   $ \hat{\text{V}} \left( \text{Ar(Si(Lo))} \right) = \left( MS _ {Ar} - MS _ {Res} \right) / 5 $ From the output, we therefore can calculate these estimates directly from the mean squares. The estimated component of variation for the residual is $MS _{Res} = 2525.7$ and for areas this is $( MS _ {Ar} – MS _ {Res} ) / 5 = (3111.8 – 2525.7) / 5 = 117.2$. Similar logic, when applied to the other terms in the analysis yields:   $ \hat{\text{V}} \left( \text{Si(Lo)} \right) = \left( MS _ {Si} - MS _ {Ar} \right) / 10 = 110.9 $   $ \hat{\text{V}} \left( \text{Lo} \right) = \left( MS _ {Lo} - MS _ {Si} \right) / 20 = 381.7 $ These estimates are all calculated automatically by the program and included in the output file in the column labeled ‘Estimate’ under the heading entitled ‘Components of variation’ (Fig. 1.29). For the species composition of molluscs in these kelp holdfast assemblages, the greatest component of variation occurred at the smallest spatial scale (the residual), followed by locations, and then areas and sites, with the latter two being comparable in size (Fig. 1.29). An important point here is that these estimates are not actual “variance components” in the traditional sense unless one is analysing a single variable and the resemblance measure used is Euclidean distance. In addition, these are obviously not the same as variance-covariance matrices used in traditional multivariate statistics either (e.g. Mardia, Kent & Bibby (1979) , Seber (1982) ), because they do not include any estimation of covariance structure at all. Rather, they are interpretable geometrically as measures of variability from a partitioning on the basis of the dissimilarity (or similarity) measure chosen. These estimates (like their univariate counterparts of variance components) will be in terms of the squared units of the dissimilarity measure chosen. Thus, in order to put these back onto the original units, PERMANOVA also calculates their square root (provided in the column labeled ‘Sq.root’ in the results file). These values are akin to a standard deviation in a traditional univariate analysis. Thus, if the value of the dissimilarity measure used has a direct interpretation (such as the Jaccard or Bray-Curtis measures, which are both percentages), then these can be examined and interpreted as well. For example, the greatest variation in molluscan composition is at the level of individual replicate holdfasts, which (according to the square root of the estimated component of variation due to the residual of 50.3) may share only around 50% of their species, even though they may be separated by just a few metres. Over and above this, holdfasts in different areas may be an additional 10-11% dissimilar in their composition, on average, and so on. Although the above design included all random factors, for which a discussion of estimating components of variation is a fairly natural one, we can also estimate the components of variation due to fixed effects. Recall that these are not measures of variance per se, but rather are sums of squared fixed effects divided by appropriate degrees of freedom (see the section on Components of variation). However, if we are interested in comparing the amount of variation that is attributable to different terms in the model, estimates of components for fixed and/or random factors are useful and are directly comparable. In fact, it is indeed these estimates of components of variation that should be used as a correct basis for comparing the relative importance of different terms in the model towards explaining overall variation ( Underwood & Petraitis (1993) ). In contrast, the raw sums of squares (whether alone or as a percentage of the total sum of squares) are not directly comparable, because different terms generally have different degrees of freedom (e.g., it would clearly be inappropriate to compare the percentage of the total sum of squares explained by a factor having only 1 degree of freedom versus some other factor that had 5 degrees of freedom). The only potentially unsettling consequence of using analogues of the ANOVA estimators to estimate components of variation is the fact that these estimates (even in the univariate case on the basis of Euclidean distance) can sometimes turn out to be negative ( Thompson (1962) , Searle, Casella & McCulloch (1992) ). This is clearly illogical and is generally accompanied by there being little or no evidence against the null hypothesis for the term in question that its component is equal to zero (i.e., a large P-value). Although there are other methods available for estimating variance components (i.e., ML, REML, Bayesian, etc., see Searle, Casella & McCulloch (1992) ), the ANOVA estimators do have the attractive quality of being unbiased34. The best solution to this issue is often to re-analyse the data after removing that term from the model (e.g., Thompson & Moore (1963) , Fletcher & Underwood (2002) ). This leads naturally to a consideration of how to remove terms from a model, also referred to in some cases as pooling (see the following section). 34 An unbiased estimator is one whose expectation is equal to the parameter it is trying to estimate. 1.29 Pooling or excluding terms For a given design file, PERMANOVA, by default, will do a partitioning according to all terms that are directly implied by the experimental design. For multi-factor designs, PERMANOVA will assume that all factors are crossed with one another, unless nesting is specified explicitly in column 2 of the design file. Although a factor cannot interact with a factor within which it is nested, factors that are crossed with one another necessarily generate interaction terms, and this full model (having all possible interactions) is generated by default. In some cases, however, one may wish to remove one or more terms from a given model. There are various reasons for wishing to remove individual terms, including:   (i) lack of evidence against the null hypothesis of that term’s component being equal to zero;   (ii) a negative estimate of that term’s component of variation;   (iii) previous studies have determined that term’s component to be zero or negligible;   (iv) hypotheses of interest require tests of models that exclude one or more particular terms. The user should be aware, however, that removing a term from a model equates with the assertion that its component of variation (that is, either S(*) or V(*) for a fixed or a random term, respectively, as the case may be) is equal to zero. By asserting that a component is equal to zero, one effectively combines, or pools, that term’s contribution (and its associated degrees of freedom) with some other term in the model. In the dialog of the PERMANOVA routine, we make a distinction between two different ways of removing a term: Excluding a term from the model – in which case the term is completely excluded and is not considered as ever having been part of the model in any form. Regardless of where the term occurs in the structure of the experimental design, excluding a term in this way is equivalent to pooling the df and SS for that term with the residual df and SS; Pooling a term – in which case the df and SS for that term is pooled with the term (or terms) which have equivalent EMS’s after that term’s component of variation is set to zero. For example, in a fully hierarchical design, this would correspond to a term being pooled with the term occurring immediately below it within the structure of the design. Complete exclusion of a term might be done, for example, in cases where we wish to construct a particular model that fully ignores those terms (e.g., in designs that lack replication, see the section on Split-plot designs). This is done by clicking on the ‘Terms…’ button in the PERMANOVA dialog. More generally, however, the removal of terms should be done with correct and appropriate pooling, where the component of variation for that term is set to zero and the EMS’s for the other terms in the model are re-evaluated. Pooling like this should be done, for example, to sequentially remove terms from a model having negative estimates for components of variation (e.g., Fletcher & Underwood (2002) ), or to remove terms having large P-values. In PERMANOVA, this is done by clicking on the ‘Pool…’ button in the PERMANOVA dialog. With respect to pooling on the basis of the reason given in (i) above, the assertion that a given term’s component is equal to zero should be made with some caution. Although a P-value > 0.05 (under the usual scientific convention) may not provide sufficient evidence to reject the null hypothesis (H$_0$), failing to reject H$_0$ is nevertheless logically very different from asserting that H$_0$ is true ( Popper (1959) , Popper (1963) )! There are differences of opinion regarding how large the P-value for a given term should be before the assertion of H$_0$ to remove that term might be justified (e.g., Hines (1996) , Janky (2000) ). “To pool or not to pool” is a decision left to the user, but we note that many practicing scientists use the rule-of-thumb suggested by Winer, Brown & Michels (1991) and Underwood (1997) that the P-value should exceed 0.25 before removing (i.e. pooling) any given term. Pooling a single term has important consequences for the construction of pseudo-F ratios, P-values and the estimation of components for the remaining terms. Thus, it is generally unwise to remove more than one term at a time (unless there are sufficient a priori reasons for doing so). The general rule suggested by Thompson & Moore (1963) and Fletcher & Underwood (2002) is, when faced with more than one term which might be removed from a given model, remove only one term at a time, beginning with the term having the smallest mean square, and at each step re-assess whether more terms should be removed or not. Fig. 1.30. Calculating the average taxonomic distinctness of molluscs for New Zealand holdfast assemblages. As an example of pooling, consider the New Zealand holdfast assemblages discussed in the previous two sections. Here, we shall focus on the analysis of a single variable – the average taxonomic distinctness (AvTD, $\Delta +$, Warwick & Clarke (1995) ) for molluscs in holdfasts. Open up the file hold.pwk and, from within this workspace, open up the aggregation file for the molluscs, called Mollusca.agg. Next, highlight the molluscs worksheet (see the section Nested design for details on obtaining this worksheet). Use the built-in tool in PRIMER to obtain the average taxonomic distinctness in each sample as a worksheet: select Analyse > Diverse and remove the $\checkmark$ in front of all of the default options except for ($\checkmark$AvTD: $\Delta +$) shown under the ‘Taxdisc’ tab & ($\checkmark$Results to worksheet) (Fig. 1.30). (Note that this also requires clicking on each of the ‘Other’, ‘Shannon’ and ‘Simpson’ tabs, in turn, in order to remove the $\checkmark$ for those options as well). Rename the resulting worksheet delta+ (which should contain only one variable, ‘Delta+’) and select Edit > Properties to give it the title: NZ holdfast AvTD for molluscs. Next, from the delta+ worksheet, calculate a Euclidean distance resemblance matrix: Analyse > Resemblance > (Analyse between •Samples) & (Measure •Euclidean distance). Do the analysis on the basis of the nested design (see the design file in Fig. 1.29) using the PERMANOVA routine with (Design worksheet: Nested design) & (Num. permutations: 9999) and all other choices as per the defaults. This analysis will result in an ANOVA partitioning (yielding values for df, SS, MS, F ratios and estimates of variance components) equivalent to that obtained using classical univariate ANOVA. (The data in this worksheet can be exported and analysed using a different statistical package to confirm this). The only difference will lie in the P-values, which of course are obtained using permutations in PERMANOVA, as opposed to using the traditional tables of the F distribution which rely on the assumption of normality. The results suggest that there is no significant variability in AvTD for molluscs among different sites (Fig. 1.31). Furthermore, the P-value for ‘Si(Lo)’ is quite large (P > 0.90) and the estimate of the variance component for ‘Si(Lo)’ is negative. Thus, using the rationale of either (i) or (ii) above, we may remove this term from the model by pooling it. Fig. 1.31. Results of PERMANOVA for the full 3-factor nested model on AvTD of molluscs. This is achieved relatively easily, by re-running the PERMANOVA routine on the basis of the same design file, but this time, click on the ‘Pool…’ button (Fig. 1.32). A new dialog box appears entitled ‘Selection’ with a list of all of the terms in the model. The user may choose which terms in the model to pool. For the present case, we wish to pool the ‘Si(Lo)’ term. In the ‘Available’ box, click on the term Si(Lo) followed by to move this into the ‘Include’ box, then ‘OK’. To understand how pooling is done in the PERMANOVA routine, consider the details of the EMS in the present example for the full model before pooling: By pooling the term ‘Si(Lo)’, we are explicitly asserting that the component V(Si(Lo)) = 0. Setting this component deliberately to zero wherever it appears yields the following: Thus, the mean square for the term ‘Si(Lo)’ and the mean square for the term ‘Ar(Si(Lo))’ are now estimating the same thing, i.e. they have the same expectation. This means that their SS and df can be added together to obtain a pooled MS, as follows: $$ MS _ {pooled} = \frac {\left( SS _ {Si(Lo)} + SS _ {Ar(Si(Lo))} \right) } { \left( df _ {Si(Lo)} + df _ {Ar(Si(Lo))} \right) } \tag{1.5} $$ The PERMANOVA output file from the analysis after pooling identifies which terms were pooled, under the heading ‘Pooled terms’ and also identifies the terms whose SS and df were combined as a consequence of this (Fig. 1.32). The new pooled term is given a unique name – it is simply called ‘Pooled’ in the present case. In the event of there being more terms to pool, potentially more than one pool of terms may occur in a single analysis. This is also catered for by the PERMANOVA routine, if required, with each pool identified by its component terms and given its own unique name. The EMS’s after pooling in the present case are: and the associated degrees of freedom for the pooled term are: $$ df _ {Pooled} = \left( df _ {Si(Lo)} + df _ {Ar(Si(Lo))} \right) = \left( 4 + 8 \right) = 12 \tag{1.6} $$ Pooling the ‘Si(Lo)’ term has resulted in all of the remaining estimated components of variation in the model being non-negative. It has also, however, substantially changed the F-ratios and tests for the remaining terms. The term ‘Lo’ is now not statistically significant (pseudo-F = 1.16, P = 0.36, Fig. 1.32), whereas before, when tested using the mean square for ‘Si(Lo)’ as the denominator, it was (pseudo-F = 19.1, P = 0.01, Fig. 1.31). This is an important point. Although pooling may result in an increase in power, caused by an increase in the denominator df for the test of a given term (here, ‘Den.df’ for the test of ‘Lo’ has gone up from 4 to 12), this is also often off-set by an increase in the denominator MS for the test after pooling (cf. $MS _ {Si(Lo)} = 0.68$, whereas $MS _{Pooled} = 11.2$), which reduces the value of pseudo-F for the test. It is usually not possible to tell a priori just how the tests for other terms in the model will be affected by pooling. Clearly, however, estimated components of variation, pseudo-F ratios and P-values will be affected by pooling. Thus, as stated previously, pooling should only be done one term at a time. From the present analysis, after pooling, there is apparently no statistically significant variability in the AvTD of molluscs among holdfasts at any of the spatial scales examined (Locations, Sites or Areas). Note that these results (Fig. 1.32), obtained by removing the ‘Si(Lo)’ term using the ‘Pool…’ button in the PERMANOVA dialog, are not the same as the results that would have been obtained if we had simply excluded ‘Si(Lo)’ from the model entirely, using the ‘Terms…’ button. In that case, the term would have been considered to be non-existent. It would not have been included in the partitioning at all and its SS and df would therefore have ended up as part of the residual. This is clearly illogical and undesirable in the present case. Removing the ‘Si(Lo)’ term should result in an increase to the df associated with the denominator MS being used to test the ‘Lo’ term (i.e., its SS and df should be combined with the ‘Ar(Si(Lo))’ term), and not merely ignored and added to the residual variation. Another possibility would be to re-cast the model with two factors: ‘Lo’ and ‘Ar(Lo)’, but we would need to be careful and make sure that the areas within any location (even those from different sites) each had a unique name in the specification of the ‘Areas’ factor levels. Fig. 1.32. Dialog to pool the ‘Si(Lo)’ term and results for the PERMANOVA analysis of the AvTD of molluscs after pooling. If pooling of a given term would result in its being combined with the residual in any event (which of course, occurs legitimately in some cases), then the results using these two approaches will be equivalent. Otherwise, however, they will not. We consider that removal of terms using the ‘Pool…’ button will be most appropriate for the majority of situations. The use of the ‘Terms…’ button to exclude terms entirely may be useful, however, to craft specific simplified models in the event that replication is lacking (see the section Designs that lack replication). The dialog available under the ‘Terms…’ button can also be used to change the order in which individual terms are fitted. Although the order in which terms are fitted is of no consequence for balanced ANOVA designs, it is important for analyses of unbalanced designs or designs including covariates using Type I (sequential) SS (see the sections Unbalanced designs and Designs with covariates). 1.30 Designs that lack replication (Plankton net study) A topic related to the issue of pooling is the issue of designs that lack replication. Familiar examples are some of the classical experimental designs, primarily from the agricultural literature, such as randomised blocks, split plots or latin squares (e.g. Mead (1988) ). For these designs the essential issue is that there is no replication of samples within cells, but rather there is only 1 sample per cell. This means that it is not possible to distinguish between variation among samples and variation due to the highest-order (most complex) interaction term. Thus, in order to proceed, the experimenter has either to assume (i) that the highest-order interaction term is zero, or (ii) that the so-called “residual” mean square in the model actually has expectation V(Res) + V(highest-order interaction). Note that, for the latter assumption to work, at least one of the factors involved in the highest-order interaction has to be random. See Gates (1995) and chapter 10 of Quinn & Keough (2002) for further discussion of these issues. From a practical perspective, for PERMANOVA to proceed with the analysis (regardless of which of the above two perspectives one chooses to take), the highest-order interaction term needs to be excluded from the analysis (see the section Pooling or excluded terms). This can either be done manually, or if the PERMANOVA routine detects that there is no within-cell replication, then it will issue a warning. If you choose to proceed by clicking ‘OK’, it will automatically exclude the highest-order interaction term from the model. If you receive this warning and you know that you do have within-cell replication, then there is a very good chance that you have mis-labeled your factor levels somehow35. An example of a two-way crossed design without replication is provided in a study by Winsor & Clarke (1940) to investigate the catch of various groups of plankton by two nets hauled horizontally, with one net being 2 metres below the other. Ten hauls were made with the pair of nets at depths of 29 and 31 meters, respectively. The experimental design is: Factor A: Position (fixed with a = 2 levels, either upper (U) or lower (L) depths). Factor B: Haul (random with b = 10 levels, labeled simply 1-10). There is only 1 value per combination of treatments, with no replication, so N = a × b = 20. This is effectively a randomised block design, where the hauls are “blocks”. The variables recorded correspond to five different groups of plankton. Standard deviations in the various groups were roughly proportional to the means, so data were transformed and are provided as logarithms of the catch numbers for each of the plankton groups. These data are located in the plank.pri file in the ‘Plankton’ folder of the ‘Examples add-on’ directory, and were provided by Snedecor (1946) . Fig. 1.33. PCA of the study of plankton from ten hauls (numbered) at either 29 m depth (upper) or 31 m depth (lower). Examination of the data (already log-transformed) reveals no zeros and that their distributions are fairly even, with no extreme values or outliers36. The variables are also on similar scales and are measured in the same units; therefore, an analysis based directly on Euclidean distances would be reasonable here – prior normalisation is not necessary. For data like these, an appropriate ordination method is principal components analysis (PCA)37. The first two principal components explained 83.6% of the total variance in the five variables (Fig. 1.33). Variability among the hauls is apparent in the diagram, but a clear difference in the plankton numbers due to the position of the nets (upper versus lower), if any, is not obvious. The PERMANOVA analysis of these data on the basis of a Euclidean distance matrix has detected significant variability among the hauls, but also has detected a significant effect of the position of the net (Fig. 1.34). Notice that the output has identified the excluded term: ‘PositionxHaul’, as an important reminder that the analysis without replication is not without an additional necessary assumption in this regard. It might seem surprising that the analysis has detected any effect of ‘Position’ at all, given the pattern seen in the PCA (Fig. 1.33). Looks can be deceiving, however. Close inspection of the plot reveals that, within almost all of the individual hauls (i.e., 1, 2, 3, 5, 7, 9, 10), the symbol for the ‘upper’ group lies to the right of the symbol for the ‘lower’ group. Only hauls 4, 6 and 8 do not conform to this pattern. We can perhaps understand the nature of this overall effect of Position by examining the averages for the ‘upper’ and ‘lower’ nets for each of the variables across all of the hauls. With the Plankton worksheet highlighted, select Tools > Average > (Samples •Averages for factor: Position) & (Variables •No averaging). The resulting worksheet shows that the average log(abundance) for all five of the plankton variables was larger for the nets towed at the shallower depth (the ‘upper’ nets) (Fig. 1.34). Fig. 1.34. Design file and PERMANOVA analysis for the two-way study that lacks replication within cells. The interaction term is, by necessity, excluded from the analysis, as it is already confounded with the residual variance. Also shown are the averages per depth for each of the 5 variables in the plankton study. Another important point here is to recognise that, had we treated the above design as if the hauls were the replicates, and ignored the variation among hauls, then we would not have detected any effect of position at all. You can check this fact by running PERMANOVA on the data using a one-way design with the factor ‘Position’ only (the result is non-significant, with P > 0.25). Thus, despite the fact that there is a consistent shift in the plankton assemblage between the upper and lower nets within each haul, the variation from haul to haul would have masked this entirely and we would have failed to detect it (as we at first did when contemplating the PCA plot), if we had not included the factor ‘Haul’ in our analysis. The advantages of “blocking” to achieve greater power to detect treatment effects have been known for a very long time (e.g., Fisher (1935) , Snedecor (1946) , Mead (1988) ), but this example shows that the phenomenon can occur equally strikingly in the analysis of multivariate data. The analysis of some other designs that lack replication have other issues, on top of the one already noted regarding the highest-order interaction being inextricably confounded with the residual. For example, the experimental design known as the latin square consists of a random allocation of t treatments to a t × t matrix of sample units, with the added constraint that there be one of every treatment in every row and one of every treatment in every column of this array. The usual model fitted to such a design partitions the sum of squares according to the following sources: rows (R), columns (C) and treatments (T). None of the potential interaction terms (R×C, R×T, C×T, R×C×T) are traditionally included in these models, because none of them can be readily unconfounded from the residual. PERMANOVA does not have separate subroutines for treating these kinds of special designs, but will not give sensible results38 unless the terms which cannot be estimated are first removed from the model. For these more complex designs lacking replication, it is up to the user to know if such interactions need to be excluded, to understand the consequences of the assumptions underlying these models if they are to be used, and to exclude the relevant terms manually, using the ‘Terms…’ button in the PERMANOVA dialog. 35 For example, you may have given the levels of factor B the names b1, b2, b3 within the first level of factor A, but then called them B1, B2, B3 within the second level of factor A, and PRIMER will not interpret these names as being the same. 36 PRIMER’s Analyse>Draftsman plot routine is very useful for visually examining the distributions and joint distributions of variables in a worksheet. 37 For more details regarding this method and its implementation in PRIMER, see chapter 4 of Clarke & Warwick (2001) and chapter 10 in Clarke & Gorley (2006) . 38 (or, at least, not the results that are traditionally provided for such designs). 1.31 Split-plot designs (Woodstock plants) Another special case of a design lacking appropriate replication is known as a split-plot design. These designs usually arise in an agricultural context, where the experimenter has applied the treatment levels for a factor (say, factor A) randomly to whole plots (usually within blocks of some kind) at a large scale, but then, within each of these whole plots, treatments for another factor (say, factor B) are applied randomly to smaller units. Thus, the whole plots are each “split” into smaller units. Although there are variations on this theme, the traditional split-plot design lacks replicates of the whole plots at the larger spatial scale (i.e., it is a randomised block design for factor A, with whole plots acting as the ‘error’). There is also commonly a lack of replication at the smaller spatial scale (i.e., the number of levels of factor B is equal to the number of sample units within each whole plot and these levels are allocated randomly and separately within each whole plot). A proposed rationale for using a split-plot design is that factors may occur naturally at different scales (e.g., Mead (1988) ). Another proposed rationale is that one may already know that factor A has important effects, and one may be willing to sacrifice information on factor A to get more precise results for factor B and the interaction A×B. Although neither of these actually provides a solid argument for ignoring the need for appropriate replication in experiments, split-plot designs do still occur from time to time in biological research and can be analysed (after some thought and with care) using PERMANOVA. For more information regarding the assumptions and potential disadvantages of split-plot designs, see Mead (1988) and Underwood (1997) . Fig. 1.35. Schematic diagram of the Woodstock split-plot design examining the potential effects of fire frequency (no burning, burning every 2 years, every 4 years or every 8 years) and the effect of grazers (by fencing some sub-plots (F), while leaving others unfenced (U)) on plant assemblages. To analyse a split-plot design in PERMANOVA, effectively two analyses must be done: one at the ‘whole-plot’ level and one at the ‘sub-plot’ level. Results of these two analyses can then be combined to form the traditional partitioning that is usually presented for such designs. An example of a split-plot design is provided by a study of the effect of fire disturbance and fencing (to exclude grazers) on the composition of plant assemblages on the central western slopes of New South Wales in south-eastern Australia ( Prober, Thiele & Hunt (2007) ). The experimental design (shown schematically in Fig. 1.35) included the following factors:   Blocks: (random with r = 4 levels).   Factor A: Fire frequency (fixed with a = 4 levels: 0 yrs, 2 yrs, 4 yrs or 8 yrs).   Whole-plots: (random and nested within Blocks and Fire frequency, unreplicated).   Factor B: Fencing (fixed with b = 2 levels: fenced and unfenced).   Sub-plots: (random and nested within all of the above, unreplicated). The two fencing treatments were randomly allocated to two sub-plots (measuring 5 m × 5 m) within each fire treatment (whole plots) and there is one of each fire treatment (4 whole plots) randomised within each block. Relative abundances (cover) of higher plant species within each sub-plot were estimated using a point-intercept technique (an 8 mm dowel placed vertically at each of 50 points on a grid across each plot). Although the design was set up at each of two locations (Woodstock and Monteagle) and data were obtained over a number of years ( Prober, Thiele & Hunt (2007) ), we consider here only data from the Woodstock location collected in 2003. We also will exclude two species: Poa sieberiana and Themeda australis, the dominant grasses, from the analysis. These have already been analysed separately in detail ( Prober, Thiele & Hunt (2007) ) and our focus here instead will be on the more subtle potential responses of subsidiary forbs and exotic species.               Table 1.2. Sources of variation and degrees of freedom for a traditional partitioning according to the split-plot design for the Woodstock experiment. Partitioning for a split-plot design is traditionally done according to Table 1.2. The upper-level factors (Blocks and Fire frequency) are tested against the whole-plot error, while the lower-level factors (Fencing and the Fire × Fencing interaction term) are tested against the sub-plot error. To do this partitioning and the necessary tests using PERMANOVA, we first focus on the top-half of Table 1.2 only. This calls for a randomised block design, but where the whole plots are effectively treated as the sample units. So, first we need to obtain distances among centroids39 for the whole plots. Open the file wsk.pri (in the ‘Woodstock’ folder of the ‘Examples add-on’ directory) containing the abundances (cover measures) for p = 117 plant species. First select all of the variables except the variables numbered 50 and 63 in the dataset (‘Poa sieb’ and ‘Themaus’, respectively). Calculate a Bray-Curtis resemblance matrix after square-root transforming the selected data and choose PERMANOVA+ > Distances among centroids… > Grouping factor: WholePlot, then click ‘OK’. This yields a new matrix of Bray-Curtis resemblances among the 16 whole plots (4 fire treatments × 4 blocks). Next, run a two-way randomised block design of Block and Fire on the resemblance matrix among centroids (‘Resem2’). As there is no replication of the whole plots, you will either have to remove the ‘Block × Fire’ interaction term manually (by clicking on the ‘Terms…’ button in the PERMANOVA dialog), or let PERMANOVA do that for you. This analysis will give results for the top half (the first four terms) of the table (Fig. 1.36). Fig. 1.36. Step one in the PERMANOVA analysis of the Woodstock split-plot design: a two-way randomised block design for whole-plots. Now, for the lower half of the table, we will want to fit the ‘Fence’ and ‘Fire × Fence’ terms, given the whole-plots, in order to get the correct sub-plot error. Go back to the original resemblance matrix among all 32 samples (‘Resem1’) and set up a PERMANOVA design file with three factors: WholePlot, Fire frequency and Fencing (Fig. 1.37). As usual, PERMANOVA attempts to construct and fit a full model, including all interaction terms implied by the structure that is specified. In the present case, it is not possible to fit all interaction terms, because of the lack of replication. Here, it is not just the lowest level in the analysis that is unreplicated, but we also lack replication at a higher level in the design. PERMANOVA will not automatically exclude the terms that would normally be excluded from a split-plot analysis, so these must be excluded manually by clicking on the ‘Terms…’ button in the PERMANOVA dialog and choosing only those terms in the model that we wish to fit (Fig. 1.37). For this design, there are a number of terms that are purposefully excluded: namely, any interactions involving whole plots as well as the main effect of Fire. (This is on top of the fact that we have also chosen to ignore any possible interactions involving Blocks.) It is important to recognise these assumptions underlying any split-plot analysis, and the requirement to manually remove terms is a good reminder of what is going on here. Of course, if we had replicated whole-plots and replicated sub-plots, we would be in a position to analyse the whole design, including all interactions, in a single PERMANOVA analysis. Fig. 1.37. Step two in the PERMANOVA analysis of the Woodstock split-plot design: a three-way analysis of sub-plots, excluding certain terms, but including whole-plots. Once we have the results from both “halves” of the split-plot analysis, we can use these to construct the complete table (Table 1.3). When combining the results of more than one analysis from a single set of data, such as this, it is a good idea to check that the partitioning has been done correctly by making sure that the sum of the individual SS add up to the total SS. This is true for the present example (Table 1.3) and will be true in general, at least for a correct partitioning of any balanced design. The analysis suggests that the fencing treatments had no significant effect on these assemblages, but fire frequency did. There was no evidence that these two factors interact significantly with one another. Spatial variation among blocks and among whole plots was also substantial (see the relative sizes of components of variation in the output file).      Table 1.3. Results from the PERMANOVA analysis of plant assemblages in response to fire frequency and fencing (removal of grazers) in the Woodstock split-plot experiment. 39 Importantly, these centroids are not calculated on the original data, they are calculated on the full set of PCO axes obtained from the resemblance matrix, in order to preserve the resemblance measure chosen as the basis of the analysis. For more details, see chapter 3. 1.32 Repeated measures (Victorian avifauna, revisited) While randomised blocks, latin squares and split-plot designs lack spatial replication, a special case of a design lacking temporal replication (and which occurs quite a lot in ecological sampling) is the repeated measures design (e.g., Gurevitch & Chester (1986) , Green (1993) ). In essence, these designs consist of individual sampling units (usually belonging to various treatments, etc.) which are repeatedly examined at several different time points. Such designs receive special attention for two reasons: (i) they do not have replication within cells, so suffer from the same issues discussed in the previous two sections and (ii) they require some additional assumptions because samples are generally considered to be non-independent through time. The first issue – lack of replication – is dealt with easily enough by PERMANOVA, as discussed above. The program (after issuing the usual warning) will simply exclude the highest-order interaction term (i.e. which will include “Time” as one of its members) and the analysis proceeds from there in the usual way. The issue of non-independence is, however, another matter. In a traditional repeated measures analysis of univariate data, the partitioning of the total sum of squares is done in the usual way, treating “Time” as a fixed factor. What is taken on as an additional assumption, however, is something known as sphericity. Sphericity is an assumption about the nature of the correlations through time for the sample units, which must be similar for the different treatments. Although a formal test of sphericity is provided by Mauchley (1940) (see Winer, Brown & Michels (1991) for an example), this approach is unfortunately rather highly susceptible to deviations from normality ( Huyhn & Mandeville (1979) ). Huynh & Feldt (1970) have demonstrated that a necessary and sufficient condition is to check the equality of the variances of the differences between levels of the repeated measures factor (“Time”) across treatments (or combinations of treatments). Thus, one calculates the differences between each pair of time points, obtains the variances of these difference values for each treatment (or treatment combination) and then checks for equality of these variances (e.g., Quinn & Keough (2002) ). If this assumption is violated, then the F ratios obtained from the analysis are no longer distributed like traditional F distributions under true null hypotheses. In this case, for traditional univariate analysis, a number of possible corrections to the degrees of freedom can be done to get a correct test (e.g., Box (1954) , Geisser & Greenhouse (1958) , Huynh & Feldt (1976) ). When dealing with multivariate response data, one might consider doing an analogous test of sphericity (of some sort) by calculating the dissimilarities (or distances) between levels of the repeated measures factor across treatments. The variances of these dissimilarities could then be compared among treatments (or treatment combinations) using, for example, PERMDISP (see chapter 2), or even using a traditional test for homogeneity of variances among groups. Recall, however, that PERMANOVA uses permutation procedures in order to generate a correct distribution of each pseudo-F statistic under a (relevant) true null hypothesis. So the only essential assumption associated with the use of PERMANOVA (whether there be repeated measures or otherwise) is the exchangeability of samples (or of appropriate residuals). It must be admitted that the correlation structure among samples through time, if any, will be effectively ignored under permutation. Thus, differences in correlation structure through time among treatments (i.e. lack of sphericity) may, therefore, produce a statistically significant result. However, we consider that differences in correlation structure through time are indicative of (at least one type of) a treatment effect, so should warrant closer inspection by the investigator in any event. Clearly, the degree to which correlation structure (in space or in time, as in repeated measures) can affect the results of permutation tests for repeated measures designs (or any other designs for that matter) warrants further study. If statistically significant results are obtained in a repeated measures analysis, the user may wish to accompany the PERMANOVA with a separate test for sphericity (in the case of univariate data) or its analogue (using dissimilarities rather than differences) for multivariate responses, in order to shed further light on the meaning of the results and appropriate inferences. An example helps to clarify these ideas. The data on Victorian avifauna, previously examined in the section Monte-Carlo P-values, actually consisted of a repeated measures design, as described by Mac Nally & Timewell (2005) . The previous analyses (Figs. 1.12 and 1.13) were based on data summed across four different observation times. However, these data are also available at the level of individual surveys done at each time in the file vicsurv.pri in the folder ‘VictAvi’ of the ‘Examples add-on’ directory. Open up the data and calculate the dissimilarity matrix among samples on the basis of the binomial deviance dissimilarity measure; re-name this matrix BinomDev. Fig. 1.38. MDS of Victorian avifauna at each of two sites within each of several states of flowering intensity (poor, medium, good or adjacent to a good site) for 4 times (given as numbers on the plot). An MDS ordination on the basis of this matrix suggests a strong effect of flowering intensity on the bird communities, with ‘good’ sites appearing to the left of the diagram, ‘poor’ or ‘medium’ sites occurring to the right, and ‘adjacent’ samples being more variable than the other treatments through time (Fig. 1.38). The repeated measures experimental design here is:     Factor A:   Treatment (fixed with a = 4 levels: poor, medium, good or adjacent).     Factor B:   Site (random, nested in Treatments with b = 2 levels, labeled simply S1-S8).     Factor C:   Time (fixed40 with c = 4 levels, labeled 1-4). After creating the design file, re-name it Repeated measures for reference. Note that there is no replication within cells, as there is only one sample per site at each time. We therefore need to exclude the highest-order interaction, ‘Site(Treatment) × Time’, from the model. See Fig. 1.39 and refer to the above section Pooling or excluding terms, if necessary, for details on how to do this. Analysis by PERMANOVA (with 9999 permutations of residuals under a reduced model) reveals significant treatment effects on the bird assemblages (Fig. 1.39), and shows that these effects are consistent through time (note that P > 0.6 for ‘TrxTi’). Pair-wise comparisons (not shown here, but you can do them easily) reveal that this effect is largely due to there being significant differences between the good sites vs the others. Having identified significant treatment effects, we may wish to examine the multivariate analogue to the test of sphericity for univariate data. Namely, we can calculate the dissimilarities between each pair of time points within each treatment. Then, we can compare the time-point differences for equality of variances. This is purely optional and is not a requirement of the repeated measures analysis when performed using PERMANOVA, but it may nevertheless shed some further light on whether the significant treatment effects detected were due only to differences in location in multivariate space (as would appear to be the case from the diagram) or whether sites might also vary in the nature of their non-independence among time points. To do this, we first need to extract the relevant dissimilarities between time-points for each sampling unit (in this case, from each site) from the full dissimilarity matrix. These are shown in Table 1.4. For example, the binomial deviance dissimilarity between times 1 and 2 for site 1 (a ‘poor’ site) is 17.85, and so on. These multivariate dissimilarity values can take on the same role as the differences among time points which are usually examined in a univariate repeated measures analysis (e.g., see Box 10.6 in Quinn & Keough (2002) )41. The estimated variances in these dissimilarities (given in the last row of Table 1.4) appear to be similar among the six paired groups (i.e. the six columns in Table 1.4). Levene’s test (using either means $F _{5,42} = 0.59$, P > 0.71 or medians $F _{5,42} = 0.32$, P > 0.90), supports the assumption of homogeneity of these variances42. This means we can fairly safely infer in this case that the treatment effects we detected were not caused by differences in dissimilarity structure within samples through time among the sites43. Fig. 1.39. PERMANOVA analysis of a repeated measures experimental design for Victorian avifauna. Before leaving the topic of repeated measures, there is one other method of analysis worth mentioning. With univariate data, one has the option of treating the individual values at each time point as separate variables, then doing a multivariate analysis among treatments. Note that this approach will only work when there is one response variable measured at different times. It is to be distinguished from the situation we have above, where there are already multivariate data (i.e., abundances of 27 species of birds) obtained at each time point. The multivariate approach to a repeated measures analysis of a single response variable is, however, straightforward to do using PERMANOVA – one simply organises the data so that measures at different time points are the “variables” and uses Euclidean distance as the basis of the analysis. Although taking a multivariate approach with univariate data (treating the time points as variables) completely avoids having to consider any notions of sphericity, it does not provide any test of the factor ‘Time’ or any tests of the interactions of ‘Time’ with (most) other factors, which would, of course, be provided by a traditional repeated measures partitioning. Green (1993) discusses other pros and cons with using the multivariate versus the univariate approach to a repeated measures analysis of a single variable (see Fig. 6 therein). Table 1.4. Binomial deviance dissimilarities between time points for each site and their estimated variances. For repeated measures designs and in other cases where there is known to be correlation structure (non-independence) among replicates, there may be other ways to analyse the data. First, one might consider rephrasing hypotheses in terms of dissimilarities between particular pairs of correlated objects, then analysing those dissimilarities in a univariate analysis. For example, Faith, Humphrey & Dostine (1991) tested the null hypothesis of no difference in the average dissimilarity between an impact and a control site from before to after the onset of a disturbance (with individual time points treated as replicates). Similarly, if one has a baseline or control against which other treatments are to be compared through time, one might consider using principal response curves (PRC, van den Brink & ter Braak (1999) ), which is simply a special form of redundancy analysis (RDA) . Distance-based redundancy analysis (dbRDA) can also be used to model changes in the community through time (or space) explicitly as a linear, quadratic or other polynomial traveling through the multivariate data cloud (e.g. Makarenkov & Legendre (2002) ), rather than treating time (or space) as an ANOVA factor (see the chapters on DISTLM and dbRDA below for more details). Larger-scale monitoring programs, which generally have many sites sampled repeatedly at many times, can also be analysed using multivariate control charts ( Anderson & Thompson (2004) ). This approach is designed to detect when (and where) individual sites deviate significantly from what would be expected, given natural temporal variability. 40 We have chosen to treat “Time” as a fixed factor here, because this is traditionally how it is treated in repeated measures experimental designs. However, the user can of course choose to treat this factor as random if this is more in line with relevant hypotheses of interest. 41 Bear in mind that the univariate differences, however, can show direction by being either positive or negative, which will affect calculated variances. In contrast, dissimilarities are always positive, so cannot show direction, per se. Therefore, this is not at all intended to be a strict test for sphericity, merely to provide a means of examining the null hypothesis of no difference in dissimilarity structure among time points for individual samples across the different treatments. See also the approach used by Clarke, Somerfield, Airoldi et al. (2006) . 42 Note the use of a traditional univariate Levene’s test here for testing homogeneity of variances. We could also have done this test using a permutational approach in the routine PERMDISP on the basis of a Euclidean distance matrix for the univariate variable produced in Table 1.1. Indeed, if the tables are used, then this is equivalent to doing a traditional Levene’s test. See the chapter on PERMDISP for more details. 43 Winer, Brown & Michels (1991) discuss how correlations may need to be examined at several different levels in more complex multi-factorial repeated measures designs (see chapter 7 therein). We do not pursue this further here. 1.33 Unbalanced designs Virtually all of the examples thus far have involved the analysis of what are known as balanced experimental designs. For these situations, there is equal replication within each level of a factor (or within each cell). Even data that lack replication have an equal number of n = 1 within each cell. However, despite a researcher’s best efforts, sometimes the design ends up being unbalanced, where there are unequal numbers of replicate samples within each factor level (or within each cell) of the design. For the one-way case (such as the Ekofisk oil-field data seen in the section One-way example), the consequences of an unbalanced design are not problematic. We can perform PERMANOVA in the usual way, with the usual partitioning of sums of squared dissimilarities. Only two consequences of unequal replication are apparent for one-way designs. First, the multiplier on the EMS for the factor of interest is no longer necessarily a whole number, as it was for the balanced case. By scrolling down the PERMANOVA results window produced in the analysis of the Ekofisk oil-field data (shown only partially in Fig. 1.10 above), we can see the multiplier for the ‘S(Di)’ component of variation in the EMS for Distance is 9.57, whereas for any of the balanced designs, the multipliers for any component in any EMS are whole numbers (e.g., see the EMS in Figs. 1.15, 1.23 or 1.29, all of which are balanced designs). For one-way designs, this does not change, however, the use of the residual MS in the denominator of the pseudo-F ratio for the test of “No differences among the groups”. The second consequence of an unbalanced design is apparent when we consider the permutations. We still randomly allocate observation units across the groups (levels of the factor), while maintaining the existing group differences in sample sizes. Each individual unit no longer has an equal chance of falling into any particular group, but instead will have a greater chance of falling into a group that has a larger sample size. However, we can still proceed easily on the basis that all possible re-arrangements of the samples by reference to the existing (albeit unbalanced) experimental design are equally likely. The more important issues facing experimenters with unbalanced designs occur when there is more than one factor in the design. In that case, the consequences are: (i) the multipliers on individual components of variation in the EMS’s are not necessarily whole numbers and these multipliers can differ for the same component when it appears in the EMS’s of different terms in the model; (ii) the main effects of factors and the interaction terms are no longer independent of one another. The latter is perhaps the most important conceptual issue in the analysis of unbalanced, as opposed to balanced, designs. This means that, like in multiple regression (see chapter 4), the order in which we choose to fit the terms matters. Fig. 1.40. Venn diagrams showing the difference between (a) a balanced and (b) an unbalanced two-way crossed design. Begin by considering a two-way crossed ANOVA design, with factors A, B and their interaction A×B. If the design is balanced, then the individual amounts of variation in the response data cloud explained by each of the terms in the model are completely independent of one another. This can be visualised using a Venn diagram (Fig. 1.40), where the total variation in the system ($SS _T$) is represented by a large circle and the residual variation, $SS _ {Res}$, is the area left over after removing all of the portions explained by the model. For a balanced design, the individual terms in the model explain separate independent portions of the total variation (Fig. 1.40a), whereas for an unbalanced design, there will be some overlap among the terms regarding the individual portions of variation that they explain (Fig. 1.40b). 1.34 Types of sums of squares (Birds from Borneo) When the design is unbalanced, there will be a number of different ways to do the partitioning, which will depend to some extent on our hypotheses and how we wish to treat the potential overlap among the terms. The different ways of doing the partitioning are called “Types” of sums of squares. More particularly, there are (at least) four types, known (perhaps unhelpfully) as Type I, II, III and IV. This terminology was initially coined by the developers of the SAS computer program (e.g., SAS Institute (1999) ), and is now in common usage. All of these types of SS produce identical results for balanced designs. Furthermore, Types II, III and IV will be identical for models with no interactions and Types III and IV will be identical if all cells are filled (i.e. if all cells have n ≥ 1). Searle (1987) provides an excellent text regarding models and hypotheses for unbalanced designs, including a comparison of the types of SS, which are briefly described below: Fig. 1.41. Schematic Venn diagrams demonstrating the conceptual differences in Types of SS for a two-way crossed unbalanced design. Type I SS. This approach may be described as sequential. Here, each term is fitted after taking into account (i.e., conditioning upon or treating as covariates) all previous terms in the model. The order of the terms listed in the design file for the analysis therefore matters. For example, in the two-way crossed design, if the order of the terms listed in the analysis were A, B and A×B, then the SS would calculate the SS for A (ignoring other terms), the SS for B (given that A is already in the model) and then the SS for A×B (given that A and B are both already in the model). This is shown diagrammatically in Fig. 1.41a. You would get a different partitioning for this same design, however, if you chose to fit factor B first and then A given B (Fig. 1.41b). Another thing to note about Type I SS is that the sum of the individual SS will add up to the total SS. The sequential analysis offered by the Type I approach may be appropriate for fully nested hierarchical models, for which there exists a natural ordering of the terms. In other cases, Type I SS may be used to explore the amount of overlap in the explained variability among terms, by changing the order of the terms in the analysis and seeing how this affects the results. Type II SS. This approach might be described as a conditional analysis. The Type II SS for a given term is defined as the reduction in the residual SS due to adding the term after all other terms have been included in the model except any terms that contain the effect being tested. In other words, main effects are not conditioned upon interaction terms that involve them. In the two-way crossed design, this amounts to fitting A given B, B given A and A×B given both A and B (Fig. 1.41c). Note that, given the explicit definition, the order in which the individual terms are fitted will not matter here. However, the individual SS in the analysis will not necessarily add up to the total SS. There will potentially be some “bits missing” as a consequence of this partitioning. (In the two-way case, the “bit missing” is the overlap of factors A and B in their explained variation). Type III SS. This approach might be described as a fully partial analysis. Every term in the model is fitted only after taking into account all other terms in the full model. Thus, for our example, we can see this amounts to fitting A given both B and A×B, B given both A and A×B, and A×B given both A and B (Fig. 1.41c). Like for Type II SS, the nature of the definition here ensures that the order in which terms are fit will not matter. However, also like Type II SS, the sum of the individual SS will not add up to the total SS, and there will be more “bits missing” (ignored regions of overlap) for the Type III case. However, complete orthogonality (independence) of all of the hypotheses is ensured using this method. Type IV SS. This approach was developed by the SAS Institute (e.g., SAS Institute (1999) ) to deal with cases of a particular kind of imbalance known as “some cells empty”, in which data from some of the cells (combinations of treatments) are actually missing entirely (i.e. n = 0 for those cells). These situations can be contrasted with situations known as “all cells filled”, which may have either unequal or equal (balanced designs) replication per cell. According to Searle (1987) , Type IV SS are for testing hypotheses determined rather arbitrarily by the SAS GLM routine, which depends not only on which cells have data in them, but also on the order in which levels of the factors happen to have been listed. As such, Type IV SS is not generally recommended and is not discussed further here. The situation of “some cells empty” is quite adequately dealt with using Type III SS. The PERMANOVA dialog box offers the user the option of using Type I, II or III SS. The default in PERMANOVA is to use Type III SS. This is primarily because most editors of journals (at least, most ecological journals) have now come to expect Type III SS to be used for unbalanced designs, simply because these will tend to be the most conservative of the three. However, there is no particular reason not to use the other types, especially if one of these is better suited to particular hypotheses of interest. For example, as already noted, a sequential analysis (Type I SS) would be quite sensible to use for an hierarchical nested design. Indeed, many statisticians would consider Type I SS to be the most sensible general approach, as no components of variation are left out (i.e. there are no “bits missing”). Type I SS also allows clarification of the relative sizes of overlapping regions, when terms are fitted in different orders. In most cases, provided the degree of imbalance in the design is modest (due, for example, to just a few missing observations here and there), the overall conclusions of the study will be little affected by this choice. Fig. 1.42. Sample sizes per cell in the unbalanced two-way layout for the Borneo birds example. A case in point is provided by an analysis of bird assemblages from Borneo, Indonesia in response to a two-way unbalanced design, as described by Cleary, Genner, Boyle et al. (2005) . A total of N = 37 sites were sampled within the Kayu Mas logging concession, close to Sangai, Central Kalimantan. The sites were cross-classified according to two factors: Logging (fixed with a = 3 levels: unlogged primary forest, forest logged in 1993/94 and forest logged in 1989/90) and Slope (fixed with b = 3 levels: lower, middle and upper). There were different numbers of sites (n) within each of the a × b = 3 × 3 = 9 cells in the design, as shown in Fig. 1.42. Within each site, spot-mapping (using calls and visual observations) was used to sample birds along each of two parallel 300 m linear transects, 50 m apart at each site. There were p = 177 bird species recorded in all and the data are located in the file born.pri, located in the ‘BorneoBirds’ folder of the ‘Examples add-on’ directory. Fig. 1.43. MDS ordination of bird assemblages from Borneo in all combinations of logging (P = primary forest, L89 = logged in 1989/90, L93 = logged in 1993/94) and slope (L = lower, M = middle and U = upper). An MDS plot of the bird communities, on the basis of log(x+1)-transformed data and Bray-Curtis similarities ( Cleary, Genner, Boyle et al. (2005) ) shows a clear effect of logging, and suggests some effects of slope as well, although these are less clear (Fig. 1.43). For illustration, four different analyses of the data were done using PERMANOVA (Fig. 1.44). First an analysis using Type I SS was done, fitting the factor of “Logging” first. Next, an analysis using Type I SS was done again, but this time fitting the factor of “Slope” first. Analyses were then done using each of Type II and Type III SS, in turn. Note that the SS for either factor using Type II SS corresponds to what is obtained if that factor is fitted second in a sequential (Type I) analysis. The Type III SS are different from all others, except for the interaction term, which in all cases was conditioned upon both of the main effects. Note also that the multipliers on individual components in each EMS differ for the different types of SS as well. This, in turn, means that the estimates of components of variation will also differ (Fig. 1.44). Fig. 1.44. PERMANOVA analyses of the two-way crossed unbalanced design for the example of Borneo birds using different Types of SS, as indicated. In summary, the analysis of an unbalanced design using different Types of SS affects the values of (i) the SS themselves (and thus the values of MS and pseudo-F) for individual terms in the model; (ii) the EMS for each term and multipliers of individual components of variation; (iii) the estimates of the sizes of components of variation. Despite all of this, for the present example, the same general conclusions would be obtained, regardless of which Type of SS we had decided to use (Fig. 1.44). There are significant differences among bird assemblages in forests having different logging histories and the slope of the site also has a significant effect. These factors did not interact with one another and logging effects were much larger than slope effects (see Fig. 1.43 and also the estimated components of variation). The extent to which different types of SS will give comparable results will depend on just how unbalanced the design is – greater imbalances will generally lead to greater overlapping regions and thus potentially greater discrepancies. Perhaps the most important point is to recognise how different choices for the type of SS in unbalanced designs correspond conceptually to different underlying hypotheses (see Table 1.5, Fig. 1.41 and Searle (1987) ).   Table 1.5. Tests done using different types of SS in a two-way crossed unbalanced design (cf. Fig. 1.41). The vertical line is to be read as “given”, thus “A | B” should be read as “factor A given factor B”. A comma should be read as “and”. 1.35 Designs with covariates (Holdfast invertebrates, revisited) A topic that is related (perhaps surprisingly) to the topic of unbalanced designs is the analysis of covariance, or ANCOVA. There are some situations where the experimenter, faced with the analysis of a set of data in response to an ANOVA-type of experimental design, would like to take into account one or more quantitative variables or covariates. The essential idea here is that the response variable(s) may be known already to have a relationship with (or to be affected by) some (more or less continuous) quantitative variable. What is of interest then is to analyse the response data cloud given this known existing relationship. That is, one may wish to perform the ANOVA partitioning only after including (i.e., fitting, conditioning upon or taking into account) the covariate(s) in the model. Unfortunately, even if the design is completely balanced, so that terms are independent of one another, this is almost certainly not the case if a covariate is added to the model44. Thus, if one or more covariates are included in an ANOVA design (to yield what is commonly referred to as an ANCOVA model), the SS for individual terms in the model are not independent of one another. This means that the Type of SS must be chosen carefully for these situations, just as for an unbalanced case. Generally, the covariate is to be fit first, with the design factors to be considered given the covariate, so there is a logical sequential order of terms. Thus, Type I SS usually makes the most sense here. However, using Type I SS does mean that the order of the fit of the design factors relative to each other will also matter, so this should be kept in mind. Thus, if Type I SS are to be used, the experimenter might also decide to re-run the analysis, swapping the order of factors in the ANOVA part of the model in order to check whether this affects any essential interpretations of the results. The order in which the individual terms are fit by PERMANOVA can be changed either by changing the order of the rows in the design file or by clicking on the ‘Terms…’ button in the PERMANOVA dialog to yield the ‘Ordered Selection’ window (e.g., Fig. 1.32), then using the up and down arrows to move the relative positions of individual terms that will be included in the sequential fit. An example of a design which might include a covariate in the model is provided by the holdfast invertebrate dataset, seen in the section Nested design. It is known that the community structure of organisms inhabiting holdfasts is affected by the volume of the holdfast habitat itself (e.g., Smith, Simpson & Cairns (1996) ). This stands to reason, given the well-known ecological phenomenon of the species-area relationship (e.g., Arrhenius (1921) , Connor & McCoy (1979) ). Larger holdfasts have not only a larger colonisable area (which might alone be dealt with adequately through a standardisation by total sample abundance), but they also have greater complexity, usually with a greater number of interstices and root-like structures (called haptera) as well. The volume of each holdfast in the study by Anderson, Diebel, Blom et al. (2005) was measured using water displacement, and should really be included as a covariate in the analysis. Fig. 1.45. Draftsman plot of environmental variables for kelp holdfasts. Open up the file hold.pri in the folder ‘HoldNZ’ of the ‘Examples add-on’ directory. Here, we shall analyse only those species (or taxa) that occurred as counts of abundances. We will not include the encrusting organisms that were recorded only as an ordinal measure from 0-3. Choose Select > Variables > Indicator levels > Indicator name: Ordinal > Levels… > (Available Y) & (Include N). Choose Tools > Duplicate and re-name the data file produced as hold.abund. We shall base the analysis on a dissimilarity measure that is a modification of the Gower measure, recently described by Anderson, Ellingsen & McArdle (2006) . This measure is interpretable as the average order-of-magnitude difference in abundance per species, where each change in order of magnitude (i.e. 1, 10, 100, 1000 on a log10 scale) is given the same weight as a change in composition from 0 to 1 (see Anderson, Ellingsen & McArdle (2006) for more details). Choose Analyse > Resemblance > (Analyse between •Samples) & (Measure •More (tab)), click on the tab labeled ‘More’ and choose (•Others: Modified Gower > Modified Gower log base: 10). Save the current workspace as hold_cov.pwk. Environmental data, including depth (in m), density of kelp plants where the holdfast was collected (per m2) and a measure of volume for each holdfast (in ml), are located in the file holdenv.pri, also located in the ‘HoldNZ’ folder. Open up this file within the hold_cov workspace just created and choose Analyse > Draftsman plot to examine the distributions of these environmental variables, focusing on the variable of ‘Volume’ in particular (Fig. 1.45). The distributions of environmental or other quantitative variables to be included as covariates in PERMANOVA should be investigated before proceeding. There are two important reasons for doing this: The model fit by PERMANOVA is linear with respect to the covariate. That is, the relationship between the covariate and the multivariate community structure as represented in the space defined by the dissimilarity measure chosen is linear. Importantly, the relationship between the covariate and the original species (or other response) variables is emphatically not linear unless Euclidean distance is used as the basis of the analysis. So, the experimenter should look (as far as possible) for outliers and other oddities in the multivariate space of the dissimilarity measure chosen (i.e., using ordination) and should also examine the distribution of the covariate as part of a general diagnostic process before proceeding. Permutation of raw data will have inflated type I error if there are outliers in the covariates ( Kennedy & Cade (1996) , Anderson & Legendre (1999) , Anderson & Robinson (2001) ). Thus, permutation of raw data is not allowed as an option if covariates are included in the PERMANOVA model. If the distribution of the covariate is skewed or bimodal, then the user may wish either to transform the variable, or to split the dataset according to any clear modalities observed. See chapter 10 of Clarke & Gorley (2006) regarding the diagnostic process for interpreting and using draftsman plots for environmental variables. The main point is that covariates used in PERMANOVA should show approximately symmetric distributions that are roughly normal, with no extreme outliers. As pointed out by Clarke & Gorley (2006) , it is not necessary to agonise over this issue! Normality is by no means an assumption of the analysis, but outliers can have a strong influence on results. A draftsman plot of the holdfast environmental variables indicates that the distribution of volumes is fairly even, with no obvious skewness or outliers (Fig. 1.45), so no transformation is necessary. Fig. 1.46. MDS of holdfast fauna with volume superimposed as bubbles. An MDS plot of the holdfast fauna with volume superimposed (as bubbles) shows a clear relationship between volume and community structure (Fig. 1.46); change in community structure from left to right across the diagram is associated with increasing volume of the holdfast. Thus, it would make sense to examine the variability in community structure at different spatial scales over and above this existing relationship with volume. Fig. 1.47. PERMANOVA dialog for the analysis of holdfast fauna in response to the fully nested design and including volume as a covariate. To proceed with the analysis, start by selecting the holdenv worksheet, highlight the single variable of Volume, choose Select >Highlighted followed by Tools >Duplicate, then re-name the resulting worksheet vol. It is important that the variable(s) to be used as covariate(s) in the model be contained in a single separate datasheet. It is also important that the labels of the samples in this data sheet match those that are contained in the resemblance matrix, although the particular order of the samples need not be the same in the two files. (PERMANOVA, like other routines in PRIMER, uses a label matching procedure to ensure correct analysis). Next, go to the Modified Gower resemblance matrix calculated from the hold.abund data sheet and create a PERMANOVA design file with three factors: Location, Site and Area, according to the fully nested design outlined in the section Nested design (e.g., Fig. 1.29). Re-name this design file Nested design. Click on the dissimilarity matrix again and choose PERMANOVA+ > PERMANOVA > (Design worksheet: Nested design) & (Covariable worksheet: vol) & (Sums of Squares •Type I SS (sequential)) & (Num. permutations: 9999) & (Permutation method •Permutation of residuals under a reduced model) & ($\checkmark$Use short names), as shown in Fig. 1.47. The results of this analysis show that there is a strong and significant effect of the covariate, i.e., a significant relationship between holdfast volume and community structure as measured by the Modified Gower dissimilarity measure (Fig. 1.48). This is not surprising, given the pattern seen in the MDS plot (Fig. 1.46). Nevertheless, even given the variation among holdfast communities due to volume, significant variability is still detected among the assemblages at each spatial scale in the design: among locations (100’s of km’s), sites (100’s of m’s) and areas (10’s of m’s) (Fig. 1.48). It is possible, in such an analysis, to include interactions between the covariate and each of the other terms in the experimental design. A significant interaction between a factor (such as Locations) and a covariate (such as volume) indicates that the nature of the relationship between the covariate and the multivariate responses differs within different levels of the factor (i.e., that the relationship between volume and holdfast community structure differs at different locations, and so on for other factors in the model). In the case of univariate data with a single factor and one covariate, as analysed in a traditional ANCOVA model, the test of the interaction term between the covariate and the factor is also known as a test for homogeneity of slopes. Sometimes ANCOVA is presented as a model without the interaction term and in these cases it is stated that homogeneity of slopes is an assumption of the analysis. However, if one includes the interaction term in the model, then clearly this type of homogeneity is no longer an assumption, as the model explicitly allows for different slopes for different levels of the factor. See Winer, Brown & Michels (1991) or Quinn & Keough (2002) for more complete discussions and examples of traditional ANCOVA models. In PERMANOVA, tick the box labeled ‘Include interactions’ in the dialog under the specification of the covariable worksheet in order to include these terms in the model. Note that, after ticking this box, you can also still click on the ‘Terms…’ button and change the order in which the terms are fitted and/or exclude particular terms from the model if you wish, whether these be interactions with the covariate or some other terms. An analysis of the holdfast data where interaction terms are included suggests there is actually no compelling reason to include any of these interactions in the model in this particular case (P > 0.18 for all interactions with the covariate, Fig. 1.48). Fig. 1.48. PERMANOVA results of the analysis of holdfast fauna, including volume as a covariate, either (a) without interactions or (b) including interactions. We re-iterate that the most important thing to remember about analyses involving covariates is that the individual terms are not independent of one another, so the Type of SS chosen for the analysis will affect the results, just as it does for unbalanced designs. What is more, the underlying complexity of this particular analysis45 (even without including the interaction terms), which is a mixed model having both random factors and a (fixed) covariate, can be appreciated by considering the EMS’s and the rather impressive gymnastics the program must go through in order to produce tests of individual terms using pseudo-F (Fig. 1.49). For more details on how these pseudo-F ratios were constructed, see the next section, Linear combinations of MS. Fig. 1.49. PERMANOVA results of the analysis of holdfast fauna, including volume as a covariate and showing the EMS, construction of pseudo-F and components of variation for each term in the model. 44 The reason non-independence is introduced becomes clearer perhaps when we consider that if the covariate is fitted first and different ranges of the continuous covariate occur in different cells, then it cannot be independent of existing factors. 45 This complexity was, unfortunately, not recognised by the authors of the original paper ( Anderson, Diebel, Blom et al. (2005) ). The analyses they presented took a “naïve” view and conditioned on the covariate without considering the effects this would have on the EMS’s and pseudo-$F$ ratios for the other factors in the model. Making mistakes is, however, part of doing science and we would be nowhere if we did not make mistakes, acknowledge our errors and learn from them! 1.36 Linear combinations of mean squares (NZ fish assemblages) Several aspects of the above analysis demonstrate its affinity with unbalanced designs. Note that: (i) the multipliers for components of variation in each EMS are not whole numbers; and (ii) the multipliers for a given component of variation are not the same in different EMS’s. As a consequence of these two things, for many of the terms in the model, there is no other single term that, alone, can provide a MS which can act as a denominator to yield a correct pseudo-F ratio. What this generally means for such cases is that some linear combination of mean squares must be used in order to construct a test of the given null hypothesis of interest. Depending on how these linear combinations are constructed, this can also mean that even the degrees of freedom used for the tests are not whole numbers (Fig. 1.49)! The good news is that the PERMANOVA routine is actually equal to this (rather horrendous) task; (i) it determines the correct EMS’s for every term; (ii) it calculates the correct linear combinations of mean squares required to construct pseudo-F ratios to test each hypothesis; and (iii) it determines the correct distributions of each pseudo-F ratio under each relevant null hypothesis using permutations, producing accurate P-values. The routine is also not bothered by non-integer degrees of freedom. The exchangeable units for a given test (e.g., Anderson & ter Braak (2003) ) are chosen using what might be called a “highest-order term” approach46. That is, consider the linear combination of mean squares required for a given test. Of the terms giving rise to those mean squares, the term in the denominator which is of highest order (excluding continuous covariates and the residual) is used to identify exchangeable units for that particular test. For example, the terms whose mean squares are included in the denominator for the pseudo-F ratio test of Location in Fig. 1.49 are: Si(Lo), Ar(Si(Lo)) and Res. The highest-order term (excluding the residual) is Ar(Si(Lo)). Therefore, all replicates within an area will be kept together as a group under permutation and the 16 different areas will be the units permuted for the test of Locations in this design. Importantly, after each permutation, the full pseudo-F ratio is constructed according to the requirements of the numerator and denominator, even if either one or both of these are linear combinations of mean squares. The need to construct linear combinations of mean squares happens not just in highly complex models with covariates, as described above. It is also (much more commonly) required for certain terms even in many balanced designs, especially multi-way designs involving more than one random factor47. For example, consider a study of temperate rocky reef fish assemblages as described by Anderson & Millar (2004) . The study consisted of surveys of fish biodiversity, where abundances of fish species were counted in 10 transects (25 m × 5 m) sampled by SCUBA divers from each of 4 sites in each of 2 habitats at each of 4 locations along the northeast coast of New Zealand. These surveys have been done each year in the austral summer. The data provided in the file fishNZ.pri (located in the folder ‘FishNZ’ in the ‘Examples add-on’ directory) are sums across the 10 transects within each site for each of p = 58 fish species from surveys conducted in each of two years: 2004 and 200548. The experimental design here is: Factor A: Year (random with a = 2 levels: 4 = 2004 and 5 = 2005). Factor B: Location (random with b = 4 levels: B = Berghan Point, H = Home Point, L = Leigh and A = Hahei). Factor C: Habitat (fixed with c = 2 levels: b = urchin-grazed ‘barrens’ and k = kelp forest). It is of interest to test the null hypothesis of no difference between the two habitats in fish assemblages. It is also (secondarily) of interest to test and quantify the variability among years and among locations in fish community structure. An analysis of the data according to the above experimental design using PERMANOVA has been done on the basis of the scaled binomial deviance dissimilarity measure (see Anderson & Millar (2004) for a description of this measure), yielding the results shown in Fig. 1.50. Fig. 1.50. Analysis of New Zealand temperate reef fish according to the three-way mixed model. Note the linear combination of mean squares needed to test the habitat main effect: ‘Ha’. What is of immediate interest to us here is the test for the main effect of ‘Habitat’ or ‘Ha’ in the output49. The EMS for this factor is: Now, to construct pseudo-F, we need to find a denominator mean square whose expectation is equal to the above when the null hypothesis that S(Ha) = 0 is true. That is, we need a denominator whose expectation is: However, there is clearly no single term that, alone, can perform this duty, because we require both the Lo×Ha and the Ye×Ha components of variation to appear here. We can see, however, that the term we seek can be obtained by constructing a linear combination of mean squares: It is desirable, however, not to include mean squares negatively, as this could generate negative pseudo-F ratios, which are not really sensible (e.g., Searle, Casella & McCulloch (1992) ). Thus, re-arranging the above (so that all mean square terms appear positively), we have: or Accordingly, if we wish to construct a pseudo-F ratio where the numerator and denominator will have the same expectation if H$_0$: S(Ha) = 0 is true, and which gets large with increases in the size of S(Ha) alone then we can use: This is precisely the pseudo-F ratio constructed by the PERMANOVA program, as stated in the results in the line for ‘Ha’ under the heading ‘Construction of Pseudo-F ratio(s) from mean squares’ (Fig. 1.50). The exchangeable units for this test, using the “highest-order term” approach, are obtained by considering all of the terms involved in the construction of pseudo-F (apart from the term being tested, which is ‘Ha’ here) and determining the one with the highest order. The terms involved are: Ye×Lo×Ha, Ye×Ha and Lo×Ha. The term with the highest order is Ye×Lo×Ha, so the exchangeable units for this test are the a × b × c = 2 × 4 × 2 = 16 cells. So, the n = 4 sites occurring within each of those cells will be permuted together as a unit and the full pseudo-F will be re-constructed after each permutation in order to test the term ‘Ha’ in this case 50. The attentive user will notice the large P-values (> 0.25) associated with a couple of the terms in the model, namely the ‘Ye×Ha’ and ‘Lo×Ha’ interaction terms. If possible, it is desirable to pool (remove) such terms, thereby simplifying the construction of pseudo-F and potentially increasing power (see the section Pooling or excluding terms). Only one term should be removed at a time, beginning with the one having the smallest mean square (e.g., Fletcher & Underwood (2002) ), as this will affect the tests and estimated components of variation for other terms in the model. By following this procedure, it is found that ‘Ye×Ha’ and ‘Lo×Ha’ can each be removed (in that order), yielding the results shown in Fig. 1.51. Note that, in this case, after pooling, no linear combinations of mean squares were required for any of the remaining tests. Fig. 1.51. PERMANOVA of the New Zealand fish assemblages after pooling of terms, showing the dialog used at the second step in the pooling procedure, when the term ‘Lo×Ha’ was pooled along with ‘Ye×Ha’, which had previously been pooled at step one. In traditional univariate ANOVA, the potential use of pooling for these kinds of situations was very important. This is because the construction and subsequent testing of traditional F ratios using linear combinations of mean squares in the numerator and denominator is fraught with difficulties, primarily because ratios of linear combinations of mean squares (called “quasi” F ratios by Quinn & Keough (2002) , see also Blackwell, Brown & Mosteller (1991) ) no longer have (known) F distributions under a true null hypothesis, even when the usual assumptions of normality, homogeneity, etc. are fulfilled (see also Searle, Casella & McCulloch (1992) ). Complicated approximations have therefore been suggested in order to obtain P-values for these cases (e.g., Satterthwaite (1946) , Gaylor & Hopper (1969) ), in the event that pooling was not possible. One of the most important advantages of the PERMANOVA routine is that it uses permutation tests to obtain P-values. Thus, as long as: (i) the test statistic is constructed correctly in the sense that it isolates the term of interest under the null hypothesis; and (ii) the permutations are done so as to create alternative realisations under a true null hypothesis by permuting appropriate exchangeable units, then the calculations of the P-values are correct and can be used for valid inference. This is true whether or not the corresponding traditional univariate test would be able to be done at all using more traditional theoretical approaches. The more general unified approach of the new PERMANOVA software caters even for these situations (such as the need for linear combinations of mean squares), where the traditional tests would be very difficult (sometimes even impossible) to formulate. Also, of course, PERMANOVA can be implemented on distance or dissimilarity matrices which have been calculated from either univariate or multivariate data. 46 The order of a term is defined here as follows: a main effect (e.g., A, B, …) is of first order, a two-way interaction (e.g., A×B) is of second order, a three-way interaction (e.g., A×B×C) is of third order, and so on. Also, a nested term, such as B(A), is of second order, while C(A×B) or C(B(A)) would both be of third order, etc. 47 The need to use linear combinations of mean squares to construct appropriate pseudo-F statistics will also arise much more commonly (or for more of the terms in the model) if the constraint that fixed effects should sum to zero across levels of random factors in mixed interactions is not applied (i.e., if one chooses to remove the $\checkmark$ in front of the option ‘Fixed effects sum to zero’ in the PERMANOVA dialog). 48 Note that these are not the same years as those analysed by Anderson & Millar (2004) , but are data from more recent surveys. 49 See the section Inference space and power regarding the logic and hypotheses underlying tests of fixed main effects even in the presence of potentially non-zero interactions with random factors. 50 Although some preliminary simulation work has indicated that this “highest-order term” approach works well in trial cases in terms of maintaining rates of type I error at chosen significance levels, a more complete study of this rather complex issue of appropriate exchangeable units for F ratios involving linear combinations of mean squares would be welcome. 1.37 Asymmetrical designs (Mediterranean molluscs) Although a previous section has been devoted to the analysis of unbalanced designs, there are some special cases of designs having missing cells which deserve extra attention. Such designs are commonly referred to as asymmetrical designs, and consist essentially of there being different numbers of levels of a nested factor within each different level of an upper-level factor. Important examples include the asymmetrical designs that can occur in studies of environmental impact. Here, there might only be a single site that is impacted, whereas there might be multiple control (or unimpacted) sites (e.g., Underwood (1992) , Underwood (1994) , Glasby (1997) ). The reason for asymmetrical designs arising frequently in the analysis of environmental impacts is that it is generally highly unlikely that an impact site of a particular type (e.g., an oil spill, a sewage outfall, the building of a particular development, etc.) will be replicated, whereas there is often no reason not to include multiple replicate control sites (at a given spatial scale) against which changes at the (purportedly) impacted site might be measured (e.g., Underwood (1992) , Underwood (1994) ). On the face of it, the experimenter might consider that such a design presents a severe case of imbalance, where not all cells are filled. In actual fact, this is not really the case, because it is only the number of levels of the nested factor that are unequal, and actually all of the terms in the partitioning will be independent of one another, just as they would be in a balanced design. An example of an asymmetrical design is provided by a study of subtidal molluscan assemblages in response to a sewage outfall in the Mediterranean ( Terlizzi, Scuderi, Fraschetti et al. (2005) ). The study area is located along the south-western coast of Apulia (Ionian Sea, south-east Italy). Sampling was undertaken in November 2002 at the outfall location and at two control or reference locations. Control locations were chosen at random from a set of eight possible locations separated by at least 2.5 km and providing comparable environmental conditions to those occurring at the outfall (in terms of slope, wave exposure and type of substrate). They were also chosen to be located on either side of the outfall, to avoid spatial pseudo-replication. At each of the three locations, three sites, separated by 80 - 100 m were randomly chosen. At each site, assemblages were sampled at a depth of 3 - 4 m on sloping rocky surfaces and n = 9 random replicates were collected (each replicate consisted of scrapings from an area measuring 20 cm × 20 cm), yielding a total of N = 81 samples. The experimental design is: Factor A: Impact versus Control (‘IvC’, fixed with a = 2 levels: I = impact and C = control). Factor B: Location (‘Loc’, random, nested in IvC with b = 2 levels nested in C and 1 level nested in I, labeled as numbers 1, 2, 3). Factor C: Site (random, nested in Loc(IvC) with c = 3 levels, labeled as numbers 1-9). A schematic diagram (Fig. 1.52) helps to clarify this design, and why it is considered asymmetrical. Fig. 1.52. Schematic diagram of the asymmetrical design for Mediterranean molluscs. The data for this design are located in the file medmoll.pri in the ‘MedMoll’ folder of the ‘Examples add-on’ directory. As in Terlizzi, Scuderi, Fraschetti et al. (2005) , to visualise patterns among sites, we may obtain an MDS plot of the averages of samples at the site level. Go to the worksheet containing the raw data and choose Tools > Average > (Samples •Averages for factor: Site), then Analyse > Resemblance > (Analyse between: •Samples) & (Measure •Bray-Curtis), followed by Analyse > MDS. Fig. 1.53. MDS of site averages for Mediterranean molluscs, where I = impact and C = control locations. There is apparent separation between the (averaged) assemblages at sites from the impact location compared to the controls on the MDS plot (Fig. 1.53), but to test this, we shall proceed with a formal PERMANOVA analysis. The difference, if any, between impact and controls must be compared with the estimated variation among control locations. Next, calculate a Bray-Curtis resemblance matrix directly from the original medmoll.pri raw data sheet (i.e., not the averaged data), then create a PERMANOVA design file according to the above experimental design, re-name the design file Asymmetric and proceed to run the analysis by choosing the following in the PERMANOVA dialog: (Design worksheet: Asymmetric) & (Test: •Main test) & (Sums of Squares •Type III (partial) ) & (Num. permutations: 9999) & (Permutation method •Permutation of residuals under a reduced model) & ($\checkmark$Do Monte Carlo tests) & ($\checkmark$Fixed effects sum to zero). For clarity in viewing these results, un-check (i.e. remove the $\checkmark$ from) the option to ‘Use short names’. The results indicate that there is significant variability among sites, but variation among control locations is not detected over and above this site-level variability (i.e. the term ‘Loc(IvC)’ is not statistically significant P > 0.12, Fig. 1.54). There is also apparently a significant difference in the structure of molluscan assemblages at the impact location compared to the control locations (P(MC) = 0.036). Note that, in the absence of any replication of outfall (impact) locations, the only basis upon which location-level variability may be measured is among the (in this case only two) control locations. We should refrain from going overboard in the extent of our inferences here – it is inappropriate to place strong importance on an approximate MC P-value which relies on asymptotic theory (i.e. its accuracy gets better as sample size increases) yet was obtained using only two control locations and hence only 1 df in the denominator. Furthermore, of course, in the absence of any data from before the outfall was built, it is not possible to infer that the difference between the impact and the controls detected here was necessarily caused by the sewage outfall. It is also impossible to know whether this difference is something that has persisted or will persist in time. What we can say is that the molluscan assemblages at the outfall at the time the data were collected were indeed distinct from those found at the two control locations in the area sampled at that time, as observed in the MDS plot (Fig. 1.53). Fig. 1.54. PERMANOVA analysis of an asymmetrical design for mediterranean molluscs. Asymmetrical designs such as this one have often previously been analysed and presented by partitioning overall location effects into two additive pieces: (i) the SS due to the contrast of the impact vs the controls and (ii) the SS due to the variability among the controls (e.g., Underwood (1994) , Glasby (1997) ). The reason for this has largely been because of the need for experimenters to utilise software to analyse these designs which did not allow for different numbers of levels of nested factors in the model. PERMANOVA, however, allows the direct analysis of each of the relevant terms of an asymmetrical design such as this, without having to run more than one analysis, and without any other special manipulations or calculations. It is not necessary for the user to specify contrasts in PERMANOVA in order to analyse an asymmetrical design. It is essential to recognise the difference between the design outlined above, which is the correct one, and the following design, which is not: Factor A: Locations (fixed or random(?) with a = 3 levels and special interest in the contrast of level 3 (impact) versus levels 1 and 2 (controls), Factor B: Sites (random, b = 3 levels, nested in Locations) Be warned! If an asymmetrical design is analysed using contrasts in PERMANOVA, then although the SS of the partitioning will be correct, the F ratios and P-values for some of the terms will almost certainly be incorrect! Note that the correct denominator MS for the test of ‘IvC’ must be at the right spatial scale, i.e. at the scale of locations, even though, in the present design (with only one impact location), our only measure of location-level variability comes from the control locations. Recall that a contrast, as a one-degree-of-freedom component partitioned from some main effect, will effectively use the same denominator as that used by the main effect (e.g., see the section Contrasts), which is not the logical choice in the present context. A clear hint to the problem underlying this approach is apparent as soon as we try to decide whether the ‘Locations’ factor should be fixed or random. It cannot be random, because the contrast of ‘IvC’ is clearly a fixed contrast of two states we are interested in. However, neither can it be fixed, because the control locations were chosen randomly. They are intended to represent a population of possible control locations and their individual levels are not of any interest in and of themselves. Once again, the rationale for analysing asymmetrical designs using contrasts in the past (and subsequently constructing the correct F tests by hand, see for example, Glasby (1997) ) was because software was not widely available to provide a direct analysis of the true design. Although an asymmetrical design might appear to be unbalanced (and it is, in the sense that the amount of information, or number of levels, used to measure variability at the scale of the nested factor is different within different levels of the upper-level factor), it does not suffer from the issue of non-independence which was described in the section on Unbalanced designs. In fact, all of the terms in an asymmetrical design such as this are completely orthogonal (independent) of one another. So, it does not matter which Type of SS the user chooses, nor in which order the terms are fitted – the same results will be obtained. (This is easily verified by choosing to re-run the above analysis using, for example, Type I SS instead.) 1.38 Environmental impacts Some further comments are appropriate here regarding experimental designs to detect environmental impact ( Green (1979) , Underwood (1991) , Underwood (1992) , Underwood (1994) ). These designs generally include measurements of a response variable of interest before and after a potential impact from one or more control site(s) and from the purportedly impacted site(s). These are referred to as “BACI” designs51 as an acronym for the two levels in each of the two major factors of the design: the temporal factor (“before”/“after”) and the spatial factor (“control”/“impact”). Importantly, a significant interaction between these two factors would (potentially) lead to inferences regarding significant impact, so the ability formally to test the interaction term(s) in such models is paramount here. With regard to extending non-parametric multivariate hypothesis-testing methods to handle such designs, Clarke (1993) conceded: “This would appear to defy development within the similarity-based framework… and must be accepted as a limitation of the current methodology, though there is clearly scope for further study here” (p. 138). Indeed, further study led to the development of PERMANOVA, which allows (under slightly less general conditions, e.g., the approach is no longer fully non-parametric) tests of interaction terms for any multi-factorial model on the basis of any resemblance measure of choice, with P-values obtained by appropriate permutation techniques. Its new implementation as an add-on to PRIMER also now allows appropriate analyses of asymmetrical designs (i.e., in the case of there being multiple control sites, but only one impact site – see the section Asymmetrical designs), as required. In essence, as PERMANOVA can be used to analyse multivariate responses to any ANOVA model, it therefore can be used readily for the analysis of assemblages in response to either BACI or beyond BACI experimental designs in environmental impact studies. A further contribution of Underwood ( Underwood (1991) , Underwood (1992) , Underwood (1994) ) in the area of experimental designs to detect environmental impacts (in addition to proposals to extend the basic model to include multiple sites and times of sampling and thus avoid pseudo-replication) was to propose the use of two-tailed F-tests to detect potential impacts on variability in the response variable, and not just to detect changes in means. Recently, Terlizzi, Anderson, Fraschetti, S. et al. (2007) have implemented these ideas for multivariate data (albeit not in the context of environmental impact, but to investigate patterns along a depth gradient). More particularly, they used bootstrapping to place confidence intervals on differences in the sizes of multivariate components of variation. Although the current PERMANOVA+ software does not implement this approach, direct measures of the components of variation for each term in the model are provided as part of the PERMANOVA output. Thus, the sizes of multivariate (or univariate) components of variation for different sub-sets of the data may be calculated directly in this way. Another approach to addressing hypotheses concerning variability in assemblage structure is to use the PERMDISP routine (see chapter 2), which can be applied either to individual replicates within groups (or cells), or can be applied to centroids (calculated from PCO axes when non-Euclidean measures are being used, see chapter 3) to compare dispersions for higher-level factors in more complex designs. 51 Much to the delight of Italian ecologists!