inst/doc/Point_Process.R

In nimbleSCR: Spatial Capture-Recapture (SCR) Methods Using 'nimble'

## ----setup, include=FALSE-----------------------------------------------------
knitr::opts_chunk$set(echo = TRUE)
knitr::opts_chunk$set(fig.width = 8, fig.height = 10) 

## ----echo = FALSE-------------------------------------------------------------
##load("C:/Users/cymi/Downloads/WolverineData.RData")
#load("~/Downloads/WolverineData.RData")
## So it works for everyone
if(Sys.info()['user'] == 'dturek') {
  baseDir <- '~/github/nimble/nimbleSCR/'                   ## Daniel
} else if(Sys.info()['user'] == 'pidu') {
  baseDir <- 'C:/Users/pidu/PACKAGES/nimbleSCR/'            ## Pierre
} else if(Sys.info()['user'] == 'cymi') {
  baseDir <- 'C:/Personal_Cloud/OneDrive/Work/nimbleSCR/'   ## Cyril
} else if(Sys.info()['user'] == 'arichbi') {
  baseDir <- 'C:/PROJECTS/nimbleSCR/'                       ## Richard
} else if(Sys.info()['user'] == 'admin') {                  ## Soumen
  baseDir <- '~/GitHubSD/nimbleSCR/' 
} else baseDir <- NULL

if(is.null(baseDir)) {
    load(system.file("extdata", "WolverineData.RData", package = "nimbleSCR"))
} else {
    load(file.path(baseDir,"nimbleSCR/inst/extdata/WolverineData.RData"))
}

## ---- warning = FALSE, message = FALSE----------------------------------------
## Load packages
library(nimble)
library(nimbleSCR)
library(basicMCMCplots)

## ---- warning = FALSE, message = FALSE----------------------------------------
## Create habitat grid
coordsHabitatGridCenter <- cbind(rep(seq(75, 5, by = -10), 10),
                                 sort(rep(seq(5, 100, by = 10), 8)))
colnames(coordsHabitatGridCenter) <- c("x","y")

## Create trap grid
coordsObsCenter <- cbind(rep(seq(15, 65, by = 10), 8),
                         sort(rep(seq(15, 85, by = 10), 6)))
colnames(coordsObsCenter) <- c("x","y")

## Plot check
plot(coordsHabitatGridCenter[,"y"] ~ coordsHabitatGridCenter[,"x"],
     xlim = c(0,80), ylim = c(0,100),
     pch = 1, cex = 1.5) 
points(coordsObsCenter[,"y"] ~ coordsObsCenter[,"x"], col="red", pch=16 ) 
par(xpd=TRUE)
legend(x = 7, y = 13,
       legend=c("Habitat window centers", "Observation window centers"),
       pt.cex = c(1.5,1),
       horiz = T,
       pch=c(1,16),
       col=c("black", "red"),
       bty = 'n')


## ---- warning = FALSE, message = FALSE----------------------------------------
## Rescale coordinates
scaledObjects <- scaleCoordsToHabitatGrid(
  coordsData = coordsObsCenter,            
  coordsHabitatGridCenter = coordsHabitatGridCenter)

## Get lower and upper cell coordinates
lowerAndUpperCoords <- getWindowCoords(
  scaledHabGridCenter = scaledObjects$coordsHabitatGridCenterScaled,
  scaledObsGridCenter = scaledObjects$coordsDataScaled,
  plot.check = F)


## ---- warning = FALSE, message = FALSE----------------------------------------
modelCode <- nimbleCode({
  ##---- SPATIAL PROCESS 
  ## Prior for AC distribution parameter
  habCoeffSlope ~ dnorm(0, sd = 10)
  
  ## Intensity of the AC distribution point process
  habIntensity[1:numHabWindows] <- exp(habCoeffSlope * habCovs[1:numHabWindows])
  sumHabIntensity <- sum(habIntensity[1:numHabWindows])
  logHabIntensity[1:numHabWindows] <- log(habIntensity[1:numHabWindows])
  logSumHabIntensity <- log(sumHabIntensity)
  
  ## AC distribution
  for(i in 1:M){
    sxy[i, 1:2] ~ dbernppAC(
      lowerCoords = lowerHabCoords[1:numHabWindows, 1:2],
      upperCoords = upperHabCoords[1:numHabWindows, 1:2],
      logIntensities = logHabIntensity[1:numHabWindows],
      logSumIntensity = logSumHabIntensity,
      habitatGrid = habitatGrid[1:numGridRows,1:numGridCols],
      numGridRows =  numGridRows,
      numGridCols = numGridCols
    )
  }
  
  ##---- DEMOGRAPHIC PROCESS
  ## Prior for data augmentation
  psi ~ dunif(0,1)
  
  ## Data augmentation
  for (i in 1:M){
    z[i] ~ dbern(psi)
  }
  
  ##---- DETECTION PROCESS
  ## Priors for detection parameters
  sigma ~ dunif(0, 50)
  detCoeffInt ~ dnorm(0, sd = 10)
  detCoeffSlope ~ dnorm(0, sd = 10)
  
  ## Intensity of the detection point process
  detIntensity[1:numObsWindows] <- exp(detCoeffInt + detCoeffSlope * detCovs[1:numObsWindows]) 
  
  ## Detection process
  for (i in 1:M){
    y[i, 1:numMaxPoints, 1:3] ~ dpoisppDetection_normal(
      lowerCoords = obsLoCoords[1:numObsWindows, 1:2],
      upperCoords = obsUpCoords[1:numObsWindows, 1:2],
      s = sxy[i, 1:2],
      sd = sigma,
      baseIntensities = detIntensity[1:numObsWindows],
      numMaxPoints = numMaxPoints,
      numWindows = numObsWindows,
      indicator = z[i]
    )
  }
  
  ##---- DERIVED QUANTITIES
  ## Number of individuals in the population
  N <- sum(z[1:M])
})

