RDA with all explanatory variables and their second-degree terms
To run this RDA we need to compute the second-degree terms of all explanatory
variables except dfs. Raw polynomial terms are highly correlated, so it is preferable
to use orthogonal polynomials, which can be computed by function poly() of
package stats. However, when a matrix is provided to it, poly() computes all
polynomial terms including the ones combining the variables, e.g. xy, x
2 y
2 , and so
on. What we need here are only the variables and their respective (orthogonal)
quadratic terms. Our function polyvars() takes care of this.
env.square <- polyvars(env2, degr = 2)
names(env.square)
spe.envsq.rda <- rda(spe.h ~ ., env.square)
R2ad <- RsquareAdj(spe.envsq.rda)$adj.r.squared
spe.envsq.fwd
adjR2thresh = R2ad)
spe.envsq.fwd
envsquare.red <- env.square[, sort(spe.envsq.fwd$order)]
(spe.envsq.fwd.rda <- rda(spe.h ~., envsquare.red))
RsquareAdj(spe.envsq.fwd.rda)
summary(spe.envsq.fwd.rda)
The result is quite different from the one obtained with the RDA based only on
the first-degree variables (see Sect. 6.3.2.6). In that analysis, three variables (ele,
oxy and bod) were selected, yielding a model with an R
2
adj ¼ 0.5401. The
polynomial RDA is far less parsimonious, having retained nine explanatory variables. As a bonus, the adjusted R
2
adj is now 0.7530, an important increase. Actually,
the nine terms retained belong to five different variables. Indeed, the linear and
quadratic terms of four variables have been retained, those of ele, oxy, slo and
amm. This means that among the species that respond to these variables, some of
them do it in a linear way (modelled by the first-degree term) and others by showing
a unimodal response (modelled by the second-degree term). A triplot of this analysis
allows a finer interpretation Fig. 6.9.
# Triplot using lc (model) site scores and scaling 2
triplot.rda(spe.envsq.fwd.rda,
site.sc = "lc",
scaling = 2,
plot.sites = FALSE,
pos.env = 1,
mult.arrow = 0.9,
mult.spe = 0.9,
mar.percent = 0
)
248
6 Canonical Ordination
To run this RDA we need to compute the second-degree terms of all explanatory
variables except dfs. Raw polynomial terms are highly correlated, so it is preferable
to use orthogonal polynomials, which can be computed by function poly() of
package stats. However, when a matrix is provided to it, poly() computes all
polynomial terms including the ones combining the variables, e.g. xy, x
2 y
2 , and so
on. What we need here are only the variables and their respective (orthogonal)
quadratic terms. Our function polyvars() takes care of this.
env.square <- polyvars(env2, degr = 2)
names(env.square)
spe.envsq.rda <- rda(spe.h ~ ., env.square)
R2ad <- RsquareAdj(spe.envsq.rda)$adj.r.squared
spe.envsq.fwd
adjR2thresh = R2ad)
spe.envsq.fwd
envsquare.red <- env.square[, sort(spe.envsq.fwd$order)]
(spe.envsq.fwd.rda <- rda(spe.h ~., envsquare.red))
RsquareAdj(spe.envsq.fwd.rda)
summary(spe.envsq.fwd.rda)
The result is quite different from the one obtained with the RDA based only on
the first-degree variables (see Sect. 6.3.2.6). In that analysis, three variables (ele,
oxy and bod) were selected, yielding a model with an R
2
adj ¼ 0.5401. The
polynomial RDA is far less parsimonious, having retained nine explanatory variables. As a bonus, the adjusted R
2
adj is now 0.7530, an important increase. Actually,
the nine terms retained belong to five different variables. Indeed, the linear and
quadratic terms of four variables have been retained, those of ele, oxy, slo and
amm. This means that among the species that respond to these variables, some of
them do it in a linear way (modelled by the first-degree term) and others by showing
a unimodal response (modelled by the second-degree term). A triplot of this analysis
allows a finer interpretation Fig. 6.9.
# Triplot using lc (model) site scores and scaling 2
triplot.rda(spe.envsq.fwd.rda,
site.sc = "lc",
scaling = 2,
plot.sites = FALSE,
pos.env = 1,
mult.arrow = 0.9,
mult.spe = 0.9,
mar.percent = 0
)
248
6 Canonical Ordination