Transcript NUMERICAL ANALYSIS OF BIOLOGICAL AND ENVIRONMENTAL …

Ordination Analysis II –
Direct Gradient Analysis
John Birks
Quantitative Methods in Palaeoecology
and Palaeoclimatology
PAGES Valdivia October 2010
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 coordinates
Principal response curves
CCA/RDA as predictive tools
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)
Environmental
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 .
Vegetational data
DUNE-MEADOW DATA
DCA
axis 2
Indirect
analysis
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=b0 + b1x1 + b2x2
Angle ()with axis 1 = arctan(b2 / b1)
Indirect
analysis
Site scores of the second DCA axis plotted against the amount of manure.
Indirect
analysis
Correlation coefficients (100  r) of the environmental
variables for the four first DCA axes for the Dune
Meadow Data
Variable
Axes
1
2
3
4
1 A1
58
24
7
9
2 moisture
76
57
7
-7
3 use
35
-21
-3
-5
6
-68
-7
-64
5 SF
22
-29
5
-60
6 BF
-28
-24
39
22
7 HF
-22
-26
-55
-14
8 NM
21
73
17
56
0.54
0.29
0.08
0.05
4 manure
Eigenvalue
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
Parameter
constant
c0
A1
c1
moisture
c2
use
c3
manure
c4
Estimate
–2.32
0.14
0.38
0.31
–0.00
s.e.
0.50
0.08
0.09
0.22
0.12
ANOVA table
d.f.
Regression
Residual 15
Total
19
s.s.
4
6.2
23.2
m.s.
17.0
0.41
1.22
F
4.25
R2 = 0.75
R2adj = 0.66
t
–4.62
1.71
4.08
1.37
–0.01
Indirect
analysis
10,6
Ey1 = b0 + b1x1 + b2x2 + ...bnxn
CA axis 1
environmental variables
x = 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)
 constrained CA
(Detrended CCA)

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.
Indirect GA
Species
PRIMARY DATA IN GRADIENT ANALYSIS
Abundances
or
+/variables
Response variables
Values
Env. vars
Direct GA
PLUS
Predictor or explanatory variables
Classes
Artificial example of unimodal response curves of five species (A-E) with respect to standardised 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
e.g.
all possible linear combinations
xj  co  c1z1 j  c2z2 j  c3z3 j  .....
zj = environmental variable at site j
c = weights
xj = resulting ‘compound’ environmental variable
CCA selects linear combination of environmental variables that maximises
dispersion of species scores, i.e. chooses the best weights (ci) 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
REF
Algorithm
1) Start with arbitrary, but unequal, site scores xi.
2) Calculate species scores by weighted averaging of site scores.
n
n
i 1
i 1
uk   y ik xi /  y ik
3) Calculate new site scores by weighted averaging of species scores.
m
m
k 1
k 1
x   y ikuk /  y ik

i
[So far, two-way weighted average algorithm of correspondence
analysis].
REF
REF
REF
REF
4) Obtain regression coefficients of site scores on the environmental variables by
weighted multiple regression.
b  Z 1RZ  Z 1Rx 
1
where
b and x* are column vectors
Z is environmental data n x (q +1)
R is n x n matrix with site totals in diagonal
5) Calculate new site scores
x  zb
or
x  b0   bz
6) Centre and standardise site scores so that:
n
 y i  xi  0 and
i 1
n
2
y
x
 i i  1
