Plotting PCA cordinates as composites of environmental parameters into NMDS ordination of species assemblages

Question

I have conducted a PCA and identified that the principal components (PC) are not driven by a single environmental parameter but are affected by several for each PC. I was then advised to retrieve the coordinates of PC1 and PC2 for each site and plot them onto my NMDS of species assemblages to visualize the potential effects of the environmental composites, instead of plotting the environmental parameters individually like usual with envfit from the vegan package in R. I now have an NMDS biplot with PC1 and PC2 as the arrows, which show there is some separation of sites from the environmental composite of PC1. The aim of this is an exploratory analysis to visualize how species assemblages are affected by the environmental parameters.

I understand the reason behind doing this, but I am not certain it makes sense statistically, nor have I been able to find any references to justify this method. I would very much appreciate some help in determining if this is an appropriate method or if there are better alternatives to this. Thank you very much! [![enter image description here][1]][1]

Answer 1

There may be no paper dealing with PCs in vector fitting over ordination. However, there is abundant literature about using PCs to replace a high number of correlated observed variables with a couple of surrogate variables. Opinions diverge, but I think the dominant modern view is: don't do this! Naturally, there are many ways of implementing this, and some ways have a better justification, but just plugging in PCs from exploratory analysis is something I wouldn't do. Ask yourself these simple questions: How do I interpret those fitted factors? Can I do it? Can I explain to myself (as a starter) what they mean? However, it should not be difficult to find someone who disagrees with this message.

A side note: Ordination plots should always have equal aspect ratio. In your plot one unit of axis is more than two times longer on dimension 2 than on dimension 1. This is particularly important in NMDS where all you try to find is the configuration, or distances among points, and these are distorted if you stretch dim 2. If you use ggplot2 you should add + coord_fixed(ratio=1) in your plotting string (and in vegan plot this is taken care automatically).