NUMERICAL ANALYSIS OF BIOLOGICAL AND ENVIRONMENTAL DATA

Transcript NUMERICAL ANALYSIS OF BIOLOGICAL AND ENVIRONMENTAL DATA

NUMERICAL ANALYSIS OF BIOLOGICAL AND ENVIRONMENTAL DATA

Lecture 7.

Direct Gradient Analysis

DIRECT GRADIENT ANALYSIS

Interpretation of ordination axes with external data Canonical or constrained ordination techniques (= direct gradient analysis) Canonical correspondence analysis (CCA) Introduction Basic terms and ordination plots Other topics in CCA Robustness Scaling and interpretation of CCA plots Example Redundancy analysis (RDA) (= constrained PCA) Scaling and interpretation of RDA plots Statistical testing of constrained ordination axes Partial constrained ordinations Partial ordinations Partitioning variance Environmental (predictor) variables and their selection Canonical correlation analysis Distance-based redundancy analysis Canonical analysis of principal co-ordinates Principal response curves Polynomial RDA and CCA CCA/RDA as predictive tools Non-linear canonical analysis of principal co ordinates Canonical Gaussian ordination Constrained additive ordination

CANODRAW

BASIS OF CLASSICAL ORDINATION INTERPRETATION AND ENVIRONMENT

We tend to assume that biological assemblages are controlled by environment, so: 1. Two sites close to each other in an indirect ordination are assumed to have similar composition, and 2. if two sites have similar composition, they are assumed to have similar environment. In addition: 3. Two sites far away from each other in ordination are assumed to have dissimilar composition, and thus 4. if two sites have different composition, they are assumed to have different environment.

J. Oksanen (2002)

DUNE-MEADOW DATA Values of environmental variables

and Ellenberg’s indicator values of species written alongside the ordered data table of the Dune Meadow Data, in which species and sites are

arranged in order of their scores on

the second DCA axis. A1: thickness of A1 horizon (cm), 9 meaning 9cm or more; moisture: moistness in five classes from 1 = dry to 5 = wet; use: 1 = hayfield, 2 = a mixture of pasture and hayfield, 3 = pasture; manure: amount applied in five classes from 0 = no manure to 5 = heavy use of manure. The meadows are classified by type of management: SF, standard farming; BF, biological farming; HF, hobby farming; NM, nature management; F, R, N refer to Ellenberg’s indicator values for moisture, acidity and nutrients, respectively .

DCA axis 2 DCA axis 1

The amount of manure written on the DCA ordination. The trend in the amount across the diagram is shown by an arrow, obtained by a multiple regression of manure on the site scores of the DCA axes. Also shown are the mean scores for the four types of management, which indicate, for example, that the nature reserves (NM) tend to lie at the top of the diagram. Ez=b 0 + b 1

1 + b 2

2 Angle (  )with axis 1 = arctan(b 2 / b 1 )

Site scores of the second DCA axis plotted against the amount of manure.

Correlation coefficients (100  r) of the environmental variables for the four first DCA axes for the Dune Meadow Data

Variable

1 A1 2 moisture 3 use 4 manure 5 SF 6 BF 7 HF 8 NM Eigenvalue

35 6 22 -28 -22 21 0.54

24 57 -21

-68

-29 -24 -26

0.29

Axes

7 -3 -7 5 39

7 -55 17 0.08

9 -7 -5 -64 -60 22 -14 56 0.05

Multiple regression of the first CA axis on four environmental variables of the dune meadow data, which shows that moisture contributes significantly to the explanation of the first axis, whereas the other variables do not.

Term

constant A1 moisture use manure

Parameter

Estimate s.e.

–2.32

0.50 0.14

0.38

0.08

0.09

0.31

–0.00

0.22

0.12

–4.62

1.71

4.08

1.37

–0.01

ANOVA table

d.f.

Regression Residual 15 Total 19

s.s.

4 6.2

23.2

m.s.

17.0

0.41

1.22

4.25

10,6 R 2 = 0.75

R 2 adj = 0.66

1 = b

0 + b 1 x 1 + b 2 x 2 + ...b

n x n

CA axis 1 environmental variables

= environmental variables

TWO-STEP APPROACH OF INDIRECT GRADIENT ANALYSIS

Standard approach to about 1985: started by D.W. Goodall in 1954 Limitations: (1) environmental variables studied may turn out to be poorly related to the first few ordination axes.

(2) may only be related to 'residual' minor directions of variation in species data.

(3) remaining variation can be substantial, especially in large data sets with many zero values.

(4) a strong relation of the environmental variables with, say, axis 5 or 6 can easily be overlooked and unnoticed.

Limitations overcome by canonical or constrained ordination techniques = multivariate direct gradient analysis.

CANONICAL ORDINATION TECHNIQUES

Ordination and regression in one technique – Cajo ter Braak 1986 Search for a weighted sum of environmental variables that fits the species best, i.e.

that gives the maximum regression sum of squares

Ordination diagram 1) patterns of variation in the species data 2) main relationships between species and each environmental variable Redundancy analysis  constrained or canonical PCA Canonical correspondence analysis (CCA)  (Detrended CCA)  constrained CA constrained DCA Axes constrained to be linear combinations of environmental variables.

In effect PCA or CA with one extra step: Do a multiple regression of site scores on the environmental variables and take as new site scores the fitted values of this regression.

Multivariate regression of Y on X.

PRIMARY DATA IN GRADIENT ANALYSIS

Abundances or +/ variables Response variables PLUS

Values Classes

Predictor or explanatory variables

Artificial example of unimodal response curves of five species (A-E) with respect to standard ised environmental variables showing different degrees of separation of the species curves moisture linear combination of moisture and phosphate CCA linear combination a: Moisture b: Linear combination of moisture and phosphate, chosen a priori c: Best linear combination of environmental variables, chosen by CCA.

Sites are shown as dots, at y = 1 if Species D is present and at y = 0 if Species D is absent

Combinations of environmental variables

e.g.

3 x moisture + 2 x phosphate all possible linear combinations

x j



c o



.....

z j

= environmental variable at site

j c

= weights

x j

= resulting ‘compound’ environmental variable CCA selects linear combination of environmental variables that maximises dispersion of species scores, i.e. chooses the best weights (

c i

) of the environmental variables.

ALTERNATING REGRESSION ALGORITHMS

- CA - DCA - CCA Algorithms for (A) Correspondence Analysis, (B) Detrended Correspondence Analysis, and (C) Canonical Correspondence Analysis, diagrammed as flowcharts. LC scores are the linear combination site scores, and WA scores are the weighted averaging scores.

REF

CANONICAL CORRESPONDENCE ANALYSIS

Algorithm

1) Start with arbitrary, but unequal, site scores x

2) Calculate species scores by weighted averaging of site scores.

u k



i n

  1

y ik x i

i n

  1

y ik

3) Calculate new site scores by weighted averaging of species scores.

x i

 

k m

  1

y ik u k

k m

  1

y ik

[So far, two-way weighted average algorithm of correspondence analysis].

REF REF REF

REF REF 4) Obtain regression coefficients of site scores on the environmental variables by weighted multiple regression.

 

  1

R x

 where

b Z R

and x* are column vectors is environmental data n x (q +1) is n x n matrix with site totals in diagonal 5) Calculate new site scores



0  

6) Centre and standardise site scores so that:

i n

  1

y i



x i

 0 and

i n

  1

y i



x i

2  1 7) Stop on convergence, i.e. when site scores are sufficiently close to site scores of previous iteration. If not, go to 2.

REF REF

CANONICAL OR CONSTRAINED CORRESPONDENCE ANALYSIS (CCA)

Ordinary correspondence analysis gives: 1. Site scores which may be regarded as reflecting the underlying gradients.

2. Species scores which may be regarded as the location of species optima in the space spanned by site scores.

Canonical or constrained correspondence analysis gives in addition: 3. Environmental scores which define the gradient space.

These optimise the interpretability of the results.

J. Oksanen (2002)

BASIC TERMS

Eigenvalue = Maximised dispersion of species scores along axis. In CCA usually smaller than in CA. If not, constraints are not useful.

Canonical coefficients = ‘Best’ weights or parameters of final regression.

Multiple correlation of regression = Species–environment correlation. Correlation between site scores that are linear combinations of the environmental variables and site scores that are WA of species scores. Multiple correlation from the regression. Can be high even with poor models. Use with care!

Species scores = WA optima of site scores, approximations to Gaussian optima along individual environmental gradients.

Site scores = Linear combinations of environmental variables (‘fitted values’ of regression) (1).

Can also be calculated as weighted averages of species scores that are themselves WA of site scores (2).

(1) LC scores are predicted or fitted values of multiple regression with constraining predictor variables 'constraints'.

(2) WA scores are weighted averages of species scores.

Generally always use (1) unless all predictor variables are 1/0 variables.

SUMMARY OF DUNE MEADOW DATA

Dune Meadow Data. Unordered table that contains 20 relevées (columns) and 30 species (rows). The right-hand column gives the abbreviation of the species names listed in the left-hand column; these abbreviations will be used throughout the book in other tables and figures. The species scores are according to the scale of van der Maarel (1979b).

Environmental data of 20 relevées from the dune meadows

Use categories: 1 = hay 2 = intermediate 3 = grazing * = mean value of variable Sample number 12 13 14 15 16 17 18 19 20 1 2 3 4 5 6 7 8 9 10 11 A1 horizon Moisture Management class type Use 4 5 5 5 5 2 1 5 5 1 1 2 2 1 1 1 5 4 2 1 5.8

6.0

9.3

11.5

5.7

4.0

4.6* 3.7

3.5

2.8

3.5

4.3

4.2

6.3

4.3

2.8

4.5

3.7

3.3

3.5

2 2 3 2 3 1 1 1 1 2 2 2 2 1 2 3 3 1 1 3 SF SF NM NM SF NM NM NM NM SF BF SF SF HF HF HF HF HF BF BF Manure class 0 3 0 0 0 0 2* 3 0 4 2 4 4 2 2 3 3 1 1 1

DCA axis 2

DCA

DCA ordination diagram of the Dune Meadow Data

DCA axis 1

Axis Two  2 = 0.29

DCA

Axis One  1 = 0.54

Correlations of environmental variables with DCA axes 1 and 2

CCA

CCA of the Dune Meadow Data. a: Ordination diagram with environmental variables represented by arrows. the c scale applies to environmental variables, the u scale to species and sites. the types of management are also shown by closed squares at the centroids of the meadows of the corresponding types of management.

DCA CCA



0.54

0.46



0.40

0.29

R axis 1 R axis 2

0.87

0.83

0.96

0.89

CANONICAL CORRESPONDENCE ANALYSIS

Canonical correspondence analysis: canonical coefficients (100 x c) and intra set correlations (100 x r) of environmental variables with the first two axes of CCA for the Dune Meadow Data. The environmental variables were standardised first to make the canonical coefficients of different environmental variables comparable. The class SF of the nominal variable 'type of management' was used as a reference class in the analysis. Variable Coefficients Axis 1 Axis 2 A1 9 Moisture 71 Use 25 Manure SF BF HF NM -7 -9 20 18 -37 -29 5 -27 16 19 92 Correlations Axis 1 Axis 2 57

21 -30 16 -37 -36 56 -17 -14 -41

-79

-70 15 -12

a b

b: Inferred ranking of the species along the variable amount of manure, based on the biplot interpretation of Part a of this figure.

BIPLOT PREDICTION OF ENVIRONMENTAL VARIABLES

• Project a site point onto environmental arrow: predict its environmental value • Exact with two constraints only • Projections are exact only in the full multi-dimensional space. Often curved when projected onto a plane Modified from J. Oksanen (2002)

REF

CLASS VARIABLES

• Class 1/0 variables usually represented as 'dummy' variables: make m - 1 indicator variables out of m levels (moisture classes 1, 2, 4, 5). Ignore class 3 • Scoring: 1 if site belongs to the class, 0 otherwise • One dummy less than levels, because one is redundant. If it does not belong to any of the m - 1 classes, it must belong to the remaining one.

• Ordered factors such as the four moisture classes can be expressed as polynomial constraints. A four-level ordered factor can be expressed in three 'dummy' variables - linear, quadratic, and cubic effects. Plot as biplot arrows. Helps to find when one can replace a multilevel factor by a single continuous variable (e.g. Moisture Linear) REF J. Oksanen (2002) REF

CCA: JOINT PLOTS AND TRIPLOTS

• You may have in a same figure • WA scores of species • WA or LC scores of sites • Biplot arrows or class centroids of environmental variables • In full space, the length of an environmental vector is 1: When projected onto ordination space • Length tells the strength of the variable • Direction shows the gradient • For every arrow, there is an equal arrow to the opposite direction, decreasing direction of the gradient • Project sample points onto a biplot arrow to get the expected value • Class variables coded as dummy variables • Plotted as class centroids • Class centroids are weighted averages • LC score shows the class centroid, WA scores show the dispersion of the centroid • With class variables only: Multiple Correspondence Analysis or Analysis of Concentration

Summary

CANOCO

Axes

Axes

Eigenvalues Species-environment correlations Cumulative percentage variance

1 2 3 4

.461 .298 .160 .134

.958 .902 .855 .889

of species data 21.8 35.9 43.5 49.8

of species-environment 37.8 62.3

75.4 86.3

relation Sum of all unconstrained eigenvalues = inertia Sum of all canonical eigenvalues = species-environment relation

Total inertia

2.115

'Fitted' species data 2.115

1.220

 1  canonical   2 eigenvalue s  100 Rules of thumb:  >0.30 strong gradient  >0.40 good niche separation of species

CONVERGENCE CRITERIA IN EIGENVALUE EXTRACTION

Oksanen & Minchin (1997) J. Vegetation Science 8, 447–454

DECORANA CANOCO

version 2 version 3.1

Tolerance

0.0001

0.00005

