How to plot Pcoa/MDS?

Question

I’m trying to plot Pcoa, but the chart is coming out pretty slick.

Before the Pcoa I did PERMANOVA using data of importance of 10 species of coastal vegetation as a function of 5 different treatments (change of salinity and water) and used the Bray-Curtis dissimilarity matrix. With this same data, I am trying to plot Pcoa. As my object is very large, I selected a part of the data. Follow the code:

library(vegan)

meus_dados <- structure(list(Treatment = c("T1", "T1", "T1", "T1", "T1", "T2", 
"T2", "T2", "T2", "T2", "T3", "T3", "T3", "T3", "T3", "T4", "T4", 
"T4", "T4", "T4", "T5", "T5", "T5", "T5", "T5"), Sal = c(20L, 
20L, 20L, 20L, 20L, 3L, 3L, 3L, 3L, 3L, 10L, 10L, 10L, 10L, 10L, 
10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L), Agua = c(6L, 
6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 2L, 2L, 2L, 2L, 2L, 15L, 
15L, 15L, 15L, 15L, 6L, 6L, 6L, 6L, 6L), Sp1 = c(0.794290128748461, 
0.676055337749016, 0.649817348325361, 0.490593956384099, 0.460140400409192, 
0.356605528960497, 0, 0.410011047234125, 0.485048880341384, 0.380487296882943, 
0.478130491326628, 0.393925420118031, 0.509315406411044, 0.548767646671349, 
0.249824697144633, 0.588139655317169, 0.458489317544165, 0.508041166651013, 
0.489174836104582, 0.508196906269527, 0.627466911428099, 0.774643095726598, 
0.543312871998383, 0.430149253292226, 0.488483347062045), Sp2 = c(0.31950289567597, 
0.658313387776009, 0.688008327519027, 0.586911420488643, 0.580577992032451, 
0.356605528960497, 0.239378990887995, 0.225481709273101, 0.242524440170692, 
0.194777553950768, 0.888178597896676, 0.996509961427144, 0.795541435184066, 
0.713158936565669, 0.802719823057339, 0.543860201049459, 0.495818132429415, 
0.455672570231077, 0.36522265070534, 0.368021150346472, 0.774093342394432, 
0, 0.38506500075798, 0.524676787317157, 0.394399669660163), Sp3 = c(0, 
0, 0, 0, 0, 0, 0.0926044346179928, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
0, 0, 0.0913056626763351, 0.145824930692855, 0, 0.127252209323957, 
0.0962662501894949, 0, 0), Sp4 = c(0, 0, 0, 0, 0, 0.257361359186398, 
0.146774556270003, 0.174456464424683, 0.187642521006036, 0.154357559497079, 
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0), Sp5 = c(0.31950289567597, 
0, 0.419594374295962, 0.316782038357181, 0.338654327178036, 0, 
0, 0, 0, 0, 0.35726949591146, 0.304782309227413, 0.440485455199368, 
0.233496302294192, 0.322893736936447, 0, 0, 0, 0, 0, 0.261044554272076, 
0.201689979912391, 0.333026509713188, 0.391728749895399, 0.254877089092083
), Sp6 = c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.214884344687665, 
0.0919770920816606, 0, 0, 0, 0, 0, 0, 0), Sp7 = c(0.31950289567597, 
0.238365530436068, 0, 0, 0, 0.43497144733486, 0.525291741706045, 
0.467329790348523, 0.485048880341384, 0.389555107901536, 0, 0, 
0, 0, 0, 0.5386670892028, 0.406489870695075, 0.350189434253656, 
0.464514166061593, 0.464110931651819, 0, 0.328942189236348, 0.248844646824885, 
0.253151032233496, 0.341099671900618), Sp8 = c(0.247201184223628, 
0.238365530436068, 0.24257994985965, 0.514142039871732, 0.465470460285241, 
0.16237693893029, 0.185208869235986, 0.174456464424683, 0.187642521006036, 
0.194777553950768, 0.276421414865236, 0.304782309227413, 0.254657703205522, 
0.233496302294192, 0.249824697144633, 0.219555369620381, 0.18509129990788, 
0.183954184163321, 0.18261132535267, 0.184010575173236, 0.337395191905393, 
0.567472525800706, 0.393484720516069, 0.400294177261723, 0.521140222285091
), Sp9 = c(0, 0, 0, 0, 0, 0.269702257697169, 0.440323668810008, 
0.410011047234125, 0.41209275713447, 0.588656150841522, 0, 0, 
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0), Sp10 = c(0, 0.188900213602838, 
0, 0, 0.15515682009508, 0, 0, 0, 0, 0, 0, 0, 0, 0.271080812174598, 
0.374737045716949, 0.109777684810191, 0.239227034735801, 0.318188460537611, 
0.315865696423145, 0.237830218279473, 0, 0, 0, 0, 0)), class = "data.frame", row.names = c(NA, 
-25L))

dist <- vegdist(meus_dados[,-c(1:3)], method = "bray")

groups <- meus_dados$Treatment
groups <- as.factor(groups)


pcoa <- cmdscale(dist)
efit <- envfit(pcoa, meus_dados[,2:3])
plot(pcoa, col = c("black", "orange", "pink", "blue", "green")[groups], pch = c(19,1, 24,5,6,7,9)[groups],
     xlim = c(-.3,0.3), ylim=c(-.3, .2),
     xlab = "PCoA 1", ylab = "PCoA 2")
abline(h = 0, v = 0, lty = 2)
plot(efit, col = "red", cex = 0.9)

And this is what I’m getting:

a large mass of colored dots. What I would like to plot the name of the 10 species instead of these points according to the 5 different treatments.

I have tried other methods such as NMDS (with the metaMDS function) for each treatment separately but gives the following error:

no non-missing arguments to min; returning Inf

Would anyone have a suggestion to plot Pcoa the way I’m thinking or any other suggestion that I can do this?

Answer 1

You can extract the position of the species with wascores. The code below plots the graph with the text at the point position. Note that the object efit now includes the three environmental variables (meus_dados[, 1:3]). In your code, also note that the values reported in the arguments xlim and ylim on the first call of the function plot hide sample units to the left of the graph.

pcoa <- cmdscale(dist)
efit <- envfit(pcoa, meus_dados[, 1:3])
ordiplot(pcoa, type = "t", cex = 0.9)
abline(h = 0, v = 0, lty = 3)
species_scores <- wascores(pcoa, meus_dados[, -c(1:3)])
plot(efit, col = "red", cex = 0.9)
text(species_scores, rownames(species_scores), cex = 0.9, col = "red")