Fit with half-normal plateau detection function and home range size estimation using nimbleSCR package

Soumen Dey

2022-11-30

In this vignette, we demonstrate how to use the nimbleSCR (Bischof et al. 2020) and NIMBLE (Valpine et al. 2017; NIMBLE Development Team 2021) packages to simulate and fit a Bayesian SCR model with a half-normal plateau detection function and derive home-range size estimates from the estimated parameters (Dey et al. (2022)). We use a method based on numerical approximation to estimate home range size from SCR models with circular detection functions.

# LOAD PACKAGES
library(coda)
library(nimble)
library(nimbleSCR)
library(basicMCMCplots)

1. Binomial SCR Model with half-normal plateau detection function

1.1 Habitat and trapping grid

For this example, we create a \(30 \times 30\) habitat grid represented by grid cell centers. We center a \(20 \times 20\) trapping grid in the habitat leaving an unsampled buffer width of 5 distance units (du) around it.

In order to use nimbleSCR’s efficient functions, we have to rescale habitat and trap coordinates to to allow a fast habitat cell ID lookup and the local evaluation (see Milleret et al. (2019) and Turek et al. (2021)
for further details). We use the ‘scaleCoordsToHabitatGrid’ function to rescale coordinates and the ‘getLocalObjects’ function to set up the objects necessary to perform the local evaluation. Special care should be taken when chosing ‘dmax’, the radius within which traps are considered, relative to \(\sigma\) (Milleret et al. 2019). Here, we are using a value \(>2 * \sigma\).

1.2 Define SCR model with half-normal plateau detection function

Here, we use a custom distribution ‘dbinomLocal_normalPlateau()’ (available in the nimbleSCR package) corresponding to the half-normal plateau detection function (Dey et al. 2022).

Note that we do not provide ‘detNums’ and ‘detIndices’ arguments in the ‘dbinomLocal_normalPlateau()’ function as we want to use this model to simulate data (see ‘?dbinomLocal_normalPlateau’ for further details). When simulating detections using this formulation of the SCR model in NIMBLE, all the information about detections (where and how many) is stored for each indiidual in the observation vector ‘y’ in the following order:

We now need to provide the maximum number of spatial recaptures that can be simulated per individual. We recommend using ‘trapLocal$numlocalindicesmax’ that defines the maximum number of traps available for detections when local evaluation is used. This will enable the simulation of as many spatial detections as allowed by the restrictions imposed by the local evaluation (defined by the ‘dmax’ argument from ‘getLocalObjects’). This means that the length of the ‘y’ observation vector for each individual is equal to the length of \(c(detNums, x, detIndices)\) and is therefore equal to \(lengthYCombined = 1+ 2 * trapLocal\$numLocalIndicesMax\).

This means that: ‘y[1]’ denotes the number of detections, ‘y[2:(trapLocal$numLocalIndicesMax+1)]’ gives the number of detections in each detector where it was detected (and ‘-1’ at the remaining indices) and ‘y[(trapLocal$numLocalIndicesMax+2):(2*trapLocal$numLocalIndicesMax+1)]’ gives the id of the traps at which detections occur (and ‘-1’ at the remaining indices).

The nimble model code of an SCR model with half-normal plateau detection function is given below.

1.3 Set parameter values to simulate

Here, we will use nimble’s functionalities to simulate a SCR data set directly from the model. For this, we need to set the values of the top-level parameters of the model:

The model formulation uses data augmentation to derive N estimates (Royle and Dorazio 2012). We therefore need to choose the total number of individuals M (detected + augmented). Here we used 500. The expected total number of individuals present in the population is ‘M * psi’.

1.4 Create data, constants and inits objects

We can now create the lists of input data and constants required by NIMBLE to build the model.

In order to simulate data from the nimble model, we need to provide simulated values as initial values.

1.5 Create NIMBLE model

we can then build the nimble SCR model object:

1.6 Simulate SCR data from the NIMBLE model

We first need to obtain the list of nodes that will be simulated. We used the ‘getDependencies’ function from NIMBLE. Using the ‘simulate’ function from NIMBLE, we will then simulate the activity center (AC) locations (‘sxy’), the state of the individual (‘z’) and SCR observation data (‘y’) given the values we provided for ‘p0’, ‘sigma’, ‘w’ and ‘psi’.

After running ‘simulate’, the simulated data are stored in the ‘model’ object. For example, we can access the simulated ‘z’ and check how many individuals were considered present:

We have simulated 257 individuals present in the population of which 162 were detected.

