What is ordisurf() doing…?

I’m writing this post for two reasons: i) someone searched on Google for the term ‘what is ordisurf doing’ and ended up on my blog, and ii) because I have been on the receiving end of reviewers comments on a paper I co-authored where they didn’t know what ordisurf() was doing either! It is hardly surprising that people who don’t know R or haven’t studied the code or the examples in the vegan documentation do not realise what ordisurf() is trying to do as there isn’t a paper in the scientific literature explaining the method. Whilst a solution to that part of the problem will have to wait until Jari, Dave Roberts and I get our acts together and write one, this post might be useful in the interim.

Before direct gradient analysis or canonical ordination was invented/used in ecology, the standard approach to analysing multivariate ecological data was to ordinate them using PCA or CA, for example, and then relate the separate, important axes of that ordination with a set of environmental variables. A multiple regression was often used to relate the two, with ordination axis score taken as the response variable and the set of environmental variables as the predictors. This is problematic for methods like nMDS that don’t have axes, where we should consider the k dimensions of the solution as a whole rather than as independent “axes” of variation, which is where ordisurf() comes in.

So we have something tangible to work with in the ensuing discussion, lets fit a simple PCA to the classic Dutch Dune Meadows data set, available in vegan, and display the resulting biplot

require(vegan)
data(dune)
data(dune.env)

dune.pca <- rda(dune)
biplot(dune.pca, scaling = 3)

For these data the main continuous, quantitative variable is A1, the thickness of the soil A1 horizon. The envfit() can be used to project a biplot arrow for this variable into the ordination

set.seed(17)
dune.ev <- envfit(dune.pca ~ A1, data = dune.env)
plot(dune.ev)

The result of the vector fitting is shown below, indicating borderline significance for the A1 horizon

R> dune.ev

***VECTORS

       PC1     PC2     r2 Pr(>r)
A1 0.98316 0.18274 0.2632  0.063 .
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
P values based on 999 permutations.

As this is an unconstrained ordination, there is no reason at all to assume that the values of the A1 horizon vary in a linear fashion across the biplot. Instead, it might be better to fit a smooth response surface of A1 values over the limits of the biplot. For that we use ordisurf().

dune.sf <- ordisurf(dune.pca ~ A1, data = dune.env, plot = FALSE, scaling = 3)
biplot(dune.pca, scaling = 3)
plot(dune.ev)
plot(dune.sf, col = "forestgreen", add = TRUE)

Combining all the various plotting elements thus far we get this figure

The fitted surface is far from linear! The object returned by ordisurf() is an augments object of class “gam” from the mgcv package, so we can use methods from that package to interrogate the result

R> summary(dune.sf)

Family: gaussian
Link function: identity

Formula:
y ~ s(x1, x2, k = knots)
<environment: 0x5df9e50>

Parametric coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept)   4.8500     0.3567    13.6  2.1e-08 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Approximate significance of smooth terms:
           edf Ref.df    F p-value
s(x1,x2) 7.583  8.621 1.94   0.147

R-sq.(adj) =  0.464   Deviance explained = 67.8%
GCV score = 4.4577  Scale est. = 2.5448    n = 20

Here we also see that there is little evidence to reject the null hypothesis. So what is ordisurf() actually doing? It is doing nothing more than fitting the following model using the gam() function

require(mgcv)
scrs <- data.frame(scores(dune.pca, display = "sites", scaling = 3))
dat <- with(dune.env, cbind(scrs, A1))
mod <- gam(A1 ~ s(PC1, PC2, k = 10), data = dat)

Line 1 loads the mgcv package. In line 2 we extract the site scores on PCA axes 1 and 2 using symmetric scaling (scaling = 3) and convert to a data frame. Then we column-bind the soil A1 horizon thickness data into object dat, which we can pass as the data object in the call to gam(). The final line of code (line 4) fits the response surface model using a model formula and a 2D smooth on the PCA axis 1 and 2 sites scores. We restrict the complexity of this smooth using k = 10 as there are only 14 unique values for A1 and the default were starting with a smooth with more degrees of freedom than unique data points.

It should now be clear that we fit a model to predict the soil A1 horizon thickness using a 2-D smooth of the PCA site scores on axes 1 and 2 as the predictor variable. This is backwards to how we might conventionally think of the problem of relating explanatory variables to ordination axes, but it is logical if you think of the model as saying “given the main pattern in the species composition described by ordination axes, how well does this pattern explain variation in response variable at the sites.”

In newer versions of vegan (>= 1.17-9) we now provide access to more of the functionality provided by Simon Wood’s mgcv package for fitting GAMs:

• you can alter the penalty used in the GCV routine that selects the degree of smoothness in the fitted smooth function (the response surface) via argument penalty. A penalty of 1.4 degrees of freedom per knot is often suggested if greater penalty on complex smooths is desired (penalty = 1.4);
• ordisurf() now accepts the method argument of gam(). Simon’s latest advice to me was that doing the smoothness selection via Marginal (Maximum) Likelihood (ML) or Restricted Maximum Likelihood (REML) gave the most reliable p-values on the smooth functions. To use this form of model fitting instead of GCV, supply method = “ML” or method = “REML” in the ordisurf() call;
• By default, gam() can penalize smooths back to linear functions/surfaces but no further. An additional penalty term can be added to the smoothness selection procedure so that smooths can be penalised all the way back to zero degrees of freedom, effectively removing those terms from the model. This is a formal means of model selection. To turn this feature on, add select = TRUE to the ordisurf() call.

As a final illustration, we compare the response surface fitted earlier with one fitted using ML smoothness selection and the extra penalty term

dune.sf2 <- ordisurf(dune.pca ~ A1, data = dune.env, plot = FALSE, scaling = 3,
                     method = "ML", select = TRUE)
biplot(dune.pca, scaling = 3)
plot(dune.ev)
plot(dune.sf2, col = "forestgreen", add = TRUE)

We get a similar pattern to before

but the surface is a lot less complex using approximately 4.5 degrees of freedom compared with approximately 8.7 in dune.sf. This doesn’t alter our interpretation of the significance of the relationship between the plant composition and A1 horizon thickness, however.

R> summary(dune.sf2)

Family: gaussian
Link function: identity

Formula:
y ~ s(x1, x2, k = knots)
<environment: 0x61d67d8>

Parametric coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept)   4.8500     0.3958   12.26 2.24e-09 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Approximate significance of smooth terms:
           edf Ref.df     F p-value
s(x1,x2) 3.519  4.999 1.883   0.156

R-sq.(adj) =   0.34   Deviance explained = 46.3%
ML score = 41.909  Scale est. = 3.1326    n = 20

It remains to be seen whether we can trust the p-values that mgcv provides for predictor data derived from an ordination. Preliminary simulations that Jari Oksanen and I have made suggest the p-values have the right Type I error rate when we use randomly generated data with no relationship to the ordination axes. However, we have only just started this work and those are the sorts of results that are best presented for peer review and not relegated to a blog post. I hope that goes some way to explaining what on Earth it is that ordisurf() does…