Is there a way to perform a nested/constrained NMDS ordination with metaMDS()? For example, I'd like to constrain my distance matrix with the metadata variables "location" and "Glacier" and I thought the code would look something like this:
ord <- metaMDS(dist ~ location*Glacier, data=metadata)
Where dist is my output from vegdist() on a standard ASV table.
As far as I know, constrained NMDS has not been invented yet (but see the second chapter in this message). Moreover, it wouldn't do what you expect it to do. Constrained metric ordination (CCA, RDA, dbRDA) work by first fitting the linear model of type y ~ location*Glacier and then submitting the fit[ted] values to ordination. So in these models all sampling units with the same location & Glacier will be identical, and will fall in the same point over each other in ordination. We get over this in CCA, RDA, dbRDA using the so-called WA scores which are the predictions for these single locations from community composition. We call these weighted average or WA scores (although for RDA & dbRDA they rather are weighted sum scores). Along these lines, you can have constrained ordination with first finding the averaged (or fit[ted]) values of your community (such as fcomm <- fitted(lm(comm ~ location*Glacier)), or fcomm <- fitted(gam(comm ~ ...) for continuous variables) and using these as input to NMDS. However, getting the scatter of original observations cannot be easily achieved (although the development version of vegan in github has a non-merged branch for adding new points in NMDS ordination).
If you really are interested, the extreme-vegan package natto in https://github.com/jarioksa/natto has function distconstrain that can directly calculate constrained dissimilarities as in vegan:::dbrda (although for non-metric and semi-metric dissimilarities these may fail). Submitting these to NMDS will give constrained NMDS with same problems as the approach outlined above: identical combinations of constraints will give identical dissimilarities. Perhaps more interestingly, you can use these to get partial NMDS (also a method that does not exist, but we have used that in one publication), or NMDS after removing the effect of these variables.