如何将向量添加到 R 中的 NMDS？

Question

我正在尝试确定盐度是否对无脊椎动物群落结构有任何影响。我已经设法制作了我的 NMDS，但每次我尝试 ordiplot 或 envfit 时都会出现数据类型不合适的错误。我可以使用其他矢量命令吗？

这是我的数据

输入（无脊椎动物）

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))

社区矩阵

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

如何更改图例标题？尝试了所有额外的命令，但出现了错误。

谢谢

Answer 1

不确定您最初的 envfit 错误是什么，但如果您只是选择盐度列然后运行它，它应该可以工作。您可以在 labs() 参数中手动更改带引号的图例标题。你不需要颜色，因为你从来没有设置颜色美学。因此，您可以将形状更改为 shape = "Site" 或任何您想调用的名称。这是您的数据的解决方法。

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.