In AndrewLJackson/SIBER: Stable Isotope Bayesian Ellipses in R

About

This example is associated with, but is not guaranteed an exact replica of the analysis contained in Sheppard et al. 2018 (Intragroup competition predicts individual foraging specialisation in a group-living mammal. Ecology Letters. doi). It is included here a demonstration of the "KAPOW" function in SIBER and not as a reproducible example of the above cited paper.

The data comprise multiple stable isotope values obtained on individual mongooses who are nested within packs. We calculate the proportional ellipse by individual in the context of its pack. An ellipse is fit to each individual within a pack; and the outline of all ellipses is then calculated as the union of all the ellipses. Finally, each individuals ellipse is calculated as a proportion of the total for that pack. Plotting is performed using ggplot and statistics generated using SIBER.

Setup

First load in the packages. This script relies heavily on the tidyverse approach to data manipulation and avoids for loops wherever possible. Note that we do not explicitly load the tidyverse package which is itself a wrapper that loads a cluster of related packages, and instead load the specific packages we require, which in this case is dplyr and purrr.

library(dplyr)
library(purrr)
library(ggplot2)
library(SIBER)

Import Data

Now load in the data:

# This loads a pre-saved object called mongoose that comprises the 
# dataframe for this analysis.
data("mongooseData")


# Ordinarily we might typically use code like this to import our data from a 
# csv file. 
# mongoose <- read.csv("mongooseFullData.csv", header = TRUE, 
#                      stringsAsFactors = FALSE)

There are lots of individuals with fewer than 4 observations which is definitely the lower limit to fit these models. There also appears to be at least one individual that appears in more than one pack, so we need to group both individual and pack and then check that there are sufficient samples.

# min sample size for individual replicates per pack.
min.n <- 4

mongoose_2 <- mongoose %>% group_by(indiv.id, pack) %>% 
  filter(n() >= min.n) %>% ungroup()

# convert pack and indiv.id to factor
mongoose_2 <- mongoose_2 %>% mutate(indiv.id = factor(indiv.id),
                                    pack = factor(pack))

# count observations 
id_pack_counts <- mongoose %>% count(pack)

knitr::kable(id_pack_counts)

p1 <- ggplot(data = mongoose_2, aes(c13, n15, color = indiv.id)) + 
  geom_point()  + 
  viridis::scale_color_viridis(discrete = TRUE, guide = FALSE) + 
  facet_wrap(~pack)

print(p1)

Visualise the packs ellipses

Split the data into a list by pack and plot the pack's collective isotopic niche as a grey shaded area, with the constituent invididual ellipses in colour along with the raw data.

# split by pack
packs <- mongoose_2 %>% split(.$pack)

# use purrr::map to apply siberKapow across each pack.
pack_boundaries <- purrr::map(packs, siberKapow, isoNames = c("c13","n15"), 
                             group = "indiv.id", pEll = 0.95)


# Define afunction to strip out the boundaries of the union of the 
# ellipses and plot them. This function returns the ggplot2 object
# but doesnt actually do the plotting which is handled afterwards.
plotBoundaries <- function(dd, ee){

  # exdtract the boundary points for each KAPOW shape.
  bdry <- data.frame(dd$bdry)

  # the plot object
  p <- ggplot(data = ee, aes(c13, n15, color = indiv.id)) + 
  geom_point()  + 
  viridis::scale_color_viridis(discrete = TRUE, guide = "legend", name = "Individual") + 
    geom_polygon(data = bdry, mapping = aes(x, y, color = NULL), alpha = 0.2) + 
    viridis::scale_fill_viridis(discrete = TRUE, guide = FALSE) + 
    theme_bw() + 
    ggtitle(paste0("Pack: ", as.character(ee$pack[1]) )) + 
    geom_polygon(data = dd$ell.coords, aes(X1, X2, group = indiv.id), 
                 alpha = 0.2, fill = NA)
  return(p)

}


# map this function over packs and return the un-printed ggplot2 objects
bndry_plots <- purrr::map2(pack_boundaries, packs, plotBoundaries)

# print them to screen / file
print(bndry_plots)

Print tables of proportions

The actual data of the proportional ellipses which is printed to the output document and also saved to *.csv file.

# KAPOW areas for each pack
total.area <- map(pack_boundaries, spatstat.geom::area)

# a function to extract ellipse coordinates, calculate areas and return
# as a vector not a list.
extractProportions <- function(x){unlist(map(x$owin.coords, spatstat.geom::area))}

# map our individual ellipse area function over packs
ellipse.areas <- map(pack_boundaries, . %>% extractProportions)

# calculate ellipses as proportions of the KAPOW for that pack by mapping
# over both the individual ellipses and pack totals and dividing.
ellipse_proportions <- map2(ellipse.areas, total.area, `/`)

# print(ellipse_proportions)

# convert to table with a nested map_df call for easier printing. 
# Probably possible to use at_depth() to simplify this, but possibly
# not as i use map_df() here.
df_proportions <- map_df(ellipse_proportions, 
                         . %>% map_df(data.frame, .id = "individual"), 
                         .id = "pack" )

# rename the ugly variable manually
df_proportions <- rename(df_proportions, Proportion = ".x..i..")

# print a nice table
knitr::kable(df_proportions, digits = 2)

# Optional code to save to csv file
# write.csv(df_proportions, file = "mongoose_kapow_niche_proportions.csv",
#           row.names = FALSE)