Reputation: 167
This might be a simple one, but I'm trying to draw ellipses around my treatments on my PCoA plot.
My data frame (sc) is:
MDS1 MDS2 Treatment
X1xF1 -0.19736183 -0.24299825 1xFlood
X1xF2 -0.17409568 -0.29727596 1xFlood
X1xF3 -0.15272444 -0.28553837 1xFlood
S1 -0.06643271 0.47049959 Start
S2 -0.15143350 0.31152966 Start
S3 -0.26156297 0.12296849 Start
X3xF1 0.29840827 0.04581617 3xFloods
X3xF2 0.50503749 -0.07011503 3xFloods
X3xF3 0.20016537 -0.05488630 3xFloods
and my code is:
ggplot(data=sc,(aes(x=MDS1,y=MDS2,colour = Treatment)))+geom_point(size=3)+
ggtitle("PCoA of samples at 'class' level(method='Bray')\n",sep=''))+
theme_bw()+guides(colour = guide_legend(override.aes = list(size=3)))+
stat_ellipse()
It plots the PCoA okay up until stat_ellipse(). I've tried it with various parameters and at best I can get one ellipse for the whole plot (although I can't seem to reproduce that now).
What I'm after is three CI ellipses for the three treatments, coloured the same as the treatments. Any help would be very appreciated!
Thanks.
Upvotes: 2
Views: 13896
Reputation: 59345
There is no stat_ellipse(...)
in the ggplot
package, so you must have retreived it from somewhere else. Care to share?? There are at least two versions that I am aware of, here, and here. Neither of these seem to work with your dataset, which is odd because both have worked with other datasets.
I finally fell back on the option of generating the ellipses externally to ggplot
, which is not that difficult really.
library(ggplot2)
library(ellipse)
centroids <- aggregate(cbind(MDS1,MDS2)~Treatment,sc,mean)
conf.rgn <- do.call(rbind,lapply(unique(sc$Treatment),function(t)
data.frame(Treatment=as.character(t),
ellipse(cov(sc[sc$Treatment==t,1:2]),
centre=as.matrix(centroids[t,2:3]),
level=0.95),
stringsAsFactors=FALSE)))
ggplot(data=sc,(aes(x=MDS1,y=MDS2,colour = Treatment)))+
geom_point(size=3)+
geom_path(data=conf.rgn)+
ggtitle(paste("PCoA of samples at 'class' level(method='Bray')\n",sep=''))+
theme_bw()+
guides(colour = guide_legend(override.aes = list(size=3)))
Upvotes: 8