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, ...)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.
- ...
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
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 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.0055
#>
#> Occurrence (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 (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 (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 (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 (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)