Wolverine Example

1. Local detector evaluation

This step restricts calculations of the detection probabilities (here using a halfnormal function) to detectors (traps) within a radius where detections are realistically possible (Cyril Milleret et al. 2019). We use the function getLocalObjects to identify the set of detectors that are within a certain distance d_max of each habitat cell center (blue points in the plot below). These reduced sets of detectors are stored in the localIndices matrix and are later used in the local evaluation of the detection model to speed up calculations. The value of d_max should be as small as possible in order to reduce computation but always large enough so that for any particular individual, the set of local traps associated with the coordinates of the activity center s include all detectors at which that individual was detected. The d_max value will therefore affect the number of columns in localIndices Here, we use d_max = 38. with the coordinates of the activity center s … We also aggregated the habitat matrix to obtain larger habitat cells (lower resolution) and obtain objects with smaller dimensions. This reduces the number of habitat cells for which we have to identify the set of detectors that are within d_max of the cell center. The goal is to create the object localIndices of the smallest dimension possible, that balances the cost of looking up relevant grid cells and reducing calculations for each grid cell.

Here, we resize the habitat matrix by a factor of 24, which corresponds to the resizeFactor argument. This means that 24x24 cells are aggregated into a single cell. The resizeFactor value will affect how many rows localIndices will be composed of.

set.seed(2)

DetectorIndex <- getLocalObjects(habitatMask = data$habitat.mx,
                               coords = data$detector.xy,
                               dmax = 38,
                               resizeFactor = 24)

constants$y.maxDet <- dim(DetectorIndex$habitatGrid)[1]
constants$x.maxDet <- dim(DetectorIndex$habitatGrid)[2]
constants$ResizeFactor <- DetectorIndex$resizeFactor
constants$n.cells <- dim(DetectorIndex$localIndices)[1]
constants$maxNBDets <- DetectorIndex$numLocalIndicesMax
data$detectorIndex <- DetectorIndex$localIndices
data$nDetectors <- DetectorIndex$numLocalIndices
data$habitatIDDet <- DetectorIndex$habitatGrid

2. Sparse representation of the observation matrix

We re-express y as a sparse representation of the detection matrix to reduce its size. In this representation, we turn the detection matrix y into three objects:

• detIndices: where each row (corresponding to one individual) contains the identification numbers of detectors at which that individual was detected.
• y: a second matrix of identical dimension, containing the number of detections of a given individual at each detector. This second matrix is necessary for modelling non-binary detections (e.g. binomial observation models)
• detNums: a vector containing the number of detectors at which each individual was detected.

ySparse <- getSparseY(x = my.jags.input$y)
data$y <- ySparse$y[,,1]  
data$yDets <- ySparse$detIndices[,,1]
data$nbDetections <- ySparse$detNums[,1]
constants$nMaxDetectors <- ySparse$maxDetNums

3. Skip unnecessary calculations

The function dbinomLocal_normal() takes the logical argument indicator that specifies whether the individual i is available (z_i = 1) for detection or not (z_i = 0). When z_i = 0, calculations of p_ij are not performed and therefore increases MCMC efficiency.

Configure and Build MCMC

We configure an MCMC algorithm to the Rmodel model object.

We assign MCMC monitors to N, σ, and p₀.

Block sampling to increase MCMC efficiency

We also remove the univariate Metropolis-Hastings samplers that were assigned by default to each dimension of the sxy variables (ACs). Instead, we add joint Metropolis-Hastings samplers (RW_block) samplers to each x and y coordinate pair sxy[i, 1:2].

conf <- configureMCMC(Rmodel, monitors = c("N", "sigma", "p0"), print = FALSE)
conf$removeSamplers("sxy")
ACnodes <- paste0("sxy[", 1:constants$n.individuals, ", 1:2]")
for(node in ACnodes) {
  conf$addSampler(target = node,
                  type = "RW_block",
                  control = list(adaptScaleOnly = TRUE),
                  silent = TRUE)
}
Rmcmc <- buildMCMC(conf)

Compile and Run MCMC

Finally, we compile both the model and MCMC objects and execute the compiled MCMC for 10 000 iterations.

Cmodel <- compileNimble(Rmodel)
Cmcmc <- compileNimble(Rmcmc, project = Rmodel)
MCMC_runtime <- system.time(
  samples <- runMCMC(Cmcmc, niter = 10000)
)