SIA plots using ggplot2

Import the data as before. Note that we get some warnings from both that “objects are masked from” various packages. This is because the packages we have just loaded have functions of the same name as those that have already been loaded (usually ones from the base R packages). This warning is telling us that, in the case of filter, when we simply call filter(), we will be using the last loaded one, i.e. from dplyr, rather than the one from stats. If you want to force R to use a particular function from a particular package you can write the long form, dplyr::filter() or stats::filter(). In many ways, it is good practice to always use this format, but in reality we are all lazy.

#install.packages("tidyverse")

library(tidyverse)

# import the data. 
# If you get an error about "No such file or directory" then you 
# need to take care about where R is currently working, and where your 
# data file is in relation to this.
mydata <- read.csv("Practical01.csv", 
                   header = TRUE, 
                   stringsAsFactors = FALSE)

# As per the summary statistics file, I am going to remove our outlier
# note the use of dplyr::filter() to avoid loading the entire 
# dplyr pacakge for just one function
mydata <- dplyr::filter(mydata, d13C < max(d13C))

# verify that our data looks correct by printing the first1 10 lines
# to screen
head(mydata)

# if the data is not a massive dataset, you might like to look  
# at it all
mydata

The graphics that are part of “base R” are not very pretty (or at least some people think so - I quite like them to be honest). Another one of Hadley Wickham’s very popular packages is ggplot2 which by default makes some very blog-friendly graphics. And of course you can change the theme (template) on them to get something more suitable for publication (and fashion seems to be changing here too).

In ggplot2, figures are created in layer, by first creating a basic layer of axes according to an “aesthetic” which is then used as a framework to add points and lines and other embellishments. If you push the layers into an object using the <- assignment as I have done here, then you will need to print() your figure to get it to render.

first.plot <- ggplot(data = mydata, 
                     mapping = aes(x = d13C, 
                                   y = d15N)) + 
  geom_point(aes(color = Taxon, shape = Taxon), size = 5) +
  ylab(expression(paste(delta^{15}, "N (\u2030)"))) +
  xlab(expression(paste(delta^{13}, "C (\u2030)"))) + 
  theme(text = element_text(size=20))

# And print our plot to screen
print(first.plot)

If you want to use the more normal plotting format for journals without gridlines, and without the light-grey background, you add theme_classic() to your plot. This is a shortcut to some default settings of the myriad options available in theme() which we used above to increase the text size in our plot.

classic.first.plot <- first.plot + theme_classic() + 
  theme(text = element_text(size=35)) + 
  coord_equal() + 
  theme(axis.ticks.length = unit(-0.2, "cm")) +
  scale_colour_viridis_d(end = 0.9)

# and print to screen
print(classic.first.plot)

# options to add to point the axis tick marks inwards
# theme(axis.ticks.length = unit(-0.1, "cm"))

Errorbar biplots

Adding the means and error bars in both directions to create the classic isotope biplot requires adding these as additional layers. We could re-write all the lines of code for the original plot above, or we can simply create a new ggplot object based on the original, and then add the layers we want, and print. We need some summary data for each of the groups to plot the intervals: we use dplyr::summarise as in the previous script today except this time we store the output into a new object which I call sbg. When we add the layers for the means and the errorbars, we need to tell the layers to use a new aesthetic mapping using the new x and y data: mapping = aes(x = mC, y = mN, ...).

# Summarise By Group (sbg)
sbg <- mydata %>% 
  group_by(Taxon) %>% 
  summarise(count = n(),
            mC = mean(d13C), 
            sdC = sd(d13C), 
            mN = mean(d15N), 
            sdN = sd(d15N) )

# make a copy of the first.plot object
# second.plot <- first.plot

# add the layers using the summary data in sbg
second.plot <- first.plot + 
  geom_errorbar(data = sbg, 
                mapping = aes(x = mC, y = mN,
                              ymin = mN - 1.96*sdN, 
                              ymax = mN + 1.96*sdN), 
                width = 0) +
  geom_errorbarh(data = sbg, 
                 mapping = aes(x = mC, y = mN,
                               xmin = mC - 1.96*sdC,
                               xmax = mC + 1.96*sdC),
                 height = 0) + 
  geom_point(data = sbg, aes(x = mC, 
                             y = mN,
                             fill = Taxon), 
             color = "black", shape = 22, size = 5,
             alpha = 0.7, show.legend = FALSE) + 
  coord_equal()

Warning: Ignoring unknown aesthetics: x

print(second.plot)

NA

Ellipses instead of errorbars

The ggplot2 function stat_ellipse allows us to easily add ellipses, of varying level which corresponds to the prediction interval. This function defaults to using the t-distribution so we will override this and specify the normal distribution as is more fitting with the SIBER approach.

# use our ellipse function to generate the ellipses for plotting

# decide how big an ellipse you want to draw
p.ell <- 0.95

# create our plot based on first.plot above
# adding the stat_ellipse() geometry. We 
# specify thee ellipse to be plotted using 
# the polygon geom, with fill and edge colour
# defined by Taxon as a grouping variable, 
# using the normal distribution and with 
# a quite high level of transparency.
ellipse.plot <- first.plot + 
  stat_ellipse(aes(group = Taxon, 
                   fill = Taxon, 
                   color = Taxon), 
               alpha = 0.25, 
               level = p.ell,
               type = "norm",
               geom = "polygon")

print(ellipse.plot)

ggsave(plot = ellipse.plot, filename = "../../images/biplot.jpeg")

Saving 9 x 6 in image

Trouble shooting

If the permil symbol ‰ is not showing correctly (and instead is printing as \(\text{\u2030}\) in your plot it is likely because your computer is not set up to use UTF-8 format character encoding. This is not a problem with your setup of R or Rstudio, but is deeper in your computer. It is fixed by changing the region or locale settings on your computer. You can access these in the system preferences area of your computer’s operating system. To check if your computer is set up to identify and interpret UTF-8 encoding, you can type sessionInfo() in the R console. You should see something like this: en_IE.UTF-8/en_IE.UTF-8/en_IE.UTF-8/C/en_IE.UTF-8/en_IE.UTF-8 under the heading “locale”. On my machine, this indicates that it is use english (en), for Ireland (IE) and the UTF-8 encoding format. If the UTF-8 format is missing from your sessionInfo() then you could try changing your operating system to a locale similar to your own region’s location and see if the UTF-8 format is supported there.