i 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
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
A1
MoistureManagement
Use
horizon
class
type
2.8
3.5
4.3
4.2
6.3
4.3
2.8
4.5
3.7
3.3
3.5
5.8
6.0
9.3
11.5
5.7
4.0
4.6*
3.7
3.5
1
1
2
2
1
1
1
5
4
2
1
4
5
5
5
5
2
1
5
5
SF
BF
SF
SF
HF
HF
HF
HF
HF
BF
BF
SF
SF
NM
NM
SF
NM
NM
NM
NM
2
2
2
2
1
2
3
3
1
1
3
2
2
3
2
3
1
1
1
1
Manure
class
4
2
4
4
2
2
3
3
1
1
1
2*
3
0
0
3
0
0
0
0
DCA
DCA
axis 2
DCA
axis 1
DCA ordination diagram of the Dune Meadow Data
DCA
Axis Two
2 = 0.29
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
1
0.54
0.46
2
0.40
0.29
R axis 1
0.87
0.96
R axis 2
0.83
0.89
CANONICAL CORRESPONDENCE ANALYSIS
Canonical correspondence analysis: canonical coefficients (100 x c) and intraset 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
A1
9
Moisture 71
Use
25
Manure
-7
SF
BF
-9
HF
18
NM
20
Correlations
Axis 2
Axis 1
-37
-29
5
-27
16
19
92
57
93
21
-30
16
-37
-36
56
Axis 2
-17
-14
-41
-79
-70
15
-12
76
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.
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)
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
CANOCO
Summary
Axes
Axes
1
2
3
4 Total inertia
Eigenvalues
.461 .298 .160 .134
2.115
Species-environment
.958 .902 .855 .889
correlations
Cumulative percentage variance
of species data
21.8 35.9 43.5 49.8
of species-environment 37.8 62.3 75.4 86.3
'Fitted' species data
relation
Sum of all unconstrained eigenvalues = inertia
2.115
Sum of all canonical eigenvalues = species-environment
1.220
relation
1  2
 canonical eigenvalues
 100
