Fitting Nmixture models in spAbundance
Jeffrey W. Doser
2023
Source:vignettes/nMixtureModels.Rmd
nMixtureModels.Rmd
Introduction
This vignette provides worked examples and explanations for fitting
singlespecies and multispecies Nmixture models in the
spAbundance
R package. We will provide step by step
examples on how to fit the following models:
 Nmixture model using
NMix()
.  Spatial Nmixture model using
spNMix()
.  Multispecies Nmixture model using
msNMix()
.  Multispecies Nmixture model with species correlations using
lfMsNMix()
.  Spatial multispecies Nmixture with species correlations using
sfMsNMix()
.
In this vignette we are only describing spAbundance
functionality to fit Nmixture models, with separate vignettes on
fitting hierarchical distance sampling models and generalized linear
mixed models. We fit all models in a Bayesian framework using custom
Markov chain Monte Carlo (MCMC) samplers written in C/C++
and called through R
’s foreign language interface. Here we
will provide a brief description of each model, with full statistical
details provided in a separate vignette. As with all model types in
spAbundance
, we will show how to perform posterior
predictive checks as a GoodnessofFit assessment, model comparison
using the Widely Applicable Information Criterion (WAIC), and
outofsample predictions using standard R helper functions (e.g.,
predict()
). Note that syntax of Nmixture models in
spAbundance
closely follows syntax for fitting occupancy
models in spOccupancy
(Doser et al.
2022), and that this vignette closely follows the documentation
on the spOccupancy
website.
To get started, we load the spAbundance
package, as well
as the coda
package, which we will use for some MCMC
summary and diagnostics. We will also use the stars
and
ggplot2
packages to create some basic plots of our results.
We then set a seed so you can reproduce the same results as we do.
Example data set: Simulated multispecies count data
As an example data set throughout this vignette, we will use a
simulated data set comprised of 6 species across 225 sites and a maximum
of 3 repeat surveys per site. Note that every site was not sampled for
all three surveys, and thus our data set is “imbalanced”. The data are
provided as part of the spAbundance
package and are loaded
with data(dataNMixSim)
. The manual page obtained using
help(dataNMixSim)
or ?dataNMixSim
contains
information on the parameter values used to generate the data set with
the simMsNMix()
function. You can look them up to have a
benchmark for comparison with the parameter estimates that we will
obtain very soon.
List of 4
$ y : int [1:6, 1:225, 1:3] 1 0 1 0 0 0 NA NA NA NA ...
$ abund.covs: num [1:225, 1:2] 0.373 0.706 0.202 1.588 0.138 ...
.. attr(*, "dimnames")=List of 2
.. ..$ : NULL
.. ..$ : chr [1:2] "abund.cov.1" "abund.factor.1"
$ det.covs :List of 2
..$ det.cov.1: num [1:225, 1:3] 1.28 NA NA NA 1.04 ...
..$ det.cov.2: num [1:225, 1:3] 2.03 NA NA NA 0.796 ...
$ coords : num [1:225, 1:2] 0 0.0714 0.1429 0.2143 0.2857 ...
.. attr(*, "dimnames")=List of 2
.. ..$ : NULL
.. ..$ : chr [1:2] "X" "Y"
The object dataNMixSim
is a list comprised of the count
data (y
), covariates on the abundance portion of the model
(abund.covs
), covariates on the detection portion of the
model (det.covs
), and the spatial coordinates of each site
(coords
) for use in spatial Nmixture models and in
plotting. This list is in the exact format required for input to
Nmixture models in spAbundance
. dataNMixSim
contains data on 6 species in the threedimensional array
y
, where the dimensions of y
correspond to
species (6), sites (225), and maximum number of replicates at any given
site (3). For singlespecies Nmixture models in Section 2 and 3, we
will only use data on one species, so we next subset the
dataNmixSim
list to only include data from the first
species in a new object data.one.sp
.
data.one.sp < dataNMixSim
data.one.sp$y < data.one.sp$y[1, , ]
table(data.one.sp$y) # Quick summary.
0 1 2 3 4 5 7 13
276 102 40 8 4 1 2 1
We see that our species is fairly rare, with most observations of the species being 0.
Singlespecies Nmixture models
Basic model description
Let \(N_j\) denote the true abundance of a species of interest at site \(j = 1, \dots, J\). We model \(N_j\) following either a Poisson or negative binomial (NB) distribution according to
\[\begin{align}\label{abundance} \begin{split} N_j &\sim \text{Poisson}(\mu_j) \text{, or, } \\ N_j &\sim \text{NB}(\mu_j, \kappa), \end{split} \end{align}\]
where \(\mu_j\) is the average abundance at site \(j\) and \(\kappa\) is a positive dispersion parameter. Smaller values of \(\kappa\) indicate overdispersion in the latent abundance values, while higher values indicate minimal overdispersion in abundance relative to a Poisson distribution. Note that as \(\kappa \rightarrow \infty\), a NB model “reverts” back to the simpler Poisson model. We model \(\mu_j\) using a log link function following
\[\begin{equation}\label{muAbund} \text{log}(\mu_j) = \boldsymbol{x}_j^\top\boldsymbol{\beta}, \end{equation}\]
where \(\boldsymbol{\beta}\) is a vector of regression coefficients for a set of covariates \(\boldsymbol{x}_j\) (including an intercept). Note that while not shown, unstructured random intercepts and slopes can be included in the equation for expected abundance. This may for instance be required for accommodating some sorts of “blocks”, such as when sites are nested in a number of different regions.
Following the standard Nmixture model of Royle (2004), we suppose observers count the number of individuals of the species of interest at each site \(j\) over a set of multiple surveys \(k = 1, \dots, K_j\), denoted as \(y_{j, k}\). Note the number of surveys \(K_j\) can vary by site (which is the case with our example simulated data set), but at least some sites must be surveyed more than once to ensure identifiability without making restrictive parameteric assumptions (Knape and KornerNievergelt 2015; Stoudt, Valpine, and Fithian 2023). We model \(y_{j, k}\) conditional on the true abundance of the species at site \(j\), \(N_j\), following
\[\begin{equation}\label{yNMix} y_{j, k} \sim \text{Binomial}(N_j, p_{j, k}), \end{equation}\]
where \(p_{j, k}\) is the probability of detecting an individual given it is present at the site. We model \(p_{j, k}\) using a logit link function in which we can allow detection probability to vary over space and/or surveys. More specifically, we have
\[\begin{equation} \label{pNMix} \text{logit}(p_{j, k}) = \boldsymbol{v}_{j, k}^\top\boldsymbol{\alpha}, \end{equation}\]
where \(\boldsymbol{\alpha}\) is a vector of the effects of a set of covariates \(\boldsymbol{v}_{j, k}\) (including an intercept). As for the model for abundance, spatially unstructured random effects can be specified for the intercepts or the slopes, e.g., for observer identity in a citizen science probject.
To complete the Bayesian specification of the model, we assign vague normal priors for the abundance (\(\boldsymbol{\beta}\)) and detection (\(\boldsymbol{\alpha}\)) regression coefficients. When fitting a model with a negative binomial distribution for abundance, we specify a uiform prior for the dispersion parameter \(\kappa\).
Fitting singlespecies Nmixture models with
NMix()
The NMix()
function fits singlespecies Nmixture
models. NMix()
has the following arguments:
NMix(abund.formula, det.formula, data, inits, priors, tuning,
n.batch, batch.length, accept.rate = 0.43, family = 'Poisson',
n.omp.threads = 1, verbose = TRUE, n.report = 100,
n.burn = round(.10 * n.batch * batch.length), n.thin = 1,
n.chains = 1, ...)
The first two arguments, abund.formula
and
det.formula
, use standard R model syntax to denote the
covariates to be included in the abundance and detection portions of the
model, respectively. Only the right hand side of the formulas are
included. Random intercepts and slopes can be included in both the
abundance and detection portions of the singlespecies Nmixture model
using lme4
syntax (Bates et al.
2015). For example, to include a random intercept for different
observers in the detection portion of the model, we would include
(1  observer)
in det.formula
, where
observer
indicates the specific observer for each data
point. The names of variables given in the formulas should correspond to
those found in data
, which is a list consisting of the
following tags: y
(count data), abund.covs
(abundance covariates), det.covs
(detection covariates).
y
is a sites x replicate matrix, abund.covs
is
a matrix or data frame with sitespecific covariate values, and
det.covs
is a list with each list element corresponding to
a covariate to include in the detection portion of the model. Covariates
on detection can vary by site and/or survey, and so these covariates may
be specified as a site by survey matrix for “observation”level
covariates (i.e., covariates that vary by site and survey) or as a
onedimensional vector for sitelevel covariates. Note the tag
offset
can also be specified to include an offset in the
model.
spAbundance
can handle “imbalanced” data sets (i.e.,
each site may not have the same number of repeat visits). This is the
case with our simulated example, where some sites have three repeat
visits, while others are only visited once. To accommodate such an
imbalanced design, the detectionnondetection data matrix y
should have a value of NA
in any site/visit combination
where there is not a survey. For example, let’s take a look at the first
10 rows of our detectionnondetection matrix in
data.one.sp
head(data.one.sp$y, 10)
[,1] [,2] [,3]
[1,] 1 NA NA
[2,] NA 2 2
[3,] NA 1 NA
[4,] NA 0 1
[5,] 0 0 0
[6,] 0 NA NA
[7,] 0 NA 0
[8,] 0 0 NA
[9,] 1 2 NA
[10,] NA 0 0
Sites 1, 3, and 6 have two NA
values and thus were only
sampled once. Sites 2, 4, 7, 8, 9, and 10 have one NA value and were
sampled twice, while site 5 was sampled on all three occasions. Such
imbalanced data sets are very common in practice. This same approach can
be used for any observationlevel covariates included in the detection
portion of the model (i.e., simply place an NA for any covariate during
site/visits when the survey was not performed). Note that the missing
values between the detectionnondetection data y
and the
detection covariates det.covs
will have to align, and if
they don’t spAbundance
will return an error. Further,
spAbundance
does not allow missing values in any sitelevel
covariates, and so if there are any missing values in sitelevel
covariates you must decide the best approach to handle this (e.g.,
impute missing values with the mean, remove sites with any missing
sitelevel covariates).
The data.one.sp
list is already in the required format
for use with the NMix()
function. Here we will model
abundance as a function of a continuous covariate
abund.cov.1
as well as a categorical variable
abund.factor.1
, which we will treat as a random effect. We
can imagine this categorical variable corresponding to a management
unit, aspect of the experimental design, or some other grouping variable
that we want to account for in our model. We model detection probability
as a function of two continuous variables that vary across sites and
replicates. We standardize all continuous covariates by using the
scale()
function in our model specification (note that
standardizing continuous covariates is highly recommended as it helps
aid convergence of the underlying MCMC algorithms):
abund.formula < ~ scale(abund.cov.1) + (1  abund.factor.1)
det.formula < ~ scale(det.cov.1) + scale(det.cov.2)
We always like to get an overview of all the data in our analysis
using the str()
function. This makes it much more clear to
us which data we put into the model and the dimensions of the given data
set.
str(data.one.sp)
List of 4
$ y : int [1:225, 1:3] 1 NA NA NA 0 0 0 0 1 NA ...
$ abund.covs: num [1:225, 1:2] 0.373 0.706 0.202 1.588 0.138 ...
.. attr(*, "dimnames")=List of 2
.. ..$ : NULL
.. ..$ : chr [1:2] "abund.cov.1" "abund.factor.1"
$ det.covs :List of 2
..$ det.cov.1: num [1:225, 1:3] 1.28 NA NA NA 1.04 ...
..$ det.cov.2: num [1:225, 1:3] 2.03 NA NA NA 0.796 ...
$ coords : num [1:225, 1:2] 0 0.0714 0.1429 0.2143 0.2857 ...
.. attr(*, "dimnames")=List of 2
.. ..$ : NULL
.. ..$ : chr [1:2] "X" "Y"
Next, we specify the initial values for the MCMC sampler in
inits
. NMix()
(and all other
spAbundance
model fitting functions) will set initial
values by default, but here we will do this explicitly, since in more
complicated cases setting initial values close to the presumed solutions
can be vital for success of an MCMCbased analysis (for instance, this
is the case when fitting distance sampling models in
spAbundance
). However, for all models described in this
vignette (in particular the nonspatial models), choice of the initial
values is largely inconsequential, with the exception being that
specifying initial values close to the presumed solutions can decrease
the amount of samples you need to run to arrive at convergence of the
MCMC chains. Thus, when first running a model in
spAbundance
, we recommend fitting the model using the
default initial values that spAbundance
provides. The
initial values that spAbundance
chooses will be reported to
the screen when setting verbose = TRUE
. After running the
model for a reasonable period, if you find the chains are taking a long
time to reach convergence, you then may wish to set the initial values
to the mean estimates of the parameters from the initial model fit, as
this will likely help reduce the amount of time you need to run the
model.
The default initial values for abundance and detection regression
coefficients (including the intercepts) are random values from a
standard normal distribution, while the default initial values for the
latent abundance effects are set to the maximum number of individuals
observed at a given site over the replicate surveys at that site. When
fitting an Nmixture model with a negative binomial distribution, the
initial value for the overdispersion parameter is drawn from the prior
distribution. Initial values are specified in a list with the following
tags: N
(latent abundance values), alpha
(detection intercept and regression coefficients), beta
(abundance intercept and regression coefficients), and
kappa
(negative binomial overdispersion parameter). Below
we set all initial values of the regression coefficients to 0, initial
values for the overdispersion parameter to 0.5, and set initial values
for N
based on the count data matrix. For the abundance
(beta
) and detection (alpha
) regression
coefficients, the initial values are passed either as a vector of length
equal to the number of estimated parameters (including an intercept, and
in the order specified in the model formula), or as a single value if
setting the same initial value for all parameters (including the
intercept). Below we take the latter approach. For the negative binomial
overdispersion parameter, the initial value is simply a single numeric
value. To specify the initial values for the latent abundance at each
site (N
), we must ensure we set the value to at least the
maximum number of individuals observed at a site on a given survey,
because we know the true abundance must be greater than or equal to the
number of individuals observed (i.e., assuming no false positives). If
the initial values for N
do not meet this criterion,
NMix()
will fail. spAbundance
will provide a
clear error message if the supplied initial values for N
are invalid. Below we use the raw count data and the
apply()
function to set the initial values to the largest
observed count at each site. For any random effects that are included in
the model, we can also specify the initial values for the random effect
variances (sigma.sq.mu
for abundance and
sigma.sq.p
for detection). By default, these will be drawn
as random values between 0.05 and 2. Here we specify the initial value
for the abundance random effect to 0.5.
# Format with explicit specification of inits for alpha and beta
# with three detection parameters and two abundance parameters
# (including the intercept).
inits < list(alpha = c(0, 0, 0),
beta = c(0, 0),
kappa = 0.5,
sigma.sq.mu = 0.5,
N = apply(data.one.sp$y, 1, max, na.rm = TRUE))
# Format with abbreviated specification of inits for alpha and beta.
inits < list(alpha = 0,
beta = 0,
kappa = 0.5,
sigma.sq.mu = 0.5,
N = apply(data.one.sp$y, 1, max, na.rm = TRUE))
We next specify the priors for the abundance and detection regression
coefficients, as well as the negative binomial overdispersion parameter.
We assume normal priors for both the detection and abundance regression
coefficients. These priors are specified in a list with tags
beta.normal
for abundance and alpha.normal
for
detection parameters (including intercepts). Each list element is then
itself a list, with the first element of the list consisting of the
hypermeans for each coefficient and the second element of the list
consisting of the hypervariances for each coefficient. Alternatively,
the hypermeans and hypervariances can be specified as a single value if
the same prior is used for all regression coefficients. By default,
spAbundance
will set the hypermeans to 0 and the
hypervariances to 100 for the abundance coefficients and 2.72 for the
detection coefficients. The variance of 2.72 corresponds to a relatively
flat prior on the probability scale (0, 1; Lunn
et al. (2013)). For the negative binomial overdispersion
parameter, we will use a uniform prior. This prior is specified as a tag
in the prior list called kappa.unif
, which should be a
vector with two values indicating the lower and upper bound of the
uniform distribution. The default prior is to set the lower bound to 0
and the upper bound to 100. Recall that lower values of
kappa
indicate substantial overdispersion and high values
of kappa
indicate minimal overdispersion. If there is
little support for overdispersion when fitting a negative binomial
model, we will likely see the estimates of kappa
be close
to the upper bound of the uniform prior distribution. For the default
prior distribution, if the estimates of kappa
are very
close to 100, this indicates little support for overdispersion in the
model, and we can likely switch to using a Poisson distribution (which
would also likely be favored by model comparison approaches). For models
with random effects in the abundance and/or detection portion of the
Nmixture model, we can also specify the prior for the random effect
variance parameter (sigma.sq.mu
for abundance and
sigma.sq.p
for detection). We assume inverseGamma priors
for these variance parameters and specify them with the tags
sigma.sq.mu.ig
and sigma.sq.p.ig
,
respectively. These priors are set as a list with two components, where
the first element is the shape parameter and the second element is the
scale parameter. The shape and scale parameters can be specified as a
single value or as vectors with length equal to the number of random
effects included in the model. The default prior distribution for random
effect variances is 0.1 for both the shape and scale parameters. Below
we use default priors for all parameters, but specify them explicitly
for clarity.
priors < list(alpha.normal = list(mean = 0, var = 2.72),
beta.normal = list(mean = 0, var = 100),
kappa.unif = c(0, 100),
sigma.sq.mu.ig = list(0.1, 0.1))
The next four arguments (tuning
, n.batch
,
batch.length
, and accept.rate
) are all related
to the specific type of MCMC sampler we use when we fit Nmixture models
in spAbundance
. The parameters in Nmixture models are all
estimated using a MetropolisHastings step, which can often be slow and
inefficient, leading to slow mixing and convergence of the MCMC chains.
To try and mitigate the slow mixing and convergence issues, we update
all parameters in Nmixture models using an algorithm called an adaptive
MetropolisHastings algorithm (see Roberts and
Rosenthal (2009) for more details on this algorithm). In this
approach, we break up the total number of MCMC samples into a set of
“batches”, where each batch has a specific number of MCMC samples. Thus,
we must specify the total number of batches (n.batch
) as
well as the number of MCMC samples each batch contains
(batch.length
) when specifying the function arguments. 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 of all model parameters is
reached. We generally recommend setting batch.length = 25
,
but in certain situations this can be increased to a larger number of
samples (e.g., 100), which can result in moderate decreases in run time.
Here we set n.batch = 1600
for a total of 40,000 MCMC
samples for each MCMC chain we run.
batch.length < 25
n.batch < 1600
# Total number of MCMC samples per chain
batch.length * n.batch
[1] 40000
Importantly, we also need to specify a target acceptance rate and
initial tuning parameters for the abundance and detection regression
coefficients (and the negative binomial overdispersion parameter and any
latent random effects if applicable). These are both features of the
adaptive algorithm we use to sample these parameters. In this adaptive
MetropolisHastings algorithm, we propose new values for the parameters
from some proposal distribution, compare them to our previous values,
and use a statistical algorithm to determine if we should accept the new
proposed value or keep the old one. The accept.rate
argument specifies the ideal proportion of times we will accept the
newly proposed values for these parameters. Roberts and Rosenthal (2009) show that if we
accept new values around 43% of the time, then this will lead to optimal
mixing and convergence of the MCMC chains. Following these
recommendations, we should strive for an algorithm that accepts new
values about 43% of the time. Thus, we recommend setting
accept.rate = 0.43
unless you have a specific reason not to
(this is the default value). The values specified in the
tuning
argument help control the initial values we will
propose for the abundance/detection coefficients and the negative
binomial overdispersion parameter. These values are supplied as input in
the form of a list with tags beta
, alpha
, and
kappa
. The initial tuning value can be any value greater
than 0, but we generally recommend starting the value out around 0.5.
These tuning values can also be thought of as tuning “variances”, as it
is these values that control the variance of the distribution we use to
generate newly proposed values for the parameters we are trying to
estimate with our MCMC algorithm. In short, the new values that we
propose for the parameters beta
, alpha
, and
kappa
come from a normal distribution with mean equal to
the current value for the given parameter and the variance equal to the
tuning parameter. Thus, the smaller this tuning parameter/variance is,
the closer our proposed values will be to the current value, and vise
versa for large values of the tuning parameter. The “ideal” value of the
tuning variance will depend on the data set, the parameter, and how much
uncertainty there is in the estimate of the parameter. This initial
tuning value that we supply is the first tuning variance that will be
used for the given parameter, and our adaptive algorithm will adjust
this tuning parameter after each batch to yield acceptance rates of
newly proposed values that are close to our target acceptance rate that
we specified in the accept.rate
argument. Information on
the acceptance rates for a few of the parameters in your model will be
displayed when setting verbose = TRUE
. After some initial
runs of the model, if you notice the final acceptance rate 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. While use of this algorithm requires
us to specify more arguments than if we didn’t “adaptively tune” our
proposal variances, this leads to much shorter run times compared to a
more simplistic approach where we do not have an “adaptive” sampling
approach, and it should thus save us time in the long haul when waiting
for these models to run. For our example here, we set the initial tuning
values to 0.5 for beta
, alpha
, and
kappa
. For models with random effects in either the
abundance or detection portions of the model, we also need to specify
tuning parameters for the latent random effect values
(beta.star
for abundance and alpha.star
for
detection). We similarly set these to 0.5.
tuning < list(beta = 0.5, alpha = 0.5, kappa = 0.5, beta.star = 0.5)
# accept.rate = 0.43 by default, so we do not specify it.
We also need to specify the length of burnin (n.burn
),
the rate at which we want to thin the posterior samples
(n.thin
), and the number of MCMC chains to run
(n.chains
). Note that currently spAbundance
runs multiple chains sequentially and does not allow chains to be run
simultaneously in parallel across multiple threads, which is something
we hope to implement in future package development. Instead, we allow
for withinchain parallelization using the n.omp.threads
argument. We can set n.omp.threads
to a number greater than
1 and smaller than the number of threads on the computer you are using.
Generally, setting n.omp.threads > 1
will not result in
decreased run times for nonspatial models in spAbundance
,
but can substantially decrease run time when fitting spatial models
(Finley, Datta, and Banerjee 2020). Here
we set n.omp.threads = 1
.
We will run the model using three chains to assess convergence using the GelmanRubin diagnostic (Rhat; Brooks and Gelman (1998)).
n.burn < 20000
n.thin < 20
n.chains < 3
We are now almost set to run the model. The family
argument is used to indicate whether we want to model abundance with a
Poisson distribution (Poisson
) or a negative binomial
distribution (NB
). Here we will start with a Poisson
distribution (the default), which we will compare to a model with a
negative binomial distribution later. The verbose
argument
is a logical value indicating whether or not MCMC sampler progress is
reported to the screen. If verbose = TRUE
, sampler progress
is reported to the screen. The argument n.report
specifies
the interval to report the MetropolisHastings sampler acceptance rate.
Note that n.report
is specified in terms of batches, not
the overall number of samples. Below we set n.report = 400
,
which will result in information on the acceptance rate and tuning
parameters every 400th batch (not sample).
We now are set to fit the model.
out < NMix(abund.formula = abund.formula,
det.formula = det.formula,
data = data.one.sp,
inits = inits,
priors = priors,
n.batch = n.batch,
batch.length = batch.length,
tuning = tuning,
n.omp.threads = 1,
n.report = 400,
family = 'Poisson',
verbose = TRUE,
n.burn = n.burn,
n.thin = n.thin,
n.chains = n.chains)

