В настоящее время у меня возникают проблемы с отображением вывода PERMANOVA в R с использованием примера набора данных, который имеет матрицу сообщества и две переменные метаданных (DIVERSITY и SITE). Выходные данные PERMANOVA предполагают, что DIVERSITY является единственным статистически значимым термином, однако в моем невоспроизводимом наборе данных результаты показывают, что оба термина, плюс термин взаимодействия (DIVERSITY * SITE), являются значимыми.
Я могу построить результаты (см. Ниже), но только с эллипсами, которые очерчивают каждую переменную отдельно. Мне нужна помощь в построении корпуса / эллипса, который четко показывает термин взаимодействия из результатов PERMANOVA (DIVERSITY * SITE). Другими словами, как я могу показать корпус / эллипс, который иллюстрирует разницу между DIVERSITY M на САЙТЕ A и DIVERSITY P на SITE B, и так далее?
Любая помощь будет принята с благодарностью!
Воспроизводимый пример:
library(vegan)
birds<-read.csv('https://raw.githubusercontent.com/collnell/lab-demo/master/bird_by_fg.csv')
birds$SITE=factor(rep(c(rep("SITEA", 4), rep("SITEB", 4), rep("SITEC", 4), rep("SITED", 4))))
bird.matrix<-as.matrix(birds[,3:9])
bird.mat<-sqrt(bird.matrix)
bird.dist<-vegdist(bird.mat, method='bray')
set.seed(36)
bird.div<-adonis2(bird.dist~DIVERSITY*SITE, data=birds, permutations = 999, method="bray", strata="PLOT")
bird.div
NMDS=metaMDS(bird.dist, k=2)
bird.perm=as.data.frame(scores(NMDS))
bird.perm$Site=birds$SITE
bird.perm$Div=birds$DIVERSITY
MDS=data.frame(MDS1=NMDS$points[,1], MDS2=NMDS$points[,2], group=birds$DIVERSITY)
veganCovEllipse<-function (cov, center = c(0, 0), scale = 1, npoints = 100)
{
theta <- (0:npoints) * 2 * pi/npoints
Circle <- cbind(cos(theta), sin(theta))
t(center + scale * t(Circle %*% chol(cov)))
}
df_ell <- data.frame()
for(g in levels(MDS$group)){
df_ell <- rbind(df_ell, cbind(as.data.frame(with(MDS[MDS$group==g,],
veganCovEllipse(cov.wt(cbind(MDS1,MDS2),wt=rep(1/length(MDS1),length(MDS1)))$cov,center=c(mean(MDS1),mean(MDS2)))))
,group=g))
}
ggplot(bird.perm, aes(x=NMDS1, y=NMDS2)) +
geom_point(aes(NMDS1, NMDS2, colour = factor(bird.perm$Site), shape = factor(bird.perm$Div)), size = 2) +
coord_fixed()+theme_classic() +
theme(panel.background = element_rect(fill = NA, colour = "black", size = 1, linetype = "solid")) +
labs(colour = "Site", shape = "Diversity") +
theme(legend.position = "right", legend.text = element_text(size = 12), legend.title = element_text(size = 12), axis.text = element_text(size = 10)) +
geom_path(data=df_ell, aes(x=MDS1, y=MDS2, colour=group), size=1, linetype=2)
Выход NMDS с переменной DIVERSITY, используемой в эллипсах. Я могу отключить переменную DIVERSITY для переменной SITE, но не могу возможность выяснить способ получения термина взаимодействия, DIVERSITY * SITE для рисования эллипсов.