In doserjef/spOccupancy: Single-Species, Multi-Species, and Integrated Spatial Occupancy Models

options(rmarkdown.html_vignette.check_title = FALSE)
knitr::opts_chunk$set(
  comment = "", eval = FALSE
)

\newcommand{\bm}{\boldsymbol} \newcommand{\pg}{\text{P{\'o}lya-Gamma}}

Introduction

This document provides statistical details on the MCMC algorithms used to fit spatially-varying coefficient (SVC) models in spOccupancy. In particular, we discuss the Gibbs samplers for the following two models presented in [@doser2024SVC]:

1. A single-species SVC occupancy model using svcPGOcc()
2. A multi-species SVC occupancy model using svcMsPGOcc()

Single-species spatially-varying coefficient occupancy model

Model description

Let $\bm{s}_j$ denote the spatial coordinates of site $j$, where $j = 1, \dots, J$. We define $z(\bm{s}_j)$ as the true presence (1) or absence (0) of the target species at site $j$ with spatial coordinates $\bm{s}_j$. We model $z(\bm{s}_j)$ as

\begin{equation}\label{z} z(\bm{s}_j) \sim \text{Bernoulli}(\psi(\bm{s}_j)), \end{equation}

where $\psi(\bm{s}_j)$ is the occupancy probability of the species at site $j$. We model $\psi(\bm{s}_j)$ according to

\begin{equation}\label{psi} \text{logit}(\psi(\bm{s}j)) = (\beta_1 + \delta_1\text{w}_1(\bm{s}_j)) + \sum{h = 2}^H\text{x}_h(\bm{s}_j){\beta_h + \delta_h\text{w}_h(\bm{s}_j)}, \end{equation}

where $\beta_1$ is an intercept, $\text{x}h(\bm{s}_j)$ is the $h$th covariate with $h = 2, \dots, H$, $\beta_h$ is the non-spatial effect of covariate $\text{x}_h(\bm{s}_j)$, and $\text{w}_1(\bm{s}_j)$ and $\text{w}_h(\bm{s}_j)$ are spatially-varying effects for the intercept and covariates, respectively. We use indicator variables $\delta_h$ for $h = 1, \dots, H$ to indicate those covariates whose effects vary spatially ($\delta_h = 1$) and those whose effects are assumed constant ($\delta_h = 0$). Note that the model reduces to a traditional single-species occupancy model when $\delta_h = 0$ for all $h$ and a spatial occupancy model [@johnson2013spatial; @doser2022spoccupancy] when $\delta_1 = 1$ and $\delta_h = 0$ for all $h > 1$. For later use, define $\tilde{H}$ as the total number of spatially-varying effects estimated in the model (i.e., $\tilde{H} = \sum{h = 1}^H\delta_h$), define $\tilde{\textbf{x}}(\bm{s}_j)$ as the $\tilde{H} \times 1$ vector of covariates at location $j$ (including an intercept if applicable), and define $\tilde{\beta}_h(\bm{s}_j) = \beta_h + \text{w}_h(\bm{s}_j)$ as the spatially-varying coefficients for those effects with $\delta_h = 1$.

Let $\mathcal{L} = {\bm{s}_1, \bm{s}_2, \dots, \bm{s}_J}$ be the set of sampled spatial locations, and define $\textbf{w}_h$ as a $J \times 1$ vector of the spatial random effects for covariate $h$ for each of the $\tilde{H}$ effects with corresponding $\delta_h = 1$. The spatially-varying effects $\textbf{w}_h$ serve as local adjustments of the covariate effects (or intercept) at each site $j$ from the overall non-spatial effect $\beta_h$, resulting in the covariate having a unique effect (i.e., $\tilde{\beta}_h(\bm{s}_j))$ on species occupancy probability at each site $j$. Following @gelfand2003spatial, we envision each $h = 1, \dots, \tilde{H}$ spatially-varying effect $\text{w}_h(\bm{s}_j)$ as a realization of a smooth latent surface ${\text{w}_h(\bm{s}) \mid \bm{s} \in \mathcal{D}}$, where $\mathcal{D}$ is the geographical domain of interest. Following standard approaches for modeling species distributions (e.g., @latimer2006building), we use Gaussian Processes (GPs) to model each of the $\tilde{H}$ smooth functions across the spatial domain. By definition, a GP model for the $h$th spatially-varying surface ${\text{w}_h(\bm{s})}$ implies that for any finite set of locations $\mathcal{L}$, the vector of random effects $\textbf{w}_h$ follows a zero-mean multivariate Gaussian distribution with a $J \times J$ covariance matrix $\bm{C}_h(\bm{s}, \bm{s}', \bm{\theta}_h)$ that is a function of the distances between any pair of site coordinates $\bm{s}$ and $\bm{s}'$ and a set of parameters ($\bm{\theta}_h$) that govern the spatial process according to a parametric covariance function. In our subsequent simulations and case study, we use an exponential covariance function such that $\bm{\theta}_h = {\sigma^2_h, \phi_h}$, where $\sigma^2_h$ is a spatial variance parameter and $\phi_h$ is a spatial decay parameter. Large values of $\sigma^2_h$ indicate large variation in the magnitude of a covariate effect across space, while values of $\sigma^2_h$ close to 0 suggest little spatial variability in the magnitude of the effect. $\phi_h$ controls the distance-dependent decay of the spatial dependence in the covariate effect and is inversely related to the spatial range, such that when $\phi_h$ is small, the covariate effect has a larger range of spatial dependence and varies more smoothly across space compared to larger values of $\phi_h$. When using an exponential correlation function, the effective spatial range, or the distance at which the spatial correlation between points drops to 0.05 [@banerjee2014hierarchical], corresponds to approximately 3 / $\phi$ (i.e., since $3 \approx -\text{log}(0.05))$.

Both frequentist and Bayesian estimation of the model defined by (\ref{z}) and (\ref{psi}) requires taking the inverse and determinant of $\tilde{H}$ dense $J \times J$ covariance matrices (i.e., $\bm{C}h(\bm{s}, \bm{s}', \bm{\theta}_h)$), which involves $O(J^3)$ computations for each of the $\tilde{H}$ spatially-varying coefficients (floating point operations or FLOPs), which quickly renders such an approach impractical for even moderately sized data sets (i.e., hundreds of spatial locations). Here we replace the GP prior for the spatially-varying coefficients with a Nearest Neighbor Gaussian Process (NNGP) prior [@datta2016hierarchical]. The NNGP is a valid GP that is based on writing the full multivariate Gaussian distribution for $\textbf{w}_h$ as a product of conditional densities, such that \begin{equation}\label{wConditional} p(\textbf{w}_h) = p(\text{w}_h(\bm{s}_1)) \cdot p(\text{w}_h(\bm{s}_2) \mid \text{w}_h(\bm{s}_1)) \cdots p(\text{w}_h(\bm{s}_J) \mid \text{w}_h(\bm{s}{J - 1}), \dots, \text{w}_h(\bm{s}_1)), \end{equation}

where $p(\cdot)$ denotes a probability density function. The NNGP prior achieves computational efficiency by replacing the conditioning sets on the right-hand side of (\ref{wConditional}) with a set of new conditioning sets, whose maximum size is determined by a pre-specified number of neighbors, $m$, where $m << J$. @datta2016hierarchical showed that $m = 15$ provides nearly identical inference to the full GP under a variety of scenarios. Let $N(\bm{s}j)$ denote the set of at most $m$ neighbors for location $\bm{s}_j$. Following @vecchia1988estimation, we set $N(\bm{s}_j)$ to be the set of at most $m$ nearest neighbors of $\bm{s}_j$ from ${\bm{s}_1, \bm{s}_2, \dots, \bm{s}{j - 1}}$ with respect to Euclidean distance. Note that this requires the set of $\mathcal{L}$ locations to have some prespecified ordering, where here we order the coordinates along the horizontal axis. Through careful construction of the neighbor sets and set of spatial locations as a directed acyclic graph, Gaussian distribution theory reveals the NNGP prior yields a new joint density for $\textbf{w}h$, denoted $\tilde{p}(\textbf{w}_h)$. Let $\textbf{w}_h(N(\bm{s}_j))$ denote the at most $m$ realizations of the $h^{\text{th}}$ NNGP at the locations in the neighbor set $N(\bm{s}_j)$. Let $C(\cdot, \bm{\theta}_h)$ denote the covariance function of the original GP from which the $h^{\text{th}}$ NNGP is derived. For any two sets $A_1$ and $A_2$, define $\text{C}{A_1, A_2}(\bm{\theta}h)$ as the covariance matrix between the observations in $A_1$ and $A_2$ for the $h^{\text{th}}$ GP. For all $h = 1, \dots, \tilde{H}$, our NNGP prior for $\textbf{w}_h$ thus takes the form \begin{equation}\label{NNGP} \tilde{p}(\textbf{w}_h) = \prod{j = 1}^J \text{Normal}(\text{w}h(\bm{s}_j) \mid \textbf{b}_h(\bm{s}_j)\textbf{w}_h(N(\bm{s}_j)), \text{f}_h(\bm{s}_j)), \end{equation} where $\textbf{b}_h(\bm{s}_j)$ is defined as \begin{equation}\label{bNNGP} \textbf{b}_h(\bm{s}_j) = \textbf{C}{\bm{s}j, N(\bm{s}_j)}(\bm{\theta}_h)\textbf{C}^{-1}{N(\bm{s}j), N(\bm{s}_j)}(\bm{\theta}_h), \end{equation} and $\text{f}_h(\bm{s}_j)$ is defined as \begin{equation}\label{fNNGP} \text{f}_h(\bm{s}_j) = \textbf{C}{\bm{s}j, \bm{s}_j}(\bm{\theta}_h) - \textbf{C}{\bm{s}j, N(\bm{s}_j)}(\bm{\theta}_h)\textbf{C}^{-1}{N(\bm{s}j), N(\bm{s}_j)}(\bm{\theta}_h)\textbf{C}{N(\bm{s}_j), \bm{s}_j}(\bm{\theta}_h). \end{equation}

To account for imperfect detection in an occupancy modeling framework, $k = 1, \dots, K(\bm{s}_j)$ sampling replicates are obtained at each site $j$ to estimate whether a nondetection of the target species is truly an absence [@mackenzie2002; @tyre2003improving]. We model the observed detection (1) or nondetection (0) of a study species at site $j$, denoted $y_k(\bm{s}_j)$, conditional on the true latent occupancy process $z(\bm{s}_j)$, following

\begin{equation}\label{yDet} y_{k}(\bm{s}j) \mid z(\bm{s}_j) \sim \text{Bernoulli}(p{k}(\bm{s}_j)z(\bm{s}_j)), \end{equation}

where $p_k(\bm{s}_j)$ 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}(\bm{s}j)) = \textbf{v}{k}(\bm{s}_j)^\top\bm{\alpha}, \end{equation}

where $\bm{\alpha}$ is a vector of regression coefficients (including an intercept) that describe the effect of site and/or observation covariates $\textbf{v}_{k}(\bm{s}_j)$ on detection. Note that following the standard occupancy model, we assume independence between the replicate surveys, conditional on the true occupancy status $z(\bm{s}_j)$ and covariates \textbf{v}$_k(\bm{s}_j)$, the true occupancy status does not change over the $K(\bm{s}_j)$ replicate surveys, and there are no false positives (i.e., if $z(\bm{s}_j) = 0$, then $y_k(\bm{s}_j) = 0$ for all $k$). We assume effects of covariates on detection probability are constant over space, but in principle spatially-varying covariate effects could be added to the detection model using the same process described above.

We assign independent Gaussian priors to all non-spatial regression coefficients ($\bm{\beta}$ and $\bm{\alpha}$), independent inverse-Gamma priors to the spatial variance parameters ($\sigma^2_h$) for each spatially-varying effect, and independent uniform priors for the spatial decay parameters ($\phi_h$).

Gibbs sampler

Update occurrence auxiliary variables ($\omega_{\beta}(\bm{s}_j)$)

We first sample the occurrence auxiliary variable ($\omega_{\beta}(\bm{s}_j)$ individually for each site $j$. Our full conditional is

\begin{equation} \omega_{\beta}(\bm{s}_j) \mid \cdot \sim \text{PG}(1, \textbf{x}(\bm{s}_j)^{\top}\bm{\beta} + \tilde{\textbf{x}}(\bm{s}_j)^\top\textbf{w}(\bm{s}_j)). \end{equation}

Update detection auxiliary variables ($\omega_{k, \alpha}(\bm{s}_j)$)

We next update the latent \pg auxiliary variable for the detection process, $\omega_{k, \alpha}(\bm{s}j)$, for each replicate $k$ at each site $j$. Note that we only need to sample $\omega{k, \alpha}(\bm{s}j)$ when $\sum{k = 1}^{K_j}y_{k}(\bm{s}_j) = 0$. Following @polson2013, we have

\begin{equation} \omega_{k, \alpha}(\bm{s}_j) \mid \cdot \sim \text{PG}(1, \textbf{v}(\bm{s}_j)^\top\bm{\alpha}). \end{equation}

Update non-spatial occurrence regression coefficients ($\bm{\beta}$)

Define $\tilde{\text{w}}(\bm{s}j) = \tilde{\text{x}}(\bm{s}_j)^\top \textbf{w}(\bm{s}_j)$, and similarly let $\tilde{\textbf{w}}$ be the $J \times 1$ vector of $\tilde{\text{w}}(\bm{s}_j)$ for all $j = 1, \dots, J$. Let $\bm{\beta} \sim N(\bm{\mu}{0, \beta}, \bm{\Sigma}{\beta})$ denote the prior distribution for $\bm{\beta}$, with $\Sigma{\beta}$ being a diagonal matrix. The \pg scheme induces a Gibbs update for the occurrence regression coefficients, which are updated at each iteration according to

\begin{equation} \label{betaSampler} \bm{\beta} \mid \cdot \sim \text{Normal}\Big([\bm{\Sigma}{\beta}^{-1} + \textbf{X}^{\top}\bm{S}{\beta}\textbf{X}]^{-1}[\textbf{X}^{\top}(\bm{z} - 0.5 \bm{1}J - \bm{S}{\beta}\tilde{\textbf{w}}) + \bm{\Sigma}{\beta}^{-1}\bm{\mu}{0, \beta}], [\bm{\Sigma}{\beta}^{-1} + \textbf{X}^{\top}\bm{S}{\beta}\textbf{X}]^{-1}\Big), \end{equation}

where $\bm{S}{\beta}$ is a diagonal $J \times J$ matrix with diagonal entries equal to the latent PG variable values ($\omega{1, \beta}, \dots, \omega_{J, \beta})$ and $\bm{z}$ is the $J \times 1$ vector of the latent occurrence values.

Update detection regression coefficients ($\bm{\alpha}$)

Similarly, let $\bm{\alpha} \sim N(\bm{\mu}{0, \alpha}, \bm{\Sigma}{\alpha})$ denote the prior distribution for $\bm{\alpha}$, with $\bm{\Sigma}_{\alpha}$ being a diagonal matrix. We sample the detection regression coefficients $\bm{\alpha}$ from

\begin{equation} \label{alphaSampler} \bm{\alpha} \mid \cdot \sim \text{Normal}\Big([\bm{\Sigma}{\alpha}^{-1} + \bm{\tilde{V}}^{\top}\bm{S}{\alpha}\bm{\tilde{V}}]^{-1}[\bm{\tilde{V}}^{\top}(\bm{\tilde{y}} - 0.5 \bm{1}{J^*}) + \bm{\Sigma}{\alpha}^{-1}\bm{\mu}{0, \alpha}], [\bm{\Sigma}{\alpha}^{-1} + \bm{\tilde{V}}^{\top}\bm{S}_{\alpha}\bm{\tilde{V}}]^{-1}\Big). \end{equation}

The detection regression coefficients $\bm{\alpha}$ are only informed by the locations where $z(\bm{s}j) = 1$, since we assume no false positive detections in the standard occupancy model. We define $J^$ as the total number of sites at the current iteration of the MCMC with $z(\bm{s}j) = 1$. $\bm{S}{\alpha}$ is a diagonal matrix with diagonal entries equal to the latent \pg variable values ($\omega_{1, 1, \alpha}, \dots, \omega_{J^, K{J^*}, \alpha}$). The matrix $\bm{\tilde{V}}$ is the matrix of detection covariates associated with the sites where $z(\bm{s}_j) = 1$. Similarly, $\bm{\tilde{y}}$ is a vector of stacked detection-nondetection data values at the entries associated with $z(\bm{s}_j) = 1$.

Update spatially-varying coefficients ($\textbf{w}(\bm{s}_j)$)

Next we update the spatially-varying coefficients $\textbf{w}(\bm{s}j)$ sequentially for each of the $j = 1, \dots, J$ spatial locations. Let $N(\bm{s}_j)$ denote the set of $m$ nearest neighbors of $\bm{s}_j$ among $\bm{s}_1, \bm{s}_2, \dots, \bm{s}{J - 1}$. For all model implementations in \texttt{spOccupancy}, we order the spatial locations by the horizontal axis. Let $\textbf{w}h(N(\bm{s}_j))$ denote the $m$ realizations of the $h^{\text{th}}$ NNGP at the locations in $N(\bm{s}_j)$, where $h$ iterates across the total $\tilde{H}$ spatially-varying coefficients. Let $C(\cdot, \bm{\theta}_h)$ denote the covariance function of the original Gaussian Process (GP) from which the $h^{\text{th}}$ NNGP is derived. For any two sets $A_1$ and $A_2$, define $\text{C}{A_1, A_2}(\bm{\theta}_h)$ as the covariance matrix between the observations in $A_1$ and $A_2$ for the $h^{\text{th}}$ GP. For $j \geq 1$, we have

\begin{equation}\label{Eq: bNNGP} \textbf{b}h(\bm{s}_j) = \textbf{C}{\bm{s}j, N(\bm{s}_j)}(\bm{\theta}_h)\textbf{C}^{-1}{N(\bm{s}_j), N(\bm{s}_j)}(\bm{\theta}_h), \end{equation}

where $\textbf{b}_h(\bm{s}_1) = \bm{0}$ for all $h = 1, \dots, \tilde{H}$. Further, we have

\begin{equation}\label{Eq: fNNGP} f_h(\bm{s}j) = \textbf{C}{\bm{s}j, \bm{s}_j}(\bm{\theta}_h) - \textbf{C}{\bm{s}j, N(\bm{s}_j)}(\bm{\theta}_h)\textbf{C}^{-1}{N(\bm{s}j), N(\bm{s}_j)}(\bm{\theta}_h)\textbf{C}{N(\bm{s}_j), \bm{s}_j}(\bm{\theta}_h). \end{equation}

For any two locations $\bm{s}1$ and $\bm{s}_2$, if $\bm{s}_1 \in N(\bm{s}_2)$ and is the $l^{\text{th}}$ member of $N(\bm{s}_2)$, then define $b_h(\bm{s}_2, \bm{s}_1)$ as the $l^{\text{th}}$ entry of $\textbf{b}_h(\bm{s}_2)$. Let $U(\bm{s}_1) = {\bm{s}_2 \in S \mid \bm{s}_1 \in N(\bm{s}_2)}$ be the collection of locations $\bm{s}_2$ for which $\bm{s}_1$ is a neighbor, where $S$ is the set of all $J$ spatial locations. For every $\bm{s}_2 \in U(\bm{s}_1)$, define $a_h(\bm{s}_2, \bm{s}_1) = \text{w}_h(\bm{s}_2) - \sum{\bm{s} \in N(\bm{s}_2), \bm{s} \neq \bm{s}_2} \text{w}_h(\bm{s})b_h(\bm{s}_2, \bm{s})$. Extending this to matrix notation, let $\bm{B}(\bm{s}_j)$ be a $\tilde{H} \times m\tilde{H}$ block matrix, with $\tilde{H} \times \tilde{H}$ diagonal blocks containing the elements of $\bm{b}_h(\bm{s}_j)$ for each of the $h = 1, \dots \tilde{H}$ spatial NNGPs for each of the specific $m$ neighbors. Let $\bm{F}(\bm{s}_j)$ be a $\tilde{H} \times \tilde{H}$ diagonal matrix with diagonal elements of $f_h(\bm{s}_j)$. Let $\bm{a}(\bm{s}, \bm{s}_j)$ contain the values $a_h(\bm{s}, \bm{s}_j)$ for each of the $h = 1, \dots, \tilde{H}$ spatially-varying coefficients. Using this notation, the full conditional for $\textbf{w}(\bm{s}_j)$ is

\begin{equation}\label{Eq: w} \begin{array}{c} \textbf{w}(\bm{s}j) \mid \cdot N{\tilde{H}}(\bm{\mu}j\bm{\Sigma}_j, \bm{\Sigma}_j) \mbox{ where, }\ \bm{\mu}_j = \bm{F}(\bm{s}_j)^{-1}\bm{B}(\bm{s}_j)\textbf{w}(N(\bm{s}_j)) + \sum{\bm{s} \in U(\bm{s}j)}\bm{B}(\bm{s}, \bm{s}_j)^\top\bm{F}(\bm{s}_j)^{-1}\bm{a}(\bm{s}, \bm{s}_j) + \ \tilde{\textbf{x}}(\bm{s}_j)\omega{\beta}(\bm{s})j((z(\bm{s}_j) - 0.5)/\omega{\beta}(\bm{s}j) - \textbf{x}(\bm{s}_j)^\top\bm{\beta}) \mbox{ and } \ \bm{\Sigma}_j = (\bm{F}(\bm{s}_j)^{-1} + \sum{\bm{s} \in U(\bm{s}j)} \bm{B}(\bm{s}, \bm{s}_j)^\top\bm{F}(\bm{s}_j)^{-1}\bm{B}(\bm{s}, \bm{s}_j) + \tilde{\textbf{x}}(\bm{s}_j)^\top\bm{S}{j, \beta}\tilde{\textbf{x}}(\bm{s})_j)^{-1}, \end{array} \end{equation}

where $\textbf{w}(N(\bm{s}j))$ is a stacked $mq \times 1$ vector of the $m$ realizations of each of the $h$ NNGPs at the locations in $N(\bm{s}_j)$ and $\bm{S}{j, \beta}$ is a $\tilde{H} \times \tilde{H}$ diagonal matrix with each diagonal entry equal to the \pg auxiliary variable at site $j$.

Update spatial decay parameters ($\bm{\phi}$)

We use a Metropolis within Gibbs step to sample $\bm{\phi}$, the spatial decay parameters for each spatially-varying coefficient. We assign uniform priors for all spatial decay parameters The full conditional posterior density for $\phi_h$ for each $h = 1, \dots, \tilde{H}$ is proportional to \begin{equation}\label{Eq: fulltheta} \begin{array}{c} p(\phi_h \mid \cdot ) \propto p_h(\phi_h)p(\textbf{w}h \mid \phi_h) \ \propto p(\phi_h) \times \prod{j = 1}^J N\left(\text{w}_h(\bm{s}_j) \mid \textbf{b}_h(\bm{s}_j)\textbf{w}_h (N(\bm{s}_j)), f_h (\bm{s}_j). \right) \end{array} \end{equation}

The same update is used for the spatial smoothness parameter $\nu_h$ if using a Matérn correlation function.

Update spatial variances ($\bm{\sigma^2}$)

The default prior for the spatial variance parameter $\sigma^2_h$ for each spatially-varying coefficient in \texttt{spOccupancy} is an inverse-Gamma distribution (i.e., $\sigma^2_h \sim IG(a_{\sigma}, b_{\sigma})$). When using an inverse-Gamma prior, we sample from the following full conditional:

\begin{equation}\label{eq:sigmaSq} \sigma^2_h \mid \cdot \sim IG(a_\sigma+J/2,\; b_\sigma+\frac 12 \sum_{j=1}^J \left(\text{w}_h(\bm{s}_j)-\textbf{b}_h(\bm{s}_j)\textbf{w}_h(N(\bm{s}_j)) \right) ^2/f^\ast_h(\bm{s}_j), \end{equation}

where $f^\ast_h(\bm{s}_j) = f_h(\bm{s}_j) / \sigma^2_h$. Note that to generate $\sigma^2_h$ for the $l$th iteration of the Gibbs sampler, we calculate $f^\ast_h(\bm{s}_j)$ using the value of $\sigma^2_h$ in the $(h - 1)$th iteration.

If using a uniform prior on the spatial variances, the update is analogous to the update for the Metropolis within Gibbs step of the spatial decay parameters in Equation \ref{Eq: fulltheta}.

Update latent occurrence values ($\bm{z}(\bm{s}_j)$)

Finally, we sample the latent occurrence states. We set $z(\bm{s}_j) = 1$ for all sites where there is at least one detection of our species of interest, and so we only need to sample $z(\bm{s}_j)$ at sites where there are no detections. Thus, for all locations with no detections, we sample $z(\bm{s}_j)$ according to

\begin{equation} z(\bm{s}j) \mid \cdot \sim \text{Bernoulli}\Bigg(\frac{\psi(\bm{s}_j) \prod{k = 1}^{K_j}(1 - p_{k}(\bm{s}j))}{1 - \psi(\bm{s}_j) + \psi(\bm{s}_j) \prod{k = 1}^{K_j}(1 - p_{k}(\bm{s}_j))}\Bigg). \end{equation}

Multi-species spatially-varying coefficient occupancy model

Model description

Now consider the case where there are multiple species of interest, $N$, that are observed during data collection. We seek to jointly model the occupancy of the $N$ species in a single model that accommodates residual correlations between species and allows for sharing of information across species via hierarchical effects. Such multi-species approaches often have improved precision and accuracy of estimates compared to single-species models [@clark2014more; @zipkin2010multi]. Given the interest in modeling large-scale patterns across space, a subset of spatial locations in $\mathcal{L}$ may fall well outside the potential range of a given species. In our subsequent modeling development, we consider the case where species-specific range maps are available as an auxiliary data source to more readily accommodate large differences in the potential locations where each species can occur.

Using similar notation to the single-species models, we model the true presence-absence state of species $i$ at site $\bm{s}_j$ following

\begin{equation}\label{zMS} z_i(\bm{s}_j) \sim \text{Bernoulli}(\psi_i(\bm{s}_j) \cdot z^\ast_i(\bm{s}_j)), \end{equation}

where $\psi_i(\bm{s}_j)$ is the occupancy probability of species $i$ at site $j$, and $z^\ast_i(\bm{s}_j)$ is a binary auxiliary data source indicating whether site $j$ is within the range of species $i$. Such data can be obtained from a variety of sources, including international databases (e.g., BirdLife International, IUCN) or expert opinion. We suggest buffering the auxiliary data range map by a suitable distance to account for potential inaccuracies in the auxiliary data. Inclusion of such auxiliary range data can drastically reduce the computational burden of the model if certain species can only exist at a subset of the spatial locations in $\mathcal{L}$ [@socolar2022biogeographic]. If auxiliary range data are not available, $z^\ast_i(\bm{s}_j)$ can be removed from (\ref{zMS}) (or equivalently, $z^\ast_i(\bm{s}_j) = 1$ for all $j$).

At sites within each species set of potentially suitable areas, we model species-specific occupancy probability as

\begin{equation}\label{psiMS} \text{logit}(\psi_i(\bm{s}j)) = (\beta{i, 1} + \delta_1\text{w}^\ast_{i, 1}(\bm{s}j)) + \sum{h = 2}^H\text{x}h(\bm{s}_j){\beta{i, h} + \delta_j\text{w}^\ast_{i, h}(\bm{s}_j)}, \end{equation}

with $\text{w}^\ast_{i, h}(\bm{s}j)$ being the spatially-varying component of the intercept ($h = 1$) or coefficient ($h > 1$) for species $i$ at site $\bm{s}_j$ and all parameters as defined before, but now spatial and non-spatial effects are unique to each species. Note we assume the same variables have spatially-varying effects across species (i.e., $\delta{i, h} = \delta_h$ for all $i$), although this could be modified to allow $\delta_h$ to vary by species by either setting values \textit{a priori} or estimating the indicator variables within the model itself. We model the non-spatial component of the $h$th regression coefficient for each species $i$ hierarchically from a common community-level distribution to share information across species [@dorazio2005; @gelfand2005modelling]. More specifically, we have

\begin{equation}\label{betaMS} \beta_{i, h} \sim \text{Normal}(\mu_{\beta_h}, \tau^2_{\beta_h}), \end{equation} where $\mu_{\beta_h}$ is the average non-spatial effect across all species in the community, and $\tau^2_{\beta_h}$ is the variability in the non-spatial effect across all species.

We seek to jointly model the species-specific spatial effects, $\textbf{w}^\ast_{i, t}$, to account for correlation in species-specific responses to covariate effects, as well as residual correlation between species after accounting for their relationships with any covariates included in the model. For a small number of species (e.g., 5), a linear model of coregionalization (LMC) framework [@gelfand2004] is a viable solution, but such an approach quickly becomes computationally intractable as the number of species in the community increases. Instead, we use a spatial factor model [@hogan2004bayesian], a dimension reduction technique that accounts for correlations in species-specific responses, while drastically reducing computational run time compared to a LMC that requires estimation of a full $N \times N$ cross-covariance matrix for each effect that is assumed to vary spatially. Here, we decompose $\text{w}^\ast_{i, h}(\bm{s}_j)$ into a linear combination of $q$ latent factors and their associated species-specific coefficients (i.e., factor loadings). Thus for each SVC in the model, we have

\begin{equation}\label{wStar} \text{w}^\ast_{i, h}(\bm{s}j) = \bm{\lambda}^\top{i, h}\textbf{w}{h}(\bm{s}_j), \end{equation} where $\bm{\lambda}{i, h}^\top$ is the $i$th row of factor loadings from the $N \times q$ loadings matrix $\bm{\Lambda}_h$, and $\textbf{w}_h(\bm{s}_j)$ is a $q \times 1$ vector of independent spatial factors at site $j$. As in the single-species SVC occupancy model, we model each of the $r = 1, \dots, q$ spatial factors for each of the $\tilde{H}$ spatially-varying effects with an NNGP prior following

\begin{equation}\label{NNGPMS} \tilde{p}(\textbf{w}{r, h}) = \prod{j = 1}^J \text{Normal}(\text{w}{r, h}(\bm{s}_j) \mid \textbf{b}{r, h}(\bm{s}j)\textbf{w}{r, h}(N(\bm{s}j)), \text{f}{r, h}(\bm{s}j)), \end{equation} with $\textbf{b}{r, h}(\bm{s}j)$ and $\text{f}{r, h}(\bm{s}_j)$ defined in (\ref{bNNGP}) and (\ref{fNNGP}), respectively.

For each SVC, we can derive an inter-species covariance matrix $\bm{\Sigma}_h = \bm{\Lambda}_h\bm{\Lambda}_h^\top$, which has rank $q << N$, and, thus, is singular. However, the inter-species covariance matrices can still be used to detect species clustering [@shirota2019sinica; @doser2023joint]. For a spatially-varying intercept, the inter-species covariance matrix provides information on the residual co-occurrence patterns between each pair of species across space after accounting for any covariates included in the model, which can be used to generate hypotheses regarding the abiotic and/or biotic drivers of residual species co-occurrence patterns (e.g., @tobler2019joint). For a spatially-varying covariate effect, a positive correlation between species in $\bm{\Sigma}_h$ indicates similar responses to an environmental covariate across space, which can be used as a model-based ordination technique to identify groups of species that respond similarly to changes in environmental variables.

Analogous to the single-species case, we model the observed detection-nondetection of each species $i$ at site $j$ during replicate survey $k$, $y_{i, k}(\bm{s}j)$ conditional on the true presence-absence of each species, $z{i}(\bm{s}_j)$, following (\ref{yDet}) and (\ref{pDet}), with all parameters now varying by species. We model the species-specific detection regression coefficients ($\bm{\alpha}_i$) hierarchically, analogous to the non-spatial occupancy regression coefficients in (\ref{betaMS}).

We assume Gaussian priors for all mean parameters and inverse-Gamma priors for variance parameters. Additional restrictions on the factors loadings matrix $\bm{\Lambda}h$ for each spatially-varying coefficient $h$ are required to ensure identifiability [@taylor2019spatial]. We fix all elements in the upper triangle to 0 and set the diagonal elements to 1. We additionally fix the spatial variance parameters $\sigma^2_h$ of each latent spatial process to 1. We assign standard Gaussian priors for the lower triangular elements in $\bm{\Lambda}$ and assign each spatial range parameter $\phi{r, h}$ an independent uniform prior.

Gibbs sampler {#sfMsPGOccSampler}

Here we describe the Gibbs sampler for fitting the spatially-varying coefficient multi-species occupancy model using svcMsPGOcc().

Update community-level occurrence coefficients ($\bm{\mu_{\beta}}$)

We first sample all community-level parameters followed by species level parameters. First we sample the community-level occurrence coefficients. Let $\bm{\mu}{\beta}$ denote the vector of all community-level occurrence means, and similarly let $\bm{T}{\beta}$ denote the variance matrix of all community-level occurrence variance parameters. Note that $\bm{T}{\beta}$ is a diagonal matrix. Let $\bm{\mu}{\beta} \sim N(\bm{\mu}{0, \beta}, \bm{\Sigma}{\beta})$ denote our prior distribution, where $\bm{\Sigma}{\beta}$ is a diagonal matrix. Note this is equivalent to assigning an independent normal prior for each coefficient. Our full conditional for the community-level regression coefficients $\bm{\mu{\beta}}$ is then

\begin{equation} \bm{\mu_{\beta}} \mid \cdot \sim N([\bm{\Sigma}{\beta}^{-1} + N\bm{T}{\beta}^{-1}]^{-1}\Big[\sum_{i = 1}^N(\bm{T}{\beta}^{-1}\bm{\beta}_i) + \bm{\Sigma}{\beta}^{-1}\bm{\mu}{0, \beta}\Big], [\bm{\Sigma}{\beta}^{-1} + N\bm{T}_{\beta}^{-1}]^{-1}). (#eq:muBeta) \end{equation}

Update community-level detection coefficients ($\bm{\mu_{\alpha}}$)

Next, we sample the community-level detection coefficients. Let $\bm{\mu}{\alpha}$ denote the vector of all community-level detection means, and similarly let $\bm{T}{\alpha}$ denote the diagonal variance matrix of all community-level detection variance parameters. Let $\bm{\mu}{\alpha} \sim N(\bm{\mu}{0, \alpha}, \bm{\Sigma}{\alpha})$ denote the prior distribution, where $\bm{\Sigma}{\alpha}$ is a diagonal matrix. Our full conditional then takes the form

\begin{equation} \bm{\mu_{\alpha}} \mid \cdot \sim N([\bm{\Sigma}{\alpha}^{-1} + N\bm{T}{\alpha}^{-1}]^{-1}\Big[\sum_{i = 1}^N(\bm{T}{\alpha}^{-1}\bm{\alpha}_i) + \bm{\Sigma}{\alpha}^{-1}\bm{\mu}{0, \alpha}\Big], [\bm{\Sigma}{\alpha}^{-1} + N\bm{T}_{\alpha}^{-1}]^{-1}). (#eq:muAlpha) \end{equation}

Update community-level occurrence variances ($\bm{\tau^2_{\beta}}$)

Let $\tau^2_{h, \beta}$ denote the community-level variance for the $h^{\text{th}}$ occurrence parameter ($h = 1, \dots, p_{\psi}$). We assign an inverse gamma prior to $\tau^2_{h, \beta}$ with shape parameter $a_{\tau_{h, \beta}}$ and scale parameter $b_{\tau_{h, \beta}}$. Our full conditional is then

\begin{equation} \tau^2_{h, \beta} \mid \cdot \sim \text{IG}(a_{\tau_{h, \beta}} + \frac{N}{2}, b_{\tau_{h, \beta}} + \frac{\sum_{i = 1}^N(\beta_{i, h} - \mu_{\beta_h})^2}{2}). (#eq:tauBeta) \end{equation}

Update community-level detection variances ($\bm{\tau^2_{\alpha}}$)

Let $\tau^2_{h, \alpha}$ denote the community-level variance for the $h^{\text{th}}$ detection parameter ($h = 1, \dots, p_{H}$). We assign an inverse gamma prior to $\tau^2_{h, \alpha}$ with shape parameter $a_{\tau_{h, \alpha}}$ and scale parameter $b_{\tau_{h, \alpha}}$. Our full conditional is then

\begin{equation} \tau^2_{h, \alpha} \mid \cdot \sim \text{IG}(a_{\tau_{h, \alpha}} + \frac{N}{2}, b_{\tau_{h, \alpha}} + \frac{\sum_{i = 1}^N(\alpha_{i, h} - \mu_{\alpha_h})^2}{2}). (#eq:tauAlpha) \end{equation}

Update species-specific occurrence auxiliary variables ($\omega_{i, \beta}(\bm{s}_j)$)

We next sample the occurrence auxiliary variable ($\omega_{i, \beta}(\bm{s}_j)$) individually for each species $i$ and site $j$. If $z_i^\ast(\bm{s}_j) = 1$ (i.e., site $\bm{s}_j$ is within the species range as defined by any auxiliary data source supplied to the model), our full conditional is

\begin{equation} \omega_{i, \beta}(\bm{s}j) \mid \cdot \sim \text{PG}(1, \bm{x}(\bm{s}_j)^{\top}\bm{\beta}_i + \tilde{\textbf{x}}(\bm{s}_j)^\top\text{w}^*{i}(\bm{s}_j)). (#eq:omegaBeta) \end{equation}

For sites and species where $z_i^\ast(\bm{s}_j) = 0$, we do not need to generate an auxiliary variable.

Update species-specific detection auxiliary variables ($\omega_{i, k, \alpha}(\bm{s}_j)$)

We next update the latent $\pg$ auxiliary variable for the detection process, $\omega_{i, k, \alpha}(\bm{s}j)$, for each replicate $k$ at each site $j$ for each species $i$. Note that we only need to sample $\omega{i, k, \alpha}(\bm{s}j)$ when $z_i(\bm{s}_j) = 1$, which can change across different MCMC iterations. For such cases when $z_i(\bm{s}_j) = 0$, set $\omega{i, k, \alpha}(\bm{s}_j) = 0$. Following @polson2013, for all $\bm{s}_j$ with $z_i(\bm{s}_j) = 1$ at the current iteration, we have

\begin{equation} \omega_{i, k, \alpha}(\bm{s}_j) \mid \cdot \sim \text{PG}(1, \bm{v}(\bm{s}_j)^\top\bm{\alpha}_i). \end{equation}

Update species-level occurrence regression coefficients ($\bm{\beta}_i$)

We update the species-level occurrence regression coefficients ($\bm{\beta}_i$), including the intercept, from the following multivariate normal full conditional

\begin{equation} \bm{\beta}i \mid \cdot \sim \text{Normal}(\textbf{L}^\ast{\beta, i}\bm{\mu}^\ast_{\beta, i}, \textbf{L}^\ast_{\beta, i}), (#eq:beta) \end{equation}

where

\begin{equation} \bm{\mu}^\ast_{\beta, i} = \bm{X}^{\top}\bm{Z}^\ast_i(\bm{z}i - 0.5 \bm{1}_J - \bm{S}{\beta}\sum_{h = 1}^{\tilde{H}}\tilde{\bm{X}}{h, J}\bm{Z}^\ast_i\textbf{W}_h^\top\bm{\lambda}{i, h}) + \bm{T}{\beta}^{-1}\bm{\mu{\beta}}, \end{equation}

and

\begin{equation} \textbf{L}^\ast_{\beta, i} = [\bm{T}{\beta}^{-1} + \bm{X}^{\top}\bm{Z}^\ast_i\bm{S}{\beta}\bm{Z}^\ast_i\bm{X}]^{-1}, \end{equation}

where $\bm{S}{\beta}$ is a diagonal $J \times J$ matrix with diagonal entries equal to the latent $\pg$ variable values for species $i$, $\bm{z}_i$ is the $J \times 1$ vector of latent occurrence values for species $i$, $\bm{1}_J$ is a $J \times 1$ vector of 1s, $\bm{Z}^\ast_i$ is a $J \times J$ diagonal matrix with diagonal elements equal to the auxiliary range binary data points for species $i$, $\tilde{\bm{X}}{h, J}$ is a $J \times J$ diagonal matrix with diagonal elements equal to the $h$th covariate value at site $j$ whose effect is assumed to vary spatially, $\textbf{W}h$ is the $q \times J$ matrix of spatial factors for the $h$th spatially-varying coefficient, and $\bm{\lambda}{i, h}$ is the $i$th row of factor loadings for the $t$th spatially-varying coefficient.

Update species-level detection regression coefficients ($\bm{\alpha}_i$)

Next, we sample the species-specific detection regression coefficients for species $i$ ($\bm{\alpha}i$). Define $K{\text{TOT}} = \sum_{j = 1}^JK_j$. Our full conditional update takes the form

\begin{equation} \bm{\alpha}i \mid \cdot \sim \text{Normal}(\textbf{L}^\ast{\alpha, i}\bm{\mu}^\ast_{\alpha, i}, \textbf{L}^\ast_{\alpha, i}), (#eq:alpha) \end{equation}

\begin{equation} \bm{\mu}^\ast_{\alpha, i} = \bm{V}^{\top}\bm{Z}i(\bm{y}_i - 0.5 \bm{1}{K_{\text{TOT}}}) + \bm{T}{\alpha}^{-1}\bm{\mu{\alpha}}, \end{equation}

and

\begin{equation} \textbf{L}^\ast_{\alpha, i} = [\bm{T}{\alpha}^{-1} + \bm{V}^{\top}\bm{Z}_i\bm{S}{\alpha}\bm{Z}_i\bm{V}]^{-1}, \end{equation}

where $\bm{S}\alpha$ is a diagonal $K{\text{TOT}} \times K_{\text{TOT}}$ matrix with diagnoal entries equal to the latent \pg variable values for species $i$, $\bm{y}i$ is the $K{\text{TOT}} \times 1$ vector of detection-nondetection values for species $i$, $\bm{1}{K{\text{TOT}}}$ is a $K_{\text{TOT}} \times 1$ vector of 1s, and $\bm{Z}i$ is a $K{\text{TOT}} \times K_{\text{TOT}}$ diagonal matrix with diagonal elements equal to the current value of the latent occurrence state $\bm{z}_i$ at the site corresponding to the given element of $\bm{y}_i$.

Update latent spatial factors for each SVC ($\textbf{w}_h(\bm{s}_j)$)

Let $N(\bm{s}j)$ denote the set of $m$ nearest neighbors of $\bm{s}_j$ among $\bm{s}_1, \bm{s}_2, \dots, \bm{s}{j - 1}$. Let $\textbf{w}{r, h}(N(\bm{s}_j))$ denote the $m$ realizations of the $r^{\text{th}}$ NNGP factor for the $h$th spatially-varying coefficient at the locations in $N(\bm{s}_j)$. Let $C(\cdot, \phi{r, h})$ denote the correlation function of the original Gaussian Process (GP) from which the NNGP is derived. For any two sets $A_1$ and $A_2$, define $\text{C}{A_1, A_2}(\phi{r, h})$ as the correlation matrix between the observations in $A_1$ and $A_2$ for the $r^{\text{th}}$ GP for the $h^{\text{th}}$ spatially-varying coefficient. For $j \geq 1$, we have

\begin{equation} \textbf{b}{r, h}(\bm{s}_j) = \textbf{C}{\bm{s}j, N(\bm{s}_j)}(\phi{r, h})\textbf{C}^{-1}{N(\bm{s}_j), N(\bm{s}_j)}(\phi{r, h}), (#eq:bNNGP) \end{equation}

where $\textbf{b}_{r, h}(\bm{s}_1) = \bm{0}$ for all $r = 1, \dots, q$ and $h = 1, \dots, \tilde{H}$. Further, we have

\begin{equation} f_{r, h}(\bm{s}j) = \textbf{C}{\bm{s}j, \bm{s}_j}(\phi{r, h}) - \textbf{C}{\bm{s}_j, N(\bm{s}_j)}(\phi{r, h})\textbf{C}^{-1}{N(\bm{s}_j), N(\bm{s}_j)}(\phi{r, h})\textbf{C}{N(\bm{s}_j), \bm{s}_j}(\phi{r, h}). (#eq:fNNGP) \end{equation}

For any two locations $\bm{s}1$ and $\bm{s}_2$, if $\bm{s}_1 \in N(\bm{s}_2)$ and is the $l^{\text{th}}$ member of $N(\bm{s}_2)$, then define $b{r, h}(\bm{s}2, \bm{s}_1)$ as the $l^{\text{th}}$ entry of $\textbf{b}{r, h}(\bm{s}2)$. Let $U(\bm{s}_1) = {\bm{s}_2 \in S \mid \bm{s}_1 \in N(\bm{s}_2)}$ be the collection of locations $\bm{s}_2$ for which $\bm{s}_1$ is a neighbor, where $S$ is the set of all $J$ spatial locations. For every $\bm{s}_2 \in U(\bm{s}_1)$, define $a{r, t}(\bm{s}2, \bm{s}_1) = \text{w}{r, h}(\bm{s}2) - \sum{\bm{s} \in N(\bm{s}2), \bm{s} \neq \bm{s}_2} \text{w}{r, h}(\bm{s})b_{r, h}(\bm{s}2, \bm{s})$. Extending this to matrix notation, let $\bm{B}_h(\bm{s}_j)$ be a $q \times mq$ block matrix, with each $q \times q$ diagonal block containing the elements of $\bm{b}{r, h}(\bm{s}j)$ for each of the $r = 1, \dots q$ spatial factors for the $h = 1, \dots, \tilde{H}$ spatially-varying coefficients for each of the specific $m$ neighbors. Let $\bm{F}_h(\bm{s}_j)$ be a $q \times q$ diagonal matrix with diagonal elements of $f{r, h}(\bm{s}j)$. Let $\bm{a}_h(\bm{s}, \bm{s}_j)$ contain the values $a{r, h}(\bm{s}, \bm{s}_j)$ for each of the $r = 1, \dots, q$ latent factors for the $h$th spatially-varying coefficient. Using this notation, the full conditional distribution for $\textbf{w}_h(\bm{s}_j)$ is

\begin{equation} \begin{array}{c} \textbf{w}h(\bm{s}_j) \mid \cdot \sim N_q(\bm{\mu}{j, h}\bm{\Sigma}{j, h}, \bm{\Sigma}{j, h}) \mbox{ where, }\ \bm{\mu}{j, h} = \bm{F}_h(\bm{s}_j)^{-1}\bm{B}_h(\bm{s}_j)\textbf{w}_h(N(\bm{s}_j)) + \sum{\bm{s} \in U(\bm{s}j)}\bm{B}_h(\bm{s}, \bm{s}_j)^\top\bm{F}_h(\bm{s}_j)^{-1}\bm{a}_h(\bm{s}, \bm{s}_j) + \ \bm{\Lambda}_h^\top\tilde{\bm{X}}{h, N}(\bm{s}j)\bm{Z}^\ast{N}(\bm{s}j)\bm{S}{N}(\bm{s}j)(\bm{S}{N}(\bm{s}j)^{-1}(\bm{z}(\bm{s}_j) - 0.5\bm{1}_N) - \bm{Z}^\ast_N(\bm{s}_j)(\bm{X}(\bm{s}_j)^\top\bm{\beta} + \sum{k \neq h}^{\tilde{H}}\tilde{\bm{X}}{k, N}(\bm{s}_j)\bm{\Lambda}_k\textbf{w}_k(\bm{s}_j))) \mbox{ and } \ \bm{\Sigma}{j, h} = (\bm{F}h(\bm{s}_j)^{-1} + \sum{\bm{s} \in U(\bm{s}j)} \bm{B}_h(\bm{s}, \bm{s}_j)^\top\bm{F}_h(\bm{s}_j)^{-1}\bm{B}_h(\bm{s}, \bm{s}_j) + \ (\bm{\Lambda}_h^\top\tilde{\bm{X}}{h, N}(\bm{s}j)\bm{Z}^\ast_N(\bm{s}_j))\bm{S}{N}(\bm{s}j)(\bm{\Lambda}_h^\top\tilde{\bm{X}}{h, N}(\bm{s}_j)\bm{Z}^\ast_N(\bm{s}_j))^\top)^{-1}, \end{array} (#eq:w) \end{equation}

where $\textbf{w}h(N(\bm{s}_j))$ is a stacked $mq \times 1$ vector of the $m$ realizations of each of the $r$ NNGPs at the locations in $N(\bm{s}_j)$, $\bm{S}{N}(\bm{s}j)$ is an $N \times N$ diagonal matrix with the $\pg$ occupancy variables for each species $i$ at site $j$ along the diagonal elements, $\tilde{\bm{X}}{h, N}(\bm{s}_j)$ is a $N \times N$ diagonal matrix with diagonal elements all equal to the value of the $h$th spatially-varying covariate at location $\bm{s}_j$, $\bm{Z}^\ast_N(\bm{s}_j)$ is a $N \times N$ diagonal matrix with diagonal elements equal to the auxiliary range binary data points at site $\bm{s}_j$ for each of the $N$ species, $\bm{X}(\bm{s}_j)^\top$ is a $N \times (NH)$ block-diagonal matrix with the $i$th diagonal block the $\bm{x}(\bm{s}_j)$ vector of $H$ covariates (including the intercept), and $\bm{\beta}$ is the $(NH) \times 1$ stacked vector of species-specific regression coefficients (including the intercept).

Update latent spatial factor loadings for each SVC ($\bm{\Lambda}_h$)

Recall we set all diagonal elements of $\bm{\Lambda}h$ to 1 and all upper triangular elements equal to 0 in order to ensure identifiability of the latent spatial factors. Given this requirement, let $q_i = \text{min}{i - 1, q}$ for $2 \leq i \leq N$, and let $\tilde{\bm{\lambda}}{i, h} = (\lambda_{i, h, 1}, \dots, \lambda_{i, h, q_i})^\top$ be the vector representing the unrestricted elements in the $i^{\text{th}}$ row of $\bm{\Lambda}h$. Define $\textbf{W}_h$ as the $J \times q$ matrix of latent spatial factors for the $h$th SVC, and let $\textbf{W}{1:i, h}$ be the first $i$ columns of $\textbf{W}h$. Using this notation, the full conditional density for $\tilde{\bm{\lambda}}{i, h}$ is $N_q(\bm{\Omega}{\tilde{\bm{\lambda}}{i, h}}\bm{\mu}{\tilde{\bm{\lambda}}{i, h}}, \bm{\Omega}{\tilde{\bm{\lambda}}{i, h}})$, where

\begin{align} \bm{\mu}{\tilde{\bm{\lambda}}{i, h}} &= \left{ \begin{matrix} (\tilde{\bm{X}}{h, J}\bm{Z}^\ast_i\textbf{W}{h, 1:(i-1)})^\top\bm{S}{i, \beta}(\bm{S}{i, \beta}^{-1}(\bm{z}i - 0.5\bm{1}_J) - \bm{Z}^\ast_i(\bm{X}\bm{\beta}_i + [\sum{k \neq h}^{\tilde{H}}\tilde{\bm{X}}{k, J}\textbf{W}_k\bm{\lambda}{k, i}] + \tilde{\bm{X}}{h, J}\dot{\textbf{w}}{h, i}))\hfill &\text{ if }&2\leq i\leq q \ (\tilde{\bm{X}}{h, J}\bm{Z}^\ast_i\textbf{W}{h})^\top\bm{S}{i, \beta}(\bm{S}{i, \beta}^{-1}(\bm{z}i - 0.5\bm{1}_J) - \bm{Z}^\ast_i(\bm{X}\bm{\beta}_i + [\sum{k \neq h}^{\tilde{H}}\tilde{\bm{X}}{k, J}\textbf{W}_k\bm{\lambda}{k, i}])) \hfill &\text{ if }&i > q \end{matrix}\right.\ \bm{\Omega}{\tilde{\bm{\lambda}}{i, h}} &= \left{ \begin{matrix} ((\tilde{\bm{X}}{h, J}\bm{Z}^\ast_i\textbf{W}{h, 1:(i-1)})^\top\bm{S}{i, \beta}(\tilde{X}{h, J}\bm{Z}^\ast_i\textbf{W}{h, 1:(i - 1)}) + I{i - 1})^{-1} \hfill &\text{ if }&2\leq i \leq q \ ((\tilde{\bm{X}}{h, J}\bm{Z}^\ast_i\textbf{W}{h})^\top\bm{S}{i, \beta}(\tilde{X}{h, J}\bm{Z}^\ast_i\textbf{W}{h}) + I{q})^{-1} \hfill &\text{ if }&i > q \end{matrix}\right., \end{align}

where $\dot{\textbf{w}}_{h, i}$ is the $i^{\text{th}}$ column of $\textbf{W}_h$, and all other variables are as defined previously.

Update spatial range parameters ($\bm{\phi}$)

We use a Metropolis within Gibbs step to sample $\bm{\phi}$. The full conditional posterior density for $\phi_{r, h}$ for each $r = 1, \dots, q$ and $h = 1, \dots, \tilde{H}$ is proportional to \begin{equation} \begin{array}{c} p(\phi_{r, h} \mid \cdot ) \propto p(\phi_{r, h})p(\textbf{w}{r, h} \mid \phi{r, h}) \ \propto p(\phi_{r, h}) \times \prod_{j = 1}^J N\left(\text{w}{r, h}(\bm{s}_j) \mid \textbf{b}{r, h}(\bm{s}j)\textbf{w}{r, h} (N(\bm{s}j)), f{r, h} (\bm{s}_j). \right) \end{array} (#eq:fullTheta) \end{equation}

We sample $\phi_{r, h}$ using a random walk Metropolis step. We use a normal proposal distribution along with a Jacobian transformation.

Update latent occurrence values ($z_{i}(\bm{s}_j)$)

Finally, we sample the latent occurrence states for each species. We set $z_{i}(\bm{s}j) = 1$ for all sites where there is at least one detection of species $i$, and so we only need to sample $z{i}(\bm{s}j)$ at sites where there are no detections. We set $z_i(\bm{s}_j) = 0$ for all sites where $z^\ast_i(\bm{s}_j) = 0$. Thus, for all locations with no detections of the species $i$ where $z^\ast_i(\bm{s}_j) = 1$, we sample $z{i}(\bm{s}_j)$ according to

\begin{equation} z_{i}(\bm{s}j) \mid \cdot \sim \text{Bernoulli}\Bigg(\frac{\psi{i}(\bm{s}j) \prod{k = 1}^{K_j}(1 - p_{i, k}(\bm{s}j))}{1 - \psi{i}(\bm{s}j) + \psi{i}(\bm{s}j) \prod{k = 1}^{K_j}(1 - p_{i, k}(\bm{s}_j))}\Bigg). (#eq:z) \end{equation}

References {-}

doserjef/spOccupancy documentation built on Feb. 28, 2025, 1:19 p.m.