Trazar la función ordiellipse del paquete vegan en el diagrama NMDS creado en ggplot2

Question

En lugar de la función de diagramación normal, estoy usando ggplot2 para crear gráficos ggplot2 . Me gustaría mostrar grupos en el gráfico ordiellipse() usando la función ordiellipse() del paquete vegan .

Ejemplo de datos:

library(vegan)
library(ggplot2)
data(dune)
# calculate distance for NMDS
sol <- metaMDS(dune)
# 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)
# draw dispersion ellipses around data points
ordiellipse(sol, MyMeta$amt, display = "sites", kind = "sd", label = T)
# 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))

¿Cómo puedo agregar ordiellipse al diagrama de ggplot2 creado con ggplot2 ?

La respuesta de Didzis Elferts a continuación funciona muy bien. ¡Gracias! Sin embargo, ahora estoy interesado en trazar el siguiente ordiellipse a la trama ggplot2 creada con ggplot2 :

ordiellipse(sol, MyMeta$amt, display = "sites", kind = "se", conf = 0.95, label = T)

Desafortunadamente, no entiendo lo suficiente sobre cómo funciona la función veganCovEllipse para poder ajustar el script yo mismo.

Answer 1

Antes que nada, agregué un grupo de columnas a su marco de datos NMDS.

 NMDS = data.frame(MDS1 = sol$points[,1], MDS2 = sol$points[,2],group=MyMeta$amt)

El segundo marco de datos contiene los valores medios de MDS1 y MDS2 para cada grupo y se usará para mostrar los nombres de los grupos en el diagrama

 NMDS.mean=aggregate(NMDS[,1:2],list(group=group),mean)

El marco de datos df_ell contiene valores para mostrar elipsis. Se calcula con la función veganCovEllipse que está oculta en el paquete vegan . Esta función se aplica a cada nivel de NMDS (grupo) y utiliza también la función cov.wt para calcular la matriz de covarianza.

 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(NMDS$group)){
df_ell <- rbind(df_ell, cbind(as.data.frame(with(NMDS[NMDS$group==g,],
veganCovEllipse(cov.wt(cbind(MDS1,MDS2),wt=rep(1/length(MDS1),length(MDS1)))$cov,center=c(mean(MDS1),mean(MDS2)))))
,group=g))
}

Ahora las elipses se trazan con la función geom_path() y annotate() usadas para trazar los nombres de los grupos.

 ggplot(data = NMDS, aes(MDS1, MDS2)) + geom_point(aes(color = group)) +
geom_path(data=df_ell, aes(x=MDS1, y=MDS2,colour=group), size=1, linetype=2)+
annotate("text",x=NMDS.mean$MDS1,y=NMDS.mean$MDS2,label=NMDS.mean$group)

La idea para el trazado de la elipse se adoptó a partir de otra pregunta de .

ACTUALIZACIÓN: solución que funciona en ambos casos

Primero, crea un marco de datos NMDS con una columna de grupo.

NMDS = data.frame(MDS1 = sol$points[,1], MDS2 = sol$points[,2],group=MyMeta$amt)

A continuación, guarde el resultado de la función ordiellipse() como algún objeto.

ord<-ordiellipse(sol, MyMeta$amt, display = "sites",
kind = "se", conf = 0.95, label = T)

El marco de datos df_ell contiene valores para mostrar elipsis. Se calcula nuevamente con la función veganCovEllipse que está oculta en el paquete vegan . Esta función se aplica a cada nivel de NMDS (grupo) y ahora usa argumentos almacenados en el objeto ord - cov , center y scale de cada nivel.

df_ell <- data.frame()
for(g in levels(NMDS$group)){
df_ell <- rbind(df_ell, cbind(as.data.frame(with(NMDS[NMDS$group==g,],
veganCovEllipse(ord[[g]]$cov,ord[[g]]$center,ord[[g]]$scale)))
,group=g))
}

El trazado se hace de la misma manera que en el ejemplo anterior. En cuanto al cálculo de las coordenadas para elipses se utiliza el objeto de ordiellipse() , esta solución funcionará con diferentes parámetros que proporcione para esta función.

ggplot(data = NMDS, aes(MDS1, MDS2)) + geom_point(aes(color = group)) +
geom_path(data=df_ell, aes(x=NMDS1, y=NMDS2,colour=group), size=1, linetype=2)