How do I add a vector onto an NMDS in R?

Question

I am trying to work out if Salinity has any effect on invertebrate community structure. I've managed to make my NMDS but every time I try ordiplot or envfit there are errors with unsuitable data type. Is there another vector command I could use?

Here is my data

dput(Invertebrates)

structure(list(X = structure(c(9L, 5L, 6L, 13L, 22L, 14L, 7L, 
8L, 19L, 21L, 11L, 23L, 12L, 1L, 17L, 10L, 3L, 4L, 15L, 2L, 20L, 
18L, 16L), .Label = c("Adelaide NS A", "Adelaide NS B", "Adelaide SS A", 
"Adelaide SS B", "Betteshanger Pond A", "Betteshanger Pond B", 
"Broad dike A", "Broad dike B", "Finglesham Brook A", "Finglesham Brook B", 
"Fowlmead Lake A", "Fowlmead lake B", "Great Mongeham A", "Great Mongeham B", 
"Ham Fen SS", "Little Downs Bridge A", "Little Downs Bridge B", 
"S3 Broad dike SS A", "S3 Broad dike SS B", "Site 6 NS A", "Site 6 NS B", 
"Site 7 SS A", "Site 7 SS B"), class = "factor"), Gammarus.pulex = c(112L, 
0L, 0L, 7L, 6L, 32L, 0L, 0L, 14L, 65L, 0L, 0L, 0L, 5L, 48L, 78L, 
8L, 4L, 1L, 3L, 58L, 24L, 10L), Ilyocoris.cimicoides = c(1L, 
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 8L, 0L, 0L, 
0L, 0L, 0L, 0L, 0L, 16L), Bithynia.tentaculata = c(3L, 0L, 0L, 
0L, 0L, 0L, 0L, 0L, 7L, 0L, 0L, 1L, 0L, 3L, 0L, 3L, 0L, 20L, 
0L, 0L, 0L, 0L, 23L), Asellus.aquaticus = c(0L, 0L, 0L, 2L, 0L, 
0L, 0L, 2L, 6L, 2L, 0L, 0L, 0L, 6L, 23L, 15L, 33L, 9L, 6L, 8L, 
2L, 50L, 46L), Sialis.lutaria = c(0L, 0L, 0L, 2L, 0L, 1L, 0L, 
0L, 0L, 2L, 0L, 0L, 0L, 2L, 0L, 1L, 9L, 2L, 3L, 0L, 0L, 0L, 0L
), Haliplus.fluviatilis = c(0L, 0L, 0L, 0L, 6L, 0L, 2L, 0L, 0L, 
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 6L, 0L, 0L, 0L), Coenagrion.pulchellum = c(0L, 
0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 2L, 0L, 2L, 0L, 0L, 
0L, 0L, 1L, 1L, 3L, 2L), Physa.fontilnalis = c(0L, 0L, 0L, 0L, 
0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 
0L, 13L, 0L), Anax.parthenope = c(0L, 0L, 0L, 0L, 0L, 0L, 0L, 
0L, 0L, 0L, 3L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L
), Corixa.punctata = c(0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 
1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 2L, 0L, 0L, 4L), Lymnaea.stagnalis = c(0L, 
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 2L, 0L, 0L, 0L, 
7L, 0L, 0L, 0L, 0L, 0L), Bithynia.leachii = c(0L, 0L, 0L, 0L, 
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 6L, 18L, 0L, 12L, 0L, 0L, 
0L, 0L, 0L, 0L), Lymnaea.truncatula = c(0L, 0L, 0L, 0L, 0L, 0L, 
0L, 0L, 0L, 0L, 0L, 0L, 0L, 4L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 
0L), Radix.palustris = c(0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 
0L, 0L, 0L, 0L, 1L, 2L, 0L, 4L, 0L, 0L, 0L, 0L, 0L, 0L), Bathyomphalus.contortus = c(0L, 
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 19L, 
14L, 0L, 0L, 0L, 0L, 4L), Gyraulus.albus = c(0L, 0L, 0L, 0L, 
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 6L, 1L, 0L, 7L, 0L, 0L, 0L, 
0L, 0L, 0L), Planorbis.planorbis = c(0L, 0L, 0L, 0L, 0L, 0L, 
0L, 0L, 0L, 0L, 0L, 0L, 0L, 4L, 0L, 0L, 4L, 1L, 0L, 12L, 0L, 
0L, 5L), Piscicola.geometra = c(0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 
0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L), 
    Dytiscus.marginalis = c(0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 
    0L, 0L, 0L, 0L, 0L, 0L, 1L, 2L, 0L, 0L, 1L, 0L, 0L, 0L), 
    Haplotaxis.gordioides = c(0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 
    0L, 0L, 0L, 0L, 0L, 0L, 0L, 2L, 0L, 0L, 0L, 0L, 0L, 0L, 0L
    ), Anisus.vortex = c(0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 
    0L, 0L, 0L, 0L, 0L, 0L, 0L, 124L, 29L, 0L, 8L, 0L, 0L, 0L
    ), Planorbis.cornatus = c(0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 
    0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 4L, 0L, 0L, 0L, 0L, 0L, 0L
    ), Radix.ovata = c(0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 
    0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 23L, 52L, 0L, 0L, 0L), Salinity = c(1, 
    4, 4, 4.5, 3, 4.5, 4, 4, 8, 3, 1, 3, 1, 2.5, 2, 1, 6, 6, 
    1, 2.5, 3, 8, 2)), class = "data.frame", row.names = c(NA, 
-23L))

