realized_density: Generate realized density surface from MCMC output

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

Generate realized density surface from MCMC output

Description

Streamlined construction of realized density surface from posterior samples of latent indicator variable (z) and activity center coordinates (s)

Usage

realized_density(
  samples,
  grid,
  ext,
  crs_,
  site,
  hab_mask = FALSE,
  discrete = FALSE,
  s_alias = "s",
  z_alias = "z"
)

Arguments

samples	Either a matrix (for a single MCMC chain) or list of posterior samples from multiple chains from MCMC sampling; possibly returned from run_classic.
grid	A matrix or array object of the the state-space grid. This is returned from grid_classic.
ext	An Extent object from the raster package. This is returned from grid_classic. Note that ext is only used when when the state-space grid is a 3-dimensional array.
crs_	The UTM coordinate reference system (EPSG code) used for your location provided as an integer (e.g., 32608 for WGS 84/UTM Zone 8N).
site	Site identifier variable included for detected and augmented individuals used as a constant in model runs.
hab_mask	Either FALSE (the default) or a matrix or arary output from mask_polygon or mask_raster functions.
discrete	Logical. Either FALSE (the default) to indicate that the MCMC output is from a 'classic' or continuous SCR model or TRUE that the MCMC output is from a 'discrete' SCR model.
s_alias	A character value used to identify the latent activity center coordinates used in the model. Default is "s".
z_alias	A character value used to identify the latent inclusion indicator used in the model. Default is "z".

Details

This function automates the construction of realized density surfaces from MCMC samples.

Value

a raster object or a list of raster objects if state-space grid is an array.

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

# 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 # rescale to kilometers
ext = as.vector(Grid$ext) # 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 = 1000, 
            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, 250, 350), s = s.st3d,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","s","z")

# get model
scr_model = get_classic(dim_y = 2, enc_dist = "binomial",sex_sigma = FALSE,
trapsClustered = 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 (exclude lengthy parameters "s" and "z")
nimSummary(out, exclude = c("s","z"), trace = TRUE)

library(tictoc)       
tic() # track time
r = realized_density(samples = out, grid = Grid$grid, 
ext = Grid$ext, crs_ = mycrs, site = NULL, hab_mask = FALSE, 
discrete = FALSE)       
toc()      

# load virdiis color pallete library      
library(viridis)
library(raster)

# make simple raster plot
label = expression("Realized density (activity centers/100 m"^2*")")
plot(r, col=viridis(100),main=label)

## End(Not run)

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