Как добавить вектор в NMDS в R? - PullRequest
Новая
       56

Как добавить вектор в NMDS в R?

0 голосов
/ 06 апреля 2020

Я пытаюсь понять, влияет ли соленость на структуру сообщества беспозвоночных. Мне удалось создать NMDS, но каждый раз, когда я пытаюсь выполнить регулярное задание или envfit, возникают ошибки с неподходящим типом данных. Есть ли еще одна векторная команда, которую я мог бы использовать?

Вот мои данные

dput (Беспозвоночные) 

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

А как мне изменить название легенды? Пробовал во всех дополнительных командах, но были ошибки.

Спасибо

1 Ответ

0 голосов
/ 10 апреля 2020

Не знаю, какой была ваша первоначальная ошибка envfit, но если вы просто выберите столбец Salinity и затем запустите его, он должен работать. Вы можете вручную изменить заголовок легенды в кавычках в аргументе labs(). Там вам не нужен цвет, потому что вы никогда не устанавливаете цвет эстетики c. Таким образом, вы можете просто изменить форму на 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.
...