run_classic.RdA wrapper function to conduct Markov Chain Monte Carlo (MCMC) sampling using nimble.
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"
)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'.
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")).
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))
# 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)
}