CANOCO

version 3.12a (‘STRICT’ CRITERIA) 0.000005

CANOCO

version 4 0.000005

Number of iterations

10 15 999 999

REPRESENTATION OF CLASS VARIABLES (1/0) IN CCA

1. Make class centroids as distinct as possible 2. Make clouds about centroids as compact as possible • Success   • LC scores are the class centroids: the expected locations, WA scores are the dispersion of the centroid • If high  , WA scores are close to LC scores • With several class variables, or together with continuous variables, the simple structure can become blurred J. Oksanen (2002)

Canonical correspondence analysis

Unimodal curves for the expected abundance response (y) of four species against an environmental gradient or variable (x). The optima, estimated by weighted averages, (u) [k=1,2,3], of three species are indicated. The curve for the species on the left is truncated and therefore appears monotonic instead of unimodal; its optimum is outside the sampled interval but, its weighted average is inside. The curves drawn are symmetric, but this is no strict requirement for CCA.

7) t-values of canonical coefficients or forward selection option in CANOCO to find minimal set of significant variables that explain data about as well as full set.

8) Can be sensitive to deviant sites, but only if there are outliers in terms of both species composition and environment. CCA usually much more robust than CA.

9) Can regard CCA as a display of the main patterns in weighted averages of each species with respect to the environmental variables.

Intermediate between CA and separate WA regressions for each species.

Separate WA regressions  NICHE.

point in q-dimensional space of environmental variables. CCA attempts to provide a low-dimensional representation of this niche.

10) ‘Dummy’ variables (e.g. group membership or classes) as environmental variables. Shows maximum separation between pre-defined groups.

11) ‘Passive’ species or samples or environmental variables. Some environmental variables active, others passive e.g. group membership – active environmental variables – passive 12) CANOCO ordination diagnostics fit of species and samples pointwise goodness of fit can be expressed either as residual distance from the ordination axis or plane, or as proportion of projection from the total chi-squared distance species tolerances, sample heterogeneity

Passive ‘fossil’ samples added into CCA of modern data

Canonical correspondence analysis (CCA) time-tracks of selected cores from the Round Loch of Glenhead; (a) K5, (b) K2, (c) K16, (d) k86, (e) K6, (f) environmental variables. Cores are presented in order of decreasing sediment accumulation rate.

13) Indicator species 14) Behaves well with simulated data.

M W Palmer (1993) Ecology 74, 2215–2230 Copes with skewed species distributions ‘noise’ in species abundance data unequal sampling designs highly intercorrelated environmental variables situations when not all environmental factors are known

Site scores along the first two axes in CCA and DCA ordinations, with varying levels of quantitative noise in species abundance. Quantitative noise was not simulated. The top set represents CCA LC scores and environmental arrows, the middle represents CCA WA scores, and the bottom represents DCA scores. Sites with equal positions along the environmental gradient 2 are connected with lines to facilitate comparisons.

Palmer, M.W. (1993) Ecology 74, 2215–2230

..continued

ROBUSTNESS OF CANONICAL CORRESPONDENCE ANALYSIS

Like all numerical techniques, CCA makes certain assumptions, most particularly that the abundance of a species is a unimodal function of position along environmental gradient. Does not have to be symmetric unimodal function.

Simulated data Palmer 1993 – CCA performs well even with highly skewed species distributions.

‘Noise’ in ecological data – errors in data collection, chance variation, site-specific factors, etc. Noise is also regarded as ‘unexplained’ or ‘residual’ variance. Regardless of cause, noise does not affect seriously CCA.

‘Noise’ in environmental data is another matter. In regression, assumed that predictor variables are measured without error. CCA is a form of regression, so noise in environmental variables can affect CCA.

Highly correlated environmental variables, e.g. soil pH and Ca. Species distributions along Ca gradient may be identical to distributions along pH gradient, even if one is ecologically unimportant. Species and object arrangement in CCA plot not upset by strong inter-correlations. CCA (like all other regression techniques) cannot tell us which is the ‘real’ important variable.

Both may be statistically significant – small amount of variation in Ca at a fixed level of pH may cause differences in species composition.

Arch – very rarely occurs in CCA. Detrended CCA generally should not be used except in special cases.

INFLUENCE OF NOISY ENVIRONMENTAL DATA ON CANONICAL CORRESPONDENCE ANALYSIS

McCune (1997) Ecology 78, 2617–2623 Simulated artificial data 10 x 10 grid. 40 species following Gaussian response model.

(1) (2) (3) (4) 2 2 10 2 + 10 environmental variables environmental variables X and Y co-ordinates with added noise TENXTEN NOISMOD (random number mean = 0, variance 17%) added to each cell random environmental variables environmental variables with added noise from NOIS1O NOISMOD random environmental variables from NOIS 10 (5) 99 NOISFULL – random environmental variables increases for axis 1 and 2.

NOISBOTH NOISFULL ‘Species-environment’ correlation increases as number of random variables Is in fact the correlation between the linear combination and WA site scores.

Poor criterion for evaluating success.

Not interpreted as measure of strength of relationship.

Monte Carlo permutation tests - NO STATISTICAL SIGNIFICANCE!

Dependence of the 'species-environment correlation,' the correlation between the LC and WA site scores, on a second matrix composed of from 1 to 99 random environmental variables. This correlation coefficient is inversely related to the degree of statistical constraint exerted by the environmental variables.

TEN x TEN NOISMOD NOISIO NOISFULL

TENXTEN NOISMOD NOISE 10 NOISBOTH NOISFULL

Monte Carlo tests

* * (99 env vars) ns * ns  1 0.77

0.65

0.20

0.65

0.12

 2 0.79

0.65

0.12

0.65

0.08

1 0.98

0.90

0.49

0.91

1.0

2 0.9

0.8

0.4

0.8

1.0

Linear combination site scores WA site scores best fit of species abundances to the environmental data best represent the assemblage structure Sensitive to noise True direct gradient analysis (multivariate regression) Aim to describe biological variation in relation to

LC scores

+ + +

WA scores

– – – environment Assemblage structure – + Which to use depends on one's aims and the nature of the data.

‘Species-environmental correlation’ better called ‘LC-WA’ correlation. Better measure of the strength of the relationship is the proportion of the variance in the species data that is explained by the environmental data. Evaluation should always be by a Monte Carlo permutation test.

LC OR WA SCORES?

M IKE P ALMER "Use LC scores, because they give the best fit with the environment and WA scores are a step from CCA towards CA." B RUCE M C C UNE "LC scores are excellent, if you have no error in constraining variables. Even with small error, LC scores can become poor, but WA scores can be good even in noisy data." LC scores are the default in CANODRAW.

Be aware of both - plot both to be sure.

J. Oksanen (2002)

DATA ORDERINGS

CCA DIAGRAM

TEN SETS OF DISTANCES TO REPRESENT, EMPHASIS ON 5, 8, AND 1 (FITTED ABUNDANCES OF SPECIES AND SITES)

Data-tables in an ecological study on species environmental relations. Primary data are the sub-table 1 of abundance values of species and the sub-tables 4 and 7 of values and class labels of quantitative and qualitative environmental variables (env. var), respectively. The primary data are input for canonical correspondence analysis (CCA). The other sub-tables contain derived (secondary) data, as the arrows indicate, named after the (dis)similarity coefficient they contain. The coefficients shown in the figure are optimal when species-environmental relations are unimodal. The CA ordination diagram represents these sub-tables, with emphasis on sub-tables 5 (weighted averages of species with respect to quantitative environmental variables), 8 (totals of species in classes of qualitative environmental variables) and 1 (with fitted, as opposed to observed, abundance values of species). The sub-tables 6, 9, and 10 contain correlations among quantitative environmental variables, means of the quantitative environmental variables in each of the classes of the qualitative variables and chi-square distances among the classes, respectively. (Chis-sq = Chi-square; Aver = Averages; Rel = Relative)

DEFAULT CCA PLOT

• Like CA biplot, but now a triplot: vectors for linear constraints.

• Classes as weighted averages or centroids.

• Most use LC scores: these are the fitted values.

• Popular to scale species relative to eigenvalues, but keep sites unscaled. Species-conditional plot.

• Sites do not display their real configuration, but their projections onto environmental vectors are the estimated values.

J. Oksanen (2002)

REF

Emphasis on

1 Species x sites 2 Species x species 3 Sites x sites

SCALING IN CCA

Hill scaling –1 SITES

Rel abundances – Turnover distances

Default scaling 2 SPECIES

Fitted abundances (rel) Chi-squared distances –

Quant env vars

4 Sites x env vars 3 5 Species x env vars – Values of env vars

Weighted averages Weighted averages

6 Env vars x env vars Effects 2 Correlations 2

Qualit env vars

7 Sites x env classes 4 Membership 1 8 Species x env classes Rel total abund 9 Env vars x env classes – 10 Env classes x env classes REF Turnover distances Membership 1

Rel total abund

Mean values of env vars –

fitted by least squares

REF 1 by centroid principle 2 change in site scores if env variable changes are one standard deviation 3 inter-set correlations 4 group centroids REF

REF Sub-tables (row numbers) that can be displayed by two differently scaled ordination diagrams in canonical correspondence analysis (CCA). Display is by the biplot rule unless noted otherwise. Hill's scaling (column 2) was the default in

CANOCO 2.1

, whereas the species-conditional biplot scaling (column 3) is the default in

CANOCO 3.1

and

. The weighted sum of squares of sites scores of an axis is equal to  /(1 and equal to 1 in scaling -1 and scaling 2, respectively. The weighted sum of squares of species scores of an axis is equal to  ) with 1/(1  )  its eigenvalue and equal to  in scaling -1 and scaling 2, respectively. If the scale unit is the same of both species and sites scores, then sites are weighted averages of species scores in scaling -1 and species are weighted averages of site scores in scaling 2. Table in italics are fitted by weighted least-squares (rel. = relative; env. = environmental; cl. = classes; - = interpretation unknown). REF Note that symmetric scaling (=3) has many optimal properties (Gabriel, 2002; ter Braak, personal communication) REF REF

Scaling 1 .

species x sites a 2 .

species x species -1: focus on sites Hill's scaling Rel. Abundances b,c Interpreta tion CENTROID UNKNOWN DISTANCE f 2: focus on species biplot scaling of CCA

Fitted rel. abund.

b  -square distances d Interpreta tion BIPLOT rule or CENTROID rule DISTANCE rule DISTANCE rule 3 .

sites x sites Quantitative env. vars 4 .

sites x env. vars g Turnover distances c,e 5 7 8 .

species x env. vars 6 .

env.vars x env. vars Qualitative env. vars .

sites x env. classes species x env. cls.

i 9 . env.vars x env. classes 10 .

env. classes x env. classes.

Weighted averages

Effects h Membership k Rel. total abund.

c,b Turnover distances c,e UNKNOWN BIPLOT ? BIPLOT CENTROID CENTROID UNKNOWN DISTANCE f Values of env.vars

Weighted averages

Correlations Membership k

Rel. total abund.

b Mean values of env. vars BIPLOT rule BIPLOT rule BIPLOT rule CENTROID rule CENTROID rule BIPLOT rule DISTANCE rule

REF REF a Site scores are linear combinations of the environmental variables. The adjective "fitted" must be deleted if site scores are proportional to the weighted average of species scores.

b The centroid principle can be applied also if sites and species scores are plotted in the same units, i in scaling -1, species that occur in a site lie around it, whereas in scaling 2, the species' distribution is centred at the species point.

c The biplot rule cannot be applied d In the definition of this coefficient, abundance must be replaced by fitted abundance values, because CCA is correspondence analysis of fitted abundance values e No explicit formula known f Chi-square distances, provided the eigenvalues of the axes are of the same magnitude g (y Environmental scores are (intra-set) correlations in scaling 2; more precisely, the coordinate of an arrow head on an axis (i.e. the score) is the weighted product-moment coefficient of the environmental variable with the axis, the weights being the abundance totals of the sites i+ ). The scores in scaling -1 are {  (1  )} ½ times those in scaling 2.

h Effect is defined as the change in site scores if the environmental variable changes one standard deviation in value (while neglecting the other variables).

i Environmental points are centroids of site points k Via centroid principle, not via biplot REF REF

INTERPRETATION OF CCA PLOTS

Centroid principle Distance principle Biplot principle (of relative abundances) Small eigenvalues, short (< 4SD) gradients – Biplot principle Large eigenvalues (> 0.40), long (> 4SD) gradients – Centroid and distance principles and some biplot principles Note that the centroid and distance principle may approximate biplot principle if gradients are short and eigenvalues small.

Differences are least important if  1  2

} Ordinal 4 classes 3 classes

CCA EXAMPLE

7 binary class variables Remove effect of seasonal variation Example data: quantitative and qualitative environmental variables (a) and qualitative covariables (b) recorded at 40 sites along two tributaries from the Hierden stream (sd: standard deviation, min: minimum, max: maximum). Aquatic macro-fauna data

Ranking environmental variables in importance by their marginal (left) and conditional (right) effects of the macrofauna in the example data-set, as obtained by forward selection. ( = fit = eigenvalue with variable j only; 

= additional fit = increase in eigenvalue; cum ( 

out by taking the month class variables as covariables.



) = cumulative total of eigenvalues 

; P = significance level of the effect, as obtained with a Monte Carlo permutation test under the null model with 199 random permutations; - additional variables tested; veg. = vegetation). Seasonal variation is partialled 2 3 Marginal effects (forward: step 1) j 1 Variable Shrubs (1/0)  1 0.25

(0.01) Source distance 0.22

EC 0.20

(0.01) (0.01) Conditional effects (forward: continued) j 1 2 3 Variable Shrubs (1/0) Source distance EC  a 0.25

0.19

(0.01) (0.01) (0.01) Cum ( 0.25

0.44

0.63

 a ) 4 Discharge 0.14

(0.03) 0.75

8 9 4 5 6 7 Discharge Total veg cover Shading Soil grain size Stream width High weedy veg 10 Cover bank veg U vs L stream 0.17

0.16

0.15

0.14

0.13

0.22

(0.01) (0.01) (0.01) (0.02) (0.05) (0.08) (0.11) (0.01) Cover emergent Cover bank veg Soil grain size 0.11

0.11

0.10

(0.10) (0. 12) (0.13) Each variable is the only env. var.

MARGINAL EFFECTS i.e. ignoring all other variables U vs L stream 0.09

(0.01)

EXTRA FIT

Change in eigenvalue if particular variable selected CONDITIONAL EFFECTS given other selected variables

Species-conditional triplot based on a canonical correspondence analysis of the example macro-invertebrate data displaying 13% of the inertia (=weighted variance) in the abundances and 69% of the variance in the weighted averages and class totals of species with respect to the environmental variables. The eigenvalue of axis 1 (horizontally) and axis 2 (vertically) are 0.35 and 0.17 respectively; the eigenvalue of the axis 3 (not displayed) is 0.13. Sites are labelled with stream code (U, L) and are ranked by distance from the source (rank number within the stream). Species (triangles) are weighted averages of site scores (circles). Quantitative environmental variables are indicated by arrows. The class variable shrub is indicated by the square points labelled Shrub and No shrub. The scale marks along the axes apply to the quantitative environmental variables; the species scores, site scores and class scores were multiplied by 0.4 to fit in the coordinate system. Only selected species are displayed which have N2>4 and a small N2-adjusted root mean square tolerance for the first two axes. The species names are abbreviated to the part in italics as follows Ceratopogonidae, Dendrocoelum lacteum, Dryops luridus, Erpobdella testacea, Glossiphonia complanata, Haliplus lineatocollis, Helodidae, Micropsectra atrofasciata, Micropsectra fusca, Micropterna sequax, Prodiamesa olivacea, Stictochironomus sp.

CANONICAL CORRESPONDENCE ANALYSIS (CCA) - A SUMMARY

• Unconstrained CA gives • Species ordination which is derived from site ordination • Site ordination which is derived from species ordination • Fitted vectors for environmental variables (indirect gradient analysis) • Constrained CA (Canonical CA) gives a direct gradient analysis • Species ordination which is derived from site ordination • Site scores which are linear combinations of environmental variables (LC scores) • Site ordination which is derived from species ordination (WA scores) so that species-environment correlation is maximised with the LC scores • Vectors of environmental variables that define the linear combination scores for sites

Gradient length estimation CA Indirectly CCA Directly Outline of ordination techniques present ed here. DCA (detrended correspondence analysis) was applied for the determina tion of the length of the gradient (LG). LG is important for choosing between ordination based on a linear or on an unimodal response model. Correspond ence analysis (CA) is not considered any further because in “microcosm experi ment discussed here LG < or = 1.5 SD units. LG < 3 SD units are considered to be typical in experimental ecotoxicology. In cases where LG < 3, ordination based on linear response models is considered to be most appropriate. PCA (principal component analysis) visualizes variation in species data in relation to best fitting theoretical variables. Environmental variables explaining this visualised variation are deduced afterwards, hence, indirectly. RDA ( redundancy analysis) visualises variation in species data directly in relation to quantified environ mental variables. Before analysis, covariables may be introduced in RDA to compensate for systematic differences in experimental units. After RDA, a permutation test can be used to examine the significance of effects.

REDUNDANCY ANALYSIS – CONSTRAINED PCA

Short (< 2SD) compositional gradients Linear or monotonic responses Reduced-rank regression PCA of y with respect to x Two-block mode C PLS PCA of instrumental variables Rao (1964) PCA RDA best hypothetical latent variable is the one that gives the smallest total residual sum of squares selects linear combination of environmental variables that gives smallest total residual sum of squares ter Braak (1994) Ecoscience 1, 127–140 Canonical community ordination Part I: Basic theory and linear methods

RDA ordination diagram of the Dune Meadow Data with environmental variables represen ted as arrows. The scale of the diagram is: 1 unit in the plot corresponds to 1 unit for the sites, to 0.067 units for the species and to 0.4 units for the environmental variables.

Redundancy analysis: canonical coefficients (100 x c) and intra-set correlations (100 x r) of environmental variables with the first two axes of RDA for the Dune Meadow Data. The environmental variables were standardized first to make the canonical coefficients of different environmental variables comparable. The class SF of the nominal variable “type of management” was used as reference class in the analysis. Variable Coefficients Axis1 Axis2 A1 Moisture Use -1 15 5 -5 9 -6 Manure SF BF HF NM -8 -10 -10 -4 16 0 -2 -13 Correlations Axis1 Axis2 54

15 -6 12 29 -26 25 -48 -40 51

86 76

-11 13

-79

PCA and RDA comparisons

Important to do the check that the environmental variables relate to the major gradients in composition detected by the PCA.

PCA RDA PCA RDA % %

Axis 1 Axis 2

29 26 Correlation 0.90

Correlation 0.95

21 17 0.82

0.89

BIPLOT INTERPRETATION

Cosine of angle  correlation Long arrows of species and environmental variables most important species unconstrained Goodness of fit  1 +  2 sum of eigenvalues constrained Euclidean distance biplot Covariance (correlation) biplot RDA covariance or correlation matrix of species RDA – constrained form of multiple regression Uses 2 (q + m) + m parameters (

env variables,

Multiple regression m (q + 1) e.g.

40 species 10 envir variables RDA 140 parameters MR 440 parameters RDA is thus reduced rank regression (RR) fitted species species)

Primary and secondary data tables in a typical community ecological study of species environment relations. Indirect methods of ordination use the tables under (a). Direct methods also use the tables under (b). The primary data are the table of abundance values and the tables of values and class labels of quantitative and qualitative environmental variables (env. var), respectively. The secondary tables are named after the (dis)similarity coefficients they contain. The appropriate coefficients must be chosen by the ecologist. The coefficients shown in the figure are optimal when species-environment relations are linear.

Tables that can be displayed by two differently scaled biplots in principal components analysis (a) and redundancy analysis (b).

The sum of squares of site scores of an axis is equal to its eigenvalue in scaling 1, and equal to 1 in scaling 2. The sum of squares of species scores of an axis is equal to 1 in scaling 1 and equal to its eigenvalue in scaling 2. Tables in bold are fitted by (weighted) least-squares.

Biplot scaling 1: focus on sites distance biplot

(a) principal components analysis species x sites sites species (b) redundancy analysis species x sites b sites b species

Quantitative env. vars.:

species x env. vars.

sites x env. vars.

d d env. vars.

Qualitative env. vars:

species x env. classes f sites x env. classes f env. classes f env. vars. x env. classes g

abundances

Euclidean distances -

fitted abundances

Euclidean distances c -

correlations

effects e

means

Euclidean distances g

2: focus on species correlation biplot abundances

correlations fitted abundances

correlations

a,c

correlations

values of env. vars correlations

means

means a

REF a Automatic if abundance is standardised by species. If abundance is only centred by species, a post-hoc rescaling of the site scores is needed so as to account for the differences in variance amongst species.

REF b Site scores are a linear combination of the environment variables instead of being a weighted sum of species abundances.

c In the definition of this coefficient, abundance must be replaced by the fitted abundance.

d Environmental scores are intraset correlations in scaling 2 and  s ½ with  s environmental variables.

times those in scaling 1 the eigenvalue of axis . In CANOCO, the scores are termed biplot scores for e Effect of the environmental variable on the ordination scores, while neglecting the other environmental variables; length of arrow is the effect size, i.e. the variance explained by the variable.

f Environmental classes are centroids of site points belonging to the class.

g membership via centroid principle, not via the biplot rules.

REF REF

Correlation biplot based on a redundancy analysis of the Dune Meadow Data displaying 43% of the variance in the abundances and 71% of the variances in the fitted abun dances. Quantitative environ ment variables are indicated by arrows. The qualitative variable Management type is indicated by the square points labelled SF, BF, HF, and NM. The displayed species are selected on the basis that more than 30% of their variance is accounted for by the diagram. Eigenvalues of the first three axes are 0.26, 0.17,and 0.07; the sum of all canonical eigenvalues is 0.61.

The scale marks along the axes apply to the species and quantitative environmental variables; the site scores and class scores were multiplied by 0.46 to fit in the coordinate system. The abbreviations are given in Jongman et al. (1987).The rule for interpreting a biplot (projection on an imaginary axis) is illustrated for the species Pla lan and sites 11 and 12.

PROPOSED NEW SCALING FOR CCA AND RDA

Gabriel, K.R. (2002) Biometrika 89, 423-436 Symmetric scaling (3) of biplots preserves the optimal fit to the species data table and preserves the (proportional) fit of at least 95% of the inter-species correlations/distances and inter-sample distances. It is a very good compromise.

Only recommended (ter Braak, pers. comm.) to deviate from symmetric scaling if the focus of study is strongly on either species (scaling 2) or on samples (scaling 1).

Data table unaffected by scaling: Species x sites Species data (PCA) Fitted species data (RDA) Relative species data (CA) Fitted relative species data (CCA) Species x environmental variables Correlations of species (RDA) Species x environmental classes Optima (WA) of species (CCA) Mean abundances of species (RDA) Relative abundances of species across classes (CCA)

Data tables with 95% preservation of proportional fit: Species x species Correlations (PCA, RDA) Chi-square distances (CA, CCA) Sites x sites Env. classes x env. classes Euclidean distances (PCA, RDA) Chi-square distances (CA, CCA) Euclidean distances (RDA) Chi-square distances (CCA) Env. variables x env. variables Sites x env. variables Correlations (RDA, CCA) Values (RDA, CCA) Sites x env. classes Env. variables x env. classes Means (RDA, CCA) Mean values of env. variables (RDA, CCA)

ALTERNATIVES TO ENVIRONMENTAL VECTORS IN CCA AND RDA

• Fitted vectors natural in constrained ordination, since these have linear constraints.

• Distant sites are different, but may be different in various ways: environmental variables may have a non-linear relation to ordination.

Contours Bubble plots GAM J. Oksanen (2002)

STATISTICAL TESTING OF CONSTRAINED ORDINATION RESULTS

Statistical significance of species-environmental relationships. Monte Carlo permutation tests.

Randomly permute the environmental data, relate to species data ‘random data set’. Calculate eigenvalue and sum of all canonical eigenvalues (trace). Repeat many times (99).

If species react to the environmental variables, observed test statistic (  1 or trace) for observed data should be larger than most (e.g. 95%) of test statistics calculated from random data. If observed value is in top 5% highest values, conclude species are significantly related to the environmental variables.

STATISTICAL SIGNIFICANCE OF CONSTRAINING VARIABLES

• CCA or RDA maximise correlation with constraining variables and eigenvalues.

• Permutation tests can be used to assess statistical significance: - Permute rows of environmental data.

- Repeat CCA or RDA with permuted data many times.

- If observed  higher than (most) permutations, it is regarded as statistically significant.

J. Oksanen (2002)

e.g.

PARTIAL CONSTRAINED ORDINATIONS (Partial CCA, RDA, etc)

pollution effects seasonal effects  COVARIABLES Eliminate (partial out) effect of covariables. Relate residual variation to pollution variables.

Replace environmental variables by their residuals obtained by regressing each pollution variable on the covariables.

Analysis is conditioned on specified variables or covariables.

These conditioning variables may typically be 'random' or background variables, and their effect is removed from the CCA or RDA based on the 'fixed' or interesting variables.

PARTIAL CCA

Natural variation due to sampling season and due to gradient from fresh to brackish water partialled out by partial CCA.

Variation due to pollution could now be assumed. Ordination diagram of a partial canonical correspond-ence analysis of diatom species (A) in dykes with as explanatory variables 24 variables-of-interest (arrows) and 2 covariables (chloride concentration and season). The diagram is symmetrically scaled [23] and shows selected species and standardized variables and, instead of individual dykes, centroids (•) of dyke clusters. The variables-of-interest shown are: BOD = biological oxygen demand, Ca = calcium, Fe = ferrous compounds, N = Kjeldahl-nitrogen, O 2 = oxygen, P = ortho-phosphate, Si= silicium compunds, WIDTH = dyke width, and soil types (CLAY, PEAT). All variables except BOD, WIDTH, CLAY and PEAT were transformed to logarithms because of their skew distribution. The diatoms shown are: Ach hun = Achnanthes hungarica,

Ach min = A. minutissima, Aph cas= Amphora castellata Giffen, Aph lyb = A. lybica, Aph ven = A. veneta, Coc pla = Cocconeis placentulata, Eun lun = Eunotia lunaris, Eun pec = E. pectinalis, Gei oli = Gomphoneis olivaceum, Gom par = Gomphonema parvulum, Mel jur = Melosira jürgensii, Nav acc = Navicula accomoda, Nav cus = N. cuspidata, Nav dis = N. diserta, Nav exi = N. exilis, Nav gre = N. gregaria, Nav per = N. permitis, Nav sem = N. seminulum, Nav sub= N. subminuscula,Nit amp = Nitzschia amphibia, Nit bre = N. bremensis v. brunsvigensis, Nit dis = N. dissipata, Nit pal = N. palea, Rho cur = Rhoicosphenia curvata.

(Adapted from H. Smit, in prep)

PARTIAL ORDINATION ANALYSIS (Partial PCA, CA, DCA)

1) 2) 3) There can be many causes of variation in ecological or other data. Not all are of major interest. In partial ordination, can ‘factor out’ influence from causes not of primary interest. Directly analogous to partial correlation or partial regression. Can have partial ordination (indirect gradient analysis) and partial constrained ordination (direct gradient analysis). Variables to be factored out are ‘COVARIABLES’ or ‘CONCOMITANT VARIABLES’. Examples are: 4) 5) 6) 7) Differences between observers.

