wolverineCode.Rmd

---
title: "Using nimbleSCR to fit point process SCR models to wolverine non-invasive genetic sampling data"
author: "Pierre Dupont"
date: "`r Sys.Date()`"
output: rmarkdown::html_vignette
bibliography: references.bib  
vignette: >
  %\VignetteIndexEntry{Wolverine SCR_pp Example}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(echo = TRUE)
knitr::opts_chunk$set(fig.width = 5, fig.height = 5) 
```

```{r, echo = FALSE}
## Load wolverine data
load("wolverineData.RData")
```

This vignette demonstrates the different wolverine SCR models presented in
"A flexible and efficient Bayesian implementation of point process models for 
spatial capture-recapture data" [@Zhang2022] using the NIMBLE [@valpine2017nimble;@nimbleSoftware2021] and nimbleSCR packages [@nimbleSCR].


### Load Libraries 
```{r, warning = FALSE, message = FALSE}
library(nimble)
library(nimbleSCR)
library(basicMCMCplots)
library(coda)
```


### Define `nimble` Models

Here, we define the first `nimble` model.
It uses i) the Bernoulli point process for modelling the distribution of individual wolverine activity centers, ii) a Poisson point process for modelling individual detections, and iii) a semi-complete data-likelihood approach [@King2016] to estimate model parameters.
```{r}
modelCode1 <- nimbleCode({  
  ##------ SPATIAL PROCESS 
  ## Intercept and slope for the log-linear model for habitat selection intensity
  habCoeffInt ~ dnorm(0, sd = 10)
  habCoeffSlope ~ dnorm(0, sd = 10)
  ## Habitat intensity for each habitat window
  habIntensity[1:numHabWindows] <- exp(habCoeffInt + habCoeffSlope * habCovs[1:numHabWindows])
  sumHabIntensity <- sum(habIntensity[1:numHabWindows])
  logHabIntensity[1:numHabWindows] <- log(habIntensity[1:numHabWindows])
  logSumHabIntensity <- log(sumHabIntensity)
  ## Activity centres of the observed individuals: a Bernoulli point process
  for(i in 1:numIdDetected){
    sxy[i,1:2] ~ dbernppAC(
      lowerCoords = habLoCoords[1:numHabWindows, 1:2],
      upperCoords = habUpCoords[1:numHabWindows, 1:2],
      logIntensities = logHabIntensity[1:numHabWindows],
      logSumIntensity = logSumHabIntensity,
      habitatGrid = habitatGrid[1:y.max,1:x.max],
      numGridRows = y.max,
      numGridCols = x.max)
  }#i
  
  ##----- DEMOGRAPHIC PROCESS
  ## Number of individuals in the population
  N ~ dpois(sumHabIntensity)
  ## Number of detected individuals
  nDetectedIndiv ~ dbin(probDetection, N)

  ##----- DETECTION PROCESS 
  ## Scale for the multivariate normal detection function
  sigma ~ dunif(0,10)
  ## Intercept and slope parameters for the log-linear model for detection intensity
  for(c in 1:numCounties){
    detCoeffInt[c] ~ dnorm(0, sd = 10)
  }#c
  for(cc in 1:numDetCovs){
    detCoeffSlope[cc] ~ dnorm(0, sd = 10)
  }#cc
  ## Baseline detection intensity for each detection window
  for(j in 1:numDetWindows){
    detIntensity[j] <- exp( detCoeffInt[detCounties[j]] +
                              detCoeffSlope[1] * detCovs[j,1] +
                              detCoeffSlope[2] * detCovs[j,2] +
                              detCoeffSlope[3] * detCovs[j,3])
  }#j
  ## Detections of the observed individuals conditional on their activity centers
  for(i in 1:numIdDetected) {
    y[i,1:(maxDetections+1),1:3] ~ dpoisppDetection_normal(
      lowerCoords = detLoCoords[1:numDetWindows,1:2],
      upperCoords = detUpCoords[1:numDetWindows,1:2],
      s = sxy[i,1:2],
      sd = sigma,
      baseIntensities = detIntensity[1:numDetWindows],
      numMaxPoints = maxDetections,
      numWindows = numDetWindows,
      indicator = 1)
  }#i

  ## The probability that an individual in the population is detected at least once
  ## i.e. one minus the void probability over all detection windows
  probDetection <- 1 - marginalVoidProbNumIntegration(
    quadNodes = quadNodes[1:maxNumNodes,1:2,1:numHabWindows],
    quadWeights = quadWeights[1:numHabWindows],
    numNodes = numNodes[1:numHabWindows],
    lowerCoords = detLoCoords[1:numDetWindows,1:2],
    upperCoords = detUpCoords[1:numDetWindows,1:2],
    sd = sigma,
    baseIntensities = detIntensity[1:numDetWindows],
    habIntensities = habIntensity[1:numHabWindows],
    sumHabIntensity = sumHabIntensity,
    numObsWindows = numDetWindows,
    numHabWindows = numHabWindows
  )
  logDetProb <- log(probDetection)
  normData ~ dnormalizer(logNormConstant = -numIdDetected*logDetProb)

})
```


Next, we can define a second `nimble` model based on the "usual" formulation of Bayesian SCR models.
In this model, individual activity center locations are modeled using a two-step process in which i) a categorical distribution is used to sample the habitat window the activity center is located in, and ii) the center of the selected habitat window is used as the location of the activity center inside the selected habitat window. The probability vector of the categorical process follows the intensity of the density point process of the first formulation. Individual detections are not modeled continuously in space, instead they are aggregated at discrete detector locations and the number of detections at each detector is modeled using a Poisson distribution.
Finally, we use the data augmentation approach [@royle2012parameter] to derive population size instead of the semi-complete data-likelihood approach.
```{r}
modelCode2 <- nimbleCode({
  ##------ SPATIAL PROCESS
  habCoeffSlope ~ dnorm(0, sd = 10)
  ## Habitat intensity for each habitat window
  habIntensity[1:numHabWindows] <- exp(habCoeffSlope * habCovs[1:numHabWindows])
  for(i in 1:M){
    sID[i] ~ dcat(habIntensity[1:numHabWindows])
  }#i

  ##----- DEMOGRAPHIC PROCESS
  psi ~ dunif(0,1)
  for(i in 1:M){
    z[i] ~ dbern(psi)
  }#i
  ## Number of individuals in the population
  N <- sum(z[1:M])

  ##----- DETECTION PROCESS
  ## Scale for the multivariate normal detection function
  sigma ~ dunif(0,10)
  ## Intercept and slope parameters for the log-linear model for detection intensity
  for(c in 1:numCounties){
    detCoeffInt[c] ~ dnorm(0, sd = 10)
  }#c
  for(cc in 1:numDetCovs){
    detCoeffSlope[cc] ~ dnorm(0, sd = 10)
  }#cc
  ## Baseline detection intensity for each detection window
  for(j in 1:numDetWindows){
    lambdaTraps[j] <- exp( detCoeffInt[detCounties[j]] +
                             detCoeffSlope[1] * detCovs[j,1] +
                             detCoeffSlope[2] * detCovs[j,2] +
                             detCoeffSlope[3] * detCovs[j,3])
  }#j

  ## Detections of the observed individuals conditional on their activity centers
  for(i in 1:M) {
    y[i,1:lengthYCombined] ~ dpoisLocal_normal(
      lambdaTraps = lambdaTraps[1:numDetWindows],
      trapCoords = detCoords[1:numDetWindows,1:2],
      s = habCoords[sID[i],1:2],
      sigma = sigma,
      localTrapsIndices = localTrapsIndices[1:numHabWindows,1:numDetWindows],
      localTrapsNum = localTrapsNum[1:numHabWindows],
      habitatGrid = habitatGrid[1:y.max,1:x.max],
      lengthYCombined = lengthYCombined,
      indicator = z[i])
  }#i
})
```


For comparison purposes, we define a third `nimble` model, which is a hybrid formulation of the two previous ones. It uses i) the Bernoulli point process for modelling the distribution of wolverine activity centers and ii) a Poisson point process to model individual detections in continuous space but uses iii) data augmentation instead of the semi-complete data-likelihood approach to estimate population size. This formulation is expected to be both a better representation of the detection process (thus more accurate than model 2) and more efficient than model 1 because of the complexity of the semi-complete data likelihood calculation.
```{r}
modelCode3 <- nimbleCode({
  ##------ SPATIAL PROCESS
  habCoeffSlope ~ dnorm(0, sd = 10)
  ## Habitat intensity for each habitat window
  habIntensity[1:numHabWindows] <- exp(habCoeffSlope * habCovs[1:numHabWindows])
  sumHabIntensity <- sum(habIntensity[1:numHabWindows])
  logHabIntensity[1:numHabWindows] <- log(habIntensity[1:numHabWindows])
  logSumHabIntensity <- log(sumHabIntensity)
  ## Activity centres of the observed individuals: a bernoulli point process
  for(i in 1:M){
    sxy[i,1:2] ~ dbernppAC(
      lowerCoords = habLoCoords[1:numHabWindows,1:2],
      upperCoords = habUpCoords[1:numHabWindows,1:2],
      logIntensities = logHabIntensity[1:numHabWindows],
      logSumIntensity = logSumHabIntensity,
      habitatGrid = habitatGrid[1:y.max,1:x.max],
      numGridRows = y.max,
      numGridCols = x.max)
  }#i

  ##----- DEMOGRAPHIC PROCESS
  psi ~ dunif(0,1)
  for(i in 1:M){
    z[i] ~ dbern(psi)
  }#i
  ## Number of individuals in the population
  N <- sum(z[1:M])

  ##----- DETECTION PROCESS
  ## Scale for the half-normal detection function
  sigma ~ dunif(0,10)
  ## Intercept and slope parameters for the log-linear model for detection intensity
  for(c in 1:numCounties){
    detCoeffInt[c] ~ dnorm(0, sd = 10)
  }#c
  for(cc in 1:numDetCovs){
    detCoeffSlope[cc] ~ dnorm(0, sd = 10)
  }#cc
  ## Baseline detection intensity for each detection window
  for(j in 1:numDetWindows){
    detIntensity[j] <- exp( detCoeffInt[detCounties[j]] +
                              detCoeffSlope[1] * detCovs[j,1] +
                              detCoeffSlope[2] * detCovs[j,2] +
                              detCoeffSlope[3] * detCovs[j,3])
  }#j

  ## Detections of the observed individuals conditional on their activity centers
  for(i in 1:M){
    y[i,1:(maxDetections+1),1:3] ~ dpoisppDetection_normal(
      lowerCoords = detLoCoords[1:numDetWindows,1:2],
      upperCoords = detUpCoords[1:numDetWindows,1:2],
      s = sxy[i,1:2],
      sd = sigma,
      baseIntensities = detIntensity[1:numDetWindows],
      numMaxPoints = maxDetections,
      numWindows = numDetWindows,
      indicator = z[i])
  }#i
})
```


We can also define a fourth `nimble` model, which uses the categorical distribution to model the distribution of wolverine activity centers, a Poisson point process to model individual detections in continuous space and data augmentation to estimate population size. This formulation is expected to be slower than version 3 because of the categorical process. 
```{r}
modelCode4 <- nimbleCode({
  ##------ SPATIAL PROCESS
  habCoeffSlope ~ dnorm(0, sd = 10)
  ## Habitat intensity for each habitat window
  habIntensity[1:numHabWindows] <- exp(habCoeffSlope * habCovs[1:numHabWindows])
  for(i in 1:M){
    sID[i] ~ dcat(habIntensity[1:numHabWindows])
  }#i

  ##----- DEMOGRAPHIC PROCESS
  psi ~ dunif(0,1)
  for(i in 1:M){
    z[i] ~ dbern(psi)
  }#i
  ## Number of individuals in the population
  N <- sum(z[1:M])

  ##----- DETECTION PROCESS
  ## Scale for the half-normal detection function
  sigma ~ dunif(0,10)
  ## Intercept and slope parameters for the log-linear model for detection intensity
  for(c in 1:numCounties){
    detCoeffInt[c] ~ dnorm(0, sd = 10)
  }#c
  for(cc in 1:numDetCovs){
    detCoeffSlope[cc] ~ dnorm(0, sd = 10)
  }#cc
  ## Baseline detection intensity for each detection window
  for(j in 1:numDetWindows){
    detIntensity[j] <- exp( detCoeffInt[detCounties[j]] +
                              detCoeffSlope[1] * detCovs[j,1] +
                              detCoeffSlope[2] * detCovs[j,2] +
                              detCoeffSlope[3] * detCovs[j,3])
  }#j

  ## Detections of the observed individuals conditional on their activity centers
  for(i in 1:M){
    y[i,1:(maxDetections+1),1:3] ~ dpoisppDetection_normal(
      lowerCoords = detLoCoords[1:numDetWindows,1:2],
      upperCoords = detUpCoords[1:numDetWindows,1:2],
      s = habCoords[sID[i],1:2],
      sd = sigma,
      baseIntensities = detIntensity[1:numDetWindows],
      numMaxPoints = maxDetections,
      numWindows = numDetWindows,
      indicator = z[i])
  }#i
})
```



### Load Wolverine Data

We load the wolverine example data directly. More details on how the data were extracted can be found in Section S1 of Appendix S2 of @Zhang2022.  
```{r, eval = FALSE}
load("wolverineData.RData")
```


## Create NIMBLE Model

Now we create `nimble` model objects, using the model code
defined above, as well as constants, data, and initial values.
```{r eval = FALSE}
Rmodel1 <- nimbleModel(modelCode1, constants1, data1, inits1)
Rmodel2 <- nimbleModel(modelCode2, constants2, data2, inits2)
Rmodel3 <- nimbleModel(modelCode3, constants3, data3, inits3)
Rmodel4 <- nimbleModel(modelCode4, constants4, data4, inits4)
```


#### Configure and Build MCMC

We configure an MCMC algorithm for each `Rmodel` model object and we assign MCMC monitors to the different parameters we want to track (e.g. $N$, $\sigma$, ...)
```{r message = FALSE, eval = FALSE}
conf1 <- configureMCMC( Rmodel1,
                       monitors = params1,
                       print = FALSE)