Preparing to run the model


Model description

Poisson Nmixture model with 225 sites.
Samples per Chain: 40000 (1600 batches of length 25)
Burnin: 20000
Thinning Rate: 20
Number of Chains: 3
Total Posterior Samples: 3000
Source compiled with OpenMP support and model fit using 1 thread(s).
Adaptive Metropolis with target acceptance rate: 43.0

Chain 1

Sampling ...
Batch: 400 of 1600, 25.00%
Parameter Acceptance Tuning
beta[1] 40.0 0.18211
beta[2] 24.0 0.16811
alpha[1] 56.0 0.26630
alpha[2] 48.0 0.27716
alpha[3] 52.0 0.30025

Batch: 800 of 1600, 50.00%
Parameter Acceptance Tuning
beta[1] 36.0 0.17150
beta[2] 48.0 0.17150
alpha[1] 40.0 0.28276
alpha[2] 32.0 0.30025
alpha[3] 28.0 0.33183

Batch: 1200 of 1600, 75.00%
Parameter Acceptance Tuning
beta[1] 44.0 0.17850
beta[2] 28.0 0.17497
alpha[1] 40.0 0.28847
alpha[2] 40.0 0.27716
alpha[3] 20.0 0.35946

Batch: 1600 of 1600, 100.00%

Chain 2

Sampling ...
Batch: 400 of 1600, 25.00%
Parameter Acceptance Tuning
beta[1] 68.0 0.16152
beta[2] 36.0 0.17497
alpha[1] 60.0 0.28847
alpha[2] 32.0 0.28847
alpha[3] 40.0 0.33183

Batch: 800 of 1600, 50.00%
Parameter Acceptance Tuning
beta[1] 36.0 0.17150
beta[2] 32.0 0.16811
alpha[1] 40.0 0.27716
alpha[2] 40.0 0.28276
alpha[3] 32.0 0.34537

Batch: 1200 of 1600, 75.00%
Parameter Acceptance Tuning
beta[1] 48.0 0.18579
beta[2] 52.0 0.17497
alpha[1] 44.0 0.30631
alpha[2] 36.0 0.29430
alpha[3] 32.0 0.35946

Batch: 1600 of 1600, 100.00%

Chain 3

Sampling ...
Batch: 400 of 1600, 25.00%
Parameter Acceptance Tuning
beta[1] 48.0 0.17850
beta[2] 40.0 0.15832
alpha[1] 24.0 0.28276
alpha[2] 32.0 0.31250
alpha[3] 44.0 0.36672

Batch: 800 of 1600, 50.00%
Parameter Acceptance Tuning
beta[1] 40.0 0.17497
beta[2] 32.0 0.17850
alpha[1] 40.0 0.27168
alpha[2] 44.0 0.30631
alpha[3] 40.0 0.35234

Batch: 1200 of 1600, 75.00%
Parameter Acceptance Tuning
beta[1] 28.0 0.16811
beta[2] 40.0 0.17497
alpha[1] 24.0 0.29430
alpha[2] 28.0 0.27716
alpha[3] 24.0 0.33853

Batch: 1600 of 1600, 100.00%
NMix()
returns a list of class NMix
with a
suite of different objects, many of them being coda::mcmc
objects of posterior samples. The “Preparing to run the model” section
will print information on default priors or initial values that are used
when they are not specified in the function call. Here we specified
everything explicitly so no information was reported.
We next use the summary()
function on the resulting
NMix()
object for a concise, informative summary of the
regression parameters and convergence of the MCMC chains.
summary(out)
Call:
NMix(abund.formula = abund.formula, det.formula = det.formula,
data = data.one.sp, inits = inits, priors = priors, tuning = tuning,
n.batch = n.batch, batch.length = batch.length, family = "Poisson",
n.omp.threads = 1, verbose = TRUE, n.report = 400, n.burn = n.burn,
n.thin = n.thin, n.chains = n.chains)
Samples per Chain: 40000
Burnin: 20000
Thinning Rate: 20
Number of Chains: 3
Total Posterior Samples: 3000
Run Time (min): 0.4257
Abundance (log scale):
Mean SD 2.5% 50% 97.5% Rhat ESS
(Intercept) 0.1468 0.2196 0.6161 0.1403 0.2626 1.0520 546
scale(abund.cov.1) 0.1157 0.0788 0.0400 0.1150 0.2716 1.0146 2753
Abundance Random Effect Variances (log scale):
Mean SD 2.5% 50% 97.5% Rhat ESS
(Intercept)abund.factor.1 0.3672 0.3502 0.0778 0.2812 1.1465 1.0638 2229
Detection (logit scale):
Mean SD 2.5% 50% 97.5% Rhat ESS
(Intercept) 0.4670 0.1839 0.1039 0.4689 0.8132 1.0001 1380
scale(det.cov.1) 0.5399 0.1359 0.8080 0.5360 0.2751 1.0027 2833
scale(det.cov.2) 0.9961 0.1713 0.6667 0.9949 1.3414 1.0027 2497
Notice that the abundance coefficients are printed on the log scale
and the detection coefficients are on the logit scale. We see that our
species is quite rare, with an average expected abundance of
approximately 0.86, exp(0.15). We also see a moderate, positive effect
of the covariate, although the 95% credible interval overlaps 0. The
random effect variance mean is approximately 0.37, indicating this is a
potentially important source of variability in abundance for the
species. Note that we can also change the quantiles that are returned in
the summary output if we desire to see additional quantiles. This is
controlled with the quantiles
argument.
Call:
NMix(abund.formula = abund.formula, det.formula = det.formula,
data = data.one.sp, inits = inits, priors = priors, tuning = tuning,
n.batch = n.batch, batch.length = batch.length, family = "Poisson",
n.omp.threads = 1, verbose = TRUE, n.report = 400, n.burn = n.burn,
n.thin = n.thin, n.chains = n.chains)
Samples per Chain: 40000
Burnin: 20000
Thinning Rate: 20
Number of Chains: 3
Total Posterior Samples: 3000
Run Time (min): 0.4257
Abundance (log scale):
Mean SD 25% 50% 75% Rhat ESS
(Intercept) 0.1468 0.2196 0.2762 0.1403 0.0125 1.0520 546
scale(abund.cov.1) 0.1157 0.0788 0.0643 0.1150 0.1701 1.0146 2753
Abundance Random Effect Variances (log scale):
Mean SD 25% 50% 75% Rhat ESS
(Intercept)abund.factor.1 0.3672 0.3502 0.1849 0.2812 0.452 1.0638 2229
Detection (logit scale):
Mean SD 25% 50% 75% Rhat ESS
(Intercept) 0.4670 0.1839 0.3433 0.4689 0.5885 1.0001 1380
scale(det.cov.1) 0.5399 0.1359 0.6335 0.5360 0.4498 1.0027 2833
scale(det.cov.2) 0.9961 0.1713 0.8806 0.9949 1.1055 1.0027 2497
The model summary also provides information on convergence of the MCMC chains in the form of the GelmanRubin diagnostic (Brooks and Gelman 1998) and the effective sample size (ESS) of the posterior samples. Here we find all Rhat values are less than 1.1. The ESS values are adequately high for all model parameters, indicating adequate mixing of the MCMC chains.
We can use the plot()
function to generate a simple
trace plot of the MCMC chains to provide additional confidence in the
convergence (or nonconvergence) of the model. The plotting
functionality for each model type in spAbundance
takes
three arguments: x
(the resulting object from fitting the
model), param
(the parameter name that you want to
display), and density
(a logical value indicating whether
to also generate a density plot in addition to the traceplot). To see
the parameter names available to use with plot()
for a
given model type, you can look at the manual page for the function,
which for models generated from NMix()
can be accessed with
?plot.NMix
.
# Abundance regression coefficients
plot(out, param = 'beta', density = FALSE)
# Detection regression coefficients
plot(out, param = 'alpha', density = FALSE)
Posterior predictive checks
The function ppcAbund()
performs a posterior predictive
check on all spAbundance
model objects as a GoodnessofFit
(GoF) assessment. The fundamental idea of GoF testing is that a good
model should generate data that closely align with the observed data. If
there are drastic differences in the true data from the data generated
under the model, our model is likely not very useful (Hobbs and Hooten 2015). In
spAbundance
, we perform posterior predictive checks using
the following approach:
 Fit the model using any of the modelfitting functions (here
NMix()
), which will generate replicated values for all observed data points from the posterior predictive distribution of the data.  Optionally bin both the actual and the replicated count data in some manner, such as by site or replicate.
 Compute a fit statistic on both the actual data and also on the modelgenerated ‘replicate data’.
 Compare the fit statistics for the true data and replicate data. If they are widely different, this suggests a lack of fit of the model to the actual data set at hand.
