Построение векторов envfit (веганский пакет) в ggplot2
Я работаю над завершением графика NMDS, который я создал в vegan и ggplot2, но не могу понять, как добавить векторы envfit для загрузки видов в график. Когда я пытаюсь это сказать "неверное графическое состояние".
Приведенный ниже пример немного изменен по сравнению с другим вопросом ( функция Plotting ordiellipse из пакета vegan на график NMDS, созданный в ggplot2), но он выразил именно тот пример, который я хотел включить, так как я использовал этот вопрос, чтобы помочь мне получить metaMDS в ggplot2 в первую очередь:
library(vegan)
library(ggplot2)
data(dune)
# calculate distance for NMDS
NMDS.log<-log(dune+1)
sol <- metaMDS(NMDS.log)
# Create meta data for grouping
MyMeta = data.frame(
sites = c(2,13,4,16,6,1,8,5,17,15,10,11,9,18,3,20,14,19,12,7),
amt = c("hi", "hi", "hi", "md", "lo", "hi", "hi", "lo", "md", "md", "lo",
"lo", "hi", "lo", "hi", "md", "md", "lo", "hi", "lo"),
row.names = "sites")
# plot NMDS using basic plot function and color points by "amt" from MyMeta
plot(sol$points, col = MyMeta$amt)
# same in ggplot2
NMDS = data.frame(MDS1 = sol$points[,1], MDS2 = sol$points[,2])
ggplot(data = NMDS, aes(MDS1, MDS2)) +
geom_point(aes(data = MyMeta, color = MyMeta$amt))
#Add species loadings
vec.sp<-envfit(sol$points, NMDS.log, perm=1000)
plot(vec.sp, p.max=0.1, col="blue")
3 ответа
Начните с добавления библиотек. Дополнительно библиотека grid является необходимым.
library(ggplot2)
library(vegan)
library(grid)
data(dune)
Выполните анализ metaMDS и сохраните результаты во фрейме данных.
NMDS.log<-log(dune+1)
sol <- metaMDS(NMDS.log)
NMDS = data.frame(MDS1 = sol$points[,1], MDS2 = sol$points[,2])
Добавьте загрузки видов и сохраните их как фрейм данных. Направления стрелок косинусов хранятся в списке vectors и матрица arrows, Чтобы получить координаты стрелок, эти значения направления должны быть умножены на квадратный корень из r2 значения, которые хранятся в vectors$r, Более простой способ - использовать функцию scores() как указано в ответе @Gavin Simpson. Затем добавьте новый столбец, содержащий species имена.
vec.sp<-envfit(sol$points, NMDS.log, perm=1000)
vec.sp.df<-as.data.frame(vec.sp$vectors$arrows*sqrt(vec.sp$vectors$r))
vec.sp.df$species<-rownames(vec.sp.df)
Стрелки добавляются с geom_segment() и названия видов с geom_text(), Для обеих задач фрейм данных vec.sp.df используется.
ggplot(data = NMDS, aes(MDS1, MDS2)) +
geom_point(aes(data = MyMeta, color = MyMeta$amt))+
geom_segment(data=vec.sp.df,aes(x=0,xend=MDS1,y=0,yend=MDS2),
arrow = arrow(length = unit(0.5, "cm")),colour="grey",inherit_aes=FALSE) +
geom_text(data=vec.sp.df,aes(x=MDS1,y=MDS2,label=species),size=5)+
coord_fixed()
Проблема с (в остальном превосходным) принятым ответом, которая объясняет, почему все векторы имеют одинаковую длину на прилагаемом рисунке [Обратите внимание, что принятый ответ теперь отредактирован для масштабирования стрелок, как я описываю ниже, чтобы избежать путаница для пользователей, сталкивающихся с вопросами и ответами], заключается в том, что то, что хранится в $vectors$arrows компонент объекта, возвращаемого envfit() являются косинусами направленных векторов. Все они имеют единичную длину, поэтому стрелки на графике @Didzis Elferts имеют одинаковую длину. Это отличается от вывода из plot(envfit(sol, NMDS.log))и возникает потому, что мы масштабируем координаты векторной стрелки по корреляции с конфигурацией ординации ("оси"). Таким образом, виды, которые показывают слабую связь с конфигурацией посвящения, получают более короткие стрелки. Масштабирование выполняется путем умножения направляющих косинусов на sqrt(r2) где r2 являются значениями, показанными в таблице вывода на печать. При добавлении векторов к существующему графику vegan также пытается масштабировать набор векторов таким образом, чтобы они заполняли доступное пространство графика при сохранении относительной длины стрелок. Как это делается, обсуждается в разделе " Подробности " ?envfit и требует использования неэкспортированной функции vegan:::ordiArrowMul(result_of_envfit),
Вот полный рабочий пример, который копирует поведение plot.envfit используя ggplot2:
library(vegan)
library(ggplot2)
library(grid)
data(dune)
# calculate distance for NMDS
NMDS.log<-log1p(dune)
set.seed(42)
sol <- metaMDS(NMDS.log)
scrs <- as.data.frame(scores(sol, display = "sites"))
scrs <- cbind(scrs, Group = c("hi","hi","hi","md","lo","hi","hi","lo","md","md",
"lo","lo","hi","lo","hi","md","md","lo","hi","lo"))
set.seed(123)
vf <- envfit(sol, NMDS.log, perm = 999)
Если мы остановимся на этом и посмотрим на vf:
> vf
***VECTORS
NMDS1 NMDS2 r2 Pr(>r)
Belper -0.78061195 -0.62501598 0.1942 0.174
Empnig -0.01315693 0.99991344 0.2501 0.054 .
Junbuf 0.22941001 -0.97332987 0.1397 0.293
Junart 0.99999981 -0.00062172 0.3647 0.022 *
Airpra -0.20995196 0.97771170 0.5376 0.002 **
Elepal 0.98959723 0.14386566 0.6634 0.001 ***
Rumace -0.87985767 -0.47523728 0.0948 0.429
.... <truncated>
Итак r2 данные используются для масштабирования значений в столбцах NMDS1 а также NMDS2, Окончательный сюжет производится с:
spp.scrs <- as.data.frame(scores(vf, display = "vectors"))
spp.scrs <- cbind(spp.scrs, Species = rownames(spp.scrs))
p <- ggplot(scrs) +
geom_point(mapping = aes(x = NMDS1, y = NMDS2, colour = Group)) +
coord_fixed() + ## need aspect ratio of 1!
geom_segment(data = spp.scrs,
aes(x = 0, xend = NMDS1, y = 0, yend = NMDS2),
arrow = arrow(length = unit(0.25, "cm")), colour = "grey") +
geom_text(data = spp.scrs, aes(x = NMDS1, y = NMDS2, label = Species),
size = 3)
Это производит:
Могу ли я добавить что-то поздно?
Envfit предоставляет значения pvalue, и иногда вы хотите просто отобразить значимые параметры (что-то веганское может сделать для вас с p.=0,05 в команде plot). Я изо всех сил пытался сделать это с ggplot2. Вот мое решение, может быть, вы найдете более элегантное?
Начиная с ответа Дидзиса сверху:
ef<-envfit(sol$points, NMDS.log, perm=1000)
ef.df<-as.data.frame(ef$vectors$arrows*sqrt(ef$vectors$r))
ef.df$species<-rownames(ef.df)
#only significant pvalues
#shortcutting ef$vectors
A <- as.list(ef$vectors)
#creating the dataframe
pvals<-as.data.frame(A$pvals)
arrows<-as.data.frame(A$arrows*sqrt(A$r))
C<-cbind(arrows, pvals)
#subset
Cred<-subset(C,pvals<0.05)
Cred <- cbind(Cred, Species = rownames(Cred))
"Cred" теперь может быть реализован в аргументе geom_segment, как обсуждалось выше.
Краткое дополнение: чтобы получить полное представление о plot.envfit функциональность в ggplot2необходимо применить коэффициент, известный как "длина стрелок полностью использует площадь участка". Я не знаю, было ли это намеренно упущено в ответах выше, поскольку это было даже специально упомянуто Гэвином? Просто извлеките требуемый коэффициент масштабирования, используяarrow_factor <- ordiArrowMul(vf) а затем вы можете применить его к обоим столбцам NMDS в spp.scrs или вы можете сделать это вручную, например
arrow_factor <- ordiArrowMul(vf)
spp.scrs <- as.data.frame(scores(vf, display = "vectors")) * arrow_factor
spp.scrs <- cbind(spp.scrs, Species = rownames(spp.scrs), Pvalues = vf$vectors$pvals, R_squared = vf$vectors$r)
# select significance similarly to `plot(vf, p.max = 0.01)`
spp.scrs <- subset(spp.scrs, Pvalues < 0.01)
# you can also add the arrow factor in here (don't do both!)
ggplot(scrs) +
geom_point(mapping = aes(x = NMDS1, y = NMDS2, colour = Group)) +
coord_fixed() + ## need aspect ratio of 1!
geom_segment(data = spp.scrs,
aes(x = 0, xend = NMDS1 * arrow_factor, y = 0, yend = NMDS2 * arrow_factor),
arrow = arrow(length = unit(0.25, "cm")), colour = "grey") +
geom_text(data = spp.scrs, aes(x = NMDS1 * arrow_factor, y = NMDS2 * arrow_factor, label = Species),
size = 3)