1.7. MCMC with NIMBLE

1.7.1 Compile the nimble model and create MCMC configuration

Here, we build the NIMBLE model again using the simulated ‘y’ as data. For simplicity, we used the simulated ‘z’ as initial values. Then we can fit the SCR model with the simulated ‘y’ data set.

## [1] -5721.569
## ===== Monitors =====
## thin = 1: N, p0, psi, sigma, w
## ===== Samplers =====
## binary sampler (500)
##   - z[]  (500 elements)
## RW_block sampler (500)
##   - sxy[]  (500 multivariate elements)
## RW sampler (4)
##   - psi
##   - sigma
##   - w
##   - p0
## [[1]]
## RW sampler: psi,  reflective: TRUE
## [[2]]
## RW sampler: sigma,  reflective: TRUE
## [[3]]
## RW sampler: w,  reflective: TRUE
## [[4]]
## RW sampler: p0,  reflective: TRUE
## [[5]]
## binary sampler: z[1],  reflective: TRUE

1.8. Computation of home range size

Estimating home-range size for a given circular detection function is equivalent to finding an estimate of the quantile \(r_\alpha\) such that \(\alpha\%\) of all movements lie within the circle of radius \(r_\alpha\) centered on activity centre \(\mathbf{\textit{s}}\). We can then define the \(\alpha\%\) home range area as \(A_\alpha\), the set of all points \(\mathbf{\textit{x}}\) in the habitat such that \(||\mathbf{\textit{x}} − \mathbf{\textit{s}}|| ≤ r_\alpha\).

While an analytical estimate can be derived for the classical half-normal detection function (see below), the half-normal plateau is a custom detection function for which an analytical estimate \(r_\alpha\) is not readily available.

The half-normal plateau detection function is written as: \[\begin{align} p_{\text{HNP}}(d) = \begin{cases} p_0,& \text{if } d < w,\\[-0.5em] p_0 \, \exp\Big{(}-\frac{(d-w)^2}{2\sigma^2}\Big{)}, & \text{if } d \geq w, \end{cases} \end{align}\] where \(p_0 \in (0,1)\), \(\sigma > 0\) and \(w \geq 0\). Here \(d\) can be interpreted as the distance between a given location and activity centre of an animal. The parameters in the function: \(p_0\) denotes the baseline detection probability, \(\sigma\) denotes the scale parameter, and \(w\) denotes the plateau length where the detection probability remains constant. For a better understanding, we also plot the detection function against distance in below.

## HALF-NORMAL PLATEAU FUNCTION
detfun.HNP <- function(D, p0, sigma, w){
  bool <- D <= w
  which_not_bool <- which(!bool)
  dens <- numeric(length = length(D))  
  dens[bool] <- 1
  adjustedD2 <- (D[which_not_bool]-w)^2
  dens[which_not_bool] <- exp(-(adjustedD2)/(2.0*sigma*sigma))
  return(p0*dens)
}
par(mfrow=c(1,1))
x <- seq(0, 15, len=10001)
plot(x, detfun.HNP(D=x, p0=0.2, sigma = 1, w = 1), type = "l", lwd=2, col = "orange",
     main = "Half-normal plateau detection function", xlab = 'd', ylab = 'Detection probability') 
lines(x, detfun.HNP(D=x, p0=0.1, sigma = 1.5, w = 3), lwd=2, lty = "solid", col = "darkcyan")
focal.vars <- c('p0', 'sigma', 'w') 
focal.values <- cbind(c(0.2, 0.1), # p0
                      c(1, 1.5), # sigma
                      c(1, 3)) # w 
legend.items <- apply(focal.values, 1, 
                      function(s){paste(paste(focal.vars,"=",s,sep=" "),collapse=", ")})
legend("topright",
       legend = legend.items,
       lty = c('solid', 'solid'), 
       bg="transparent",
       col=c('orange', "darkcyan"),
       cex = 0.75,
       bty = 'n')

The function ‘getHomeRangeArea()’ uses numerical approximation (to find the root) and estimate of \(95\%\) home range radius and area for a set of circular detection functions available in nimbleSCR, viz., half-normal (detFun = 0), half-normal plateau (detFun = 1), exponential (detFun = 2), asymmetric logistic (detFun = 3), bimodal (detFun = 4), and donut (detFun = 5).

The function ‘getHomeRangeArea()’ has two parts: a setup and and a run part. If we only want to calculate HR radius for a single set of parameter values, there is no need to compile. But for multiple sets of parameter values (e.g., an mcmc sample) one should compile the function before implementing the run part to reduce the computation time. ‘getHomeRangeArea()’ returns the estimates of home range radius and area for a given set of arguments with respect to a particular detection function using bisection algorithm:

The run part of the function can be run without specifying any arguments. For more information on the detection functions, please refer to Dey et al. (2022). We provide some examples to compute home range radius and area in below.

Next we obtain estimate \(95\%\) home range radius and area for half-normal plateau detection function (detFun = 1). Note that, we only have to provide the parameter values for ‘sigma’ and ‘w’ (i.e., ‘p0’ is not needed) to compute home range radius and area.

prob <- 0.95
paramnames.hr <- c("HRradius", "HRarea")
params <- c(sigma, w)
names(params) <- c("sigma", "w")
HRAnim <- getHomeRangeArea( x = params, detFun = 1, prob = prob, d = 6, 
                      xlim = c(0, dim(habitatMask)[2]), ylim = c(0, dim(habitatMask)[1]),
                      nBreaks = 800, tol = 1E-5, nIter = 2000)

# Different values of argument "detFun"
# 0 = Half-normal, 1 = Half-normal plateau, 2 = Exponential,
# 3 = Aysmmetric logistic, 4 = Bimodal, 5 = Donut.
HR.hnp <- c(HRAnim$run())
names(HR.hnp) <- paramnames.hr
print(HR.hnp)
##  HRradius    HRarea 
##  3.546511 39.514144

Simulated values of \(95\%\) home range radius and area are 3.55 du and 39.51 du\(^2\), respectively.

We can also obtain estimates of home range size and its associated uncertainty from the MCMC chains of the parameters \(p_0, \sigma\), and \(w\).

myNimbleOutput <- as.mcmc.list(myNimbleOutput) 
if(is.list(myNimbleOutput)){posteriorSamples <- do.call(rbind, myNimbleOutput)} 
if(!is.list(myNimbleOutput)){posteriorSamples <- as.matrix(myNimbleOutput)} 
theseIter <- round(seq(1, nrow(posteriorSamples), by = 10))
posteriorSamples <- posteriorSamples[theseIter,names(params)]
HRAnim.mat <- getHomeRangeArea(x = posteriorSamples, detFun = 1, prob = prob, d = 6, 
                         xlim = c(0, dim(habitatMask)[2]), ylim = c(0, dim(habitatMask)[1]),
                         nBreaks = 800, tol = 1E-5, nIter = 2000)

cHRAnim.arr <- compileNimble(HRAnim.mat, resetFunctions = TRUE)

HRA.Runtime <- system.time(
  HR.chain <- cHRAnim.arr$run()
)
print(HRA.Runtime)
##    user  system elapsed 
##  49.364   0.096  49.603
dimnames(HR.chain)[[2]] <- paramnames.hr

From the obtained sample ‘HR.chain’ of \(95\%\) home range radius and area, we can obtain summary statistics of these two metrics.

HRest <- do.call(rbind, lapply(c(1:2), function(j){
  c(mean(HR.chain[,j], na.rm = T), sd(HR.chain[,j], na.rm = T), quantile(HR.chain[,j], probs = c(0.025, 0.975), na.rm = T))
}))
dimnames(HRest) <- list(paramnames.hr, c("MEAN", "SD", "2.5%", "97.5%"))

cat("Numerical estimates using MCMC samples: \n", sep = "")
## Numerical estimates using MCMC samples:
print(HRest)
##               MEAN         SD      2.5%     97.5%
## HRradius  3.558721 0.09277004  3.394315  3.756545
## HRarea   39.813718 2.08249689 36.195473 44.332986

Based on the above analysis, the mean estimate of the \(95\%\) home range radius is 3.56 du and the \(95\%\) home range area is 39.81 du\(^2\).

For visual representation of the distributions of \(95\%\) home range radius and area, we also get the traceplots and density of these metrics.

HR.mcmc.sample <- as.mcmc.list(lapply(1:nchains, function(x) {
  niterPerChain <- nrow(HR.chain) / nchains
  theseRows <- ((x - 1) * niterPerChain + 1):(x * niterPerChain)
  this.MCMC.chain <- mcmc(HR.chain[theseRows, ])
  return(this.MCMC.chain)
})) 
names(HR.mcmc.sample) <- paste0("chain", 1:nchains)
chainsPlot(HR.mcmc.sample, line = c(HR.hnp))

