In NorskRegnesentral/pupR: Harp and hooded seal pup production estimation
knitr::opts_chunk$set(
collapse = TRUE,
comment = "#>",
fig.width = 7,
fig.height = 5
)
library(pupR)
Overview
The purpose of this vignette is to show how to use the pupR package for estimating the pup abundance of harp and hooded seal populations based on aerial photos. It will guide you through how to prepare and read the data, and how to estimate the pup abundance. Also there is a function for finding a correction factor to correct the estimated abundance for pups which have not been born yet or pups that have left ice based on multiple staging data. A complete analysis of a survey caried out on harp and hooded seals in 2012 in the West Ice (along the coast of Greenland) is included to demonstrate the use of this package. In the appendix each method for estimating the abundance and the method for estimating the correction factor is described.
To load the pupR package type
library(pupR)
Data used and how to load them
To estimate the pup abundance of harp and hooded seals we take aerial photos over the whelping patch. The number of pups are counted on each photo either manually or through machine learning techniques. The whelping patch is normally a large area so it is hard to cover the complete patch with aerial photos. Instead strip transects with a certain distance between tracs are flown. There is also no complete coverage of each transect by aerial photos. So we use different estimation techniques (described in this vignette) to estimate the total number of pups born from the sampled aera.
The whelping season lasts for several weeks so at the time of conducting the aerial survey you might risk that some seals has already matured enough to have left the ice and some seals might not have been born yet. This means that the estimate you have obtained from the aerial photos is underestimated. However due to knowledghe of the length of various age stages (some stages lasts for hours, and others for days) one can estimate a correction factor to scale up the estimated obtained from aerial photos to include missing seals.
So, to get an estimate of the total abundance of seal pups you need two data sets. One from the aerial photos, and another where you over a period of time inspect the whelping patch and count the number of seals in the various age stage groups.
Demo data
When installing the pupR package a demo data set is included so that one can explore how to use the package.
To load the demo data:
data("pupsDemoData")
The demo data is from a survey carried out in the West Ice in 2012. When loading the demo data you get a list containing three elements. One element is a list containing the counts from the aerial photos, dataList. The other two elements, called HarpStages and HoodedStages contains the staging data for the harp and hooded seals, respectievely. The dataList list contains the pup counts for each photo for both the harp and the hooded population, along with lots of other information such as photo ID, which transect a given photo is etc. See table 1 for a complete overview of the variables in the data list.
names(demoData)
The data files
All data that is described below will be loaded and used to estimate the pup abundance of the harp and hooded seals from aerial strip transect data. Having the data on a strip transect format is only required by one of the estimation methods. The others are more general and can deal with different sampling formats.
The data set is divided in two parts. One part contains a file where the following information is needed in separate columns:
- Reader ID
- Photo ID
- Which transect the photo is on (not required)
- Latitude of photo
- Longitude of photo
- Altitude of the photo. Along with some known camera parameters this is used to find the size of the area covered by each photo.
- Area covered by the photo (this is only added if the size of the area covered is known but not required. Using the altitude where the photo is taken from along with the parameters this is calculated)
- Number of pups on a given photo. One separate column for each reader. Some photos are read twice and a "true" number of pups is agreed upon by the readers. This is used to correct for readers errors.
- Comments, this column you can add various comments. If some photos for some reason are to be removed then write Remove in this column.
Important: The names of the various colums have been hard-coded in order to avoid trouble when loading the data. This might be changed in a future release of the package.
The other part of the data set is the staging data. The staging data is a series of staging surveys for a period of time where we find the proportion of pups in various pre-specified stages (with known duration). The staging data can either be obtained from flying low and count the number of seals in each stage, or having land based counts. These staging data is used to estimate a correction factor used to correct for pups that has not been born yet or pups that has left the ice. For harp seal pups this factor is genenrally close to 1 as the duration of the various stages is quite long so that the likelihood of catching all the pups when carrying out the photo survey is quite high. The duration various stages of the hooded pups are much shorter so it is likely to miss some pups at the time of the photo survey. Due to this the correction factor often turns out to be significant.
To load the data set we use the loadAndPrepareData function. To read the data set the you need to specify the name of the survey (survey), the name of the data file (fname), wether it is harp or hooded seals or both (population), and various camera parameters used to calculate the length and the width (and thus the area coverered) of the photos (CamParamLength, CamParamWidth, and CamParam) based on the altitude the photos are taken from. These camera parameters are unique to the camera used.
This function not only loads the data but does some preparations of the data to make it ready for the estimation routines.
dataList = loadAndPrepareData(survey = "WestIce2012",
fname = "WestIce2012.csv",
population = c("harp","hood"))
The function writes out information of the proportion of seals
Next we correct the pup counts for readers error using the correctReading function, i.e.,
dataList = demoData$dataList
harpStages = demoData$HarpStages
hoodedStages = demoData$HoodedStages
correctedList = correctReading(dataList = dataList)
Estimating the number of pups
In this section we will show how to estimate the number of pups based on the pup counts from the aerial photos.
The Kingsley method
We first present a method developed by Kingsley (ref). This method requires that the data is on a strip transect format. We use the kingsley function to estimate the pup abundance and we often estimate the pup abundance based on both the uncorrected pup counts and the corrected pup counts.
# The harp seal population
#Estimate the number of pups using the uncorrected pup counts
EstHarp = kingsley(dataList$xycord$x,
dataList$data$HarpFinalCounts,
dataList$data$Area,
dataList$data$Transect)
#Estimate the number of pups using readers corrected pup counts
EstHarpCorr = kingsley(correctedList$xycord$x,
correctedList$data$CountsHarpCorr,
correctedList$data$Area,
correctedList$data$Transect)
# The hooded seal population
#Estimate the number of pups using the uncorrected pup counts
EstHooded = kingsley(dataList$xycord$x,
dataList$data$HoodedFinalCounts,
dataList$data$Area,
dataList$data$Transect)
#Estimate the number of pups using readers corrected pup counts
EstHoodedCorr = kingsley(correctedList$xycord$x,
correctedList$data$CountsHoodedCorr,
correctedList$data$Area,
correctedList$data$Transect)
The output of the kingsley function is a data frame containing the estimated number of pups N, the standard error SE, and the coefficient of variation CV.
When correcting for readers errors one estimate an uncertainty of the readers correction. This is obtained from either VmeasHarp or VmeasHooded in the output list from the correctReading function. Thus, both the Kingsley method and correcting for readers errors has an uncertainty, and to find the total uncertainty of the estimated pup abundance one need to combine these uncertainties, i.e.,
#Total uncertainty estimate of pup estimation + uncertainty of reader errors
VarHarpCorr = EstHarpCorr$SE^2+correctedList$VmeasHarp
VarHoodedCorr = EstHoodedCorr$SE^2+correctedList$VmeasHooded
GAM based pup abundance
In this section we present a method based on Generalized Additive Models (GAMs) to estimate the pup density over the whelping patch. From this we obtain both an estimate and a uncertainty of the pup abundance. For details of the GAM method see the Appendix.
To estimate the pup abundance using GAMs we use the GAMestimate function. As we did for the Kingsley method we estimate the pup abundance using both the uncorrected and the corrected pup counts. We can also specify for wich population to estimate the pup abundance, or both. In the example below we estimate the pup abundance for both populations.
In addition we have to specify for which distribution to use in the model. At the moment one can choose between the Poisson- or the Negative binomial- distributions, depending on the sparsity of the data. In this example we use the negative binomial distribution.
To estimate the pup abundance using uncorrected pup counts:
EstGAM = GAMestimate(harpcounts =
dataList$data$HarpFinalCounts,
hoodedcounts =
dataList$data$HoodedFinalCounts,
area = dataList$data$Area,
xycord = dataList$xycord,
transect = dataList$data$Transect,
distr = "negbin")
To estimate the pup abundance using corrected pup counts:
EstGAMCorr = GAMestimate(harpcounts =
correctedList$data$CountsHarpCorr,
hoodedcounts =correctedList$data$CountsHoodedCorr,
area = correctedList$data$Area,
xycord = correctedList$xycord,
transect = correctedList$data$Transect,
distr = "negbin")
We could do the same with a Poisson distribution, and then check the goodnes-off-fit for the different models. If you do that it will show that the Negative Binomial distribution is a better model due to over-dispersion.
The total uncertainty estimate of the pup abundance using GAMs using corrected pup counts is:
VarHarpGAMCorr = EstGAMCorr$SEHarpGAM^2+correctedList$VmeasHarp
VarHoodedGAMCorr = EstGAMCorr$SEHoodedGAM^2+correctedList$VmeasHooded
Estimating the birth distribution - the correction factor
In this section we will demonstrate how you can estimate the correction factor to correct for pups which have not been born yet or pups that have left the ice. The correction factor is estimated from a different data set than the aerial photos. The data needed is staging data. These data is obtained by flying over the whelping patch for a period of time and count the number of pups in various pre-defined age/growt stages with known length, see Appendix for model details. This applies for both populations, but each population are divided in different stage classes and the length of the various stages are different. After the pups reach a certain mature stage it is assumed that the pup leaves the ice, and is thus unavailable for counting if the aerial photo survey are carried out after the pups have left the ice. Based on these data and the known length of the various stages one can estimate the date of the peak of the whelping season and through a model developed by (ref) one can estimate the number of pups on the ice at any date in the whelping season.
To estimate the birth distribution along with the correction factor for the date of the photographic survey we use the birthDist function, i.e.,
#Estimating the birth distribution for the harp seals
bdistHarp = birthDist(data = harpStages,
harpLengthStages = c(2.4,4.42,11.39),
harpKappa = 12.4,
population = "harp",
datePhoto = 28)
#Estimating the birth distribution for the hooded seals
bdistHooded = birthDist(data = hoodedStages,
hoodedLengthStages = c(2,1,4),
hoodedKappa = 8.6,
population = "hooded",
datePhoto = 28)
Note that we have to estimate the birth distribution separately for each population.
As we now have estimated the birth distribution and obtained the correction factor for the date when the aerial survey was carried out we can obtain a pup abundance estimated adusted for pups that have not been born yet and pups that have left the ice. Earlier we found a total uncertainty of the estimated pup abundance (not adjusted for missing seals) when using the readers corrected pup counts by combining the uncertainty in the estimation procedure and the uncertainty of the readers correction. We need to do the same for the final pup abundance estimate which adjusts for missing seals based on the readers corrected data.
We have combined all the estimates in a list:
QHarp = bdistHarp$PropIce # Correction factor for harp seals
SEQHarp = bdistHarp$PropIceSD # Uncertainty of the harp seal correction factor
QHooded = bdistHooded$PropIce # Correction factor for harp seals
SEQHooded = bdistHooded$PropIceSD # Uncertainty of the harp seal correction factor
# Correcting for missing seals
estimates = list()
# Harp seal population
# Estimated pup abundance based on raw pup counts
estimates$EstHarpFinal = EstHarp$N/QHarp
# Estimated pup abundance based on readers corrected pup counts
estimates$EstHarpFinalCorr = EstHarpCorr$N/QHarp
#---------------------------------------------------
# Hooded seal population
# Estimated pup abundance based on raw pup counts
estimates$EstHoodedFinal = EstHooded$N/QHooded
# Estimated pup abundance based on readers corrected pup counts
estimates$EstHoodedFinalCorr = EstHoodedCorr$N/QHooded
# Calculating the total uncertainty
estimates$VarTotalHarp = (1/QHarp)^2*EstHarp$SE^2+(EstHarp$N/QHarp^2)*SEQHarp^2
estimates$VarTotalHooded = (1/QHooded)^2*EstHooded$SE^2+(EstHooded$N/QHooded^2)*SEQHooded^2
estimates$VarTotalHarpCorr = (1/QHarp)^2*VarHarpCorr+(EstHarpCorr$N/QHarp^2)*SEQHarp^2
estimates$VarTotalHoodedCorr = (1/QHooded)^2*VarHoodedCorr+(EstHoodedCorr$N/QHooded^2)*SEQHooded^2
We do the same for the GAM estimates:
estimatesGAM = list()
estimatesGAM$EstHarpGAMFinal = EstGAM$NHarpGAM/QHarp
estimatesGAM$EstHarpGAMFinalCorr = EstGAMCorr$NHarpGAM/QHarp
estimatesGAM$EstHoodedGAMFinal = EstGAM$NHoodedGAM/QHooded
estimatesGAM$EstHoodedGAMFinalCorr = EstGAMCorr$NHoodedGAM/QHooded
estimatesGAM$VarTotalHarpGAM = (1/QHarp)^2*EstGAM$VarHarpGAM+(EstGAM$NHarpGAM/QHarp^2)*SEQHarp^2
estimatesGAM$VarTotalHarpGAMCorr = (1/QHarp)^2*VarHarpGAMCorr+(EstGAMCorr$NHarpGAM/QHarp^2)*SEQHarp^2
estimatesGAM$VarTotalHoodedGAM = (1/QHooded)^2*EstGAM$VarHoodedGAM+(EstGAM$NHoodedGAM/QHooded^2)*SEQHooded^2
estimatesGAM$VarTotalHoodedGAMCorr = (1/QHooded)^2*VarHoodedGAMCorr+(EstGAMCorr$NHoodedGAM/QHooded^2)*SEQHooded^2
In summary we now have both uncorrected and readers corrected pup counts (for both populations), we have pup abundance estimates (both populations, and based on both corrected and uncorrected pup counts), we got a correction factor for pups not born yet and pups that have already left the ice, and we got a total uncertainty in the final pup abundance estimate (based on the uncertainty of the readers corrections and the uncertainty of the correction factor from the birth distribution).
With the list we can print out the final results to the screen.
For the Kingley method:
#Print pup production estimates based on the Kingsley method to screen
printKingsley2screen(estimates)
For the GAM method:
#Print pup production estimates based on the GAM method to screen
printGAM2screen(estimatesGAM)
Appendix
Complete script for estimating the pup production
In this section we present a complete analysis of the 2012 harp and hooded survey carried out in the West Ice (along the east coast of Greenland).
# Loading the pupR package
library(pupR)
########################
# Loading the data
########################
# If using the demo data set - loading the demo data set
data("pupsDemoData")
# Break up the elements of the list in the demo data
dataList = demoData$dataList
harpStages = demoData$HarpStages
hoodedStages = demoData$HoodedStages
# If using the complete data set - loading the data
dataList = loadAndPrepareData(survey = "WestIce2012",
fname = "WestIce2012.csv",
population = c("harp","hood"))
# Correct for readers errors
correctedList = correctReading(dataList = dataList)
###################################
# Abundance estimation
###################################
# The Kingsley method
#*********************
#
# The harp seal population
#----------------------------
#
# Using the uncorrected pup counts
EstHarp = kingsley(dataList$xycord$x,
dataList$data$HarpFinalCounts,
dataList$data$Area,
dataList$data$Transect)
# Using readers corrected pup counts
EstHarpCorr = kingsley(correctedList$xycord$x,
correctedList$data$CountsHarpCorr,
correctedList$data$Area,
correctedList$data$Transect)
# The hooded seal population
#------------------------------
#
# Using the uncorrected pup counts
EstHooded = kingsley(dataList$xycord$x,
dataList$data$HoodedFinalCounts,
dataList$data$Area,
dataList$data$Transect)
# Using readers corrected pup counts
EstHoodedCorr = kingsley(correctedList$xycord$x,
correctedList$data$CountsHoodedCorr,
correctedList$data$Area,
correctedList$data$Transect)
# Calculating the total uncertainty estimate of the abundance estimate
VarHarpCorr = EstHarpCorr$SE^2+correctedList$VmeasHarp
VarHoodedCorr = EstHoodedCorr$SE^2+correctedList$VmeasHooded
# The GAM method
#*****************
#
# Both the harp and the hooded populationS
# Using the pup counts not corrected for readers errors
EstGAM = GAMestimate(harpcounts =
dataList$data$HarpFinalCounts,
hoodedcounts =
dataList$data$HoodedFinalCounts,
area = dataList$data$Area,
xycord = dataList$xycord,
transect = dataList$data$Transect,
distr = "negbin")
# Using readers corrected pup counts
EstGAMCorr = GAMestimate(harpcounts =
correctedList$data$CountsHarpCorr,
hoodedcounts =correctedList$data$CountsHoodedCorr,
area = correctedList$data$Area,
xycord = correctedList$xycord,
transect = correctedList$data$Transect,
distr = "negbin")
# Calculating the total uncertainty
VarHarpGAMCorr = EstGAMCorr$SEHarpGAM^2+correctedList$VmeasHarp
VarHoodedGAMCorr = EstGAMCorr$SEHoodedGAM^2+correctedList$VmeasHooded
# Correct for seals not born yet or seals which has left the ice
#*****************************************************************
#Estimating the birth distribution for the harp seals
bdistHarp = birthDist(data = harpStages,
harpLengthStages = c(2.4,4.42,11.39),
harpKappa = 12.4,
population = "harp",
datePhoto = 28)
#Estimating the birth distribution for the hooded seals
bdistHooded = birthDist(data = hoodedStages,
hoodedLengthStages = c(2,1,4),
hoodedKappa = 8.6,
population = "hooded",
datePhoto = 28)
# Extracting the correction factor and the uncertainty
QHarp = bdistHarp$PropIce # Correction factor for harp seals
SEQHarp = bdistHarp$PropIceSD # Uncertainty of the harp seal correction factor
QHooded = bdistHooded$PropIce # Correction factor for harp seals
SEQHooded = bdistHooded$PropIceSD # Uncertainty of the harp seal correction factor
# Correcting for missing seals
estimates = list()
# Correcting the Kingsley method estimates
#------------------------------------------
# Harp seal population
# Estimated pup abundance based on raw pup counts
estimates$EstHarpFinal = EstHarp$N/QHarp
# Estimated pup abundance based on readers corrected pup counts
estimates$EstHarpFinalCorr = EstHarpCorr$N/QHarp
# Hooded seal population
# Estimated pup abundance based on raw pup counts
estimates$EstHoodedFinal = EstHooded$N/QHooded
# Estimated pup abundance based on readers corrected pup counts
estimates$EstHoodedFinalCorr = EstHoodedCorr$N/QHooded
# Calculating the total uncertainty
estimates$VarTotalHarp = (1/QHarp)^2*EstHarp$SE^2+(EstHarp$N/QHarp^2)*SEQHarp^2
estimates$VarTotalHooded = (1/QHooded)^2*EstHooded$SE^2+(EstHooded$N/QHooded^2)*SEQHooded^2
estimates$VarTotalHarpCorr = (1/QHarp)^2*VarHarpCorr+(EstHarpCorr$N/QHarp^2)*SEQHarp^2
estimates$VarTotalHoodedCorr = (1/QHooded)^2*VarHoodedCorr+(EstHoodedCorr$N/QHooded^2)*SEQHooded^2
# Correcting the GAM estimates
#-------------------------------
estimatesGAM = list()
estimatesGAM$EstHarpGAMFinal = EstGAM$NHarpGAM/QHarp
estimatesGAM$EstHarpGAMFinalCorr = EstGAMCorr$NHarpGAM/QHarp
estimatesGAM$EstHoodedGAMFinal = EstGAM$NHoodedGAM/QHooded
estimatesGAM$EstHoodedGAMFinalCorr = EstGAMCorr$NHoodedGAM/QHooded
estimatesGAM$VarTotalHarpGAM = (1/QHarp)^2*EstGAM$VarHarpGAM+(EstGAM$NHarpGAM/QHarp^2)*SEQHarp^2
estimatesGAM$VarTotalHarpGAMCorr = (1/QHarp)^2*VarHarpGAMCorr+(EstGAMCorr$NHarpGAM/QHarp^2)*SEQHarp^2
estimatesGAM$VarTotalHoodedGAM = (1/QHooded)^2*EstGAM$VarHoodedGAM+(EstGAM$NHoodedGAM/QHooded^2)*SEQHooded^2
estimatesGAM$VarTotalHoodedGAMCorr = (1/QHooded)^2*VarHoodedGAMCorr+(EstGAMCorr$NHoodedGAM/QHooded^2)*SEQHooded^2
# Print the results to screen
#******************************
#
# Pup production estimates based on the Kingsley method:
printKingsley2screen(estimates)
# Pup production estimates based on the GAM method
printGAM2screen(estimatesGAM)
Pup production estimation methods
The photographic surveys are based on a systematic sampling design with a single random start and a sampling unit of transects of variable length. The sampled data are pup counts from aerial photos. In this section we will review the Kingsley method and the GAM method.
The Kingsley method
The estimated number of pups on the ice at the time of survey may be written as [@Salberg2008; @Øigård2010]:
$$\hat{N}=T\sum_{j=1}^{J}W_jx_j$$
where $W_j=l_j/A_j$, $A_j$ is the area covered of all photographs on transect $j$, $l_j$ is the length of transect $j$, $J$ is the total number of transects, and $x_j=\sum_{k=1}^{P_{ij}}\hat{n}_{j,k}$ is the sum of the corrected counts on transect $j$. The number of photos on the $j$th transect is $P_j$ and $T$ is the spacing between transects in the survey. This estimator takes into account changes in transect width along transects and between transects due to changes in flight altitude.
The estimates of error variance $V^s$, based on serial differences between transects were calculated as [@Salberg2008]:
$$V^s=\frac{TJ}{2(J-1)}\Biggl(T-\frac{\sum_{j=1}^{J}A_j}{\sum_{j=1}^{J}l_j}\Biggr)\sum_{j=1}^{J}l_j(W_jx_j-W_{j+1}x_{j+1})^2$$
This estimator assumes that the mean is constant between two neighboring transects. For the seal pup data this assumption is often not valid due to clustered data, and we will have an unwanted contribution from the difference between the transect count mean values which causes an overestimate of the variance of the pup production estimate [@Cochran1977]. However, if the seals are homogenously spread over a large area this assumption is fine.
The variance associated with mis-classification of pups, i.e., readers errors, for the whole survey is then [@Salberg2008]:
$$\begin{aligned} V^{meas}=T^2\Biggl[\sum_{j=1}^{J}W^2_jP_j\sigma^2+\Biggl(\sum_{j=1}^{J}W_jP_j\Biggr)^2var(\alpha)+2cov(\alpha,\beta)\Biggl(\sum_{j=1}^{J}W_jP_j\Biggr) \ \Biggl(\sum_{j=1}^{J}Wj\sum_{k=1}^{P_j}n_{j,k}\Biggr)+var(\beta)\Biggl(\sum_{j=1}^{J}Wj\sum_{k=1}^{P_j}n_{j,k}\Biggr)^2\Biggr] \end{aligned}$$
If the intercept term is not statistically significant on a specified level it could be dropped from the linear correction model. The variance expression is then simplified to
$$V^{meas}=T^2\Biggl[\sum_{j=1}^{J}W^2_jP_j\sigma^2+\Biggl(\sum_{j=1}^{J}W_jP_j\Biggr)^2\Biggr]$$
To obtain the total sampling variance of the survey, the variance associated with the mis-identification corrections $V^{meas}$ was added to the sampling variance $V^s$, i.e.:
$$V=V^s+V^{meas}$$
The GAM method
In the survey analysis developed in Salberg et al. (2009) the data was analyzed using spatial modeling methods based on Generalized Additive Models. Even thought we are dealing with count data, a Poisson error distribution would not be appropriate because the seal count data are over-dispersed. Overdispersed data are often a consequence of the population being clustered (McCullagh and Nelder, 1989). We therefore assume that the data are negative binomial distributed, which has been applied previous to model over-dispersed count data (Thurston et al., 2000; Hinde and Demétrio, 1998; Gardner et al., 1995; Augustin et al., 1998).
The counted number of pups in the $k$th photograph is
$$n_k = A_kd_k$$
where $A_k$ is the area covered by the $k$th photograph, and $d_k$ is the density of pups in photograph $k$. We assume that $n_k$ is negative binomial distributed with mean $\mu_k$, and shape parameter $\kappa$. The purpose of the GAM is to model the pup density over the patch as a function of spatial location ${\bf x} = [x_1, x_2]$. Using the logarithmic link function $g(\mu_k) = \eta_k = \log(\mu_k)$, the pup density is modeled as
$$\mu_k=e^{\eta_k}=\exp\left[\log\left(A_k\right)+S\left(x_k\right)\right],$$
where $\log(A_k)$ is an offset variable, and $S(\cdot)$ is a smoothing function of the spatial covariates, and ${\bf x}k$ is the spatial location at the $k$th sampling point. Note that the estimated expected seal density is then
$$d_k=\exp\left(S\left(x_k\right)\right].$$
The smoothing function $S({\bf x})$ is modeled using a thin plate regression spline, and the degree of smoothness is determined using generalized cross validation (Wood, 2003). Once the model has been chosen, we may predict the seal density at any location in the patch. Hence, the GAM provides us with a smoothed expected seal density surface over the entire survey area. To estimate the total pup production in the patch, we numerically integrate the expected seal density over space (Augustin et al., 1998; Borchers et al., 1997; Salberg et al., 2009)
$$\hat{N}{GAM}=\sum_{\text{fine grid } j}\hat{N}{GAM}^j$$
where $\hat{N}{GAM}^j$ is the estimated number of pups in the $j$th fine grid are at spatial location $x_j$. The method used to calculate the variance of the pup production follows exactly the procedures described in detail in Salberg et al. (2009), and has the following compact form
$$V_i^{GAM}=\mu^T_gX_g\text{cov}(\beta)X_g^T\mu_g,$$
where ${\bf\mu}_g$ is a vector collecting all fine grid pup density estimates, ${\bf X}_g$ is the matrix that linearly maps the estimated parameter vector $\beta$ to the smoothed expected seal density surface. We refer to the original manuscript for a detailed explanation of the various quantities (Salberg et al., 2009).
Estimating the birth distribution
This was done by using information on the proportion of pups in seven distinct age-dependent stages. These easily recognizable descriptive age categories were based on pelage colour and body condition, overall appearance, and muscular coordination, as described for the northwest Atlantic harp seals by [@Stewart1980]:
- Newborn: Pup still wet, bright yellow colour often present. Often associated with wet placentas and blood stained snow.
- Yellowcoat: Pup dry, yellow amniotic stain still persistent on pelt. The pup is lean and moving awkwardly.
- Thin whitecoat: Amniotic stain faded, pup with visible neck and often conical in shape, pelage white.
- Fat whitecoat: Visibly fatter, neck not visible, cylindrical in shape, pelage still white.
- Greycoat: Darker juvenile pelt beginning to grow in under the white lanugo giving a grey cast to the pelt, “salt-and-pepper”-look in later stages.
- Ragged-jackets: Lanugo shed in patches, at least a handful from torso (nose, tail and flippers do not count).
- Beaters: Fully moulted pups (a handful of lanugo may remain).
Prior to the survey, classifications of pup stages were standardized among observers to ensure consistency. To determine the proportion of pups in each stage on a given day, random samples of pups were obtained by flying a series of transects over the patch. Pups were classified from the helicopter hovering just above the animals. The spacing between transects depended on the size of the actual patch.
A similar procedure was followed for hooded seals where information on the proportion of pups in each of five distinct age-dependent stages was used to assess the temporal distribution of births.
These arbitrary, but easily recognizable age categories were based on pelage colour and body condition, overall appearance, and muscular coordination, as described for northwest Atlantic hooded seals by @Bowen1987 and @Stenson1988, and used in the previous surveys in the Greenland Sea [@Salberg2008; @Øigård2010]:
- Unborn: Parturient females.
- Newborn: Skin in loose folds along flanks, fur saturated to wet, entire pelage with yellowish hue, awkward body movements. Mother present. Often associated with wet placentas and blood stained snow.
- Thin blueback: Pup dry, ventrum white, neck well defined, trunk conical in shape. Mother present. Mainly 1-2 days old.
- Fat blueback: Ventrum white, neck not distinguishable, trunk fusiform in shape. Mother present. Mainly 2-4 days old.
- Solitary blueback: As in fat blueback, but mother not present. Mainly 4 days or older.
The temporal distribution of births for both harp and hooded seals was estimated using the method developed in @Reed1968 and adapted for modelling the birth distribution for harp and hooded seals in @Bowen1987, and @Myers1989. The life cycles of the seals were assumed to be divided into $k$ identifiable age-dependent stages $S_1,.....,S_k$. Birth takes place into state $S_1$ and the pup then progresses in succession through states $S_1,S_2,....$ until it attains maturity when reaching state $S_k$. All pups reaching state $S_k$ eventually die in that state, either from hunting or natural causes [@Reed1968]. We assumed that for both seal populations the birth rate could be adequately described by a continuous function of time, $m_1(t)$ which denoted the temporal distribution of births. The distribution of births over time was assumed to be a normal distribution with mean value $\mu_1$ and standard deviation $\sigma_1$.
The various development stages are denoted by the subscript $j$, and a pup passes from stage $j$ to stage $j+1$. The stage durations are specified in terms of transition intensity functions $\phi_j(t)$, which is the probability that an animal passes from stage $j$ to $j+1$ in the interval $[\tau,\tau+\Delta t]$ and has survived. Here $\tau$ is the time spent in stage $j$. The stage duration was assumed to be a semi-Markov process, i.e. the transition intensities depend only on the current stage and the time so far spent in that stage [@Bowen1987]. The rate at which pups enter the stage $j$ at time $t$ were denoted by $m_j(t)$ and given by a recurrence relationship @Myers1989:
$$m_j(t)={\int_0^\infty}m_{j-1}(t-\tau)\phi_{j-1}(\tau)d\tau\quad j=1,..,k$$
The proportion of pups that will be observed on the ice in stage $j$ at time $t$ is [@Bowen1987; @Myers1989]:
$$n_j(t)={\int_0^\infty}m_{j-1}(t-\tau)\Biggl(1-{\int_0^\tau}\phi(s)ds\Biggr)d\tau$$
This equation assumes no pup morality during these stages and that all pups on the ice are visible. In @Bowen1987, -@Bowen2007 and @Myers1989 the transition intensity functions $\phi_j(t)$ were assumed to follow a Gamma distribution with shape parameter $\kappa_j$ and scale parameter $\rho_j$ for stage $j$. The product between the shape parameter and the scale parameter, $\rho_j\kappa_j$, gives the mean duration of stage $j$.
The numbers of individuals observed to be of stage $j$ at time $t_i$ were denoted $S_{ij}$. The $S_{ij}$'s were obtained by taking a random sample of the pup abundance and determining the stage of each individual. The predicted proportions of each stage present at time $t_i$, $P_{ij}$, are calculated as in @Myers1989, i.e. by estimating the parameters $\hat{\mu}_1$ and $\hat{\sigma}_1$ of the birth distribution. The proportion of pups on the ice at time $t$ was estimated using [@Salberg2008; @Øigård2010]:
$$Q(t)={\sum_{j=1}^k}\eta_j(t)$$
The estimated variance of the proportion of pups on the ice at a given time was estimated by simulating from the proportion of pups in the various stages obtained from the staging by simulating from a multinomial distribution with $k$ stages [@Salberg2008].
To correct for pups still not born, and pups that had left the ice at the time of the photographic survey, the estimated numbers of pups on the ice at the time of the survey were corrected by:
$$\hat{N}^{corr}=\frac{\hat{N}}{\hat{Q}}$$
where $\hat{Q}$ is the estimated proportion of pups visible on the photographs at the time of the survey.
The estimates of $N_i$ and $Q$ are independent and therefore the error variance of the estimated total number of pups born in the patch $\hat{N}^{corr}$ may be obtained using the $\delta$-method [e.g. @Casella1990]:
$$V^{corr}=\Biggl(\frac{1}{Q}\Biggr)^2V\Biggl(\frac{N}{Q^2}\Biggr)V^q$$
where $V^q$ is the estimated variance of $\hat{Q}$.
To conform to the procedure used in 2012, we used the following binning of the various stages of the harp seal pups: stage 1 = Newborn/Yellow, stage 2 = Thin white, and stage 3 = Fat white/Greycoat.
The following binning of the various stages of the hooded seal pups was: stage 1 = Newborn and Thin, stage 2 = Fat, and stage 3 = Solitary.
References
NorskRegnesentral/pupR documentation built on Sept. 17, 2020, 5:52 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
knitr::opts_chunk$set( collapse = TRUE, comment = "#>", fig.width = 7, fig.height = 5 )
library(pupR)
Overview
The purpose of this vignette is to show how to use the pupR package for estimating the pup abundance of harp and hooded seal populations based on aerial photos. It will guide you through how to prepare and read the data, and how to estimate the pup abundance. Also there is a function for finding a correction factor to correct the estimated abundance for pups which have not been born yet or pups that have left ice based on multiple staging data. A complete analysis of a survey caried out on harp and hooded seals in 2012 in the West Ice (along the coast of Greenland) is included to demonstrate the use of this package. In the appendix each method for estimating the abundance and the method for estimating the correction factor is described.
To load the pupR package type
library(pupR)
Data used and how to load them
To estimate the pup abundance of harp and hooded seals we take aerial photos over the whelping patch. The number of pups are counted on each photo either manually or through machine learning techniques. The whelping patch is normally a large area so it is hard to cover the complete patch with aerial photos. Instead strip transects with a certain distance between tracs are flown. There is also no complete coverage of each transect by aerial photos. So we use different estimation techniques (described in this vignette) to estimate the total number of pups born from the sampled aera.
The whelping season lasts for several weeks so at the time of conducting the aerial survey you might risk that some seals has already matured enough to have left the ice and some seals might not have been born yet. This means that the estimate you have obtained from the aerial photos is underestimated. However due to knowledghe of the length of various age stages (some stages lasts for hours, and others for days) one can estimate a correction factor to scale up the estimated obtained from aerial photos to include missing seals.
So, to get an estimate of the total abundance of seal pups you need two data sets. One from the aerial photos, and another where you over a period of time inspect the whelping patch and count the number of seals in the various age stage groups.
Demo data
When installing the pupR package a demo data set is included so that one can explore how to use the package.
To load the demo data:
data("pupsDemoData")
The demo data is from a survey carried out in the West Ice in 2012. When loading the demo data you get a list containing three elements. One element is a list containing the counts from the aerial photos, dataList. The other two elements, called HarpStages and HoodedStages contains the staging data for the harp and hooded seals, respectievely. The dataList list contains the pup counts for each photo for both the harp and the hooded population, along with lots of other information such as photo ID, which transect a given photo is etc. See table 1 for a complete overview of the variables in the data list.
names(demoData)
The data files
All data that is described below will be loaded and used to estimate the pup abundance of the harp and hooded seals from aerial strip transect data. Having the data on a strip transect format is only required by one of the estimation methods. The others are more general and can deal with different sampling formats.
The data set is divided in two parts. One part contains a file where the following information is needed in separate columns:
- Reader ID
- Photo ID
- Which transect the photo is on (not required)
- Latitude of photo
- Longitude of photo
- Altitude of the photo. Along with some known camera parameters this is used to find the size of the area covered by each photo.
- Area covered by the photo (this is only added if the size of the area covered is known but not required. Using the altitude where the photo is taken from along with the parameters this is calculated)
- Number of pups on a given photo. One separate column for each reader. Some photos are read twice and a "true" number of pups is agreed upon by the readers. This is used to correct for readers errors.
- Comments, this column you can add various comments. If some photos for some reason are to be removed then write Remove in this column.
Important: The names of the various colums have been hard-coded in order to avoid trouble when loading the data. This might be changed in a future release of the package.
The other part of the data set is the staging data. The staging data is a series of staging surveys for a period of time where we find the proportion of pups in various pre-specified stages (with known duration). The staging data can either be obtained from flying low and count the number of seals in each stage, or having land based counts. These staging data is used to estimate a correction factor used to correct for pups that has not been born yet or pups that has left the ice. For harp seal pups this factor is genenrally close to 1 as the duration of the various stages is quite long so that the likelihood of catching all the pups when carrying out the photo survey is quite high. The duration various stages of the hooded pups are much shorter so it is likely to miss some pups at the time of the photo survey. Due to this the correction factor often turns out to be significant.
To load the data set we use the loadAndPrepareData function. To read the data set the you need to specify the name of the survey (survey), the name of the data file (fname), wether it is harp or hooded seals or both (population), and various camera parameters used to calculate the length and the width (and thus the area coverered) of the photos (CamParamLength, CamParamWidth, and CamParam) based on the altitude the photos are taken from. These camera parameters are unique to the camera used.
This function not only loads the data but does some preparations of the data to make it ready for the estimation routines.
dataList = loadAndPrepareData(survey = "WestIce2012", fname = "WestIce2012.csv", population = c("harp","hood"))
The function writes out information of the proportion of seals
Next we correct the pup counts for readers error using the correctReading function, i.e.,
dataList = demoData$dataList harpStages = demoData$HarpStages hoodedStages = demoData$HoodedStages
correctedList = correctReading(dataList = dataList)
Estimating the number of pups
In this section we will show how to estimate the number of pups based on the pup counts from the aerial photos.
The Kingsley method
We first present a method developed by Kingsley (ref). This method requires that the data is on a strip transect format. We use the kingsley function to estimate the pup abundance and we often estimate the pup abundance based on both the uncorrected pup counts and the corrected pup counts.
# The harp seal population #Estimate the number of pups using the uncorrected pup counts EstHarp = kingsley(dataList$xycord$x, dataList$data$HarpFinalCounts, dataList$data$Area, dataList$data$Transect) #Estimate the number of pups using readers corrected pup counts EstHarpCorr = kingsley(correctedList$xycord$x, correctedList$data$CountsHarpCorr, correctedList$data$Area, correctedList$data$Transect) # The hooded seal population #Estimate the number of pups using the uncorrected pup counts EstHooded = kingsley(dataList$xycord$x, dataList$data$HoodedFinalCounts, dataList$data$Area, dataList$data$Transect) #Estimate the number of pups using readers corrected pup counts EstHoodedCorr = kingsley(correctedList$xycord$x, correctedList$data$CountsHoodedCorr, correctedList$data$Area, correctedList$data$Transect)
The output of the kingsley function is a data frame containing the estimated number of pups N, the standard error SE, and the coefficient of variation CV.
When correcting for readers errors one estimate an uncertainty of the readers correction. This is obtained from either VmeasHarp or VmeasHooded in the output list from the correctReading function. Thus, both the Kingsley method and correcting for readers errors has an uncertainty, and to find the total uncertainty of the estimated pup abundance one need to combine these uncertainties, i.e.,
#Total uncertainty estimate of pup estimation + uncertainty of reader errors VarHarpCorr = EstHarpCorr$SE^2+correctedList$VmeasHarp VarHoodedCorr = EstHoodedCorr$SE^2+correctedList$VmeasHooded
GAM based pup abundance
In this section we present a method based on Generalized Additive Models (GAMs) to estimate the pup density over the whelping patch. From this we obtain both an estimate and a uncertainty of the pup abundance. For details of the GAM method see the Appendix.
To estimate the pup abundance using GAMs we use the GAMestimate function. As we did for the Kingsley method we estimate the pup abundance using both the uncorrected and the corrected pup counts. We can also specify for wich population to estimate the pup abundance, or both. In the example below we estimate the pup abundance for both populations.
In addition we have to specify for which distribution to use in the model. At the moment one can choose between the Poisson- or the Negative binomial- distributions, depending on the sparsity of the data. In this example we use the negative binomial distribution.
To estimate the pup abundance using uncorrected pup counts:
EstGAM = GAMestimate(harpcounts = dataList$data$HarpFinalCounts, hoodedcounts = dataList$data$HoodedFinalCounts, area = dataList$data$Area, xycord = dataList$xycord, transect = dataList$data$Transect, distr = "negbin")
To estimate the pup abundance using corrected pup counts:
EstGAMCorr = GAMestimate(harpcounts = correctedList$data$CountsHarpCorr, hoodedcounts =correctedList$data$CountsHoodedCorr, area = correctedList$data$Area, xycord = correctedList$xycord, transect = correctedList$data$Transect, distr = "negbin")
We could do the same with a Poisson distribution, and then check the goodnes-off-fit for the different models. If you do that it will show that the Negative Binomial distribution is a better model due to over-dispersion.
The total uncertainty estimate of the pup abundance using GAMs using corrected pup counts is:
VarHarpGAMCorr = EstGAMCorr$SEHarpGAM^2+correctedList$VmeasHarp VarHoodedGAMCorr = EstGAMCorr$SEHoodedGAM^2+correctedList$VmeasHooded
Estimating the birth distribution - the correction factor
In this section we will demonstrate how you can estimate the correction factor to correct for pups which have not been born yet or pups that have left the ice. The correction factor is estimated from a different data set than the aerial photos. The data needed is staging data. These data is obtained by flying over the whelping patch for a period of time and count the number of pups in various pre-defined age/growt stages with known length, see Appendix for model details. This applies for both populations, but each population are divided in different stage classes and the length of the various stages are different. After the pups reach a certain mature stage it is assumed that the pup leaves the ice, and is thus unavailable for counting if the aerial photo survey are carried out after the pups have left the ice. Based on these data and the known length of the various stages one can estimate the date of the peak of the whelping season and through a model developed by (ref) one can estimate the number of pups on the ice at any date in the whelping season.
To estimate the birth distribution along with the correction factor for the date of the photographic survey we use the birthDist function, i.e.,
#Estimating the birth distribution for the harp seals bdistHarp = birthDist(data = harpStages, harpLengthStages = c(2.4,4.42,11.39), harpKappa = 12.4, population = "harp", datePhoto = 28) #Estimating the birth distribution for the hooded seals bdistHooded = birthDist(data = hoodedStages, hoodedLengthStages = c(2,1,4), hoodedKappa = 8.6, population = "hooded", datePhoto = 28)
Note that we have to estimate the birth distribution separately for each population.
As we now have estimated the birth distribution and obtained the correction factor for the date when the aerial survey was carried out we can obtain a pup abundance estimated adusted for pups that have not been born yet and pups that have left the ice. Earlier we found a total uncertainty of the estimated pup abundance (not adjusted for missing seals) when using the readers corrected pup counts by combining the uncertainty in the estimation procedure and the uncertainty of the readers correction. We need to do the same for the final pup abundance estimate which adjusts for missing seals based on the readers corrected data.
We have combined all the estimates in a list:
QHarp = bdistHarp$PropIce # Correction factor for harp seals SEQHarp = bdistHarp$PropIceSD # Uncertainty of the harp seal correction factor QHooded = bdistHooded$PropIce # Correction factor for harp seals SEQHooded = bdistHooded$PropIceSD # Uncertainty of the harp seal correction factor # Correcting for missing seals estimates = list() # Harp seal population # Estimated pup abundance based on raw pup counts estimates$EstHarpFinal = EstHarp$N/QHarp # Estimated pup abundance based on readers corrected pup counts estimates$EstHarpFinalCorr = EstHarpCorr$N/QHarp #--------------------------------------------------- # Hooded seal population # Estimated pup abundance based on raw pup counts estimates$EstHoodedFinal = EstHooded$N/QHooded # Estimated pup abundance based on readers corrected pup counts estimates$EstHoodedFinalCorr = EstHoodedCorr$N/QHooded # Calculating the total uncertainty estimates$VarTotalHarp = (1/QHarp)^2*EstHarp$SE^2+(EstHarp$N/QHarp^2)*SEQHarp^2 estimates$VarTotalHooded = (1/QHooded)^2*EstHooded$SE^2+(EstHooded$N/QHooded^2)*SEQHooded^2 estimates$VarTotalHarpCorr = (1/QHarp)^2*VarHarpCorr+(EstHarpCorr$N/QHarp^2)*SEQHarp^2 estimates$VarTotalHoodedCorr = (1/QHooded)^2*VarHoodedCorr+(EstHoodedCorr$N/QHooded^2)*SEQHooded^2
We do the same for the GAM estimates:
estimatesGAM = list() estimatesGAM$EstHarpGAMFinal = EstGAM$NHarpGAM/QHarp estimatesGAM$EstHarpGAMFinalCorr = EstGAMCorr$NHarpGAM/QHarp estimatesGAM$EstHoodedGAMFinal = EstGAM$NHoodedGAM/QHooded estimatesGAM$EstHoodedGAMFinalCorr = EstGAMCorr$NHoodedGAM/QHooded estimatesGAM$VarTotalHarpGAM = (1/QHarp)^2*EstGAM$VarHarpGAM+(EstGAM$NHarpGAM/QHarp^2)*SEQHarp^2 estimatesGAM$VarTotalHarpGAMCorr = (1/QHarp)^2*VarHarpGAMCorr+(EstGAMCorr$NHarpGAM/QHarp^2)*SEQHarp^2 estimatesGAM$VarTotalHoodedGAM = (1/QHooded)^2*EstGAM$VarHoodedGAM+(EstGAM$NHoodedGAM/QHooded^2)*SEQHooded^2 estimatesGAM$VarTotalHoodedGAMCorr = (1/QHooded)^2*VarHoodedGAMCorr+(EstGAMCorr$NHoodedGAM/QHooded^2)*SEQHooded^2
In summary we now have both uncorrected and readers corrected pup counts (for both populations), we have pup abundance estimates (both populations, and based on both corrected and uncorrected pup counts), we got a correction factor for pups not born yet and pups that have already left the ice, and we got a total uncertainty in the final pup abundance estimate (based on the uncertainty of the readers corrections and the uncertainty of the correction factor from the birth distribution).
With the list we can print out the final results to the screen.
For the Kingley method:
#Print pup production estimates based on the Kingsley method to screen printKingsley2screen(estimates)
For the GAM method:
#Print pup production estimates based on the GAM method to screen printGAM2screen(estimatesGAM)
Appendix
Complete script for estimating the pup production
In this section we present a complete analysis of the 2012 harp and hooded survey carried out in the West Ice (along the east coast of Greenland).
# Loading the pupR package library(pupR) ######################## # Loading the data ######################## # If using the demo data set - loading the demo data set data("pupsDemoData") # Break up the elements of the list in the demo data dataList = demoData$dataList harpStages = demoData$HarpStages hoodedStages = demoData$HoodedStages # If using the complete data set - loading the data dataList = loadAndPrepareData(survey = "WestIce2012", fname = "WestIce2012.csv", population = c("harp","hood")) # Correct for readers errors correctedList = correctReading(dataList = dataList) ################################### # Abundance estimation ################################### # The Kingsley method #********************* # # The harp seal population #---------------------------- # # Using the uncorrected pup counts EstHarp = kingsley(dataList$xycord$x, dataList$data$HarpFinalCounts, dataList$data$Area, dataList$data$Transect) # Using readers corrected pup counts EstHarpCorr = kingsley(correctedList$xycord$x, correctedList$data$CountsHarpCorr, correctedList$data$Area, correctedList$data$Transect) # The hooded seal population #------------------------------ # # Using the uncorrected pup counts EstHooded = kingsley(dataList$xycord$x, dataList$data$HoodedFinalCounts, dataList$data$Area, dataList$data$Transect) # Using readers corrected pup counts EstHoodedCorr = kingsley(correctedList$xycord$x, correctedList$data$CountsHoodedCorr, correctedList$data$Area, correctedList$data$Transect) # Calculating the total uncertainty estimate of the abundance estimate VarHarpCorr = EstHarpCorr$SE^2+correctedList$VmeasHarp VarHoodedCorr = EstHoodedCorr$SE^2+correctedList$VmeasHooded # The GAM method #***************** # # Both the harp and the hooded populationS # Using the pup counts not corrected for readers errors EstGAM = GAMestimate(harpcounts = dataList$data$HarpFinalCounts, hoodedcounts = dataList$data$HoodedFinalCounts, area = dataList$data$Area, xycord = dataList$xycord, transect = dataList$data$Transect, distr = "negbin") # Using readers corrected pup counts EstGAMCorr = GAMestimate(harpcounts = correctedList$data$CountsHarpCorr, hoodedcounts =correctedList$data$CountsHoodedCorr, area = correctedList$data$Area, xycord = correctedList$xycord, transect = correctedList$data$Transect, distr = "negbin") # Calculating the total uncertainty VarHarpGAMCorr = EstGAMCorr$SEHarpGAM^2+correctedList$VmeasHarp VarHoodedGAMCorr = EstGAMCorr$SEHoodedGAM^2+correctedList$VmeasHooded # Correct for seals not born yet or seals which has left the ice #***************************************************************** #Estimating the birth distribution for the harp seals bdistHarp = birthDist(data = harpStages, harpLengthStages = c(2.4,4.42,11.39), harpKappa = 12.4, population = "harp", datePhoto = 28) #Estimating the birth distribution for the hooded seals bdistHooded = birthDist(data = hoodedStages, hoodedLengthStages = c(2,1,4), hoodedKappa = 8.6, population = "hooded", datePhoto = 28) # Extracting the correction factor and the uncertainty QHarp = bdistHarp$PropIce # Correction factor for harp seals SEQHarp = bdistHarp$PropIceSD # Uncertainty of the harp seal correction factor QHooded = bdistHooded$PropIce # Correction factor for harp seals SEQHooded = bdistHooded$PropIceSD # Uncertainty of the harp seal correction factor # Correcting for missing seals estimates = list() # Correcting the Kingsley method estimates #------------------------------------------ # Harp seal population # Estimated pup abundance based on raw pup counts estimates$EstHarpFinal = EstHarp$N/QHarp # Estimated pup abundance based on readers corrected pup counts estimates$EstHarpFinalCorr = EstHarpCorr$N/QHarp # Hooded seal population # Estimated pup abundance based on raw pup counts estimates$EstHoodedFinal = EstHooded$N/QHooded # Estimated pup abundance based on readers corrected pup counts estimates$EstHoodedFinalCorr = EstHoodedCorr$N/QHooded # Calculating the total uncertainty estimates$VarTotalHarp = (1/QHarp)^2*EstHarp$SE^2+(EstHarp$N/QHarp^2)*SEQHarp^2 estimates$VarTotalHooded = (1/QHooded)^2*EstHooded$SE^2+(EstHooded$N/QHooded^2)*SEQHooded^2 estimates$VarTotalHarpCorr = (1/QHarp)^2*VarHarpCorr+(EstHarpCorr$N/QHarp^2)*SEQHarp^2 estimates$VarTotalHoodedCorr = (1/QHooded)^2*VarHoodedCorr+(EstHoodedCorr$N/QHooded^2)*SEQHooded^2 # Correcting the GAM estimates #------------------------------- estimatesGAM = list() estimatesGAM$EstHarpGAMFinal = EstGAM$NHarpGAM/QHarp estimatesGAM$EstHarpGAMFinalCorr = EstGAMCorr$NHarpGAM/QHarp estimatesGAM$EstHoodedGAMFinal = EstGAM$NHoodedGAM/QHooded estimatesGAM$EstHoodedGAMFinalCorr = EstGAMCorr$NHoodedGAM/QHooded estimatesGAM$VarTotalHarpGAM = (1/QHarp)^2*EstGAM$VarHarpGAM+(EstGAM$NHarpGAM/QHarp^2)*SEQHarp^2 estimatesGAM$VarTotalHarpGAMCorr = (1/QHarp)^2*VarHarpGAMCorr+(EstGAMCorr$NHarpGAM/QHarp^2)*SEQHarp^2 estimatesGAM$VarTotalHoodedGAM = (1/QHooded)^2*EstGAM$VarHoodedGAM+(EstGAM$NHoodedGAM/QHooded^2)*SEQHooded^2 estimatesGAM$VarTotalHoodedGAMCorr = (1/QHooded)^2*VarHoodedGAMCorr+(EstGAMCorr$NHoodedGAM/QHooded^2)*SEQHooded^2 # Print the results to screen #****************************** # # Pup production estimates based on the Kingsley method: printKingsley2screen(estimates) # Pup production estimates based on the GAM method printGAM2screen(estimatesGAM)
Pup production estimation methods
The photographic surveys are based on a systematic sampling design with a single random start and a sampling unit of transects of variable length. The sampled data are pup counts from aerial photos. In this section we will review the Kingsley method and the GAM method.
The Kingsley method
The estimated number of pups on the ice at the time of survey may be written as [@Salberg2008; @Øigård2010]:
$$\hat{N}=T\sum_{j=1}^{J}W_jx_j$$ where $W_j=l_j/A_j$, $A_j$ is the area covered of all photographs on transect $j$, $l_j$ is the length of transect $j$, $J$ is the total number of transects, and $x_j=\sum_{k=1}^{P_{ij}}\hat{n}_{j,k}$ is the sum of the corrected counts on transect $j$. The number of photos on the $j$th transect is $P_j$ and $T$ is the spacing between transects in the survey. This estimator takes into account changes in transect width along transects and between transects due to changes in flight altitude. The estimates of error variance $V^s$, based on serial differences between transects were calculated as [@Salberg2008]:
$$V^s=\frac{TJ}{2(J-1)}\Biggl(T-\frac{\sum_{j=1}^{J}A_j}{\sum_{j=1}^{J}l_j}\Biggr)\sum_{j=1}^{J}l_j(W_jx_j-W_{j+1}x_{j+1})^2$$ This estimator assumes that the mean is constant between two neighboring transects. For the seal pup data this assumption is often not valid due to clustered data, and we will have an unwanted contribution from the difference between the transect count mean values which causes an overestimate of the variance of the pup production estimate [@Cochran1977]. However, if the seals are homogenously spread over a large area this assumption is fine.
The variance associated with mis-classification of pups, i.e., readers errors, for the whole survey is then [@Salberg2008]:
$$\begin{aligned} V^{meas}=T^2\Biggl[\sum_{j=1}^{J}W^2_jP_j\sigma^2+\Biggl(\sum_{j=1}^{J}W_jP_j\Biggr)^2var(\alpha)+2cov(\alpha,\beta)\Biggl(\sum_{j=1}^{J}W_jP_j\Biggr) \ \Biggl(\sum_{j=1}^{J}Wj\sum_{k=1}^{P_j}n_{j,k}\Biggr)+var(\beta)\Biggl(\sum_{j=1}^{J}Wj\sum_{k=1}^{P_j}n_{j,k}\Biggr)^2\Biggr] \end{aligned}$$
If the intercept term is not statistically significant on a specified level it could be dropped from the linear correction model. The variance expression is then simplified to
$$V^{meas}=T^2\Biggl[\sum_{j=1}^{J}W^2_jP_j\sigma^2+\Biggl(\sum_{j=1}^{J}W_jP_j\Biggr)^2\Biggr]$$
To obtain the total sampling variance of the survey, the variance associated with the mis-identification corrections $V^{meas}$ was added to the sampling variance $V^s$, i.e.:
$$V=V^s+V^{meas}$$
The GAM method
In the survey analysis developed in Salberg et al. (2009) the data was analyzed using spatial modeling methods based on Generalized Additive Models. Even thought we are dealing with count data, a Poisson error distribution would not be appropriate because the seal count data are over-dispersed. Overdispersed data are often a consequence of the population being clustered (McCullagh and Nelder, 1989). We therefore assume that the data are negative binomial distributed, which has been applied previous to model over-dispersed count data (Thurston et al., 2000; Hinde and Demétrio, 1998; Gardner et al., 1995; Augustin et al., 1998).
The counted number of pups in the $k$th photograph is $$n_k = A_kd_k$$ where $A_k$ is the area covered by the $k$th photograph, and $d_k$ is the density of pups in photograph $k$. We assume that $n_k$ is negative binomial distributed with mean $\mu_k$, and shape parameter $\kappa$. The purpose of the GAM is to model the pup density over the patch as a function of spatial location ${\bf x} = [x_1, x_2]$. Using the logarithmic link function $g(\mu_k) = \eta_k = \log(\mu_k)$, the pup density is modeled as $$\mu_k=e^{\eta_k}=\exp\left[\log\left(A_k\right)+S\left(x_k\right)\right],$$ where $\log(A_k)$ is an offset variable, and $S(\cdot)$ is a smoothing function of the spatial covariates, and ${\bf x}k$ is the spatial location at the $k$th sampling point. Note that the estimated expected seal density is then $$d_k=\exp\left(S\left(x_k\right)\right].$$ The smoothing function $S({\bf x})$ is modeled using a thin plate regression spline, and the degree of smoothness is determined using generalized cross validation (Wood, 2003). Once the model has been chosen, we may predict the seal density at any location in the patch. Hence, the GAM provides us with a smoothed expected seal density surface over the entire survey area. To estimate the total pup production in the patch, we numerically integrate the expected seal density over space (Augustin et al., 1998; Borchers et al., 1997; Salberg et al., 2009) $$\hat{N}{GAM}=\sum_{\text{fine grid } j}\hat{N}{GAM}^j$$ where $\hat{N}{GAM}^j$ is the estimated number of pups in the $j$th fine grid are at spatial location $x_j$. The method used to calculate the variance of the pup production follows exactly the procedures described in detail in Salberg et al. (2009), and has the following compact form $$V_i^{GAM}=\mu^T_gX_g\text{cov}(\beta)X_g^T\mu_g,$$ where ${\bf\mu}_g$ is a vector collecting all fine grid pup density estimates, ${\bf X}_g$ is the matrix that linearly maps the estimated parameter vector $\beta$ to the smoothed expected seal density surface. We refer to the original manuscript for a detailed explanation of the various quantities (Salberg et al., 2009).
Estimating the birth distribution
This was done by using information on the proportion of pups in seven distinct age-dependent stages. These easily recognizable descriptive age categories were based on pelage colour and body condition, overall appearance, and muscular coordination, as described for the northwest Atlantic harp seals by [@Stewart1980]:
- Newborn: Pup still wet, bright yellow colour often present. Often associated with wet placentas and blood stained snow.
- Yellowcoat: Pup dry, yellow amniotic stain still persistent on pelt. The pup is lean and moving awkwardly.
- Thin whitecoat: Amniotic stain faded, pup with visible neck and often conical in shape, pelage white.
- Fat whitecoat: Visibly fatter, neck not visible, cylindrical in shape, pelage still white.
- Greycoat: Darker juvenile pelt beginning to grow in under the white lanugo giving a grey cast to the pelt, “salt-and-pepper”-look in later stages.
- Ragged-jackets: Lanugo shed in patches, at least a handful from torso (nose, tail and flippers do not count).
- Beaters: Fully moulted pups (a handful of lanugo may remain).
Prior to the survey, classifications of pup stages were standardized among observers to ensure consistency. To determine the proportion of pups in each stage on a given day, random samples of pups were obtained by flying a series of transects over the patch. Pups were classified from the helicopter hovering just above the animals. The spacing between transects depended on the size of the actual patch.
A similar procedure was followed for hooded seals where information on the proportion of pups in each of five distinct age-dependent stages was used to assess the temporal distribution of births. These arbitrary, but easily recognizable age categories were based on pelage colour and body condition, overall appearance, and muscular coordination, as described for northwest Atlantic hooded seals by @Bowen1987 and @Stenson1988, and used in the previous surveys in the Greenland Sea [@Salberg2008; @Øigård2010]:
- Unborn: Parturient females.
- Newborn: Skin in loose folds along flanks, fur saturated to wet, entire pelage with yellowish hue, awkward body movements. Mother present. Often associated with wet placentas and blood stained snow.
- Thin blueback: Pup dry, ventrum white, neck well defined, trunk conical in shape. Mother present. Mainly 1-2 days old.
- Fat blueback: Ventrum white, neck not distinguishable, trunk fusiform in shape. Mother present. Mainly 2-4 days old.
- Solitary blueback: As in fat blueback, but mother not present. Mainly 4 days or older.
The temporal distribution of births for both harp and hooded seals was estimated using the method developed in @Reed1968 and adapted for modelling the birth distribution for harp and hooded seals in @Bowen1987, and @Myers1989. The life cycles of the seals were assumed to be divided into $k$ identifiable age-dependent stages $S_1,.....,S_k$. Birth takes place into state $S_1$ and the pup then progresses in succession through states $S_1,S_2,....$ until it attains maturity when reaching state $S_k$. All pups reaching state $S_k$ eventually die in that state, either from hunting or natural causes [@Reed1968]. We assumed that for both seal populations the birth rate could be adequately described by a continuous function of time, $m_1(t)$ which denoted the temporal distribution of births. The distribution of births over time was assumed to be a normal distribution with mean value $\mu_1$ and standard deviation $\sigma_1$.
The various development stages are denoted by the subscript $j$, and a pup passes from stage $j$ to stage $j+1$. The stage durations are specified in terms of transition intensity functions $\phi_j(t)$, which is the probability that an animal passes from stage $j$ to $j+1$ in the interval $[\tau,\tau+\Delta t]$ and has survived. Here $\tau$ is the time spent in stage $j$. The stage duration was assumed to be a semi-Markov process, i.e. the transition intensities depend only on the current stage and the time so far spent in that stage [@Bowen1987]. The rate at which pups enter the stage $j$ at time $t$ were denoted by $m_j(t)$ and given by a recurrence relationship @Myers1989:
$$m_j(t)={\int_0^\infty}m_{j-1}(t-\tau)\phi_{j-1}(\tau)d\tau\quad j=1,..,k$$
The proportion of pups that will be observed on the ice in stage $j$ at time $t$ is [@Bowen1987; @Myers1989]:
$$n_j(t)={\int_0^\infty}m_{j-1}(t-\tau)\Biggl(1-{\int_0^\tau}\phi(s)ds\Biggr)d\tau$$
This equation assumes no pup morality during these stages and that all pups on the ice are visible. In @Bowen1987, -@Bowen2007 and @Myers1989 the transition intensity functions $\phi_j(t)$ were assumed to follow a Gamma distribution with shape parameter $\kappa_j$ and scale parameter $\rho_j$ for stage $j$. The product between the shape parameter and the scale parameter, $\rho_j\kappa_j$, gives the mean duration of stage $j$. The numbers of individuals observed to be of stage $j$ at time $t_i$ were denoted $S_{ij}$. The $S_{ij}$'s were obtained by taking a random sample of the pup abundance and determining the stage of each individual. The predicted proportions of each stage present at time $t_i$, $P_{ij}$, are calculated as in @Myers1989, i.e. by estimating the parameters $\hat{\mu}_1$ and $\hat{\sigma}_1$ of the birth distribution. The proportion of pups on the ice at time $t$ was estimated using [@Salberg2008; @Øigård2010]:
$$Q(t)={\sum_{j=1}^k}\eta_j(t)$$
The estimated variance of the proportion of pups on the ice at a given time was estimated by simulating from the proportion of pups in the various stages obtained from the staging by simulating from a multinomial distribution with $k$ stages [@Salberg2008].
To correct for pups still not born, and pups that had left the ice at the time of the photographic survey, the estimated numbers of pups on the ice at the time of the survey were corrected by:
$$\hat{N}^{corr}=\frac{\hat{N}}{\hat{Q}}$$
where $\hat{Q}$ is the estimated proportion of pups visible on the photographs at the time of the survey.
The estimates of $N_i$ and $Q$ are independent and therefore the error variance of the estimated total number of pups born in the patch $\hat{N}^{corr}$ may be obtained using the $\delta$-method [e.g. @Casella1990]:
$$V^{corr}=\Biggl(\frac{1}{Q}\Biggr)^2V\Biggl(\frac{N}{Q^2}\Biggr)V^q$$
where $V^q$ is the estimated variance of $\hat{Q}$.
To conform to the procedure used in 2012, we used the following binning of the various stages of the harp seal pups: stage 1 = Newborn/Yellow, stage 2 = Thin white, and stage 3 = Fat white/Greycoat.
The following binning of the various stages of the hooded seal pups was: stage 1 = Newborn and Thin, stage 2 = Fat, and stage 3 = Solitary.
References
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.