## ---- warning = FALSE, message = FALSE----------------------------------------
sigma <- 1
psi <- 0.6
detCoeffInt <- 0.1
detCoeffSlope <- 0.5
habCoeffSlope <- -1.5

## ---- warning = FALSE, message = FALSE----------------------------------------
M <- 150

## ---- warning = FALSE, message = FALSE----------------------------------------
numMaxPoints <- 19 + 1

## ---- warning = FALSE, message = FALSE----------------------------------------
detCovs <- runif(dim(lowerAndUpperCoords$lowerObsCoords)[1],-1,1)
habCovs <- runif(dim(lowerAndUpperCoords$lowerHabCoords)[1],-1,1)

## ---- warning = FALSE, message = FALSE----------------------------------------
nimConstants <- list( M = M,
                      numObsWindows = dim(lowerAndUpperCoords$lowerObsCoords)[1],
                      numMaxPoints = numMaxPoints,
                      numHabWindows = dim(lowerAndUpperCoords$upperHabCoords)[1],
                      habitatGrid = lowerAndUpperCoords$habitatGrid,
                      numGridRows = dim(lowerAndUpperCoords$habitatGrid)[1],
                      numGridCols = dim(lowerAndUpperCoords$habitatGrid)[2])

nimData <- list( obsLoCoords = lowerAndUpperCoords$lowerObsCoords,
                 obsUpCoords = lowerAndUpperCoords$upperObsCoords,
                 lowerHabCoords = lowerAndUpperCoords$lowerHabCoords,
                 upperHabCoords = lowerAndUpperCoords$upperHabCoords,
                 detCovs = detCovs,
                 habCovs = habCovs)

## ---- warning = FALSE, message = FALSE----------------------------------------
nimInits <- list( psi = psi,
                  sigma = sigma,
                  detCoeffInt = detCoeffInt,
                  detCoeffSlope = detCoeffSlope,
                  habCoeffSlope = habCoeffSlope)

## ---- warning = FALSE, message = FALSE----------------------------------------
model <- nimbleModel( code = modelCode,
                      constants = nimConstants,
                      data = nimData,
                      inits = nimInits,
                      check = F,
                      calculate = F)

## ---- warning = FALSE, message = FALSE----------------------------------------
nodesToSim <- model$getDependencies(names(nimInits), self = F)
set.seed(1)
model$simulate(nodesToSim, includeData = FALSE)

## ---- warning = FALSE, message = FALSE----------------------------------------
N <- sum(model$z)

## ---- warning = FALSE, message = FALSE----------------------------------------
i = 7
## Number of detections for individual i
model$y[i,1,1]

## Plot of the habitat and trap grids 
plot( scaledObjects$coordsHabitatGridCenterScaled[,"y"] ~ scaledObjects$coordsHabitatGridCenterScaled[,"x"],
      pch = 1, cex = 0.5)
rect( xleft = lowerAndUpperCoords$lowerHabCoords[,1] ,
     ybottom = lowerAndUpperCoords$lowerHabCoords[,2] ,
     xright = lowerAndUpperCoords$upperHabCoords[,1],
     ytop = lowerAndUpperCoords$upperHabCoords[,2],
     col = adjustcolor("red", alpha.f = 0.4),
     border = "red")

rect( xleft = lowerAndUpperCoords$lowerObsCoords[,1] ,
     ybottom = lowerAndUpperCoords$lowerObsCoords[,2] ,
     xright = lowerAndUpperCoords$upperObsCoords[,1],
     ytop = lowerAndUpperCoords$upperObsCoords[,2],
     col = adjustcolor("blue",alpha.f = 0.4),
     border = "blue")

