Is there a way to extract and/or translate summary.ordisurf output to view the relationship with each ordination axis in R?

Question

I have a question regarding plotting and interpreting the relationship between continuous environmental variables to an NMDS ordination of species abundances using R.

• This is a programming follow-up question to a question on Cross-Validated (stats.stackexchange).

In R's vegan package, function envfit() conveniently provides estimated coefficients for drafting a vector in 2-dimensional ordination space. This is convenient because we get a quantitative way to compare a given environmental variable to both axes of the ordination separately.

However, this really only works if the environmental variables have a linear relationship with the axes -- something again discussed in this CV post. When the variables have a non-linear relationship to each ordination axis, another vegan function, ordisurf() can be used.

But instead of providing a direct quantitative relationship of the variable with each axis, ordisurf() instead simply provides summary output for a gam model fit overall.

Question: Is there a way to extract and/or translate the summary.ordisurf output into the context of individual estimated coefficients for each ordination axis separately?

• Vegan coauthor, Gavin Simpson, encouraged me to ask here on SO for an answer.

---

Below is example code (mostly borrowed from Gavin). ordisurf code:

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

## fit NMDS using Bray-Curtis dissimilarity (default)
set.seed(12)
sol <- metaMDS(dune)

## NMDS plot
plot(sol)
## Fit and add the 2d surface
sol.s <- ordisurf(sol ~ A1, data = dune.env, method = "REML", 
                  select = TRUE)
## look at the fitted model
summary(sol.s)

This produces:

> summary(sol.s)

Family: gaussian 
Link function: identity 

Formula:
y ~ s(x1, x2, k = knots[1], bs = bs[1])
<environment: 0x2fb78a0>

Parametric coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   4.8500     0.4105   11.81 9.65e-10 ***
---
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) 1.591      9 0.863  0.0203 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

R-sq.(adj) =   0.29   Deviance explained =   35%
REML score = 41.587  Scale est. = 3.3706    n = 20

whereas envfit() would produce something separating the relationships with both axes (i.e., "NMDS1" and "NMDS2"):

envfit(sol ~ A1, data = dune.env)


***VECTORS

     NMDS1   NMDS2     r2 Pr(>r)  
A1 0.96473 0.26323 0.3649   0.02 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Permutation: free
Number of permutations: 999