inst/doc/WolverinePointProcess.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 = 5, fig.height = 5) 

## ---- echo = FALSE------------------------------------------------------------
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'] == 'weizhang'){               
    baseDir <- '/Users/weizhang/Documents/GitHub/nimbleSCR/' ## Wei 
} else baseDir <- NULL

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

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

## -----------------------------------------------------------------------------
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)

})

## -----------------------------------------------------------------------------
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
})

## -----------------------------------------------------------------------------
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
})

## -----------------------------------------------------------------------------
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
})

## ---- eval = FALSE------------------------------------------------------------
#  load("WolverinePointProcess_Data.RData")

## ----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)

## ----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)

## ----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)

## ----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"))

## ----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"))
##}

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

## -----------------------------------------------------------------------------
round(MCMC_samples1_summary, 2)
round(MCMC_samples2_summary, 2)
round(MCMC_samples3_summary, 2)
round(MCMC_samples4_summary, 2)

## ----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"))

## -----------------------------------------------------------------------------
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")] 

## -----------------------------------------------------------------------------
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)