## Plot the activity center of individual i
points( model$sxy[i, 2] ~ model$sxy[i, 1],
        col = "orange", pch = 16)

## Plot detections of individual i
dets <- model$y[i,2:model$y[i,1,1], ]
points( dets[,2] ~ dets[,1],
        col = "green", pch = 16)

par(xpd = TRUE)
legend(x = -1, y = 13,
       legend = c("Habitat windows",
                  "Observation windows",
                  "Simulated AC",
                  "Detections"),
       pt.cex = c(1,1),
       horiz = T,
       pch = c(16, 16, 16, 16),
       col = c("red", "blue", "orange", "green"),
       bty = 'n')

## ---- warning = FALSE, message = FALSE----------------------------------------
nimData1 <- nimData
nimData1$y <- model$y
nimInits1 <- nimInits
nimInits1$z <- model$z
nimInits1$sxy <- model$sxy

## Create and compile the NIMBLE model
model <- nimbleModel( code = modelCode,
                      constants = nimConstants,
                      data = nimData1,
                      inits = nimInits1,
                      check = F,
                      calculate = F)

## Check the initial log-likelihood 
model$calculate()

## ----eval = FALSE-------------------------------------------------------------
#  cmodel <- compileNimble(model)

## ---- warning = FALSE, message = FALSE----------------------------------------
MCMCconf <- configureMCMC(model = model,
                          monitors  = c("N","sigma","psi","detCoeffInt",
                                        "detCoeffSlope","habCoeffSlope"),
                          control = list(reflective = TRUE),
                          thin = 10)
MCMC <- buildMCMC(MCMCconf)

## ----eval = FALSE-------------------------------------------------------------
#  cMCMC <- compileNimble(MCMC, project = model, resetFunctions = TRUE)
#  
#  ## Run MCMC
#  MCMCRuntime <- system.time(samples <- runMCMC( mcmc = cMCMC,
#                                                        nburnin = 500,
#                                                        niter = 10000,
#                                                        nchains = 3,
#                                                        samplesAsCodaMCMC = TRUE))

## ----eval = FALSE, echo = FALSE-----------------------------------------------
#  save(samples, MCMCRuntime, file = file.path(baseDir,"nimbleSCR/inst/extdata/pointProcess_samples.RData"))

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

## ---- warning = FALSE, message = FALSE----------------------------------------
## Print runtime
MCMCRuntime

## Traceplots and density plots for the tracked parameters
chainsPlot(samples, line = c(N, detCoeffInt, detCoeffSlope, habCoeffSlope, N/M, sigma))

## ---- warning = FALSE, message = FALSE----------------------------------------
modelCodeSemiCompleteLikelihood <- nimbleCode({
  #----- SPATIAL PROCESS
  ## Priors
  habCoeffInt ~ dnorm(0, sd = 10)
  habCoeffSlope ~ dnorm(0, sd = 10)

  ## Intensity of the AC distribution point process
  habIntensity[1:numHabWindows] <- exp(habCoeffInt + habCoeffSlope * habCovs[1:numHabWindows])
  sumHabIntensity <- sum(habIntensity[1:numHabWindows])
  logHabIntensity[1:numHabWindows] <- log(habIntensity[1:numHabWindows])
  logSumHabIntensity <- log(sum(habIntensity[1:numHabWindows] ))

  ## AC distribution
  for(i in 1:nDetected){
    sxy[i, 1:2] ~ dbernppAC(
      lowerCoords = lowerHabCoords[1:numHabWindows, 1:2],
      upperCoords = upperHabCoords[1:numHabWindows, 1:2],
      logIntensities = logHabIntensity[1:numHabWindows],
      logSumIntensity = logSumHabIntensity,
      habitatGrid = habitatGrid[1:numGridRows,1:numGridCols],
      numGridRows =  numGridRows,
      numGridCols = numGridCols
    )
  }

  ##---- DEMOGRAPHIC PROCESS
  ## Number of individuals in the population
  N ~ dpois(sumHabIntensity)
  ## Number of detected individuals
  nDetectedIndiv ~ dbin(probDetection, N)

  ##---- DETECTION PROCESS
  ## Probability that an individual in the population is detected at least once
  ## i.e. 1 - void probability over all detection windows
  probDetection <- 1 - marginalVoidProbNumIntegration(
    quadNodes = quadNodes[1:nNodes, 1:2, 1:numHabWindows],
    quadWeights = quadWeights[1:numHabWindows],
    numNodes = numNodes[1:numHabWindows],
    lowerCoords = obsLoCoords[1:numObsWindows, 1:2],
    upperCoords = obsUpCoords[1:numObsWindows, 1:2],
    sd = sigma,
    baseIntensities = detIntensity[1:numObsWindows],
    habIntensities = habIntensity[1:numHabWindows],
    sumHabIntensity = sumHabIntensity,
    numObsWindows = numObsWindows,
    numHabWindows = numHabWindows
  )

  ## Priors for detection parameters
  sigma ~ dunif(0, 50)
  detCoeffInt ~ dnorm(0, sd = 10)
  detCoeffSlope ~ dnorm(0, sd = 10)

  ## Intensity of the detection point process
  detIntensity[1:numObsWindows] <- exp(detCoeffInt + detCoeffSlope * detCovs[1:numObsWindows])
  ## Detection process
  ## Note that this conditions on the fact that individuals are detected (at least once)
  ## So, at the bottom of this model code we deduct log(probDetection) from the log-likelihood
  ## function for each individual
  for (i in 1:nDetected){
    y[i, 1:numMaxPoints, 1:3] ~ dpoisppDetection_normal(
      lowerCoords = obsLoCoords[1:numObsWindows, 1:2],
      upperCoords = obsUpCoords[1:numObsWindows, 1:2],
      s = sxy[i, 1:2],
      sd = sigma,
      baseIntensities = detIntensity[1:numObsWindows],
      numMaxPoints = numMaxPoints,
      numWindows = numObsWindows,
      indicator = 1
    )
  }
  ## Normalization: normData can be any scalar in the data provided when building the model
  ## The dnormalizer is a custom distribution defined for efficiency, where the input data
  ## does not matter. It makes it possible to use the general dpoippDetection_normal function
  ## when either data augmentation or the SCDL is employed
  logDetProb <- log(probDetection)
  normData ~ dnormalizer(logNormConstant = -nDetected * logDetProb)
})