The numerical estimate of home range radius and area can be improved further by tuning the function arguments: the number of bins (‘nBreaks’), the tolerance (‘tol’) and the uppler limit of the number of iterations for numerical approximation using bi-section algorithm (‘nIter’). However, tuning of these parameters will also affect the computation time of getHomeRangeArea().

2. Calculating home range radius and area with other detection functions

We can also fit SCR models with half-normal or exponential detection function using custom distributions such as ‘dbinomLocal_normal()’ or ‘dbinomLocal_exp()’, respectively. An example of SCR model fitting using half-normal detection and other efficient nimbleSCR features is given here. Fitting of SCR model using exponential detection function can be performed following similar steps (see ‘?dbinomLocal_exp’ for further details). Below, we provide examples to compute home range radius and area for the other detection functions available (Dey et al. 2022).

Half-normal (detFun = 0)

The half-normal detection function is written as: \[\begin{align} p_{\text{HN}}(d) = p_0 \, \exp\Big{(}-\frac{d^2}{2\sigma^2}\Big{)}, \end{align}\] where \(p_{0} \in (0,1)\), \(\sigma > 0\). Here \(p_0\) denotes the baseline detection probability and \(\sigma\) denotes the scale parameter. For a better understanding, we also plot the detection function against distance.

Next we obtain the \(95\%\) home range radius and area for half-normal detection function. Note that we only have to provide the parameter values for ‘sigma’ (i.e., ‘p0’ is not needed) to compute home range radius and area.

##  HRradius    HRarea 
##  4.897252 75.345054

Exponential (detFun = 2)

The exponential detection function is written as: \[\begin{align} p_{\text{EXP}}(d) = p_0 \, \exp\Big{(}-\text{rate} \, * \, d\Big{)}, \end{align}\] where \(p_{0} \in (0,1)\), \(\text{rate} > 0\). Here \(p_0\) denotes the baseline detection probability and \(\text{rate}\) denotes the rate parameter. For a better understanding, we also plot the detection function against distance.

Next we obtain the \(95\%\) home range radius and area for exponential detection function. Again, we only have to provide the parameter values for ‘rate’ (i.e., ‘p0’ is not needed) to compute home range radius and area.

##   HRradius     HRarea 
##   9.389359 276.963004

Asymmetric logistic (detFun = 3)

The asymmetric logistic detection function is written as: \[\begin{align} p_{\text{AL}}(d) = p_0 \, \Big{[} 1 + (d/\sigma)^{\alpha_a} g(d) \, + (d/\sigma)^{\alpha_a \alpha_b} (1-g(d)) \Big{]}^{-1}, \end{align}\] where \(g(d) =\big{\{} 1 + (d/\sigma)^{\nu} \big{\}}^{-1}\), \(\nu = \frac{2|\alpha_a|\alpha_b}{1+\alpha_b}\), \(p_{0} \in (0,1)\), \(\sigma > 0\), \(\alpha_a \in \mathbb{R}\) and \(\alpha_b > 0\). Here \(p_0\) denotes the baseline detection probability, \(\alpha_a\) denotes the first curvature parameter, \(\alpha_b\) denotes the second curvature parameter and \(\sigma\) denotes the distnace from the activity centre where the asymmetric logistic detection function takes the value \(p_0/2\). For a better understanding, we also plot the detection function against distance.

Next we obtain the \(95\%\) home range radius and area for asymmetric logistic detection function. Note that we only have to provide the parameter values for ‘sigma’, ‘alpha.a’ and ‘alpha.b’ (i.e., ‘p0’ is not needed) to compute home range radius and area.

##  HRradius    HRarea 
##  4.257453 56.944218

Bimodal (detFun = 4)

The bimodal detection function is written as: \[\begin{align} p_{\text{BI}}(d) = p_{0a} \, \exp\Big{(}-\frac{d^2}{2\sigma_{a}^2}\Big{)} \, + p_{0b} \, \exp\Big{(}-\frac{(d-w)^2}{2\sigma_{b}^2}\Big{)} \end{align}\] where \(p_{0a} \in (0,1)\), \(p_{0b} \in (0,1)\), \(\sigma_{a} > 0\), \(\sigma_{b} > 0\) and \(w \geq 0\). Here \(p_{0a}\) denotes the baseline detection probability of the first peak, \(p_{0b}\) denotes the baseline detection probability of the second peak, \(\sigma_{a}\) denotes the scale parameter for the left tail, \(\sigma_{b}\) denotes the scale parameter for the right tail and \(w\) denotes the distance between the two peaks. For a better understanding, we also plot the detection function against distance.

