CCA species coordinates relative to ordination surface

Question

I want to display species coordinates relative to an environmental variable using ordination surface (ordisurf in vegan package). In order to do this, I guess I have take the scaling of species scores into consideration. plot.cca function in vegan package seems to scale species proportional to eigenvalues (see example below) and leave site scores unscaled:

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

mod <- cca(varespec ~ Al + P + K, varechem)
summary(mod, display = "species")

# Scaling 2 for species and site scores
# * Species are scaled proportional to eigenvalues
# * Sites are unscaled: weighted dispersion equal on all dimensions

plot(mod, type = "n")
points(mod, display = "sites", pch = 19, cex = 0.1)
ordisurf(mod, varechem$Al, add = TRUE)
text(mod, display = "species")

Species scores seem aligned correctly. Now if I specify another scaling:

scal <- "sites" 
plot(mod, type = "n", scaling = scal)
points(mod, display = "sites", pch = 19, cex = 0.1, scalin = scal)
ordisurf(mod, varechem$Al, add = TRUE, scaling = scal)
text(mod, display = "species", scaling = scal, cex = 0.7)

Species scores get higher values than site scores. The ordination surface, however, seems linked to the site scores.

Related questions:

1. Does the default scaling plot species scores such that they are correctly related to the ordination surface model?
2. (Why) can(not) I specify other scaling such as "sites" (row scores) and still get the species scores correctly aligned related to ordination surface if I consequently specify this for all layers (as above)?

Answer 1

The scaling argument scales one set of score to the other. You're just fiddling with the relative positions of one set versus another.

The ordisurf() surface is inherently constrained to the convex hull (or very slightly beyond it) of the site scores — in whatever scaling you choose to apply — because the model it is fitting is:

$$\hat{y}_i = f(x_{1i}, x_{2i})$$

where $x_{j}$ are the site scores on the two axes of the ordination to which you are fitting the surface.

We don't extrapolate at all in ordisurf() (IIRC we chop the grid of points we predict at to create the contours off if they lie beyond the convex hull of the site scores in the 2-d ordination space being fitted.)

There is a new basis in mgcv (as of early 2017) which allows for extrapolation in a reasonable manner — the b-spline basis. It's a bit more complicated to set up as you need to carefully control where beyond the data the penalty should apply etc and then configure the knots accordingly. But you could try that if you want to extrapolate away from the fitted surface. My gut feeling is that it probably isn't worth the effort.

I guess that explains your second question, so as to the first, the relationship between the surface and the species scores is the same as the relationship between the species scores and the site scores in the scaling you have used, plus some extra uncertainty. You'd have to look up what the relationship is as I don't recall exactly what the properties are in the scaling used and would have to look it up myself :-)