## ---- warning = FALSE, message = FALSE----------------------------------------
idDetected <- which(nimData1$y[,1,1] > 0)
## Subset data to detected individuals only
nimData1$y <- nimData1$y[idDetected,,]
## Provide the number of detected individuals as constant
nimConstants$nDetected <- length(idDetected)
## With this model, we also need to provide the number of detected individuals as data for the estimation of population size.
nimData1$nDetectedIndiv <- length(idDetected)
## As mentioned above, "normData" can take any value.
nimData1$normData <- 1

## We also provide initial values for the new parameters that need to be estimated
nimInits1$N <- 100
nimInits1$habCoeffInt <- 0.5
nimInits1$sxy <- nimInits1$sxy[idDetected,]

## ---- warning = FALSE, message = FALSE----------------------------------------
## Number of equal subintervals for each dimension of a grid cell
nPtsPerDim <- 2
## Number of points to use for the numerical integration for each grid cell
nNodes <- nPtsPerDim^2
## Generate midpoint nodes coordinates for numerical integration using the "getMidPointNodes" function
nodesRes <- getMidPointNodes( nimData1$lowerHabCoords,
                           nimData1$upperHabCoords,
                           nPtsPerDim)
## Add this info to the data and constant objects
nimData1$quadNodes <- nodesRes$quadNodes
nimData1$quadWeights <- nodesRes$quadWeights
nimData1$numNodes <- rep(nNodes,dim(nimData1$lowerHabCoords)[1])
nimConstants$nNodes <- dim(nodesRes$quadNodes)[1]

## ---- eval = TRUE, warning = FALSE, message = FALSE---------------------------
model <- nimbleModel(code = modelCodeSemiCompleteLikelihood,
                     constants = nimConstants,
                     data = nimData1,
                     inits = nimInits1,
                     check = F,
                     calculate = F)

model$calculate()

## ----eval = FALSE-------------------------------------------------------------
#  cmodel <- compileNimble(model)
#  cmodel$calculate()

## -----------------------------------------------------------------------------
MCMCconf <- configureMCMC(model = model,
                          monitors =  c("N","sigma","probDetection","habCoeffInt", "detCoeffInt","detCoeffSlope","habCoeffSlope"),
                          control = list(reflective = TRUE),
                          thin = 10)

MCMC <- buildMCMC(MCMCconf)

## ---- eval = FALSE, warning = FALSE, message = FALSE--------------------------
#  cMCMC <- compileNimble(MCMC, project = model, resetFunctions = TRUE)
#  
#  ## Run MCMC
#  MCMCRuntime1 <- system.time(samples1 <- runMCMC( mcmc = cMCMC,
#                                                       nburnin = 500,
#                                                       niter = 10000,
#                                                       nchains = 3,
#                                                       samplesAsCodaMCMC = TRUE))

## ----eval = FALSE, echo = FALSE-----------------------------------------------
#  save(samples1, MCMCRuntime1, file = file.path(baseDir,"nimbleSCR/inst/extdata/pointProcess_samples1.RData"))

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

## ---- warning = FALSE, message = FALSE----------------------------------------
MCMCRuntime1

## Plot check
chainsPlot(samples1, line = c(N, detCoeffInt, detCoeffSlope, NA, habCoeffSlope, NA, sigma))