Next we obtain the \(95\%\) home range radius and area for bimodal detection function. Here, we have to provide all the parameter values for ‘p0.a’, ‘sigma.a’, ‘p0.b’, ‘sigma.b’ and ‘w’ to compute home range radius and area.

##  HRradius    HRarea 
##  3.973037 49.590105

Donut (detFun = 5)

The donut detection function is written as: \[\begin{align} p_{\text{DN}}(d) = \begin{cases} p_0 \, \exp\Big{(}-\frac{(d-w)^2}{2\sigma_{a}^2}\Big{)},& \text{if } d < w,\\[-0.5em] p_0 \, \exp\Big{(}-\frac{(d-w)^2}{2\sigma_{b}^2}\Big{)}, & \text{if } d \geq w, \end{cases} \end{align}\] where \(p_0 \in (0,1)\), \(\sigma_{a} > 0\), \(\sigma_{b} > 0\) and \(w \geq 0\). Here \(p_0\) denotes the baseline detection probability, \(\sigma_{a}\) denotes the scale parameter for the left tail, \(\sigma_{b}\) denotes the scale parameter for the right tail and \(w\) denotes the distance between the home-range center and the “donut” centre. For a better understanding, we also plot the detection function against distance.

Next we obtain the \(95\%\) home range radius and area for donut detection function. Note that we only have to provide the parameter values for ‘sigma.a’, ‘sigma.b’ and ‘w’ (i.e., ‘p0’ is not needed) to compute home range radius or area.

##  HRradius    HRarea 
##  3.119946 30.580466

REMARK

The analytical estimate of home range radius for the half-normal function is: \(\hat{r}_\alpha=\sigma \sqrt{q(\alpha,2)}\), where \(q(\alpha,2)\) is the \(\alpha\%\) quantile of a chi-square distribution with 2 degrees of freedom (Royle et al. 2014).

sigma <- 2
q<-qchisq(prob,2)
radius<-sigma*sqrt(q)
area<-pi*(radius^2)
HR.true<-c(radius,area)
names(HR.true) <- paramnames.hr
cat("Analytical estimates: \n", sep = "")
## Analytical estimates:
print(HR.true)
##  HRradius    HRarea 
##  4.895494 75.290964

Note that, the estimates of \(95\%\) home range radius and area obtained from our method using numerically approximation (4.9, 75.35) are very close to the above analytical estimate (4.9, 75.29).

REFERENCES

Bischof, Richard, Daniel Turek, Cyril Milleret, Torbjørn Ergon, Pierre Dupont, and Perry de Valpine. 2020. NimbleSCR: Spatial Capture-Recapture (Scr) Methods Using ’Nimble’. https://CRAN.R-project.org/package=nimbleSCR.

Dey, Soumen, Richard Bischof, Pierre P. A. Dupont, and Cyril Milleret. 2022. “Does the Punishment Fit the Crime? Consequences and Diagnosis of Misspecified Detection Functions in Bayesian Spatial Capture–Recapture Modeling.” Ecology and Evolution 12 (2): e8600.

Milleret, Cyril, Pierre Dupont, Christophe Bonenfant, Henrik Brøseth, Øystein Flagstad, Chris Sutherland, and Richard Bischof. 2019. “A Local Evaluation of the Individual State-Space to Scale up Bayesian Spatial Capture–Recapture.” Ecology and Evolution 9 (1): 352–63.

NIMBLE Development Team. 2021. NIMBLE: MCMC, Particle Filtering, and Programmable Hierarchical Modeling (Verison 0.12.10. https://doi.org/10.5281/zenodo.5562925.

Royle, J. Andrew, Richard B. Chandler, Rahel Sollmann, and Beth Gardner. 2014. Spatial Capture–Recapture. Academic Press.

Royle, J. Andrew, and Robert M. Dorazio. 2012. “Parameter-Expanded Data Augmentation for Bayesian Analysis of Capture–Recapture Models.” Journal of Ornithology 152 (2): 521–37.

Turek, Daniel, Cyril Milleret, Torbjørn Ergon, Henrik Brøseth, Pierre Dupont, Richard Bischof, and Perry De Valpine. 2021. “Efficient Estimation of Large-Scale Spatial Capture–Recapture Models.” Ecosphere 12 (2): e03385.

Valpine, Perry de, Daniel Turek, Christopher J. Paciorek, Clifford Anderson-Bergman, Duncan Temple Lang, and Rastislav Bodik. 2017. “Programming with Models: Writing Statistical Algorithms for General Model Structures with NIMBLE.” Journal of Computational and Graphical Statistics 26 (2): 403–13.