Generate realized density surface from MCMC output — realized

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

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".

Value

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

Details

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

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