run_classic: Run spatial capture-recapture models in nimble using parallel...

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

Run 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_classic(
  model,
  data,
  constants,
  inits,
  params,
  dimensions = NULL,
  nimFun = NULL,
  niter = 1000,
  nburnin = 100,
  thin = 1,
  nchains = 1,
  parallel = FALSE,
  RNGseed = NULL,
  s_alias = "s"
)

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'.
s_alias	A character value used to identify the latent activity center coordinates used in the model. Default is "s". Note that the length of s_alias must be either 1 (e.g., "s") or 2 (e.g., c("s","su")).

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

# add some random noise to locations
traps <- traps + runif(prod(dim(traps)),-20,20) 

mysigma = 300 # simulate sigma of 300 m
mycrs = 32608 # EPSG for WGS 84 / UTM zone 8N
pixelWidth = 100 # store pixelWidth

Grid = grid_classic(X = traps, crs_ = mycrs, buff = 3*mysigma, 
res = pixelWidth) 

# create polygon to use as a 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 = 1, N = 200, K = 4, base_encounter = 0.15, 
enc_dist = "binomial",hab_mask = hab_mask, setSeed = 50)

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

# rescale inputs
rescale_list = rescale_classic(X = traps, ext = Grid$ext, s.st = s.st3d, 
hab_mask = hab_mask)

# store rescaled extent
ext = rescale_list$ext

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

# add rescaled traps
data$X = rescale_list$X

# prepare constants
constants = list(M = 500,n0 = nrow(data$y),J=dim(data$y)[2], 
                 K=dim(data3d$y)[3],x_lower = ext[1], x_upper = ext[2], 
                 y_lower = ext[3], y_upper = ext[4],sigma_upper = 1000, 
                 A = sum(hab_mask)*(pixelWidth)^2,pixelWidth=pixelWidth)

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

# add hab_mask and OK for habitat check
data$hab_mask = hab_mask
data$OK = rep(1,constants$M)

# get initial activity center starting values
s.st3d = rescale_list$s.st

# define all initial values
inits = list(sigma = runif(1, 250, 350), s = s.st3d,psi=runif(1,0.2,0.3),
          p0 = runif(1, 0.05, 0.15),   pOK = data$OK, 
          z=c(rep(NA,constants$n0),rep(0,constants$M-constants$n0)))

# parameters to monitor
params = c("sigma","psi","p0","N","D")

# get model
scr_model = get_classic(dim_y = 2, enc_dist = "binomial",sex_sigma = FALSE,
hab_mask=TRUE)

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

# summarize output
samples = do.call(rbind, out)
par(mfrow=c(1,1))
hist(samples[,which(dimnames(out[[1]])[[2]]=="N")], xlab = "Abundance", 
xlim = c(0,500), main="")
abline(v=200, col="red") # add line for simulated abundance

# not run
#nimSummary(out)

## End(Not run)

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