0

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

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
theforestecologist
  • 4,667
  • 5
  • 54
  • 91
  • I do not quite understand what you need. All information about the fit is in the returned gam object, and that is what you need to handle. However, for 2D model it is a thinplate spline gam which does not have such a simple interpretation of marginal responses. Neither should you want to have such an interpretation as it makes no sense (it is not rotation invariant). For fitted vectors (envfit) it is just a technical tool to help displaying linear trend surface on the plane. The info for axes is not rotation invariant, but the fitted plane is invariant w.r.t. configuration on 2D plane. – Jari Oksanen Jan 19 '23 at 08:24
  • I'd misunderstood what you were looking for - @JariOksanen is quite right when he says you shouldn't be trying to decompose this into per axis components with NMDS. You could do that with PCA, CA, etc. but I thought you want to get values of the surface (i.e. predict from the GAM at a grid of points over the ordination). You can use `isotropic = FALSE` to allow for a more marginal interpretation (it would then use a tensor product smooth) but even with a TPRS you still evaluate the function over axis 1 for fixed value of axis 2. BUt really, don't do this for an NMDS – Gavin Simpson Jan 20 '23 at 10:15
  • Hmm. Thanks to both of you for this feedback. My follow-up: *is* there an appropriate way to achieve my goal? Goal: Ito correlate environmental variables with each axis of an NMDS ordination so as to try to explain the distribution of points in NMDS space. – theforestecologist Jan 21 '23 at 01:34

0 Answers0