conf2 <- configureMCMC( Rmodel2,
                       monitors = params2,
                       print = FALSE)
conf3 <- configureMCMC( Rmodel3,
                       monitors = params3,
                       print = FALSE)
conf4 <- configureMCMC( Rmodel4,
                       monitors = params4,
                       print = FALSE)
Rmcmc1 <- buildMCMC(conf1)
Rmcmc2 <- buildMCMC(conf2)
Rmcmc3 <- buildMCMC(conf3)
Rmcmc4 <- buildMCMC(conf4)
```



#### Compile and Run MCMC

Finally, we compile the models and MCMC objects and
execute the compiled MCMC algorithms for 4 chains of 11000 iterations each.
```{r eval = FALSE}
Cmodel1 <- compileNimble(Rmodel1)
Cmcmc1 <- compileNimble(Rmcmc1, project = Rmodel1)
MCMC_runtime1 <- system.time(
  MCMC_samples1 <- runMCMC(Cmcmc1,
                           niter = 11000,
                           nburnin = 1000,
                           nchains = 4))

Cmodel2 <- compileNimble(Rmodel2)
Cmcmc2 <- compileNimble(Rmcmc2, project = Rmodel2)
MCMC_runtime2 <- system.time(
  MCMC_samples2 <- runMCMC(Cmcmc2,  
                           niter = 11000,
                           nburnin = 1000,
                           nchains = 4))