community matrix

com<-Invertebrates[,2:24]
m_com<-as.matrix(com)

NMDS

nmds=metaMDS(m_com,distance="euclidean")
    data.scores=as.data.frame(scores(nmds))
    xx = ggplot(data.scores, aes(x = NMDS1, y = NMDS2)) + 
    geom_point(size = 4, aes( shape = Invertebrates[,1]))+ 
    theme(axis.text.y = element_text(colour = "black", size = 12, face = "bold"), 
    axis.text.x = element_text(colour = "black", face = "bold", size = 12), 
    legend.text = element_text(size = 12, face ="bold", colour ="black"), 
    legend.position = "right", axis.title.y = element_text(face = "bold", size = 14), 
    axis.title.x = element_text(face = "bold", size = 14, colour = "black"), 
    legend.title = element_text(size = 14, colour = "black", face = "bold"), 
    panel.background = element_blank(), panel.border = element_rect(colour = "black", fill = NA, size = 1.2),
    legend.key=element_blank()) + 
    labs(x = "NMDS1", colour = Invertebrates[,1], y = "NMDS2", shape = Invertebrates[,1])+
    scale_shape_manual(values=LETTERS[1:23])
    xx

And how do I change the legend title? Tried in all the extra commands but there were errors.

Thanks

Answer 1

Not sure what your original envfit error was, but if you just select the Salinity column and then run it, it should work. You can manually change the legend title in quotes in the labs() argument. You don't need colour there because you never set a colour aesthetic. So you can just change shape to shape = "Site" or whatever you want to call it. Here's the work-through with your data.

library(vegan)
library(ggplot2)
library(dplyr)

# Community Matrix
com <- Invertebrates[,2:24]
m_com <- as.matrix(com)

# NMDS
nmds = metaMDS(m_com, distance = "euclidean", trymax = 100)
data.scores = as.data.frame(scores(nmds))
nmds$stress
stressplot(nmds)
# Note, I set trymax = 100 because at the default 20 it was not converging

# envfit. First make dataframe for environmental variables
env <- select(Invertebrates, Salinity)
ef <- envfit(nmds, env, permu = 999)
ef
# Note that salinity is NOT significantly correlated with community structure

# Default vegan plot with vector
plot(nmds)
plot(ef)

# Edit dataframe
str(Invertebrates)
names(Invertebrates)[1] <- "Site"
Invertebrates$NMDS1 <- scores(nmds)[,1]
Invertebrates$NMDS2 <- scores(nmds)[,2]

# Make dataframe with vector to add to ggplot
# See https://stackoverflow.com/questions/14711470/plotting-envfit-vectors-vegan-package-in-ggplot2 for more discussion
vec.df <- as.data.frame(ef$vectors$arrows*sqrt(ef$vectors$r))
vec.df$variables <- rownames(vec.df)
# Note this is set up to work if you have more than one environmental variable


# ggplot
xx = ggplot(Invertebrates, aes(x = NMDS1, y = NMDS2)) + 
  geom_point(size = 4, aes(shape = Site)) + 
  geom_segment(data = vec.df,
               aes(x = 0, xend = NMDS1, y = 0, yend = NMDS2),
               arrow = arrow(length = unit(0.5, "cm")),
               colour="blue",
               inherit.aes = FALSE) + 
  geom_text(data = vec.df,
            aes(x = NMDS1, y=NMDS2, label = variables),
            size=4) +
  labs(x = "NMDS1", y = "NMDS2", shape = "Site") +
  scale_shape_manual(values=LETTERS[1:23]) +
  theme(axis.text = element_text(colour = "black", face = "bold", size = 12),
        axis.title = element_text(face = "bold", size = 14, colour = "black"), 
        legend.text = element_text(size = 8, face ="bold", colour ="black"), 
        legend.position = "right", axis.title.y = element_text(face = "bold", size = 14), 
        legend.title = element_text(size = 14, colour = "black", face = "bold"),
        legend.key=element_blank(),
        panel.background = element_blank(), 
        panel.border = element_rect(colour = "black", fill = NA, size = 1.2))
xx
# Note that the vector in the default vegan plot is longer because it goes through an additional scale to fill the graphical space while here it is just based on the correlation, which there was none, so it is a very short vector. I would recommend not putting the vector in this case because the correlation is not significant.