################################################################################################################################
## Supporting information for the article: Meyer et al. (2017) "Marine mammal population decline linked to obscured by-catch" ##
## This file includes R code to calculate the Petersen estimates as published in ##
## Childerhouse S, Muller C, Burns T, French R, Kay E (2016) ##
## "BPM-16-Final Report for CSP Project NZ sea lion ground component 2015-16 v1.3" ##
## (Tech. Rep., Blue Planet Marine). ##
## Standard deviation of pup production in 1997 at Sandy Bay was calculated from the reported ##
## standard error and replicate size; for the year 1995 at Sandy Bay and the years 1995 and 1997 ##
## at Dundas Islands raw data or information on sample size is not available. ##
################################################################################################################################
rm(list=ls())
## -- Load packages ---
library(data.table)
## -- Load data from mark-recapture study ---
setwd(" ... ") # set working directory
dt <- fread("Dataset_S1.csv", header = TRUE)
## -- Calculate Petersen estimate (Pi) and add number of dead pups to get full pup production estimate
## Pi = { [ ((M+1) * (Ci+1))/ (Ri + 1) ] - 1 } + No_of_dead_pups
## M: no. of pups marked in breeding season
## Ci: no. of pups examined in recapture sample i (i.e. marked and unmarked)
## Ri: no. of marked pups in recapture sample i
dt[, Pi := ( ( (No_of_marked_individuals+1) * (No_of_marked_individuals_recaptured + No_of_unmarked_individuals_captured + 1)
/ (No_of_marked_individuals_recaptured + 1) ) - 1) + No_of_dead_pups, by = c("Subpopulation", "Year", "Sample")]
## -- Pi statistics (i.e. mean, sd, se, 95% confidence intervals)
dt[, Pi_mean := round(mean(Pi), 0), by = c("Subpopulation", "Year")]
dt[, Pi_sd := round(sd(Pi), 0), by = c("Subpopulation", "Year")]
dt[, Pi_se := round(Pi_sd / sqrt(max(Sample)), 0), by = c("Subpopulation", "Year")]
dt[, Pi_lower := Pi_mean - (Pi_se * 1.96), by = c("Subpopulation", "Year")]
dt[, Pi_upper := Pi_mean + (Pi_se * 1.96), by = c("Subpopulation", "Year")]
## -- Prepare table for AR(1) model
dt_mean <- dt[ , .(unique(Pi_mean), unique(Pi_sd), unique(Pi_se), max(Sample)), by = c("Subpopulation", "Year")]
names(dt_mean) <- c("Subpopulation", "Year", "Pi_mean", "Pi_sd", "Pi_se", "Sample_size")
setkey(dt_mean, Subpopulation, Year)
## -- Add means and standard errors reported for 1995 and 1997 for sub-populations at Sandy Bay and Dundas Island
dt_mean[Subpopulation=="Sandy_Bay" & Year==1995, Pi_mean := 467] # Table 1 in Gales, N. J., & Fletcher, D. J. (1999). Abundance, distribution and status of the New Zealand sea lion, Phocarctos hookeri. Wildlife Research, 26(1), 35-52.
dt_mean[Subpopulation=="Sandy_Bay" & Year==1995, Pi_se := 4] # Table 1 in Gales, N. J., & Fletcher, D. J. (1999). Abundance, distribution and status of the New Zealand sea lion, Phocarctos hookeri. Wildlife Research, 26(1), 35-52.
dt_mean[Subpopulation=="Sandy_Bay" & Year==1996, Pi_mean := 455] # Table 1 in Gales, N. J., & Fletcher, D. J. (1999). Abundance, distribution and status of the New Zealand sea lion, Phocarctos hookeri. Wildlife Research, 26(1), 35-52.
dt_mean[Subpopulation=="Sandy_Bay" & Year==1996, Pi_se := 3] # Table 1 in Gales, N. J., & Fletcher, D. J. (1999). Abundance, distribution and status of the New Zealand sea lion, Phocarctos hookeri. Wildlife Research, 26(1), 35-52.
dt_mean[Subpopulation=="Sandy_Bay" & Year==1997, Pi_mean := 509] # Appendix 1 in Childerhouse S, Muller C, Burns T, French R, Kay E (2016) "BPM-16-Final Report for CSP Project NZ sea lion ground component 2015-16 v1.3" (Tech. Rep., Blue Planet Marine).
dt_mean[Subpopulation=="Sandy_Bay" & Year==1997, Pi_se := 9] # Field report for pupping season in 1997
dt_mean[Subpopulation=="Dundas_Island" & Year==1995, Pi_mean := 1837] # Table 1 in Gales, N. J., & Fletcher, D. J. (1999). Abundance, distribution and status of the New Zealand sea lion, Phocarctos hookeri. Wildlife Research, 26(1), 35-52.
dt_mean[Subpopulation=="Dundas_Island" & Year==1995, Pi_se := 20] # Table 1 in Gales, N. J., & Fletcher, D. J. (1999). Abundance, distribution and status of the New Zealand sea lion, Phocarctos hookeri. Wildlife Research, 26(1), 35-52.
dt_mean[Subpopulation=="Dundas_Island" & Year==1996, Pi_mean := 2017] # Table 1 in Gales, N. J., & Fletcher, D. J. (1999). Abundance, distribution and status of the New Zealand sea lion, Phocarctos hookeri. Wildlife Research, 26(1), 35-52.
dt_mean[Subpopulation=="Dundas_Island" & Year==1996, Pi_se := 22] # Table 1 in Gales, N. J., & Fletcher, D. J. (1999). Abundance, distribution and status of the New Zealand sea lion, Phocarctos hookeri. Wildlife Research, 26(1), 35-52.
dt_mean[Subpopulation=="Dundas_Island" & Year==1997, Pi_mean := 2260] # Appendix 1 in Childerhouse S, Muller C, Burns T, French R, Kay E (2016) "BPM-16-Final Report for CSP Project NZ sea lion ground component 2015-16 v1.3" (Tech. Rep., Blue Planet Marine).
dt_mean[Subpopulation=="Dundas_Island" & Year==1997, Pi_se := 24] # Field report for pupping season in 1997
## -- Add sample size in year 1997 for Sandy Bay (known from field report)
dt_mean[Subpopulation=="Sandy_Bay" & Year==1997, Sample_size := 9]
## -- Calculate SD at Sandy Bay in 1997
dt_mean[Subpopulation=="Sandy_Bay" & Year==1997, Pi_sd := Pi_se * sqrt(Sample_size)]
## -- only keep information required for AR(1) model (i.e. year, mean estimate, and sd)
dt_mean[, Pi_se := NULL]
dt_mean[, Sample_size := NULL]
setkey(dt_mean, Subpopulation, Year)
write.csv(dt_mean, file = "Dataset_S2.csv")