Just a quick post while I'm waiting for some code to finish running...
If you use stable isotope analysis and R, like I do, you've probably come across the excellent package 'siar'. This is really great for extracting isotope-related stats, like layman metrics or bayesian ellipses. My only niggle is that I'm a ggplot2 advocate and I find that R base graphics tricky and sometimes inconsistent (I also find it much easier to customise the plot with ggplot2). It's not a big deal but I decided I'd prefer to be able to plot isotopic data, along with the standard ellipses that get exported from siar, in ggplot2. Here is the resultant code
The solution is currently a bit clunky. Andrew Jackson (one of the people behind siar) will be looking at implementing a more sophisticated solution into siar, accessible through cran. In the meantime, if you want the code, here it is...
#first you need to subset your data manually by each species (or whatever group you're using)
sp1 <- subset(data, data$species=="sp1")
sp2 <- subset(data, data$species=="sp2")
sp3 <- subset(data, data$species=="sp3")
#next calculate standard.ellipse for each subset
SE1 <- standard.ellipse(sp1$d13C, sp1$d15N)
SE2 <- standard.ellipse(sp2$d13C, sp2$d15N)
SE3 <- standard.ellipse(sp3$d13C, sp3$d15N)
#create name variable for each group (so ggplot2 knows how to colour it)
#the name needs to be the same as your original data frame that contains your isotopic data
sp1_ <- rep("sp1", length(SE1$xSEAc))
sp2_ <- rep("sp2", length(SE2$xSEAc))
sp3_ <- rep("sp3", length(SE3$xSEAc))
#create new data frame with your names and ellipse outputs
species <- c(sp1_,sp2_,sp3_,sp4_)
x <- c(SE1$xSEAc,SE2$xSEAc,SE3$xSEAc)
y <- c(SE1$ySEAc,SE2$ySEAc,SE3$ySEAc
df_SE <- data.frame(x,y,species)
plot(df_SE$x, df_SE$y)
#When you plot SE1$xSEAc vs SE1$ySEAc you should see the ellipse(s) form as a series of points
#This gives you your SEAc (sample size corrected area)
#If you also want hull total area (TA), implement the next bit too
#I got the find_hull function from the stack overflow page linked below
find_hull <- function(df) df[chull(df$d13C, df$d15N), ]
hulls <- ddply(data, "species", find_hull)
#Plot your data. Example below, customise as necessary
plot <- ggplot(data, aes(x=d13C, y=d15N, color=species)) +
theme_bw() +
geom_polygon(data = hulls, alpha = 0.1, linetype=3) +
geom_point(size=3) +
ylab(expression(paste(delta^{15}, "N (\u2030)"))) +
xlab(expression(paste(delta^{13}, "C (\u2030)"))) +
geom_path(data=df_SE, aes(x=x, y=y, color=species), linetype=1, size=2)
plot
R basic plot (left) and new ggplot2 plot (right)
I prefer the ggplot2 version for a few reasons:
- Most/ all of my other plots are ggplot2, so figures in my papers can be more consistent
- I think you can do more with factors, colour scales, legends etc. (maybe I'm just not that good with R base)
