# Matrix of fitted values (eq. 11.10)
Yhat <- Xcr %*% B
# Matrix of residuals
Yres <- Yc - Yhat
## 3. PCA on fitted values
# Covariance matrix (eq. 11.12)
S <- cov(Yhat)
# Eigenvalue decomposition
eigenS <- eigen(S)
# How many canonical axes?
kc <- length(which(eigenS$values > 0.00000001))
# Eigenvalues of canonical axes
ev <- eigenS$values[1 : kc]
# Total variance (inertia) of the centred matrix Yc
trace = sum(diag(cov(Yc)))
# Orthonormal eigenvectors (contributions of response
# variables, scaling 1)
U <- eigenS$vectors[, 1 : kc]
row.names(U) <- colnames(Y)
# Site scores (vegan's wa scores, scaling 1; eq.11.17)
F <- Yc %*% U
row.names(F) <- row.names(Y)
# Site constraints (vegan's 'lc' scores, scaling 1;
# eq. 11.18)
Z <- Yhat %*% U
row.names(Z) <- row.names(Y)
## 2. Computation of the multivariate linear regression
# Matrix of regression coefficients (eq. 11.9)
B <- solve(t(Xcr) %*% Xcr) %*% t(Xcr) %*% Yc
# Species-environment correlations
corXZ <- cor(X, Z)
# Diagonal matrix of weights
D <- diag(sqrt(ev / trace))
# Canonical coefficients (eq. 11.19)
CC <- B %*% U
row.names(CC) <- colnames(X)
254
6 Canonical Ordination