Rules of thumb:
>0.30
strong gradient
>0.40
good niche separation of species
OTHER CCA TOPICS
1) Environmental variables
continuous
–
biplot arrows
classes
–
centroid (weighted average) of sites belonging to that class
2) CA approximates ML solution of Gaussian model
CCA approximates ML solution of Gaussian model if CA axis is close to the linear combination of environmental variables. [Johnson & Altman (1999) Environmetrics 10, 39-52]
In CCA species compositional data are explained through a Gaussian unimodal response
model in which the explanatory variable is a linear combination of environmental
variables.
3) CCA – very robust, major assumption is that response model is UNIMODAL.
(Tolerances, maxima, and location of optima can be violated - see Johnson & Altman
1999)
4) Constraints become less and less strict the more environmental variables there are. If
q, number of environmental variables ≥ number of samples -1, no real constraints and
CCA = CA.
5) Arch effect may crop up. Detrending (by polynomials) DCCA. Useful for estimating
gradient lengths (use segments).
6) Arch effect can often be removed by dropping superfluous environmental variables,
especially those highly correlated with the arched axis.
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  point in q-dimensional space of environmental variables.
NICHE.
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
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
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)
2
2
(3)
(4)
10
2
+
10
environmental variables
environmental variables
X and Y co-ordinates
TENxTEN
with added noise
NOISMOD
(random number mean = 0, variance 17%)
added to each cell
random environmental variables
NOIS1O
environmental variables with added noise from
NOISMOD
random environmental variables from
NOIS 10
NOISBOTH
NOISFULL
(5) 99 random environmental variables
NOISFULL –
‘Species-environment’ correlation increases as number of random variables
increases for axis 1 and 2.
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!
TEN x TEN
NOISMOD
NOISIO
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.
NOISFULL
Monte Carlo tests 1
TENxTEN
*
0.77
NOISMOD
*
0.65
NOISE 10
ns
0.20
NOISBOTH
*
0.65
NOISFULL
ns
0.12
(99 env vars)
Linear combination site
scores
WA site scores
2
0.79
0.65
0.12
0.65
0.08
r1
0.98
0.90
0.49
0.91
1.0
r2
0.9
0.8
0.4
0.8
1.0
best fit of species abundances to
the environmental data
best represent the assemblage
structure
LC scores WA scores
+
–
+
–
Sensitive to noise
True direct gradient analysis
(multivariate regression)
Aim to describe biological
+
–
variation in relation to
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?
MIKE PALMER
"Use LC scores, because they give
the best fit with the environment
and WA scores are a step from CCA
towards CA."
BRUCE MCCUNE
"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)
SCALING IN CCA
Hill scaling
–1
Emphasis on
1 Species x sites
SITES
3 Sites x sites
2
SPECIES
Rel abundances
2 Species x species
Default scaling
–
Turnover distances
Fitted abundances (rel)
Chi-squared distances
–
Quant env vars
4 Sites x env
vars3
–
Values of env vars
5 Species x env vars
Weighted averages Weighted averages
6 Env vars x env vars
Effects2
Correlations2
fitted by least
squares
1 by
centroid
principle
2 change
in site
scores if env
variable changes
are one standard
deviation
Qualit env vars
7 Sites x env classes4 Membership1
Membership1
8 Species x env classes Rel total abund
Rel total abund
9 Env vars x env classes
Mean values of env vars
10 Env classes x env
classes
–
Turnover distances
–
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 4. The weighted sum of
squares of sites scores of an axis is equal to /(1-) with  its eigenvalue
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 1/(1-) 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).
Note that symmetric scaling (=3) has many optimal properties
(Gabriel, 2002; ter Braak, personal communication)
REF
REF
-1: focus on sites
Hill's scaling
Interpreta- 2: focus on species
tion
biplot scaling of CCA
1. species x sitesa
Rel. Abundancesb,c
CENTROID
Fitted rel. abund.b
BIPLOT rule or
CENTROID rule
2. species x species
-
UNKNOWN
-square
distancesd
DISTANCE rule
3. sites x sites
Turnover
distancesc,e
DISTANCE
f
DISTANCE rule
Scaling
Interpretation
Quantitative env. vars
4. sites x env. varsg
-
UNKNOWN
Values of env.vars
BIPLOT rule
5. species x env. vars
Weighted
averages
BIPLOT
Weighted averages
BIPLOT rule
6. env.vars x env. vars
Effectsh
? BIPLOT
Correlations
BIPLOT rule
7. sites x env. classesi
Membershipk
CENTROID
Membershipk
CENTROID rule
8. species x env. cls.
Rel. total
abund.c,b
CENTROID
Rel. total abund.b
CENTROID rule
9. env.vars x env.
classes
-
UNKNOWN
Mean values of
env. vars
BIPLOT rule
10. env. classes x env.
Turnover
DISTANCE
f
DISTANCE rule
Qualitative env. vars
REF
REF
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.
a
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
In the definition of this coefficient, abundance must be replaced by fitted abundance
values, because CCA is correspondence analysis of fitted abundance values
d
e
No explicit formula known
f
Chi-square distances, provided the eigenvalues of the axes are of the same magnitude
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
(yi+). The scores in scaling -1 are {(1-)}½ times those in scaling 2.
g
Effect is defined as the change in site scores if the environmental variable changes one
standard deviation in value (while neglecting the other variables).
h
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
CCA EXAMPLE
} Ordinal
4 classes
3 classes
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. (1 = fit = eigenvalue with
variable j only; a = additional fit = increase in eigenvalue; cum (a) = cumulative total of eigenvalues a; 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
out by taking the month class variables as covariables.
Marginal effects (forward: step 1)
Conditional effects (forward: continued)
j
Variable
1
P
j
Variable
a
P
Cum (a)
1
Shrubs (1/0)
0.25
(0.01)
1
Shrubs (1/0)
0.25
(0.01)
0.25
2
Source distance
0.22
(0.01)
2
Source distance
0.19
(0.01)
0.44
3
EC
0.20
(0.01)
3
EC
0.19
(0.01)
0.63
4
Discharge
0.17
(0.01)
4
Discharge
0.14
(0.03)
0.75
5
Total veg cover
0.16
(0.01)
6
Shading
0.15
(0.01)
-
Cover emergent
0.11
(0.10)
-
7
Soil grain size
0.14
(0.02)
-
Cover bank veg
0.11
(0. 12) -
8
Stream width
0.14
(0.05)
-
Soil grain size
0.10
(0.13)
-
9
High weedy veg
0.14
(0.08)
10
Cover bank veg
0.13
(0.11)
-
U vs L stream
0.22
(0.01)
-
U vs L stream
0.09
(0.01)
-
EXTRA FIT
Each variable is the only env. var.
Change in eigenvalue if particular variable selected
MARGINAL EFFECTS i.e. ignoring all
other variables
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
Outline of ordination techniques presented here. DCA (detrended correspondence
analysis) was applied for the determinaGradient length
tion of the length of the gradient (LG).
estimation
LG is important for choosing between
ordination based on a linear or on an
unimodal response model. Correspondence analysis (CA) is not considered any
CA
further because in “microcosm experiIndirectly
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
CCA
on linear response models is considered
Directly
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 environmental 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
PCA RDA -
Rao (1964)
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 represented 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
-1
-5
Moisture 15
9
Use
5
-6
Manure
-8
16
SF
BF
-10
0
HF
-10
-2
NM
-4
-13
Correlations
Axis1
Axis2
54
-6
92
12
15
29
-26
86
25
76
-48
-11
-40
13
51
-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.
Axis 1 Axis 2
PCA
%
29
21
RDA
%
26
17
PCA
Correlation
0.90
0.82
RDA
Correlation
0.95
0.89
BIPLOT INTERPRETATION
Cosine of angle  correlation
Long arrows of species and environmental variables most important
Goodness of fit
1 + 2
sum of
eigenvalues
unconstrained
species
constrained
fitted species
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 (q env variables, m species)
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)
Primary and secondary data tables in a typical community ecological study of speciesenvironment 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
abundances
sites
Euclidean distances
species
(b) redundancy analysis
species x sitesb
fitted abundances
sitesb
Euclidean distancesc
species
Quantitative env. vars.:
species x env. vars. d
correlations
sites x env. vars. d
env. vars. d
effectse
Qualitative env. vars:
species x env. classesf
means
g
sites x env. classesf
env. classesf
Euclidean distances
env. vars. x env. classes
-
2: focus on species
correlation biplot
abundances
correlationsa
fitted abundances
correlationsa,c
correlations
values of env. vars
correlations
means
g
means
REF
REF
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.
a
Site scores are a linear combination of the environment variables instead of being a
weighted sum of species abundances.
b
c
In the definition of this coefficient, abundance must be replaced by the fitted abundance.
Environmental scores are intraset correlations in scaling 2 and s½ times those in scaling 1
with s the eigenvalue of axis . In CANOCO, the scores are termed biplot scores for
environmental variables.
d
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.
e
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 abundances. Quantitative environment 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)
Optima (WA) of species (CCA)
Species x environmental classes
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
Euclidean distances (PCA, RDA)
Chi-square distances (CA, CCA)
Env. classes x env. classes
Euclidean distances (RDA)
Chi-square distances (CCA)
Env. variables x env. variables
Correlations (RDA, CCA)
Sites x env. variables
Values (RDA, CCA)
Sites x env. classes
Means (RDA, CCA)
Env. variables x env. classes
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)
PARTIAL CONSTRAINED ORDINATIONS
(Partial CCA, RDA, etc)
e.g.
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,
O2 = oxygen, P = ortho-phosphate, Si= siliciumcompunds, 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)
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:
1)
Differences between observers.
2)
Time of observation.
3)
Between-plot variation when interest is temporal trends within repeatedly sampled
plots.
4)
Uninteresting gradients, e.g. elevation when interest is on grazing effects.
5)
Temporal or spatial dependence, e.g. stratigraphical depth, transect position, x and y
co-ordinates. Help remove autocorrelation and make objects more independent.
6)
Collecting habitat – outflow, shore, lake centre.
7)
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