Cmodel3 <- compileNimble(Rmodel3)
Cmcmc3 <- compileNimble(Rmcmc3, project = Rmodel3)
MCMC_runtime3 <- system.time(
  MCMC_samples3 <- runMCMC(Cmcmc3,  
                           niter = 11000,
                           nburnin = 1000,
                           nchains = 4))

Cmodel4 <- compileNimble(Rmodel4)
Cmcmc4 <- compileNimble(Rmcmc4, project = Rmodel4)
MCMC_runtime4 <- system.time(
  MCMC_samples4 <- runMCMC(Cmcmc4,  
                           niter = 11000,
                           nburnin = 1000,
                           nchains = 4))

MCMC_samples1_summary <- samplesSummary(do.call(rbind,MCMC_samples1))
MCMC_samples2_summary <- samplesSummary(do.call(rbind,MCMC_samples2))
MCMC_samples3_summary <- samplesSummary(do.call(rbind,MCMC_samples3))
MCMC_samples4_summary <- samplesSummary(do.call(rbind,MCMC_samples4))

MCMC_samples1_ess <- effectiveSize(MCMC_samples1)
MCMC_samples2_ess <- effectiveSize(MCMC_samples2)
MCMC_samples3_ess <- effectiveSize(MCMC_samples3)
MCMC_samples4_ess <- effectiveSize(MCMC_samples4)
```

```{r echo = FALSE, eval = FALSE}
##save(MCMC_samples1, MCMC_samples2, MCMC_samples3, MCMC_samples4,
##     MCMC_runtime1, MCMC_runtime2, MCMC_runtime3, MCMC_runtime4,
##     file = file.path(baseDir,
##                      "nimbleSCR/inst/extdata/WolverinePointProcess_samples.RData"))