Time of observation.

Between-plot variation when interest is temporal trends within repeatedly sampled plots.

Uninteresting gradients, e.g. elevation when interest is on grazing effects.

Temporal or spatial dependence, e.g. stratigraphical depth, transect position, x and y co-ordinates. Help remove autocorrelation and make objects more independent.

Collecting habitat – outflow, shore, lake centre.

Everything – partial out effects of all factors to see residual variation in data. Given ecological knowledge of sites and/or species, can try to interpret residual variation. May indicate environmental variables not measured, may be largely random, etc.

PARTIAL ORDINATIONS

e.g. partial out the effects of some covariables prior to indirect gradient analysis within-plot change between-plot differences PRIMARY INTEREST NOT OF INTEREST Partial plot identity, ordination of residual variation, i.e. within-plot change.

e.g. Swaine & Greig-Smith (1980) Bakker et al. (1990) J Ecol 68, 33–41 J Plankton Research 12, 947–972

COVARIABLES IN CCA AND RDA

Background variables or 'covariables' Partial CCA Partial RDA Vegetation Environmental variables or 'constraints' CCA RDA Vegetation (residual) CA DCA PCA Vegetation (residual) "Nuisance" variables or other background factors can be removed before studying interesting factors. Partial CCA or partial RDA.

Permutation tests are for environmental variables only.

Residual variation can be analysed at any level. Can partition the variance.

Final residual shows what you cannot explain with available environmental variables.

Interpretation of final residual based on other correlates and/or ecological knowledge.

PARTITIONING VARIATION

ANOVA Two-way ANOVA   total SS = regression SS + residual SS between group (factor 1) + between treatments (factor 2) + interactions + error component Ecology 73, 1045–1055 Borcard et al. (1992) Variance or variation decomposition into 4 components Important to consider groups of environmental variables relevant at same level of ecological relevance (e.g. micro-scale, species-level, assemblage level, etc.).

Variation = variance in RDA Variation = inertia in CCA = chi-square statistic of data divided by the data’s total = sum of all eigenvalues of CA

Total inertia = total variance 1.164

Sum canonical eigenvalues = 0.663

Explained variance  Unexplained variance = T – E What of explained variance component?

57% 57% 43% Soil variables (pH, Ca, LOI) Land-use variables (e.g. grazing, mowing) Not independent Do CCA/RDA using 1) Soil variables only 2) Land-use variables only 3) Partial analysis Soil 4) Partial analysis Land-use a) Soil variation independent of land-use (3) b) Land-use structured (covarying) soil variation (1–3) c) Land-use independent of soil (4)   canonical eigenvalues 0.521

canonical eigenvalues 0.503

Land-use covariables Soil covariables 0.160

0.361

0.142

0.160

0.142

13.7% 31% 12.2% Total explained variance d) Unexplained 56.9% 43.1%

a b c d

unique covariance unique unexplained

CANOCO

VARIATION PARTITIONING OR DECOMPOSITION WITH 3 OR MORE SETS OF PREDICTOR (EXPLANATORY) VARIABLES

Qinghong & Bråkenheim, (1995) Water, Air and Soil Pollution 85, 1587–1592 Three sets of predictors – Series of RDA and partial RDA Climate (C), Geography (G) and Deposition of Pollutants (D)

Predictors

G+C+D D G+C G+C -

Covariables

G+C D D Joint effect D  G+C=0.784-0.132=0.679-0.027=0.652

C D+G G+D G+D C C Joint effect C  D+G=0.737-0.106=0.706-0.074=0.631

Sum of canonical

 0.811

0.027

0.811

0.784

0.132

0.679

0.811

0.106

0.706

0.074

0.737

0.812

0.811

Predictors

G D+C D+G -

Covariables

D+C G G Joint effect G  D+C=0.777-0.228=0.538-0.034=0.549

Sum of canonical

 0.034

0.777

0.228

0.538

0.811

All predictors Pure deposition Pure climate Pure geography

Canonical eigenvalues

0.811

0.027

0.106

0.034

PD PC PG Joint G + C Joint G + D 0.132

0.074

Joint D + C 0.228

Unexplained variance 1 – 0.811 = 0.189

C D

PD PC CD CDG DG CG PG

Covariance terms CD DG CG CDG

CD + DG + CDG = CD + CG + CDG = DG + CG + CDG = 0.652

0.631

0.549

PD + PC + CD = 0.027 + 0.106 + CD = 0.777 – 0.549 = 0.228

PD PC (D+C) (DG + CG + CDG) PD + PG + DG = 0.027 + 0.034 + DG = 0.706 – 0.631 = 0.074

PD PG (G+D) (CD + CG + CDG) PC + PG + CG = 0.106 + 0.034 + CG = 0.784 – 0.652 = 0.132

PC PG (G+C) (CD + DG + CDG)  CD = 0.095

DG = 0.013

CG = –0.008

CDG = 0.652 – 0.013 – 0.095

= 0.631 – (–0.008) – 0.095

= 0.544

= 0.054 – (–0.008) – 0.013

= 0.544

Total explained variance 0.811 consists of: Common climate + deposition Common deposition + geography 0.095

0.013

Common climate + geography 0.008

Common climate + geography + deposition 0.544

Unique climate PC 0.106

Unique geography PG 0.034

Unique deposition PD 0.027

Unexplained variance 0.189

See also Qinghong Liu (1997) – Environmetrics 8, 75–85 Anderson & Gribble (1998) – Australian J. Ecology 23, 158-167 Total variation: 1) random variation 2) unique variation from a specific predictor variable or set of predictor variables 3) common variation contributed by all predictor variables considered together and in all possible combinations Usually only interpretable with 2 or 3 'subsets' of predictors.

In CCA and RDA, the constraints are linear. If levels of the environmental variables are not uncorrelated (orthogonal), may find negative 'components of variation'.

REF REF

'NEGATIVE' VARIANCES

In variance partitioning, the groups of predictor variables used should be non-linearly independent for unbiased partitioning or decomposition.

If the groups of variables have polynomial dependencies, some of the variance components may be negative. Negative variances are, in theory, impossible.

High-order dependencies commonly arise with high numbers of variables and high number of groups of variables.

Beware of inter-relationships between predictor variables and between groups of predictors. Problem common to all regression based techniques, including (partial) CCA or RDA.

REF Careful model selection (minimal adequate model) is essential for many purposes, including variance partitioning.

REF

UNBIASED ESTIMATES OF VARIANCE COMPONENTS

Peres-Neto et al. 2006 Ecology 87: 2614-2615 Legendre 2007 Journal of Plant Ecology (in press) Variation in response variable y or Y = Variation explained by X [a] [b] [c] Variation explained by Z Unexplained (residual) variation = [d] Rectangle = 100% of the variation in

(response variable(s)) Fraction

is the intersection (not interaction) or covariance of the amounts of variation explained by linear model of

and

Partial linear regression

y ~X|Z – partial linear regression of response variable

on predictor variables

(

variables) whilst controlling for the linear effect of

containing

covariables

Partial linear regression y = x + covariables z Partial R 2 y~X|Z R 2 y~X|Z = SS (fitted values of y~X|Z) (SS (fitted values) + SS (residuals)) = [a] ([a] + [d]) F statistic used to test significance of partial R 2 takes into account number of covariables

; in ordinary multiple regression

= 0 F = (R 2 partial m) ((1 – R 2 partial ) (n – 1 – m – q)) F = ([a] m) ([d] (n – 1 – m – q))

= number of predictor variables , n = number of objects Also compute y~Z|X as well as y~X|Z

Unadjusted and adjusted coefficients of determination R 2 R 2 = Regression SS Total SS = 1 – Residual SS Total SS

(A)

R 2 adj = 1 – Residual mean square = 1 – (1 – R 2 ) Total df Total mean square Residual df

(B)

R 2 adj takes into account numbers of degrees of freedom associated with numerator and denominator in equation A In ordinary multiple regression total df = (n – 1) residuals df = (n – m - 1) where

and = number of observations and

= number of explanatory variables

Partial canonical analysis

In RDA, canonical total SS.

R 2 is ratio of the sum of each response variable’s regression (or fitted values) SS to the sum of all response variables’ Significance of F tested by permutation with ( m x p ) and p (n – m – 1) degrees of freedom where

= number of response variables Can calculate R 2 adjusted (equation B) to produce unbiased estimates of the real contributions of variables in

to the explanation of the response matrix

2 ( Y |

) adj  1 



 1

 1  1 

  1  1  1

n m

 1 ( 1 

) where

= number of predictor variables in matrix

Variation partitioning of

explanatory variables with respect to

and

, two sets of Canonical analysis Y~X Y~W Y~(X,W) R 2 R 2 2 of Y~X of Y~W R 2 adj * and fractions [a+b]=R a 2 of Y~X [b+c]=R a 2 of Y~W R 2 of Y~(X,W) [a+b+c]=R a 2 [a]=[a+b]-[b] Can test significance Yes Yes Yes Yes [b]=[a+b]+[b+c]-[a+b-c] No [c]=[b+c]-[b] Residuals = [d]=1-[a+b+c] Yes No * Calculated as 1-(1-R 2 )(Total df/Residual df) = 1 – Residual mean square/Total mean square varpart in R vegan

Significance of fractions a and c tested by ‘permutation of the residuals’ option in CANOCO Fraction [a]

  [

]

  [

] (

 1 



)  Fraction [c]

  [

]

 ([

] (

 1 



)) where

m q

= number of explanatory variables in predictor set

Significance of fractions [b] (covariance) and [d] cannot be assessed statistically (residuals)

Oribatid mite data: 35 species, 70 samples, 5 habitat variables, and spatial descriptors Fractions Unadjusted Adjusted Difference Habitat [a] Shared habitat and space [b] 0.136 * 0.391

0.091 * 0.345

ns Spatial [c] Residual [d] 0.232 * 0.241

0.101 *

0.462

ns varpart in R vegan

Representation of components Fraction [a] dot maps of ‘bubble’ lots Fraction [b] Fraction [c] Fraction [d] interpolation maps interpolation maps (e.g. kriging) maps either dot or interpolated Can be expanded to 3 or 4 tables of explanatory variables. Gets very complex!

Applicable to CCA but algebra even more complicated.

ENVIRONMENTAL CONSTRAINTS AND CURVATURE IN ORDINATIONS

• Curvature often cured because axes are forced to be linear combination of environmental variables (constraints).

• High number of constraints = no constraint.

• Absolute limit: number of constraints = min (M, N) - 1, but release from the constraints can begin much earlier.

• Reduce environmental variables so that only the important remain: heuristic value better than statistical criteria.

• Reduces multicollinearity as well.

J. Oksanen (2002)

Classification of gradient analysis techniques by type of problem, response model and method of estimation

Method of estimation

Type of problem

Linear Least Squares Maximum Likelihood Unimodal Weighted Averaging

Regression Calibration Ordination Constrained ordination Partial ordination 2 1 Multiple regression Gaussian regression Linear calibration 'inverse regression' Principal components analysis (PCA) Redundancy analysis (RDA) 4 Partial component analysis Gaussian calibration Gaussian ordination Gaussian canonical ordination Partial Gaussian ordination Weighted averaging of site scores Weighted averaging of species scores (WA) Correspondence analysis (CA); detrended correspondence analysis (DCA) Canonical correspondence analysis (CCA); detrended CCA Partial correspondence analysis; partial DCA Partial constrained ordination 3 Partial redundancy analysis Partial Gaussian canonical ordination Partial canonical correspondence analysis; partial detrended CCA 1 = constrained multivariate regression 2 = ordination after regression on covariables 3 = constrained ordination after regression on covariables = constrained partial multivariate regression 4 = 'reduced rank regression' = “PCA of y with respect to x”

ENVIRONMENTAL VARIABLES IN CONSTRAINED ORDINATIONS

1) Choice can greatly influence the results. Fewer the environmental variables, the more constrained the ordination is.

2) Possible to have one only – can evaluate its explanatory power.

3) Can always remove superfluous variables if they are confusing or difficult to interpret. Can often remove large number without any marked effect. Remember post-hoc removal of variables is not valid in a hypothesis-testing analysis.

4) Linear combinations – environmental variables cannot be linear combinations of other variables. If a variable is a linear combination of other variables, singular matrix results, leads to analogous process of dividing by zero.

Examples: total cations, Ca, Mg, Na, K, etc. Delete total cations - % clay, % silt, % sand dummy variables (granite or limestone or basalt) 5) Transformation of environmental data – how do we scale environmental variables in such a way that vegetation ‘perceives’ the environment? Need educated guesses.

Log transformation usually sensible – 1 unit difference in N or P is probably more important at low concentrations than at high concentrations.

As statistical significance in CANOCO is assessed by randomisation tests, no need to transform data to fulfil statistical assumptions.

Transformations useful to dampen influence of outliers. Environmental data automatically standardised in RDA and CCA.

6) Dummy variables – factors such as bedrock type, land-use history, management, etc, usually described by categorical or class variables. 1 if belongs to class, 0 if it does not. For every categorical variable with K categories, only need K – 1 dummy variables e.g.

Plot 1 2 3 4 5 6

Granite

1 0 1 0 0 0

Limestone Basalt Gabbro

0 1 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 7) Circular data – some variables are circular (e.g. aspect) and large values are very close to small values. Aspect – transform to trigonometric functions.

northness = cosine (aspect) eastness = sine (aspect) Northness will be near 1 if aspect is generally northward and –1 if southward. Close to 0 if west or east.

Alternatively for aspect southness = 180 - |aspect - 180| (S = 180, N =0) westness = |180 - |aspect - 270|| (W = 180, E = 0) Day of year – usually not a problem unless dealing with sampling over whole year. Can create ‘winterness’ and ‘springness’ variables as for aspect.

