KAPOW! example

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")
## Warning in data("mongooseData"): data set 'mongooseData' not found
# 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)
pack n
11 128
14A 3
14B 30
17 36
19 42
1B 118
1H 112
2 140
21 42
24 12
4B 29
4E 8
7A 83

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)
## Warning: `group_by_()` was deprecated in dplyr 0.7.0.
## ℹ Please use `group_by()` instead.
## ℹ See vignette('programming') for more help
## ℹ The deprecated feature was likely used in the SIBER package.
##   Please report the issue to the authors.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
# 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)
## Warning: The `guide` argument in `scale_*()` cannot be `FALSE`. This was deprecated in
## ggplot2 3.3.4.
## ℹ Please use "none" instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.