#=====================
#= Bayesian ellypses =
#=====================

#Adapted from https://gist.github.com/AndrewLJackson/b16fb216d8a9a96c20a4a979ec97e7b0
# See ggplot2-isotopes.Rmd

library(SIBER)

library(ggplot2)
library(dplyr)

  
######## Plot ellipses with area only #####

# upload data

mydata <- read.csv("Narwhals_ref_towtusked.csv", header=T)

mydata4 <- mutate(mydata, group = factor(group), community = factor(community))

# plot the data
first.plot <- ggplot(data = mydata4, aes(iso1, iso2)) +
  geom_point(aes(color = group, shape = community), size = 2,
             stroke = ifelse(mydata4$community == "circle" | mydata4$community == "anomalous", 2, 1)) +    
   ylab(expression(paste(delta^{15}, "N (\u2030)")))+
  xlab(expression(paste(delta^{13}, "C (\u2030)"))) + 
  theme(text = element_text(size=15))+
  scale_colour_manual(values=c("black","purple","darkolivegreen4","darkcyan","mediumpurple"),drop=FALSE)+
  scale_fill_manual(values=c("black","purple","darkolivegreen4","darkcyan","mediumpurple"),drop=FALSE)+
  scale_shape_manual(values=c(21, 17))

# add the ellipses and plot
p.ell <- 0.40 

ellipse.plot <- first.plot +  theme_classic() +
  stat_ellipse(aes(group = interaction(community, group), 
                   fill = group, 
                   color = group), 
               alpha = 0.2, 
               level = p.ell,
               type = "norm",
               geom = "polygon")

print(ellipse.plot)

# change plotting options
ellipse.plot2 <-ellipse.plot +  theme_classic() +
  theme(axis.text.x = element_text(size=16)) +
  theme(axis.text.y = element_text(size=16)) +
  theme(axis.title.x = element_text(size=17)) +
  theme(axis.title.y = element_text(size=17))


# 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)
# here I name the mean of iso1 as mC for meanCarbon, and similarly for 
# sdC, and mN, sdN in the produced table.
sbg <- mydata4 %>% group_by(group) %>% 
  summarise(count = length(group),
            mC = mean(iso1), 
            sdC = sd(iso1), 
            mN = mean(iso2), 
            sdN = sd(iso2) )

# add the means and error bars to the plots
second.plot <- ellipse.plot2 + 
  geom_point(data = sbg, aes(mC, mN, group = group), 
             color = c("black","purple","darkolivegreen4","darkcyan","mediumpurple"), shape = 22, size = 5,
            alpha = 1) + # to play with the transparence
  geom_errorbar(data = sbg, 
                mapping = aes(x = mC,
                          #    ymin = mN - 1.96*sdN, # I change it to use SD only
                          ymin = mN - sdN,
                              ymax = mN + sdN), 
                color = c("black","purple","darkolivegreen4","darkcyan","mediumpurple"),
                width = 0, inherit.aes = FALSE)+
  geom_errorbarh(data = sbg, 
                 mapping = aes(y = mN,
                               xmin = mC - sdC,
                               xmax = mC + sdC),
                 color = c("black","purple","darkolivegreen4","darkcyan","mediumpurple"),
                 height = 0, inherit.aes = FALSE)
print(second.plot)

# Add the anonalous narwhals for which there is too little data to plot any ellipses
new <- read.csv("Narwhals_other_anomalous.csv", header=T)

third.plot <- second.plot + 
  geom_point(data = new, aes(color = group, shape = community), size = 2, 
             stroke = ifelse(new$community == "anomalous", 2, 1)) +
  scale_colour_manual(values=c("black","purple","darkolivegreen4","darkcyan","mediumpurple"), drop=FALSE) +
  scale_shape_manual(values=c(21, 17, 8))

print(third.plot)

# We then edit the plot on powerpoint to change symbols and legend