total SS = regression SS + residual SS
Two-way ANOVA

between group (factor 1)
+ between treatments (factor 2)
+ interactions
+ error component
Borcard et al. (1992)
Ecology 73, 1045–1055
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, assemblagelevel, 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
57%
57%
43%
What of explained variance component?
Soil variables (pH, Ca, LOI)
Land-use variables (e.g. grazing, mowing)
Not independent
Do CCA/RDA using
Soil variables only
canonical eigenvalues
Land-use variables only
canonical eigenvalues
Partial analysis Soil
Land-use covariables
Partial analysis Land-use
Soil covariables
Soil variation independent of land-use (3)
0.160
Land-use structured (covarying) soil variation (1–3)
0.361
Land-use independent of soil (4)
0.142
Total explained variance
d) Unexplained
1)
2)
3)
4)
a)
b)
c)
a
unique
b
c
unique
covariance
0.521
0.503
0.160
0.142
13.7%
31%
12.2%
56.9%
43.1%
d
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 –
Climate (C), Geography (G) and Deposition of Pollutants (D)
Series of RDA and partial RDA
Predictors
Covariables
G+C+D
D
G+C
G+C
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.811
0.679
0.106
0.812
0.706
0.074
0.811
0.737
Sum of canonical 
0.034
0.811
0.777
0.228
0.811
0.538
Predictors
Covariables
G
D+G
D+C
D+C
G
G
Joint effect
GD+C=0.777-0.228=0.538-0.034=0.549
Canonical eigenvalues
All predictors
0.811
Pure deposition
0.027
PD
Pure climate
0.106
PC
Pure geography
0.034
PG
Joint G + C
0.132
Joint G + D
0.074
Joint D + C
0.228
Unexplained variance 1 – 0.811 = 0.189
D
Covariance
terms
PD
CD
CD
DG
CG
CDG
DG
CDG
PC
C
CG
PG
G
CD + DG + CDG =
0.652
CD + CG + CDG =
0.631
DG + CG + CDG =
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
CD = 0.095
CDG
PG
DG = 0.013
(G+C)
(CD +
DG +
CDG)
CG = –0.008
= 0.652 – 0.013 – 0.095
= 0.544
= 0.631 – (–0.008) – 0.095
= 0.544
= 0.054 – (–0.008) – 0.013
= 0.544
Total explained variance 0.811 consists of:
Common
Common
Common
Common
climate + deposition
deposition + geography
climate + geography
climate + geography + deposition
0.095
0.013
0.008
0.544
Unique climate PC
Unique geography PG
Unique deposition PD
Unexplained variance
0.106
0.034
0.027
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'.
'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 regressionbased techniques, including (partial) CCA or RDA.
Careful model selection (minimal adequate model) is essential
for many purposes, including variance partitioning.
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
Regression
Multiple regression
Gaussian regression
Weighted averaging of site
scores
Calibration
Linear calibration
'inverse regression'
Gaussian
calibration
Weighted averaging of species
scores (WA)
Ordination
Principal
components
analysis (PCA)
Correspondence analysis (CA);
Gaussian ordination detrended correspondence
analysis (DCA)
Constrained ordination1
Redundancy
analysis (RDA)4
Gaussian canonical
ordination
Canonical correspondence
analysis (CCA); detrended CCA
Partial ordination2
Partial component
analysis
Partial Gaussian
ordination
Partial correspondence
analysis; partial DCA
Partial constrained
ordination3
Partial redundancy
analysis
Partial Gaussian
canonical
ordination
Partial canonical
correspondence analysis;
partial detrended CCA
1=
Unimodal Weighted Averaging
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
0
1
0
1
0
0
Basalt Gabbro
0
0
0
0
0
0
0
0
1
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 (pH2) (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 CONSTRAINED 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 speciesenvironmental 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 (cj) are proportional to their
VIF.
number of predictors
Var c j   VIF residualvariance/n  q  1
number of samples
VIF is related to the (partial) multiple correlation coefficient Rj between variable j and
the other environmental variables.
VIF 
1
1  R 2j
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.
Useful for finding minimal set of variables.
Not unique, e.g. pH and
Ca (and other variables).
VIF
VIF
pH
123.8
6.4
Ca
2.6
–
Mg
8.1
29.3
'AIC' FOR MODEL SELECTION IN CCA AND RDA
Jari Oksanen (2004) VEGAN 1.7-6
R
deviance.cca
deviance.rda
Find statistics in CCA and RDA that resemble deviance and assess an AIClike 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.
AIC -
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 L ˆ  2k
where k = number of free parameters in the model
Lˆ = model maximum likelihood
From estimated residual sum of squares (RSS) in CCA
h
( sum of all eigenvalues   tracei )
i 1
where h = number of predictor variables in model,
 
log L(ˆ)   12 n log( ˆ 2 )
where
log = loge
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
 2K(K  1) 
AICc  AIC  
 n  K  1 
As in GLM, interested in differences in AIC between models
i = AICci – min AICc
Data
– 5 cruises (DEM-1 – DEM-5)
- 8 models
Godinez-Dominguez & Freire, 2003
CCA permutation tests
27 p < 0.05 (global test)
40 models
20 p < 0.05 (first CCA axis)
AICc 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.
CANONICAL CORRELATION ANALYSIS – CANCOR
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
x  b y 
1
k
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. m must be < n – q
CANCOR biplot
Differs from RDA also in error component.
van der Meer (1991) J. Expl. Mar. Biol. Ecol. 148, 105–120
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

no effect
CCA, RDA affects results

no effect
CANCOR

no effect
no effect
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
2
d 'ij  (dij  2c1 )0.5 for i  j
where c1 is equal to absolute value of largest negative eigenvalue of matrix used in PCoA 1
aij  
1 2
dij
2
D
ij  aij  ai  a j  a..
1
Use all principal co-ordinate sample scores (n - 1 or m, whichever is less) as RESPONSE
(species) data in RDA. Use dummy variables for experimental design as predictors in X in
RDA.
Now under framework of RDA and battery of permutation tests, can analyse structured
experiments but WHOLE ASSEMBLAGE (cf. MANOVA but where m >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)
Test of one factor
in a single-factor model
Matrix Y
(replicates x
principal coordinates)
Redundancy analysis (RDA)
F# statistic
Test of interaction term
in multifactorial model
Matrix Y
(replicates x
principal coordinates)
Partial redundancy analysis
(partial RDA) F# statistic
Correction for
negative eigenvalues
Matrix X
(dummy variables
for the factor)
Test of F# by permutation
Matrix X
(dummy
variables
for the
interaction)
Matrix XC
(dummy
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
Multivariate RDA statistic
Total sum of squares
sum of all eigenvalues of Y
Treatment sum of squares = SSTr
trace = sum of all canonical eigenvalues of Y on X
Treatment degrees of freedom = dfTr
q
Residual sum of squares = SSRes
rss = sum of all eigenvalues – trace
Residual degrees of freedom = dfRes
nT – q – 1
Treatment mean square = SStr/dfTr = MSTr
trace/q
Residual mean square SSRes/dfRes = MSRes
rss/(nT – q - 1)
F
MSTr
MSRes
F# 
trace/ q
rss /(nT  q  1)
Tr = treatment, q = number of linearly independent dummy variables,
nT = 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
Y = ZB + XC + E
random error
responses
covariables
B & C - unknown but fixed regression coefficients
(Note Z = covariables and X = predictors here)
To test H0 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 Er of regression Y on Z (covariables)
REDUCED MODEL OR NULL MODEL
CANOCO 3 & 4
4. Permute residuals Ef of regression Y on Z and X (covariables and predictors)
FULL MODEL
CANOCO 3 & 4
1&2
1
2
3&4
DESIGN-BASED PERMUTATIONS
Wrong type I error, low power
OK but no basis for testing interaction effects
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 on X and Z. Full model. Gives lower type II error, but only
slightly so.
(If no covariables Y = XC + E, does not matter if samples in Y or X are permuted. CANOCO
permutes X)
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- RDA ANO- MANOVA
SIM TEL
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
*
*
*
*
*
*
*
*
*
*
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. X may contain the codes of an ANOVA
model (design matrix) or it may contain one or more predictor variables (e.g.,
environmental variable) of interest.
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 X on Y on the basis of any distance
measure of choice and uses permutations of the observations to assess
statistical significance.
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 Qt, number of eigenvectors. Select the
lowest misclassification error (in the case of groups) or the minimum residual sum of
squares (in the case of quantitative variables in X). Also Qt must not exceed m or n and is
chosen so that the proportion of variance explained by the first Qt 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 Qt axis model (trace and first
eigenvalue).
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)
Distance
measure
preserved
Relationship of ordination axes
with original variables
Takes into account
correlation structure
Redundancy Analysis (RDA)
Euclidean
distance
Linear with X, linear with
^
fitted values, Y
= X(X'X)-1 X'Y
... among variables in X,
but not among variables in
Y
Canonical Correspondence
Analysis (CCA)
Chi-square
distance
Linear with X, approx unimodal
with ^Y, linear with fitted
values, Y*
... among variables in X,
but not among variables in
Y
Canonical Correlation Analysis
(CCorA, COR)
Mahalanobis
distance
Linear with X, linear with Y
... among variables in X,
and among variables in Y
Canonical Discriminant
Analysis (CDA; Canonical
Variate Analysis CVA;
Discriminant Function
Analysis, DFA)
Mahalanobis
distance
Linear with X, linear with Y
... among variables in X,
and among variables in Y
Canonical Analysis of Principal
Coordinates (CAP; Generalized
Discriminant Analysis)
Any chosen
distance or
dissimilarity
Linear with X, linear with Qt;
unknown with Y (depends on
distance measure)
... among variables in X,
and among principal
coordinates Qt
CRITERION FOR DRAWING ORDINATION AXES
RDA
•
Finds axis of maximum correlation between Y and some linear
combination of variables in X (i.e., multivariate regression of Y on
^
X, followed by PCA on fitted values, Y).
CCA
•
Same as RDA, but Y are transformed to Y* and weights (square roots
of row sums) are used in multiple regression.
CCorA
•
Finds linear combination of variables in Y and X that are maximally
correlated with one another.
CVA
•
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.
CAP
•
Finds linear combination of axes in Qt and in X that are maximally
correlated, or (if X 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
Yotk + bk cdt + 
Yd(i)tk
=
Yd(i)tk
= abundance counts of taxon k at time t in replicate i of treatment d
Yotk
= mean abundance of taxon k in controls (o) at time t
cdt
= principal response of treatment d at time t (PRC)
bk
= weight of species k with respect to cdt

= error term with mean of zero and variance 2k
d(i)tk
d(i)tk
where
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 replicated), 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 k)
allow species interpretation. Higher the weight, more the actual species response is likely
to follow the PRC pattern, because the response pattern = bk cdt. 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 phytoplankton 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 (bk) can be
interpreted as the affinity of
the taxon with the principal
response curves.
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 betweenmonth variation; this is displayed on the horizontal axis. 57% of all variance could be allocated to betweensite 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 (bk) can be interpreted as the affinity of the variables with the
Principal Response Curves (cdt).
See biggest differences for the two upstream sites, with lower NO 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
in the upstream sites due to pollution.
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.
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
Ym and Xm  Modelm
where subscript m = modern
'Palaeo' CCA/RDA – modelling and reconstruction
Ym and Xm  Modelm

Yo and Modelm  Xo
where subscript o = fossil or
'palaeo' data
(Lecture 8 on Environmental Reconstruction)
Predictive CCA/RDA – modelling and prediction
Ym and Xm  Modelm
Xf and Modelm  Yf
where subscript f = future
(predicted) data
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).
Model performance
Axis 1
Axis 2
Species variance
CA
0.41
0.21
34.3%
CCA
0.28
0.11
21.8%
Total inertia 1.79
Constrained inertia 0.68
38% explained by topography
Species fell into five groups of kappa and other performance values
Kappa > 0.5
e.g. Carex curvula, Veronica alpina
Kappa > 0.4
e.g. Androsace alpina, Oreochloa disticha
Kappa > 0.3
e.g. Saxifraga oppositifolia, Primula glutinosa
Kappa > 0.25
e.g. Ranunculus glacialis, Cerastium uniflorum
Kappa < 0.1
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
SOFTWARE FOR CONSTRAINED ORDINATIONS
CANOCO + CANODRAW
[CANCOR, CAP]
canonical
correlation analysis
constrained principal
co-ordinates analysis
DISTPCOA
distance-based redundancy
analysis via principal coordinates analysis
R(VEGAN) (CCA, RDA)
R (Non-linear CAP)
CANODRAW 4
Pie symbols plot
CANODRAW 4
CANODRAW 4
Response curves fitted
using GAM
Attribute plot
Sample diagram with
principal response curves
Biplot with environmental
variables & sites
T-value biplot
Isolines in RDA
ordination diagram
1987
1987
2002
2003
Mark Hill
Cajo ter Braak
Marti Anderson
Petr Šmilauer