# save(MCMC_samples1_summary, MCMC_samples2_summary, MCMC_samples3_summary, MCMC_samples4_summary,
#      MCMC_samples1_ess, MCMC_samples2_ess, MCMC_samples3_ess, MCMC_samples4_ess,
#      MCMC_runtime1, MCMC_runtime2, MCMC_runtime3, MCMC_runtime4,
#      file = file.path(baseDir,
#                       "nimbleSCR/inst/extdata/WolverinePointProcess_summary.RData"))
```

```{r echo = FALSE}
##if(is.null(baseDir)) {
##    load(system.file("extdata",
##                     "WolverinePointProcess_samples.RData",
##                     package = "nimbleSCR"))
##} else {
##    load(file.path(baseDir,"nimbleSCR/inst/extdata/WolverinePointProcess_samples.RData"))
##}
load("wolverineResults.RData")
```


## Results
We can then look at the summary statistics of posterior distributions for different models:
```{r }
round(MCMC_samples1_summary, 2)
round(MCMC_samples2_summary, 2)
round(MCMC_samples3_summary, 2)
round(MCMC_samples4_summary, 2)
```

<!--
And examine traceplots and posterior distributions for $N$ and $\sigma$:
-->
```{r echo = FALSE, eval = FALSE}
chainsPlot(MCMC_samples1, var = c("N","sigma"))
chainsPlot(MCMC_samples2, var = c("N","sigma"))
chainsPlot(MCMC_samples3, var = c("N","sigma"))
chainsPlot(MCMC_samples4, var = c("N","sigma"))
```

Next, we can check the posterior effective sample size (ESS) resulting from our
40000 posterior samples for the population size estimates ($N$):  
```{r }
round(MCMC_samples1_ess,2)[c("N","sigma")] 
round(MCMC_samples2_ess,2)[c("N","sigma")] 
round(MCMC_samples3_ess,2)[c("N","sigma")] 
round(MCMC_samples4_ess,2)[c("N","sigma")] 
```
We can also calculate the MCMC efficiency; this corresponds to the rate of generating effectively independent posterior samples, per second of MCMC runtime:
```{r }
round(MCMC_samples1_ess,2)["N"]/round(MCMC_runtime1[3],1)
round(MCMC_samples2_ess,2)["N"]/round(MCMC_runtime2[3],1)
round(MCMC_samples3_ess,2)["N"]/round(MCMC_runtime3[3],1)
round(MCMC_samples4_ess,2)["N"]/round(MCMC_runtime4[3],1)
```


## Conclusion
We can see that model 2 and 4 seem to overestimate individual space use ($\sigma$), and underestimate $N$, as a consequence of the coarse resolutions used for the habitat and detectors and the aggregation of activity centers and detections to the centers of the habitat and detector windows, respectively.

We can also see that the semi-likelihood approach is much slower than the data-augmentation for comparable models and results (Runtime : `r round(MCMC_runtime1[3]/3600,1)`hrs for model 1 against `r round(MCMC_runtime3[3]/3600,1)`hrs for model 3).

However, when looking at the efficiency, model 2 seems to perform better, thanks to the faster detection model using discrete locations (Efficiency : `r round(MCMC_samples2_ess,2)["N"]/round(MCMC_runtime2[3],1)`ESS/sec) (... but at the cost of poor population size estimates).

Overall, model 3, combining the detection point process and data augmentation seems to be the best compromise with accurate population size estimates and relatively good efficiency (Efficiency : `r round(MCMC_samples3_ess,2)["N"]/round(MCMC_runtime3[3],1)`ESS/sec) 



## References