run_discrete: Run discrete spatial capture-recapture models in 'nimble'...

In sitkensis22/localSCR: Supporting functions for advanced Bayesian spatial capture-recapture modeling using 'nimble'

Run discrete spatial capture-recapture models in 'nimble' using parallel processing

Description

A wrapper function to conduct Markov Chain Monte Carlo (MCMC) sampling using 'nimble'.

Usage

run_discrete(
  model,
  data,
  constants,
  inits,
  params,
  dimensions = NULL,
  nimFun = NULL,
  niter = 1000,
  nburnin = 100,
  thin = 1,
  nchains = 1,
  parallel = FALSE,
  RNGseed = NULL
)

Arguments

model	The nimbleCode used to define model in nimble package, possibly generated from get_classic.
data	A list of data inputs needed to run SCR models in nimble.
constants	A list of constants needed to run SCR models in nimble.
inits	Starting values for stochastic parameters to begin MCMC sampling.
params	A vector of character strings that define the parameters to trace in the MCMC simulation.
dimensions	A named list of dimensions for variables. Only needed for variables used with empty indices in model code that are not provided in constants or data.
nimFun	A list of nimbleFunction(s) to be used in the model when using parallel = TRUE.
niter	An integer value of the total number of MCMC iterations to run per chain.
nburnin	An integer value of the number of MCMC iterations to discard as 'burnin'.
thin	An integer value of the amount of thinning of the chain. For example, thin=2 would retain every other MCMC sample.
nchains	An integer value for the number of MCMC chains to run
parallel	A logical value indicating whether MCMC chains shoud be run in parallel processing. Default is FALSE.
RNGseed	An integer value specifying the random number generating seed using in parallel processing. Also used when parallel = FALSE in setSeed within runMCMC function in 'nimble'. This ensures that the MCMC samples will be the same during each run using the same data, etc. Default is NULL.

Details

This function provides a wrapper to easily run Bayesian SCR models using nimble.

Value

A list of MCMC samples for each parameter traced with length equal to the number of chains run.

Author(s)

Daniel Eacker

Examples

## Not run: 
# simulate a single trap array with random positional noise
x <- seq(-800, 800, length.out = 5)
y <- seq(-800, 800, length.out = 5)
traps <- as.matrix(expand.grid(x = x, y = y))
set.seed(200)
traps <- traps + runif(prod(dim(traps)),-20,20) 

mysigma = c(220, 300) # simulate sex-specific scaling parameter
mycrs = 32608 # EPSG for WGS 84 / UTM zone 8N
pixelWidth = 100 # store pixelWidth or grid resolution

# create state-space
Grid = grid_classic(X = traps, crs_ = mycrs, buff = 3*mysigma, res = pixelWidth)

# create polygon for mask
library(sf)
poly = st_sfc(st_polygon(x=list(matrix(c(-1765,-1765,1730,-1650,1600,1650,0,1350,
           -800,1700,-1850,1000,-1765,-1765),ncol=2, byrow=TRUE))), crs =  mycrs)

# create habitat mask
hab_mask = mask_polygon(poly = poly, grid = Grid$grid, crs_ = mycrs, 
prev_mask = NULL)

# simulate data for uniform state-space and habitat mask
data3d = sim_classic(X = traps, ext = Grid$ext, crs_ = mycrs, sigma_ = mysigma, 
             prop_sex = 0.6, N = 200, K = 4, base_encounter = 0.15, enc_dist = "binomial",
             hab_mask = hab_mask, setSeed = 100)

# total augmented population size 
M = 400

# get initial activity center starting values
s.st = initialize_classic(y=data3d$y, M=M, X=traps, ext=Grid$ext, 
                           hab_mask = hab_mask)

# convert traps and starting locations to discrete format
d_list = discretize_classic(X=traps, grid = Grid$grid, s.st = s.st, 
         crs_ = mycrs, hab_mask = hab_mask)

# inspect discrete data list
str(d_list)

# prepare data
data = list(y=data3d$y)
data$y = data$y[which(apply(data$y, 1, sum)!=0),,] # remove augmented records 
# covert to 2d by summing over individuals and traps
data$y = apply(data$y, c(1,2), sum) 

# add discretized traps
data$X = d_list$X100 # convert units to km

# add grid
data$grid = d_list$grid/1000 # convert units to km

# prepare constants (note get density in activity center/100 m2 rather 
# than activity centers/m2)
constants = list(M = M,n0 = nrow(data$y),J=dim(data$y)[2], K=dim(data3d$y)[3],
nPix=d_list$nPix,pixArea = 0.1^2,
sigma_upper = 1, A = (sum(hab_mask)*(pixelWidth/100)^2))

# augment sex
data$sex = c(data3d$sex,rep(NA,constants$M-length(data3d$sex)))

# add z and zeros vector data for latent inclusion indicator
data$z = c(rep(1,constants$n0),rep(NA,constants$M - constants$n0))
data$zeros =  c(rep(NA,constants$n0),rep(0,constants$M - constants$n0))

# define all initial values 
inits = list(sigma = runif(2, 0.250, 0.350), s = d_list$s.st,alpha0 = 3,
  p0 = runif(1, 0.05, 0.15),psi_sex = runif(1,0.5, 0.7),
        sex=c(rep(NA,constants$n0),rbinom(constants$M-constants$n0,1,0.6)),
        z=c(rep(NA,constants$n0),rep(0,constants$M-constants$n0)))

# parameters to monitor
params = c("sigma","psi","p0","N","D","EN","psi_sex","alpha0","s","z")

   
# get model
discrete_model = get_discrete(type="marked",dim_y = 2, enc_dist = "binomial",
                 sex_sigma = TRUE,trapsClustered=FALSE)

# run model
library(tictoc)
tic() # track time elapsed
out = run_discrete(model = discrete_model, data=data, constants=constants,
       inits=inits, params = params,niter = 10000, nburnin=1000, 
       thin=1, nchains=2, parallel=TRUE,  RNGseed = 500)
toc()

nimSummary(out, exclude = c("s","z"))

## End(Not run)

sitkensis22/localSCR documentation built on May 15, 2022, 5:26 p.m.