Differences between vegan::rda and calibrate::rda

[marchtaylor]

Hello,

I'm trying to understand the rda scores returned by summary.rda. Specifically, I have noticed that there are differences between the constraining (e.g. environmental) matrix variable scores between vegan::rda and calibrate::rda.

In the example below, I show there differences visually for the first 3 RDA axes scores from both functions. Basically, calibrate seems to do a simple projection of the constraining data onto its site scores in order to derive the environmental variable scores (e.g. summary()$biplot in vegan). I haven't dealt with all of the scaling issues in the comparison, but I'm more concerned when the scores are not linearly related to each other. They are in fact linearly related for all other outputs (species scores, env. site scores, comm. site scores).

If you would be so kind as to help me understand these differences. I am especially interested in the variable scores of both matrices as I need to project these scores back onto another PCA, which was used to derive the RDA input matrices.

Kind regards

Example:

library(vegan)
library(calibrate)
data(varespec)
data(varechem)

RDAveg <- vegan::rda(X = varespec[,1:10], Y = varechem, scale = FALSE) # number of species reduced
RDAcal <- calibrate::rda(Y = varespec[,1:10], X = varechem, scaling = 0)


# extract scores from vegan::rda ------------------------------------------
S <- summary(RDAveg, axes = 5, scaling = 1)
# orthonormal environmental scores (columns of X)
RDAveg$A <- S$biplot
# orthonormal species scores (columns of Y)
RDAveg$B <- S$species
# orthonormal site scores (rows of X)
RDAveg$V <- S$constraints
# orthonormal site scores (rows of Y)
RDAveg$W <- S$sites
 

# extract scores from calibrate::rda ------------------------------------------
# orthonormal environmental scores (columns of X)
RDAcal$A <- RDAcal$Gxs # Xa <- scale(varechem, scale = F); solve(t(Xa) %*% Xa) %*% t(Xa) %*% RDAcal$Fs
# orthonormal species scores (columns of Y)
RDAcal$B <- RDAcal$Gys
# orthonormal site scores (rows of X)
RDAcal$V <- RDAcal$Fs
# orthonormal site scores (rows of Y)
RDAcal$W <- scale(varespec[,1:10], center = T, scale = F) %*% RDAcal$Gys


# plot --------------------------------------------------------------------
op <- par(mfcol = c(2,2), oma = c(0,0,1.5,0), mar = c(3,3,1.5,0.5), mgp = c(2,0.5,0))
plot(c(t(RDAveg$A[,1:3])), c(t(RDAcal$A[,1:3])), col = 1:3, xlab = "vegan", ylab = "calibrate", main = "A: env. scores")
plot(c(t(RDAveg$B[,1:3])), c(t(RDAcal$B[,1:3])), col = 1:3, xlab = "vegan", ylab = "calibrate", main = "B: spp. scores")
plot(c(t(RDAveg$V[,1:3])), c(t(RDAcal$V[,1:3])), col = 1:3, xlab = "vegan", ylab = "calibrate", main = "V: site score (env.)")
plot(c(t(RDAveg$W[,1:3])), c(t(RDAcal$W[,1:3])), col = 1:3, xlab = "vegan", ylab = "calibrate", main = "W: site scores (spp.)")
mtext(text = "Comparison of scaled RDA biplot scores (axes 1:3)", line = 0.5, outer = T, adj = 0, font = 2)
par(op)

Output: