Spatially-varying coefficient models in spOccupancy
Jeffrey W. Doser
2022
Source:vignettes/svcUnivariateHTML.Rmd
svcUnivariateHTML.RmdIntroduction
This vignette details spOccupancy functionality introduced in v0.5.0 to fit single-species spatially-varying coefficient (SVC) models. When fitting models across large spatial domains, it is increasingly likely that the effect of spatial and/or temporal covariates on species occurrence varies across space. Ignoring nonstationarity can result in reduced predictive performance and biased inference on the effect of covariates at different spatial regions (Finley 2011; Rollinson et al. 2021). SVC occupancy models are a flexible extension of spatial occupancy models that allows for not only the intercept to vary across space, but also allows the effects of the covariates themselves to vary spatially, resulting in each spatial location having a unique effect of the covariate (Pease, Pacifici, and Kays 2022). These models can generate key insights into how environmental factors differentially influence species across its range, how temporal trends vary across different parts of a species range, and the relative importance of different covariate effects at different parts of the species range.
Here we introduce four functions for fitting SVC models in spOccupancy. The function svcPGOcc() fits a single-season SVC occupancy model, and is an extension to the basic spatial occupancy model fit by spPGOcc(). The function svcTPGOcc() fits a multi-season SVC occupancy model, which serves as an extension to the spatio-temporal occupancy model fit by stPGOcc() where replicated detection-nondetection data are collected over multiple primary time periods (i.e., breeding seasons, years). Additionally, we include two functions for fitting SVC generalized linear models (GLMs) where we ignore imperfect detection: svcPGBinom() fits a single-season SVC GLM and svcTPGBinom() fits a multi-species GLM. In these functions, we also allow for modeling binomial data instead of the usual binary detection-nondetection data, which may be applicable for certain species or scenarios where imperfect detection may not be as large of an issue (e.g., modeling tree species distributions using a nested plot/subplot design).
As per usual, we use Pólya-Gamma data augmentation to yield a computationally efficient Gibbs sampler for GLM-type models with a logit link function (Polson, Scott, and Windle 2013), and use Nearest Neighbor Gaussian Processes (Datta et al. 2016) in all SVC models to greatly reduce the computational burden encountered when fitting models with spatial random effects formulated as a Gaussian process.
In this vignette, we will walk through each of the four single-species SVC models in spOccupancy, detailing how to fit the models, do model comparison using the Widely Available Information Criterion (WAIC) and k-fold cross-validation, as well as make predictions across an area of interest to generate maps of the spatially-varying coefficients. We will work with simulated data sets and will walk through how to simulate data sets using the spOccupancy functions simOcc(), simTOcc(), simBinom(), and simTBinom(). We will not go into explicit detail for some of the model-fitting function arguments (e.g., priors, initial values) as the syntax is nearly identical to other spOccupancy functions, so we encourage you to first work through the spOccupancy introductory vignette if you are not at least vaguely familiar with spOccupancy syntax.
Below, we first load the spOccupancy package, as well as the ggplot2 package to create some basic plots of our results. We also set a seed so you can reproduce our results.
library(spOccupancy)
library(ggplot2)
set.seed(300)Spatially-varying coefficients occupancy model
Basic model description
Let \(\boldsymbol{s}_j\) denote the spatial coordinates of site \(j\), where \(j = 1, \dots, J\). We define \(z(\boldsymbol{s}_j)\) as the true presence (1) or absence (0) at site \(j\) with spatial coordinates \(\boldsymbol{s}_j\). We model \(z(\boldsymbol{s}_j)\) as
\[\begin{equation}\label{z} z(\boldsymbol{s}_j) \sim \text{Bernoulli}(\psi(\boldsymbol{s}_j)), \end{equation}\]
where \(\psi(\boldsymbol{s}_j)\) is the occurrence probability at site \(j\). We model \(\psi(\boldsymbol{s}_j)\) according to
\[\begin{equation}\label{psi} \text{logit}(\psi(\boldsymbol{s}_j)) = \textbf{x}(\boldsymbol{s}_j)\boldsymbol{\beta} + \tilde{\textbf{x}}(\boldsymbol{s}_j)\textbf{w}(\boldsymbol{s}_j), \end{equation}\]
where \(\boldsymbol{\beta}\) is a vector of \(q\) regression coefficients (including an intercept) that describe the non-spatial effects of covariates \(\textbf{x}(\boldsymbol{s}_j)\), and \(\textbf{w}(\boldsymbol{s}_j)\) is a vector of \(\tilde{q}\) spatially-varying effects of covariates \(\tilde{\textbf{x}}(\boldsymbol{s}_j)\). Note that \(\tilde{\textbf{x}}(\boldsymbol{s}_j)\) may be identical to \(\textbf{x}(\boldsymbol{s}_j)\) if all covariates effects are assumed to vary spatially, or a subset of \(\textbf{x}(\boldsymbol{s}_j)\) if some effects are assumed to be constant across space. The model reduces to a traditional single-species occupancy model when all covariate effects are assumed constant across space and a spatial occupancy model (Johnson et al. 2013; Doser et al. 2022) when only the intercept is assumed to vary across space.
The spatially-varying effects \(\textbf{w}(\boldsymbol{s}_j)\) serve as local adjustments of the covariate effects at site \(j\) from the overall non-spatial effects \(\boldsymbol{\beta}\), resulting in the covariate having a unique effect on species occurrence at each site \(j\). Following Gelfand et al. (2003), we model each \(r = 1, \dots, \tilde{q}\) spatially-varying effect \(\text{w}_r(\boldsymbol{s}_j)\) using a zero-mean spatial Gaussian process. More specifically, we have
\[\begin{equation}\label{w} \text{$\text{w}_r(\boldsymbol{s})$} \sim N(\boldsymbol{0}, \boldsymbol{C}_r(\boldsymbol{s}, \boldsymbol{s}', \boldsymbol{\theta}_r)), \end{equation}\]
where \(\boldsymbol{C}_r(\boldsymbol{s}, \boldsymbol{s}', \boldsymbol{\theta}_r)\) is a \(J \times J\) covariance matrix that is a function of the distances between any pair of site coordinates \(\boldsymbol{s}\) and \(\boldsymbol{s}'\) and a set of parameters (\(\boldsymbol{\theta}_r\)) that govern the spatial process according to a spatial correlation function. Our associated software implementation in supports four correlation functions: exponential, spherical, Gaussian, and Matern (Banerjee, Carlin, and Gelfand 2003). For the exponential, spherical, and Gaussian correlation functions, \(\boldsymbol{\theta}_r = \{\sigma^2_r, \phi_r\}\), where \(\sigma^2_r\) is a spatial variance parameter and \(\phi_r\) is a spatial decay parameter. Large values of \(\sigma^2_r\) indicate large variation in the magnitude of a covariate effect across space, while values of \(\sigma^2_r\) close to 0 suggests little spatial variability in the magnitude of the effect. \(\phi_r\) controls the range of the spatial dependence in the covariate effect and is inversely related to the spatial range, such that when \(\phi_r\) is small, the covariate effect has a larger range of spatial dependence and varies more smoothly across space compared to larger values of \(\phi_r\). The Matern correlation function has an additional smoothness parameter \(\nu_r\), which provides further flexibility in the smoothness and decay of the spatial process. To avoid the computational burden associated with fitting the full Gaussian process model, we use NNGPs as a computationally efficient and statistically robust alternative (Datta et al. 2016). See the supplemental material in the spOccupancy manuscript and Datta et al. (2016) for more details on NNGPs and their implementations in occupancy models.
To account for imperfect detection in an occupancy modeling framework (MacKenzie et al. 2002, @tyre2003improving), \(k = 1, \dots, K(\boldsymbol{s}_j)\) sampling replicates are obtained at each site \(j\). We model the observed detection (1) or nondetection (0) of a study species at site \(j\), denoted \(y_k(\boldsymbol{s}_j)\), conditional on the true latent occupancy process \(z(\boldsymbol{s}_j)\) following
\[\begin{equation}\label{yDet} y_{k}(\boldsymbol{s}_j) \mid z(\boldsymbol{s}_j) \sim \text{Bernoulli}(p_{k}(\boldsymbol{s}_j)z(\boldsymbol{s}_j)), \end{equation}\]
where \(p_k(\boldsymbol{s})\) is the probability of detecting the species at site \(j\) during replicate \(k\) given the species is truly present at the site. We model detection probability as a function of site and/or observation-level covariates according to
\[\begin{equation}\label{pDet} \text{logit}(p_{k}(\boldsymbol{s}_j)) = \boldsymbol{v}_{k}(\boldsymbol{s}_j)\boldsymbol{\alpha}, \end{equation}\]
where \(\boldsymbol{\alpha}\) is a vector of regression coefficients (including an intercept) that describe the effect of site and/or observation covariates \(\boldsymbol{v}_{k}(\boldsymbol{s}_j)\) on detection.
To complete the Bayesian specification of the model, we assign Gaussian priors to the non-spatial occurrence and detection regression coefficients, an inverse-Gamma or uniform prior to the spatial variance parameters, and a uniform prior for the spatial decay (and smoothness if using a Matern correlation function) parameters.
Simulating data with simOcc()
The function simOcc() simulates single-species detection-nondetection data. simOcc() has the following arguments.
simOcc(J.x, J.y, n.rep, beta, alpha, psi.RE = list(),
p.RE = list(), sp = FALSE, svc.cols = 1, cov.model,
sigma.sq, phi, nu, x.positive = FALSE, ...)J.x and J.y indicate the number of spatial locations to simulate data along a horizontal and vertical axis, respectively, such that J.x * J.y is the total number of sites (i.e., J). n.rep is a numeric vector of length J that indicates the number of replicates at each of the J sites (denoted as K in the previous model description). beta and alpha are numeric vectors containing the intercept and any regression coefficient parameters for the occurrence and detection portions of the occupancy model, respectively. psi.RE and p.RE are lists that are used to specify random intercepts on occurrence and detection, respectively. These are only specified when we want to simulate data with unstructured random intercepts. Each list should be comprised of two tags: levels, a vector that specifies the number of levels for each random effect included in the model, and sigma.sq.psi or sigma.sq.p, which specify the variances of the random effects for each random effect included in the model. sp is a logical value indicating whether to simulate data with at least one spatial Gaussian process for either the intercept or some of the occupancy covariate effects. svc.cols is a numeric vector indicating which of the covariates (including the intercept) are generated with spatially-varying effects. By default, svc.cols = 1, which corresponds to a spatially-varying intercept when sp = TRUE (i.e., a spatial occupancy model). cov.model specifies the covariance function used to model the spatial dependence structure, with supported values of exponential, matern, spherical, and gaussian. Finally, sigma.sq is the spatial variance parameter, phi is the spatial range parameter, and nu is the spatial smoothness parameter (only applicable when cov.model = 'matern'). Note that sigma.sq, phi, and nu should have the same length as the number of spatially-varying effects specified in svc.cols. Lastly, the x.positive argument indicates whether or not the occupancy covariates should be simulated to be only positive from a Uniform(0, 1) distribution (TRUE) or both positive and negative and simulated from a Normal(0, 1) distribution (FALSE).
Below we simulate data across 1600 sites with anywhere between 1-4 replicates at a given site, a single covariate effect on occurrence, and a single covariate effect on detection. We assume both the occupancy intercept and the effect of the covariate vary across space, so we set svc.cols = c(1, 2). We use a spherical correlation function. We do not include any unstructured random effects on occurrence or detection.
J.x <- 40
J.y <- 40
# Total number of sites
(J <- J.x * J.y)[1] 1600
# Number of replicates at each site
n.rep <- sample(4:4, J, replace = TRUE)
# Intercept and covariate effect on occurrence
# Note these are the non-spatial effects.
beta <- c(-0.5, -0.2)
# Intercept and covariate effect on detection
alpha <- c(0.9, -0.3)
# No unstructured random intercept on occurrence
psi.RE <- list()
# No unstructured random intercept on detection
p.RE <- list()
# Spatial range for intercept and covariate effect
phi <- c(3 / .8, 3 / .7)
# Spatial variance for intercept and covariate effect
sigma.sq <- c(1, 0.5)
# Simulate the occupancy covariate from a Normal(0, 1) distribution
x.positive <- FALSE
# Spatially-varying coefficient columns
svc.cols <- c(1, 2)
# Simulate the data
dat <- simOcc(J.x = J.x, J.y = J.y, n.rep = n.rep, beta = beta,
alpha = alpha, psi.RE = psi.RE, p.RE = p.RE,
sp = TRUE, sigma.sq = sigma.sq, phi = phi,
cov.model = 'spherical', svc.cols = svc.cols,
x.positive = x.positive)Next, let’s explore the simulated data a bit before we move on (plotting code adapted from Hooten and Hefley (2019)).
str(dat)List of 13
$ X : num [1:1600, 1:2] 1 1 1 1 1 1 1 1 1 1 ...
$ X.p : num [1:1600, 1:4, 1:2] 1 1 1 1 1 1 1 1 1 1 ...
$ coords : num [1:1600, 1:2] 0 0.0256 0.0513 0.0769 0.1026 ...
..- attr(*, "dimnames")=List of 2
.. ..$ : NULL
.. ..$ : chr [1:2] "Var1" "Var2"
$ w : num [1:1600, 1:2] -1.184 -1.456 -0.633 -0.124 -0.206 ...
$ psi : num [1:1600, 1] 0.144 0.186 0.246 0.36 0.344 ...
$ z : int [1:1600] 0 0 0 0 1 1 0 0 0 0 ...
$ p : num [1:1600, 1:4] 0.672 0.738 0.581 0.737 0.697 ...
$ y : int [1:1600, 1:4] 0 0 0 0 1 1 0 0 0 0 ...
$ X.p.re : logi NA
$ X.re : logi NA
$ X.w : num [1:1600, 1:2] 1 1 1 1 1 1 1 1 1 1 ...
$ alpha.star: logi NA
$ beta.star : logi NAThe simulated data object consists of the following objects: X (the occurrence design matrix), X.p (the detection design matrix), coords (the spatial coordinates of each site), w (the latent spatial process for any covariates (and intercept) whose effects vary across space), psi (occurrence probability), z (the latent occupancy status), y (the detection-nondetection data), X.w (the design matrix for the spatially-varying coefficients), X.p.re (the detection random effect levels for each site), X.re (the occurrence random effect levels for each site), alpha.star (the detection random effects for each level of the random effect), beta.star (the occurrence random effects for each level of the random effect). Note because we did not include any unstructured effects on detection or occurrence, the objects associated with the unstructured random effects all have a value of NA.
# Detection-nondetection data
y <- dat$y
# Occurrence design matrix for fixed effects
X <- dat$X
# Detection design matrix for fixed effets
X.p <- dat$X.p
# Occurrence values
psi <- dat$psi
# Spatial coordinates
coords <- dat$coords
# Spatially-varying intercept and covariate effects
w <- dat$w
# Simple plot of the occurrence probability across space.
# Dark points indicate high occurrence.
plot(coords, type = "n", xlab = "", ylab = "", asp = TRUE,
main = "Simulated Occurrence", bty = 'n')
points(coords, pch=15, cex = 2.1, col = rgb(0,0,0,(psi-min(psi))/diff(range(psi))))
We see there is clear variation in occurrence probability across the simulated spatial region.
We can also visualize the spatially-varying intercept and spatially-varying covariate effect. We do that by adding the spatial component of the intercept and covariate effect (stored in the w matrix) to the non-spatial components of the intercept and covariate effect (stored in beta), and then visualizing using the same code as before
# Intercept
int.effect <- beta[1] + w[, 1]
cov.effect <- beta[2] + w[, 2]
# Dark points indicate more positive effects, white points indicate more negative effects.
plot(coords, type = "n", xlab = "", ylab = "", asp = TRUE, main = "Intercept",
bty = 'n')
points(coords, pch=15, cex = 2.1,
col = rgb(0,0,0,(int.effect-min(int.effect))/diff(range(int.effect))))
plot(coords, type = "n", xlab = "", ylab = "", asp = TRUE, main = "Covariate Effect",
bty = 'n')
points(coords, pch=15, cex = 2.1,
col = rgb(0,0,0,(cov.effect-min(cov.effect))/diff(range(cov.effect))))
Note the spatial intercept corresponds fairly closely with the map of occurrence probability, which makes sense.
The final step before we can fit the model is to package up the data in a list for use in spOccupancy model fitting functions. This requires creating a list that consists of the detection-nondetection data (y), occurrence covariates (occ.covs), detection covariates (det.covs), and coordinates (coords). See the introductory vignette for more details. For our example here (and throughout the vignette), we will fit the model to 75% of the data points (1200 locations) and subsequently predict at the remaining 400 values to show the predictive ability of the model.
# Subset data for prediction.
# Split into fitting and prediction data set
pred.indx <- sample(1:J, round(J * .25), replace = FALSE)
y.fit <- y[-pred.indx, ]
y.pred <- y[pred.indx, ]
X.fit <- X[-pred.indx, ]
X.pred <- X[pred.indx, ]
X.p.fit <- X.p[-pred.indx, , ]
X.p.pred <- X.p[pred.indx, , ]
coords.fit <- coords[-pred.indx, ]
coords.pred <- coords[pred.indx, ]
psi.fit <- psi[-pred.indx]
psi.pred <- psi[pred.indx]
w.fit <- w[-pred.indx, ]
w.pred <- w[pred.indx, ]
# Package all data into a list
# Occurrence covariates
occ.covs <- X.fit[, 2, drop = FALSE]
colnames(occ.covs) <- c('occ.cov.1')
# Detection covariates
det.covs <- list(det.cov.1 = X.p.fit[, , 2])
# Package into a list for spOccupancy
data.list <- list(y = y.fit,
occ.covs = occ.covs,
det.covs = det.covs,
coords = coords.fit)
# Take a look at the data structure.
str(data.list)List of 4
$ y : int [1:1200, 1:4] 0 0 0 1 1 0 0 0 0 1 ...
$ occ.covs: num [1:1200, 1] 0.536 1.828 0.125 -0.364 1.433 ...
..- attr(*, "dimnames")=List of 2
.. ..$ : NULL
.. ..$ : chr "occ.cov.1"
$ det.covs:List of 1
..$ det.cov.1: num [1:1200, 1:4] 0.612 -0.446 1.916 0.226 -0.813 ...
$ coords : num [1:1200, 1:2] 0 0.0256 0.0513 0.1026 0.1282 ...
..- attr(*, "dimnames")=List of 2
.. ..$ : NULL
.. ..$ : chr [1:2] "Var1" "Var2"Fitting spatially-varying coefficients occupancy models with svcPGOcc()
The function svcPGOcc() fits single-season SVC occupancy models. svcPGOcc() has the following arguments:
svcPGOcc(occ.formula, det.formula, data, inits, priors,
tuning, svc.cols = 1, cov.model = "exponential", NNGP = TRUE,
n.neighbors = 15, search.type = "cb", n.batch,
batch.length, accept.rate = 0.43,
n.omp.threads = 1, verbose = TRUE, n.report = 100,
n.burn = round(.10 * n.batch * batch.length),
n.thin = 1, n.chains = 1, k.fold, k.fold.threads = 1,
k.fold.seed = 100, k.fold.only = FALSE, ...)The arguments to svcPGOcc() are identical as those for a spatial occupancy model fit using spPGOcc(), with the addition of the svc.cols argument to specify the SVCs in the model. occ.formula and det.formula contain the R model formulas for the occurrence and detection portions of the occupancy model. Here we fit the model with a single covariate on both occupancy and detection
occ.formula <- ~ occ.cov.1
det.formula <- ~ det.cov.1The svc.cols argument is used to specify the covariates whose effects are estimated as SVCs. svc.cols can either be a numeric indexing vector with integer numbers corresponding to the order in which you specified covariates in the occ.formula argument, or can be a character vector with the names of the covariates specified in occ.formula. Note that by default, svc.cols = 1, which is equivalent to fitting a spatial occupancy model with spPGOcc(). This clearly shows how SVC occupancy models are a simple extension of a spatial occupancy model. A spatial occupancy model is simply an SVC occupancy model where the only spatially varying “covariate” is the intercept. If specifying svc.cols as a character vector, use '(Intercept)' to specify a spatially-varying intercept. Here we set svc.cols to include a spatially-varying intercept and spatially-varying effect of the occurrence covariate. We also specify the cov.model argument to indicate we will use a spherical correlation function (note spOccupancy uses the same correlation function for all SVCs). Usually, our default correlation function is the exponential, which is what we use throughout most of the other vignettes, but here we use a spherical correlation function to show how spOccupancy can handle these other functions.
svc.cols <- c(1, 2)
# OR
# svc.cols <- c('(Intercept)', 'occ.cov.1')
cov.model <- 'spherical'We next specify the initial values, which are specified in a list analogous to a spatial occupancy model using spPGOcc, with the only difference being that if you supply initial values for the spatial random effects w, these must be specified as a two-dimensional matrix with rows corresponding to SVC and column corresponding to site.
dist.data <- dist(data.list$coords)
inits.list <- list(alpha = 0, beta = 0, sigma.sq = 0.5,
phi = 3 / mean(dist.data),
z = apply(data.list$y, 1, max, na.rm = TRUE),
w = matrix(0, length(svc.cols), nrow(data.list$y)))We next specify the priors to use for all parameters in the model (alternatively, we could not specify priors and simply use the default values svcPGOcc() provides). We will use the default normal priors for the occurrence (beta) and detection (alpha) regression coefficients. For the spatial decay parameter (phi), we specify a uniform prior with bounds based on the maximum and minimum inter-site distances. Our default prior for phi is to set the lower bound to 3 / max and upper bound to 3 / min, where max and min are the maximum and minimum inter-site distances, respectively. This results in a prior that states the effective spatial range is anywhere between the maximum distance between sites and the smallest distance between sites. Lastly, we specify an inverse-Gamma prior for sigma.sq. Following Banerjee, Carlin, and Gelfand (2003), we generally will set the scale parameter of the inverse-Gamma to 2 and the shape parameter to our buest guess of the spatial variance. We could also specify a uniform prior for the spatial variance parameter. For binary data, very large values of sigma.sq can result in undesirable and unrealistic values of the spatial random effects on the logit scale, and so a uniform prior can be used to restrict sigma.sq to some maximum value (e.g., 5) that is reasonable on the logit scale (Wright et al. 2021).
priors.list <- list(alpha.normal = list(mean = 0, var = 2.72),
beta.normal = list(mean = 0, var = 2.72),
sigma.sq.ig = list(a = 2, b = 0.5),
phi.unif = list(a = 3 / max(dist.data),
b = 3 / min(dist.data)))The next three arguments (n.batch, batch.length, and accept.rate) are all related to the Adaptive MCMC sampler used when we fit the model. Updates for all parameters with a uniform prior (in this case the spatial decay parameter phi and the spatial variance parameter sigma.sq) require the use of a Metropolis-Hastings algorithm. We implement an adaptive Metropolis-Hastings algorithm as discussed in Roberts and Rosenthal (2009). This algorithm adjusts the tuning values for each parameter that requires a Metropolis-Hastings update within the sampler itself. This process results in a more efficient sampler than if we were to fix the tuning parameters prior to fitting the model. The parameter accept.rate is the target acceptance rate for each parameter, and the algorithm will adjust the tuning parameters to hover around this value. The default value is 0.43, which we suggest leaving as is unless you have a good reason to change it. The tuning parameters are updated after a single “batch”. We break up the total number of MCMC samples into a set of “batches”, where each batch has a specific number of samples. We must specify both the total number of batches (n.batch) as well as the number of MCMC samples each batch contains (batch.length). Thus, the total number of MCMC samples is n.batch * batch.length. Typically, we set batch.length = 25 and then play around with n.batch until convergence is reached. Here we set n.batch = 800 for a total of 20000 MCMC samples. We will additionally specify a burn-in period of length 10000 and a thinning rate of 10. We run the model for 3 chains, ultimately resulting in 3000 posterior samples. Importantly, we also need to specify an initial value for the tuning parameters for the spatial decay parameter, spatial variance parameter if using a uniform prior for sigma.sq, and the smoothness parameter (if cov.model = 'matern'). These values are supplied as input in the form of a list with tags phi, sigma.sq, and nu. The initial tuning value can be any value greater than 0, but we recommend starting the value out around 0.5. After some initial runs of the model, if you notice the final acceptance rate of a parameter is much larger or smaller than the target acceptance rate (accept.rate), you can then change the initial tuning value to get closer to the target rate. Here we set the initial tuning value for phi to 0.2 after some initial exploratory runs of the model.
batch.length <- 25
n.batch <- 800
n.burn <- 10000
n.thin <- 10
n.chains <- 1
tuning.list <- list(phi = 0.2)We are now ready to run the model. We set the verbose argument equal to TRUE and the n.report argument to 100 to report progress on the MCMC chain after every 100th batch. Additionally, we fit the model with an NNGP (NNGP = TRUE) using 5 neighbors (n.neighbors = 5). See the supplemental material in Doser et al. (2022) for more information on choosing the number of neighbors in the NNGP approximation.
n.omp.threads <- 1
verbose <- TRUE
n.report <- 100 # Report progress at every 100th batch.
# Approx. run time: 4.2 min
out.svc <- svcPGOcc(occ.formula = occ.formula,
det.formula = det.formula,
data = data.list,
inits = inits.list,
n.batch = n.batch,
batch.length = batch.length,
priors = priors.list,
svc.cols = svc.cols,
cov.model = cov.model,
NNGP = TRUE,
n.neighbors = 5,
tuning = tuning.list,
n.report = n.report,
n.burn = n.burn,
n.thin = n.thin,
n.chains = n.chains)----------------------------------------
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 1200 sites.
Samples per chain: 20000 (800 batches of length 25)
Burn-in: 10000
Thinning Rate: 10
Number of Chains: 1
Total Posterior Samples: 1000
Number of spatially-varying coefficients: 2
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: 100 of 800, 12.50%
Coefficient Parameter Acceptance Tuning
1 phi 52.0 0.18279
2 phi 52.0 0.21025
-------------------------------------------------
Batch: 200 of 800, 25.00%
Coefficient Parameter Acceptance Tuning
1 phi 44.0 0.17917
2 phi 24.0 0.15268
-------------------------------------------------
Batch: 300 of 800, 37.50%
Coefficient Parameter Acceptance Tuning
1 phi 40.0 0.15268
2 phi 20.0 0.17562
-------------------------------------------------
Batch: 400 of 800, 50.00%
Coefficient Parameter Acceptance Tuning
1 phi 40.0 0.14378
2 phi 64.0 0.16539
-------------------------------------------------
Batch: 500 of 800, 62.50%
Coefficient Parameter Acceptance Tuning
1 phi 52.0 0.14378
2 phi 40.0 0.18648
-------------------------------------------------
Batch: 600 of 800, 75.00%
Coefficient Parameter Acceptance Tuning
1 phi 32.0 0.13815
2 phi 32.0 0.18648
-------------------------------------------------
Batch: 700 of 800, 87.50%
Coefficient Parameter Acceptance Tuning
1 phi 44.0 0.15268
2 phi 44.0 0.28381
-------------------------------------------------
Batch: 800 of 800, 100.00%We can take a look at the model results using the summary() function and compare them to the true values we used to simulate the data.
summary(out.svc)
Call:
svcPGOcc(occ.formula = occ.formula, det.formula = det.formula,
data = data.list, inits = inits.list, priors = priors.list,
tuning = tuning.list, svc.cols = svc.cols, cov.model = cov.model,
NNGP = TRUE, n.neighbors = 5, n.batch = n.batch, batch.length = batch.length,
n.report = n.report, n.burn = n.burn, n.thin = n.thin, n.chains = n.chains)
Samples per Chain: 20000
Burn-in: 10000
Thinning Rate: 10
Number of Chains: 1
Total Posterior Samples: 1000
Run Time (min): 8.0548
Occurrence (logit scale):
Mean SD 2.5% 50% 97.5% Rhat ESS
(Intercept) -0.4291 0.1995 -0.7904 -0.4434 -0.0070 NA 91
occ.cov.1 -0.1082 0.2259 -0.5427 -0.1135 0.4487 NA 16
Detection (logit scale):
Mean SD 2.5% 50% 97.5% Rhat ESS
(Intercept) 0.8751 0.0701 0.7384 0.8734 1.0106 NA 1000
det.cov.1 -0.2621 0.0688 -0.3968 -0.2606 -0.1325 NA 1000
Spatial Covariance:
Mean SD 2.5% 50% 97.5% Rhat ESS
sigma.sq-(Intercept) 1.5021 0.4981 0.8141 1.4217 2.7294 NA 41
sigma.sq-occ.cov.1 0.5076 0.1856 0.2048 0.4835 0.9883 NA 50
phi-(Intercept) 6.9408 1.4537 4.1067 6.8634 9.8749 NA 52
phi-occ.cov.1 5.2571 2.0523 2.2056 5.1009 9.8257 NA 18
# True values
beta[1] -0.5 -0.2
alpha[1] 0.9 -0.3
sigma.sq[1] 1.0 0.5
phi[1] 3.750000 4.285714Because we only ran one chain of the model, we see the Rhat values are reported as NA. For a complete analysis, we would run the model for multiple chains, make sure the Rhat values are less than 1.1, and also ensure the effective sample sizes are adequately large. Here, the ESS values are somewhat low for the occurrence parameters and spatial covariance parameters, but we will continue interpreting the results for our exploratory purposes here. We see our model does a good job of recovering the true occurrence and detection regression coefficient values. The spatial variance parameters are also quite close to the estimated values. The spatial decay parameter value for the intercept is a bit larger than the simulated value. The spatial decay parameters are only weakly identifiable (i.e., there is very little information to estimate them), and thus estimating their true values can be a difficult task, in particular when fitting a model with multiple SVCs. Generally, we do not attempt to interpret the spatial decay parameters when fitting spatially-explicit occupancy models. Instead, we will often interpret the actual estimated spatial process values at each location, which are of particular interest in spatially-varying coefficient models.
Next, let’s take a look at the resulting objects contained in the out.svc list.
names(out.svc) [1] "rhat" "beta.samples" "alpha.samples" "theta.samples"
[5] "coords" "z.samples" "X" "X.re"
[9] "X.w" "w.samples" "psi.samples" "like.samples"
[13] "X.p" "X.p.re" "y" "ESS"
[17] "call" "n.samples" "n.neighbors" "cov.model.indx"
[21] "svc.cols" "type" "n.post" "n.thin"
[25] "n.burn" "n.chains" "pRE" "psiRE"
[29] "run.time" The resulting model object contains a variety of things, most of which are just used in subsequent functions for posterior predictive checks, prediction, and summarization. The objects that end in “samples” are the posterior MCMC samples for the different objects. See ?svcPGOcc for more information.
To extract the estimates of the spatially-varying coefficients at each of the spatial locations in the data set used to fit the model, we need to combine the non-spatial component of the coefficient (contained in out.svc$beta.samples) and the spatial component of the coefficient (contained in out.svc$w.samples). Recall that in an SVC occupancy model, the total effect of a covariate at any given location is the sum of the non-spatial effect and the adjustment of the effect at that specific location. We provide the function getSVCSamples() to extract the SVCs at each location.
svc.samples <- getSVCSamples(out.svc)
str(svc.samples)List of 2
$ (Intercept): 'mcmc' num [1:1000, 1:1200] -1.497 -0.158 -0.62 -0.127 -0.355 ...
..- attr(*, "mcpar")= num [1:3] 1 1000 1
$ occ.cov.1 : 'mcmc' num [1:1000, 1:1200] 0.6315 -0.4277 -0.7302 -0.6828 0.0796 ...
..- attr(*, "mcpar")= num [1:3] 1 1000 1The resulting object, here called svc.samples, is a list with each component corresponding to a matrix of the MCMC samples of each spatially-varying coefficient estimated in the model, with rows corresponding to MCMC sample and column corresponding to site.
Below we plot the true SVCs for the covariate at the 1200 locations used for fitting the model compared to the mean estimates from our model.
# Get true covariate values at the locations used to fit the mdoel
cov.effect.fit <- beta[2] + w.fit[, 2]
# Get mean values of the SVC for the covariate
svc.cov.mean <- apply(svc.samples$occ.cov.1, 2, mean)
# Dark points indicate more positive effects, white points
# indicate more negative effects.
plot(coords.fit, type = "n", xlab = "", ylab = "", asp = TRUE,
main = "True values", bty = 'n')
points(coords.fit, pch=15, cex = 2.1,
col = rgb(0,0,0,(cov.effect.fit-min(cov.effect.fit))/diff(range(cov.effect.fit))))
plot(coords.fit, type = "n", xlab = "", ylab = "", asp = TRUE,
main = "Estimated values", bty = 'n')
points(coords.fit, pch=15, cex = 2.1,
col = rgb(0,0,0,(svc.cov.mean-min(svc.cov.mean))/diff(range(svc.cov.mean))))
Note that because we held out 400 random values across the study area, some of the white squares in the above image correspond to locations where we did not have any data (we will subsequently predict at these locations). We see our estimates align pretty closely with the true values used to simulate the data. Our mean estimates are smoother than the values used to generate the data. This is what we would expect, as the true values are a single instance of a simulated spatial process, whereas the mean values we have plotted average across individual instances to generate a more smoothed estimate. Overall, the model seems to accurately identify locations of low and high effects of the covariate.
Posterior Predictive Checks
The spOccupancy function ppcOcc() performs a posterior predictive check for all spOccupancy model objects as an assessment of Goodness of Fit (GoF). The key idea of GoF testing is that a good model should generate data that closely align with the observed data. If there are large differences in the observed data from the model-generated data, our model is likely not very useful (Hooten and Hobbs 2015). We can use the ppcOcc() and summary() functions to generate a Bayesian p-value as a quick assessment of model fit. A Bayesian p-value that hovers around 0.5 indicates adequate model fit, while values less than 0.1 or greater than 0.9 suggest our model does not fit the data well. See the introductory spOccupancy vignette and help page for ppcOcc() for more details. Below we perform a posterior predictive check with the Freeman-Tukey statistic, grouping the data by individual sites.
Call:
ppcOcc(object = out.svc, fit.stat = "freeman-tukey", group = 1)
Samples per Chain: 20000
Burn-in: 10000
Thinning Rate: 10
Number of Chains: 1
Total Posterior Samples: 1000
Bayesian p-value: 0.496
Fit statistic: freeman-tukey We see our Bayesian p-value is very close to the optimal 0.5, suggesting adequate model fit (it would be a bit concerning if we didn’t see a good model fit here since we are using simulated data!).
Model Selection using WAIC
The spOccupancy function waicOcc() calculates the Widely Applicable Information Criterion (WAIC) for all spOccupancy fitted model objects. The WAIC is a fully Bayesian information criterion that is adequate for comparing a set of hierarchical models and selecting the best-performing model for final analysis (see the introductory spOccupancy vignette for more details). Smaller values of WAIC indicate a better performing model. Below, we fit a spatial occupancy model without the SVC of the covariate effect using the spOcc() function (we could equivalently do this by using svcPGOcc() and setting svc.cols = 1).
# Using default priors and initial values.
# Approx. run time: 4.2 min
out.sp <- spPGOcc(occ.formula = occ.formula,
det.formula = det.formula,
data = data.list,
n.batch = n.batch,
batch.length = batch.length,
cov.model = cov.model,
NNGP = TRUE,
n.neighbors = 5,
tuning = tuning.list,
n.report = n.report,
n.burn = n.burn,
n.thin = n.thin,
n.chains = n.chains)----------------------------------------
Preparing to run the model
----------------------------------------No prior specified for beta.normal.
Setting prior mean to 0 and prior variance to 2.72No prior specified for alpha.normal.
Setting prior mean to 0 and prior variance to 2.72No prior specified for phi.unif.
Setting uniform bounds based on the range of observed spatial coordinates.No prior specified for sigma.sq.
Using an inverse-Gamma prior with the shape parameter set to 2 and scale parameter to 1.z.inits is not specified in initial values.
Setting initial values based on observed databeta is not specified in initial values.
Setting initial values to random values from the prior distributionalpha is not specified in initial values.
Setting initial values to random values from the prior distributionphi is not specified in initial values.
Setting initial value to random value from the prior distributionsigma.sq is not specified in initial values.
Setting initial value to random value from the prior distributionw is not specified in initial values.
Setting initial value to 0----------------------------------------
Building the neighbor list
----------------------------------------
----------------------------------------
Building the neighbors of neighbors list
----------------------------------------
----------------------------------------
Model description
----------------------------------------
NNGP Occupancy model with Polya-Gamma latent
variable fit with 1200 sites.
Samples per chain: 20000 (800 batches of length 25)
Burn-in: 10000
Thinning Rate: 10
Number of Chains: 1
Total Posterior Samples: 1000
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: 100 of 800, 12.50%
Parameter Acceptance Tuning
phi 44.0 0.15576
-------------------------------------------------
Batch: 200 of 800, 25.00%
Parameter Acceptance Tuning
phi 48.0 0.16212
-------------------------------------------------
Batch: 300 of 800, 37.50%
Parameter Acceptance Tuning
phi 32.0 0.16212
-------------------------------------------------
Batch: 400 of 800, 50.00%
Parameter Acceptance Tuning
phi 40.0 0.14378
-------------------------------------------------
Batch: 500 of 800, 62.50%
Parameter Acceptance Tuning
phi 60.0 0.15891
-------------------------------------------------
Batch: 600 of 800, 75.00%
Parameter Acceptance Tuning
phi 36.0 0.16873
-------------------------------------------------
Batch: 700 of 800, 87.50%
Parameter Acceptance Tuning
phi 40.0 0.16873
-------------------------------------------------
Batch: 800 of 800, 100.00%
# WAIC for the SVC model
waicOcc(out.svc) elpd pD WAIC
-1185.427 151.770 2674.394
# WAIC for the spatial occupancy model
waicOcc(out.sp) elpd pD WAIC
-1232.0892 118.5872 2701.3528 We see the WAIC for the SVC occupancy model is lower than the WAIC for the spatial occupancy model, indicating the SVC model has improved model performance.
Prediction
Finally, we can use the predict() function with all spOccupancy model-fitting functions to generate a series of posterior predictive samples at new locations (as well as already sampled locations), given a set of covariates and their spatial locations. Note that we can predict both new occupancy values as well as new detection values.
Below we predict occupancy probability at the 400 locations we held out when fitting the model. The predict() function for svcPGOcc() objects requires the model object, the design matrix of covariates at the new locations (including the intercept if specified in the model), and the spatial coordinates of the new locations. Below we predict across the 400 “new” locations and plot them in comparison to the true values we used to simulate the data.
# Take a look at X.pred, the design matrix for the prediction locations
head(X.pred) [,1] [,2]
[1,] 1 1.1569942
[2,] 1 -0.8225662
[3,] 1 -1.2718065
[4,] 1 2.6700613
[5,] 1 -0.4581872
[6,] 1 0.1231113
# Predict occupancy at the 400 new sites
out.pred <- predict(out.svc, X.pred, coords.pred)----------------------------------------
Prediction description
----------------------------------------
NNGP Occupancy model with Polya-Gamma latent
variable fit with 1200 observations.
Number of covariates 2 (including intercept if specified).
Number of spatially-varying covariates 2 (including intercept if specified).
Using the spherical spatial correlation model.
Using 5 nearest neighbors.
Number of MCMC samples 1000.
Predicting at 400 non-sampled locations.
Source compiled with OpenMP support and model fit using 1 threads.
-------------------------------------------------
Predicting
-------------------------------------------------
Location: 100 of 400, 25.00%
Location: 200 of 400, 50.00%
Location: 300 of 400, 75.00%
Location: 400 of 400, 100.00%
Generating latent occupancy state
# Use the getSVCSamples() function to extract the SVC values
# at the prediction locations
svc.pred.samples <- getSVCSamples(out.svc, pred.object = out.pred)
# True covariate effect values at new locations
cov.effect.pred <- beta[2] + w.pred[, 2]
# Get mean values of the SVC for the covariate
svc.cov.pred.mean <- apply(svc.pred.samples$occ.cov.1, 2, mean)
# Dark points indicate more positive effects, white points indicate more
# negative effects.
plot(coords.pred, type = "n", xlab = "", ylab = "", asp = TRUE,
main = "True values", bty = 'n')
points(coords.pred, pch=15, cex = 2.1,
col = rgb(0,0,0,(cov.effect.pred-min(cov.effect.pred))/
diff(range(cov.effect.pred))))
plot(coords.pred, type = "n", xlab = "", ylab = "", asp = TRUE,
main = "Estimated values", bty = 'n')
points(coords.pred, pch=15, cex = 2.1,
col = rgb(0,0,0,(svc.cov.pred.mean-min(svc.cov.pred.mean))/
diff(range(svc.cov.pred.mean))))
We see a pretty close correspondence between the true values and the predicted values.
Spatially-varying coefficients generalized linear model
Accounting for imperfect detection explicitly in an occupancy modeling framework may not be feasible given the data. In particular, single-visit (or nonreplicated) detection-nondetection data are a very common data source due to the extra expenses of doing multiple surveys at a given location. In this case, we cannot separately estimate imperfect detection from occupancy probability in an occupancy modeling framework without making strict assumptions (Lele, Moreno, and Bayne 2012), and so instead we might fit a more standard generalized linear model (GLM).
Additionally, occupancy models are much less frequently used for modeling plant species distributions than animal distributions. While imperfect detection is still prevalent in plant distribution studies (Chen et al. 2013), we may wish to model plant distribution data using a standard generalized linear modeling framework, in particular if we do not have replicate surveys available in our data set. If our inference focuses solely on the effects of covariates on where a species occurs and not the overall estimates of “occupancy probability”, this may be a feasible approach. However, all inferences made from such models that ignore imperfect detection should carefully consider how detection probability may vary across space and/or time, and what consequences this has on the interpretation of the results. For example, if a spatially-varying covariate influences both the probability a species occurs at a location and the probability the species is detected, one would interpret the effect of the covariate as the effect on the confounded process of detection/occupancy probability.
Basic model description
As before, let \(y_k(\boldsymbol{s}_j)\) be the detection (1) or non-detection (0) of our species of interest at site \(j\) with spatial coordinates \(\boldsymbol{s}_j\). When not explicitly accounting for imperfect detection in an occupancy model, we will model the detection-nondetection data directly in a SVC GLM framework. Specifically, we have
\[\begin{equation}\label{yNoDet} y^*(\boldsymbol{s}_j) \sim \text{Binomial}(K(\boldsymbol{s}_j), \psi(\boldsymbol{s}_j)), \end{equation}\]
where \(y^*(\boldsymbol{s}_j) = \sum_{k = 1}^{K(\boldsymbol{s}_j)}y_{k}(\boldsymbol{s}_j)\), \(K(\boldsymbol{s}_j)\) is the number of replicate surveys at site \(j\), and \(\boldsymbol{\psi}(\boldsymbol{s}_j)\) is the occurrence probability at site \(j\). When only one replicate survey is available at each site, the Binomial likelihood in the previous equation reduces to a Bernoulli likelihood. We model \(\psi(\boldsymbol{s}_j)\) as before. Note that we can incorporate spatially-varying covariates in the model of \(\psi(\boldsymbol{s}_j)\) that may influence detection probability of the species, but our estimates of \(\psi(\boldsymbol{s}_j)\) are still interpreted as relative occurrence probabilities.
Simulating data with simBinom()
The function simBinom() simulates single-species detection-nondetection data for which detection is assumed perfect. simBinom() has the following arguments, which are very similar to those we saw previously with simOcc().
simBinom(J.x, J.y, weights, beta, psi.RE = list(),
sp = FALSE, svc.cols = 1, cov.model, sigma.sq, phi, nu,
x.positive = FALSE, ...)J.x and J.y indicate the number of spatial locations to simulate data along a horizontal and vertical axis, respectively, such that J.x * J.y is the total number of sites (i.e., J). weights is a numeric vector of length J that indicates the number of Bernoulli trials (replicates) at each of the J sites (denoted as K in the previous model description). beta is a numeric vector containing the intercept and any regression coefficient parameters for the model, respectively. psi.RE is a list that is used to specify unstructured random intercepts, respectively. These are only specified when we want to simulate data with random intercepts. All other arguments are the same as those we saw simOcc().
We simulate data across 1600 sites where we assume there is only a single replicate survey at a given site and a single covariate effect on relative occurrence. We assume both the intercept and the effect of the covariate vary across space, so we set svc.cols = c(1, 2). We use an exponential correlation function. We do not include any unstructured random effects on occurrence or detection.
# Set seed again to get the same data set
set.seed(488)
J.x <- 40
J.y <- 40
# Total number of sites
(J <- J.x * J.y)[1] 1600
# Number of replicates at each site
weights <- rep(1, J)
# Intercept and covariate effect
# Note these are the non-spatial effects.
beta <- c(-0.5, -0.2)
# No unstructured random intercepts
psi.RE <- list()
# Spatial range for intercept and covariate effect
phi <- c(3 / .8, 3 / .7)
# Spatial variance for intercept and covariate effect
sigma.sq <- c(1, 0.5)
# Simulate the covariate from a Normal(0, 1) distribution
x.positive <- FALSE
# Spatially-varying coefficient columns
svc.cols <- c(1, 2)
# Simulate the data
dat <- simBinom(J.x = J.x, J.y = J.y, weights = weights, beta = beta,
psi.RE = psi.RE, sp = TRUE, sigma.sq = sigma.sq, phi = phi,
cov.model = 'exponential', svc.cols = svc.cols,
x.positive = x.positive)Next, let’s explore the simulated data a bit before we move on (plotting code adapted from Hooten and Hefley (2019)).
str(dat)List of 8
$ X : num [1:1600, 1:2] 1 1 1 1 1 1 1 1 1 1 ...
$ coords : num [1:1600, 1:2] 0 0.0256 0.0513 0.0769 0.1026 ...
..- attr(*, "dimnames")=List of 2
.. ..$ : NULL
.. ..$ : chr [1:2] "Var1" "Var2"
$ w : num [1:1600, 1:2] 0.61481 0.00417 -0.26095 -0.54444 -0.9741 ...
$ psi : num [1:1600, 1] 0.539 0.362 0.402 0.338 0.329 ...
$ y : int [1:1600] 0 1 0 1 0 0 0 0 1 0 ...
$ X.re : logi NA
$ X.w : num [1:1600, 1:2] 1 1 1 1 1 1 1 1 1 1 ...
$ beta.star: logi NAThe simulated data object consists of the following objects: X (the design matrix), coords (the spatial coordinates of each site), w (the latent spatial process for any covariates (and intercept) whose effects vary across space), psi (relative occurrence probability), y (the detection-nondetection data), X.w (the design matrix for the spatially-varying coefficients), X.re (the occurrence random effect levels for each site), beta.star (the random effects for each level of the unstructured random effect). Note because we did not include any unstructured effects, the objects associated with the unstructured random effects all have a value of NA.
# Detection-nondetection data
y <- dat$y
# Design matrix
X <- dat$X
# Occurrence values
psi <- dat$psi
# Spatial coordinates
coords <- dat$coords
# Spatially-varying intercept and covariate effects
w <- dat$w
# Simple plot of the relative occurrence probability across space.
# Dark points indicate high occurrence.
plot(coords, type = "n", xlab = "", ylab = "", asp = TRUE,
main = "Simulated Occurrence", bty = 'n')
points(coords, pch=15, cex = 2.1, col = rgb(0,0,0,(psi-min(psi))/diff(range(psi))))
Lastly, we package the data up in a list for use in spOccupancy model fitting functions. For SVC GLMs, this requires creating a list that consists of the detection-nondetection data (y), covariates (covs), the binomial weights aka the number of replicates (weights), and coordinates (coords). Note that the covariates here are stored in a single matrix object covs rather than split apart into the occurrence (occ.covs) and detection (det.covs) covariates as we do for occupancy models. As we did for the SVC occupancy model, we will fit the model to 75% of the data points (1200 locations) and subsequently predict at the remaining 400 values to show the predictive ability of the model.
# Subset data for prediction.
# Split into fitting and prediction data set
pred.indx <- sample(1:J, round(J * .25), replace = FALSE)
y.fit <- y[-pred.indx]
y.pred <- y[pred.indx]
X.fit <- X[-pred.indx, ]
X.pred <- X[pred.indx, ]
coords.fit <- coords[-pred.indx, ]
coords.pred <- coords[pred.indx, ]
psi.fit <- psi[-pred.indx]
psi.pred <- psi[pred.indx]
w.fit <- w[-pred.indx, ]
w.pred <- w[pred.indx, ]
weights.fit <- weights[-pred.indx]
weights.pred <- weights[pred.indx]
# Package all data into a list
# Covariates
covs <- X.fit[, 2, drop = FALSE]
colnames(covs) <- c('cov.1')
# Package into a list for spOccupancy
data.list <- list(y = y.fit,
covs = covs,
coords = coords.fit,
weights = weights.fit)
# Take a look at the data structure.
str(data.list)List of 4
$ y : int [1:1200] 0 1 0 0 0 0 1 1 0 1 ...
$ covs : num [1:1200, 1] -0.178 0.791 0.733 -1.758 1.847 ...
..- attr(*, "dimnames")=List of 2
.. ..$ : NULL
.. ..$ : chr "cov.1"
$ coords : num [1:1200, 1:2] 0 0.0256 0.0513 0.1026 0.1282 ...
..- attr(*, "dimnames")=List of 2
.. ..$ : NULL
.. ..$ : chr [1:2] "Var1" "Var2"
$ weights: num [1:1200] 1 1 1 1 1 1 1 1 1 1 ...Fitting spatially-varying coefficients generalized linear models with svcPGBinom()
The function svcPGBinom() fits single-season SVC GLMs. svcPGBinom() has the following arguments, which are nearly identical to those we saw with svcPGOcc():
svcPGBinom(formula, data, inits, priors, tuning, svc.cols = 1,
cov.model = "exponential", NNGP = TRUE,
n.neighbors = 15, search.type = "cb", n.batch,
batch.length, accept.rate = 0.43,
n.omp.threads = 1, verbose = TRUE, n.report = 100,
n.burn = round(.10 * n.batch * batch.length),
n.thin = 1, n.chains = 1, k.fold, k.fold.threads = 1,
k.fold.seed = 100, k.fold.only = FALSE, ...)The only difference between the arguments of svcPGBinom() and svcPGOcc() is that now we only specify a single formula, formula, since we do not separately model detection probability from occurrence probability. Below we specify the formula, including the single covariate we simulated the data with.
formula <- ~ cov.1As before, the svc.cols argument specifies the covariates whose effects are estimated as SVCs. Here we set svc.cols to include a spatially-varying intercept and spatially-varying effect of the covariate. We also specify the cov.model argument to indicate we will use an exponential correlation function.
svc.cols <- c(1, 2)
cov.model <- 'exponential'We next specify the initial values, which is exactly analogous to svcPGOcc(), but we don’t have any initial values for the detection parameters or the latent occupancy values, since we don’t estimate those in an SVC GLM.
dist.data <- dist(data.list$coords)
inits.list <- list(beta = 0, sigma.sq = 0.5,
phi = 3 / mean(dist.data),
w = matrix(0, length(svc.cols), length(data.list$y)))We next specify the priors to use for all parameters in the model, which we set to be the same as those we used for svcPGOcc().
priors.list <- list(beta.normal = list(mean = 0, var = 2.72),
sigma.sq.ig = list(a = 2, b = 0.5),
phi.unif = list(a = 3 / max(dist.data),
b = 3 / min(dist.data)))Finally, we specify the MCMC criteria exactly as we saw previously and then we are all set to run the model.
batch.length <- 25
n.batch <- 800
n.burn <- 10000
n.thin <- 10
n.chains <- 1
tuning.list <- list(phi = 0.2, sigma.sq = 0.2)
n.omp.threads <- 1
verbose <- TRUE
n.report <- 100 # Report progress at every 100th batch.
# Approx. run time: 3.7 min
out.svc <- svcPGBinom(formula = formula,
data = data.list,
inits = inits.list,
n.batch = n.batch,
batch.length = batch.length,
priors = priors.list,
svc.cols = svc.cols,
cov.model = cov.model,
NNGP = TRUE,
n.neighbors = 5,
tuning = tuning.list,
n.report = n.report,
n.burn = n.burn,
n.thin = n.thin,
n.chains = n.chains)----------------------------------------
Preparing to run the model
----------------------------------------
----------------------------------------
Building the neighbor list
----------------------------------------
----------------------------------------
Building the neighbors of neighbors list
----------------------------------------
----------------------------------------
Model description
----------------------------------------
Spatial NNGP Binomial model with Polya-Gamma latent
variable fit with 1200 sites.
Samples per chain: 20000 (800 batches of length 25)
Burn-in: 10000
Thinning Rate: 10
Number of Chains: 1
Total Posterior Samples: 1000
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: 100 of 800, 12.50%
Coefficient Parameter Acceptance Tuning
1 phi 36.0 0.18648
2 phi 24.0 0.22326
-------------------------------------------------
Batch: 200 of 800, 25.00%
Coefficient Parameter Acceptance Tuning
1 phi 40.0 0.16539
2 phi 44.0 0.16873
-------------------------------------------------
Batch: 300 of 800, 37.50%
Coefficient Parameter Acceptance Tuning
1 phi 52.0 0.17562
2 phi 44.0 0.17917
-------------------------------------------------
Batch: 400 of 800, 50.00%
Coefficient Parameter Acceptance Tuning
1 phi 24.0 0.22777
2 phi 28.0 0.17214
-------------------------------------------------
Batch: 500 of 800, 62.50%
Coefficient Parameter Acceptance Tuning
1 phi 76.0 0.20201
2 phi 32.0 0.17214
-------------------------------------------------
Batch: 600 of 800, 75.00%
Coefficient Parameter Acceptance Tuning
1 phi 44.0 0.15891
2 phi 68.0 0.20609
-------------------------------------------------
Batch: 700 of 800, 87.50%
Coefficient Parameter Acceptance Tuning
1 phi 48.0 0.19025
2 phi 36.0 0.19409
-------------------------------------------------
Batch: 800 of 800, 100.00%
# Compare to values used to generate the data
summary(out.svc)
Call:
svcPGBinom(formula = formula, data = data.list, inits = inits.list,
priors = priors.list, tuning = tuning.list, svc.cols = svc.cols,
cov.model = cov.model, NNGP = TRUE, n.neighbors = 5, n.batch = n.batch,
batch.length = batch.length, n.report = n.report, n.burn = n.burn,
n.thin = n.thin, n.chains = n.chains)
Samples per Chain: 20000
Burn-in: 10000
Thinning Rate: 10
Number of Chains: 1
Total Posterior Samples: 1000
Run Time (min): 7.0394
Occurrence (logit scale):
Mean SD 2.5% 50% 97.5% Rhat ESS
(Intercept) -0.2728 0.1473 -0.5770 -0.2596 0.0057 NA 34
cov.1 0.0831 0.1470 -0.2261 0.0913 0.3611 NA 22
Spatial Covariance:
Mean SD 2.5% 50% 97.5% Rhat ESS
sigma.sq-(Intercept) 0.3826 0.1631 0.1793 0.3473 0.8226 NA 31
sigma.sq-cov.1 0.2842 0.1788 0.0846 0.2438 0.8155 NA 14
phi-(Intercept) 6.8393 2.6174 2.8451 6.4441 13.8589 NA 27
phi-cov.1 14.9416 15.3815 2.7489 9.4919 66.0294 NA 10
# True values
beta[1] -0.5 -0.2
sigma.sq[1] 1.0 0.5
phi[1] 3.750000 4.285714Let’s next extract the full SVC values using the getSVCSamples() function and then compare the estimated values to those used to generate the model. This time, we won’t make a map of the predicted values, but rather will just look at a scatter plot of fitted values vs. true values.
# Intercept ---------------------------------------------------------------
svc.samples <- getSVCSamples(out.svc)
int.quants <- apply(svc.samples[["(Intercept)"]], 2,
quantile, probs = c(0.025, 0.5, 0.975))
svc.true.fit <- beta + w.fit
plot(svc.true.fit[, 1], int.quants[2, ], pch = 19,
ylim = c(min(int.quants[1, ]), max(int.quants[3, ])),
xlab = "True", ylab = "Fit", main = 'Intercept')
abline(0, 1)
arrows(svc.true.fit[, 1], int.quants[2, ], svc.true.fit[, 1],
col = 'gray', int.quants[1, ], length = 0.02, angle = 90)
arrows(svc.true.fit[, 1], int.quants[1, ], svc.true.fit[, 1],
col = 'gray', int.quants[3, ], length = 0.02, angle = 90)
points(svc.true.fit[, 1], int.quants[2, ], pch = 19)
# Covariate ---------------------------
cov.quants <- apply(svc.samples[["cov.1"]], 2, quantile,
probs = c(0.025, 0.5, 0.975))
svc.true.fit <- beta + w.fit
plot(svc.true.fit[, 2], cov.quants[2, ], pch = 19,
ylim = c(min(cov.quants[1, ]), max(cov.quants[3, ])),
xlab = "True", ylab = "Fit", main = 'Covariate')
abline(0, 1)
arrows(svc.true.fit[, 2], cov.quants[2, ], svc.true.fit[, 2],
col = 'gray', cov.quants[1, ], length = 0.02, angle = 90)
arrows(svc.true.fit[, 2], cov.quants[1, ], svc.true.fit[, 2],
col = 'gray', cov.quants[3, ], length = 0.02, angle = 90)
points(svc.true.fit[, 2], cov.quants[2, ], pch = 19)
Here we see our estimated values capture the patterns in the simulated values. Notice that the estimated values do not fall exactly on the one to one line, and rather the estimated values are all seemingly closer to zero than the true values, in particular for the sites that have a large magnitude of the effect. This is a common phenomenon in spatial models that use binary data in a GLM framework as a result of the lack of information in binary data as well as the link function (e.g., logit) that transforms binary data to the real scale. Generally, as the number of spatial locations increase and the spatial range increases (i.e., the spatial decay parameter \(\phi\) gets smaller), we will tend to see the values be more closely approximated by the model. Regardless, we generally see good coverage of the 95% credible intervals for the true values (ideally we want 95% of the gray lines to overlap the black one to one line), and the estimated means capture the pattern in the effects across space.
Model selection using WAIC
We can do model selection using WAIC as we saw with svcPGOcc(). Below we fit a spatial GLM that assumes the covariate effect is stationary across the study area. We then compare the two models with WAIC.
out.glm <- svcPGBinom(formula = formula,
data = data.list,
n.batch = n.batch,
batch.length = batch.length,
svc.cols = 1,
cov.model = cov.model,
NNGP = TRUE,
n.neighbors = 5,
tuning = tuning.list,
n.report = n.report,
n.burn = n.burn,
n.thin = n.thin,
n.chains = n.chains)----------------------------------------
Preparing to run the model
----------------------------------------No prior specified for beta.normal.
Setting prior mean to 0 and prior variance to 2.72No prior specified for phi.unif.
Setting uniform bounds based on the range of observed spatial coordinates.No prior specified for sigma.sq.
Using an inverse-Gamma prior with the shape parameter set to 2 and scale parameter to 1.beta is not specified in initial values.
Setting initial values to random values from the prior distributionphi is not specified in initial values.
Setting initial value to random values from the prior distributionsigma.sq is not specified in initial values.
Setting initial values to random values from the prior distributionw is not specified in initial values.
Setting initial value to 0----------------------------------------
Building the neighbor list
----------------------------------------
----------------------------------------
Building the neighbors of neighbors list
----------------------------------------
----------------------------------------
Model description
----------------------------------------
Spatial NNGP Binomial model with Polya-Gamma latent
variable fit with 1200 sites.
Samples per chain: 20000 (800 batches of length 25)
Burn-in: 10000
Thinning Rate: 10
Number of Chains: 1
Total Posterior Samples: 1000
Number of spatially-varying coefficients: 1
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: 100 of 800, 12.50%
Coefficient Parameter Acceptance Tuning
1 phi 64.0 0.19801
-------------------------------------------------
Batch: 200 of 800, 25.00%
Coefficient Parameter Acceptance Tuning
1 phi 24.0 0.17917
-------------------------------------------------
Batch: 300 of 800, 37.50%
Coefficient Parameter Acceptance Tuning
1 phi 52.0 0.15268
-------------------------------------------------
Batch: 400 of 800, 50.00%
Coefficient Parameter Acceptance Tuning
1 phi 28.0 0.18279
-------------------------------------------------
Batch: 500 of 800, 62.50%
Coefficient Parameter Acceptance Tuning
1 phi 52.0 0.16539
-------------------------------------------------
Batch: 600 of 800, 75.00%
Coefficient Parameter Acceptance Tuning
1 phi 44.0 0.17214
-------------------------------------------------
Batch: 700 of 800, 87.50%
Coefficient Parameter Acceptance Tuning
1 phi 32.0 0.19025
-------------------------------------------------
Batch: 800 of 800, 100.00%
# SVC model
waicOcc(out.svc) elpd pD WAIC
-715.54045 80.42124 1591.92337
# Spatially-varying intercept model
waicOcc(out.glm) elpd pD WAIC
-735.14310 67.59059 1605.46739 As expected, the SVC GLM outperforms the spatial GLM that assumes a constant covariate effect.
Prediction
Finally, we can use the predict() function just as we saw with svcPGOcc() to predict relative occurrence and the effects of the covariates at new (and old) locations.
# Take a look at X.pred, the design matrix for the prediction locations
head(X.pred) [,1] [,2]
[1,] 1 -1.3248935
[2,] 1 0.4733202
[3,] 1 -0.4613708
[4,] 1 0.6230023
[5,] 1 0.2933181
[6,] 1 -0.5235239
# Predict occupancy at the 400 new sites
out.pred <- predict(out.svc, X.pred, coords.pred, weights.pred)----------------------------------------
Prediction description
----------------------------------------
NNGP Occupancy model with Polya-Gamma latent
variable fit with 1200 observations.
Number of covariates 2 (including intercept if specified).
Number of spatially-varying covariates 2 (including intercept if specified).
Using the exponential spatial correlation model.
Using 5 nearest neighbors.
Number of MCMC samples 1000.
Predicting at 400 non-sampled locations.
Source compiled with OpenMP support and model fit using 1 threads.
-------------------------------------------------
Predicting
-------------------------------------------------
Location: 100 of 400, 25.00%
Location: 200 of 400, 50.00%
Location: 300 of 400, 75.00%
Location: 400 of 400, 100.00%
Generating latent occupancy state
# Use the getSVCSamples() function to extract the SVC values
# at the prediction locations
svc.pred.samples <- getSVCSamples(out.svc, pred.object = out.pred)
# True covariate effect values at new locations
cov.effect.pred <- beta[2] + w.pred[, 2]
# Get median and 95% CIs of the SVC for the covariate effect
cov.pred.quants <- apply(svc.pred.samples[["cov.1"]], 2,
quantile, probs = c(0.025, 0.5, 0.975))
plot(cov.effect.pred, cov.pred.quants[2, ], pch = 19,
ylim = c(min(cov.pred.quants[1, ]), max(cov.pred.quants[3, ])),
xlab = "True", ylab = "Fit", main = 'Covariate')
abline(0, 1)
arrows(cov.effect.pred, cov.pred.quants[2, ], cov.effect.pred,
col = 'gray', cov.pred.quants[1, ],
length = 0.02, angle = 90)
arrows(cov.effect.pred, cov.pred.quants[1, ], cov.effect.pred,
col = 'gray', cov.pred.quants[3, ],
length = 0.02, angle = 90)
points(cov.effect.pred, cov.pred.quants[2, ], pch = 19)
As with the fitted values, we see the predicted values seem to capture the pattern in the true estimates.
Spatially-varying coefficients multi-season occupancy model
Consider the case where detection-nondetection data are collected across multiple seasons (e.g., years, breeding periods) during which the true occupancy status can change. Such spatio-temporal data can be used to understand occurrence trends over time (e.g., Isaac et al. (2014)), as well as the environmental variables that drive spatio-temporal shifts in species distributions (e.g., Rushing et al. (2020)). As an extension to the SVC occupancy model, here we present a SVC multi-season occupancy model that estimates spatially-varying effects (as described previously) of covariates and accounts for temporal autocorrelation using temporal random effects. This model is an extension of the multi-season spatial (stPGOcc()) and non-spatial (tPGOcc()) occupancy models introduced in spOccupancy v0.4.0. See the vignette on multiseason models for additional background information on modeling spatio-temporal patterns in occupancy.
Basic model description
Let \(y_{k, t}(\boldsymbol{s}_j)\) denote the detection or nondetection of the species of interest at site \(j\) during replicate survey \(k\) during time period \(t\), with \(t = 1, \dots, T\). We model \(y_{k, t}(\boldsymbol{s}_j)\) conditional on the true occupancy status, \(z_t(\boldsymbol{s}_j)\) of the species at site \(j\) during time \(t\) according to
\[\begin{equation}\label{zTime} z_{t}(\boldsymbol{s}_j) \sim \text{Bernoulli}(\psi_t(\boldsymbol{s}_j)), \end{equation}\]
where \(\psi_t(\boldsymbol{s}_j)\) is the occurrence probability at site \(j\) during primary time period \(t\). We model \(\psi_{t}(\boldsymbol{s}_j)\) following
\[\begin{equation}\label{psiTime} \text{logit}(\psi_t(\boldsymbol{s}_j)) = \textbf{x}_t(\boldsymbol{s}_j)\boldsymbol{\beta} + \tilde{\textbf{x}}_t(\boldsymbol{s}_j)\textbf{w}(\boldsymbol{s}_j) + \eta_t, \end{equation}\]
where \(\eta_t\) is a random temporal effect and the other parameters are defined as before. Note that the covariates can now vary across space and/or time period. Because we assume the non-spatial (\(\boldsymbol{\beta})\) and spatially-varying (\(\textbf{w}(\boldsymbol{s}_j)\)) coefficients are constant over the \(T\) primary time periods, they represent the average covariate effects across the temporal extent of the data. We model \(\eta_t\) as either an unstructured random effect (i.e., \(\eta_t \sim \text{Normal}(0, \sigma^2_{T})\)) or using a first-order autoregressive (i.e., AR(1)) covariance structure in which we estimate a temporal variance parameter, \(\sigma^2_T\), and a temporal correlation parameter, \(\rho\).
The data \(y_{k, t}(\boldsymbol{s}_j)\) are modeled conditional on the true occupancy status \(z_t(\boldsymbol{s}_j)\) analogous to the SVC occupancy model, with detection probability now allowed to vary across site, replicates, and/or primary time periods.
Simulating data with simTOcc()
The function simTOcc() simulates single-species multi-season detection-nondetection data. simTOcc() has the following arguments, which again closely resemble all other data simulation functions in spOccupancy.
simTOcc(J.x, J.y, n.time, n.rep, beta, alpha, sp.only = 0, trend = TRUE,
psi.RE = list(), p.RE = list(), sp = FALSE, svc.cols = 1, cov.model,
sigma.sq, phi, nu, ar1 = FALSE, rho, sigma.sq.t, x.positive = FALSE,...)J.x and J.y indicate the number of spatial locations to simulate data along a horizontal and vertical axis, respectively, such that J.x * J.y is the total number of sites (i.e., J). n.time indicates the number of time periods over which data are collected. n.rep is a numeric matrix with rows corresponding to sites and columns corresponding to time periods where each element corresponds to the number of replicate surveys performed within a given time period at a given site (denoted as K in the previous model description). Note that throughout we will refer to “time period” as the additional temporal replication in multi-season models (often called primary time period) and will use the term “replicate” to refer to the multiple visits obtained at a given site during a given time period (often called secondary time periods). beta and alpha are numeric vectors containing the intercept and any regression coefficient parameters for the occurrence and detection portions of the occupancy model, respectively. sp.only is a numeric vector specifying which occurrence covariates should only vary over space and not over time. trend is a logical value indicating whether or not a trend parameter is included in the model (the trend is assumed to be the second covariate listed in beta). psi.RE and p.RE are lists that are used to specify random intercepts on occurrence and detection, respectively, which are analogous to previous data simulation functions. sp is a logical value indicating whether to simulate data with at least one spatial Gaussian process for either the intercept or some of the occupancy covariate effects. svc.cols is a numeric vector indicating which of the covariates (including the intercept) are generated with spatially-varying effects. cov.model specifies the covariance function used to model the spatial dependence structure. Finally, sigma.sq is the spatial variance parameter, phi is the spatial range parameter, and nu is the spatial smoothness parameter (only applicable when cov.model = 'matern'). Note that sigma.sq, phi, and nu should have the same length as the number of spatially-varying effects specified in svc.cols. Lastly, the x.positive argument indicates whether or not the occupancy covariates should be simulated to be only positive from a Uniform(0, 1) distribution (TRUE) or both positive and negative and simulated from a Normal(0, 1) distribution (FALSE).
For this simulation, we will attempt to estimate a spatially-varying trend parameter. Such an analysis can provide useful information for monitoring programs by allowing for prioritization of locations where occupancy trends are declining. Additionally, it can often help generate hypotheses as to what factors are driving the heterogenous trends across space. Below we simulate data across 400 sites and ten time periods with anywhere between 1-4 replicates at a given time-period, a spatially-varying intercept, a spatially-varying trend on occupancy, and a single covariate effect on detection. Note that we vary the number of time periods each site is sampled. We assume both the occupancy intercept and the trend vary across space, so we set svc.cols = c(1, 2). We use an exponential correlation function. We do not include any unstructured random effects on occurrence or detection.
set.seed(10)
J.x <- 20
J.y <- 20
# Total number of sites
(J <- J.x * J.y)[1] 400
# Total number of time periods at each site
n.time <- sample(10, J, replace = TRUE)
n.time.max <- max(n.time)
# Number of replicates at each site/time period combination
n.rep <- matrix(NA, J, max(n.time))
for (j in 1:J) {
n.rep[j, 1:n.time[j]] <- sample(4, n.time[j], replace = TRUE)
}
sp.only <- 0
# Simulate a trend parameter
trend <- TRUE
# Intercept and trend on occurrence
beta <- c(-0.5, -0.2)
# Intercept and covariate effect on detection
alpha <- c(0.9, -0.3)
# No unstructured random intercept on occurrence
psi.RE <- list()
# No unstructured random intercept on detection
p.RE <- list()
# Spatial range for intercept and trend
phi <- c(3 / .8, 3 / .7)
# Spatial variance for intercept and trend
sigma.sq <- c(1, 0.5)
# Simulate data with an AR1 temporal process
ar1 <- TRUE
# Temporal variance of AR1 random effect
sigma.sq.t <- 0.5
# Temporal correlation of AR1 random effect
rho <- 0.25
# Simulate the occupancy covariate from a Normal(0, 1) distribution
x.positive <- FALSE
# Spatially-varying coefficient columns
svc.cols <- c(1, 2)
# Simulate the data
dat <- simTOcc(J.x = J.x, J.y = J.y, n.time = n.time, n.rep = n.rep, beta = beta,
alpha = alpha, sp.only = sp.only, trend = trend,
psi.RE = psi.RE, p.RE = p.RE, sp = TRUE,
sigma.sq = sigma.sq, phi = phi,
cov.model = 'exponential', svc.cols = svc.cols,
x.positive = x.positive, ar1 = ar1, sigma.sq.t = sigma.sq.t,
rho = rho)First let’s take a look at the true occupancy values averaged across all spatial locations to see if there appears to be an overall trend over time
str(dat)List of 13
$ X : num [1:400, 1:10, 1:2] 1 1 1 1 1 1 1 1 1 1 ...
$ X.p : num [1:400, 1:10, 1:4, 1:2] 1 1 1 1 1 1 1 1 1 1 ...
$ coords : num [1:400, 1:2] 0 0.0526 0.1053 0.1579 0.2105 ...
..- attr(*, "dimnames")=List of 2
.. ..$ : NULL
.. ..$ : chr [1:2] "Var1" "Var2"
$ psi : num [1:400, 1:10] 0.535 0.395 0.229 0.465 0.91 ...
$ z : int [1:400, 1:10] 1 1 1 1 1 1 0 1 0 0 ...
$ p : num [1:400, 1:10, 1:4] 0.658 0.722 0.568 0.741 0.664 ...
$ y : int [1:400, 1:10, 1:4] 1 1 0 0 1 1 0 1 0 0 ...
$ w : num [1:400, 1:2] -0.188959 -0.000185 -1.212923 -0.730247 0.348285 ...
$ X.p.re : logi NA
$ X.re : logi NA
$ alpha.star: logi NA
$ beta.star : logi NA
$ eta : num [1:10, 1] 0.194 0.425 0.587 0.161 1.795 ...
# Plot average occupancy probability
plot(apply(dat$psi, 2, mean), pch = 19, xlab = 'Time Period',
ylab = 'Average Occupancy Probability',
cex = 1.5, frame = FALSE, ylim = c(0, 1))
The simulated data object consists of the same objects as we saw previously with simOcc(). However, because of the additional time period dimension, all of the objects (except for spatial coordinates and spatially-varying effects) now have an additional dimension of time period. For example, the occurrence design matrix X is now a three-dimensional array with dimensions corresponding to site, time period, and covariate. We see on average occupancy is decreasing over the simulated ten year period, which makes sense given the average value we used to simulate the data (-0.2).
Below we extract the data we need to fit the model and create a plot of the spatially-varying trend across the study region
# Detection-nondetection data
y <- dat$y
# Occurrence design matrix for fixed effects
X <- dat$X
# Detection design matrix for fixed effets
X.p <- dat$X.p
# Occurrence values
psi <- dat$psi
# Spatial coordinates
coords <- dat$coords
# Spatially-varying intercept and covariate effects
w <- dat$w
cov.effect <- beta[2] + w[, 2]
plot.dat <- data.frame(x = coords[, 1],
y = coords[, 2],
cov.effect = cov.effect)
ggplot(plot.dat, aes(x = x, y = y, fill = cov.effect)) +
geom_raster() +
scale_fill_gradient2(midpoint = 0, low = '#B2182B', mid = 'white',
high = '#2166AC', na.value = NA) +
theme_bw()
We see heterogenous trends across the study area with some locations showing a large decrease in occurrence (red) and others showing a large increase (blue).
We finally package up the data into the data list format for multi-season occupancy models. As before, we will fit the model with 75% of the spatial locations and then predict at the remaining 25% of the locations. With multi-season occupancy models, the occurrence covariates can now vary across both space and time, and so we specify occ.covs as a list rather than as a matrix/data frame as we saw with svcPGOcc(). See the multi-season occupancy model vignette for additional details.
# Subset data for prediction.
# Split into fitting and prediction data set
pred.indx <- sample(1:J, round(J * .25), replace = FALSE)
y.fit <- y[-pred.indx, , ]
y.pred <- y[pred.indx, , ]
X.fit <- X[-pred.indx, , ]
X.pred <- X[pred.indx, , ]
X.p.fit <- X.p[-pred.indx, , , ]
X.p.pred <- X.p[pred.indx, , , ]
coords.fit <- coords[-pred.indx, ]
coords.pred <- coords[pred.indx, ]
psi.fit <- psi[-pred.indx, ]
psi.pred <- psi[pred.indx, ]
w.fit <- w[-pred.indx, ]
w.pred <- w[pred.indx, ]
# Package all data into a list
# Occurrence covariates
occ.covs <- list(trend = X.fit[, , 2])
# Detection covariates
det.covs <- list(det.cov.1 = X.p.fit[, , , 2])
# Package into a list for spOccupancy
data.list <- list(y = y.fit,
occ.covs = occ.covs,
det.covs = det.covs,
coords = coords.fit)
# Take a look at the data structure.
str(data.list)List of 4
$ y : int [1:300, 1:10, 1:4] 1 1 0 1 0 1 0 0 1 0 ...
$ occ.covs:List of 1
..$ trend: num [1:300, 1:10] -1.57 -1.57 -1.57 -1.57 -1.57 ...
$ det.covs:List of 1
..$ det.cov.1: num [1:300, 1:10, 1:4] 0.818 -0.177 -0.501 0.353 0.701 ...
$ coords : num [1:300, 1:2] 0 0.0526 0.1579 0.2632 0.3158 ...
..- attr(*, "dimnames")=List of 2
.. ..$ : NULL
.. ..$ : chr [1:2] "Var1" "Var2"Fitting spatially-varying coefficients multi-season occupancy models with svcPGOcc()
The function svcTPGOcc() fits spatially-varyig coefficients multi-season occupancy models. It’s arguments are exactly identical to svcPGOcc(), with the only addition being the ar1 argument.
svcTPGOcc(occ.formula, det.formula, data, inits, priors,
tuning, svc.cols = 1, cov.model = 'exponential', NNGP = TRUE,
n.neighbors = 15, search.type = 'cb', n.batch,
batch.length, accept.rate = 0.43, n.omp.threads = 1,
verbose = TRUE, ar1 = FALSE, n.report = 100,
n.burn = round(.10 * n.batch * batch.length),
n.thin = 1, n.chains = 1, k.fold, k.fold.threads = 1,
k.fold.seed = 100, k.fold.only = FALSE, ...)Given the similarity in the syntax for svcTPGOcc() and svcPGOcc(), we won’t go into all that much detail on them here. For additional details on the input parameters related to the temporal parameters, see the multi-season occupancy model vignette.
The occ.formula and det.formula are specified in the same manner as we saw previously. The data argument contains the list of the formatted data, which we discussed in the preceding section for multi-season models.
data <- data.list
str(data)List of 4
$ y : int [1:300, 1:10, 1:4] 1 1 0 1 0 1 0 0 1 0 ...
$ occ.covs:List of 1
..$ trend: num [1:300, 1:10] -1.57 -1.57 -1.57 -1.57 -1.57 ...
$ det.covs:List of 1
..$ det.cov.1: num [1:300, 1:10, 1:4] 0.818 -0.177 -0.501 0.353 0.701 ...
$ coords : num [1:300, 1:2] 0 0.0526 0.1579 0.2632 0.3158 ...
..- attr(*, "dimnames")=List of 2
.. ..$ : NULL
.. ..$ : chr [1:2] "Var1" "Var2"
occ.formula <- ~ trend
det.formula <- ~ det.cov.1Next we specify the initial values for all model parameters. We will fit the model with an AR(1) temporal covariance structure, which includes two additional parameters: sigma.sq.t (a temporal variance parameter) and rho (a temporal correlation parameter). We specifiy initial values for z to 1 if the species was ever observed at the given site/time period combination and 0 otherwise.
z.inits <- apply(data$y, c(1, 2), function(a) as.numeric(sum(a, na.rm = TRUE) > 0))
inits <- list(beta = 0, alpha = 0, z = z.inits,
sigma.sq = 1, phi = 3 / 0.5,
sigma.sq.t = 1.5, rho = 0.2)We specify priors in the priors argument just as we saw with svcPGOcc(). We use an inverse-Gamma prior for the temporal variance sigma.sq.t and a uniform prior for the temporal correlation parameter rho.
priors <- list(beta.normal = list(mean = 0, var = 2.72),
alpha.normal = list(mean = 0, var = 2.72),
sigma.sq.t.ig = c(2, 0.5),
rho.unif = c(-1, 1),
sigma.sq.ig = list(a = 2, b = 1),
phi.unif = list(a = 3 / 1, b = 3 / .1))Next we specify parameters associated with the spatial random effects. In particular, we set cov.model = 'exponential' to use an exponential spatial correlation function, n.neighbors = 5 to use an NNGP with 5 nearest neighbors, and svc.cols = c(1, 2) to specify a spatially-varying intercept and trend. We also specify the ar1 argument to indicate we will use an AR(1) temporal covariance structure.
cov.model <- 'exponential'
svc.cols <- c(1, 2)
n.neighbors <- 5
ar1 <- TRUEFinally, we set the number of MCMC batches, batch length, the amount of burn-in, and our thinning rate.
n.batch <- 600
batch.length <- 25
# Total number of samples
n.batch * batch.length[1] 15000
n.burn <- 10000
n.thin <- 20We now run the model with svcTPGOcc() and take a look at a summary of the results using summary().
# Approx. run time: ~ 1.1 min
out.svc.trend <- svcTPGOcc(occ.formula = occ.formula,
det.formula = det.formula,
data = data,
inits = inits,
priors = priors,
cov.model = cov.model,
svc.cols = svc.cols,
n.neighbors = n.neighbors,
n.batch = n.batch,
batch.length = batch.length,
verbose = TRUE,
ar1 = ar1,
n.report = 200,
n.burn = n.burn,
n.thin = n.thin,
n.chains = 1)----------------------------------------
Preparing the data
----------------------------------------
----------------------------------------
Building the neighbor list
----------------------------------------
----------------------------------------
Building the neighbors of neighbors list
----------------------------------------
----------------------------------------
Model description
----------------------------------------
Spatial NNGP Trend Occupancy Model with Polya-Gamma latent
variable fit with 300 sites and 10 years.
Samples per chain: 15000 (600 batches of length 25)
Burn-in: 10000
Thinning Rate: 20
Number of Chains: 1
Total Posterior Samples: 250
Number of spatially-varying coefficients: 2
Using the exponential spatial correlation model.
Using 5 nearest neighbors.
Using an AR(1) temporal autocorrelation matrix.
Source compiled with OpenMP support and model fit using 1 thread(s).
Adaptive Metropolis with target acceptance rate: 43.0
----------------------------------------
Chain 1
----------------------------------------
Sampling ...
Batch: 200 of 600, 33.33%
Coefficient Parameter Acceptance Tuning
1 phi 40.0 0.66365
2 phi 48.0 0.57695
1 rho 40.0 1.11628
-------------------------------------------------
Batch: 400 of 600, 66.67%
Coefficient Parameter Acceptance Tuning
1 phi 40.0 0.63763
2 phi 44.0 0.53259
1 rho 48.0 1.13883
-------------------------------------------------
Batch: 600 of 600, 100.00%
summary(out.svc.trend)
Call:
svcTPGOcc(occ.formula = occ.formula, det.formula = det.formula,
data = data, inits = inits, priors = priors, svc.cols = svc.cols,
cov.model = cov.model, n.neighbors = n.neighbors, n.batch = n.batch,
batch.length = batch.length, verbose = TRUE, ar1 = ar1, n.report = 200,
n.burn = n.burn, n.thin = n.thin, n.chains = 1)
Samples per Chain: 15000
Burn-in: 10000
Thinning Rate: 20
Number of Chains: 1
Total Posterior Samples: 250
Run Time (min): 1.9582
Occurrence (logit scale):
Mean SD 2.5% 50% 97.5% Rhat ESS
(Intercept) -0.3418 0.3825 -1.0919 -0.3266 0.4991 NA 53
trend -0.7806 0.3669 -1.3587 -0.8054 0.0165 NA 24
Detection (logit scale):
Mean SD 2.5% 50% 97.5% Rhat ESS
(Intercept) 0.9469 0.0564 0.8448 0.9492 1.0501 NA 250
det.cov.1 -0.3115 0.0515 -0.4147 -0.3087 -0.1962 NA 250
Spatio-temporal Covariance:
Mean SD 2.5% 50% 97.5% Rhat ESS
sigma.sq-(Intercept) 0.6124 0.2408 0.2654 0.5746 1.1159 NA 86
sigma.sq-trend 0.7903 0.3047 0.3939 0.7391 1.5609 NA 75
phi-(Intercept) 5.8699 2.2306 3.2610 5.2182 11.8320 NA 45
phi-trend 7.9204 3.1375 3.5599 7.5867 15.1393 NA 52
sigma.sq.t 0.6199 0.3944 0.2203 0.5300 1.3731 NA 250
rho 0.0847 0.2361 -0.3542 0.0803 0.5332 NA 250
# Compare estimates to true values
beta[1] -0.5 -0.2
alpha[1] 0.9 -0.3
sigma.sq[1] 1.0 0.5
phi[1] 3.750000 4.285714
rho[1] 0.25
sigma.sq.t[1] 0.5Next, we extract the full SVC values using the getSVCSamples() function and then compare the estimated values to those used to generate the model.
# Intercept ---------------------------------------------------------------
svc.samples <- getSVCSamples(out.svc.trend)
int.quants <- apply(svc.samples[["(Intercept)"]], 2, quantile,
probs = c(0.025, 0.5, 0.975))
svc.true.fit <- beta + w.fit
plot(svc.true.fit[, 1], int.quants[2, ], pch = 19,
ylim = c(min(int.quants[1, ]), max(int.quants[3, ])),
xlab = "True", ylab = "Fit", main = 'Intercept')
abline(0, 1)
arrows(svc.true.fit[, 1], int.quants[2, ], svc.true.fit[, 1], col = 'gray',
int.quants[1, ], length = 0.02, angle = 90)
arrows(svc.true.fit[, 1], int.quants[1, ], svc.true.fit[, 1], col = 'gray',
int.quants[3, ], length = 0.02, angle = 90)
points(svc.true.fit[, 1], int.quants[2, ], pch = 19)
# Trend -------------------------------------------------------------------
trend.quants <- apply(svc.samples[["trend"]], 2, quantile,
probs = c(0.025, 0.5, 0.975))
plot(svc.true.fit[, 2], trend.quants[2, ], pch = 19,
ylim = c(min(trend.quants[1, ]), max(trend.quants[3, ])),
xlab = "True", ylab = "Fit", main = 'Spatially-Varying Trends')
abline(0, 1)
arrows(svc.true.fit[, 2], trend.quants[2, ], svc.true.fit[, 2], col = 'gray',
trend.quants[1, ], length = 0.02, angle = 90)
arrows(svc.true.fit[, 2], trend.quants[1, ], svc.true.fit[, 2], col = 'gray',
trend.quants[3, ], length = 0.02, angle = 90)
points(svc.true.fit[, 2], trend.quants[2, ], pch = 19)
We see our model does a good job of recovering the true parameter values for both the spatially-varying intercept and spatially-varying trend.
Posterior Predictive Checks
Next we run a posterior predictive check using a Freeman-Tukey statistic and grouping by site.
ppc.out <- ppcOcc(out.svc.trend, fit.stat = 'freeman-tukey', group = 1)Currently on time period 1 out of 10Currently on time period 2 out of 10Currently on time period 3 out of 10Currently on time period 4 out of 10Currently on time period 5 out of 10Currently on time period 6 out of 10Currently on time period 7 out of 10Currently on time period 8 out of 10Currently on time period 9 out of 10Currently on time period 10 out of 10
summary(ppc.out)
Call:
ppcOcc(object = out.svc.trend, fit.stat = "freeman-tukey", group = 1)
Samples per Chain: 15000
Burn-in: 10000
Thinning Rate: 20
Number of Chains: 1
Total Posterior Samples: 250
----------------------------------------
All time periods combined
----------------------------------------
Bayesian p-value: 0.5016
----------------------------------------
Individual time periods
----------------------------------------
Time Period 1 Bayesian p-value: 0.456
Time Period 2 Bayesian p-value: 0.456
Time Period 3 Bayesian p-value: 0.492
Time Period 4 Bayesian p-value: 0.404
Time Period 5 Bayesian p-value: 0.732
Time Period 6 Bayesian p-value: 0.64
Time Period 7 Bayesian p-value: 0.392
Time Period 8 Bayesian p-value: 0.46
Time Period 9 Bayesian p-value: 0.42
Time Period 10 Bayesian p-value: 0.564
Fit statistic: freeman-tukey As expected, we see all Bayesian p-values are quite close to the optimal 0.5.
Model Comparison using WAIC
Here we use WAIC to compare to a model that assumes the temporal trend is constant over space using the stPGOcc() function (or equivalently, we could use svcTPGOcc() with svc.cols = 1).
priors.stPGOcc <- list(beta.normal = priors$beta.normal,
alpha.normal = priors$alpha.normal,
sigma.sq.t.ig = priors$sigma.sq.t.ig,
rho.unif = priors$rho.unif,
phi.unif = c(3 / 1, 3 / 0.1),
sigma.sq.ig = c(2, 1))
out.stPGOcc <- stPGOcc(occ.formula = occ.formula,
det.formula = det.formula,
data = data,
inits = inits,
priors = priors.stPGOcc,
cov.model = cov.model,
n.neighbors = n.neighbors,
n.batch = n.batch,
batch.length = batch.length,
verbose = TRUE,
ar1 = ar1,
n.report = 200,
n.burn = n.burn,
n.thin = n.thin,
n.chains = 1)----------------------------------------
Preparing the data
----------------------------------------
----------------------------------------
Building the neighbor list
----------------------------------------
----------------------------------------
Building the neighbors of neighbors list
----------------------------------------
----------------------------------------
Model description
----------------------------------------
Spatial NNGP Multi-season Occupancy Model with Polya-Gamma latent
variable fit with 300 sites and 10 primary time periods.
Samples per chain: 15000 (600 batches of length 25)
Burn-in: 10000
Thinning Rate: 20
Number of Chains: 1
Total Posterior Samples: 250
Using the exponential spatial correlation model.
Using 5 nearest neighbors.
Using an AR(1) temporal autocorrelation matrix.
Source compiled with OpenMP support and model fit using 1 thread(s).
Adaptive Metropolis with target acceptance rate: 43.0
----------------------------------------
Chain 1
----------------------------------------
Sampling ...
Batch: 200 of 600, 33.33%
Parameter Acceptance Tuning
phi 44.0 0.66365
rho 52.0 1.09417
-------------------------------------------------
Batch: 400 of 600, 66.67%
Parameter Acceptance Tuning
phi 36.0 0.69073
rho 36.0 1.07251
-------------------------------------------------
Batch: 600 of 600, 100.00%
# SVC multi-season occupancy model
waicOcc(out.svc.trend) elpd pD WAIC
-1987.9099 118.2237 4212.2672
# Spatial multi-season occupancy model
waicOcc(out.stPGOcc) elpd pD WAIC
-2068.97911 67.75701 4273.47224 The lower WAIC value lends support for the model with the spatially-varying trend compared to a model where the trend is assumed constant over space.
Prediction
Finally, we can use the predict() function to predict relative occurrence and the effects of the covariates at new (and old) locations. The t.cols argument specifies the time periods over which you want to generate predictions. Here we will predict occupancy at the hold out locations for all ten time periods.
# Take a look at X.pred, the design matrix for the prediction locations
head(X.pred), , 1
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
[1,] 1 1 1 1 1 1 1 1 1 1
[2,] 1 1 1 1 1 1 1 1 1 1
[3,] 1 1 1 1 1 1 1 1 1 1
[4,] 1 1 1 1 1 1 1 1 1 1
[5,] 1 1 1 1 1 1 1 1 1 1
[6,] 1 1 1 1 1 1 1 1 1 1
, , 2
[,1] [,2] [,3] [,4] [,5] [,6] [,7]
[1,] -1.566503 -1.218391 -0.8702795 -0.5221677 -0.1740559 0.1740559 0.5221677
[2,] -1.566503 -1.218391 -0.8702795 -0.5221677 -0.1740559 0.1740559 0.5221677
[3,] -1.566503 -1.218391 -0.8702795 -0.5221677 -0.1740559 0.1740559 0.5221677
[4,] -1.566503 -1.218391 -0.8702795 -0.5221677 -0.1740559 0.1740559 0.5221677
[5,] -1.566503 -1.218391 -0.8702795 -0.5221677 -0.1740559 0.1740559 0.5221677
[6,] -1.566503 -1.218391 -0.8702795 -0.5221677 -0.1740559 0.1740559 0.5221677
[,8] [,9] [,10]
[1,] 0.8702795 1.218391 1.566503
[2,] 0.8702795 1.218391 1.566503
[3,] 0.8702795 1.218391 1.566503
[4,] 0.8702795 1.218391 1.566503
[5,] 0.8702795 1.218391 1.566503
[6,] 0.8702795 1.218391 1.566503
# Predict occupancy at the 400 new sites
out.pred <- predict(out.svc.trend, X.pred, coords.pred, t.cols = 1:n.time.max)----------------------------------------
Prediction description
----------------------------------------
Spatial NNGP Multi-season Occupancy model with Polya-Gamma latent
variable fit with 300 observations and 10 years.
Number of fixed covariates 2 (including intercept if specified).
Using the exponential spatial correlation model.
Using 5 nearest neighbors.
Number of MCMC samples 250.
Predicting at 100 non-sampled locations.
Source compiled with OpenMP support and model fit using 1 threads.
-------------------------------------------------
Predicting
-------------------------------------------------
Location: 100 of 100, 100.00%
Generating latent occupancy state
# Use the getSVCSamples() function to extract the SVC values
# at the prediction locations
svc.pred.samples <- getSVCSamples(out.svc.trend, pred.object = out.pred)
# True covariate effect values at new locations
trend.pred <- beta[2] + w.pred[, 2]
# Get meian and 95% CIs of the SVC for the trendariate effect
trend.pred.quants <- apply(svc.pred.samples[["trend"]], 2, quantile,
probs = c(0.025, 0.5, 0.975))
plot(trend.pred, trend.pred.quants[2, ], pch = 19,
ylim = c(min(trend.pred.quants[1, ]), max(trend.pred.quants[3, ])),
xlab = "True", ylab = "Fit", main = 'Spatially-varying trend')
abline(0, 1)
arrows(trend.pred, trend.pred.quants[2, ], trend.pred, col = 'gray',
trend.pred.quants[1, ], length = 0.02, angle = 90)
arrows(trend.pred, trend.pred.quants[1, ], trend.pred, col = 'gray',
trend.pred.quants[3, ], length = 0.02, angle = 90)
points(trend.pred, trend.pred.quants[2, ], pch = 19)
Spatially-varying coefficients multi-season generalized linear model
Basic model description
To model spatio-temporal patterns in occurrence when ignoring imperfect detection, we directly model the detection-nondetection data analogous to the SVC GLM, with all data and parameters now varying over time. The occurrence status of the species at site \(j\) during primary time period \(t\) (\(\psi_t(\boldsymbol{s}_j)\)) is analogous to the SVC multi-season occupancy model.
Simulating data with simTBinom()
The function simTBinom() simulates single-species multi-season detection-nondetection data assuming perfect detection. simTBinom() has the following arguments, which again closely resemble all other data simulation functions in spOccupancy.
simTBinom(J.x, J.y, n.time, weights, beta, sp.only = 0,
trend = TRUE, psi.RE = list(), sp = FALSE,
cov.model, sigma.sq, phi, nu, svc.cols = 1,
ar1 = FALSE, rho, sigma.sq.t, x.positive = FALSE, ...) All parameters are the same as those described previously for simTBinom() and simTOcc(). Below we simulate data from 400 sites over 10 years, with uneven sampling over the 10 years across the sites. As with the SVC multi-season occupancy model, our goal for this simulation is to assess a spatially-varying trend.
set.seed(500)
J.x <- 20
J.y <- 20
# Total number of sites
(J <- J.x * J.y)[1] 400
# Total number of time periods at each site
n.time <- sample(10, J, replace = TRUE)
n.time.max <- max(n.time)
# Number of Binomial replicates at each site/time period combination.
weights <- matrix(NA, J, max(n.time))
for (j in 1:J) {
weights[j, 1:n.time[j]] <- sample(4, n.time[j], replace = TRUE)
}
sp.only <- 0
# Simulate a trend parameter
trend <- TRUE
# Intercept and trend on relative occurrence
beta <- c(-0.5, -0.2)
# No unstructured random intercept on relative occurrence
psi.RE <- list()
# Spatial range for intercept and trend
phi <- c(3 / .8, 3 / .7)
# Spatial variance for intercept and trend
sigma.sq <- c(1, 0.5)
# Simulate data with an AR1 temporal process
ar1 <- TRUE
# Temporal variance of AR1 random effect
sigma.sq.t <- 0.5
# Temporal correlation of AR1 random effect
rho <- 0.25
# Spatially-varying coefficient columns
svc.cols <- c(1, 2)
# Simulate the data
dat <- simTBinom(J.x = J.x, J.y = J.y, n.time = n.time, weights = weights,
beta = beta, sp.only = sp.only, trend = trend, psi.RE = psi.RE,
sp = TRUE, sigma.sq = sigma.sq, phi = phi,
cov.model = 'exponential', svc.cols = svc.cols,
ar1 = ar1, sigma.sq.t = sigma.sq.t, rho = rho)
str(dat)List of 8
$ X : num [1:400, 1:10, 1:2] 1 1 1 1 1 1 1 1 1 1 ...
$ coords : num [1:400, 1:2] 0 0.0526 0.1053 0.1579 0.2105 ...
..- attr(*, "dimnames")=List of 2
.. ..$ : NULL
.. ..$ : chr [1:2] "Var1" "Var2"
$ psi : num [1:400, 1:10] 0.25286 0.06884 0.0829 0.03187 0.00752 ...
$ y : int [1:400, 1:10] 2 0 1 0 0 0 0 0 0 0 ...
$ w : num [1:400, 1:2] -1.246 -1.52 -1.033 -0.798 -1.627 ...
$ X.re : logi NA
$ beta.star: logi NA
$ eta : num [1:10, 1] -0.3429 0.0803 -0.5192 -0.1028 -0.0615 ...
# Plot average occupancy probability
plot(apply(dat$psi, 2, mean), pch = 19, xlab = 'Time Period',
ylab = 'Average Occupancy Probability',
cex = 1.5, frame = FALSE, ylim = c(0, 1))
We see a small decline in relative occurrence probability over the simulated ten years.
Below we extract the data we need to fit the model and create a plot of the spatially-varying trend across the study region
# Detection-nondetection data
y <- dat$y
# Design matrix for fixed effects
X <- dat$X
# Occurrence values
psi <- dat$psi
# Spatial coordinates
coords <- dat$coords
# Spatially-varying intercept and covariate effects
w <- dat$w
trend <- beta[2] + w[, 2]
plot.dat <- data.frame(x = coords[, 1],
y = coords[, 2],
trend = trend)
ggplot(plot.dat, aes(x = x, y = y, fill = trend)) +
geom_raster() +
scale_fill_gradient2(midpoint = 0, low = '#B2182B', mid = 'white', high = '#2166AC',
na.value = NA) +
theme_bw() +
labs(x = 'Easting', y = 'Northing', fill = 'Trend')
Lastly, we package up the data in a list required for fitting the multi-season spatially-varying coefficient GLM. We create a list with the following objects: y (the detection-nondetection data specified as a two-dimensional matrix), covs (the covariates specified as a list), weights (a matrix of binomial weights indicating the number of Bernoulli trials), and coords (a matrix of spatial coordinates). As before, we subset the data for prediction to only fit the model with 75% of the sites.
# Subset data for prediction.
# Split into fitting and prediction data set
pred.indx <- sample(1:J, round(J * .25), replace = FALSE)
y.fit <- y[-pred.indx, ]
y.pred <- y[pred.indx, ]
X.fit <- X[-pred.indx, , ]
X.pred <- X[pred.indx, , ]
coords.fit <- coords[-pred.indx, ]
coords.pred <- coords[pred.indx, ]
psi.fit <- psi[-pred.indx, ]
psi.pred <- psi[pred.indx, ]
w.fit <- w[-pred.indx, ]
w.pred <- w[pred.indx, ]
weights.fit <- weights[-pred.indx, ]
weights.pred <- weights[pred.indx, ]
# Package all data into a list
# Occurrence covariates
covs <- list(trend = X.fit[, , 2])
# Package into a list for spOccupancy
data.list <- list(y = y.fit,
covs = covs,
weights = weights.fit,
coords = coords.fit)
# Take a look at the data structure.
str(data.list)List of 4
$ y : int [1:300, 1:10] 2 0 1 0 0 0 0 0 0 0 ...
$ covs :List of 1
..$ trend: num [1:300, 1:10] -1.57 -1.57 -1.57 -1.57 -1.57 ...
$ weights: int [1:300, 1:10] 4 1 1 3 4 3 4 3 2 3 ...
$ coords : num [1:300, 1:2] 0 0.0526 0.1053 0.1579 0.2105 ...
..- attr(*, "dimnames")=List of 2
.. ..$ : NULL
.. ..$ : chr [1:2] "Var1" "Var2"Fitting spatially-varying coefficients generalized linear models with svcTPGBinom()
The function svcTPGBinom() fits spatially-varyig coefficients multi-season generalized linear models.
svcTPGBinom(formula, data, inits, priors,
tuning, svc.cols = 1, cov.model = 'exponential', NNGP = TRUE,
n.neighbors = 15, search.type = 'cb', n.batch,
batch.length, accept.rate = 0.43, n.omp.threads = 1,
verbose = TRUE, ar1 = FALSE, n.report = 100,
n.burn = round(.10 * n.batch * batch.length),
n.thin = 1, n.chains = 1, k.fold, k.fold.threads = 1,
k.fold.seed = 100, k.fold.only = FALSE, ...)There are no new arguments or details for fitting models with svcTPGBinom(), so we will go ahead and fit the model after specifying all the model arguments.
# Formatted data list
data <- data.list
str(data)List of 4
$ y : int [1:300, 1:10] 2 0 1 0 0 0 0 0 0 0 ...
$ covs :List of 1
..$ trend: num [1:300, 1:10] -1.57 -1.57 -1.57 -1.57 -1.57 ...
$ weights: int [1:300, 1:10] 4 1 1 3 4 3 4 3 2 3 ...
$ coords : num [1:300, 1:2] 0 0.0526 0.1053 0.1579 0.2105 ...
..- attr(*, "dimnames")=List of 2
.. ..$ : NULL
.. ..$ : chr [1:2] "Var1" "Var2"
# Model formula
formula <- ~ trend
# Initial values
inits <- list(beta = 0, sigma.sq = 1, phi = 3 / 0.5,
sigma.sq.t = 1.5, rho = 0.2)
# Priors
priors <- list(beta.normal = list(mean = 0, var = 2.72),
sigma.sq.t.ig = c(2, 0.5),
rho.unif = c(-1, 1),
sigma.sq.ig = list(a = 2, b = 1),
phi.unif = list(a = 3 / 1, b = 3 / .1))
# Exponential covariance function
cov.model <- 'exponential'
# Spatially-varying intercept and trend
svc.cols <- c(1, 2)
# Fit model with a NNGP with 5 neighbors
n.neighbors <- 5
# Include an AR(1) temporal covariance structure.
ar1 <- TRUE
# MCMC parameters
n.batch <- 600
batch.length <- 25
# Total number of samples
n.batch * batch.length[1] 15000
n.burn <- 10000
n.thin <- 20
# Approx. run time: ~ 1.5 min
out.svc.trend <- svcTPGBinom(formula = formula,
data = data,
inits = inits,
priors = priors,
cov.model = cov.model,
svc.cols = svc.cols,
n.neighbors = n.neighbors,
n.batch = n.batch,
batch.length = batch.length,
verbose = TRUE,
ar1 = ar1,
n.report = 200,
n.burn = n.burn,
n.thin = n.thin,
n.chains = 1)----------------------------------------
Preparing the data
----------------------------------------
----------------------------------------
Building the neighbor list
----------------------------------------
----------------------------------------
Building the neighbors of neighbors list
----------------------------------------
----------------------------------------
Model description
----------------------------------------
Spatial NNGP Multi-season Binomial Model with Polya-Gamma latent
variable fit with 300 sites and 10 years.
Samples per chain: 15000 (600 batches of length 25)
Burn-in: 10000
Thinning Rate: 20
Number of Chains: 1
Total Posterior Samples: 250
Number of spatially-varying coefficients: 2
Using the exponential spatial correlation model.
Using 5 nearest neighbors.
Using an AR(1) temporal autocorrelation matrix.
Source compiled with OpenMP support and model fit using 1 thread(s).
Adaptive Metropolis with target acceptance rate: 43.0
----------------------------------------
Chain 1
----------------------------------------
Sampling ...
Batch: 200 of 600, 33.33%
Coefficient Parameter Acceptance Tuning
1 phi 28.0 0.58860
2 phi 60.0 0.70469
1 rho 32.0 1.28403
-------------------------------------------------
Batch: 400 of 600, 66.67%
Coefficient Parameter Acceptance Tuning
1 phi 40.0 0.53259
2 phi 40.0 0.81058
1 rho 28.0 1.30996
-------------------------------------------------
Batch: 600 of 600, 100.00%
summary(out.svc.trend)
Call:
svcTPGBinom(formula = formula, data = data, inits = inits, priors = priors,
svc.cols = svc.cols, cov.model = cov.model, n.neighbors = n.neighbors,
n.batch = n.batch, batch.length = batch.length, verbose = TRUE,
ar1 = ar1, n.report = 200, n.burn = n.burn, n.thin = n.thin,
n.chains = 1)
Samples per Chain: 15000
Burn-in: 10000
Thinning Rate: 20
Number of Chains: 1
Total Posterior Samples: 250
Run Time (min): 1.8373
Occurrence (logit scale):
Mean SD 2.5% 50% 97.5% Rhat ESS
(Intercept) -0.1490 0.2990 -0.7855 -0.131 0.4739 NA 32
trend -0.2644 0.3101 -0.7975 -0.279 0.3572 NA 20
Spatio-temporal Covariance:
Mean SD 2.5% 50% 97.5% Rhat ESS
sigma.sq-(Intercept) 0.7342 0.1810 0.4590 0.7109 1.1274 NA 141
sigma.sq-trend 0.7398 0.2080 0.4181 0.7101 1.1573 NA 120
phi-(Intercept) 7.6085 2.4949 3.7415 7.2722 13.0823 NA 160
phi-trend 5.0191 1.4382 3.0930 4.7755 8.6837 NA 81
sigma.sq.t 0.2524 0.1365 0.0930 0.2154 0.6464 NA 250
rho -0.0135 0.2508 -0.5345 -0.0059 0.4355 NA 250
# Compare estimates to true values
beta[1] -0.5 -0.2
alpha[1] 0.9 -0.3
sigma.sq[1] 1.0 0.5
phi[1] 3.750000 4.285714
rho[1] 0.25
sigma.sq.t[1] 0.5Let’s extract the spatially-varying intercept and spatially-varying trend estimates from our model using getSVCSamples() and compare them to the true values.
# Intercept ---------------------------------------------------------------
svc.samples <- getSVCSamples(out.svc.trend)
int.quants <- apply(svc.samples[["(Intercept)"]], 2, quantile,
probs = c(0.025, 0.5, 0.975))
svc.true.fit <- beta + w.fit
plot(svc.true.fit[, 1], int.quants[2, ], pch = 19,
ylim = c(min(int.quants[1, ]), max(int.quants[3, ])),
xlab = "True", ylab = "Fit", main = 'Intercept')
abline(0, 1)
arrows(svc.true.fit[, 1], int.quants[2, ], svc.true.fit[, 1], col = 'gray',
int.quants[1, ], length = 0.02, angle = 90)
arrows(svc.true.fit[, 1], int.quants[1, ], svc.true.fit[, 1], col = 'gray',
int.quants[3, ], length = 0.02, angle = 90)
points(svc.true.fit[, 1], int.quants[2, ], pch = 19)
# Trend -------------------------------------------------------------------
trend.quants <- apply(svc.samples[["trend"]], 2, quantile,
probs = c(0.025, 0.5, 0.975))
plot(svc.true.fit[, 2], trend.quants[2, ], pch = 19,
ylim = c(min(trend.quants[1, ]), max(trend.quants[3, ])),
xlab = "True", ylab = "Fit", main = 'Spatially-Varying Trends')
abline(0, 1)
arrows(svc.true.fit[, 2], trend.quants[2, ], svc.true.fit[, 2], col = 'gray',
trend.quants[1, ], length = 0.02, angle = 90)
arrows(svc.true.fit[, 2], trend.quants[1, ], svc.true.fit[, 2], col = 'gray',
trend.quants[3, ], length = 0.02, angle = 90)
points(svc.true.fit[, 2], trend.quants[2, ], pch = 19)
We see our model does a good job of recovering the true parameter values for both the spatially-varying intercept and spatially-varying trend.
Model Comparison using WAIC
Here we use WAIC to compare our full SVC model to a model that assumes the temporal trend is constant over space by setting svc.cols = 1.
out.trend.const <- svcTPGBinom(formula = formula,
data = data,
inits = inits,
priors = priors,
svc.cols = 1,
cov.model = cov.model,
n.neighbors = n.neighbors,
n.batch = n.batch,
batch.length = batch.length,
verbose = TRUE,
ar1 = ar1,
n.report = 200,
n.burn = n.burn,
n.thin = n.thin,
n.chains = 1)----------------------------------------
Preparing the data
----------------------------------------
----------------------------------------
Building the neighbor list
----------------------------------------
----------------------------------------
Building the neighbors of neighbors list
----------------------------------------
----------------------------------------
Model description
----------------------------------------
Spatial NNGP Multi-season Binomial Model with Polya-Gamma latent
variable fit with 300 sites and 10 years.
Samples per chain: 15000 (600 batches of length 25)
Burn-in: 10000
Thinning Rate: 20
Number of Chains: 1
Total Posterior Samples: 250
Number of spatially-varying coefficients: 1
Using the exponential spatial correlation model.
Using 5 nearest neighbors.
Using an AR(1) temporal autocorrelation matrix.
Source compiled with OpenMP support and model fit using 1 thread(s).
Adaptive Metropolis with target acceptance rate: 43.0
----------------------------------------
Chain 1
----------------------------------------
Sampling ...
Batch: 200 of 600, 33.33%
Coefficient Parameter Acceptance Tuning
1 phi 48.0 0.58860
1 rho 36.0 1.18530
-------------------------------------------------
Batch: 400 of 600, 66.67%
Coefficient Parameter Acceptance Tuning
1 phi 44.0 0.58860
1 rho 36.0 1.30996
-------------------------------------------------
Batch: 600 of 600, 100.00%
# SVC multi-season GLM
waicOcc(out.svc.trend) elpd pD WAIC
-1517.713 198.032 3431.490
# Spatial multi-season GLM
waicOcc(out.trend.const) elpd pD WAIC
-1674.1340 162.3747 3673.0174 Prediction
Finally, we use the predict() function to predict relative occurrence and the effects of the covariates at new (and old) locations over the simulated ten tme periods.
# Take a look at X.pred, the design matrix for the prediction locations
head(X.pred), , 1
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
[1,] 1 1 1 1 1 1 1 1 1 1
[2,] 1 1 1 1 1 1 1 1 1 1
[3,] 1 1 1 1 1 1 1 1 1 1
[4,] 1 1 1 1 1 1 1 1 1 1
[5,] 1 1 1 1 1 1 1 1 1 1
[6,] 1 1 1 1 1 1 1 1 1 1
, , 2
[,1] [,2] [,3] [,4] [,5] [,6] [,7]
[1,] -1.566503 -1.218391 -0.8702795 -0.5221677 -0.1740559 0.1740559 0.5221677
[2,] -1.566503 -1.218391 -0.8702795 -0.5221677 -0.1740559 0.1740559 0.5221677
[3,] -1.566503 -1.218391 -0.8702795 -0.5221677 -0.1740559 0.1740559 0.5221677
[4,] -1.566503 -1.218391 -0.8702795 -0.5221677 -0.1740559 0.1740559 0.5221677
[5,] -1.566503 -1.218391 -0.8702795 -0.5221677 -0.1740559 0.1740559 0.5221677
[6,] -1.566503 -1.218391 -0.8702795 -0.5221677 -0.1740559 0.1740559 0.5221677
[,8] [,9] [,10]
[1,] 0.8702795 1.218391 1.566503
[2,] 0.8702795 1.218391 1.566503
[3,] 0.8702795 1.218391 1.566503
[4,] 0.8702795 1.218391 1.566503
[5,] 0.8702795 1.218391 1.566503
[6,] 0.8702795 1.218391 1.566503
# Predict occupancy at the 400 new sites
out.pred <- predict(out.svc.trend, X.pred, coords.pred, t.cols = 1:n.time.max,
weights = weights.pred)----------------------------------------
Prediction description
----------------------------------------
Spatial NNGP Multi-season Occupancy model with Polya-Gamma latent
variable fit with 300 observations and 10 years.
Number of fixed covariates 2 (including intercept if specified).
Using the exponential spatial correlation model.
Using 5 nearest neighbors.
Number of MCMC samples 250.
Predicting at 100 non-sampled locations.
Source compiled with OpenMP support and model fit using 1 threads.
-------------------------------------------------
Predicting
-------------------------------------------------
Location: 100 of 100, 100.00%
Generating latent occupancy state
# Use the getSVCSamples() function to extract the SVC values
# at the prediction locations
svc.pred.samples <- getSVCSamples(out.svc.trend, pred.object = out.pred)
# True covariate effect values at new locations
trend.pred <- beta[2] + w.pred[, 2]
# Get meian and 95% CIs of the SVC for the trendariate effect
trend.pred.quants <- apply(svc.pred.samples[["trend"]], 2, quantile,
probs = c(0.025, 0.5, 0.975))
plot(trend.pred, trend.pred.quants[2, ], pch = 19,
ylim = c(min(trend.pred.quants[1, ]), max(trend.pred.quants[3, ])),
xlab = "True", ylab = "Fit", main = 'Spatially-varying trend')
abline(0, 1)
arrows(trend.pred, trend.pred.quants[2, ], trend.pred, col = 'gray',
trend.pred.quants[1, ], length = 0.02, angle = 90)
arrows(trend.pred, trend.pred.quants[1, ], trend.pred, col = 'gray',
trend.pred.quants[3, ], length = 0.02, angle = 90)
points(trend.pred, trend.pred.quants[2, ], pch = 19)
References
Banerjee, Sudipto, Bradley P Carlin, and Alan E Gelfand. 2003. Hierarchical Modeling and Analysis for Spatial Data. Chapman; Hall/CRC.
Chen, Guoke, Marc Kéry, Matthias Plattner, Keping Ma, and Beth Gardner. 2013. “Imperfect Detection Is the Rule Rather Than the Exception in Plant Distribution Studies.” Journal of Ecology 101 (1): 183–91.
Datta, Abhirup, Sudipto Banerjee, Andrew O Finley, and Alan E Gelfand. 2016. “Hierarchical Nearest-Neighbor Gaussian Process Models for Large Geostatistical Datasets.” Journal of the American Statistical Association 111 (514): 800–812.
Doser, Jeffrey W, Andrew O Finley, Marc Kéry, and Elise F Zipkin. 2022. “spOccupancy: An R package for single-species, multi-species, and integrated spatial occupancy models.” Methods in Ecology and Evolution 13 (8): 1670–8.
Finley, Andrew O. 2011. “Comparing Spatially-Varying Coefficients Models for Analysis of Ecological Data with Non-Stationary and Anisotropic Residual Dependence.” Methods in Ecology and Evolution 2 (2): 143–54.
Gelfand, Alan E, Hyon-Jung Kim, CF Sirmans, and Sudipto Banerjee. 2003. “Spatial Modeling with Spatially Varying Coefficient Processes.” Journal of the American Statistical Association 98 (462): 387–96.
Hooten, Mevin B, and Trevor J Hefley. 2019. Bringing Bayesian Models to Life. CRC Press.
Hooten, Mevin B, and N Thompson Hobbs. 2015. “A Guide to Bayesian Model Selection for Ecologists.” Ecological Monographs 85 (1): 3–28.
Isaac, Nick JB, Arco J van Strien, Tom A August, Marnix P de Zeeuw, and David B Roy. 2014. “Statistics for Citizen Science: Extracting Signals of Change from Noisy Ecological Data.” Methods in Ecology and Evolution 5 (10): 1052–60.
Johnson, Devin S, Paul B Conn, Mevin B Hooten, Justina C Ray, and Bruce A Pond. 2013. “Spatial Occupancy Models for Large Data Sets.” Ecology 94 (4): 801–8.
Lele, Subhash R, Monica Moreno, and Erin Bayne. 2012. “Dealing with Detection Error in Site Occupancy Surveys: What Can We Do with a Single Survey?” Journal of Plant Ecology 5 (1): 22–31.
MacKenzie, Darryl I, James D Nichols, Gideon B Lachman, Sam Droege, J Andrew Royle, and Catherine A Langtimm. 2002. “Estimating Site Occupancy Rates When Detection Probabilities Are Less Than One.” Ecology 83 (8): 2248–55.
Pease, Brent S, Krishna Pacifici, and Roland Kays. 2022. “Exploring Spatial Nonstationarity for Four Mammal Species Reveals Regional Variation in Environmental Relationships.” Ecosphere 13 (8): e4166.
Polson, Nicholas G, James G Scott, and Jesse Windle. 2013. “Bayesian Inference for Logistic Models Using Pólya–Gamma Latent Variables.” Journal of the American Statistical Association 108 (504): 1339–49.
Roberts, Gareth O, and Jeffrey S Rosenthal. 2009. “Examples of Adaptive Mcmc.” Journal of Computational and Graphical Statistics 18 (2): 349–67.
Rollinson, Christine R, Andrew O Finley, M Ross Alexander, Sudipto Banerjee, Kelly-Ann Dixon Hamil, Lauren E Koenig, Dexter Henry Locke, et al. 2021. “Working Across Space and Time: Nonstationarity in Ecological Research and Application.” Frontiers in Ecology and the Environment 19 (1): 66–72.
Rushing, Clark S, J Andrew Royle, David J Ziolkowski Jr, and Keith L Pardieck. 2020. “Migratory Behavior and Winter Geography Drive Differential Range Shifts of Eastern Birds in Response to Recent Climate Change.” Proceedings of the National Academy of Sciences 117 (23): 12897–12903.
Tyre, Andrew J, Brigitte Tenhumberg, Scott A Field, Darren Niejalke, Kirsten Parris, and Hugh P Possingham. 2003. “Improving Precision and Reducing Bias in Biological Surveys: Estimating False-Negative Error Rates.” Ecological Applications 13 (6): 1790–1801.
Wright, Wilson J, Kathryn M Irvine, Thomas J Rodhouse, and Andrea R Litt. 2021. “Spatial Gaussian Processes Improve Multi-Species Occupancy Models When Range Boundaries Are Uncertain and Nonoverlapping.” Ecology and Evolution.