Return to site

#Rstats - Make nice siar plots in ggplot2

#stable isotopes #Rstats #ggplot2

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:

  1. Most/ all of my other plots are ggplot2, so figures in my papers can be more consistent
  2. I think you can do more with factors, colour scales, legends etc. (maybe I'm just not that good with R base)
All Posts
×

Almost done…

We just sent you an email. Please click the link in the email to confirm your subscription!

OKSubscriptions powered by Strikingly