Can I correlate the 'species' with the assemblage variation?

[tjquimpo]

I'm using capscale to discern variability in my three study sites, with benthic cover as the environmental variables. I have the already plotted the 'sites' in an ordination plot, and I have added vector strengths of my environmental variables.

My question is, can I also correlate 'species' with the assemblage variation? This is what I mean for correlating species with assemblage variation (kindly see the image from Lindfield et al. 2016: mesophotic depths as refuge areas for fisher-targeted species on coral reefs).

In this paper, they used PRIMER v6 to correlate species with the assemlage variation. I have looked for a similar procedure in vegan, but I could not find any topic on this.

Any help would be greatly appreciated.

Thank you, Tim

[jarioksa]

People do so. Actually, we do so even in capscale if you call the function with the community data and calculate the dissimilarities within capscale.

Can you do so? I don't know. The correlations of species are correct if you have used Euclidean distances. In that case your results will be identical with rda. The correlations are calculated like the dissimilarities were Euclidean, and if they are not, the correlation are not exactly correct. This may not be absolutely wrong either -- it depends on the kind of dissimilarities you have used. People do so, though.

We are planning to phase out capscale and replace it with dbrda. In dbrda we do not automatically estimate those correlations, but in capscale we do -- but that is going out.

[tjquimpo]

Thank you for your reply Dr. Oksanen. I'm wondering if you could possibly share a script wherein I could replicate the above graph? I looked at your two papers on vegan ordination methods (Oksanen 2012: constrained ordination tutorial with R; Oksanen 2015: multivariate analysis of ecological communities in R: vegan tutorial), and was sadly unable to find this topic.

I wrote this reproducible example using your data, and was able to generate the attached file. I just need to fit the species that contributed to the dissimilarity (after which the graph will be similar to Lindfield et al. 2016).

# load data

data(varechem)
data(varespec)

# subset only the elements
elem = subset(varechem, select = c(N,P,K,Ca,Mg,S,Al,Fe,Mn,Zn,Mo))

# run dbrada

dbrda = capscale(varespec ~ N + P + K + Ca + Mg + S + Al + Fe + Mn + Zn +
           Mo, elem, distance = 'bray')

# vector fit

vf = envfit(dbrda, elem, permutations = 999)

# isolate significant values

elem2 = subset(elem, select = -c(K,S,Zn,Mo))
vf = envfit(dbrda, elem2, permutations = 999)

# plot dbrda (assuming this is the best model)
plot(dbrda, type='n')
points(dbrda, pch = 16, col = 'blue')
plot(vf, col = 'red')

[jarioksa]

Starting from your # vector fit comment:

## you should fit models to the LC scores, hence display = "lc"
vf <- envfit(dbrda, elem, permutations = 999, display="lc")
sf <- envfit(dbrda, varespec, display = "lc") # species correlations
plot(dbrda, type = "n")
## automatic selection of items with the p.max argument
plot(vf, col = "red", p.max = 0.05)
plot(sf, col = "blue", p.max = 0.05)

This adds species as correlations and hence shows relative changes. The standard output showed the absolute changes and hence abundant species were further away from the origin than rare species; with correlations you do not have such a difference. To display the species points as correlations, you can set argument correlation = TRUE in the standard plot, and it seems that you need to set manually the scaling constant to really see the species as our heuristic fails here ( @gavinsimpson : I think we need to tweak the setting of const when we have correlation = TRUE):

plot(dbrda, correlation = TRUE, const = c(100,1))

This standard plot should be consistent with your design version.

[tjquimpo]

Thank you for the answer Dr. Oksanen. I tried it out and this is exactly what I was looking for.