run_discrete.RdA wrapper function to conduct Markov Chain Monte Carlo (MCMC) sampling using 'nimble'.
run_discrete(
model,
data,
constants,
inits,
params,
dimensions = NULL,
nimFun = NULL,
niter = 1000,
nburnin = 100,
thin = 1,
nchains = 1,
parallel = FALSE,
RNGseed = NULL
)The nimbleCode used to define model in nimble package,
possibly generated from get_classic.
A list of data inputs needed to run SCR models in nimble.
A list of constants needed to run SCR models in nimble.
Starting values for stochastic parameters to begin MCMC sampling.
A vector of character strings that define the parameters to trace in the MCMC simulation.
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.
A list of nimbleFunction(s) to be used in the model when using
parallel = TRUE.
An integer value of the total number of MCMC iterations to run per chain.
An integer value of the number of MCMC iterations to discard as 'burnin'.
An integer value of the amount of thinning of the chain. For
example, thin=2 would retain every other MCMC sample.
An integer value for the number of MCMC chains to run
A logical value indicating whether MCMC chains shoud be run
in parallel processing. Default is FALSE.
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.
A list of MCMC samples for each parameter traced with length equal to the number of chains run.
This function provides a wrapper to easily run Bayesian SCR models
using nimble.
if (FALSE) {
# 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"))
}