8) Vegetation-derived variables – maximum height, total biomass, total cover, light penetration, % open ground can all be ‘environmental’ variables. Such variables SHOULD NOT BE USED in hypothesis testing, as danger of circular reasoning.

9) Interaction terms – e.g. elevation * precipitation. Easy to implement, difficult to interpret. If elevation and precipitation interact to influence species composition, easy to make term but the ecological meaning of where in environmental space the stands or species are is unclear. Huge number of possibilities N variables  ½ N (N – 1) possible interactions. 5 variables  10 interactions.

AVOID quadratic terms [e.g. pH * pH (pH 2 ) (cf. multiple regression and polynomial terms)]. Can create an ARCH effect or warpage of ordination space.

Try to avoid interaction terms except in clearly defined hypothesis-testing studies where the null hypothesis is that ‘variables c and d do not interact together to influence the species composition’.

For interaction to be significant, eigenvalue 1 of the analysis with product term should be considerably greater than  1 when there is no product term and the t-value associated with the product term should be greater than 2 in absolute value.

Avoid product variables to avoid ‘data dredging’.

SELECTING ENVIRONMENTAL VARIABLES IN CON STRAINED ORDINATION ANALYSIS (e.g. CCA, RDA)

1) The fewer the environmental variables, stronger the constraints are.

2) With q  (number of samples – 1) environmental variables, the analysis is unconstrained.

3) Small numbers of environmental variables may remove any arch effect.

4) Want to try to find MINIMAL ADEQUATE SET of environmental variables that explain the species data about as well as the FULL SET.

5) Automatic selection (e.g. forward selection) can be dangerous: a) Several sets can be almost equally good. Automatic selection finds one but may not be the best.

b) Selection order may change the result and important variables may not be selected.

c) Small changes in the data can change the selected variables. Difficult to draw reliable conclusions about relative importance of variables. Omission of a variable does not mean it is not ecologically important.

6) If you are lucky, there may only be one minimal adequate model but do not assume that there is only one such model.

7) How do we go about finding a minimal adequate model or set of environmental variables?

1) Start with all explanatory variables in the analysis, FULL MODEL. Consider sum of canonical eigenvalues (amount of explained variance), eigenvalues and species environmental correlations.

2) Try to simplify full model by deleting variables but not reducing the model performance. May be impossible to remove variables without some loss of information.

Deletion criteria:

a) Deletion on external criteria – variables not relevant.

b) Deletion on correlation structure – variables may be highly correlated (e.g. pH, Ca, Mg, CEC). Any one could be used as a proxy for the others. Best to choose the one that is likely to be the most direct cause of vegetation response.

Can do a PCA of environmental variables to explore correlation structure of variables.

c) Interpretability – variables with short arrows.

d) Non-significant – delete any that are non-significant (p > 0.05) in analysis with one environmental variable only in CCA or RDA.

e) Ecological importance f) Stepwise analysis – forward selection, add one variable at a time until no other variables ‘significantly’ explain residual variation in species data. 3) Final selection must be based on ecological and statistical criteria. The purpose of numerical data analysis is 'INSIGHT', not complex statistics!

WHAT IS DONE?

1) CCA (or RDA) is performed on each variable separately and marginal effects are listed in order.

2) Select the variable with largest marginal effect (= eigenvalue) and test its statistical significance by unrestricted Monte Carlo permutation tests and 999 permutations. Accept if p < 0.05.

3) This variable is now used as a covariable and the variables are now listed in order of their conditional effects (i.e. variance explained when allowing for effects of variable one selected). Evaluate its statistical significance and apply Bonferroni-type correction for simultaneous multiple tests, namely  1 =  /t where t =number of tests.

For  = 0.05 t = 1,  1 = 0.05

t = 2,  1 = 0.025

t = 3.  1 = 0.0166

t = 4,  1 = 0.0125

With 999 permutations (i.e. p of 0.001 can be evaluated), becomes very slow. Required if you are to properly evaluate the Monte Carlo permutation probabilities.

These tests do not control for overall Type I error. In practical terms this means that too many variables will be judged ‘significant’.

Alternatively, can stop when the increase in fit when including a variable is less than 1.0% (EXTRA FIT).

MINIMAL ADEQUATE MODEL IN CCA

13 environmental variables 3 environmental variables J. Oksanen (2002)

OTHER PROBLEMS

1) Selection of categorical variables coded as dummy variables. Suppose there are 5 bedrock types but only ‘granite’ is selected by forward selection. Should you select the other variables as well? If you consider the different bedrock types to be independent, the answer is NO. If you consider there to be one categorical variable (bedrock) with five states, the answer is YES.

2) Last two remaining variables within a category will always have identical fit because they contain identical information (if it is not z then it must be y). Does not matter which you choose. Select the commoner category.

3) No guarantee that forward selection (or any other stepwise procedure) will result in ‘best’ set of variables. Only way is perform constrained ordinations for every conceivable combination of variables. Currently impossible with current technology.

4) Accept that minimal adequate model is one possible solution only.

5) For exploratory, descriptive studies, do not be reluctant to use a priori ecological knowledge.

VARIANCE INFLATION FACTORS (VIF)

Variance of estimated regression (= canonical) coefficients (c

) are proportional to their VIF.

number of predictors Var  VIF   residual variance  



 1  number of samples VIF is related to the (partial) multiple correlation coefficient

R j

the other environmental variables. VIF  1  1

between variable

and Not unique, e.g. pH and Ca (and other variables).

If VIF > 20, that variable is almost perfectly correlated with other variables and has no unique contribution to the regression equation. Regression (= canonical) coefficient unstable, not worth considering.

pH Ca VIF

123.8

2.6

VIF

6.4

– Useful for finding minimal set of variables.

8.1

29.3

'AIC' FOR MODEL SELECTION IN CCA AND RDA

AIC Jari Oksanen (2004) VEGAN 1.7-6 R deviance.cca

deviance.rda

Find statistics in CCA and RDA that resemble deviance and assess an AIC like statistic as in regression model building.

Deviance of CCA = chi-square of the residual data matrix after fitting the constraints.

Deviance of RDA = average residual variance per species.

Can be used to help select between possible models in CCA or RDA.

index of fit that takes account of the parsimony of the model by penalising for the number of parameters.

smaller the values, better the fit.

here equals the residual deviance + 2x number of regression (canonical) coefficients fitted.

STAGES IN 'AIC' MODEL SELECTION IN CCA AND RDA

1. Define a null model into which variables are sequentially added in order of their statistical importance. Null model is unconstrained PCA or CA.

2. Now do stepping by either a forward selection of a backward elimination of the predictor variables. Need to define an upper and lower scope for the stepping to occur within.

Forward selection – lower scope = null model (no predictors) - upper scope = full model (all predictors included) Backward elimination – lower scope = full model - upper scope = null model

3. At each step, the effect of adding or deleting a variable is evaluated in terms of the AIC criterion. Low AIC values are to be preferred.

4. If a lower AIC can be achieved by adding or deleting a variable at a stage, then this predictor variable is added/deleted.

5. Useful to use both backward elimination and forward selection at each step. Start with full model, eliminate first variable, then the next, try to add either variable back into the model, and so on.

6. After the final model is derived (lowest AIC), can test this model to see if the effects of the constraining predictor variables are statistically significant. Use Monte Carlo permutation test under the reduced model.

'AIC' MODEL SELECTION

Godínez-Domínguez & Freire 2003 Marine Ecology Progress Series 253, 17-24 Rather than a statistical test of one null hypothesis, AIC provides a methodology for selecting an a priori set of alternative hypotheses.

1. Definition of set of a priori models 2. Statistical fit of models to data (e.g. CCA) 3. Selection of 'best' model – Akaike Information Criterion (AIC) AIC = 

2 log

  

2

where

= number of free parameters in the model

 ˆ = model maximum likelihood

From estimated residual sum of squares (RSS) in CCA

(



sum of all eigenvalue s



i h

  1

trace

)

where h = number of predictor variables in model ,

log

 

1 2

log(

)

where log = log e n = sample size, and 2 = RSS/n

To avoid bias in AIC due to links between sample size and number of parameters, corrected AIC is AIC c  AIC    2

(

K n



  1 ) 1   As in GLM, interested in differences in AIC between models  i = AIC ci – min AIC c

Data – 5 cruises (DEM-1 – DEM-5) - 8 models Godinez-Dominguez & Freire, 2003

CCA permutation tests 40 models 27 p < 0.05 (global test) 20 p < 0.05 (first CCA axis) AIC c approach to find most parsimonious model

Can determine not only the 'best' model but rank the different underlying hypotheses according to AIC criteria of parsimony. Spatial models 2 and 5, namely depth stratification but no difference between sheltered and exposed stations, are most appropriate for these data.

REF

CANONICAL CORRELATION ANALYSIS – CANCOR

REF Standard linear technique for relating two sets of variables.

Similar to RDA – assumes linear response model.

Selects canonical coefficients for species and environmental variables to MAXIMISE species – environmental correlation  canonical correlation   1 RDA species scores are simply weighted sums of site scores 

1  

b k y ki

 CANCOR species scores are b parameters estimated by multiple regression of site scores on species variables  number of species << number of sites In fact number of species + number of environmental variables must be smaller than number of sites i.e.

must be <

–

CANCOR biplot Differs from RDA also in error component.

van der Meer (1991) J. Expl. Mar. Biol. Ecol. 148, 105–120 REF

REF REF RDA uncorrelated independent errors with equal variances (least-squares technique).

CANCOR correlated normal errors (maximum likelihood technique). Is in realty, a GLM.

 residual correlations between errors are additional parameters in CANCOR. Many species. Cannot estimate them reliably from data from few sites.

Generalised variance minimised in CANCOR = product of eigenvalues of matrix ‘volume’ of hyperellipsoid.

Total variance minimised in RDA = sum of diagonal elements = sum of eigenvalues.

Linear transformation Non-linear transformation Linear transformation of species data of species data of environmental data

CA, PCA affects results CCA, RDA affects results CANCOR REF no effect    no effect no effect no effect REF

COMPARISON WITH CANONICAL CORRELATION ANALYSIS

McKecknie et al. (1975) Genetics 81, 571–594 Genetic variability in 21 colonies of butterfly Euphydryas editha in relation to 11 environmental variables in California and Oregon.

Gene frequencies of phosphoglucose-isomerase (Pgi) determined by electrophoresis.

What are relations between Pgi frequencies and environment?

Canonical correlation analysis

Can corr  1  0 .

77 0.88

 2  0 .

56 0.75

 3  0 .

17 0.41

Not significant Small sample size Leaving aside significance, first CC contrast between percentage of 0.4 mobility genes and other genes and precipitation and low max temperatures.

i.e. 0.4 mobility gene lacking in colonies from areas of high precipitation and low max temperatures.

Frequency of 1.0 mobility genes high in colonies at high altitudes and high precipitation and low temperatures.

Environmental variables and phosphoglucose isomerase Pg and Oregon. i gene frequencies for colonies of the butterfly Euphdryas editha in California Data source: McKechnie et al. (1975). Environmental variables have been rounded to integers for simplicity. The original data were for 21 colonies. For the present example, five colonies with small samples for the estimation of gene frequencies have been excluded so as to make all of the estimates about equally reliable.

Colonies of Euphydryas editha in California and Oregon

Pg i Pg i Pg i Pg m

Canonical correspondence analysis

 1 = 0.16

0.4

1 1.0 4 1.16 5 1.4

6 0.6

2  2 = 0.07 Pg i 0.8 3 opposite to latitude, i.e. southerly high max temp, high min temp ?

high ppt n

Ordination diagram based on canonical correspondence analysis of allele frequencies in the butterfly Euphdryas editha in the 21 colonies ( alleles (   ) with respect to 11 environmental variables (data from McKechnie et al. 1975). The ) from the three loci Pgm,Pgi and Hk, are numbered in order of increasing mobility class. Infrequent alleles are displayed by a dot. The environmental variables (arrows) are: altitude; latitude; annual precipitation (Annuprec), annual maximum and minimum temperature (Max-temp, Min-temp); highest and lowest temperature in posi-diapause (H Posdia,L-Posdia), the adult phase (H-Adul, L-Adul) and the pre diapause (H-predia, L-predia). The abbreviations of colony names follow Manly (1985). (  = 0 0.16,  = 0.07, scaling  = ½).

CONSTRAINED LINEAR ORDINATION (PCA FRAMEWORK)

• Canonical Correlation Analysis (CANCOR) • Continuous environmental variables and vegetation • Can be computed only if number of sites > number of species + number of env. vars +1 • Redundancy analysis (RDA) • As CANCOR but assumes that error variance constant for all plant species • Technically possible to estimate in vegetation data, unlike CANCOR • Canonical Variates Analysis (CVA) or Discriminant Analysis – see lecture 9 • Predict class membership using continuous variables • For instance, pre-determined vegetation type using vegetation data • LC score shows the centroids, Weighted sum scores show the dispersion and overlap

DISTANCE-BASED REDUNDANCY ANALYSIS

DISTPCOA

Pierre Legendre & Marti Anderson (1999) Ecol. Monogr. 69, 1-24.

RDA but with any distance coefficient

RDA Euclidean distance Absolute abundances Quantity dominated CCA chi-square metric Relative abundances Shape/composition dominated Does it matter?

Total biomass or cover and species composition Varying e.g. ridge  snow bed gradient Other dissimilarities Bray & Curtis Jaccard +/ Gower mixed data non-Euclidean non-Euclidean non-Euclidean semi-metric semi-metric semi-metric

Basic idea

Reduce sample x sample DC matrix (any DC) to principal co-ordinates (principal co-ordinates analysis, classical scaling, metric scaling – Torgerson, Gower) but with correction for negative eigenvalues to preserve distances.

PCoA – embeds the Euclidean part of DC matrix, rest are negative eigenvalues for which no real axes exist. These correspond to variation in distance matrix, which cannot be represented in Euclidean space. If only use positive eigenvalues, RDA gives a biased estimate of the fraction of variance in original data.

