I am trying to group data points of an NMDS in ggplot by adding ellipses using the ordiellipse
function following the helpful advice from this post. However, although I don't get any error messages or warnings, calculating the ellipse data produces an empty data frame.
The dataset is available here, and my code is as following:
library(vegan)
library(ggplot2)
setwd("C:")
veg_matrix <- read.csv("161019_vegetation_matrix.csv",header=T); veg_matrix[2:123] <- lapply(veg_matrix[2:123],as.character); veg_matrix[2:123] <- lapply(veg_matrix[2:123],as.numeric)
rownames(veg_matrix) <- veg_matrix[,1]; veg_matrix <- veg_matrix[,-1] # remove non-numeric rownames (=plot codes) from dataset
speciescomp_NMDS = metaMDS(veg_matrix, k=2, trymax=100, distance="raup", na.rm=T)
plot(speciescomp_NMDS,display="sites",type="n")
# add grouping data
plot_scores <- as.data.frame(scores(speciescomp_NMDS))
plot_scores$plot <- rownames(plot_scores);plot_scores <- cbind(plot_scores,data.frame(matrix(unlist(strsplit(plot_scores$plot,"_")),nrow=24,byrow=T)))[,-3]; colnames(plot_scores)[c(3,4)] <- c("summit","aspect")
plot_scores$group <- ""
plot_scores$group[plot_scores$summit=="BUF"|plot_scores$summit=="SES"] <- "low"
plot_scores$group[plot_scores$summit=="CHA"|plot_scores$summit=="MIN"] <- "intermediate"
plot_scores$group[plot_scores$summit=="CUO"|plot_scores$summit=="GAJ"] <- "high"
plot_scores$group <- transform(plot_scores, as.factor(plot_scores$group))
# compute ellipse data
ord <- ordiellipse(speciescomp_NMDS,plot_scores$group,display = "sites", kind = "sd", conf = .95, label=T)
f_ellipse <- 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_level <- data.frame()
for(g in levels(plot_scores$group)){
if(is.null(ord[[g]])) next
df_ell_level <- rbind(df_ell_level,
cbind(as.data.frame(with(plot_scores[plot_scores$group==g,],
f_ellipse(ord[[g]]$cov,ord[[g]]$center,ord[[g]]$scale)))
,level=group))
}
Any ideas on how to handle this one will be warmly appreciated!
First convert group
to a factor
plot_scores <- transform(plot_scores, group = as.factor(group))
If we now fix your loop to be:
df_ell_level <- data.frame()
for(g in levels(plot_scores$group)){
if(is.null(ord[[g]])) next
df_ell_level <-
rbind(df_ell_level,
cbind(as.data.frame(f_ellipse(ord[[g]]$cov, ord[[g]]$center, ord[[g]]$scale)),
level=g)
)
}
Then I get a data frame with 301 rows. But I've cut out a lot of your unnecessary code from the call, and you really don't want to be rbind()
ing data frames together in a loop like this.