Extract spatially-varying coefficient MCMC samples

Function for extracting the full spatially-varying coefficient MCMC samples from an spOccupancy model object.

Usage

getSVCSamples(object, pred.object, ...)

Arguments

object

an object of class svcPGOcc, svcPGBinom, svcTPGOcc, svcTPGBinom.

pred.object

a prediction object from a spatially-varying coefficient model fit using spOccupancy. Should be of class predict.svcPGOcc, predict.svcPGBinom, predict.svcTPGOcc, or predict.svcTPGBinom. If specified, SVC samples are extracted at the prediction locations.

...

currently no additional arguments

Author

Jeffrey W. Doser doserjef@msu.edu,

Value

A list of coda::mcmc objects of the spatially-varying coefficient MCMC samples for all spatially-varying coefficients estimated in the model (including the intercept if specified). Note these values correspond to the sum of the estimated spatial and non-spatial effect to give the overall effect of the covariate at each location. Each element of the list is a two-dimensional matrix where dimensions correspond to MCMC sample and site. If pred.object is specified, values are returned for the prediction locations instead of the sampled locations.

Examples

set.seed(400)
# Simulate Data -----------------------------------------------------------
J.x <- 8
J.y <- 8
J <- J.x * J.y
n.rep <- sample(2:4, J, replace = TRUE)
beta <- c(0.5, 2)
p.occ <- length(beta)
alpha <- c(0, 1)
p.det <- length(alpha)
phi <- c(3 / .6, 3 / .8)
sigma.sq <- c(1.2, 0.7)
svc.cols <- c(1, 2)
dat <- simOcc(J.x = J.x, J.y = J.y, n.rep = n.rep, beta = beta, alpha = alpha, 
              sigma.sq = sigma.sq, phi = phi, sp = TRUE, cov.model = 'exponential', 
              svc.cols = svc.cols)
# Detection-nondetection data
y <- dat$y
# Occupancy covariates
X <- dat$X
# Detection covarites
X.p <- dat$X.p
# Spatial coordinates
coords <- dat$coords

# Package all data into a list
occ.covs <- X[, -1, drop = FALSE]
colnames(occ.covs) <- c('occ.cov')
det.covs <- list(det.cov.1 = X.p[, , 2])
data.list <- list(y = y, 
                  occ.covs = occ.covs, 
                  det.covs = det.covs, 
                  coords = coords)

# Number of batches
n.batch <- 10
# Batch length
batch.length <- 25
n.iter <- n.batch * batch.length
# Priors 
prior.list <- list(beta.normal = list(mean = 0, var = 2.72), 
                   alpha.normal = list(mean = 0, var = 2.72),
                   sigma.sq.ig = list(a = 2, b = 1), 
                   phi.unif = list(a = 3/1, b = 3/.1)) 
# Initial values
inits.list <- list(alpha = 0, beta = 0,
                   phi = 3 / .5, 
                   sigma.sq = 2,
                   w = matrix(0, nrow = length(svc.cols), ncol = nrow(X)),
                   z = apply(y, 1, max, na.rm = TRUE))
# Tuning
tuning.list <- list(phi = 1) 

out <- svcPGOcc(occ.formula = ~ occ.cov, 
                det.formula = ~ det.cov.1, 
                data = data.list, 
                inits = inits.list, 
                n.batch = n.batch, 
                batch.length = batch.length, 
                accept.rate = 0.43, 
                priors = prior.list,
                cov.model = 'exponential', 
                svc.cols = c(1, 2),
                tuning = tuning.list, 
                n.omp.threads = 1, 
                verbose = TRUE, 
                NNGP = TRUE, 
                n.neighbors = 5, 
                search.type = 'cb', 
                n.report = 10, 
                n.burn = 50, 
                n.thin = 1)
#> ----------------------------------------
#> 	Preparing to run the model
#> ----------------------------------------
#> ----------------------------------------
#> 	Building the neighbor list
#> ----------------------------------------
#> ----------------------------------------
#> Building the neighbors of neighbors list
#> ----------------------------------------
#> ----------------------------------------
#> 	Model description
#> ----------------------------------------
#> NNGP Occupancy model with Polya-Gamma latent
#> variable fit with 64 sites.
#> 
#> Samples per chain: 250 (10 batches of length 25)
#> Burn-in: 50 
#> Thinning Rate: 1 
#> Number of Chains: 1 
#> Total Posterior Samples: 200 
#> 
#> Number of spatially-varying coefficients: 2 
#> Using the exponential spatial correlation model.
#> 
#> Using 5 nearest neighbors.
#> 
#> Source compiled with OpenMP support and model fit using 1 thread(s).
#> 
#> Adaptive Metropolis with target acceptance rate: 43.0
#> ----------------------------------------
#> 	Chain 1
#> ----------------------------------------
#> Sampling ... 
#> Batch: 10 of 10, 100.00%

svc.samples <- getSVCSamples(out)
str(svc.samples)
#> List of 2
#>  $ (Intercept): 'mcmc' num [1:200, 1:64] 0.612 -1.485 1.695 2.008 0.631 ...
#>   ..- attr(*, "mcpar")= num [1:3] 1 200 1
#>  $ occ.cov    : 'mcmc' num [1:200, 1:64] 2.95 2 1.47 2.27 3.89 ...
#>   ..- attr(*, "mcpar")= num [1:3] 1 200 1