Correction for negative eigenvalues

 (

d ij

2  2

1 ) 0 .

for i



where c 1 is equal to absolute value of largest negative eigenvalue of matrix used in PCoA  1 

a ij ij

  

a ij

1 2 

d ij

a i



a j



D  1 Use all principal co-ordinate sample scores ( n - 1 or

, whichever is less) as RESPONSE (species) data in RDA. Use dummy variables for experimental design as predictors in

in RDA.

Now under framework of RDA and battery of permutation tests, can analyse structured experiments but WHOLE ASSEMBLAGE (cf. MANOVA but where

> N ).

Now can test null hypothesis (as in MANOVA) that assemblages from different treatments are no more different than would be expected due to random chance at a given level of probability. BUT unlike non-parametric tests (ANOSIM, Mantel tests), can test for interactions between factors in multivariate data but using any DC (not only Euclidean as in ANOVA / MANOVA). Using permutation tests means we do not have to worry about multivariate normality or homogeneity of covariance matrices within groups, or abundance of zero values as in ecological data.

DISTPCoA

www.umontreal.ca

Raw data (replicates x species) Distance matrix (Bray-Curtis, etc) Principal coordinate analysis (PCoA) Correction for negative eigenvalues Test of one factor in a single-factor model Matrix Y (replicates x principal coordinates) Redundancy analysis (RDA) F # statistic Matrix X (dummy variables for the factor) Test of F # by permutation Test of interaction term in multifactorial model Matrix Y (replicates x principal coordinates) Partial redundancy analysis (partial RDA) F # statistic Matrix X (dummy variables for the interaction) Matrix X (dummy C variables for the main effects) Test of F # by permutation under the full model

REF REF Correspondence between the various components of the univariate F-statistics and the multivariate RDA statistics in the one-factor case.

Univariate ANOVA Total sum of squares Treatment sum of squares = SS Tr Treatment degrees of freedom = df Tr Residual sum of squares = SS Res Residual degrees of freedom = df Res Treatment mean square = SS tr /df Tr = MS Tr Residual mean square SS Res /df Res = MS Res

 MS Tr MS Res Multivariate RDA statistic sum of all eigenvalues of Y trace = sum of all canonical eigenvalues of Y on X q rss = sum of all eigenvalues – trace n T – q – 1 trace/q rss/(n T – q - 1)

#  rss trace /(

n T

 /

q q

 1 ) Tr n T = treatment, q = number of linearly independent dummy variables, = number of replicates REF REF

Placing non-metric distances into Euclidean space first, then use ANOVA/MANOVA within RDA with permutation tests-.

Builds on ANOVA as form of multiple regression with orthogonal dummy variables as predictors.

MATCH between ANOVA statistics and RDA statistics.

PERMUTATION TESTS in RDA (also CCA) in CANOCO.

What is shuffled?

predictors responses Y = Z B + X C + E covariables random error B & C - unknown but fixed regression coefficients (Note Z = covariables and X = predictors here)

To test H 0 C = 0 (i.e. the effect of X on Y ) 1. Permute rows of Y 2. Permute rows of X (env. data)

CANOCO 2

3. Permute residuals E r of regression Y on REDUCED MODEL OR NULL MODEL Z (covariables)

CANOCO 3 & 4

4. Permute residuals E f FULL MODEL of regression Y on Z and X (covariables and predictors)

CANOCO 3 & 4

1 & 2 1 2 DESIGN-BASED PERMUTATIONS Wrong type I error, low power OK but no basis for testing interaction effects 3 & 4 MODEL-BASED PERMUTATIONS 3. Permute residuals of Y on Z (covariables) Default in CANOCO 3 & 4 Reduced model – maintains type I error in small data sets. Without covariables, gives exact Monte Carlo significance level. Retains structure in X and Z .

4. Permute residuals of Y slightly so.

on X and Z . Full model. Gives lower type II error, but only (If no covariables Y = XC + E , does not matter if samples in Y permutes X ) or X are permuted.

CANOCO

In DISTPCOA, do RDA with Y as principal co-ordinates scores, X defines dummy variables to code for interaction terms, and Z defines dummy variables for main effects (covariables) if interested in interactions. Can determine components of variation attributable to individual factors and interaction terms as in a linear model for multivariate data BUT using any DC that integrates both quantities and composition.

Can test the significance of individual terms and interaction terms for any complex multi-factorial experimental design. Cannot be applied to unbalanced data. If unbalanced because of missing or lost values, use missing data replacements (if other replicates).

Distance-based RDA offers special advantages to ecological researchers not shared by any other single multivariate method. These are: Shares characteristics with:

MAN OVA

1 The researcher has the flexibility to choose an appropriate dissimilarity measure, including those with semi-metric qualities, such as the Bray-Curtis measure 2 PCoA puts the information on dissimilarities among the replicates into a Euclidean framework which can then be assessed using linear models 3 A correction for negative eigenvalues in the PCoA, if needed, can be done such that probabilities obtained by a permutation test using the RDA F # statistic are unaffected (correction method 1) 4 By using the multiple regression approach to analysis of variance, with dummy variables coding for the experimental design, RDA can be used to determine the components of variation attributable to individual factors and interaction terms in a linear model for multivariate data 5 Multivariate test statistics for any term on a linear model can be calculated, with regard to analogous univariate expected mean squares 6 Statistical tests of multivariate hypothesis using RDA statistics are based on permutations, meaning that there is no assumption of multinormality of response variables in the analysis. Also, there are no restrictions to the number of variables that can be included in RDA 7 Permutations of residuals using the method of ter Braak (1992) allows the permutation test to be structured precisely to the hypothesis and the full linear model of the design under consideration 8 The significance of multivariate interaction terms can be tested * * *

RDA

* * * * *

ANO SIM

* *

MAN TEL

* *

DISTANCE-BASED MULTIVARIATE ANALYSIS FOR A LINEAR MODEL

McArdle, B.H. & Anderson, M.J. (2001) Ecology 82; 290-297 DISTLM, DISTLMforward - www.stat.auckland.ac.nz/~mja DISTLM - multivariate multiple regression of any symmetric distance matrix with or without forward selection of individual predictors or sets of predictors and associated permutation tests.

Y response variables ( n x m ) X predictor variables ( n x q ) (1/0 or continuous variables) Performs a non-parametric test of the multivariate null hypothesis of no relationship between Y and X on the basis of any distance measure of choice, using permutations of the observations. environmental variable) of interest.

X may contain the codes of an ANOVA model (design matrix) or it may contain one or more predictor variables (e.g., Like Legendre and Anderson's (1999) distance-based redundancy analysis but with no correction for negative eigenvalues. Shown theoretically that partitioning the variability in X according to a design matrix or model can be achieved directly from the distance matrix itself, even if the distance measure is semi-metric (e.g., Bray-Curtis distance).

CANONICAL ANALYSIS OF PRINCIPAL CO-ORDINATES

Anderson, M.J. & Willis, T.J. (2003) Ecology 84, 511-525 CAP – www.stat.auckland.ac.nz/~mja CAP - canonical analysis of principal co-ordinates based on any symmetric distance matrix including permutation tests.

Y response variables ( n x m ) X predictor variables ( n x q ) (1/0 or continuous variables) Performs canonical analysis of effects of statistical significance. X on Y on the basis of any distance measure of choice and uses permutations of the observations to assess If X contains 1/0 coding of an ANOVA model (design matrix), result is a generalised discriminant analysis. If X contains one or more quantitative predictor variables, result is a generalised canonical correlation analysis.

Output 1. Eigenvalues and eigenvectors from the PCOORD. Can use latter to plot an indirect ordination of data.

2. Canonical correlations and squared canonical correlations.

3. Canonical axis scores.

4. Correlations of original variables ( Y ) with canonical axes.

5. Diagnostics to help determine appropriate value for squares (in the case of quantitative variables in X Q

, number of eigenvectors. Select the lowest misclassification error (in the case of groups) or the minimum residual sum of ). Also Q

must not exceed

and is chosen so that the proportion of variance explained by the first Q

axes is more that 60% and less than 100% of the total variation in the original dissimilarity matrix.

6. In the case of groups, table of results for 'leave-one-out' classification of individual observations to groups.

7. Results of permutation test to test statistical significance of Q

eigenvalue).

