Summarize MCMC chain output from nimble

Summarizes lists of MCMC chain output from nimble

nimSummary(d, trace = FALSE, plot_all = FALSE, exclude = NULL, digits = 3)

Arguments

d

A list of MCMC samples from each chain returned from run_classic.

trace

A logical value indicating whether or not traces of MCMC samples should be plotted. Default is FALSE.

plot_all

A logical value indicating whether or not all parameters should be included in plots. This assumes that some parameters are excluded in the summary table (i.e., exclude != NULL). Default is FALSE.

exclude

Either NULL or a scalar or vector containing parameter(s) to exclude from summary. Note that high dimensional parameters (e.g., s[1, 1, 1]) can be excluded by calling exclude = "s". Default is NULL.

digits

An integer value indicating how many digits the output should be rounded to.

Value

A dataframe of summary statistics for MCMC samples.

Details

This function summarizes Bayesian SCR models run using nimble including mean and quantiles of samples, as well as effective sample size and Rhat statistics. Note that f0 is the proportion of samples that are greater than zero. Also, at least 2 chains must be run to use this function.

Author

Daniel Eacker

Examples

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

# create grid and extent
Grid = grid_classic(X = traps, crs_ = mycrs, buff = 3*mysigma, 
res = pixelWidth)

# simulate encounter data
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 = FALSE, setSeed = 200)

# prepare data
data = list(y=data3d$y)
data$y = data$y[which(apply(data$y, 1, sum)!=0),,] # remove augmented records
data$y = apply(data$y, c(1,2), sum) # covert to 2d by summing over occasions
data$X = traps/1000 # rescale to kilometers
ext = as.vector(Grid$ext)/1000 # recale to kilometers

# 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 = 1, 
           A = prod(ext[2]-ext[1],ext[4]-ext[3]))

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

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

# define all initial values
inits = list(sigma = runif(1, 0.250, 0.350), s = s.st3d/1000,
             psi=runif(1,0.2,0.3), p0 = runif(1, 0.05, 0.15),
             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)

# run model
library(tictoc)
tic() # track time elapsed
out = run_classic(model = scr_model, data=data, constants=constants, 
                  inits=inits, params = params,niter = 5000, nburnin=1000,
                  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

# summarize output
nimSummary(out, trace=TRUE, plot_all=TRUE)
}