To perform a posterior predictive check, we send the resulting
NMix
model object as input to the ppcAbund()
function, along with a fit statistic (fit.stat
) and a
numeric value indicating how to group, or bin, the data
(group
). Currently supported fit statistics include the
FreemanTukey statistic and the ChiSquared statistic
(freemantukey
or chisquared
, respectively,
Kéry and Royle (2016)). Currently,
ppcAbund()
allows the user to group the data by row (site;
group = 1
), column (replicate; group = 2
), or
not at all (group = 0
). ppcAbund()
will then
return a set of posterior samples for the fit statistic (or discrepancy
measure) using the actual data (fit.y
) and model generated
replicate data set (fit.y.rep
), summed across all data
points in the chosen manner. For example, when setting
group = 1
, spAbundance
will first sum all of
the count values at a given site across all replicates at that site,
then calculate the fit statistic using the sitelevel sums. When setting
group = 0
, spAbundance
calculates the fit
statistic directly on the count value associated with each site and
replicate. We generally recommend performing a posterior predictive
check using multiple forms of grouping, as they may reveal (or fail to
reveal) different inadequacies of the model for the specific data set at
hand (Kéry and Royle 2016). Throughout
this vignette, we will display different types of posterior predictive
checks using different combinations of the fit statistic and grouping
approach.
The resulting values from a call to ppcAbund()
can be
used with the summary()
function to generate a Bayesian
pvalue, which is the probability, under the fitted model, to obtain a
value of the fit statistic that is more extreme (i.e., larger) than the
one observed, i.e., for the actual data. Bayesian pvalues are sensitive
to individual values, so we may also want to explore the discrepancy
measures for each (potentially “grouped”) data point.
ppcAbund()
returns a matrix of posterior quantiles for the
fit statistic for both the observed (fit.y.group.quants
)
and model generated, replicate data
(fit.y.rep.group.quants
) for each “grouped” data point.
We next perform a posterior predictive check using the FreemanTukey
statistic grouping the data by sites. We summarize the posterior
predictive check with the summary()
function, which reports
a Bayesian pvalue. A Bayesian pvalue 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 (Hobbs and Hooten 2015). As always with a
simulationbased analysis using MCMC, you will get numerically slightly
different values.
Call:
ppcAbund(object = out, fit.stat = "freemantukey", group = 1)
Samples per Chain: 40000
Burnin: 20000
Thinning Rate: 20
Number of Chains: 3
Total Posterior Samples: 3000
Bayesian pvalue: 0
Fit statistic: freemantukey
The Bayesian pvalue here is 0, indicating inadequate model fit. In
particular, a value close to 0 indicates there is more variability in
the observed data points than the replicate data points generated by our
model. This could be a result of (1) missing sources of variability in
true abundance; (2) missing sources of variability in detection
probability; or (3) missing sources of variability in both abundance and
detection. We will see later in the vignette that our low Bayesian
pvalue in this case is a result of additional spatial variability in
abundance. See the introductory
spOccupancy
vignette for ways to further explore
resulting objects from posterior predictive checks.
A brief note on generation of replicate data for posterior predictive checks
For Nmixture models, the ppcAbund()
function contains
an additional argument type
, which is used to indicate the
specific form of replicate data value to generate from the given model.
This argument can take two values: marginal
and
conditional
. For type = 'conditional'
, the
replicate data values are generated conditional on the estimated values
of the latent abundance random effects N
, while for
type = 'marginal'
, the replicate data values are not
generated conditional on the latent effects N
. When
calculating the replicate data, we generate a different replicate data
set for each MCMC sample \(l\), such
that we can generate a full distribution of replicate data values (i.e.,
the posterior predictive distribution of the data). More specifically,
for the conditional replicate data values, the value at site \(j\) and replicate \(k\) for MCMC sample \(l\), denoted as \(y^{(l)}_{\text{rep}, j, k}\) is calculated
according to
\[\begin{equation} y^{(l)}_{\text{rep}, j, k} \sim \text{Binomial}(N^{(l)}_j, p^{(l)}_{j, k}), \end{equation}\]
where \(N^{(l)}_j\) is the estimated abundance at site \(j\) for MCMC sample \(l\) and \(p^{(l)}_{j, k}\) is the probability of detecting an individual at site \(j\) during replicate \(k\) for MCMC sample \(l\). The reason we refer to these as “conditional” replicate data values is because the values of \(N^{(l)}_j\), which are estimated directly when fitting the model, are actually conditional on the observed data values. Recall that when supplying initial values for \(N_j\), we had to ensure the initial values were greater than or equal to the observed data values. This showcases the conditional nature of the \(N^{(l)}_j\) estimates (which form the posterior distribution for \(N_j\)) when fitting the model, as the smallest value \(N^{(l)}_j\) can take at any given iteration of the MCMC algorithm \(l\) is the largest observed number of individuals at that site across all \(K_j\) replicates. Thus, replicate values generated in this manner are in a sense conditional on the observed data values given the constraints on \(N^{(l)}_j\).
Instead, we can imagine calculating replicate data in a slightly different way that eliminates the conditional nature of the \(N^{(l)}_j\) estimates. The second approach, which we refer to as the “marginal” approach, generates replicate data by first predicting a value of latent abundance at site \(j\) using the expected abundance at site \(j\) at MCMC sample \(l\) (i.e., \(\mu^{(l)}_j\)), and then subsequently generating a replicate data point for each replicate \(k\) at site \(j\). More specifically, we generate marginal replicate data values according to
\[\begin{equation} \begin{split} N^{(l)}_{\text{rep}, j} &\sim \text{Poisson}(\mu^{(l)}_j), \\ y^{(l)}_{\text{rep}, j, k} &\sim \text{Binomial}(N^{(l)}_{\text{rep}, j}, p^{(l)}_{j, k}). \end{split} \end{equation}\]
Note that the Poisson distribution would be replaced with a negative binomial distribution if that is used to fit the model. By calculating replicate values in this manner, the \(N^{(l)}_{\text{rep}, j}\) values are no longer required to be at least as large as the maximum number of individuals ever observed at a site.
We performed a small simulation study to assess the ability of the
two approaches to detect unmodeled spatial heterogeneity in abundance.
Very briefly, we simulated 100 data sets where count data were generated
at 225 sites with a maximum of three replicate surveys at any given
site. We simulated abundance as a function of two covariates and
additional residual spatial autocorrelation. We then fit an Nmixture
model to each of the data sets using NMix()
in which we did
not account for the residual spatial autocorrelation. Finally, we
performed a posterior predictive check and calculated a Bayesian pvalue
using both the marginal and conditional approach. For the posterior
predictive check, we used the FreemanTukey statistic and grouped the
data by site. Ideally, our posterior predictive check should indicate to
us that there is additional variation in the data that we are not
accounting for. The full script used to generate the data is available
for download here. We
load the conditional and marginal Bayesian pvalues from the simulation
study below and look at their values.
# Load the results from online
load(url("https://www.jeffdoser.com/files/misc/nmixsim/nmixsimresults.rda"))
# Conditional Bayesian pvalues
round(conditional.bps, 2)
[1] 0.00 0.02 0.05 0.01 0.03 0.00 0.31 0.63 0.07 0.00 0.00 0.01 0.04 0.00 0.01
[16] 0.00 0.00 0.00 0.01 0.01 0.03 0.00 0.00 0.02 0.00 0.21 0.00 0.00 0.09 0.00
[31] 0.07 0.00 0.44 0.04 0.17 0.00 0.00 0.00 0.00 0.00 0.05 0.01 0.01 0.00 0.00
[46] 0.01 0.00 0.00 0.36 0.03 0.00 0.00 0.01 0.01 0.08 0.00 0.00 0.00 0.02 0.01
[61] 0.00 0.04 0.00 0.13 0.00 0.00 0.01 0.06 0.06 0.00 0.00 0.00 0.00 0.04 0.00
[76] 0.00 0.24 0.00 0.00 0.02 0.00 0.00 0.00 0.00 0.02 0.01 0.00 0.00 0.00 0.52
[91] 0.04 0.38 0.01 0.00 0.01 0.01 0.08 0.00 0.00 0.02
# Marginal Bayesian pvalues
round(marginal.bps, 2)
[1] 0.00 0.00 0.00 0.00 0.00 0.00 0.45 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
[16] 0.00 0.00 0.00 0.00 0.00 0.01 0.00 0.00 0.00 0.00 0.05 0.00 0.00 0.00 0.00
[31] 0.00 0.00 0.23 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
[46] 0.00 0.00 0.00 0.03 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
[61] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.01 0.00 0.00 0.00 0.00 0.00 0.00 0.00
[76] 0.00 0.02 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.33
[91] 0.00 0.01 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
# How many conditional values are not rejected (e.g., > 0.1 and < 0.9)
sum(conditional.bps > 0.1 & conditional.bps < 0.9)
[1] 10
# How many marginal values are not rejected (e.g., > 0.1 and < 0.9)
sum(marginal.bps > 0.1 & marginal.bps < 0.9)
[1] 3
Here we see the marginal approach performs slightly better at detecting the unmodeled spatial heterogeneity compared to the conditional Bayesian pvalues. In particular, we find support for inadequate model fit in 97 out of 100 simulations for the marginal approach, while we find support for inadequate model fit in 90 of the simulations using the conditional approach. Because we simulated the data with additional spatial variation in abundance that we are not accounting for, our two approaches are doing a good job of detecting the inadequate performance of the nonspatial Nmixture model, with the marginal approach performing slightly better.
This concept of what we called “marginal” and “conditional” replicate values is not new to the statistical ecology literature. In particular, Conn et al. (2018) discuss this concept (and many more aspects of GoodnessofFit checking in hierarchical Bayesian models), and give a nice visual summary of this idea in their Figure 4 in the context of hierarchical spatial models. Marshall and Spiegelhalter (2003) discuss a similar concept in the context of disease mapping, in which they refer to what we call a “marginal” approach as a “mixed” posterior predictive check and advocate for its use. We believe additional exploration of posterior predictive checks and other associated GoodnessofFit tests in hierarchical models like Nmixture models, hierarchical distance sampling models, and occupancy models is an important avenue of future research. Here we provide options to perform posterior predictive checks using both types (with the default being marginal), and encourage users to explore multiple forms of performing posterior predictive checks to increase one’s comfort that the fitted model is an adequate representation of the datagenerating process. Our small and simple simulation study showed similar performance of the two approaches with the marginal approach performing slightly better than the conditional approach, but additional, more indepth simulation studies are needed to better assess their performance. These results align with Marshall and Spiegelhalter (2003), who found the “marginal” approach to outperform the more traditional “conditional” approach in the context of a spatial regression model. Further exploring this potential pattern for Nmixture, hierarchical distance sampling, and occupancytype models is an important area of future work. This discussion applies to all Nmixture models discussed in this vignette.
Model selection using WAIC
Posterior predictive checks allow us to assess how well our model fits the data, but they are not very useful if we want to compare multiple competing models and ultimately select a final model based on some criterion. Bayesian model selection is very much a constantly changing field. See Hooten and Hobbs (2015) for an accessible overview of Bayesian model selection for ecologists.
For Bayesian hierarchical models like Nmixture models, the most common Bayesian model selection criterion, the deviance information criterion or DIC, is not applicable (Hooten and Hobbs 2015). Instead, the Widely Applicable Information Criterion (Watanabe 2010) is often recommended to compare a set of models and select the bestperforming model for final analysis.
The WAIC is calculated for all spAbundance
model objects
using the function waicAbund()
. We calculate the WAIC
as
\[ \text{WAIC} = 2 \times (\text{elppd}  \text{pD}), \]
where elppd is the expected log pointwise predictive density and pD is the effective number of parameters. We calculate elppd by calculating the likelihood for each posterior sample, taking the mean of these likelihood values, taking the log of the mean of the likelihood values, and summing these values across all sites. We calculate the effective number of parameters by calculating the variance of the log likelihood for each site taken over all posterior samples, and then summing these values across all sites.
We calculate the WAIC using waicAbund()
for our model
below (as always, note some slight differences with your solutions due
to Monte Carlo error).
waicAbund(out)
N.max not specified. Setting upper index of integration of N to 10 plus
the largest estimated abundance value at each site in object$N.samples
elpd pD WAIC
369.16550 22.74705 783.82511
Note the somewhat cryptic message that is displayed to the screen
when you run the previous line. When calculating WAIC for Nmixture
models (and hierarchical
distance sampling models), we need to integrate out the latent
abundance values, which requires setting an upper bound to the potential
value of the latent abundance values \(N_j\) at each spatial location. By default,
waicAbund()
will set that upper bound to the largest
abundance value at each site plus 10 (as indicated by the message). This
upper bound can be controlled further with the N.max
argument in waicAbund()
. See the help page for
waicAbund
for details.
Now let’s do a bit of model comparison. We will fit the same model as
before, except now we will use a negative binomial distribution. We fit
the model below, setting verbose = FALSE
to hide the
messages printed by the model fitting functions.
# Approx run time: .5 minute
out.2 < NMix(abund.formula = ~ abund.cov.1 + (1  abund.factor.1),
det.formula = det.formula, data = data.one.sp,
inits = inits, priors = priors, n.batch = n.batch,
batch.length = batch.length, tuning = tuning,
n.omp.threads = 1, family = 'NB', verbose = FALSE,
n.burn = n.burn, n.thin = n.thin, n.chains = n.chains)
# Compare models with WAIC
# Poisson
waicAbund(out)
N.max not specified. Setting upper index of integration of N to 10 plus
the largest estimated abundance value at each site in object$N.samples
elpd pD WAIC
369.16550 22.74705 783.82511
# Negative binomial
waicAbund(out.2)
N.max not specified. Setting upper index of integration of N to 10 plus
the largest estimated abundance value at each site in object$N.samples
elpd pD WAIC
355.28485 15.51822 741.60615
From this simple model comparison exercise, we see the negative binomial model outperforms the Poisson model, indicating there is substantial overdispersion in latent abundance that is not captured by the covariate/random effect included in the model.
Prediction
All model objects from a call to spAbundance
modelfitting functions can be used with predict()
to
generate a series of posterior predictive samples at new locations,
given the values of all covariates used in the model fitting process.
Given that we standardized the covariate values when we fit the model,
we need to standardize the covariate values for prediction using the
exact same values of the mean and standard deviation of the covariate
values used to fit the data.
We envision the predict()
functions will be used most
often for two purposes: (1) generating marginal effects plots; and (2)
predicting across some region of interest. The approach to do the
prediction is the same for both purposes. In this vignette, we will use
the predict()
function to generate marginal effects plots
for the effect of the covariate on abundance. See the hierarchical
distance sampling vignette for examples of how to do prediction with
the purpose of creating a map of abudance across space (the syntax is
exactly the same for Nmixture models).
First, we create a set of evenlyspaced values of the covariate across the range of observed values, and subsequently standardize them by the mean and standard deviation of the covariate used to fit the model. Note that we could of course predict across some other range of values that are interesting for a given study, of course subject to the caveat that any extrapolation beyond the observed range of a covariate may be risky.
cov.pred.vals < seq(min(dataNMixSim$abund.covs[, 'abund.cov.1']),
max(dataNMixSim$abund.covs[, 'abund.cov.1']),
length.out = 100)
# Scale predicted values by mean and sd used to fit the model
cov.pred.vals.s < (cov.pred.vals  mean(dataNMixSim$abund.covs[, 'abund.cov.1'])) /
sd(dataNMixSim$abund.covs[, 'abund.cov.1'])
For NMix()
, the predict()
function takes
four arguments:

object
: theNMix
fitted model object. 
X.0
: a matrix or data frame consisting of the design matrix for the prediction locations (which must include an intercept if our model contained one). 
ignore.RE
: a logical value indicating whether or not to remove random effects from the predicted values. By default, this is set toFALSE
, and so prediction will include the random effects. 
type
: a quoted keyword indicating whether we want to predictabundance
ordetection
. This is by default set toabundance
.
Here, we only seek to predict to visualize the relationship between
the covariate and abundance, and so we will predict without the random
effect and set ignore.RE = TRUE
. Thus, our design matrix
that we form below only consists of the intercept and covariate value.
If we were to predict with the random effect, the random effect values
at the prediction locations would be included as another column in
X.0
.
X.0 < as.matrix(data.frame(intercept = 1, abund.cov.1 = cov.pred.vals.s))
out.pred < predict(out.2, X.0, ignore.RE = TRUE)
str(out.pred)
List of 3
$ mu.0.samples: 'mcmc' num [1:3000, 1:100] 0.735 0.596 0.859 0.678 0.556 ...
.. attr(*, "mcpar")= num [1:3] 1 3000 1
$ N.0.samples : 'mcmc' num [1:3000, 1:100] 0 0 3 0 2 0 1 1 2 0 ...
.. attr(*, "mcpar")= num [1:3] 1 3000 1
$ call : language predict.NMix(object = out.2, X.0 = X.0, ignore.RE = TRUE)
 attr(*, "class")= chr "predict.NMix"
The resulting object consists of posterior predictive samples for the
expected abundances (mu.0.samples
) and latent abundance
values (N.0.samples
). The beauty of the Bayesian paradigm,
and the MCMC computing machinery, is that these predictions all have
fully propagated uncertainty. Below, we produce a plot showing the
relationship between expected abundance and the covariate, with the 95%
credible interval shown in grey and the posterior median shown with the
black line.
# Get the lower bound, median, and 95% credible interval
mu.0.quants < apply(out.pred$mu.0.samples, 2, quantile,
prob = c(0.025, 0.5, 0.975))
mu.plot.dat < data.frame(mu.med = mu.0.quants[2, ],
mu.low = mu.0.quants[1, ],
mu.high = mu.0.quants[3, ],
abund.cov.1 = cov.pred.vals)
ggplot(mu.plot.dat, aes(x = abund.cov.1, y = mu.med)) +
geom_ribbon(aes(ymin = mu.low, ymax = mu.high), fill = 'grey70', alpha = 0.5) +
geom_line() +
theme_bw() +
labs(x = 'Covariate', y = 'Expected abundance')
Here we see a positive relationship between expected abundance and the covariate, but there is substantial uncertainty in this relationship.
Singlespecies spatial Nmixture models
Basic model description
When working across large spatial domains, accounting for residual spatial autocorrelation in species distributions can often improve predictive performance of a model, leading to more accurate predictions of species abundance patterns (Guélat and Kéry 2018). We here extend the basic singlespecies Nmixture model to incorporate a spatial random effect that accounts for unexplained spatial variation in species abundance across a region of interest. Let \(\boldsymbol{s}_j\) denote the geographical coordinates of site \(j\) for \(j = 1, \dots, J\). In all spatiallyexplicit models, we include \(\boldsymbol{s}_j\) directly in the notation of spatiallyindexed variables to indicate the model is spatiallyexplicit. More specifically, the expected abundance at site \(j\) with coordinates \(\boldsymbol{s}_j\), \(\mu(\boldsymbol{s}_j)\), now takes the form
\[\begin{equation} \text{log}(\mu(\boldsymbol{s}_j) = \boldsymbol{x}(\boldsymbol{s}_j)^{\top}\boldsymbol{\beta} + \text{w}(\boldsymbol{s}_j), \end{equation}\]
where \(\text{w}(\boldsymbol{s}_j)\) is a spatial random effect modeled with a Nearest Neighbor Gaussian Process (NNGP; Datta et al. (2016)). More specifically, we have
\[\begin{equation} \textbf{w}(\boldsymbol{s}) \sim N(\boldsymbol{0}, \boldsymbol{\tilde{\Sigma}}(\boldsymbol{s}, \boldsymbol{s}', \boldsymbol{\theta})), \end{equation}\]
where \(\boldsymbol{\tilde{\Sigma}}(\boldsymbol{s}, \boldsymbol{s}', \boldsymbol{\theta})\) is the NNGPderived spatial covariance matrix that originates from the full \(J \times J\) covariance matrix \(\boldsymbol{\Sigma}(\boldsymbol{s}, \boldsymbol{s}', \boldsymbol{\theta})\) 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})\) that govern the spatial process. The vector \(\boldsymbol{\theta}\) is equal to \(\boldsymbol{\theta} = \{\sigma^2, \phi, \nu\}\), where \(\sigma^2\) is a spatial variance parameter, \(\phi\) is a spatial decay parameter, and \(\nu\) is a spatial smoothness parameter. \(\nu\) is only specified when using a Matern correlation function. The detection portion of the Nmixture model remains unchanged from the nonspatial Nmixture model. The NNGP is a computationally efficient alternative to working with a full Gaussian process model, which is notoriously slow for even moderately large data sets. See Datta et al. (2016) and Finley et al. (2019) for complete statistical details on the NNGP.
Fitting singlespecies spatial Nmixture models with
spNMix()
We will fit the same Nmixture model that we fit previously using
NMix()
, but we will now make the model spatiallyexplicit
by incorporating a spatial process with spNMix()
. The
spNMix()
function fits singlespecies spatial Nmixture
models.
spNMix(abund.formula, det.formula, data, inits, priors, tuning,
cov.model = 'exponential', NNGP = TRUE,
n.neighbors = 15, search.type = 'cb',
n.batch, batch.length, accept.rate = 0.43, family = 'Poisson',
n.omp.threads = 1, verbose = TRUE, n.report = 100,
n.burn = round(.10 * n.batch * batch.length), n.thin = 1,
n.chains = 1, ...)
The arguments to spNMix()
are very similar to those we
saw with NMix()
, with a few additional components. The
abundance (abund.formula
) and detection
(det.formula
) formulas, as well as the list of data
(data
), take the same form as we saw in
NMix()
, with random slopes and intercepts allowed in both
the abundance and detection models. Notice the coords
matrix in the data.one.sp
list of data. We did not use this
for NMix()
, but specifying the spatial coordinates in
data
is necessary for all spatially explicit models in
spAbundance
.
abund.formula < ~ scale(abund.cov.1) + (1  abund.factor.1)
det.formula < ~ scale(det.cov.1) + scale(det.cov.2)
str(data.one.sp) # coords is required for spNMix()
List of 4
$ y : int [1:225, 1:3] 1 NA NA NA 0 0 0 0 1 NA ...
$ abund.covs: num [1:225, 1:2] 0.373 0.706 0.202 1.588 0.138 ...
.. attr(*, "dimnames")=List of 2
.. ..$ : NULL
.. ..$ : chr [1:2] "abund.cov.1" "abund.factor.1"
$ det.covs :List of 2
..$ det.cov.1: num [1:225, 1:3] 1.28 NA NA NA 1.04 ...
..$ det.cov.2: num [1:225, 1:3] 2.03 NA NA NA 0.796 ...
$ coords : num [1:225, 1:2] 0 0.0714 0.1429 0.2143 0.2857 ...
.. attr(*, "dimnames")=List of 2
.. ..$ : NULL
.. ..$ : chr [1:2] "X" "Y"
The initial values (inits
) are again specified in a
list. Valid tags for initial values now additionally include the
parameters associated with the spatial random effects. These include:
sigma.sq
(spatial variance parameter), phi
(spatial decay parameter), w
(the latent spatial random
effects at each site), and nu
(spatial smoothness
parameter), where the latter is only specified if adopting a Matern
covariance function (i.e., cov.model = 'matern'
).
spAbundance
supports four spatial covariance models
(exponential
, spherical
,
gaussian
, and matern
), which are specified in
the cov.model
argument. Throughout this vignette, we will
use an exponential covariance model, which we often use as our default
covariance model when fitting spatiallyexplicit models and is commonly
used throughout ecology. To determine which covariance function to use,
we can fit models with the different covariance functions and compare
them using WAIC to select the best performing function. We will note
that the Matern covariance function has the additional spatial
smoothness parameter \(\nu\) and thus
can often be more flexible than the other functions. However, because we
need to estimate an additional parameter, this also tends to require
more data (i.e., a larger number of sites) than the other covariance
functions, and so we encourage use of the three simpler functions if
your data set is small. We note that model estimates are generally
fairly robust to the different covariance functions, although certain
functions may provide substantially better estimates depending on the
specific form of the underlying spatial autocorrelation in the data. For
example, the Gaussian covariance function is often useful for accounting
for spatial autocorrelation that is very smooth (i.e., long range
spatial dependence). See Chapter 2 in Banerjee,
Carlin, and Gelfand (2003) for a more thorough discussion of
these functions and their mathematical properties.
The default initial values for phi
, and nu
are set to random values from the prior distribution, while the default
initial value for sigma.sq
is set to a random value between
0.05 and 3. In all spatiallyexplicit models described in this vignette,
the spatial decay parameter phi
is often the most sensitive
to initial values. In general, the spatial decay parameter will often
have poor mixing and take longer to converge than the rest of the
parameters in the model, so specifying an initial value that is
reasonably close to the resulting value can really help decrease run
times for complicated models. As an initial value for the spatial decay
parameter phi
, we compute the mean distance between points
in our coordinates matrix and then set it equal to 3 divided by this
mean distance. When using an exponential covariance function, \(\frac{3}{\phi}\) is the effective range, or
the distance at which the residual spatial correlation between two sites
drops to 0.05 (Banerjee, Carlin, and Gelfand
2003). Thus our initial guess for this effective range is the
average distance between sites across the simulated region. As with all
other parameters, we generally recommend using the default initial
values for an initial model run, and if the model is taking a very long
time to converge you can rerun the model with initial values based on
the posterior means of estimated parameters from the initial model fit.
For the spatial variance parameter sigma.sq
, we set the
initial value to 1. This corresponds to a moderate amount of spatial
variance. Further, we set the initial values of the latent spatial
random effects at each site to 0. The initial values for these random
effects has an extremely small influence on the model results, so we
generally recommend setting their initial values to 0 as we have done
here (this is also the default). However, if you are running your model
for a very long time and are seeing very slow convergence of the MCMC
chains, setting the initial values of the spatial random effects to the
mean estimates from a previous run of the model could help reach
convergence faster.
# Pairwise distances between all sites
dist.mat < dist(data.one.sp$coords)
# Exponential covariance model
cov.model < 'exponential'
# Specify list of inits
inits < list(alpha = 0,
beta = 0,
kappa = 0.5,
sigma.sq.mu = 0.5,
N = apply(data.one.sp$y, 1, max, na.rm = TRUE),
sigma.sq = 1,
phi = 3 / mean(dist.mat),
w = rep(0, nrow(data.one.sp$y)))
The parameter NNGP
is a logical value that specifies
whether we want to use an NNGP to fit the model. Currently, only
NNGP = TRUE
is supported, but we may eventually add
functionality to fit full Gaussian Process models. The arguments
n.neighbors
and search.type
specify the number
of neighbors used in the NNGP and the nearest neighbor search algorithm,
respectively, to use for the NNGP model. Generally, the default values
of these arguments will be adequate. Datta et al.
(2016) showed that setting n.neighbors = 15
is
usually sufficient, although for certain data sets a good approximation
can be achieved with as few as 5 neighbors, which could substantially
decrease run time for the model. We generally recommend leaving
search.type = "cb"
, as this results in a fast code book
nearest neighbor search algorithm. However, details on when you may want
to change this are described in Finley, Datta,
and Banerjee (2020). We will run an NNGP model using the default
value for search.type
and setting
n.neighbors = 15
(both the defaults).
NNGP < TRUE
n.neighbors < 15
search.type < 'cb'
Priors are again specified in a list in the argument
priors
. We follow standard recommendations for prior
distributions from the spatial statistics literature (Banerjee, Carlin, and Gelfand 2003). We assume
an inverse gamma prior for the spatial variance parameter
sigma.sq
(the tag of which is sigma.sq.ig
),
and uniform priors for the spatial decay parameter phi
and
smoothness parameter nu
(if using the Matern correlation
function), with the associated tags phi.unif
and
nu.unif
. The hyperparameters of the inverse Gamma are
passed as a vector of length two, with the first and second elements
corresponding to the shape and scale, respectively. The lower and upper
bounds of the uniform distribution are passed as a twoelement vector
for the uniform priors. We also allow users to restrict the spatial
variance further by specifying a uniform prior (with the tag
sigma.sq.unif
), which can potentially be useful to place a
more informative prior on the spatial parameters. Generally, we use an
inverseGamma prior.
Note that the priors for the spatial parameters in a
spatiallyexplicit model must be at least weakly informative for the
model to converge (Banerjee, Carlin, and Gelfand
2003). For the inverseGamma prior on the spatial variance, we
typically set the shape parameter to 2 and the scale parameter equal to
our best guess of the spatial variance. The default prior hyperparameter
values for the spatial variance \(\sigma^2\) are a shape parameter of 2 and a
scale parameter of 1. This weakly informative prior suggests a prior
mean of 1 for the spatial variance, which is a moderately small amount
of spatial variation. Here we will use this default prior. For the
spatial decay parameter, our default approach is to set the lower and
upper bounds of the uniform prior based on the minimum and maximum
distances between sites in the data. More specifically, by default we
set the lower bound to 3 / max
and the upper bound to
3 / min
, where min
and max
are
the minimum and maximum distances between sites in the data set,
respectively. This equates to a vague prior that states the spatial
autocorrelation in the data could only exist between sites that are very
close together, or could span across the entire observed study area. If
additional information is known on the extent of the spatial
autocorrelation in the data, you may place more restrictive bounds on
the uniform prior, which would reduce the amount of time needed for
adequate mixing and convergence of the MCMC chains. Here we use this
default approach, but will explicitly set the values for
transparency.
min.dist < min(dist.mat)
max.dist < max(dist.mat)
priors < list(alpha.normal = list(mean = 0, var = 2.72),
beta.normal = list(mean = 0, var = 100),
kappa.unif = c(0, 100),
sigma.sq.mu.ig = list(0.1, 0.1),
sigma.sq.ig = c(2, 1),
phi.unif = c(3 / max.dist, 3 / min.dist))
We again split our MCMC algorithm up into a set of batches and use an
adaptive sampler to adaptively tune the variances that we propose new
values from. We specify the initial tuning values again in the
tuning
argument, and now need to add phi
and
w
to the parameters that must be tuned. Note that we do not
need to add sigma.sq
, as this parameter can be sampled with
a more efficient approach (i.e., it’s full conditional distribution is
available in closed form).
tuning < list(beta = 0.5, alpha = 0.5, kappa = 0.5, beta.star = 0.5,
w = 0.5, phi = 0.5)
The argument n.omp.threads
specifies the number of
threads to use for withinchain parallelization, while
verbose
specifies whether or not to print the progress of
the sampler. As before, the argument n.report
specifies the
interval to report the MetropolisHastings sampler acceptance rate.
Below we set n.report = 400
, which will result in
information on the acceptance rate and tuning parameters every 400th
batch.
verbose < TRUE
batch.length < 25
n.batch < 1600
# Total number of MCMC samples per chain
batch.length * n.batch
[1] 40000
n.report < 400
n.omp.threads < 1
We will use the same amount of burnin and thinning as we did with
the nonspatial model, and we’ll also first fit a model with a Poisson
distribution for abundance. We next fit the model and summarize the
results using the summary()
function.
n.burn < 20000
n.thin < 20
n.chains < 3
# Approx run time: 3.5 min
out.sp < spNMix(abund.formula = abund.formula,
det.formula = det.formula,
data = data.one.sp,
inits = inits,
priors = priors,
n.batch = n.batch,
batch.length = batch.length,
tuning = tuning,
cov.model = cov.model,
NNGP = NNGP,
n.neighbors = n.neighbors,
search.type = search.type,
n.omp.threads = n.omp.threads,
n.report = n.report,
family = 'Poisson',
verbose = TRUE,
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 Poisson Nmixture model with 225 sites.
Samples per Chain: 40000 (1600 batches of length 25)
Burnin: 20000
Thinning Rate: 20
Number of Chains: 3
Total Posterior Samples: 3000
Using the exponential spatial correlation model.
Using 15 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: 400 of 1600, 25.00%
Parameter Acceptance Tuning
beta[1] 36.0 0.15832
beta[2] 44.0 0.15832
alpha[1] 32.0 0.26102
alpha[2] 32.0 0.27168
alpha[3] 32.0 0.32525
phi 12.0 0.52564

Batch: 800 of 1600, 50.00%
Parameter Acceptance Tuning
beta[1] 48.0 0.18579
beta[2] 48.0 0.16811
alpha[1] 36.0 0.28276
alpha[2] 44.0 0.26630
alpha[3] 16.0 0.34537
phi 36.0 0.59265

Batch: 1200 of 1600, 75.00%
Parameter Acceptance Tuning
beta[1] 40.0 0.17497
beta[2] 52.0 0.15518
alpha[1] 40.0 0.26102
alpha[2] 48.0 0.27168
alpha[3] 36.0 0.33853
phi 60.0 0.55814

Batch: 1600 of 1600, 100.00%

Chain 2

Sampling ...
Batch: 400 of 1600, 25.00%
Parameter Acceptance Tuning
beta[1] 64.0 0.18211
beta[2] 32.0 0.16811
alpha[1] 40.0 0.28276
alpha[2] 48.0 0.25585
alpha[3] 28.0 0.33853
phi 24.0 0.53625

Batch: 800 of 1600, 50.00%
Parameter Acceptance Tuning
beta[1] 36.0 0.16478
beta[2] 40.0 0.15211
alpha[1] 24.0 0.26630
alpha[2] 40.0 0.27716
alpha[3] 36.0 0.33853
phi 32.0 0.55814

Batch: 1200 of 1600, 75.00%
Parameter Acceptance Tuning
beta[1] 36.0 0.15832
beta[2] 24.0 0.16152
alpha[1] 64.0 0.28276
alpha[2] 48.0 0.28847
alpha[3] 44.0 0.31881
phi 16.0 0.60462

Batch: 1600 of 1600, 100.00%

Chain 3

Sampling ...
Batch: 400 of 1600, 25.00%
Parameter Acceptance Tuning
beta[1] 40.0 0.17497
beta[2] 36.0 0.15211
alpha[1] 44.0 0.30025
alpha[2] 64.0 0.28276
alpha[3] 48.0 0.35234
phi 28.0 0.64201

Batch: 800 of 1600, 50.00%
Parameter Acceptance Tuning
beta[1] 52.0 0.16152
beta[2] 48.0 0.15832
alpha[1] 48.0 0.27168
alpha[2] 20.0 0.28276
alpha[3] 44.0 0.37413
phi 28.0 0.66821

Batch: 1200 of 1600, 75.00%
Parameter Acceptance Tuning
beta[1] 28.0 0.16478
beta[2] 68.0 0.16152
alpha[1] 52.0 0.26102
alpha[2] 52.0 0.27716
alpha[3] 44.0 0.37413
phi 36.0 0.52564

Batch: 1600 of 1600, 100.00%
summary(out.sp)
Call:
spNMix(abund.formula = abund.formula, det.formula = det.formula,
data = data.one.sp, inits = inits, priors = priors, tuning = tuning,
cov.model = cov.model, NNGP = NNGP, n.neighbors = n.neighbors,
search.type = search.type, n.batch = n.batch, batch.length = batch.length,
family = "Poisson", n.omp.threads = n.omp.threads, verbose = TRUE,
n.report = n.report, n.burn = n.burn, n.thin = n.thin, n.chains = n.chains)
Samples per Chain: 40000
Burnin: 20000
Thinning Rate: 20
Number of Chains: 3
Total Posterior Samples: 3000
Run Time (min): 3.2288
Abundance (log scale):
Mean SD 2.5% 50% 97.5% Rhat ESS
(Intercept) 0.4751 0.2852 1.0373 0.4704 0.0704 1.0034 320
scale(abund.cov.1) 0.0877 0.1028 0.1172 0.0858 0.2856 1.0040 2201
Abundance Random Effect Variances (log scale):
Mean SD 2.5% 50% 97.5% Rhat ESS
(Intercept)abund.factor.1 0.3049 0.2694 0.0576 0.2299 0.9987 1.0051 2288
Detection (logit scale):
Mean SD 2.5% 50% 97.5% Rhat ESS
(Intercept) 0.3141 0.2066 0.0813 0.3180 0.7122 1.0202 836
scale(det.cov.1) 0.4902 0.1394 0.7657 0.4854 0.2261 1.0005 2400
scale(det.cov.2) 0.9302 0.1670 0.6205 0.9270 1.2619 1.0040 2202
Spatial Covariance:
Mean SD 2.5% 50% 97.5% Rhat ESS
sigma.sq 0.8130 0.2362 0.4363 0.7827 1.3424 1.0115 869
phi 13.9162 6.4785 5.0179 12.6938 31.5196 1.0114 380
Looking at the model summary we see adequate convergence of all model
parameters. The summary()
output looks the same as what we
saw previously, with the additional section titled “Spatial Covariance”.
There we see the estimate of the spatial variance
(sigma.sq
) and spatial decay (phi
) parameters.
Interpretation of the spatial variance parameter can follow
interpretation of variance parameters in “regular” (i.e., unstructured)
random effect variances: when the variance is close to 0, that indicates
little support for inclusion of the spatial random effect. When the
variance is large, it indicates substantial support for the spatial
random effect. What indicates “large” vs. “small” isn’t necessarily
straightforward. Here the estimate is about 0.81, which is fairly large
given that the value is on the log scale. We can also look at the
magnitude of the spatial random effect estimates themselves to give an
indication as to how much residual spatial autocorrelation there is in
the abundance estimates. The posterior samples for the spatial random
effects are stored in the w.samples
tag of the resulting
model fit list. Here we calculate the means of the spatial random
effects and plot a histogram of their values.
We see a decent amount of large magnitude random effect values, both above and below zero. If there was limited support for spatial autocorrelation, these values would all be very close to zero.
Posterior predictive checks
Posterior predictive checks proceed exactly before using the
ppcAbund()
function.
Call:
ppcAbund(object = out.sp, fit.stat = "freemantukey", group = 1)
Samples per Chain: 40000
Burnin: 20000
Thinning Rate: 20
Number of Chains: 3
Total Posterior Samples: 3000
Bayesian pvalue: 0.5097
Fit statistic: freemantukey
Here we see a striking contrast to the Bayesian pvalue from the
nonspatial Poisson Nmixture model, which was essentially 0. Here, our
estimate is very close to 0.5, indicating our model is adequately
representing the variability in the data with the addition of the
spatial random effect in abundance. If you take a look at the manual
page for dataNMixSim
, you can see how we simulated this
data set and why the above result makes sense.
Model selection using WAIC
We next compare the spatial Poisson Nmixture model to the best
performing nonspatial Nmixture model, which used a negative binomial
distribution (stored in out.2
).
# Negative binomial nonspatial model
waicAbund(out.2)
N.max not specified. Setting upper index of integration of N to 10 plus
the largest estimated abundance value at each site in object$N.samples
elpd pD WAIC
355.28485 15.51822 741.60615
# Poisson spatial model
waicAbund(out.sp)
N.max not specified. Setting upper index of integration of N to 10 plus
the largest estimated abundance value at each site in object$N.samples
elpd pD WAIC
295.77726 52.14108 695.83667
Here we see the spatial Poisson model outperforms the negative
binomial nonspatial model. Both of these models can be viewed as a
“standard” Nmixture model that accounts for overdispersion in two
different ways: the negative binomial model has a parameter
kappa
that controls the amount of overdispersion, which is
not assumed to have any spatial structure, while the spatial Poisson
model can be thought of as a Poissonlognormal model, where the
lognormal random effects are assigned a spatial structure. Here, the
latter model does a more adequate job of accounting for overdispersion
in the latent abundance values. We could also fit a spatial negative
binomial model, which essentially includes two mechanisms to account for
overdispersion. We leave this to the interested reader to fit this model
and see how it compares for this data set. We will note that spatial
negative binomial models can often be quite difficult to successfully
fit. This is particularly the case when there is no medium to large
range spatial autocorrelation in the abundance values. When this is the
case, the spatial random effects will be very difficult to identify from
the negative binomial overdispersion parameter, and we may seem some
unidentifiability occurring between the negative binomial overdispersion
parameter and the parameters controlling the spatial dependence. In such
a case, a spatial Poisson model or a negative binomial nonspatial model
will usually suffice. Alternatively, we could set an informative prior
on the spatial range to only allow the spatial random effects to explain
longrange (i.e., broadscale) spatial autocorrelation (if any exists)
and let the negative binomial overdispersion parameter account for
shortrange (i.e., finescale) overdispersion. We could then compare
that model to simpler versions of the model using WAIC.
Prediction
We can similarly predict across a region of interest using the
predict()
function as we saw with the nonspatial Nmixture
model. Here we again generate a marginal effects plot for the effect of
the covariate on abundance. When generating a marginal effects plot, we
set the include.sp
argument to FALSE
, which
means that we will generate predictions only using the covariates, their
effects, and any unstructured random effects included in the model. We
only recommend doing this if the goal is to produce a marginal effects
plot, as the inclusion of the spatial random effects often drastically
improves model predictive performance when interest lies in predicting
at new regions. See ?predict.spNMix
for an example of
predicting abundance at a set of new locations.
# Recall what the prediction design matrix looks like
str(X.0)
num [1:100, 1:2] 1 1 1 1 1 1 1 1 1 1 ...
 attr(*, "dimnames")=List of 2
..$ : NULL
..$ : chr [1:2] "intercept" "abund.cov.1"
List of 3
$ mu.0.samples: 'mcmc' num [1:3000, 1:100] 0.715 0.643 0.962 1.523 1.01 ...
.. attr(*, "mcpar")= num [1:3] 1 3000 1
$ N.0.samples : 'mcmc' int [1:3000, 1:100] 1 0 1 2 3 1 1 0 1 1 ...
.. attr(*, "mcpar")= num [1:3] 1 3000 1
$ call : language predict.spNMix(object = out.sp, X.0 = X.0, ignore.RE = TRUE, include.sp = FALSE)
 attr(*, "class")= chr "predict.spNMix"
# Get the lower bound, median, and 95% credible interval
mu.0.quants < apply(out.pred$mu.0.samples, 2, quantile,
prob = c(0.025, 0.5, 0.975))
mu.plot.dat < data.frame(mu.med = mu.0.quants[2, ],
mu.low = mu.0.quants[1, ],
mu.high = mu.0.quants[3, ],
abund.cov.1 = cov.pred.vals)
ggplot(mu.plot.dat, aes(x = abund.cov.1, y = mu.med)) +
geom_ribbon(aes(ymin = mu.low, ymax = mu.high), fill = 'grey70', alpha = 0.5) +
geom_line() +
theme_bw() +
labs(x = 'Covariate', y = 'Expected abundance')
Multispecies Nmixture models
Basic model description
Now consider the case where repeated count data, \(y_{i, j, k}\), are collected for multiple species \(i = 1, \dots, I\) at each survey location \(j\) during survey \(k\). We are now interested in estimating the abundance of each species \(i\) at each location \(j\), denoted as \(N_{i, j}\). We model \(N_{i, j}\) analogous to the singlespecies Nmixture model, with expected abundance now varying by species and site according to
\[\begin{equation}\label{mumsNMix} \text{log}(\mu_{i, j}) = \boldsymbol{x}_j^\top\boldsymbol{\beta}_i, \end{equation}\]
where \(\boldsymbol{\beta}_i\) are the speciesspecific effects of covariates \(\boldsymbol{x}_j\) (including an intercept) . When \(N_i(\boldsymbol{s}_j)\) is modeled using a negative binomial distribution, we estimate a separate dispersion parameter \(\kappa_i\) for each species. We model \(\boldsymbol{\beta}_i\) as random effects arising from a common, communitylevel normal distribution, which leads to increased precision of speciesspecific effects compared to singlespecies models (Yamaura et al. 2012). For example, the speciesspecific abundance intercept \(\beta0_i\) is modeled according to
\[\begin{equation}\label{betaComm} \beta0_i \sim \text{Normal}(\mu_{\beta0}, \tau^2_{\beta0}), \end{equation}\]
where \(\mu_{\beta0}\) is the communitylevel abundance intercept, and \(\tau^2_{\beta0}\) is the variance of the intercept across all \(I\) species. The observation portion of the multispecies Nmixture model is identical to the singlespecies model, with all parameters indexed by species, and the speciesspecific coefficients \(\boldsymbol{\alpha}_i\) are modeled hierarchically, analogous to the speciesspecific abundance coefficients just described.
We assign normal priors to the communitylevel abundance (\(\boldsymbol{\mu}_{\beta}\)) and detection (\(\boldsymbol{\mu}_{\alpha}\)) mean parameters and inverseGamma priors to the communitylevel variance parameters (\(\boldsymbol{\tau^2}_{\beta}\) and \(\boldsymbol{\tau^2}_{\alpha}\)). We give independent uniform priors to each of the speciesspecific negative binomial dispersion parameters \(\kappa_i\).
Fitting multispecies Nmixture models with
msNMix()
spAbundance
uses nearly identical syntax for fitting
multispecies Nmixture models as it does for singlespecies models and
provides the same functionality for posterior predictive checks, model
assessment and selection using WAIC, and prediction. The
msNMix()
function fits the multispecies Nmixture model
first introduced by Yamaura et al. (2012).
msNMix()
has the following syntax
msNMix(abund.formula, det.formula, data, inits, priors,
tuning, n.batch, batch.length, accept.rate = 0.43,
family = 'Poisson', n.omp.threads = 1, verbose = TRUE, n.report = 100,
n.burn = round(.10 * n.samples), n.thin = 1,
n.chains = 1, ...)
Notice these are the exact same arguments we saw with
NMix()
. We will again use our simulated data set contained
in dataNMixSim
, but now we will use data from all 6 species
in the simulated community. Below we reload the dataNMixSim
data set to get a fresh copy.
data(dataNMixSim)
We will model abundance for all species as a function of the
continuous covariate and a random effect of the categorical variable, as
well as the two continuous covariates for detection probability. For
multispecies models, the multispecies detectionnondetection data
y
is now a threedimensional array with dimensions
corresponding to species, sites, and replicates. This is how the data
are provided in the dataNMixSim
object, so we don’t need to
do any additional preparations. For guidance on preparing raw data into
such a threedimensional array format, please see this
vignette on the spOccupancy
website, which while
providing guidance for fitting multispecies models in
spOccupancy
, the format for spAbundance
is
exactly the same and all the same guidance applies here.
abund.ms.formula < ~ scale(abund.cov.1) + (1  abund.factor.1)
det.ms.formula < ~ scale(det.cov.1) + scale(det.cov.2)
str(dataNMixSim)
List of 4
$ y : int [1:6, 1:225, 1:3] 1 0 1 0 0 0 NA NA NA NA ...
$ abund.covs: num [1:225, 1:2] 0.373 0.706 0.202 1.588 0.138 ...
.. attr(*, "dimnames")=List of 2
.. ..$ : NULL
.. ..$ : chr [1:2] "abund.cov.1" "abund.factor.1"
$ det.covs :List of 2
..$ det.cov.1: num [1:225, 1:3] 1.28 NA NA NA 1.04 ...
..$ det.cov.2: num [1:225, 1:3] 2.03 NA NA NA 0.796 ...
$ coords : num [1:225, 1:2] 0 0.0714 0.1429 0.2143 0.2857 ...
.. attr(*, "dimnames")=List of 2
.. ..$ : NULL
.. ..$ : chr [1:2] "X" "Y"
Next we specify the initial values in inits
. For
multispecies Nmixture models, we supply initial values for
communitylevel and specieslevel parameters. In msNMix()
,
we will supply initial values for the following parameters:
alpha.comm
(communitylevel detection coefficients),
beta.comm
(communitylevel abundance coefficients),
alpha
(specieslevel detection coefficients),
beta
(specieslevel abundance coefficients),
tau.sq.beta
(communitylevel abundance variance
parameters), tau.sq.alpha
(communitylevel detection
variance parameters), N
(latent abundance variables for all
species), kappa
(specieslevel negative binomial
overdispersion parameters), sigma.sq.mu
(random effect
variances for abundance), and sigma.sq.p
(random effect
variances for detection). These are all specified in a single list.
Initial values for communitylevel parameters (including the random
effect variances) are either vectors of length corresponding to the
number of communitylevel detection or abundance parameters in the model
(including the intercepts) or a single value if all parameters are
assigned the same initial values. Initial values for species level
parameters are either matrices with the number of rows indicating the
number of species, and each column corresponding to a different
regression parameter, or a single value if the same initial value is
used for all species and parameters. The initial values for the latent
abundance matrix are specified as a matrix with \(I\) rows corresponding to the number of
species and \(J\) columns corresponding
to the number of sites.
# Number of species
n.sp < dim(dataNMixSim$y)[1]
ms.inits < list(alpha.comm = 0,
beta.comm = 0,
beta = 0,
alpha = 0,
tau.sq.beta = 1,
kappa = 1,
tau.sq.alpha = 1,
sigma.sq.mu = 0.5,
N = apply(dataNMixSim$y, c(1, 2), max, na.rm = TRUE))
In multispecies models, we specify priors on the communitylevel
coefficients (or hyperparameters) rather than the specieslevel effects,
with the exception that we still assign uniform priors for the negative
binomial overdispersion parameter. For nonspatial models, these priors
are specified with the following tags: beta.comm.normal
(normal prior on the communitylevel abundance mean effects),
alpha.comm.normal
(normal prior on the communitylevel
detection mean effects), tau.sq.beta.ig
(inverseGamma
prior on the communitylevel abundance variance parameters),
tau.sq.alpha.ig
(inverseGamma prior on the communitylevel
detection variance parameters), sigma.sq.mu.ig
(inverseGamma prior on the abundance random effect variances),
sigma.sq.p.ig
(inverseGamma prior on the detection random
effect variances), kappa.unif
(uniform prior on the
speciesspecific overdispersion parameters). For all parameters except
the speciesspecific overdispersion parameters, each tag consists of a
list with elements corresponding to the mean and variance for normal
priors and scale and shape for inverseGamma priors. Values can be
specified individually for each parameter or as a single value if the
same prior is assigned to all parameters of a given type. For the
speciesspecific overdispersion parameters, the prior is specified as a
list with two elements representing the lower and upper bound of the
uniform distribution, where each of the elements can be a single value
if the same prior is used for all species, or a vector of values for
each species if specifying a different prior for each species.
By default, we set the prior hyperparameter values for the
communitylevel means to a mean of 0 and a variance of 100 for abundance
and 2.72 for detection, which results in a vague prior on the
probability scale as we discussed for the singlespecies Nmixture
model. Below we specify these priors explicitly. For the communitylevel
variance parameters, by default we set the scale and shape parameters to
0.1 following the recommendations of (Lunn et al.
2013), which results in a weakly informative prior on the
communitylevel variances. This may lead to shrinkage of the
communitylevel variance towards zero under certain circumstances so
that all species will have fairly similar values for the
speciesspecific covariate effect (Gelman
2006), although we have found multispecies models to be
relatively robust to this prior specification. We are in the process of
implementing halfCauchy and halft priors for the variance parameters
in all multispecies models in spAbundance
(and
spOccupancy
), which generally result in less shrinkage
toward zero than the inverseGamma (Gelman
2006). For the negative binomial overdispersion parameters, we by
default use a lower bound of 0 and upper bound of 100 for the uniform
priors for each species, as we did in the singlespecies (univariate)
models. Below we explicitly set these default prior values for our
example.
ms.priors < list(beta.comm.normal = list(mean = 0, var = 100),
alpha.comm.normal = list(mean = 0, var = 2.72),
tau.sq.beta.ig = list(a = 0.1, b = 0.1),
tau.sq.alpha.ig = list(a = 0.1, b = 0.1),
sigma.sq.mu.ig = list(a = 0.1, b = 0.1),
kappa.unif = list(a = 0, b = 100))
All that’s now left to do is specify the number of threads to use
(n.omp.threads
), the number of MCMC samples (which we do by
specifying the number of batches n.batch
and batch length
batch.length
), the initial tuning variances
(tuning
), the amount of samples to discard as burnin
(n.burn
), the thinning rate (n.thin
), and
arguments to control the display of sampler progress
(verbose
, n.report
). Note for the tuning
variances, we do not need to specify initial tuning values for any of
the communitylevel parameters, as those parameters can be sampled with
an efficient Gibbs update. We will also run the model with a Poisson
distribution for abundance, which later we will shortly compare to a
negative binomial.
# Specify initial tuning values
ms.tuning < list(beta = 0.3, alpha = 0.3, beta.star = 0.5, kappa = 0.5)
# Approx. run time: 2.5 min
out.ms < msNMix(abund.formula = abund.ms.formula,
det.formula = det.ms.formula,
data = dataNMixSim,
inits = ms.inits,
n.batch = 1600,
tuning = ms.tuning,
batch.length = 25,
priors = ms.priors,
n.omp.threads = 1,
verbose = TRUE,
n.report = 400,
n.burn = 20000,
n.thin = 20,
n.chains = 3)

Preparing to run the model


Model description

Multispecies Poisson NMixture model fit with 225 sites and 6 species.
Samples per Chain: 40000
Burnin: 20000
Thinning Rate: 20
Number of Chains: 3
Total Posterior Samples: 3000
Source compiled with OpenMP support and model fit using 1 thread(s).
Adaptive Metropolis with target acceptance rate: 43.0

Chain 1

Sampling ...
Batch: 400 of 1600, 25.00%
Species Parameter Acceptance Tuning
1 beta[1] 40.0 0.16966
1 alpha[1] 40.0 0.25821
2 beta[1] 36.0 0.15047
2 alpha[1] 24.0 0.25310
3 beta[1] 28.0 0.13081
3 alpha[1] 24.0 0.21568
4 beta[1] 40.0 0.30914
4 alpha[1] 52.0 0.45205
5 beta[1] 44.0 0.10290
5 alpha[1] 44.0 0.18015
6 beta[1] 64.0 0.12569
6 alpha[1] 36.0 0.22003

Batch: 800 of 1600, 50.00%
Species Parameter Acceptance Tuning
1 beta[1] 40.0 0.17308
1 alpha[1] 32.0 0.29113
2 beta[1] 32.0 0.15047
2 alpha[1] 28.0 0.24318
3 beta[1] 36.0 0.11602
3 alpha[1] 72.0 0.21141
4 beta[1] 20.0 0.31538
4 alpha[1] 44.0 0.48000
5 beta[1] 36.0 0.11147
5 alpha[1] 40.0 0.18750
6 beta[1] 40.0 0.12320
6 alpha[1] 40.0 0.23836

Batch: 1200 of 1600, 75.00%
Species Parameter Acceptance Tuning
1 beta[1] 40.0 0.17308
1 alpha[1] 44.0 0.29113
2 beta[1] 60.0 0.15351
2 alpha[1] 44.0 0.27972
3 beta[1] 40.0 0.12320
3 alpha[1] 28.0 0.23836
4 beta[1] 48.0 0.34165
4 alpha[1] 28.0 0.45205
5 beta[1] 36.0 0.10710
5 alpha[1] 32.0 0.19910
6 beta[1] 48.0 0.12076
6 alpha[1] 44.0 0.23836

Batch: 1600 of 1600, 100.00%

Chain 2

Sampling ...
Batch: 400 of 1600, 25.00%
Species Parameter Acceptance Tuning
1 beta[1] 40.0 0.15978
1 alpha[1] 40.0 0.29113
2 beta[1] 48.0 0.14457
2 alpha[1] 52.0 0.24809
3 beta[1] 52.0 0.11372
3 alpha[1] 24.0 0.21568
4 beta[1] 40.0 0.29701
4 alpha[1] 40.0 0.50968
5 beta[1] 28.0 0.10086
5 alpha[1] 44.0 0.18015
6 beta[1] 32.0 0.12076
6 alpha[1] 56.0 0.24809

Batch: 800 of 1600, 50.00%
Species Parameter Acceptance Tuning
1 beta[1] 48.0 0.17308
1 alpha[1] 56.0 0.27418
2 beta[1] 48.0 0.15351
2 alpha[1] 52.0 0.27418
3 beta[1] 48.0 0.12076
3 alpha[1] 40.0 0.21568
4 beta[1] 48.0 0.31538
4 alpha[1] 52.0 0.48000
5 beta[1] 60.0 0.10498
5 alpha[1] 52.0 0.19129
6 beta[1] 28.0 0.11837
6 alpha[1] 36.0 0.23836

Batch: 1200 of 1600, 75.00%
Species Parameter Acceptance Tuning
1 beta[1] 48.0 0.17658
1 alpha[1] 60.0 0.26343
2 beta[1] 44.0 0.15351
2 alpha[1] 48.0 0.24809
3 beta[1] 48.0 0.11837
3 alpha[1] 44.0 0.23364
4 beta[1] 40.0 0.31538
4 alpha[1] 48.0 0.45205
5 beta[1] 56.0 0.10927
5 alpha[1] 52.0 0.19515
6 beta[1] 48.0 0.12822
6 alpha[1] 40.0 0.22448

Batch: 1600 of 1600, 100.00%

Chain 3

Sampling ...
Batch: 400 of 1600, 25.00%
Species Parameter Acceptance Tuning
1 beta[1] 56.0 0.16630
1 alpha[1] 44.0 0.28537
2 beta[1] 28.0 0.15351
2 alpha[1] 32.0 0.25310
3 beta[1] 56.0 0.12076
3 alpha[1] 28.0 0.22003
4 beta[1] 24.0 0.33488
4 alpha[1] 44.0 0.50968
5 beta[1] 40.0 0.10290
5 alpha[1] 40.0 0.19129
6 beta[1] 44.0 0.13081
6 alpha[1] 44.0 0.22901

Batch: 800 of 1600, 50.00%
Species Parameter Acceptance Tuning
1 beta[1] 56.0 0.16966
1 alpha[1] 48.0 0.27418
2 beta[1] 40.0 0.14749
2 alpha[1] 48.0 0.25821
3 beta[1] 40.0 0.12569
3 alpha[1] 60.0 0.23364
4 beta[1] 48.0 0.31538
4 alpha[1] 52.0 0.44309
5 beta[1] 40.0 0.10710
5 alpha[1] 48.0 0.19515
6 beta[1] 40.0 0.13081
6 alpha[1] 52.0 0.22003

Batch: 1200 of 1600, 75.00%
Species Parameter Acceptance Tuning
1 beta[1] 32.0 0.17658
1 alpha[1] 48.0 0.25310
2 beta[1] 48.0 0.13346
2 alpha[1] 32.0 0.24809
3 beta[1] 40.0 0.10710
3 alpha[1] 40.0 0.22901
4 beta[1] 76.0 0.32175
4 alpha[1] 48.0 0.49959
5 beta[1] 40.0 0.10710
5 alpha[1] 44.0 0.19515
6 beta[1] 48.0 0.12076
6 alpha[1] 40.0 0.24318

Batch: 1600 of 1600, 100.00%
The resulting object out.ms
is a list of class
msNMix
consisting primarily of posterior samples of all
community and specieslevel parameters, as well as some additional
objects that are used for summaries, prediction, and model fit
evaluation. We can display a nice summary of these results using the
summary()
function as before. For multispecies objects,
when using summary we need to specify the level of parameters we want to
summarize. We do this using the argument level
, which takes
values community
, species
, or
both
to print results for communitylevel parameters,
specieslevel parameters, or all parameters. By default,
level
is set to both
to display both species
and communitylevel parameters.
summary(out.ms)
Call:
msNMix(abund.formula = abund.ms.formula, det.formula = det.ms.formula,
data = dataNMixSim, inits = ms.inits, priors = ms.priors,
tuning = ms.tuning, n.batch = 1600, batch.length = 25, n.omp.threads = 1,
verbose = TRUE, n.report = 400, n.burn = 20000, n.thin = 20,
n.chains = 3)
Samples per Chain: 40000
Burnin: 20000
Thinning Rate: 20
Number of Chains: 3
Total Posterior Samples: 3000
Run Time (min): 2.4987

Community Level

Abundance Means (log scale):
Mean SD 2.5% 50% 97.5% Rhat ESS
(Intercept) 0.3867 0.4782 1.3236 0.3758 0.5399 1.0009 2439
scale(abund.cov.1) 0.1353 0.5066 1.1359 0.1392 0.9583 1.0007 3212
Abundance Variances (log scale):
Mean SD 2.5% 50% 97.5% Rhat ESS
(Intercept) 1.2619 2.4750 0.2251 0.8321 4.6964 1.1792 3000
scale(abund.cov.1) 1.5388 1.8009 0.3465 1.0757 5.6716 1.0039 3000
Abundance Random Effect Variances (log scale):
Mean SD 2.5% 50% 97.5% Rhat ESS
(Intercept)abund.factor.1 0.6758 0.1574 0.4275 0.6574 1.0426 1.0034 1423
Detection Means (logit scale):
Mean SD 2.5% 50% 97.5% Rhat ESS
(Intercept) 0.5682 0.1859 0.1799 0.5682 0.9245 1.0015 2400
scale(det.cov.1) 0.1100 0.5416 0.9836 0.1060 1.1858 0.9997 2614
scale(det.cov.2) 0.9698 0.5128 0.1450 0.9919 1.9410 1.0015 3000
Detection Variances (logit scale):
Mean SD 2.5% 50% 97.5% Rhat ESS
(Intercept) 0.1702 0.2013 0.0318 0.1157 0.6221 1.0160 2320
scale(det.cov.1) 2.1259 2.1278 0.5435 1.5253 7.7415 1.0074 3188
scale(det.cov.2) 1.8679 1.9805 0.3887 1.2942 6.9378 1.0039 3000

Species Level

Abundance (log scale):
Mean SD 2.5% 50% 97.5% Rhat ESS
(Intercept)sp1 0.2086 0.2625 0.7255 0.2131 0.3059 1.0043 501
(Intercept)sp2 0.0168 0.2640 0.5179 0.0146 0.5018 1.0142 311
(Intercept)sp3 0.1312 0.2577 0.4049 0.1353 0.6226 1.0378 236
(Intercept)sp4 1.9118 0.3648 2.6323 1.9147 1.1928 1.0021 971
(Intercept)sp5 0.3397 0.2888 0.2583 0.3546 0.8691 1.0875 137
(Intercept)sp6 0.6901 0.2789 1.2225 0.6967 0.1459 1.0082 263
scale(abund.cov.1)sp1 0.1173 0.0816 0.0421 0.1178 0.2759 1.0062 2744
scale(abund.cov.1)sp2 0.5972 0.0731 0.7412 0.5963 0.4532 0.9996 2714
scale(abund.cov.1)sp3 0.2619 0.0559 0.1500 0.2608 0.3748 1.0074 2712
scale(abund.cov.1)sp4 1.1474 0.1869 1.5369 1.1390 0.7949 1.0008 1965
scale(abund.cov.1)sp5 0.9804 0.0578 1.0931 0.9794 0.8670 1.0035 1597
scale(abund.cov.1)sp6 1.4246 0.0633 1.2979 1.4259 1.5509 1.0043 1442
Detection (logit scale):
Mean SD 2.5% 50% 97.5% Rhat ESS
(Intercept)sp1 0.4861 0.1627 0.1615 0.4848 0.8003 1.0137 1637
(Intercept)sp2 0.8167 0.1771 0.4846 0.8138 1.1693 1.0015 1451
(Intercept)sp3 0.2626 0.1543 0.0484 0.2655 0.5563 1.0010 1342
(Intercept)sp4 0.6428 0.2706 0.1216 0.6396 1.1957 1.0011 1923
(Intercept)sp5 0.7079 0.1402 0.4405 0.7081 0.9987 1.0009 1224
(Intercept)sp6 0.5224 0.1700 0.1747 0.5247 0.8492 1.0010 1189
scale(det.cov.1)sp1 0.5343 0.1351 0.8155 0.5316 0.2786 1.0030 2600
scale(det.cov.1)sp2 0.0078 0.1297 0.2558 0.0102 0.2458 1.0013 2621
scale(det.cov.1)sp3 0.7548 0.1257 0.5187 0.7545 1.0042 1.0095 2687
scale(det.cov.1)sp4 1.0724 0.2626 0.5941 1.0656 1.6245 1.0052 3000
scale(det.cov.1)sp5 1.8104 0.1473 2.1086 1.8090 1.5354 1.0010 1312
scale(det.cov.1)sp6 1.3690 0.1620 1.0602 1.3693 1.6920 1.0051 1660
scale(det.cov.2)sp1 1.0053 0.1674 0.6791 1.0019 1.3441 1.0021 2534
scale(det.cov.2)sp2 1.1255 0.1621 0.8155 1.1228 1.4495 1.0122 2402
scale(det.cov.2)sp3 1.6414 0.1527 1.3548 1.6379 1.9530 1.0004 2360
scale(det.cov.2)sp4 0.8272 0.3621 1.5895 0.8104 0.1528 1.0010 2583
scale(det.cov.2)sp5 1.1793 0.1176 0.9523 1.1774 1.4133 1.0038 1832
scale(det.cov.2)sp6 2.3820 0.2351 1.9330 2.3796 2.8453 1.0003 1689
# Or more explicitly
# summary(out.ms, level = 'both')
We see adequate convergence and mixing of all model parameters, with ESS values fairly large and Rhat values less than 1.1. Looking at the communitylevel variance parameters, we can see there is variation in average abundance of each species (i.e., the intercept) as well as the covariate effect. Looking at the speciesspecific covariate effects, we can see the effects range from highly positive to highly negative.
Posterior predictive checks
As with singlespecies models, we can use the ppcAbund()
function to perform posterior predictive checks, and summarize the check
with a Bayesian pvalue using the summary()
function. The
summary()
function again requires the level
argument to specify if you want an overall Bayesian pvalue for the
entire community (level = 'community'
), each individual
species (level = 'species'
), or for both
(level = 'both'
). By default, we set `level = ‘both’. This
time, we use a chisquared test statistic, grouping the data by
replicates (again, we just do this to display the different ways of
performing posterior predictive checks, and for a full analysis we would
explore multiple fit statistics and ways of grouping the data).
ppc.ms.out < ppcAbund(out.ms, fit.stat = 'chisquared', group = 2)
Currently on species 1 out of 6
Currently on species 2 out of 6
Currently on species 3 out of 6
Currently on species 4 out of 6
Currently on species 5 out of 6
Currently on species 6 out of 6
summary(ppc.ms.out)
Call:
ppcAbund(object = out.ms, fit.stat = "chisquared", group = 2)
Samples per Chain: 40000
Burnin: 20000
Thinning Rate: 20
Number of Chains: 3
Total Posterior Samples: 3000

Community Level

Bayesian pvalue: 0.2651

Species Level

sp1 Bayesian pvalue: 0.3497
sp2 Bayesian pvalue: 0.0897
sp3 Bayesian pvalue: 3e04
sp4 Bayesian pvalue: 0.459
sp5 Bayesian pvalue: 0.268
sp6 Bayesian pvalue: 0.4237
Fit statistic: chisquared
Here we see mixed results of the posterior predictive check, with the third species having a very small Bayesian pvalue, suggesting the model does not fit very well for that species, which as with the singlespecies case, is likely related to additional spatial variation in abundance that we are not accounting for.
Model selection using WAIC
We can compute the WAIC for comparison with alternative models using
the waicAbund()
function. For multispecies models, we can
calculate the WAIC for the entire data set, or we can calculate it
separately for each species. This is done by using the logical
by.sp
argument. Note that the WAIC for the entire data set
is simply the sum of all the WAIC values for the individual responses.
Below we calculate the WAIC by species for the multispecies model and
compare the WAIC for species 1 to the corresponding WAIC for the
nonspatial singlespecies model we fit for that species.
# Multispecies Poisson Nmixture model
waicAbund(out.ms, by.sp = TRUE)
N.max not specified. Setting upper index of integration of N to 10 plus
the largest estimated abundance value at each site in object$N.samples
Currently on species 1 out of 6
Currently on species 2 out of 6
Currently on species 3 out of 6
Currently on species 4 out of 6
Currently on species 5 out of 6
Currently on species 6 out of 6
elpd pD WAIC
1 368.2380 23.62859 783.7331
2 466.3158 54.39663 1041.4248
3 564.2021 71.35770 1271.1197
4 118.1010 12.09634 260.3946
5 532.5067 37.83546 1140.6842
6 295.3850 27.93759 646.6451
# Singlespecies Poisson Nmixture model for species 1
waicAbund(out)
N.max not specified. Setting upper index of integration of N to 10 plus
the largest estimated abundance value at each site in object$N.samples
elpd pD WAIC
369.16550 22.74705 783.82511
Here we see nearly identical WAIC values for species 1 across a singlespecies model vs. a multispecies model. Next, we fit a multispecies negative binomial Nmixture model for comparison with the Poisson multispecies model
# Approx. run time: 2.5 min
out.ms.nb < msNMix(abund.formula = abund.ms.formula,
det.formula = det.ms.formula,
data = dataNMixSim,
inits = ms.inits,
n.batch = 1600,
tuning = ms.tuning,
batch.length = 25,
priors = ms.priors,
n.omp.threads = 1,
verbose = FALSE,
family = 'NB',
n.burn = 20000,
n.thin = 20,
n.chains = 3)
# Negative binomial
waicAbund(out.ms.nb, by.sp = TRUE)
N.max not specified. Setting upper index of integration of N to 10 plus
the largest estimated abundance value at each site in object$N.samples
Currently on species 1 out of 6
Currently on species 2 out of 6
Currently on species 3 out of 6
Currently on species 4 out of 6
Currently on species 5 out of 6
Currently on species 6 out of 6
elpd pD WAIC
1 354.1900 15.69356 739.7672
2 390.3124 17.27289 815.1707
3 399.1375 13.82221 825.9193
4 118.4489 11.56989 260.0377
5 480.0320 13.47496 987.0139
6 276.1147 11.86407 575.9575
# Poisson
waicAbund(out.ms, by.sp = TRUE)
N.max not specified. Setting upper index of integration of N to 10 plus
the largest estimated abundance value at each site in object$N.samples
Currently on species 1 out of 6
Currently on species 2 out of 6
Currently on species 3 out of 6
Currently on species 4 out of 6
Currently on species 5 out of 6
Currently on species 6 out of 6
elpd pD WAIC
1 368.2380 23.62859 783.7331
2 466.3158 54.39663 1041.4248
3 564.2021 71.35770 1271.1197
4 118.1010 12.09634 260.3946
5 532.5067 37.83546 1140.6842
6 295.3850 27.93759 646.6451
Here we see the negative binomial model provides substantial improvements for most species in the data set.
Prediction
Prediction proceeds exactly as we have seen previously with the singlespecies nonspatial Nmixture model. Below we generate a marginal effects plot for the effect of the covariate on each of the six species in the data set. Because we saw substantial improvements in WAIC for the first two species in the data set, we will use the negative binomial model
List of 3
$ mu.0.samples: num [1:3000, 1:6, 1:100] 0.39 0.61 0.514 0.49 0.793 ...
$ N.0.samples : num [1:3000, 1:6, 1:100] 0 0 0 1 0 1 1 0 0 0 ...
$ call : language predict.msNMix(object = out.ms.nb, X.0 = X.0, ignore.RE = TRUE)
 attr(*, "class")= chr "predict.msNMix"
Notice that mu.0.samples
(expected abundances) and
N.0.samples
(mean abundances) are now threedimensional
arrays, with dimensions corresponding to posterior sample, species, and
site. Now, we make a marginal effect plot for each species.
mu.0.quants < apply(out.pred$mu.0.samples, c(2, 3), quantile,
prob = c(0.025, 0.5, 0.975))
# Put stuff together for plotting with ggplot2
mu.0.plot.df < data.frame(mu.med = c(mu.0.quants[2, , ]),
mu.low = c(mu.0.quants[1, , ]),
mu.high = c(mu.0.quants[3, , ]),
sp = rep(1:n.sp, nrow(X.0)),
abund.cov.1 = rep(cov.pred.vals, each = n.sp))
ggplot(mu.0.plot.df, aes(x = abund.cov.1, y = mu.med)) +
geom_ribbon(aes(ymin = mu.low, ymax = mu.high), fill = 'grey70') +
geom_line() +
theme_bw() +
facet_wrap(vars(sp), scales = 'free_y') +
labs(x = 'Covariate', y = 'Expected abundance')
Latent factor multispecies Nmixture models
Basic model description
The latent factor multispecies Nmixture model is identical to the previously described multispecies Nmixture model except we include an additional component in the model to accommodate residual correlations between species. This type of model is often referred to as an abundancebased joint species distribution model (JSDM) in the statistical ecology literature (Warton et al. 2015). The previously described multispecies Nmixture model can be viewed as a simplified version of the latent factor multispecies Nmixture model, where we assume there are no residual correlations between species. This model was described in Section 8.5.4 of Kéry and Royle (2021), which they implemented in JAGS.
The model is identical to the previously described multispecies Nmixture model, with the only addition being a speciesspecific random effect at each site added to the equation for expected abundance. More specifically, we model speciesspecific abundance as
\[\begin{equation} \text{log}(\mu_{i, j}) = \boldsymbol{x}_j^\top\boldsymbol{\beta}_i + \text{w}^\ast_{i, j}. \end{equation}\]
The speciesspecific random effect \(\text{w}^\ast_{i, j}\) is used to account for residual correlations between species by assuming that correlation amongst the species can be explained by speciesspecific effects of a set of \(q\) latent variables. More specifically, we use a factor modeling approach (Lopes and West 2004) to account for species residual correlations in a computationally efficient manner (Tobler et al. 2019). This approach is ideal for large groups of species, where estimating a full \(I \times I\) covariance matrix as done in Mimnagh et al. (2022) would be computationally intractable (and/or require massive amounts of data). Specifically, we decompose \(\text{w}^\ast_{i, j}\) into a linear combination of \(q\) latent variables (i.e., factors) and their associated speciesspecific coefficients (i.e., factor loadings) according to
\[\begin{equation}\label{factorModel} \text{w}^\ast_{i, j} = \boldsymbol{\lambda}_i^\top\text{w}_j, \end{equation}\]
where \(\boldsymbol{\lambda}_i^\top\) is the \(i\) row of factor loadings from an \(I \times q\) loadings matrix \(\boldsymbol{\Lambda}\), and \(\text{w}_j\) is a \(q \times 1\) vector of independent latent factors at site \(j\). By setting \(q << N\), we achieve dimension reduction to efficiently model communities with a large number of species (TaylorRodriguez et al. 2019; Tobler et al. 2019; Doser, Finley, and Banerjee 2023). The approach accounts for residual species correlations via their speciesspecific responses to the \(q\) factors. We model each latent factor as a standard normal random variable. To ensure identifiability of the latent factors from the latent factor loadings, we fix the upper triangle of the factor loadings matrix to 0 and the diagonal elements to 1. We assign standard normal priors to the lower triangular elements of the factor loadings matrix. All other priors are identical to the multispecies Nmixture model previously discussed.
The constraints on the factor loadings matrix ensures identifiability
of the factor loadings from the latent factors, but this also results in
important practical considerations when fitting these models (e.g.,
ordering of the species, initial values). This
vignette on the spOccupancy
website discusses these
(and other) considerations. All the advice applied to factor models fit
with detectionnondetection data in spOccupancy
also
applies to factor models fit in spAbundance
.
Fitting latent factor multispecies Nmixture models with
lfMsNMix()
The function lfMsNMix()
fits latent factor multispecies
Nmixture models in spAbundance
. The arguments are
identical to those in msNMix()
with one additional argument
that specifies the number of factors we wish to use in the model
(n.factors
):
lfMsNMix(abund.formula, det.formula, data, inits, priors,
tuning, n.factors, n.batch, batch.length, accept.rate = 0.43,
family = 'Poisson', n.omp.threads = 1, verbose = TRUE, n.report = 100,
n.burn = round(.10 * n.samples), n.thin = 1,
n.chains = 1, ...)
For guidance on choosing the number of latent factors, see this
vignette on the spOccupancy
website. Here we will fit a
model with 3 latent factors.
n.factors < 3
There are only a few slight differences in how we go about fitting a
multispecies Nmixture model with latent factors compared to one
without latent factors. The data
format as well as
abund.formula
and det.formula
remain the same
as before.
abund.ms.formula < ~ scale(abund.cov.1) + (1  abund.factor.1)
det.ms.formula < ~ scale(det.cov.1) + scale(det.cov.2)
str(dataNMixSim)
List of 4
$ y : int [1:6, 1:225, 1:3] 1 0 1 0 0 0 NA NA NA NA ...
$ abund.covs: num [1:225, 1:2] 0.373 0.706 0.202 1.588 0.138 ...
.. attr(*, "dimnames")=List of 2
.. ..$ : NULL
.. ..$ : chr [1:2] "abund.cov.1" "abund.factor.1"
$ det.covs :List of 2
..$ det.cov.1: num [1:225, 1:3] 1.28 NA NA NA 1.04 ...
..$ det.cov.2: num [1:225, 1:3] 2.03 NA NA NA 0.796 ...
$ coords : num [1:225, 1:2] 0 0.0714 0.1429 0.2143 0.2857 ...
.. attr(*, "dimnames")=List of 2
.. ..$ : NULL
.. ..$ : chr [1:2] "X" "Y"
Initial values are specified as with msNMix()
, with the
exception that we can now specify initial values for the latent factor
loadings matrix lambda
and the latent factors
w
. Initial values for the speciesspecific factor loadings
(lambda
) are specified as a numeric matrix with \(I\) rows and \(q\) columns, where \(I\) is the number of species and \(q\) is the number of latent factors used in
the model. The diagonal elements of the matrix must be 1, and values in
the upper triangle must be set to 0 to ensure identifiability of the
latent factors. Note that the default initial values for the lower
triangle of the factor loadings matrix is 0. We can also specify the
initial values for the latent factors. Below we set these to 0 (which is
the default).
# Number of species
n.sp < dim(dataNMixSim$y)[1]
# Initiate all lambda initial values to 0.
lambda.inits < matrix(0, n.sp, n.factors)
# Set diagonal elements to 1
diag(lambda.inits) < 1
# Set lower triangular elements to random values from a standard normal distribution
lambda.inits[lower.tri(lambda.inits)] < rnorm(sum(lower.tri(lambda.inits)))
# Check it out.
lambda.inits
[,1] [,2] [,3]
[1,] 1.0000000 0.0000000 0.0000000
[2,] 0.4466988 1.0000000 0.0000000
[3,] 1.7183729 1.2129439 1.0000000
[4,] 1.9870724 0.3601976 0.4133971
[5,] 1.3758751 0.1253551 1.1161277
[6,] 0.2866810 0.8837145 0.7653206
w.inits < matrix(0, n.factors, ncol(dataNMixSim$y))
ms.inits < list(alpha.comm = 0,
beta.comm = 0,
beta = 0,
alpha = 0,
tau.sq.beta = 1,
lambda = lambda.inits,
kappa = 1,
w = w.inits,
tau.sq.alpha = 1,
sigma.sq.mu = 0.5,
N = apply(dataNMixSim$y, c(1, 2), max, na.rm = TRUE))
Priors are specified exactly the same as we saw with
msNMix()
. We do not need to explicitly specify the prior
for the factor loadings, as we require the prior for the
lowertriangular values to be Normal(0, 1).
ms.priors < list(beta.comm.normal = list(mean = 0, var = 2.72),
alpha.comm.normal = list(mean = 0, var = 2.72),
tau.sq.beta.ig = list(a = 0.1, b = 0.1),
tau.sq.alpha.ig = list(a = 0.1, b = 0.1),
sigma.sq.mu.ig = list(a = 0.1, b = 0.1),
kappa.unif = list(a = 0, b = 100))
We finish up by specifying the tuning values, where now we specify
the initial tuning variance for the factor loadings lambda
as well as the latent factors w
. We then run the model for
40,000 iterations with a burnin of 20,000 samples and a thinning rate
of 20, for each of 3 chains, yielding a total of 3000 posterior samples.
We will fit the model with a Poisson distribution for abundance.
# Specify initial tuning values
ms.tuning < list(beta = 0.3, alpha = 0.3, beta.star = 0.5, kappa = 0.5, w = 0.5, lambda = 0.5)
# Approx. run time: ~4.5 min
out.lf.ms < lfMsNMix(abund.formula = abund.ms.formula,
det.formula = det.ms.formula,
data = dataNMixSim,
inits = ms.inits,
n.batch = 1600,
n.factors = n.factors,
tuning = ms.tuning,
batch.length = 25,
priors = ms.priors,
n.omp.threads = 1,
verbose = TRUE,
n.report = 400,
n.burn = 20000,
n.thin = 20,
n.chains = 3)

Preparing to run the model


Model description

Latent Factor Multispecies Poisson NMixture model fit with 225 sites and 6 species.
Samples per Chain: 40000
Burnin: 20000
Thinning Rate: 20
Number of Chains: 3
Total Posterior Samples: 3000
Using 3 latent factors.
Source compiled with OpenMP support and model fit using 1 thread(s).
Adaptive Metropolis with target acceptance rate: 43.0

Chain 1

Sampling ...
Batch: 400 of 1600, 25.00%
Species Parameter Acceptance Tuning
1 beta[1] 48.0 0.17308
1 alpha[1] 24.0 0.24318
2 beta[1] 36.0 0.13890
2 alpha[1] 40.0 0.21568
3 beta[1] 48.0 0.11372
3 alpha[1] 32.0 0.21568
4 beta[1] 36.0 0.32175
4 alpha[1] 52.0 0.42572
5 beta[1] 44.0 0.10290
5 alpha[1] 40.0 0.17658
6 beta[1] 44.0 0.12320
6 alpha[1] 44.0 0.25821

Batch: 800 of 1600, 50.00%
Species Parameter Acceptance Tuning
1 beta[1] 24.0 0.16630
1 alpha[1] 48.0 0.24318
2 beta[1] 52.0 0.13890
2 alpha[1] 36.0 0.23364
3 beta[1] 32.0 0.11372
3 alpha[1] 52.0 0.21568
4 beta[1] 36.0 0.30302
4 alpha[1] 44.0 0.46118
5 beta[1] 48.0 0.10290
5 alpha[1] 56.0 0.18750
6 beta[1] 44.0 0.13346
6 alpha[1] 28.0 0.24318

Batch: 1200 of 1600, 75.00%
Species Parameter Acceptance Tuning
1 beta[1] 36.0 0.16301
1 alpha[1] 44.0 0.26343
2 beta[1] 64.0 0.13081
2 alpha[1] 68.0 0.22901
3 beta[1] 48.0 0.11372
3 alpha[1] 40.0 0.22003
4 beta[1] 28.0 0.30302
4 alpha[1] 48.0 0.47049
5 beta[1] 52.0 0.10290
5 alpha[1] 44.0 0.19515
6 beta[1] 24.0 0.12569
6 alpha[1] 36.0 0.21568

Batch: 1600 of 1600, 100.00%

Chain 2

Sampling ...
Batch: 400 of 1600, 25.00%
Species Parameter Acceptance Tuning
1 beta[1] 48.0 0.17308
1 alpha[1] 60.0 0.24318
2 beta[1] 44.0 0.13615
2 alpha[1] 32.0 0.21568
3 beta[1] 32.0 0.12076
3 alpha[1] 40.0 0.21568
4 beta[1] 40.0 0.30914
4 alpha[1] 28.0 0.47049
5 beta[1] 28.0 0.10086
5 alpha[1] 60.0 0.16966
6 beta[1] 36.0 0.11372
6 alpha[1] 56.0 0.22901

Batch: 800 of 1600, 50.00%
Species Parameter Acceptance Tuning
1 beta[1] 36.0 0.17658
1 alpha[1] 60.0 0.24318
2 beta[1] 44.0 0.14171
2 alpha[1] 40.0 0.22901
3 beta[1] 40.0 0.11602
3 alpha[1] 32.0 0.22003
4 beta[1] 64.0 0.32175
4 alpha[1] 32.0 0.48000
5 beta[1] 40.0 0.09691
5 alpha[1] 36.0 0.18015
6 beta[1] 40.0 0.11837
6 alpha[1] 40.0 0.24318

Batch: 1200 of 1600, 75.00%
Species Parameter Acceptance Tuning
1 beta[1] 28.0 0.16630
1 alpha[1] 48.0 0.26875
2 beta[1] 40.0 0.14171
2 alpha[1] 44.0 0.22003
3 beta[1] 52.0 0.11837
3 alpha[1] 32.0 0.22448
4 beta[1] 48.0 0.33488
4 alpha[1] 36.0 0.44309
5 beta[1] 56.0 0.10290
5 alpha[1] 32.0 0.17658
6 beta[1] 36.0 0.12076
6 alpha[1] 36.0 0.21568

Batch: 1600 of 1600, 100.00%

Chain 3

Sampling ...
Batch: 400 of 1600, 25.00%
Species Parameter Acceptance Tuning
1 beta[1] 52.0 0.16966
1 alpha[1] 32.0 0.26343
2 beta[1] 36.0 0.13890
2 alpha[1] 40.0 0.22901
3 beta[1] 52.0 0.12320
3 alpha[1] 48.0 0.22003
4 beta[1] 24.0 0.33488
4 alpha[1] 24.0 0.44309
5 beta[1] 40.0 0.09691
5 alpha[1] 48.0 0.19129
6 beta[1] 32.0 0.11602
6 alpha[1] 20.0 0.24318

Batch: 800 of 1600, 50.00%
Species Parameter Acceptance Tuning
1 beta[1] 40.0 0.16966
1 alpha[1] 40.0 0.25821
2 beta[1] 28.0 0.14457
2 alpha[1] 36.0 0.23836
3 beta[1] 36.0 0.11372
3 alpha[1] 40.0 0.23364
4 beta[1] 44.0 0.29701
4 alpha[1] 40.0 0.47049
5 beta[1] 28.0 0.09499
5 alpha[1] 48.0 0.18379
6 beta[1] 52.0 0.13081
6 alpha[1] 44.0 0.24809

Batch: 1200 of 1600, 75.00%
Species Parameter Acceptance Tuning
1 beta[1] 44.0 0.16966
1 alpha[1] 56.0 0.24809
2 beta[1] 28.0 0.14171
2 alpha[1] 56.0 0.21568
3 beta[1] 36.0 0.12320
3 alpha[1] 60.0 0.23364
4 beta[1] 40.0 0.32825
4 alpha[1] 44.0 0.48000
5 beta[1] 48.0 0.10086
5 alpha[1] 44.0 0.17658
6 beta[1] 36.0 0.11372
6 alpha[1] 36.0 0.23836

Batch: 1600 of 1600, 100.00%
summary(out.lf.ms)
Call:
lfMsNMix(abund.formula = abund.ms.formula, det.formula = det.ms.formula,
data = dataNMixSim, inits = ms.inits, priors = ms.priors,
tuning = ms.tuning, n.factors = n.factors, n.batch = 1600,
batch.length = 25, n.omp.threads = 1, verbose = TRUE, n.report = 400,
n.burn = 20000, n.thin = 20, n.chains = 3)
Samples per Chain: 40000
Burnin: 20000
Thinning Rate: 20
Number of Chains: 3
Total Posterior Samples: 3000
Run Time (min): 4.319

Community Level

Abundance Means (log scale):
Mean SD 2.5% 50% 97.5% Rhat ESS
(Intercept) 0.8825 0.3938 1.6756 0.8837 0.0535 1.0037 2173
scale(abund.cov.1) 0.1509 0.4788 1.0818 0.1584 0.8438 1.0072 3000
Abundance Variances (log scale):
Mean SD 2.5% 50% 97.5% Rhat ESS
(Intercept) 0.9806 1.1568 0.1689 0.6804 3.8457 1.0197 2411
scale(abund.cov.1) 1.6222 1.6289 0.4050 1.1464 5.6940 1.0013 3000
Abundance Random Effect Variances (log scale):
Mean SD 2.5% 50% 97.5% Rhat ESS
(Intercept)abund.factor.1 0.4554 0.1332 0.2529 0.4354 0.775 1.0022 809
Detection Means (logit scale):
Mean SD 2.5% 50% 97.5% Rhat ESS
(Intercept) 0.3889 0.1798 0.0425 0.3882 0.7502 1.0068 1632
scale(det.cov.1) 0.0832 0.4880 0.8937 0.0831 1.0988 0.9998 2135
scale(det.cov.2) 0.9462 0.5132 0.1472 0.9713 1.9332 1.0008 3000
Detection Variances (logit scale):
Mean SD 2.5% 50% 97.5% Rhat ESS
(Intercept) 0.1485 0.1752 0.0289 0.0983 0.5856 1.0063 2178
scale(det.cov.1) 1.7171 1.7784 0.4326 1.2601 5.6998 1.0208 3000
scale(det.cov.2) 1.8326 1.7815 0.3998 1.3052 6.6865 1.0107 3000

Species Level

Abundance (log scale):
Mean SD 2.5% 50% 97.5% Rhat ESS
(Intercept)sp1 0.6103 0.2461 1.1064 0.6098 0.1426 1.0216 542
(Intercept)sp2 0.7792 0.2611 1.3163 0.7715 0.2763 1.0327 359
(Intercept)sp3 0.7493 0.2516 1.2543 0.7542 0.2437 1.0234 316
(Intercept)sp4 2.1805 0.3800 2.9564 2.1779 1.4495 1.0135 906
(Intercept)sp5 0.0434 0.2493 0.5470 0.0278 0.4239 1.0645 211
(Intercept)sp6 1.2122 0.2702 1.7537 1.2104 0.6697 1.0482 275
scale(abund.cov.1)sp1 0.0911 0.1133 0.1349 0.0895 0.3161 1.0016 1872
scale(abund.cov.1)sp2 0.6498 0.1459 0.9364 0.6466 0.3712 1.0004 867
scale(abund.cov.1)sp3 0.2263 0.1504 0.0691 0.2292 0.5233 1.0314 544
scale(abund.cov.1)sp4 1.1878 0.2230 1.6583 1.1738 0.7742 1.0030 1374
scale(abund.cov.1)sp5 0.9623 0.1049 1.1701 0.9612 0.7581 1.0237 598
scale(abund.cov.1)sp6 1.5217 0.1423 1.2484 1.5162 1.8015 1.0051 313
Detection (logit scale):
Mean SD 2.5% 50% 97.5% Rhat ESS
(Intercept)sp1 0.3321 0.1829 0.0310 0.3364 0.6792 1.0004 1183
(Intercept)sp2 0.3788 0.2012 0.0137 0.3778 0.7676 1.0067 693
(Intercept)sp3 0.1727 0.1684 0.1728 0.1766 0.4935 1.0013 1101
(Intercept)sp4 0.5234 0.2651 0.0289 0.5128 1.0854 1.0037 1615
(Intercept)sp5 0.4695 0.1562 0.1618 0.4693 0.7693 1.0112 750
(Intercept)sp6 0.4967 0.1800 0.1488 0.4963 0.8401 1.0054 944
scale(det.cov.1)sp1 0.4734 0.1362 0.7518 0.4739 0.2130 1.0000 2471
scale(det.cov.1)sp2 0.0926 0.1415 0.3590 0.0951 0.1851 1.0067 1699
scale(det.cov.1)sp3 0.4350 0.1302 0.1888 0.4320 0.6992 1.0007 1770
scale(det.cov.1)sp4 1.0416 0.2600 0.5745 1.0309 1.5820 1.0008 2527
scale(det.cov.1)sp5 1.6153 0.1495 1.9110 1.6167 1.3151 1.0031 1271
scale(det.cov.1)sp6 1.2655 0.1803 0.9228 1.2595 1.6293 1.0091 1319
scale(det.cov.2)sp1 0.9304 0.1705 0.6081 0.9285 1.2792 1.0039 1970
scale(det.cov.2)sp2 0.9190 0.1750 0.5854 0.9133 1.2746 1.0082 1415
scale(det.cov.2)sp3 1.6974 0.1808 1.3632 1.6932 2.0722 1.0059 1268
scale(det.cov.2)sp4 0.8048 0.3666 1.5727 0.7855 0.1357 0.9999 2341
scale(det.cov.2)sp5 1.0518 0.1196 0.8219 1.0503 1.2946 1.0022 1325
scale(det.cov.2)sp6 2.3880 0.2568 1.8873 2.3819 2.8921 1.0061 1083
Note that the factor loadings themselves are not shown in the
summary()
output, but are available in the
lambda.samples
portion of the model list object.
# Rhats for lambda (the factor loadings)
out.lf.ms$rhat$lambda.lower.tri
[1] 1.017919 1.002436 1.003892 1.006558 1.011546 1.001850 1.010589 1.016527
[9] 1.022762 1.004086 1.007158 1.002827
# ESS for lambda
out.lf.ms$ESS$lambda
sp11 sp21 sp31 sp41 sp51 sp61 sp12 sp22
0.0000 261.2692 280.1641 1233.4555 327.5301 352.5995 0.0000 0.0000
sp32 sp42 sp52 sp62 sp13 sp23 sp33 sp43
236.3318 1070.1062 600.0428 272.8063 0.0000 0.0000 0.0000 901.9948
sp53 sp63
369.5350 520.2929
# Posterior quantiles for the latent factor loadings
summary(out.lf.ms$lambda.samples)$quantiles
2.5% 25% 50% 75% 97.5%
sp11 1.00000000 1.00000000 1.000000000 1.00000000 1.0000000
sp21 1.45020705 1.20217628 1.078588400 0.96892345 0.7621898
sp31 0.27796007 0.02699291 0.165861381 0.31021834 0.5751812
sp41 1.09700719 0.80531318 0.663312970 0.52961592 0.2660978
sp51 0.27722303 0.10164874 0.014920187 0.07385355 0.2669811
sp61 0.37972903 0.15107336 0.027611659 0.09431470 0.3500602
sp12 0.00000000 0.00000000 0.000000000 0.00000000 0.0000000
sp22 1.00000000 1.00000000 1.000000000 1.00000000 1.0000000
sp32 1.64391555 1.35223224 1.218405507 1.09120533 0.8597993
sp42 0.46494581 0.14087466 0.008538031 0.16599481 0.4846738
sp52 1.17796425 1.02918400 0.951172842 0.88042895 0.7406470
sp62 0.37193630 0.12395744 0.003766319 0.12976906 0.3571754
sp13 0.00000000 0.00000000 0.000000000 0.00000000 0.0000000
sp23 0.00000000 0.00000000 0.000000000 0.00000000 0.0000000
sp33 1.00000000 1.00000000 1.000000000 1.00000000 1.0000000
sp43 0.55155951 0.20084753 0.002504226 0.18554121 0.5627040
sp53 0.03197952 0.20272816 0.287836738 0.36823027 0.5302082
sp63 0.73614401 0.89519090 0.994576255 1.09378040 1.3018179
Notice the Rhat values are only reported for the elements of the
factor loadings matrix that are estimated and not fixed at any specific
value, while the ESS values are reported for all elements of the factor
loadings matrix, and take value 0 for those parameters that are fixed
for identifiability purposes. We can inspect the latent factor loadings
as well as the latent factors (stored in
out.lf.ms$w.samples
) to provide information on any
groupings that arise from the species in the modeled community. See
Appendix S1 in Doser, Finley, and Banerjee
(2023) for a discussion on using factor models as a modelbased
ordination technique.
Posterior predictive checks
We again use the ppcAbund()
function to perform a
posterior predictive check of our model
ppc.out.lf.ms < ppcAbund(out.lf.ms, fit.stat = 'chisquared', group = 2)
Currently on species 1 out of 6
Currently on species 2 out of 6
Currently on species 3 out of 6
Currently on species 4 out of 6
Currently on species 5 out of 6
Currently on species 6 out of 6
# Summarize with a Bayesian pvalue
summary(ppc.out.lf.ms)
Call:
ppcAbund(object = out.lf.ms, fit.stat = "chisquared", group = 2)
Samples per Chain: 40000
Burnin: 20000
Thinning Rate: 20
Number of Chains: 3
Total Posterior Samples: 3000

Community Level

Bayesian pvalue: 0.4874

Species Level

sp1 Bayesian pvalue: 0.441
sp2 Bayesian pvalue: 0.5447
sp3 Bayesian pvalue: 0.4443
sp4 Bayesian pvalue: 0.4447
sp5 Bayesian pvalue: 0.581
sp6 Bayesian pvalue: 0.469
Fit statistic: chisquared
Model selection using WAIC
We can use waicAbund()
to calculate the WAIC for
comparison to other models. Below, we compare the regular multispecies
Nmixture model to the latent factor multispecies Nmixture model.
# With latent factors
waicAbund(out.lf.ms, by.sp = TRUE)
N.max not specified. Setting upper index of integration of N to 10 plus
the largest estimated abundance value at each site in object$N.samples
Currently on species 1 out of 6
Currently on species 2 out of 6
Currently on species 3 out of 6
Currently on species 4 out of 6
Currently on species 5 out of 6
Currently on species 6 out of 6
elpd pD WAIC
1 288.5915 48.64827 674.4796
2 279.1837 55.95906 670.2855
3 281.2069 51.24437 664.9025
4 105.9096 19.12716 250.0736
5 379.0836 53.68815 865.5435
6 215.6488 44.41150 520.1205
# Without latent factors
waicAbund(out.ms, by.sp = TRUE)
N.max not specified. Setting upper index of integration of N to 10 plus
the largest estimated abundance value at each site in object$N.samples
Currently on species 1 out of 6
Currently on species 2 out of 6
Currently on species 3 out of 6
Currently on species 4 out of 6
Currently on species 5 out of 6
Currently on species 6 out of 6
elpd pD WAIC
1 368.2380 23.62859 783.7331
2 466.3158 54.39663 1041.4248
3 564.2021 71.35770 1271.1197
4 118.1010 12.09634 260.3946
5 532.5067 37.83546 1140.6842
6 295.3850 27.93759 646.6451
Prediction
Prediction proceeds exactly as before with the multispecies
Nmixture model, so we can simply copy the code from above and rerun it
with our latent factor model. The only difference is that we also
specify the include.w = FALSE
argument, which allows us to
visualize the expected abundance relationships without the sitelevel
random effects, analogous to the include.sp
argument that
we saw with the singlespecies spatial model.
out.pred < predict(out.lf.ms, X.0, ignore.RE = TRUE, include.w = FALSE)
# Look at the resulting object
str(out.pred)
List of 3
$ mu.0.samples: num [1:3000, 1:6, 1:100] 0.706 0.373 0.499 0.423 0.404 ...
$ N.0.samples : int [1:3000, 1:6, 1:100] 0 0 1 0 0 0 3 1 0 1 ...
$ call : language predict.lfMsNMix(object = out.lf.ms, X.0 = X.0, ignore.RE = TRUE, include.w = FALSE)
 attr(*, "class")= chr "predict.lfMsNMix"
mu.0.quants < apply(out.pred$mu.0.samples, c(2, 3), quantile,
prob = c(0.025, 0.5, 0.975))
# Put stuff together for plotting with ggplot2
mu.0.plot.df < data.frame(mu.med = c(mu.0.quants[2, , ]),
mu.low = c(mu.0.quants[1, , ]),
mu.high = c(mu.0.quants[3, , ]),
sp = rep(1:n.sp, nrow(X.0)),
abund.cov.1 = rep(cov.pred.vals, each = n.sp))
ggplot(mu.0.plot.df, aes(x = abund.cov.1, y = mu.med)) +
geom_ribbon(aes(ymin = mu.low, ymax = mu.high), fill = 'grey70') +
geom_line() +
theme_bw() +
facet_wrap(vars(sp), scales = 'free_y') +
labs(x = 'Covariate', y = 'Expected abundance')
Spatial factor multispecies Nmixture models
Basic model description
Our final, and most complex, Nmixture model that we fit in
spAbundance
is a multispecies spatiallyexplicit Nmixture
model. This model is nearly identical to the latent factor multispecies
Nmixture model, except we give a spatial structure to the latent
factors instead of assuming they are independent from each other. By
modeling the latent factors with a spatial structure, we will often see
such a model have improved predictive performance relative to a latent
factor multispecies Nmixture model (Doser,
Finley, and Banerjee 2023). This model again is an
abundancebased JSDM, but now it simultaneously accounts for imperfect
detection, residual correlations between species, and spatial
autocorrelation. This is a direct extension of the spatial factor
multispecies occupancy model presented in Doser,
Finley, and Banerjee (2023).
The model is identical to the previously described latent factor multispecies Nmixture model with the exception that the latent factors are assumed to have a spatial structure to them. More specifically, our model for speciesspecific abundance is again
\[\begin{equation} \text{log}(\mu_{i}(\boldsymbol{s}_j)) = \boldsymbol{x}(\boldsymbol{s}_j)^\top\boldsymbol{\beta}_i + \text{w}^\ast_{i}(\boldsymbol{s}_j). \end{equation}\]
Note again that for spatial models we use the notation \(\boldsymbol{s}_j\) to make it clear that the model is spatiallyexplicit and relies upon the coordinates (\(\boldsymbol{s}_j\)) at each site \(j\) to estimate this spatial pattern. As with the latent factor model, we decompose \(\text{w}^\ast_i(\boldsymbol{s}_j)\) into a linear combination of \(q\) latent variables (i.e., factors) and their associated speciesspecific coefficients (i.e., factor loadings) according to
\[\begin{equation} \text{w}^\ast_{i}(\boldsymbol{s}_j) = \boldsymbol{\lambda}_i^\top\textbf{w}(\boldsymbol{s}_j). \end{equation}\]
Now, instead of modeling \(\textbf{w}(\boldsymbol{s}_j)\) as independent normal latent variables, we assume \(\textbf{w}(\boldsymbol{s}_j)\) arise from spatial processes, allowing us to account for residual correlation in speciesspecific abundance. More specifically, each spatial factor \(\textbf{w}_r\) for each \(r = 1, \dots, q\) is modeled using a Nearest Neighbor Gaussian Process (Datta et al. 2016), i.e.,
\[\begin{equation} \textbf{w}_r \sim \text{Normal}(\boldsymbol{0}, \tilde{\boldsymbol{C}}_r(\boldsymbol{\theta}_r)), \end{equation}\]
where \(\tilde{\boldsymbol{C}}_r(\boldsymbol{\theta}_r)\) is the NNGPderived correlation matrix for the \(r^{\text{th}}\) spatial factor. Note that the spatial variance parameter for each spatial factor is assumed to be 1 for identifiability purposes. The vector \(\boldsymbol{\theta}_r\) consists of parameters governing the spatial process for each spatial factor according to some spatial correlation function, as we saw with the singlespecies Nmixture model. Thus, for the exponential, spherical, and Gaussian correlation functions, \(\boldsymbol{\theta}_r\) includes a spatial decay parameter \(\phi_r\), while the Matern correlation function includes an additional spatial smoothness parameter, \(\nu_r\).
We assume the same priors and identifiability constraints as the latent factor multispecies Nmixture model. We assign a uniform prior for the spatial decay parameter, \(\phi_r\), and the spatial smoothness parameters, \(\nu_r\), if using a Matern correlation function.
Notice that this spatial factor modeling approach is the only
approach we implement in spAbundance
for modeling
multispecies datasets while accounting for spatial autocorrelation.
While we could envision fitting a separate spatial random effect for
each species in our multispecies data set (as we implemented in
spOccupancy
in the spMsPGOcc
function), we
prefer (and recommend) using the spatial factor modeling approach as (1)
it is far more computationally efficient than fitting a separate spatial
effect for each species; (2) it explicitly acknowledges dependence
between species; (3) estimating a separate spatial effect for each
species would be very difficult to do with multiple rare species in the
data set; and (4) even if the species are independent (i.e., there are
no residual correlations), the spatial factor modeling approach performs
extremely similarly to a model with a separate spatial process for each
species (Doser, Finley, and Banerjee
2023), while still being substantially faster.
Fitting spatial factor multispecies Nmixture models with
sfMsNMix()
The function sfMsNMix()
fits spatial factor
multispecies Nmixture models in spAbundance
. The
arguments are very similar to lfMsNMix()
and
spNMix()
.
sfMsNMix(abund.formula, det.formula, data, inits, priors,
tuning, cov.model = 'exponential', NNGP = TRUE, n.neighbors = 15,
search.type = 'cb', n.factors, n.batch, batch.length, accept.rate = 0.43,
family = 'Poisson', n.omp.threads = 1, verbose = TRUE, n.report = 100,
n.burn = round(.10 * n.batch * batch.length), n.thin = 1,
n.chains = 1, ...)
We will again fit the model with three spatial factors, and will use the same covariates/random effects on abundance and detection as we have done throughout the vignette
n.factors < 3
abund.ms.formula < ~ scale(abund.cov.1) + (1  abund.factor.1)
det.ms.formula < ~ scale(det.cov.1) + scale(det.cov.2)
# Just a reminder of the data format for fitting multispecies Nmixture models.
str(dataNMixSim)
List of 4
$ y : int [1:6, 1:225, 1:3] 1 0 1 0 0 0 NA NA NA NA ...
$ abund.covs: num [1:225, 1:2] 0.373 0.706 0.202 1.588 0.138 ...
.. attr(*, "dimnames")=List of 2
.. ..$ : NULL
.. ..$ : chr [1:2] "abund.cov.1" "abund.factor.1"
$ det.covs :List of 2
..$ det.cov.1: num [1:225, 1:3] 1.28 NA NA NA 1.04 ...
..$ det.cov.2: num [1:225, 1:3] 2.03 NA NA NA 0.796 ...
$ coords : num [1:225, 1:2] 0 0.0714 0.1429 0.2143 0.2857 ...
.. attr(*, "dimnames")=List of 2
.. ..$ : NULL
.. ..$ : chr [1:2] "X" "Y"
Initial values are identical to what we saw with
lfMsNMix()
, but we will also specify the initial values for
the spatial decay parameters for each spatial factor (we use an
exponential correlation model so we do not need to specify any initial
values for the nu
, the spatial smoothness parameter used
when cov.model = 'matern'
).
# Number of species
n.sp < dim(dataNMixSim$y)[1]
# Initiate all lambda initial values to 0.
lambda.inits < matrix(0, n.sp, n.factors)
# Set diagonal elements to 1
diag(lambda.inits) < 1
# Set lower triangular elements to random values from a standard normal distribution
lambda.inits[lower.tri(lambda.inits)] < rnorm(sum(lower.tri(lambda.inits)))
# Check it out.
lambda.inits
[,1] [,2] [,3]
[1,] 1.00000000 0.00000000 0.0000000
[2,] 1.71182094 1.00000000 0.0000000
[3,] 1.02105136 1.30685972 1.0000000
[4,] 2.07933362 1.77353032 1.1023205
[5,] 0.03136039 0.04679833 0.6785791
[6,] 1.05228942 1.62190825 0.8646384
w.inits < matrix(0, n.factors, ncol(dataNMixSim$y))
# Pairwise distances between all sites
dist.mat < dist(dataNMixSim$coords)
# Exponential covariance model
cov.model < 'exponential'
ms.inits < list(alpha.comm = 0,
beta.comm = 0,
beta = 0,
alpha = 0,
tau.sq.beta = 1,
lambda = lambda.inits,
kappa = 1,
phi = 3 / mean(dist.mat),
w = w.inits,
tau.sq.alpha = 1,
sigma.sq.mu = 0.5,
N = apply(dataNMixSim$y, c(1, 2), max, na.rm = TRUE))
Priors are the same as with the latent factor multispecies Nmixture
model, where we also add in our default prior for the spatial decay
parameters, which allows the effective spatial range of each spatial
factor to range from the maximum intersite distance to the minimum
intersite distance. Notice the prior for phi
is now
specified as a list as opposed to an atomic vector as we did for the
singlespecies case. If desired, the list format allows you to easily
specify a different prior for each of the spatial factors. See
?sfMsNMix
for more details.
max.dist < max(dist.mat)
min.dist < min(dist.mat)
ms.priors < list(beta.comm.normal = list(mean = 0, var = 2.72),
alpha.comm.normal = list(mean = 0, var = 2.72),
tau.sq.beta.ig = list(a = 0.1, b = 0.1),
tau.sq.alpha.ig = list(a = 0.1, b = 0.1),
sigma.sq.mu.ig = list(a = 0.1, b = 0.1),
kappa.unif = list(a = 0, b = 100),
phi.unif = list(3 / max.dist, 3 / min.dist))
We lastly specify the tuning values and run the model for 40,000 iterations with a burnin of 20,000 samples and a thinning rate of 20, for each of 3 chains, yielding a total of 3000 posterior samples. As before, we will use a Poisson distribution and 15 nearest neighbors.
# Specify initial tuning values
ms.tuning < list(beta = 0.3, alpha = 0.3, beta.star = 0.5, kappa = 0.5,
w = 0.5, lambda = 0.5, phi = 0.5)
# Approx. run time: ~13 min
out.sf.ms < sfMsNMix(abund.formula = abund.ms.formula,
det.formula = det.ms.formula,
data = dataNMixSim,
inits = ms.inits,
n.batch = 1600,
n.factors = n.factors,
tuning = ms.tuning,
batch.length = 25,
priors = ms.priors,
n.neighbors = 15,
n.omp.threads = 1,
verbose = TRUE,
n.report = 400,
n.burn = 20000,
n.thin = 20,
n.chains = 3)

Preparing to run the model


Building the neighbor list


Building the neighbors of neighbors list


Model description

Spatial Factor NNGP Multispecies Poisson NMixture model fit with 225 sites and 6 species.
Samples per Chain: 40000
Burnin: 20000
Thinning Rate: 20
Number of Chains: 3
Total Posterior Samples: 3000
Using the exponential spatial correlation model.
Using 3 latent spatial factors.
Using 15 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: 400 of 1600, 25.00%
Species Parameter Acceptance Tuning
1 beta[1] 24.0 0.15978
1 alpha[1] 44.0 0.24809
2 beta[1] 64.0 0.14171
2 alpha[1] 36.0 0.22448
3 beta[1] 48.0 0.12076
3 alpha[1] 28.0 0.22003
4 beta[1] 56.0 0.31538
4 alpha[1] 28.0 0.42572
5 beta[1] 32.0 0.10086
5 alpha[1] 40.0 0.19129
6 beta[1] 48.0 0.12569
6 alpha[1] 40.0 0.23364
1 phi 32.0 0.49502
2 phi 32.0 0.66821
3 phi 32.0 0.47561

Batch: 800 of 1600, 50.00%
Species Parameter Acceptance Tuning
1 beta[1] 44.0 0.16630
1 alpha[1] 32.0 0.24809
2 beta[1] 44.0 0.13890
2 alpha[1] 52.0 0.21141
3 beta[1] 48.0 0.10927
3 alpha[1] 68.0 0.20312
4 beta[1] 36.0 0.30914
4 alpha[1] 48.0 0.45205
5 beta[1] 32.0 0.09691
5 alpha[1] 36.0 0.18379
6 beta[1] 44.0 0.12076
6 alpha[1] 36.0 0.23364
1 phi 48.0 0.53625
2 phi 28.0 0.68171
3 phi 52.0 0.46620

Batch: 1200 of 1600, 75.00%
Species Parameter Acceptance Tuning
1 beta[1] 16.0 0.16966
1 alpha[1] 40.0 0.25310
2 beta[1] 36.0 0.11837
2 alpha[1] 40.0 0.22003
3 beta[1] 28.0 0.11602
3 alpha[1] 40.0 0.22448
4 beta[1] 36.0 0.33488
4 alpha[1] 48.0 0.46118
5 beta[1] 40.0 0.09887
5 alpha[1] 36.0 0.17658
6 beta[1] 64.0 0.12569
6 alpha[1] 64.0 0.24318
1 phi 48.0 0.59265
2 phi 32.0 0.68171
3 phi 44.0 0.49502

Batch: 1600 of 1600, 100.00%

Chain 2

Sampling ...
Batch: 400 of 1600, 25.00%
Species Parameter Acceptance Tuning
1 beta[1] 52.0 0.16966
1 alpha[1] 52.0 0.26875
2 beta[1] 52.0 0.13081
2 alpha[1] 40.0 0.22448
3 beta[1] 56.0 0.11837
3 alpha[1] 52.0 0.20722
4 beta[1] 28.0 0.32825
4 alpha[1] 20.0 0.46118
5 beta[1] 52.0 0.09691
5 alpha[1] 72.0 0.19129
6 beta[1] 20.0 0.12320
6 alpha[1] 16.0 0.21568
1 phi 52.0 0.55814
2 phi 48.0 0.69548
3 phi 64.0 0.45697

Batch: 800 of 1600, 50.00%
Species Parameter Acceptance Tuning
1 beta[1] 28.0 0.16966
1 alpha[1] 40.0 0.25310
2 beta[1] 48.0 0.13615
2 alpha[1] 32.0 0.20722
3 beta[1] 44.0 0.12320
3 alpha[1] 60.0 0.21141
4 beta[1] 24.0 0.33488
4 alpha[1] 56.0 0.45205
5 beta[1] 36.0 0.10498
5 alpha[1] 44.0 0.16966
6 beta[1] 52.0 0.13346
6 alpha[1] 24.0 0.23836
1 phi 44.0 0.52564
2 phi 48.0 0.65498
3 phi 48.0 0.47561

Batch: 1200 of 1600, 75.00%
Species Parameter Acceptance Tuning
1 beta[1] 48.0 0.15351
1 alpha[1] 32.0 0.24809
2 beta[1] 40.0 0.13346
2 alpha[1] 56.0 0.22901
3 beta[1] 48.0 0.11837
3 alpha[1] 28.0 0.20722
4 beta[1] 48.0 0.30914
4 alpha[1] 48.0 0.45205
5 beta[1] 36.0 0.10927
5 alpha[1] 32.0 0.17658
6 beta[1] 24.0 0.13346
6 alpha[1] 36.0 0.21141
1 phi 28.0 0.59265
2 phi 36.0 0.69548
3 phi 60.0 0.52564

Batch: 1600 of 1600, 100.00%

Chain 3

Sampling ...
Batch: 400 of 1600, 25.00%
Species Parameter Acceptance Tuning
1 beta[1] 56.0 0.15978
1 alpha[1] 40.0 0.26343
2 beta[1] 40.0 0.13615
2 alpha[1] 52.0 0.21141
3 beta[1] 20.0 0.11602
3 alpha[1] 32.0 0.20722
4 beta[1] 32.0 0.30302
4 alpha[1] 36.0 0.47049
5 beta[1] 32.0 0.09887
5 alpha[1] 48.0 0.17308
6 beta[1] 72.0 0.12569
6 alpha[1] 40.0 0.23836
1 phi 36.0 0.69548
2 phi 52.0 0.55814
3 phi 28.0 1.03754

Batch: 800 of 1600, 50.00%
Species Parameter Acceptance Tuning
1 beta[1] 40.0 0.15978
1 alpha[1] 52.0 0.27418
2 beta[1] 36.0 0.14171
2 alpha[1] 48.0 0.22901
3 beta[1] 48.0 0.12320
3 alpha[1] 44.0 0.23836
4 beta[1] 32.0 0.29113
4 alpha[1] 40.0 0.46118
5 beta[1] 32.0 0.10086
5 alpha[1] 36.0 0.18750
6 beta[1] 48.0 0.13890
6 alpha[1] 40.0 0.22003
1 phi 52.0 0.59265
2 phi 36.0 0.70953
3 phi 44.0 1.01700

Batch: 1200 of 1600, 75.00%
Species Parameter Acceptance Tuning
1 beta[1] 36.0 0.17658
1 alpha[1] 36.0 0.24809
2 beta[1] 28.0 0.13081
2 alpha[1] 52.0 0.21141
3 beta[1] 60.0 0.11147
3 alpha[1] 48.0 0.22448
4 beta[1] 40.0 0.32825
4 alpha[1] 48.0 0.48969
5 beta[1] 32.0 0.09887
5 alpha[1] 44.0 0.18015
6 beta[1] 48.0 0.13615
6 alpha[1] 40.0 0.22901
1 phi 56.0 0.62930
2 phi 56.0 0.68171
3 phi 64.0 0.93881

Batch: 1600 of 1600, 100.00%
summary(out.sf.ms)
Call:
sfMsNMix(abund.formula = abund.ms.formula, det.formula = det.ms.formula,
data = dataNMixSim, inits = ms.inits, priors = ms.priors,
tuning = ms.tuning, n.neighbors = 15, n.factors = n.factors,
n.batch = 1600, batch.length = 25, n.omp.threads = 1, verbose = TRUE,
n.report = 400, n.burn = 20000, n.thin = 20, n.chains = 3)
Samples per Chain: 40000
Burnin: 20000
Thinning Rate: 20
Number of Chains: 3
Total Posterior Samples: 3000
Run Time (min): 13.2365

Community Level

Abundance Means (log scale):
Mean SD 2.5% 50% 97.5% Rhat ESS
(Intercept) 0.9343 0.4255 1.7281 0.9423 0.0466 1.0051 838
scale(abund.cov.1) 0.1342 0.4585 1.0300 0.1450 0.8202 1.0022 3000
Abundance Variances (log scale):
Mean SD 2.5% 50% 97.5% Rhat ESS
(Intercept) 1.0600 2.8494 0.1754 0.6861 3.9406 1.3365 1819
scale(abund.cov.1) 1.4854 1.4044 0.3510 1.0904 5.1941 1.0010 3000
Abundance Random Effect Variances (log scale):
Mean SD 2.5% 50% 97.5% Rhat ESS
(Intercept)abund.factor.1 0.4757 0.1302 0.2729 0.4593 0.7851 1.0386 950
Detection Means (logit scale):
Mean SD 2.5% 50% 97.5% Rhat ESS
(Intercept) 0.3864 0.1820 0.0318 0.3857 0.7504 1.0115 1822
scale(det.cov.1) 0.0892 0.5032 0.9121 0.0922 1.0757 1.0011 3000
scale(det.cov.2) 0.9268 0.5178 0.1819 0.9603 1.8848 1.0035 3608
Detection Variances (logit scale):
Mean SD 2.5% 50% 97.5% Rhat ESS
(Intercept) 0.1446 0.1868 0.0258 0.0950 0.5407 1.0068 2675
scale(det.cov.1) 1.7729 1.6139 0.4210 1.2867 6.0963 1.0231 2697
scale(det.cov.2) 1.7281 1.8076 0.3667 1.2113 6.4946 1.0038 2823

Species Level

Abundance (log scale):
Mean SD 2.5% 50% 97.5% Rhat ESS
(Intercept)sp1 0.6241 0.2730 1.1502 0.6256 0.0598 1.0071 433
(Intercept)sp2 0.7179 0.3571 1.3657 0.7298 0.0583 1.0431 113
(Intercept)sp3 0.9324 0.4650 1.9352 0.9168 0.0717 1.0722 57
(Intercept)sp4 2.1753 0.3847 2.9390 2.1663 1.4555 1.0075 799
(Intercept)sp5 0.1160 0.3368 0.8362 0.1020 0.5234 1.0527 105
(Intercept)sp6 1.3345 0.3976 2.1784 1.3099 0.6095 1.0125 124
scale(abund.cov.1)sp1 0.0875 0.1075 0.1250 0.0878 0.3053 1.0183 1413
scale(abund.cov.1)sp2 0.6506 0.1350 0.9079 0.6502 0.3923 1.0007 980
scale(abund.cov.1)sp3 0.2507 0.1197 0.0204 0.2501 0.4865 1.0246 920
scale(abund.cov.1)sp4 1.2147 0.2193 1.6649 1.2065 0.8095 1.0059 1396
scale(abund.cov.1)sp5 0.9317 0.0926 1.1191 0.9305 0.7539 1.0083 965
scale(abund.cov.1)sp6 1.4397 0.1222 1.2001 1.4381 1.6868 1.0543 459
Detection (logit scale):
Mean SD 2.5% 50% 97.5% Rhat ESS
(Intercept)sp1 0.3219 0.1845 0.0775 0.3225 0.6723 1.0103 1049
(Intercept)sp2 0.3575 0.2066 0.0765 0.3647 0.7545 1.0028 760
(Intercept)sp3 0.1776 0.1681 0.1671 0.1823 0.5054 1.0189 805
(Intercept)sp4 0.5308 0.2705 0.0301 0.5201 1.0877 1.0260 1679
(Intercept)sp5 0.4642 0.1564 0.1446 0.4676 0.7662 1.0236 891
(Intercept)sp6 0.4952 0.1789 0.1478 0.4942 0.8619 1.0273 985
scale(det.cov.1)sp1 0.4703 0.1385 0.7379 0.4707 0.2039 1.0013 2320
scale(det.cov.1)sp2 0.0889 0.1390 0.3608 0.0913 0.1706 1.0171 1428
scale(det.cov.1)sp3 0.4532 0.1319 0.1973 0.4502 0.7131 1.0029 1780
scale(det.cov.1)sp4 1.0334 0.2660 0.5437 1.0201 1.5782 1.0056 2637
scale(det.cov.1)sp5 1.6132 0.1482 1.9107 1.6040 1.3382 1.0120 1486
scale(det.cov.1)sp6 1.2433 0.1745 0.9113 1.2396 1.6013 1.0048 1485
scale(det.cov.2)sp1 0.9434 0.1669 0.6170 0.9423 1.2647 1.0119 2127
scale(det.cov.2)sp2 0.9025 0.1780 0.5557 0.8980 1.2580 1.0176 1362
scale(det.cov.2)sp3 1.6538 0.1808 1.3076 1.6526 2.0223 1.0114 1472
scale(det.cov.2)sp4 0.7864 0.3481 1.5126 0.7709 0.1586 1.0019 2330
scale(det.cov.2)sp5 1.0453 0.1219 0.8089 1.0455 1.2848 1.0122 1632
scale(det.cov.2)sp6 2.2877 0.2514 1.8327 2.2807 2.7917 1.0061 1312

Spatial Covariance

Mean SD 2.5% 50% 97.5% Rhat ESS
phi1 14.5129 4.7160 7.1550 13.8661 25.3105 1.6126 131
phi2 12.3673 7.5036 2.4742 12.2804 28.5995 4.7612 13
phi3 11.5633 9.0965 2.9467 7.5924 36.4446 10.3050 14
Posterior predictive checks
As before, we can use the ppcAbund()
function to perform
a posterior predictive check of our model
ppc.out.sf.ms < ppcAbund(out.sf.ms, fit.stat = 'freemantukey', group = 1)
Currently on species 1 out of 6
Currently on species 2 out of 6
Currently on species 3 out of 6
Currently on species 4 out of 6
Currently on species 5 out of 6
Currently on species 6 out of 6
# Summarize with a Bayesian pvalue
summary(ppc.out.sf.ms)
Call:
ppcAbund(object = out.sf.ms, fit.stat = "freemantukey", group = 1)
Samples per Chain: 40000
Burnin: 20000
Thinning Rate: 20
Number of Chains: 3
Total Posterior Samples: 3000

Community Level

Bayesian pvalue: 0.6259

Species Level

sp1 Bayesian pvalue: 0.6587
sp2 Bayesian pvalue: 0.696
sp3 Bayesian pvalue: 0.7397
sp4 Bayesian pvalue: 0.5527
sp5 Bayesian pvalue: 0.759
sp6 Bayesian pvalue: 0.3497
Fit statistic: freemantukey
Model selection using WAIC
We use waicAbund()
to calculate the WAIC for all species
in the data set, and compare the WAIC to that obtained using the
nonspatial latent factor multispecies Nmixture model.
# Nonspatial latent factor model
waicAbund(out.lf.ms, by.sp = FALSE)
N.max not specified. Setting upper index of integration of N to 10 plus
the largest estimated abundance value at each site in object$N.samples
Currently on species 1 out of 6
Currently on species 2 out of 6
Currently on species 3 out of 6
Currently on species 4 out of 6
Currently on species 5 out of 6
Currently on species 6 out of 6
elpd pD WAIC
1549.6241 273.0785 3645.4052
# Spatial factor model
waicAbund(out.sf.ms, by.sp = FALSE)
N.max not specified. Setting upper index of integration of N to 10 plus
the largest estimated abundance value at each site in object$N.samples
Currently on species 1 out of 6
Currently on species 2 out of 6
Currently on species 3 out of 6
Currently on species 4 out of 6
Currently on species 5 out of 6
Currently on species 6 out of 6
elpd pD WAIC
1553.1935 254.9933 3616.3734
Prediction
We can again use the predict()
function to predict
abundance at a set of new locations to generate a species distribution
map or to generate a marginal effects plot. Here, we again predict
abundance for the purpose of generating a marginal effects plot, and so
we set the argument include.sp = FALSE
to ignore the
spatial random effects. See the distance sampling vignette for an
example of how to generate a map of abundance across some region of
interest.
# Recall what the prediction design matrix looks like
str(X.0)
num [1:100, 1:2] 1 1 1 1 1 1 1 1 1 1 ...
 attr(*, "dimnames")=List of 2
..$ : NULL
..$ : chr [1:2] "intercept" "abund.cov.1"
List of 4
$ mu.0.samples: num [1:3000, 1:6, 1:100] 0.358 0.242 0.338 0.352 0.42 ...
$ N.0.samples : int [1:3000, 1:6, 1:100] 0 0 0 0 0 0 0 1 0 1 ...
$ call : language predict.sfMsNMix(object = out.sf.ms, X.0 = X.0, ignore.RE = TRUE, include.sp = FALSE)
$ object.class: chr "sfMsNMix"
 attr(*, "class")= chr "predict.sfMsNMix"
mu.0.quants < apply(out.pred$mu.0.samples, c(2, 3), quantile,
prob = c(0.025, 0.5, 0.975))
# Put stuff together for plotting with ggplot2
mu.0.plot.df < data.frame(mu.med = c(mu.0.quants[2, , ]),
mu.low = c(mu.0.quants[1, , ]),
mu.high = c(mu.0.quants[3, , ]),
sp = rep(1:n.sp, nrow(X.0)),
abund.cov.1 = rep(cov.pred.vals, each = n.sp))
ggplot(mu.0.plot.df, aes(x = abund.cov.1, y = mu.med)) +
geom_ribbon(aes(ymin = mu.low, ymax = mu.high), fill = 'grey70') +
geom_line() +
theme_bw() +
facet_wrap(vars(sp), scales = 'free_y') +
labs(x = 'Covariate', y = 'Expected abundance')