axis model (trace and first 8. Scores to construct constrained ordination diagram to compare with unconstrained ordination diagram.

Very good at highlighting and testing for group differences (e.g. sampling times) as CAP finds axes that maximise separation between groups.

With quantitative predictors, CAP finds axes that maximise correlation with predictor variables.

'AIC' for model selection deviance-capscale Jari Oksanen VEGAN 1.7-6 R

Extensions by Jari Oksanen (capscale in Vegan R) 1. Axes are weighted by their corresponding eigenvalues so that the ordination distances are best approximations of the original dissimilarities.

2. Uses all axes with positive eigenvalues. Guarantees that the results are the best approximation of the original dissimilarities.

3. Adds species scores as weighted sums of the (residual) species data.

4. Negative eigenvalues are harmless and can be ignored. Often most sensible to use dissimilarity coefficients that do not have negative eigenvalues. Square-root transformation of the species data prior to calculating dissimilarities can drastically reduce the number of negative eigenvalues.

Note that CAP with Euclidean distance is identical to RDA in sample, species, and biplot scores (except for possible reversal of sign).

Possible uses of canonical analysis of principal co-ordinates 1. As in CCA or RDA with biological and environmental data.

2. Fit models to data with rare or unusual samples or species that may upset CCA.

3. Analyse many environmental variables in relation to external (e.g. geographical, geological, topographical) constraints with Monte Carlo permutation tests. In other words, do a multivariate linear regression but not have to worry about the data meeting the assumptions of least-squares estimation and models.

Examples: Willis & Anderson 2003 Marine Ecology Progress Series 257: 209-221 (cryptic reef fish assemblages) Edgar et al. 2003 Journal of Biogeography31: 1107-1124 (shallow reef fish and invertebrate assemblages)

SUMMARY OF CONSTRAINED ORDINATION METHODS

Methods of constrained ordination relating response variables, Y (species abundance variables) with predictor variables, X (such as quantitative environmental variables or qualitative variables that identify factors or groups as in ANOVA).

Name of methods (acronyms, synonyms) Redundancy Analysis (RDA) Canonical Correspondence Analysis (CCA) Distance measure preserved Euclidean distance Chi-square distance Canonical Correlation Analysis (CCorA, COR) Canonical Discriminant Analysis (CDA; Canonical Variate Analysis CVA; Discriminant Function Analysis, DFA) Canonical Analysis of Principal Coordinates (CAP; Generalized Discriminant Analysis) Mahalanobis distance Mahalanobis distance Any chosen distance or dissimilarity Relationship of ordination axes with original variables Linear with X, linear with Q distance measure) t ; unknown with Y (depends on Takes into account correlation structure Linear with X, linear with fitted values, Y = X(X'X) -1 X'Y Linear with X, approx unimodal with Y, linear with fitted values, Y* Linear with X, linear with Y Linear with X, linear with Y ... among variables in X, but not among variables in Y ... among variables in X, but not among variables in Y ... among variables in X, and among variables in Y ... among variables in X, and among variables in Y ... among variables in X, and among principal coordinates Q t

CRITERION FOR DRAWING ORDINATION AXES

RDA CCA CCorA CVA CAP • Finds axis of maximum correlation between Y and some linear combination of variables in X X (i.e., multivariate regression of Y , followed by PCA on fitted values, Y ^ ).

on • Same as RDA, but Y are transformed to Y* and weights (square roots of row sums) are used in multiple regression.

• Finds linear combination of variables in Y correlated with one another.

and X that are maximally • Finds axis that maximises differences among group locations. Same as CCorA when X contains group identifiers. Equivalent analysis is regression of X on Y , provided X contains orthogonal contrast vectors.

• Finds linear combination of axes in Q t correlated, or (if X and in X that are maximally contains group identifiers) finds axis in PCO space that maximises differences among group locations.

PRINCIPAL RESPONSE CURVES (PRC)

van der Brink, P. & ter Braak, C.J.F. (1999) Environmental Toxicology & Chemistry 18, 138-148 van der Brink, P. & ter Braak, C.J.F. (1998) Aquatic Ecology 32, 163-178 PRC is a means of analysing repeated measurement designs and of testing and displaying optimal treatment effects that change across time.

Based on RDA (= reduced rank regression) that is adjusted for changes across time in the control treatment. Allows focus on time-dependent treatment effects. Plot resulting principal component against time in PRC diagram.

Developed in ecotoxicology; also used in repeated measures in experimental ecology and in descriptive ecology where spatial replication is substituted for temporal replication.

Highlights differences in measurement end-points betweeen treatments and the reference control.

PRC MODEL

d(i)tk

otk



d(i)tk

= Y

otk

+ b

+ 

d(i)tk

= mean abundance of taxon = weight of species

with respect to

c dt

where = abundance counts of taxon

at time

in replicate

of treatment

d k

in controls ( o ) at time = principal response of treatment

at time

(PRC) = error term with mean of zero and variance  2 k

Modelling the abundance of particular species as a sum of three terms, mean abundance in control, a treatment effect, and an error term.

Data input - species data (often log transformed) for different treatments at different times - predictor variables of dummy variables (1/0) to indicate all combinations of treatment and sampling time ('indicator variables') - covariables of dummy variables to indicate sampling time Do partial RDA with responses, predictors, and covariables and delete the predictor variables that represent the control. This ensures that the treatment effects are expressed as deviations from the control.

PRC MODEL (continued)

Simple example - three treatments: C = control, L = low dosage (not rep licated), H = high dosage sampled at four times (W0, W1, W2, W3), six species.

PRC MODEL (continued)

* * * * deleted in the RDA

PRC PLOTS

One curve for each treatment expressed as deviation from the control. Species weights (b

) allow species interpretation. Higher the weight, more the actual species response is likely to follow the PRC pattern, because the response pattern = b

. Taxa with high negative weight are inferred to show opposite pattern. Taxa with near zero weight show no response.

Significance of PRC can be tested by Monte Carlo permutation of the whole time series within each treatment.

Can use the second RDA axis to generate a second PRC diagram to rank 2 model.

PRC PLOTS (continued)

See maximum deviation from control after 4 weeks, maximum effect larger for 150  g/l treatment than for 50  g/l treatment.

Chlamydomonas has high negative weight and this had highest abundances in high doses after treatment began.

Principal response curves resulting from the analysis of the example data set, indicating the effects of the herbicide linuron on the phyto plankton community. Of all variance, 47% could be attributed to sampling date, and is on the horizontal axis. Of all variance, 30% could be attributed to treatment. Of the variance explained by treatment, 23% is displayed on the vertical axis. The lines represent the course of the treatment levels in time. The species weight (b response curves.

) can be interpreted as the affinity of the taxon with the principal

PRINCIPAL RESPONSE CURVES AND ANALYSIS OF MONITORING DATA

PRC usually used with experimental data. Can be used with (bio)monitoring data.

Samples at several dates at several sites of a river, some upstream of a sewage treatment plant (STP) (300 m, 100 m), in the STP outlet, and some downstream (100 m, 1 km). 795 samples, 5 sites, 1994-2002.

PRC using sampling month as covariable, product of sampling month and site as explanatory variables. Used STP outlet as the reference site.

Van der Brink, P. et al. (2003) Austr. J. Ecotoxicology 9: 141-156

Principal Response Curves indicating the effects of the outlet of a sewage treatment plant on some monthly averages of physico-chemical characteristics of a river. Of all variance, 24% could be attributed to between month variation; this is displayed on the horizontal axis. 57% of all variance could be allocated to between site differences, the remaining 19% to within-month variation. Of the between-site variation, 58% is displayed on the vertical axis. The lines represent the course of the sites in time with respect to the outlet. The weight of the physico-chemical variable (b

) can be interpreted as the affinity of the variables with the Principal Response Curves (c

See biggest differences for the two upstream sites, with lower NO in the upstream sites due to pollution.

x , total N, conductivity, salinity, total P, and temperature and higher values of turbidity and faecal coliforms. STP outlet leads to increases in N, P, temperature, etc. Downstream values decrease but are not as low as upstream sites. STP successfully reduces faecal coliforms as their values are higher

PRINCIPAL RESPONSE CURVES – A SUMMARY

Filters out mean abundance patterns across time in the control. Focuses on deviation between treatment and control. PRC displays major patterns in those deviations and provides good summary of response curves of individual taxa.

PRC helps to highlight 'signal' from 'noise' in ecological data in replicated experimental studies.

Simplified RDA - simplified by representing the time trajectory for the controls as a horizontal line and taking the control as the reference to which other treatments are compared.

PRC gives simple representation of how treatment effects develop over time at the assemblage level.

ECOLOGICAL APPLICATIONS OF PRINCIPAL RESPONSE CURVES

Vandvik, V. 2004. J. Ecology 92: 86-96 – revegetation of gaps in replicated successional series 0, 10, and 40 years after abandonment.

Vandvik, V et al. 2005. J. Applied Ecology 42: 139-149 – post-fire succession in heathlands.

Heegaard, E. & Vandvik, V. 2004. Oecologia 139: 459-466 – space-for-time substitution to analyse compositional differences in local ridge-snowbed gradients in 100 m altitudinal bands from 1140 to 1550 m.

REF

POLYNOMIAL RDA AND POLYNOMIAL CCA

Makarenkov & Legendre (2002) Ecology 83, 1146-1161 Instead of using multiple linear regression to model the relationship between each response variable y in Y and the explanatory variables X , use polynomial regression Y = b 1 x 1 + b 2 x 1 2 + b 3 x 2 + b 4 x 2 2 ... ... b p x p 2 Aim is to develop a multivariate reduced-rank regression model that produces a noticeable increase in the amount of variation in linear RDA and CCA.

Y accounted for by the model compared to standard REF

POLYRDACCA

www.fas.umontreal.ca/biol/legendre REF

REF REF REF RDA - can be viewed as a series of multiple linear regressions followed by an eigenvalue decomposition of the table of fitted values CCA approximate unimodal responses of species to environmental gradients through a chi-square transformation of the species data. Assumes linear relationships between Y and X .

Resulting fitted values are no longer a linear combination of the explanatory variables but a polynomial combination of variables, e.g.

b 1 x 1 + b 2 x 1 2 + b 3 x 2 + b 4 x 2 2 + b 5 x 1 x 2 etc.

Tries to use polynomial relationships to model non-linear relationships only when these cannot be modelled well by simple linear models.

Need to reduce the number of predictor variables so as to avoid too many terms and 'overfitting'. Overfitting occurs if number of regression terms > n-1 where n = number of observations.

Polynomial algorithm proceeds by successively reducing the matrix of explanatory variables while increasing the R 2 for the response variable by introducing polynomial terms.

REF

Display of results

1. Plot individual terms by biplot arrows. Can easily get too many arrows. May only plot those with a correlation greater than some pre-specified value.

2. Represent each explanatory variable by a single arrow corresponding to the multiple correlation of x and x 2 with the ordination axes.

REF REF

Tests of significance

1. Test statistical significance of linear model by permutation tests as in

CANOCO

2. 3. Test statistical significance of polynomial model by related permutation test.

Test if the difference in variance accounted for between the polynomial and the linear models is significant and hence if the simpler linear model is the minimal adequate model.

REF

REF Linear RDA Polynomial RDA REF P = 0.001 (999 permutations) REF Notes: Matrix Y: hunting spider species 1-12. Matrix X: water content, reflection of soil surface. Either set of biplot scores can be used to represent the environmental variables in biplots. Users of the program may also request the matrix of regression coefficients B of the multiple linear regressions of Y on X (if classical linear RDA or CCA is computed) or the polynomial coefficients for each response variable y of Y (if polynomial RDA or CCA is used). The program may also carry out permutation tests of the significance of the linear and polynomial models, as well as the significance of the difference in variance accounted for between the two models.

Significance of difference in explained variance p = 0.002 (999 permutations). Use the polynomial RDA model.

REF

REF Linear RDA 35.4% Polynomial RDA 53.7% REF REF Full arrows = biplot scores of environmental variables in polynomial: Dashed arrows = biplot scores based on multiple correlations of (water, water 2 ) and (reflection, relection 2 ) with the axes.

REF

CCA/RDA AS PREDICTIVE TOOLS

Prediction is important challenge in environmental science.

• Given environmental shift, how will species respond?

• Given environmental data only (e.g. satellite image data), what biotic assemblages could be expected?

Conventional CCA/RDA – description and modelling Y m and X m  Model m where subscript 'Palaeo' CCA/RDA – modelling and reconstruction m = modern Y m Y o and X m  and Model m Model m  X o where subscript o 'palaeo' data = fossil or (Lecture 8 on Environmental Reconstruction) Predictive CCA/RDA – modelling and prediction Y m X f and and X m  Model m Model m  Y f where subscript (predicted) data f = future

Gottfried et al. 1998. Arctic and Alpine Research 30: 207-211 Schrankogel (3497 m) Tyrol, eastern central Alps.

1000 1x1 m plots between 2830 and 3100 m, ecotonal transition between alpine zone (vegetation cover >50%) and nival zone (vegetation cover <10%).

Vegetation data- species +/- and relative abundance of 19 species Environmental data – Digital Elevation Model (DEM) with pixel size of plots in GIS ARC/INFO - gives 17 topographic descriptors at 10 resolutions plus altitude.

Gottfried et al. 1998

CCA with forward selection to give 37 predictor variables Calculated CCA sample scores for 650,000 cells of DEM area as weighted linear combination of environmental variables times the canonical coefficients.

For each predicted environment of each cell, estimated which of the 1000 plots it is closest to CCA space.

Gottfried et al. 1999

Extrapolate vegetation data from those plots to the 650,000 cells to predict species distributions, vegetation types, species richness, etc.

To evaluate predictions, did 10-fold cross-validation, namely model with 90% of the plots, predict with left-out 10%, and repeat 10 times. Compare predictions with actual observed data. Also calculated Cohen's kappa statistic between observed and predicted distributions (0 = uncorrelated, 1 = perfect match).

Total inertia 1.79

Model performance CA CCA Axis 1 0.41

0.28

Axis 2 0.21

0.11

Constrained inertia 0.68

Species variance 34.3% 21.8% 38% explained by topography Species fell into five groups of kappa and other performance values Kappa > 0.5

Kappa > 0.4

Kappa > 0.3

Kappa > 0.25

Kappa < 0.1

e.g. Carex curvula, Veronica alpina e.g. Androsace alpina, Oreochloa disticha e.g. Saxifraga oppositifolia, Primula glutinosa e.g. Ranunculus glacialis, Cerastium uniflorum e.g Poa laxa

Predicted distributions Gottfried et al. 1998

Predicted vegetation types Gottfried et al. 1998

Best predictors are topographic measures of roughness and curvature rather than simple elevation, slope, or exposure.

Modelled richness patterns decline with altitude but with a maximum richness at the alpine-nival ecotone.

What might happen under future climate warming of 1ºC or 2ºC?

Gottfried et al. 1999 Diversity and Distributions 5: 241-251.

Calculated altitudinal lapse rate using 33 temperature loggers.

Assume that the altitudinal limits are determined by temperature. Knowing the temperature lapse rate, predict species distributions for +1ºC and +2ºC temperature increases.

Done by increasing the values of the altitude variable in the environmental data for the sample plots corresponding to the lapse rate. Repeated the CCA/GIS interpolation/mapping.

Assumes that species growing at lower altitudes and hence warmer situations will occur in a future warmer climate in the same topographical conditions Predicted distribution patterns at +0º, +1º and +2º Gottfried et al. 1999

Predicted distribution patterns at +0º, +1º and +2º Gottfried et al. 1999

Predictions - 19 species modelled, about 2 will become extinct will be a reduction in genetic diversity as some 'topographical forms' will be lost Dirnböck et al. 2003 Applied Vegetation Science 6: 85-96.

Same approach but with topographic descriptors and infra-red spectral data, and 3 nearest neighbours to predict rather than 1.

Predictive mapping of 17 vegetation types between 1600 m and 2277 m on Hochschwab, eastern Alps.

Dirnböck et al. 2003

69.4% accuracy, Cohen kappa of 0.64

Topography good predictors of different alpine grasslands Infra-red spectra good predictors of different pioneer vegetation types Unexplained variation - land-use history, soil variation especially nutrients like N, P, and K

CANONICAL CORRESPONDENCE ANALYSIS

Builds on:

1) Weighted averaging of indicator species and extends WA to the simultaneous analysis of many species and many environmental variables.

2) Reciprocal averaging (= correspondence analysis) by adding the statistical methodology of regression. General framework of estimation and statistical testing of the effects of explanatory variables on biological communities.

Major Uses:

1. Identify environmental gradients in ecological data-sets.

2. In palaeoecology, used as a preliminary to determine what variables influence present-day community compositions well enough to warrant palaeoenvironmental reconstruction.

3. Add 'fossil' samples into modern 'environmental' space.

4. Study seasonal and spatial and temporal variation in communities and how this variation can be explained by environmental variation. Variance can be decomposed into seasonal, temporal, spatial, environmental and random components.

5. Niche analysis – niche-space partitioning where species probability or abundance is unimodal function of environment.

6. Impact studies.

7. Predictive studies.

8. Experimental data analysis.

Powerful alternative to multivariate analysis of variance MANOVA.

e.g. analysis of BACI (before-after-control-impact) studies with and without replication of the impacted site e.g. repeated measurement designs e.g. experimental plot (= block) designs e.g. split-plot designs

See

ter Braak C.J.F. & Verdonschot P.F.M (1995) Aquatic Sciences 57, 255–289 Canonical correspondence analysis and related multivariate methods in aquatic ecology

NON-LINEAR CANONICAL ANALYSIS OF PRINCIPAL CO-ORDINATES

Millar, Anderson & Zunan 2005 Ecology 86: 2245-2251 All the methods discussed so far (RDA, CCA, Canonical correlation analysis, distance-based RDA, CAP, non-linear RDA or CCA) are intrinsically linear with respect to the environmental variables and are fitting a linear model to the distances of dissimilarities implicit in the particular method (e.g. Euclidean distances in RDA, chi squared distance in CCA, any distance measure in distance-based RDA or CAP).

In some ecological situations, the relationship between community structure and environmental variables may be intrinsically NON LINEAR.

e.g. community composition at different distances from a pollution source. At what point along the gradient does the impact become negligible? Environmental variable (distance from pollution source) has non-linear, exponential effect.

e.g. community change through time. After disturbance, succession occurs, with invasion and replacement of species with time. At what point in time does community composition stabilise?

Non-linear canonical analysis of principal co-ordinates (NCAP) of

relation to

in 1. Principal co-ordinates analysis of

to obtain

co-ordinates using any distance measure.

orthogonal principal 2. Non-linear canonical correlation analysis of

maximise correlation of

in relation to with respect to a non-linear

transformation of X – to Use non-linear link function (e.g. exponential) just as we do when we extend a traditional univariate linear model to a GLM through the use of a link function. Non-linear optimisation procedures as in GLM.

3. Determine the smallest number of principal co-ordinates using BIC in the NCAP context BIC  log  1 

2   

log   



where

= number of parameters,

determination = constant,

2 = coefficient of

4. Test null hypothesis of no association between

by randomisation tests of rows of

and

and 5. Test if the non-linear gradient fitted is preferred to a linear gradient. Can view the non-linear model as the 'full' model and the linear model as the 'reduced' model. Use a pseudo-F statistic to test the null hypothesis of a linear gradient and derive its significance by permutation of rows of

6. Determine the approximate 95% confidence intervals of the regression coefficients by bootstrapping.

Unlike non-linear canonical correlation analysis or non-linear polynomial RDA or CCA, no limit to the number of variables in

as the analysis is based on a distance measure of choice.

No explicit assumptions about distributions in the biological data (e.g. counts, presence/absence) as randomisation tests and bootstrapping are used to test formally the statistical significance of the observed non-linear gradient.

Virtually any non-linear form of gradient can be used along with any reasonable distance measure.

Of considerable potential in the analysis of multivariate non linear ecological systems.

CANONICAL GAUSSIAN ORDINATION Something New but a bit Heavy!

CCA can be estimated in two ways (1) weighted averaging estimation – heuristic, rapid, computationally simple, approximation to (2) maximum likelihood (ML) estimation – 'ideal' but theoretically rigorous and computationally difficult CANOCO implements the weighted averaging approach. ML canonical Gaussian ordination (CGO) not computationally possible until very recently.

Yee (2004) Ecological Monographs 74: 685-701 Quadratic reduced-rank vector generalised linear models (QRR-VGLMs) as a way of implementing CGO.

Basis is RR-VGLMs where by reduced-rank regression multivariate regression (many responses and many predictors) is implemented in the framework of GLM with specified link functions, error distributions, etc.

QRR-VGLMs are simply RR-VGLMs but with the addition of a quadratic form to each linear predictor, so that bell-shaped responses can be modelled as functions of environmental variables or gradients.

VGLM package in R

Hunting spider data – 12 species, 28 sites

QRR-VGLM – assume equal or unequal species tolerances

QRR-VGLM results very similar to CCA results

Advantages of canonical Gaussian ordination (CGO) 1. CGO has clear statistical assumptions 2. CGO provides estimates of maxima, optima, and tolerances of species

3. CGO avoids complications, heuristics, and statistical simplifications and inefficiencies implicit in weighted averaging and detrending.

4. CGO is built on firm statistical theory.

5. CGO is more flexible than CCA.

CCA makes four assumptions (equal tolerances, equal maxima, homogeneously distributed data, site scores homogeneously distributed over the gradients – assumptions 1-3 imply a species-packing model). These assumptions cannot simultaneously hold with real ecological data.

CGO does not make any of these assumptions though the equal-tolerance assumption makes the ordination plot more easily interpretable.

6. Residual diagnostic plots are available. (In theory possible in CCA.)

Disadvantages of CGO are: 1. Statistical assumption of symmetric unimodal Gaussian responses may be unrealistic in practice – skewed or multimodal responses.

2. Sensitive to outliers, sparse data, multicollinearity, high leverage points, as in other GLMs.

3. The global log-likelihood function that is maximised may have several maxima, and convergence to a local solution can never be ruled out.

4. Good initial starting values can be difficult to obtain, especially as the number of species and environmental variables increases.

5. Difficult to establish statistically the rank R of the data.

6. Numerically very intensive. CPU power continues to increase exponentially, so the size of data that can be handled will also increase.

QRR-VGLM advantages: 1. Basic framework is very broad, so error distributions such as binomial, negative binomial, and zero-inflated Poisson can be handled, along with compositional data. As Yee (2006) say "CGO operates on various types of species data such as presence/absence data and Poisson counts; CCA treats all types of data the same, which, to purists, is a bit like using a baseball bat to play golf, squash, and tennis." 2. Can accommodate linear combinations of non-linear functions of explanatory variables. Can include factors and functions such as orthogonal polynomials and B-splines. Polynomial CCA is ad hoc and too limited.

3. Any rank R can be fitted in theory.

4. Approximate statistical significance of regression parameters can be tested.

5. Statistical hypothesis testing in terms of deviance is possible.

6. Overdispersion and underdispersion in the data can be detected. Not possible in CCA (Overdispersion - variance in data exceeds nominal variance under an assumed model: underdispersion is opposite).

QRR-VGLM disadvantages 1. Requires large amounts of memory as large model matrices are constructed to perform the fitting.

2. Specialist software needed – VGAM package in R.

3. Convergence problems may occur, especially with 'dirty' real data.

4. Much skill and statistical understanding are needed to use QRR-VGLMs.

5. 'Exact' standard errors for all estimated parameters currently not possible.

6. Currently computationally possible for rank 2 with 500 sites x 5 species x 5 predictors or 60 sites x 10 species x 10 predictors. Computational needs increase very rapidly with number of species and with the rank.

Important message, not mentioned by Yee (2004, 2006) is that the CCA solution is very similar to the CGO QRR-VGLM (and the CAO RR-VGAM) model results.

Reassuring for all CCA users!

Remarkable just how robust CCA is and how good weighted average estimation is.

Next step: to move from GLM framework of specified response models to GAM framework with smooth data-driven functions i.e. RR-VGAMs.

Yee (2006) Ecology 97: 203-213

CONSTRAINED ADDITIVE ORDINATION (CAO)

CAO, unlike CA, DCA, CCA, CGO, etc, does NOT assume any underlying species response model (e.g. unimodal, symmetric, etc).

Computes optimal gradients and flexible GAM-like response curves. Can see response curves as they really are in relation to the major gradients.

With one gradient, CAO is a generalisation of constrained quadratic (or Gaussian) ordination (CQO or CGO).

Replaces symmetric bell-shaped responses in CQO (=CGO) with completely flexible smooth curves estimated using smoothers.

CAO models are GAMs fitted to a very small number of latent variables.

Data-driven rather than model-driven, so CAO allows the data "to speak for themselves". No assumptions.

Can now think of ordination methods based on modern regression procedures as unconstrained or constrained (=canonical) and linear, quadratic, or smooth depending on assumed underlying responses.

VGAM package in R

Assumptions of CCA, CQO (=CGO), and CAO

Equal tolerances Equal maxima Species optima uniformly distributed Site scores uniformly distributed Symmetric bell-shaped responses CCA

+ + + + +

CQO

+/ +/ +

CAO

Ordination (O) procedures based on modern regression procedures (Yee 2006) Shape Unconstrained (U) Constrained (C) Linear (L) Quadratic (Q) Smooth (GAM) Unconstrained linear ordination e.g. unconstrained vector generalised linear model (U-VGLM) Unconstrained quadratic ordination (UQO) e.g. quadratic unconstrained VGLM (QU-VGLM), Gaussian ordination Unconstrained additive ordination (UAO) e.g. U VGAM, prinicipal curves Constrained LO e.g. reduced rank regression (=RDA), generalised reduced rank regression, RR-VGLM Constrained quadratic ordination (CQO=CGO) ordination e.g. QRR-VGLM, Gaussian logit Constrained additive ordination (CAO) curves e.g. RR VGAM, constrained principal VGLM = vector generalised linear model; VGAM = vector generalised additive model; RR = reduced rank

Yee (2006) proposes that ordination, prediction (estimating

and calibrating (estimating

from

from x ) can all be performed from the ), same model.

Can, in theory, have Ordination CLO CQO (= CGO) CAO Prediction CLP CQP CAP Calibration CLC CQC CAC

Reduced-rank vector-based generalised additive models (RR-VGLM) Approximates each surface in R-dimensional latent variable space by an additive model.

Interpretable and avoids the 'curse of dimensionality' when R becomes large. Splines or LOESS smoothers are unusable if R is large and are uninterpretable if R > 2.

Technical problems arise in how much flexibility to allow in the GAM smoother, the effective non-linear degrees of freedom.

Smooth curve with zero non-linear degrees of freedom is linear, so CAO = CLO.

Values between 1.5 and 2.5 non-linear degrees of freedom usually necessary.

Hunting spider data: 12 species 28 sites 6 environmental variables Each species has three non-linear degrees of freedom in rank-1 CAO model Clear unimodal response curves Similar to De'ath's (1999) principal curves for the same data Yee (2006)

Yee (2006) CAO model. Very similar to CQO (and CCA) results CQO (= CGO) preferred to CAO as the final model as it is simpler and parameters such as species tolerances are well defined (as in GLM compared to GAM).

CAO best used as an exploratory tool and to help CQO work better.

RR-VGLMs more powerful than ordinary GAMs because they search for optimal latent variables before applying the additive model.

GAMs are a special case of RR-VGAMs because GAMs smooth one environmental variable at a time whereas each RR-VGAM smooth is on a composite of the environmental variables after searching for the optimal combination or weighting.

At present, software only handles rank-1 RR-VGAMs.

Cannot handle interactions in X.

Are "fragile with dirty data" (Yee 2006).

SOFTWARE FOR CONSTRAINED ORDINATIONS

canonical correlation analysis distance-based redundancy analysis via principal co ordinates analysis

CANOCO + CANODRAW

[CANCOR, CAP]

DISTPCOA POLYRDACCA

constrained principal co-ordinates analysis polynomial RDA and CCA R(VEGAN) (CCA, RDA) R (Non-linear CAP) R(VGLM) (CGO) R(VGAM) CAO)

