Function for prediction at new locations for single-species integrated spatial occupancy models
predict.spIntPGOcc.RdThe function predict collects posterior predictive samples for a set of new locations given an object of class `spIntPGOcc`.
Usage
# S3 method for spIntPGOcc
predict(object, X.0, coords.0, n.omp.threads = 1, verbose = TRUE,
n.report = 100, ignore.RE = FALSE, type = 'occupancy', ...)Arguments
- object
an object of class
spIntPGOcc.- X.0
the design matrix for prediction locations. This should include a column of 1s for the intercept. Covariates should have the same column names as those used when fitting the model with
spIntPGOcc.- coords.0
the spatial coordinates corresponding to
X.0. Note thatspOccupancyassumes coordinates are specified in a projected coordinate system.- n.omp.threads
a positive integer indicating the number of threads to use for SMP parallel processing. The package must be compiled for OpenMP support. For most Intel-based machines, we recommend setting
n.omp.threadsup to the number of hyperthreaded cores. Note,n.omp.threads> 1 might not work on some systems.- verbose
if
TRUE, model specification and progress of the sampler is printed to the screen. Otherwise, nothing is printed to the screen.- n.report
the interval to report sampling progress.
- ignore.RE
logical value that specifies whether or not to remove random occurrence (or detection if
type = 'detection') effects from the subsequent predictions. IfTRUE, random effects will be included. IfFALSE, random effects will be set to 0 and predictions will only be generated from the fixed effects.- type
a quoted keyword indicating what type of prediction to produce. Valid keywords are 'occupancy' to predict latent occupancy probability and latent occupancy values (this is the default), or 'detection' to predict detection probability given new values of detection covariates. Note that prediction of detection probability is not currently supported for integrated models.
- ...
currently no additional arguments
Author
Jeffrey W. Doser doserjef@msu.edu,
Andrew O. Finley finleya@msu.edu
Value
An object of class predict.spIntPGOcc that is a list comprised of:
- psi.0.samples
a
codaobject of posterior predictive samples for the latent occurrence probability values.- z.0.samples
a
codaobject of posterior predictive samples for the latent occurrence values.
The return object will include additional objects used for standard extractor functions.
Examples
set.seed(400)
# Simulate Data -----------------------------------------------------------
# Number of locations in each direction. This is the total region of interest
# where some sites may or may not have a data source.
J.x <- 8
J.y <- 8
J.all <- J.x * J.y
# Number of data sources.
n.data <- 4
# Sites for each data source.
J.obs <- sample(ceiling(0.2 * J.all):ceiling(0.5 * J.all), n.data, replace = TRUE)
# Replicates for each data source.
n.rep <- list()
for (i in 1:n.data) {
n.rep[[i]] <- sample(1:4, size = J.obs[i], replace = TRUE)
}
# Occupancy covariates
beta <- c(0.5, 0.5)
p.occ <- length(beta)
# Detection covariates
alpha <- list()
alpha[[1]] <- runif(2, 0, 1)
alpha[[2]] <- runif(3, 0, 1)
alpha[[3]] <- runif(2, -1, 1)
alpha[[4]] <- runif(4, -1, 1)
p.det.long <- sapply(alpha, length)
p.det <- sum(p.det.long)
sigma.sq <- 2
phi <- 3 / .5
sp <- TRUE
# Simulate occupancy data.
dat <- simIntOcc(n.data = n.data, J.x = J.x, J.y = J.y, J.obs = J.obs,
n.rep = n.rep, beta = beta, alpha = alpha, sp = sp,
phi = phi, sigma.sq = sigma.sq, cov.model = 'spherical')
y <- dat$y
X <- dat$X.obs
X.p <- dat$X.p
sites <- dat$sites
X.0 <- dat$X.pred
psi.0 <- dat$psi.pred
coords <- as.matrix(dat$coords.obs)
coords.0 <- as.matrix(dat$coords.pred)
# Package all data into a list
occ.covs <- X[, 2, drop = FALSE]
colnames(occ.covs) <- c('occ.cov')
det.covs <- list()
# Add covariates one by one
det.covs[[1]] <- list(det.cov.1.1 = X.p[[1]][, , 2])
det.covs[[2]] <- list(det.cov.2.1 = X.p[[2]][, , 2],
det.cov.2.2 = X.p[[2]][, , 3])
det.covs[[3]] <- list(det.cov.3.1 = X.p[[3]][, , 2])
det.covs[[4]] <- list(det.cov.4.1 = X.p[[4]][, , 2],
det.cov.4.2 = X.p[[4]][, , 3],
det.cov.4.3 = X.p[[4]][, , 4])
data.list <- list(y = y,
occ.covs = occ.covs,
det.covs = det.covs,
sites = sites,
coords = coords)
J <- length(dat$z.obs)
# Initial values
inits.list <- list(alpha = list(0, 0, 0, 0),
beta = 0,
phi = 3 / .5,
sigma.sq = 2,
w = rep(0, J),
z = rep(1, J))
# Priors
prior.list <- list(beta.normal = list(mean = 0, var = 2.72),
alpha.normal = list(mean = list(0, 0, 0, 0),
var = list(2.72, 2.72, 2.72, 2.72)),
phi.unif = c(3/1, 3/.1),
sigma.sq.ig = c(2, 2))
# Tuning
tuning.list <- list(phi = 1)
# Number of batches
n.batch <- 40
# Batch length
batch.length <- 25
# Note that this is just a test case and more iterations/chains may need to
# be run to ensure convergence.
out <- spIntPGOcc(occ.formula = ~ occ.cov,
det.formula = list(f.1 = ~ det.cov.1.1,
f.2 = ~ det.cov.2.1 + det.cov.2.2,
f.3 = ~ det.cov.3.1,
f.4 = ~ det.cov.4.1 + det.cov.4.2 + det.cov.4.3),
data = data.list,
inits = inits.list,
n.batch = n.batch,
batch.length = batch.length,
accept.rate = 0.43,
priors = prior.list,
cov.model = "spherical",
tuning = tuning.list,
n.omp.threads = 1,
verbose = TRUE,
NNGP = TRUE,
n.neighbors = 5,
search.type = 'cb',
n.report = 10,
n.burn = 500,
n.thin = 1)
#> ----------------------------------------
#> Preparing to run the model
#> ----------------------------------------
#> ----------------------------------------
#> Building the neighbor list
#> ----------------------------------------
#> ----------------------------------------
#> Building the neighbors of neighbors list
#> ----------------------------------------
#> ----------------------------------------
#> Model description
#> ----------------------------------------
#> NNGP Spatial Integrated Occupancy Model with Polya-Gamma latent
#> variable fit with 54 sites.
#>
#> Integrating 4 occupancy data sets.
#>
#> Samples per chain: 1000 (40 batches of length 25)
#> Burn-in: 500
#> Thinning Rate: 1
#> Number of Chains: 1
#> Total Posterior Samples: 500
#>
#> Using the spherical 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 40, 25.00%
#> Parameter Acceptance Tuning
#> phi 80.0 1.11628
#> -------------------------------------------------
#> Batch: 20 of 40, 50.00%
#> Parameter Acceptance Tuning
#> phi 80.0 1.23368
#> -------------------------------------------------
#> Batch: 30 of 40, 75.00%
#> Parameter Acceptance Tuning
#> phi 64.0 1.36343
#> -------------------------------------------------
#> Batch: 40 of 40, 100.00%
summary(out)
#>
#> Call:
#> spIntPGOcc(occ.formula = ~occ.cov, det.formula = list(f.1 = ~det.cov.1.1,
#> f.2 = ~det.cov.2.1 + det.cov.2.2, f.3 = ~det.cov.3.1, f.4 = ~det.cov.4.1 +
#> det.cov.4.2 + det.cov.4.3), data = data.list, inits = inits.list,
#> priors = prior.list, tuning = tuning.list, cov.model = "spherical",
#> NNGP = TRUE, n.neighbors = 5, search.type = "cb", n.batch = n.batch,
#> batch.length = batch.length, accept.rate = 0.43, n.omp.threads = 1,
#> verbose = TRUE, n.report = 10, n.burn = 500, n.thin = 1)
#>
#> Samples per Chain: 1000
#> Burn-in: 500
#> Thinning Rate: 1
#> Number of Chains: 1
#> Total Posterior Samples: 500
#> Run Time (min): 0.0034
#>
#> ----------------------------------------
#> Occurrence
#> ----------------------------------------
#> Fixed Effects (logit scale):
#> Mean SD 2.5% 50% 97.5% Rhat ESS
#> (Intercept) 0.2798 0.4286 -0.581 0.2637 1.1541 NA 83
#> occ.cov 0.6526 0.4421 -0.173 0.6305 1.5651 NA 161
#>
#> ----------------------------------------
#> Data source 1 Detection
#> ----------------------------------------
#> Fixed Effects (logit scale):
#> Mean SD 2.5% 50% 97.5% Rhat ESS
#> (Intercept) 0.7223 0.6498 -0.5366 0.7176 1.9761 NA 358
#> det.cov.1.1 0.7057 0.6615 -0.5103 0.6900 1.9990 NA 311
#>
#> ----------------------------------------
#> Data source 2 Detection
#> ----------------------------------------
#> Fixed Effects (logit scale):
#> Mean SD 2.5% 50% 97.5% Rhat ESS
#> (Intercept) 0.6186 0.4118 -0.1791 0.6143 1.4116 NA 248
#> det.cov.2.1 0.4807 0.3914 -0.2465 0.4852 1.2507 NA 345
#> det.cov.2.2 0.3613 0.5043 -0.5803 0.3474 1.3753 NA 354
#>
#> ----------------------------------------
#> Data source 3 Detection
#> ----------------------------------------
#> Fixed Effects (logit scale):
#> Mean SD 2.5% 50% 97.5% Rhat ESS
#> (Intercept) -0.1296 0.3552 -0.8130 -0.1452 0.5566 NA 317
#> det.cov.3.1 1.0606 0.4334 0.2623 1.0484 1.9994 NA 273
#>
#> ----------------------------------------
#> Data source 4 Detection
#> ----------------------------------------
#> Fixed Effects (logit scale):
#> Mean SD 2.5% 50% 97.5% Rhat ESS
#> (Intercept) -0.6703 0.5981 -1.8455 -0.6418 0.5291 NA 266
#> det.cov.4.1 1.0340 0.7884 -0.3737 0.9873 2.6922 NA 183
#> det.cov.4.2 -0.2900 0.6176 -1.6070 -0.2510 0.8572 NA 327
#> det.cov.4.3 -1.7656 0.7798 -3.3158 -1.7166 -0.4113 NA 146
#>
#> ----------------------------------------
#> Spatial Covariance
#> ----------------------------------------
#> Mean SD 2.5% 50% 97.5% Rhat ESS
#> sigma.sq 2.2455 2.1599 0.3612 1.5353 7.4709 NA 13
#> phi 18.0089 7.3563 5.0786 17.7260 29.5405 NA 47
# Predict at new locations ------------------------------------------------
out.pred <- predict(out, X.0, coords.0, verbose = FALSE)