Pie symbols plot

CANODRAW 4

Response curves fitted using GAM Sample diagram with principal response curves T-value biplot Attribute plot Biplot with environmental variables & sites Isolines in RDA ordination diagram

1987 2002 2003 1987

Mark Hill Thomas Yee Cajo ter Braak Marti Anderson Petr Šmilauer

NUMERICAL ANALYSIS OF BIOLOGICAL AND ENVIRONMENTAL DATA

Transcript NUMERICAL ANALYSIS OF BIOLOGICAL AND ENVIRONMENTAL DATA

NUMERICAL ANALYSIS OF BIOLOGICAL AND ENVIRONMENTAL DATA

Lecture 7.

Direct Gradient Analysis

DIRECT GRADIENT ANALYSIS

BASIS OF CLASSICAL ORDINATION INTERPRETATION AND ENVIRONMENT

TWO-STEP APPROACH OF INDIRECT GRADIENT ANALYSIS

CANONICAL ORDINATION TECHNIQUES

PRIMARY DATA IN GRADIENT ANALYSIS

.....

ALTERNATING REGRESSION ALGORITHMS

CANONICAL CORRESPONDENCE ANALYSIS

CANONICAL OR CONSTRAINED CORRESPONDENCE ANALYSIS (CCA)

BASIC TERMS

SUMMARY OF DUNE MEADOW DATA

DCA

CCA

CANONICAL CORRESPONDENCE ANALYSIS

BIPLOT PREDICTION OF ENVIRONMENTAL VARIABLES

CLASS VARIABLES

CCA: JOINT PLOTS AND TRIPLOTS

CANOCO

CONVERGENCE CRITERIA IN EIGENVALUE EXTRACTION

OTHER CCA TOPICS

REPRESENTATION OF CLASS VARIABLES (1/0) IN CCA

ROBUSTNESS OF CANONICAL CORRESPONDENCE ANALYSIS

INFLUENCE OF NOISY ENVIRONMENTAL DATA ON CANONICAL CORRESPONDENCE ANALYSIS

LC OR WA SCORES?

DATA ORDERINGS

CCA DIAGRAM

DEFAULT CCA PLOT

SCALING IN CCA

INTERPRETATION OF CCA PLOTS

CCA EXAMPLE

CANONICAL CORRESPONDENCE ANALYSIS (CCA) - A SUMMARY

REDUNDANCY ANALYSIS – CONSTRAINED PCA

PCA and RDA comparisons

BIPLOT INTERPRETATION

PROPOSED NEW SCALING FOR CCA AND RDA

ALTERNATIVES TO ENVIRONMENTAL VECTORS IN CCA AND RDA

STATISTICAL TESTING OF CONSTRAINED ORDINATION RESULTS

STATISTICAL SIGNIFICANCE OF CONSTRAINING VARIABLES

PARTIAL CONSTRAINED ORDINATIONS (Partial CCA, RDA, etc)

PARTIAL CCA

PARTIAL ORDINATION ANALYSIS (Partial PCA, CA, DCA)

PARTIAL ORDINATIONS

COVARIABLES IN CCA AND RDA

PARTITIONING VARIATION

VARIATION PARTITIONING OR DECOMPOSITION WITH 3 OR MORE SETS OF PREDICTOR (EXPLANATORY) VARIABLES

'NEGATIVE' VARIANCES

UNBIASED ESTIMATES OF VARIANCE COMPONENTS

Partial linear regression

Partial canonical analysis

ENVIRONMENTAL CONSTRAINTS AND CURVATURE IN ORDINATIONS

Classification of gradient analysis techniques by type of problem, response model and method of estimation

ENVIRONMENTAL VARIABLES IN CONSTRAINED ORDINATIONS

SELECTING ENVIRONMENTAL VARIABLES IN CON STRAINED ORDINATION ANALYSIS (e.g. CCA, RDA)

MINIMAL ADEQUATE MODEL IN CCA

VARIANCE INFLATION FACTORS (VIF)

'AIC' FOR MODEL SELECTION IN CCA AND RDA

STAGES IN 'AIC' MODEL SELECTION IN CCA AND RDA

'AIC' MODEL SELECTION

2 log

2

(

sum of all eigenvalue s

trace

)

log

1 2

log(

)

Can determine not only the 'best' model but rank the different underlying hypotheses according to AIC criteria of parsimony. Spatial models 2 and 5, namely depth stratification but no difference between sheltered and exposed stations, are most appropriate for these data.

CANONICAL CORRELATION ANALYSIS – CANCOR

COMPARISON WITH CANONICAL CORRELATION ANALYSIS

CONSTRAINED LINEAR ORDINATION (PCA FRAMEWORK)

DISTANCE-BASED REDUNDANCY ANALYSIS

DISTANCE-BASED MULTIVARIATE ANALYSIS FOR A LINEAR MODEL

CANONICAL ANALYSIS OF PRINCIPAL CO-ORDINATES