# Measurement of $\omega$ meson production in pp and
p–Pb collisions at $\sqrt{s_{\scalebox{0.5}{NN}}}=5.02\leavevmode\nobreak\
$TeV with ALICE

Nicolas Strangmann

###### Abstract

The ALICE experiment at the LHC investigates the properties of the hot and
dense nuclear matter created in heavy-ion collisions. By comparing the
particle production in pp and p–Pb collisions, possible nuclear initial state
effects can be isolated. Measurements of the $\omega$ meson
$p_{\text{T}}$-spectra in pp and p–Pb collisions not only allow for a
determination of the nuclear modification factor $R_{\text{pPb}}$, but also
provide insight into the fragmentation process and serve as vital input for
decay background simulations for direct photons.

In this contribution, measurements of the $\omega$ meson production in pp and
p–Pb collisions at $\sqrt{s_{\scalebox{0.5}{NN}}}=5.02\leavevmode\nobreak\
$TeV are presented. This includes the signal extraction and various
corrections of the $\omega$ meson yields, leading to their production cross
sections and the first measured nuclear modification factor $R_{\text{pPb}}$
of the $\omega$ meson at LHC energies.

##  1 Motivation

Measurements of particle production at the LHC can help to constrain parton
distribution functions as well as fragmentation functions. While the
production of pseudoscalar mesons ($\pi$, $\eta$, $K$) has been extensively
studied in recent years [1, 2, 3, 4, 5], the production of vector mesons like
the $\omega$ has been investigated less [6], hampering theoretical models to
better describe the production of vector mesons [7]. Moreover, measuring the
production of $\omega$ mesons in pp and p–Pb collisions, one can probe Cold
Nuclear Matter (CNM) effects. Modifications due to this cold medium can be
quantified by the nuclear modification factor $R_{\text{pPb}}$, the ratio of
the production cross section in p–Pb collisions to that in pp collisions
scaled by the number of nucleons within the Pb-nucleus:

$\displaystyle
R_{\text{pPb}}=\frac{1}{A_{\text{pPb}}}\frac{\text{d}^{2}\sigma_{\text{pPb}}/\text{d}p_{\text{T}}\text{d}y}{\text{d}^{2}\sigma_{\text{pp}}/\text{d}p_{\text{T}}\text{d}y}.$
(1)

Finally, a precise measurement of the $\omega$ meson production serves as
vital input for direct photon measurements, as the $\omega$ meson is the third
largest contributor to the background of the photon measurement [8], and
therefore this measurement of $\omega$ mesons reduces the systematic
uncertainties of direct photon analyses.

##  2 Experimental setup and data sets

ALICE as a multipurpose experiment is very well suited for the measurement of
$\omega$ mesons. The $\omega$ mesons decay very close to the collision vertex
($c\tau=23\leavevmode\nobreak\ \text{f}m$) [9] and can be reconstructed in
ALICE via the dominant decay channel into three pions
($\omega\rightarrow\pi^{+}\pi^{-}\pi^{0}$). Charged pions can be tracked in
the inner tracking system (ITS) and the time projection chamber (TPC). The
neutral pions decay almost exclusively (99 %) into two photons [9]. These are
either measured with the electromagnetic calorimeter (EMCal) or via the
reconstruction of $e^{+}e^{-}$ pairs in the TPC originating from their
conversion in the detector material with the so-called photon conversion
method (PCM). The data presented here was collected with a minimum bias
trigger requiring a coincident signal in both VZERO detectors on either side
of the experiment. In November 2016, ALICE recorded about 600 million p–Pb
collisions, and in November 2017 about 900 million pp collisions at
$\sqrt{s_{\scalebox{0.5}{NN}}}=5.02\leavevmode\nobreak\ $TeV. The integrated
luminosity inspected in this analysis in pp (p–Pb) collisions amounts to
$\mathcal{L}=18.5\leavevmode\nobreak\ $n$b^{-1}$ ($0.3\leavevmode\nobreak\
$n$b^{-1}$).

##  3 Analysis

###  3.1 Pion reconstruction

The reconstruction of the $\omega$ meson relies on a pure sample of pions to
increase the significance of the $\omega$ signal and therefore the
$p_{\text{T}}$-reach of the measurement. Charged pions are selected from all
reconstructed tracks based on their specific energy loss in the TPC
(n$\sigma<3$). The decay photons measured in the EMCal and those reconstructed
using PCM are combined to reconstruct neutral pions. This results in the three
neutral pion reconstruction methods EMCal, PCM, and PCM-EMCal, in case one of
the photons converts. For the subsequent $\omega$ meson reconstruction only
photon pairs with a reconstructed invariant mass close to the neutral pion
mass are used.

###  3.2 $\omega$ meson reconstruction

Combining all $\pi^{+}$, $\pi^{-}$ and $\pi^{0}$ mesons within each event,
$\omega$ meson candidates are reconstructed. Figure 1 shows the resulting
number of $\omega$ meson candidates as a function of their reconstructed
invariant mass in selected $p_{\text{T}}$ intervals in p–Pb collisions for the
three neutral pion reconstruction methods.

(a) PCM

(b) PCM-EMCal

(c) EMCal

Figure 1: Distributions of the invariant mass of reconstructed $\omega$ meson
candidates in p–Pb collisions using the three different $\pi^{0}$
reconstruction methods.

The combinatorial background can be described with a third order polynomial
function. The $\omega$ meson signal, as shown in red in figure 1, is obtained
by subtracting the estimated background. The mass and width of the
reconstructed $\omega$ mesons are extracted by parameterizing the signal with
a Gaussian. The values obtained are shown in figure 2 for data (full markers)
as well as for simulations (open markers). PYTHIA Monash 8.2 and DPMJET
combined with GEANT3 were used to simulate the detector response. While the
reconstructed mass is in agreement with literature [9], the width is broadened
by the detector resolutions. The extracted properties of the $\omega$ meson
are well described by the Monte Carlo simulations. The raw $\omega$ meson
yield is extracted by integration of the signal around the reconstructed
$\omega$ mass.

(a) Mass and width in pp collisions

(b) Mass and width in p–Pb collisions

Figure 2: Mass and width of the $\omega$ meson signal extracted using the
different $\pi^{0}$ reconstruction methods. Full markers represent data, open
markers the simulation, and the dotted lines show published values [9].

###  3.3 Corrections and combination of methods

The Lorentz invariant production cross section is calculated from the
extracted raw $\omega$ meson yield $N^{\omega}$ by applying corrections, as
summarized in the following equation:

$\displaystyle
E\frac{\text{d}^{3}\sigma_{\omega}}{\text{d}^{3}p}=\frac{1}{\mathcal{L}_{\text{int}}}\frac{1}{2\pi
p_{\text{T}}}\frac{1}{A\epsilon\mathcal{B}}\frac{N^{\omega}}{\Delta
p_{\text{T}}\Delta y}.$ (2)

Besides normalization by the Luminosity $\mathcal{L}_{\text{int}}$ and
$2\pi\Delta y$, this includes the branching ratio
$\mathcal{B}(\omega\rightarrow\pi^{+}\pi^{-}\gamma\gamma)\approx 89\%$ as well
as the detector acceptance $A$ and reconstruction efficiency $\epsilon$. Using
equations 2 and 1, the cross sections in pp and p–Pb as well as the
$R_{\text{pPb}}$ were calculated individually for each $\pi^{0}$
reconstruction method. To account for the difference in the investigated
rapidity regions in the two collision systems (see figure 4(a)), the pp
spectrum used for the $R_{\text{pPb}}$ calculation was corrected by a factor
of about $0.99$ based on the ratio of $\omega$ spectra generated with PYTHIA
Monash 8.2 for both kinematic regions. For the three reconstruction methods of
the $\pi^{0}$ and thus $\omega$, the systematic uncertainties are estimated
from variations of the applied pion selection criteria and extraction
parameters. Figure 3 shows the ratio of the extracted invariant cross sections
to a Tsallis parametrization in pp and p–Pb. A good agreement is observed
between these methods within the statistical (bars) and systematic (boxes)
uncertainties. The cross sections and nuclear modification factors extracted
using the different methods are combined based on their uncertainties and
correlations, using the Best Linear Unbiased Estimate (BLUE) method [10].
While PCM allows for an extraction of the $\omega$ down to low $p_{\text{T}}$
thanks to its good energy resolution, the photon conversion probability of
around 8% restricts the range at higher transverse momenta, where the EMCal
measurement dominates the measurement.

(a) Ratios in pp collisions

(b) Ratios in p–Pb collisions

Figure 3: Ratios of the invariant cross sections extracted using the three
different $\pi^{0}$ reconstruction methods to a combined Tsallis
parametrization.

##  4 Results

#### Production cross section

The extracted production cross sections of $\omega$ mesons for pp and p–Pb
collisions at $\sqrt{s_{\scalebox{0.5}{NN}}}=5.02\leavevmode\nobreak\ $TeV are
shown in figure 4(a). The pp cross section covers the transverse momentum
range $1.8\leq p_{\text{T}}<16\,\text{GeV}/c$ within $|y|<0.85$. In p–Pb
collisions, the cross section is extracted for transverse momenta $2.2\leq
p_{\text{T}}<20\,\text{GeV}/c$ for $-1.315<y<0.385$. The lower panel in figure
4(a) shows the ratios of the data as well as various MC predictions to the
respective Tsallis parametrization. While EPOS LHC describes the spectrum in
p–Pb well, it overpredicts the production in pp by up to 100 %. PYTHIA Monash
8.2 and DPMJET describe the shape of the spectra in pp and p–Pb collisions,
respectively, but over-/underestimate them by 30-40 %.

####  $\omega/\pi^{0}$ ratio

The ratio of $\omega$ mesons to neutral pions as function of $p_{\text{T}}$ at
$\sqrt{s_{\scalebox{0.5}{NN}}}=5.02\leavevmode\nobreak\ $TeV is shown in
figure 4(b) for pp (orange) and p–Pb (violet) collisions. The result of a
measurement of $\omega$ production at $\sqrt{s}=13\leavevmode\nobreak\ $TeV is
shown in green. The agreement between the three measurements suggest an
independence of the $\omega/\pi^{0}$ ratio of the collision system and energy
within the uncertainties.

#### Nuclear modification factor

Figure 4(c) shows the first measurement of the nuclear modification factor
$R_{\text{pPb}}$ of the $\omega$ meson at LHC energies. As $R_{\text{pPb}}$ is
consistent with unity, no nuclear modification is visible over the measured
$p_{\text{T}}$ range. The measurement of $\pi^{0}$ $R_{\text{pPb}}$ [2] and
$\omega$ $R_{\text{dAu}}$ at $\sqrt{s_{\scalebox{0.5}{NN}}}$ = 200 GeV are
both in agreement with the nuclear modification factor presented here within
uncertainties.

(a) Measured production cross section of the $\omega$ meson in pp
and p–Pb collisions and selected MC generator predictions

(b) $\omega/\pi^{0}$ ratio measurements

(c) $\omega$ $R_{\text{pPb}}$ compared to $\pi^{0}$ $R_{\text{pPb}}$ [2] and
$\omega$ $R_{\text{dAu}}$
at $\sqrt{s_{\scalebox{0.5}{NN}}}$ = 200 GeV [11]

Figure 4: Measurements of $\omega$ meson production at
$\sqrt{s_{\scalebox{0.5}{NN}}}=5.02\leavevmode\nobreak\ $TeV.

##  5 Conclusion

The $\omega$ meson production cross section was measured in pp and p–Pb
collisions at $\sqrt{s_{\scalebox{0.5}{NN}}}=5.02\leavevmode\nobreak\ $TeV.
An agreement between the corresponding $\omega/\pi^{0}$ ratios in pp and p–Pb
collisions and a previous measurement at $\sqrt{s}=13\leavevmode\nobreak\ $TeV
suggest the production ratio to be independent of the collision system and
energy. A first measurement of the nuclear modification factor of the $\omega$
mesons presented in this contribution shows no visible nuclear modification
over the measured $p_{\text{T}}$ range.

This contribution was supported by BMBF and the Helmholtz Association.

## References

  * [1] ALICE Collaboration. 2012. Neutral pion and $\eta$ meson production in proton-proton collisions at $\sqrt{s}=0.9$ TeV and $\sqrt{s}=7$ TeV. Phys. Lett. B 717 (2012) 162
  * [2] ALICE Collaboration. 2018. Neutral pion and $\eta$ meson production in p–Pb collisions at $\sqrt{s_{\scalebox{0.5}{NN}}}\leavevmode\nobreak\ =\leavevmode\nobreak\ 5.02\leavevmode\nobreak\ $TeV. Eur. Phys. J. C (2018) 78:624
  * [3] ALICE Collaboration. 2018. $\pi^{0}$ and $\eta$ meson production in proton-proton collisions at $\sqrt{s}=8$ TeV. Eur. Phys. J. C 78 (2018) 263
  * [4] ALICE Collaboration. 2015. Measurement of pion, kaon and proton production in proton-proton collisions at $\sqrt{s}\leavevmode\nobreak\ =\leavevmode\nobreak\ 7\leavevmode\nobreak\ $TeV. EPJC 75 (2015) 226
  * [5] LHCb Collaboration. 2022. Nuclear modification factor of neutral pions in the forward and backward regions in pPb collisions. https://arxiv.org/abs/2204.10608
  * [6] ALICE Collaboration. 2020. Production of $\omega$ mesons in pp collisions at $\sqrt{s}\leavevmode\nobreak\ =7\leavevmode\nobreak\ $TeV. Eur. Phys. J. C 80 (2020) 1130
  * [7] H. Saveetha and D. Indumathi. 2017. Fragmentation of $\omega$ and $\phi$ mesons in $e^{+}e^{-}$ and pp collisions at NLO. International Journal of Modern Physics A 32.33 (2017), p. 1750199.
  * [8] ALICE Collaboration. 2016. Direct photon production in Pb-Pb collisions at $\sqrt{s_{\scalebox{0.5}{NN}}}\leavevmode\nobreak\ =\leavevmode\nobreak\ 2.76\leavevmode\nobreak\ $TeV. Phys. Lett. B 754 (2016) 235-248
  * [9] P.A. Zyla et al. (Particle Data Group), Prog. Theor. Exp. Phys. 2020, 083C01 (2020)
  * [10] Luca Lista. 2017. Combination of measurements and the BLUE method. EPJ Web of Conferences 137 (2017). Ed. by Y. Foka, N. Brambilla, and V. Kovalenko, p. 11006
  * [11] PHENIX Collaboration. 2007. Production of $\omega$ mesons at large transverse momenta in p+p and d+Au collisions at $\sqrt{s_{\scalebox{0.5}{NN}}}\leavevmode\nobreak\ =\leavevmode\nobreak\ 200\leavevmode\nobreak\ $GeV. Phys. Rev. C 75(5):051902(R)


# A Framework of Zero-Inflated Bayesian Negative Binomial Regression Models
For Spatiotemporal Data

Qing He  Hsin-Hsiung Huang

###### Abstract

Spatiotemporal data analysis with massive zeros is widely used in many areas
such as epidemiology and public health. We use a Bayesian framework to fit
zero-inflated negative binomial models and employ a set of latent variables
from Pólya-Gamma distributions to derive an efficient Gibbs sampler. The
proposed model accommodates varying spatial and temporal random effects
through Gaussian process priors, which have both the simplicity and
flexibility in modeling nonlinear relationships through a covariance function.
To conquer the computation bottleneck that GPs may suffer when the sample size
is large, we adopt the nearest-neighbor GP approach that approximates the
covariance matrix using local experts. For the simulation study, we adopt
multiple settings with varying sizes of spatial locations to evaluate the
performance of the proposed model such as spatial and temporal random effects
estimation and compare the result to other methods. We also apply the proposed
model to the COVID-19 death counts in the state of Florida, USA from 3/25/2020
through 7/29/2020 to examine relationships between social vulnerability and
COVID-19 deaths.

###### keywords:

nearest-neighbor Gaussian process , negative binomial model , Pólya-Gamma
scheme , social vulnerability , spatiotemporal effects
MSC: 62F15, 62C10, 91B72

††journal: Journal of Statistical Planning and Inference

[inst1]organization=Department of Statistics and Data Science,addressline=
University of Central Florida, city=Orlando, postcode=32816, state=Florida,
country=USA

##  1 Introduction

Zero-inflated models have been widely used for handling count data with
excessive zeros (Lewsey and Thomson, 2004, Cheung, 2002). By construction,
zero-inflated models assume that the zeros come from one of the cases (Greene,
1994): 1) structural zeros corresponding to individuals who are not at risk
for an event, and therefore have no opportunity for a positive count and 2)
random zeros which correspond to a latent class of individuals who are at risk
for an event but nevertheless have an observed response of zero. A zero-
inflated negative binomial (ZINB) model is a popular choice for modeling zero-
inflated data as it gives more reliable parameter estimates when the nonzero
counts are over-dispersed compared to other models like the zero-inflated
Poisson model (Yau et al., 2003). Bayesian approaches to fitting zero-inflated
models have gained attention recently (Ghosh et al., 2006, Neelon et al.,
2010, Zuur et al., 2012). There is an abundance of advantages of Bayesian
zero-inflated models including using prior beliefs with Bayesian inference
instead of deriving asymptotically approximate distributions for estimation of
parameters (Gelman et al., 1995). Bayesian approaches produce tractable
inference for ZINB models and have been implemented through software such as R
(R Core Team, 2021), Stan (Stan Development Team, 2021), and WinBUGS (Lunn et
al., 2013).

It is important to account for spatial and temporal structures in areas such
as epidemiology and public health. Within the Bayesian framework, it is
natural to build a Bayesian hierarchical model with a prior distribution for
spatial and temporal effects (Kang and Cressie, 2011). Distance-based
exponential or Matérn covariance functions are commonly adopted to specify
spatial correlations in geostatistical data, while conditionally
autoregressive models are often used for temporal effects (Cressie and Wikle,
2015, Banerjee et al., 2003, Zhang et al., 2022) through a hierarchical
modeling framework (Cressie and Wikle, 2015). Neelon et al. (2022) proposed a
cubic B-spline model to fit both temporal (overall time trend) and spatial
(county-specific trend) random effects in a Bayesian framework, which depends
on the number and locations of the interior knots and specification of the
degree of basis functions.

There is limited research with the focus on the Bayesian ZINB model with the
consideration of flexible spatial and temporal effects (Wang et al., 2016, Gu
et al., 2020). Neelon (2019) proposed a Bayesian ZINB model with spatial and
temporal random effects but assumed that temporal effects changed in a fixed
increment, which may not reflect the randomness nature of temporal effects. In
this study, we focus on proposing a flexible approach to model spatial and
temporal effects in a Bayesian ZINB framework. A promising option is using a
Gaussian process (GP) prior which characterize spatial or temporal effects
through a kernel function and enables tractable nonparametric Bayesian
inference. A GP prior is a probability distribution over infinite number of
possible functions which make it a powerful tool to model nonlinear patterns.
GP accounts for quantification of uncertainty, which relies on the predictive
conditional distribution of a new spatial location given a set of observed
locations (Cousin et al., 2016). GP can be viewed as a spline in a reproducing
kernel Hilbert space with the reproducing kernel of a covariance function
(Kimeldorf and Wahba, 1970, Wahba, 1990). It has become a popular modeling
tool in multivariate and geostatistical settings. For example, Diana et al.
(2021) proposed a Bayesian hierarchical occupancy model using GPs within a
logistic regression framework for spatial and temporal random effects.
However, the covariance matrix used in GP models may hinder its implementation
for large datasets due to computational issues. Recent research has developed
various approximation techniques to overcome the computational limitations
(Smola and Bartlett, 2001, Quinonero-Candela and Rasmussen, 2005, Snelson and
Ghahramani, 2007, Lázaro-Gredilla et al., 2010). There are two broad classes
of these techniques, i.e., global approximation methods and local
approximation methods (Liu et al., 2020). For global approximation methods,
sparsity of the covariance matrix is achieved through global distillation
including using a subset of the training data, exploiting the sparse structure
of the matrix, and low-rank models (Chalupka et al., 2013, Quinonero-Candela
and Rasmussen, 2005, Titsias, 2009). For example, low-rank models embed the
original process into a lower-dimensional subspace. However, it usually
requires a large subspace when used for large spatial datasets. Local
approximation methods are based on multiple sets of local experts to improve
the scalability based on the concept of divide-and-conquer (Liu et al., 2020,
Gramacy, 2016). This class of approximations has the capability of capturing
nonstationary features from the localized experts. Depending on the partition
of the input space, there are inductive local experts which employs static
partitioning (Vasudevan et al., 2009) and transductive local experts which
employ a dynamic partition (Datta et al., 2016a, b). As an example of
transductive local experts, Datta et al. (2016a) proposed a nearest-neighbor
Gaussian process (NNGP) that uses conditional independence given information
from neighboring points for large geostatistical datasets. It has well-defined
sparse precision matrices for its finite-dimensional realization and provides
fully process-based inference on underlying spatial processes.

In this paper, we adopt the NNGP model (Datta et al., 2016a) to model the
spatial and temporal effects because of its ease of use and well-defined
sparse precision matrix. The proposed Bayesian zero-inflated negative binomial
(ZINB) model uses nearest-neighbor Gaussian process priors to fit spatial and
temporal random effects that account for the over-dispersed and zero-inflated
count response by incorporating the spatial and temporal covariance matrices.
It allows us to obtain more reliable and stable inferences since NNGP yields a
lower-dimensional spatial and temporal covariance matrices and parameter
space, therefore requires much smaller matrix inversions. The propose
framework is a scalable hierarchical ZINB model with flexible spatial and
temporal effects through GP. First, the hierarchical model is based on work in
Neelon (2019) to estimate regression parameters and the spatial and temporal
effects. A key component to the parameter estimator is to introduce a set of
latent variables through the Pólya-Gamma (PG) data augmentation, an underlying
latent variable methodology (Pillow and Scott, 2012, Polson et al., 2013).
With the PG scheme, posterior inference is highly efficient via Gibbs sampling
by updating from conditional distributions with conjugate priors and explicit
conditional normal distributions. To illustrate Gibbs sampling for generating
a sequence of samples, we use an example of a pair of random variables
$(\mathcal{X},\mathcal{Y})$. A sequence of random variables

$\mathcal{X}^{\prime}_{0},\mathcal{Y}^{\prime}_{0},\mathcal{X}^{\prime}_{1},\mathcal{Y}^{\prime}_{1},\ldots,\mathcal{X}^{\prime}_{t},\mathcal{Y}^{\prime}_{t},\ldots,$

is generated iteratively as follows: 1. set $\mathcal{X}^{\prime}_{0}$ at some
initial value; 2. sample the rest from $\mathcal{Y}^{\prime}_{j}\sim
f(y\mid\mathcal{X}^{\prime}_{j}=x^{\prime}_{j})$ and
$\mathcal{X}^{\prime}_{j+1}\sim
f(x\mid\mathcal{Y}^{\prime}_{j}=y^{\prime}_{j})$ alternately (Casella and
George, 1992). It has been shown that the PG method yields superior
performance compared to other Bayesian methods in the context of structural
equation models with logistic regression for binary variables (Kim et al.,
2018) and the Bayesian ZINB (Neelon, 2019). Therefore we propose a more
flexible structure for the spatial and temporal random effects through the GP
prior in order to achieve fully model-based Bayesian inference for the spatial
and temporal effects. More specifically, to facilitate the computation
process, we incorporate the framework with a conjugate latent NNGP model
(Datta et al., 2016a, b, Zhang et al., 2021b), which is scalable for massive
spatial datasets on modest computing environments (Datta et al., 2016a, b,
Zhang et al., 2021b).

##  2 Bayesian Inference for ZINB Model Using NNPG

###  2.1 The Pólya-Gamma Scheme

Let $y_{i}$ as the outcome of the $i$th observation, which may be a count
$y_{i}=0,1,2,\ldots$ or a binary indicator $y_{i}=0,1$. We aim to model
$Y_{i}$ and a set of covariates $X_{i}=[1,x_{i1},\ldots,x_{iP}]^{T}$ through
$\mbox{E}[Y_{i}\mid X_{i}]=g^{-1}(X_{i}\beta),$ where
$\beta=[\beta_{0},\ldots,\beta_{P}]^{T}$ are the regression coefficients and
$g$ is the canonical link function in generalized linear models (GLMs)
(McCullagh and Nelder, 2019). We utilize the efficient PG scheme for
estimating the regression coefficients in GLMs. An advantage of the PG scheme
is that the posterior distributions of the parameters of interest enable us to
make use of Gibbs sampling approach within logistic regression and NB
regression frameworks (Polson et al., 2013) through an augmented PG
distributed random variable $\omega$. Generally, assuming a prior distribution
for regression coefficient $\beta\sim\mathcal{N}(b,\Sigma_{0})$, the full
conditional distributions induced by PG for the Gibbs sampler are:

$\displaystyle\omega_{i}\mid\beta\sim\mathcal{PG}(n_{i},X_{i}\beta),\mbox{ and
}\beta\mid\Omega,z\sim\mathcal{N}(\mu,\Sigma),$

where

$\displaystyle\Sigma=(\Sigma_{0}^{-1}+X^{T}\Omega X)^{-1},\mbox{ and
}\mu=\Sigma\left(\Sigma_{0}^{-1}b+X^{T}\Omega z\right).$

Here $\mathcal{PG}$ denotes a PG distribution,
$\Omega=\text{diag}(\omega_{1},\ldots,\omega_{n})$ and $z$ is the latent
variable with $z_{i}=\frac{y_{i}-n_{i}/2}{\omega_{i}}$. Specifically, we use
the PG scheme for estimating the regression coefficients in the Bernoulli and
NB parts of the ZINB regression model, respectively.

In the binomial logistic regression model (logit model), let $n_{i}=1\;\forall
i$ and $y_{i}\sim\text{Binom}(1,\frac{1}{1+e^{-\eta_{i}}})$ where
$\eta_{i}=x_{i}^{T}\beta$ . Following Polson et al. (2013), the full
conditional for $\beta$ in a binomial regression is

$p(\beta\mid
y,r,\omega)\propto\pi(\beta)\exp\left[-\frac{1}{2}(z-X\beta)^{T}\Omega(z-X\beta)\right],$

where $\pi(\beta)$ is the prior distribution, $z$ is an $n\times 1$ vector
with $z_{i}=\frac{y_{i}-1/2}{\omega_{i}}$ and $\Omega=$ diag$(\omega_{i})$. It
is readily seen that $z\mid\beta,\Omega\sim\mathcal{N}(X\beta,\Omega)$, which
leads to the convenient Gibbs sampler for $\beta$. In the negative binomial
regression model (Pillow and Scott, 2012), $n_{i}=y_{i}+r$ and $r$ can be
obtained using either the Metropolis–Hastings algorithm in which a uniform
prior is applied with positive candidate values of $r$ are drawn from a zero-
truncated normal proposal centered at the current value of $r$, or the two-
stage Gibbs sampling (Zhou and Carin, 2013, Dadaneh et al., 2018). Consider
the following model for a count response $Y_{i}$,

$\displaystyle p(Y_{i}=y_{i}\mid r,\beta)$
$\displaystyle\stackrel{{\scriptstyle
d}}{{=}}\frac{\Gamma(y_{i}+r)}{\Gamma(r)y_{i}!}(1-\psi_{i})^{r}(\psi_{i})^{y_{i}},\;r>0,\;\mbox{
where }$ $\displaystyle\psi_{i}$
$\displaystyle=\frac{\exp(x^{T}_{i}\beta)}{1+\exp(x^{T}_{i}\beta)}=\frac{\exp(\eta_{i})}{1+\exp(\eta_{i})},$

where the NB probability parameter $\psi_{i}$ is parameterized using the expit
(inverse-logit) function, which allows us to apply the same properties of the
Pólya-Gamma density as in the logistic case.

The mean and variance of $Y_{i}$ are

$\displaystyle\mbox{E}(Y_{i}\mid r,\beta)$
$\displaystyle=\frac{\psi_{i}}{1-\psi_{i}}=r\exp(\eta_{i})=\mu_{i},$
$\displaystyle\mbox{Var}(Y_{i}\mid r,\beta)$
$\displaystyle=\frac{r\psi_{i}}{(1-\psi_{i})^{2}}=r\exp(\eta_{i})[1+\exp(\eta_{i})]=\mu_{i}(1+\mu_{i}/r).$

The parameter $r$ accounts for overdispersion in the data. As $r\to\infty$,
the counts become increasingly dispersed and the NB distribution converges to
the Poisson distribution. The above parameterization also leads to a Gaussian
linear model where $z\mid\beta,\Omega\sim\mathcal{N}(X\beta,\Omega)$ where $z$
is an $n\times 1$ vector with $z_{i}=\frac{y_{i}-r}{2\omega_{i}}$,
$\omega_{i}\sim\mathcal{PG}(y_{i}+r,X_{i}\beta)$ and $\Omega=$
diag$(\omega_{i})$.

##  3 Zero-inflated Negative Binomial Model

###  3.1 The Proposed Spatiotemporal ZINB-NNGP Model

We are interested in a spatiotemporal random intercept model on data with
excessive zeros. Assume observations are collected at a fixed collection of
distinct locations denoted by $\mathbf{S}=\\{1,\ldots,S\\}$ and across $T$
distinct time segments (e.g., weeks, months, seasons and years) denoted by
$\mathbf{T}=\\{1,\ldots,T\\}$. We index the observations by $j=1,\ldots N$,
where $N$ is the total sample size. For a sampling unit $\\{s,t\\}$ where
$s\in\mathbf{S}$ and $t\in\mathbf{T}$, we denote its number of observations by
$n_{s,t}$ ($n_{s,t}\geq 0$). In a special case, $n_{s,t}=1,\forall
s\in\mathbf{S},t\in\mathbf{T}$, that is there is exactly one observation for
each sampling unit $\\{s,t\\}$ and $N=S\times T$. In the proposed framework,
it does not require the sampling units to have equal sample sizes or at least
one observation. Sparsity and unequal sizes are considered in the simulation
study.

Let response $Y_{j}$ be the observed count of the $j$th observation,
$\phi_{j}$ be the probability of belonging to the at-risk group, and
$\psi_{j}$ be the success rate for observations in the at-risk group. The
Bayesian ZINB model (Neelon, 2019) with a latent at-risk indicator variable
$W_{j}$ is introduced as follows:

$Y_{j}\sim(1-\phi_{j})1_{(W_{j}=0\land
Y_{j}=0)}+\phi_{j}\mathcal{NB}(\mu_{j},r)1_{(W_{j}=1)},\;j=1,\ldots,N,$

where $\mu_{j}$ is the mean of the NB distribution and $r$ is the
overdispersion parameter. A value of $1$ for $W_{j}$ means that the $j$th
observation belongs to the at-risk group, and a value of $0$ means the “not
at-risk” group. Let $s_{j}$ and $t_{j}$ be its spatial location and time
point, respectively. We model the at-risk probability $\phi_{j}$ using a logit
model and $Y_{j}\mid W_{j}=1$ using a negative binomial regression model:

$\displaystyle Pr(Y_{j}=0)$
$\displaystyle=(1-\phi_{j})+\phi_{j}(1-\psi_{j})^{r},$
$\displaystyle\mbox{logit}(\phi_{j})$
$\displaystyle=\mbox{logit}[Pr(W_{j}=1\mid\alpha,a_{s_{j}},b_{t_{j}},\epsilon_{11s_{j}},\epsilon_{12t_{j}})]$
$\displaystyle=x^{T}_{j}\alpha+a_{s_{j}}+b_{t_{j}}+\epsilon_{11s_{j}}+\epsilon_{11t_{j}}$
$\displaystyle=\eta_{1j},$ (1) $\displaystyle p(y_{j}\mid
r,\beta,c_{s_{j}},d_{t_{j}},\epsilon_{21s_{j}},\epsilon_{22t_{j}},W_{j}=1)$
$\displaystyle\stackrel{{\scriptstyle
d}}{{=}}\frac{\Gamma(y_{j}+r)}{\Gamma(r)y_{j}!}(1-\psi_{j})^{r}\psi_{j}^{y_{j}},\;\forall
j\;\mbox{ s.t. }W_{j}=1,\;$ $\displaystyle\psi_{j}$
$\displaystyle=\frac{\exp(x^{T}_{j}\beta+c_{s_{j}}+d_{t_{j}}+\epsilon_{21s_{j}}+\epsilon_{22t_{j}})}{1+\exp(x^{T}_{j}\beta+c_{s_{j}}+d_{t_{j}}+\epsilon_{21s_{j}}+\epsilon_{22t_{j}})}$
$\displaystyle=\frac{\exp(\eta_{2j})}{1+\exp(\eta_{2j})},$ (2)

where $a_{s_{j}}$ and $\epsilon_{11{s_{j}}}$ are the spatial random effects
and spatial random noise in the logit model, $b_{t_{j}}$ and
$\epsilon_{12{t_{j}}}$ are the temporal random effects and temporal random
noise in the logit model, $c_{s_{j}}$ and $\epsilon_{21{s_{j}}}$ are the
spatial random effects and spatial random noise in the negative binomial
regression model, and $d_{t}$ and $\epsilon_{22{t_{j}}}$ are the temporal
random effects and temporal random noise in the negative binomial model. The
two models are often referred as the binary component (the logit model) and
count component (the negative binomial model), which we will use extensively
in the paper. The proposed model incorporates additive spatial and temporal
random effects as intercepts in the binary and count component. Compared with
the existing Bayesian ZINB model (Neelon, 2019), the proposed model is more
flexible through the use of GP priors on the spatial and temporal random
effects. To elaborate the GP priors, first we define the spatial random
effects corresponding to the sampling locations $\mathbf{S}=\\{1,\ldots,S\\}$
as

$\displaystyle\vec{a}$ $\displaystyle=(a_{1},\ldots,a_{S})^{T},\mbox{in the
binary component and }$ $\displaystyle\vec{c}$
$\displaystyle=(c_{1},\ldots,c_{S})^{T},\mbox{in the count component}.$

Similarly, we define the temporal random effects of temporal points
$\mathbf{T}=\\{1,\ldots,T\\}$ as

$\displaystyle\vec{b}$ $\displaystyle=(b_{1},\ldots,b_{T})^{T},\mbox{in the
binary component and}$ $\displaystyle\vec{d}$
$\displaystyle=(d_{1},\ldots,d_{T})^{T},\mbox{in the count component}.$

A list of descriptions of the notations is given in Table 1.

Notation |  Description
---|---
$S$ |  the number of spatial locations
$T$ |  the number of temporal points
$\mathbf{S}=\\{1,2,\ldots,S\\}$ |  the list of spatial locations
$\mathbf{T}=\\{1,2,\ldots,T\\}$ |  the list of temporal points
$N$ |  the number of observations
$Y_{j},j=1,\ldots,N$ |  the response of the $j$th observation
$s_{j}\in\mathbf{S}$ |  the spatial location of the $j$th observation
$t_{j}\in\mathbf{T}$ |  the temporal point of the $j$th observation
$\phi_{j}$ |  the probability of belonging to the at-risk group
$\psi_{j}$ |  the success rate for observations in the at-risk group
$W_{j}$ |  the indicator for belonging to the at-risk group
$\mu_{j}$ |  the mean of negative binomial distribution
$r$ |  the dispersion parameter of negative binomial distribution
$a_{s_{j}}$ |  the spatial random effect in the logit model
$\epsilon_{11s_{j}}$ |  the spatial random noise in the logit model
$b_{t_{j}}$ |  the temporal random effect in the logit model
$\epsilon_{12t_{j}}$ |  the temporal random noise in the logit model
$c_{s_{j}}$ |  the spatial random effect in the negative binomial model
$\epsilon_{21s_{j}}$ |  the spatial random noise in the negative binomial model
$d_{t_{j}}$ |  the temporal random effect in the negative binomial model
$\epsilon_{22t_{j}}$ |  the temporal random noise in the negative binomial model
Table 1: The description of the notations used in the proposed model.

To account for correlations between spatial locations or temporal points, we
employ GPs for both the spatial and temporal effects. Specifically, we assume
that $\vec{a}=(a_{1},\ldots,a_{S})^{T}$ in the binary component is distributed
according to a GP with parameters $(\sigma_{11},l_{11})$. Let $h_{s}$ be the
specific location information of site $s$ (e.g., latitude and longitude) which
can be used to compute distance between different locations. It corresponds to
assuming that
$(a_{1},\ldots,a_{S})^{T}\sim\mathcal{N}(0,K_{\sigma_{11},l_{11}}(h_{1},\ldots,h_{S}))$,
where
$K_{\sigma_{11},l_{11}}(h_{1},\ldots,h_{S})_{i,j}=\sigma_{11}^{2}e^{-\frac{|h_{i}-h_{j}|^{2}}{l_{11}^{2}}},i,j=1,\ldots,S$.
Parameter $\sigma_{11}$ determines the overall variability of the GP, while
parameter $l_{11}$ measures the correlation between location $i$ and $j$.
Similarly, the autocorrelated temporal random effect
$\vec{b}=(b_{1},\ldots,b_{T})^{T}$ in the count component is distributed
according to a GP with parameters $(\sigma_{12},l_{12})$. Let $w_{t}$ be the
index of temporal point $t$. For example, if the data are at monthly level
spanning from January to December, then the corresponding indices could be
$(1,\ldots,12)$. The GP prior can be written as
$(b_{1},\ldots,b_{T})^{T}\sim\mathcal{N}(0,K_{\sigma_{12},l_{12}}(w_{1},\ldots,w_{T}))$.

Let $V_{1}$ and $V_{2}$ be the $N\times S$ and $N\times T$ design matrices for
the spatial and temporal random effects for the observed data. The following
hierarchical structure completes the definition of the binary component with
spatial random effect $\vec{a}=(a_{1},\ldots,a_{S})^{T}$, temporal random
effect $\vec{b}=(b_{1},\ldots,b_{T})^{T}$, spatial random noise
$\vec{\epsilon}_{11}=(\epsilon_{111},\ldots,\epsilon_{11S})^{T}$, and temporal
random noise $\vec{\epsilon}_{12}=(\epsilon_{121},\ldots,\epsilon_{12T})^{T}$
including the prior distributions of relevant parameters. The hierarchical
structure for the negative binomial regression with random effects
$\vec{c}=(c_{1},\ldots,c_{S})^{T}$ , $\vec{d}=(d_{1},\ldots,d_{T})^{T}$,
$\vec{\epsilon}_{21}=(\epsilon_{211},\ldots,\epsilon_{21S})^{T}$, and
$\vec{\epsilon}_{22}=(\epsilon_{221},\ldots,\epsilon_{22T})^{T}$ are in the
same fashion.

$\displaystyle(a_{1},\ldots,a_{S})^{T}$
$\displaystyle\sim\mathcal{N}(0,K_{\sigma_{11},l_{11}}(h_{1},\ldots,h_{S})),$
$\displaystyle\sigma^{2}_{11}$
$\displaystyle\sim\mathcal{IG}(a_{\sigma_{1}},b_{\sigma_{1}}),\;l_{11}\sim
Gamma(a_{l_{1}},b_{l_{1}}),$ $\displaystyle(b_{1},\ldots,b_{T})^{T}$
$\displaystyle\sim\mathcal{N}(0,K_{\sigma_{12},l_{12}}(w_{1},\ldots,w_{T})),$
$\displaystyle\sigma^{2}_{12}$
$\displaystyle\sim\mathcal{IG}(a_{\sigma_{2}},b_{\sigma_{2}}),\;l_{12}\sim
Gamma(a_{l_{2}},b_{l_{2}}),$ $\displaystyle\epsilon_{11s}$
$\displaystyle\stackrel{{\scriptstyle
iid}}{{\sim}}\mathcal{N}(0,\sigma^{2}_{\epsilon_{11}}),\;\sigma^{2}_{\epsilon_{11}}\sim\mathcal{IG}(a_{\epsilon},b_{\epsilon}),s=1,\ldots,S,$
$\displaystyle\epsilon_{12t}$ $\displaystyle\stackrel{{\scriptstyle
iid}}{{\sim}}\mathcal{N}(0,\sigma^{2}_{\epsilon_{12}}),\;\sigma^{2}_{\epsilon_{12}}\sim\mathcal{IG}(a_{\epsilon},b_{\epsilon}),t=1,\ldots,T.$

####  3.1.1 Latent Nearest-Neighbor Gaussian Process

In a spatial regression model at location $s\in\mathbf{S}$ and time
$t\in\mathbf{T}$ in a spatiotemporal domain $\mathcal{D}$, assume that

$\displaystyle
y(s,t)=\mu+w_{1}(s)+w_{2}(t)+\epsilon_{1}(s)+\epsilon_{2}(t),\;(s,t)\in\mathcal{D},$

where $\mu$ is the mean function that is itself modeled as a linear
combination of known covariates, $w_{1}(s)$ follows a latent GP with mean zero
and a positive-definite cross-covariance matrix $C_{\phi}(s,s^{\prime})$, and
$\epsilon_{1}(s)$ follows a normal distribution with mean zero and variance
$\sigma^{2}$. In the same fashion, we can define the temporal random effect
$w_{2}(t)$ and the temporal random noise $\epsilon_{2}(t)$. This type of
regression model is called a latent GP model since the GP is used as a prior
in the latent process $w(s)$ and $w(t)$, while a response GP model imposes a
GP on the outcome $y(s,t)$. In the proposed Bayesian spatiotemporal ZINB model
(2), we use the latent GP model with both spatial and temporal random effects.

In NNGP (Datta et al., 2016b), the joint density of $(w_{1},\cdots,w_{S})$ was
approximated by a series of conditional densities of size at most $m$, where
$m\leq S$ and $S$ is the size of distinct spatial locations. It was shown that
the approximation was essentially a multivariate Gaussian density with
covariance matrix $\tilde{C}_{\phi}(\mathbf{S},\mathbf{S})$. Using the spatial
random effects $w_{1}(\mathbf{S})$ described above as an example, by selecting
a list of nearest neighbors with size $m$ for each location in $\mathbf{S}$,
the Gaussian density $p(w_{1}(\mathbf{S}))$ can be approximated as follows
(Datta et al., 2016b, Zhang et al., 2021b):

$\displaystyle
p(w_{1}(\mathbf{S}))=N(w_{1}(\mathbf{S})\mid\mathbf{0},C_{\phi}(\mathbf{S},\mathbf{S}))\approx
N(w_{1}(\mathbf{S})\mid\mathbf{0},\tilde{C}_{\phi}(\mathbf{S},\mathbf{S})).$
(3)

The approximation $\tilde{C}_{\phi}(\mathbf{S},\mathbf{S})$ is computationally
efficient because its inverse $\tilde{C}_{\phi}(\mathbf{S},\mathbf{S})^{-1}$
is sparse and can be written as:

$\tilde{C}_{\phi}(\mathbf{S},\mathbf{S})^{-1}=(I-A_{S})^{T}D^{-1}_{S}(I-A_{S}),$

where $A_{S}$ is a lower triangular matrix and $D_{S}$ is a diagonal matrix.
To compute ${A}_{S}$, suppose $N_{m}(i)$ be the set of column indices of the
$m$ nearest neighbors that contain nonzero entries in the $i$th row of
${A}_{S}$, where $m=1,\ldots,S-1$ is a hyperparameter controlling the size of
nearest neighbors. Let ${A}_{S}=[{a}_{1}:\cdots:{a}_{S}]^{T}$ and
$D_{S}=\text{diag}(d_{1},d_{2},\ldots,d_{S})$. The first row of ${A}_{S}$ has
all elements equal to 0 and $d_{1}=C_{\phi}(1,1)$. For $i=2,\ldots,S$, we
obtain the nonzero entries at column positions indexed by $N_{m}(i)$ in
$A_{S}$ and the diagonal elements in $D_{S}$ as follows:

$\displaystyle a_{i}(N_{m}(i))$
$\displaystyle=C_{\phi}(i,N_{m}(i))C_{\phi}(N_{m}(i),N_{m}(i))^{-1}\mbox{ and
}$ (4) $\displaystyle d_{i}$
$\displaystyle=C_{\phi}(i,i)-C_{\phi}(i,N_{m}(i))C_{\phi}(N_{m}(i),N_{m}(i))^{-1}C_{\phi}(N_{m}(i),i).$
(5)

####  3.1.2 Spatio-temporal correlations via NNGP in ZINB

In this section, we describe how to embed spatiotemporal correlations through
NNGP in the proposed Bayesian ZINB framework. As we assume GP priors for both
the spatial and temporal random effects in the proposed model, the latent NNGP
can be used for one random process or both when needed. Since the latent NNGP
is designed for a spatial process initially, we use the spatial random effects
to demonstrate the application of latent NNGP for easier understanding. When
the number of spatial locations $S$ is large, the computation bottleneck in
the conjugate Gaussian regression model lies in computing
$K_{\sigma_{11},l_{11}}^{-1}$ and $K_{\sigma_{21},l_{21}}^{-1}$ when updating
$\sigma_{11}$ and $\sigma_{21}$. To facilitate it, we can find a sparse
alternative for the Cholesky decomposition of $K_{\sigma_{11},l_{11}}^{-1}$
and $K_{\sigma_{21},l_{21}}^{-1}$. We use the covariance matrix for the
spatial random effects from the binary component as an example to illustrate
the methodology. First, we rewrite the covariance matrix
$C_{\phi}(\mathbf{S},\mathbf{S})=K_{\sigma_{11},l_{11}}(h_{1},\ldots,h_{S})=\sigma_{11}\rho_{l_{11}}$.
Let $\tilde{\rho}_{l_{11}}$ be the NNGP approximation of $\rho_{l_{11}}$ with
its inverse $\tilde{\rho}^{-1}_{l_{11}}=(I-A_{S})^{T}D^{-1}_{S}(I-A_{S})$, as
described in Section 3.1.1. Nonzero entries in the $i$th row where
$i=2,\ldots,N$ of $A_{S}$ and $D_{S}$ in (4) and (5) become

$\displaystyle a_{i}(N_{m}(i))$
$\displaystyle=\rho_{l_{11}}(i,N_{m}(i))(\rho_{l_{11}}(N_{m}(i),N_{m}(i)))^{-1}$
(6)

and

$\displaystyle d_{i}$
$\displaystyle=1-a_{i}(N_{m}(i))\rho_{l_{11}}(N_{m}(s_{i}),i).$ (7)

In this same fashion, we model the temporal correlation. The latent NNGP is
applied to facilitate to computation bottleneck when updating the
hyperparameters for the spatial and temporal random effects. Hence, when the
spatial size and/or temporal size are large, the latent NNGP can be used.
Otherwise, one may just compute the inverse of the covariance matrix. The
step-by-step posterior sampling procedure is described in the appendix.

##  4 Simulated and Real Data Analyses

In this section, we apply the proposed Bayesian ZINB-NNGP model to both
simulated and COVID-19 fatality data. For the simulation data, a total number
of $N$ observations are generated from $S$ locations across $T$ time points.
First, we consider exactly one repetition for each sampling unit $(s,t)$, we
vary the spatial sizes to show if the proposed model can be generalized well
to a large spatial dimension. We set $S=200$ and $T=20$ in simulation 1 and
set $S=500$ and $T=20$ in simulation 2. Second, we relax the constraint on the
number of repetitions for each sampling unit by simulating a dynamic number of
repetitions for each sampling unit from a Poisson distribution. This is shown
in simulation 3.

For the COVID-19 analysis, the daily level data on the total number of
COVID-19 deaths in Florida counties in the early stage from 3/25/2020 through
7/29/2020 are analyzed. There are 69 counties ($S=69$) and 127 days ($T=127$)
in the dataset. The reported death count of county $s$ at day $t$, denoted by
$y_{s,t}$, is the response variable. We apply the proposed model to analyze
the association of COVID-19 death at-risk and count rates with social
vulnerability, other sociodemographic characteristics (i.e., health insurance
coverage, urbanicity), population health care resources (i.e., primary care
physicians), population health measures (i.e., life expectancy, obesity), and
population density.

###  4.1 Simulation Data

####  4.1.1 Simulation 1: a moderate spatial dimension and one repetition in
each sampling unit

In simulation 1, we have simulated a dataset with $S=200$ and $T=20$, which
results in a total sample size of $N=4,000$. For the fixed effect, we sample
$x\sim\mathcal{N}(0,1)$, set the intercept and slope in the binary component
at $\alpha_{0}=-0.25$ and $\alpha_{1}=0.25$, and the intercept and slope in
the count component as $\beta_{0}=0.5$ and $\beta_{1}=-0.25$. The spatial
locations are generated from a unit square. The latent NNGP approximation is
used when updating $\sigma_{11}$ and $\sigma_{21}$ for the spatial random
effects. For the size of nearest neighbors used in the latent NNGP model, we
have tested different values from 8 to 15 and found that the results are
relatively stable. In the results reported for both the simulations and COVID
analysis, we set the size of nearest neighbors as 13. The posterior means and
95% credible intervals as well as the true values of the parameters and
hyperparameters are presented in Table 2. We also use the Bayesian
spatiotemporal ZINB model with random effects proposed by Neelon (2019) and
ZINB regression models from R package ‘pscl’ (Jackman et al., 2015) for
comparison. Due to the difference in model structure, both comparative models
do not have results for the hyperparameters used in the proposed framework.

The results in Table 2 show that the 95% credible intervals of the proposed
model have covered the true values for all the parameters (i.e.,
$\alpha_{0},\alpha_{1},\beta_{0},\beta_{1}$), while the intervals of the
Bayesian ZINB model by Neelon (2019) do not cover the true values of
$\alpha_{0}$ and $\beta_{0}$ and the intervals of the traditional ZINB do not
cover the true values for parameters $\alpha_{0}$, $\alpha_{1}$ and
$\beta_{0}$. Figure 1 presents the map of the true and posterior mean spatial
random effects in the binary component (in the top panel) and the count
component (in the bottom panel). In both components, the fitted spatial
pattern closely mirrors the true distribution, suggesting the proposed model
recovers the underlying spatial pattern in the data. Figure 2 presents the
true and fitted temporal effects for each temporal point. The predicted
temporal patterns capture the trend of the true patterns in general.

Para. |  True Value |  Proposed Model |  Bayesian ZINB by Neelon (2019) |  Traditional ZINB
---|---|---|---|---
$\alpha_{0}$ |  -0.25 |  -0.07 (-1.68, 1.57) |  0.25 (-0.05, 0.68) |  0.62(0.43, 0.81)
$\alpha_{1}$ |  0.25 |  0.17 (0.06, 0.29) |  0.20 (0.07, 0.35) |  -0.24(-0.32, -0.16)
$\beta_{0}$ |  0.50 |  0.82 (-0.72, 2.47) |  0.92 (0.69, 1.14) |  -0.45(-0.91, 0.02)
$\beta_{1}$ |  -0.25 |  -0.20 (-0.27, -0.12) |  -0.20 (-0.28, -0.12) |  -0.22(-0.37, -0.07)
$l_{11}$ |  0.35 |  1.05 (0.64, 1.71) |  |
$\sigma_{11}$ |  0.50 |  1.81 (1.40, 2.31) |  |
$l_{12}$ |  1.00 |  3.45 (0.45, 4.92) |  |
$\sigma_{12}$ |  0.20 |  0.77 (0.43, 1.19) |  |
$l_{21}$ |  0.35 |  0.32 (0.25, 1.44) |  |
$\sigma_{21}$ |  0.50 |  1.21 (0.87, 2.15) |  |
$l_{22}$ |  1.00 |  2.55 (0.15, 4.63) |  |
$\sigma_{22}$ |  0.20 |  0.81 (0.42, 1.28) |  |
$\sigma_{\epsilon_{11}}$ |  0.05 |  0.05 (0.03, 0.08) |  |
$\sigma_{\epsilon_{12}}$ |  0.05 |  0.05 (0.03, 0.08) |  |
$\sigma_{\epsilon_{21}}$ |  0.05 |  0.05 (0.03, 0.08) |  |
$\sigma_{\epsilon_{22}}$ |  0.05 |  0.05 (0.03, 0.08) |  |
$r$ |  1.00 |  0.95 (0.63, 1.27) |  0.73 (0.52, 1.14) |
Table 2: Parameter estimates and 95% credible intervals for the proposed
model, the spatiotemporal ZINB model by Neelon (2019), and the traditional
ZINB model in simulation 1. Note that the comparative models have different
model assumptions and therefore do not have posterior results for
hyperparameters used in the paper. Figure 1: The true and fitted spatial
effects in simulation 1 using the proposed model. Each dot represents a
location and the color represent the scale of the spatial random effects.
Figure 2: The true and fitted temporal effects in simulation 1 using the
proposed model.

####  4.1.2 Simulation 2: a large spatial dimension and one repetition in each
sampling unit

In simulation 2, we use the same setting as simulation 1 and increase the
spatial dimension from $200$ to $500$, which results in a total sample size of
$N=10,000$. The latent NNGP approximation is used to update $\sigma_{11}$ and
$\sigma_{21}$ for the spatial random effects. The true and fitted results are
shown in Table 3. The 95% credible intervals of the proposed model have
covered the true values for all the parameters, while the intervals of the
traditional ZINB model do not cover the true values for most parameters. We
find that the Bayesian ZINB algorithm proposed by Neelon (2019) breaks in the
first few iterations due to singularity issues which are caused by the high
spatial dimension.

The true and posterior mean spatial and temporal random effects of simulation
2 are shown in Figure 3 and 4. Both show that the fitted random effects are
consistent with the true random effects.

Para. |  True Value |  Proposed Model |  Traditional ZINB
---|---|---|---
$\alpha_{0}$ |  -0.25 |  0.00 (-1.94, 1.84) |  0.32(0.14, 0.50)
$\alpha_{1}$ |  0.25 |  0.20 (0.13, 0.28) |  -0.23(-0.28, -0.17)
$\beta_{0}$ |  0.50 |  0.07 (-2.05, 2.14) |  -0.20(-0.59, 0.20)
$\beta_{1}$ |  -0.25 |  -0.23 (-0.28, -0.17) |  -0.23(-0.34, -0.12)
$l_{11}$ |  0.35 |  3.02 (2.13, 3.54) |
$\sigma_{11}$ |  0.50 |  2.87 (2.34, 3.45) |
$l_{12}$ |  1.00 |  3.92 (2.37, 4.94) |
$\sigma_{12}$ |  0.20 |  0.77 (0.44, 1.22) |
$l_{21}$ |  0.35 |  0.39 (0.20, 2.89) |
$\sigma_{21}$ |  0.50 |  1.74 (1.22, 3.41) |
$l_{22}$ |  1.00 |  3.05 (1.16, 4.79) |
$\sigma_{22}$ |  0.20 |  0.83 (0.43, 1.33) |
$\sigma_{\epsilon_{11}}$ |  0.05 |  0.05 (0.03, 0.07) |
$\sigma_{\epsilon_{12}}$ |  0.05 |  0.05 (0.03, 0.08) |
$\sigma_{\epsilon_{21}}$ |  0.05 |  0.05 (0.03, 0.08) |
$\sigma_{\epsilon_{22}}$ |  0.05 |  0.05 (0.03, 0.08) |
$r$ |  1.00 |  1.08 (0.80, 1.29) |
Table 3: Parameter estimates and 95% credible intervals for the proposed model
and the traditional ZINB model in simulation 2. Figure 3: The true and fitted
spatial effects in simulation 2 using the proposed model. Figure 4: The true
and fitted temporal effects in simulation 2 using the proposed model.

####  4.1.3 Simulation 3: dynamic number of repetitions in each sampling unit

In simulation 3, we run the simulation with a dynamic number of repetitions in
each sampling unit. We set $S=200$ and $T=20$ like simulation 1. For each
sampling unit $(s,t)$ the number of repetitions $n_{s,t}$ is randomly sampled
from a Poisson distribution with mean $2$, which results in a total sample
size of $N=7,867$. The posterior estimates and $95\%$ credible regions are
presented in Table 4. As the spatial size is relatively small and there are
repetitions in each sampling unit in simulation 3, the Bayesian ZINB (Neelon,
2019) is comparable to the proposed model. Figure 5 and 6 show the true and
fitted spatial and temporal random effects using the proposed model. The
fitted random effects are better than those in simulation 1 and 2 where there
is only one repetition per sampling unit.

Para. |  True Value |  Proposed Method |  Bayesian ZINB by Neelon (2019) |  Traditional ZINB
---|---|---|---|---
$\alpha_{0}$ |  -0.25 |  -0.27 (-2.00, 1.44) |  -0.34 (-0.48, -0.18) |  0.45 (0.30, 0.60)
$\alpha_{1}$ |  0.25 |  0.18 (0.11, 0.26) |  0.19 (0.11, 0.26) |  -0.20 (-0.26, -0.15)
$\beta_{0}$ |  0.50 |  0.31 (-1.68, 2.25) |  0.38 (0.2, 0.55) |  -0.05 (-0.34, 0.25)
$\beta_{1}$ |  -0.25 |  -0.20 (-0.26, -0.15) |  -0.20 (-0.25, -0.15) |  -0.21(-0.31, -0.12)
$l_{11}$ |  0.35 |  1.32(0.98, 1.86) |  |
$\sigma_{11}$ |  0.50 |  1.95 (1.53, 2.44) |  |
$l_{12}$ |  1.00 |  3.91 (2.04, 4.96) |  |
$\sigma_{12}$ |  0.20 |  0.74 (0.43, 1.14) |  |
$l_{21}$ |  0.35 |  0.39 (0.26, 1.77) |  |
$\sigma_{21}$ |  0.50 |  1.37 (0.86, 2.59) |  |
$l_{22}$ |  1.00 |  3.29 (0.99, 4.88) |  |
$\sigma_{22}$ |  0.20 |  0.80 (0.47, 1.24) |  |
$\sigma_{\epsilon_{11}}$ |  0.05 |  0.05 (0.03, 0.08) |  |
$\sigma_{\epsilon_{12}}$ |  0.05 |  0.05 (0.03, 0.08) |  |
$\sigma_{\epsilon_{21}}$ |  0.05 |  0.05 (0.03, 0.07) |  |
$\sigma_{\epsilon_{22}}$ |  0.05 |  0.05 (0.03, 0.08) |  |
$r$ |  1.00 |  1.19 (0.93, 1.46) |  1.13 (0.90, 0.55) |
Table 4: Parameter estimates and 95% credible intervals in simulation 3 using
the proposed model and the spatiotemporal ZINB model. Figure 5: The true and
fitted spatial effects in simulation 3 using the proposed model. Figure 6: The
true and fitted temporal effects in simulation 3 using the proposed model.

###  4.2 COVID-19 fatality rates modeling

COVID-19 has impacted populations around the world, with the fatality rate
varying dramatically across counties of Florida (Karmakar et al., 2021,
Khedhiri, 2021, Zhang et al., 2021a). Using statistical methods to identify
factors which are associated with COVID-19 fatality/mortality rates is one of
the most important topics. In the early stage of COVID-19, many counties were
at low risk of having positive cases and did not have reported deaths. Hence,
the proposed ZINB-NNGP model is a good option to model the mortality rate with
excessive zeros. In this section, we apply the proposed model to understand
the association of COVID-19 mortality at-risk and count rates with social
vulnerability, adjusting other sociodemographic characteristics (i.e., health
insurance coverage, urbanicity), population health care resources (i.e.,
primary care physicians), population health measures (i.e., life expectancy,
obesity), and population density. We are not able to include the results of
the Bayesian ZINB by Neelon (2019) for the COVID analysis because this
Bayesian ZINB algorithm cracks after few iterations due to a singularity issue
in high-dimensional spatial/temporal covariance matrix inversions, which also
happened in simulation 2.

We consider the daily COVID-19 deaths in Florida at the county level reported
from 3/25/2020 through 7/29/2020. More specifically, because of the variation
in the schedule on which deaths are reported by days of week, we analyze the
7-day average deaths from the New York Times GitHub repository
(https://github.com/nytimes/covid-19-data). The final analytic sample
comprises 8,351 observations (67 counties across 127 study days). The county
level deaths in Florida during the selected time period are shown in the left
plot of Figure 11. There were approximately 74% counties with zero deaths
reported in total, 14% counties with an average of one death per day and 12%
counties with 2 ore more deaths per day.

Social vulnerability is measured by the social vulnerability index (SVI),
which is a national and state-specific county ranking system and designed to
assist public health officials in identifying communities in need of support
and resources in an event of pandemic (Centers for Disease Control and
Prevention, 2018). The SVI has four different themes, namely, socioeconomic
status, household composition, race/ethnicity/language, and
housing/transportation, which comprise the overall SVI. We use the overall SVI
in the analysis and categorize it to three levels: high ($>75$ percentile),
moderate (between 25 percentile and 75 percentile), and low ($<25$
percentile). A higher level of overall SVI indicates more vulnerability, as
suggested in the literature (Hughes et al., 2021, Bruckhaus et al., 2022). We
also adjust the following county-level characteristics in the model: other
sociodemographic characteristics (i.e., health insurance coverage,
urbanicity); population health care resources (i.e., primary care physicians);
population health measures (i.e., life expectancy, obesity); and population
density. The variables used in the analysis and the data sources are listed in
Table 5 (Dasgupta et al., 2020).

A map of the overall SVI among Florida counties (panel (a)) and a time series
plot of the reported daily average deaths for each SVI category (panel (b))
are shown in Figure 7. In panel (a), it is shown that most counties have
moderate overall SVI (green color in the map). In panel (b), the counties with
moderate overall SVI have more deaths than the high and low group. It is
different from the finding that counties with high SVI had more COVID-19
deaths in other research focused on the U.S. or other locations (Karmakar et
al., 2021, Neelon et al., 2022). To further understand how the COVID-19 deaths
spread across the three SVI categories, we show the map of the daily average
deaths in four different months in Figure 8. Counties at a higher risk such as
Miami-Dade and Hillsborough have more COVID-19 deaths and they have moderate
overall SVI scores.

Variable |  Year |  Source
---|---|---
COVID-19 deaths |  2020 |  The New York Times
Social Vulnerability Index and component measures |  |
Overall Social Vulnerability Index |  2014-2018 |  Centers for Disease Control and Prevention
Socioeconomic status index |  2014-2018 |  Centers for Disease Control and Prevention
Poverty rate |  2014-2018 |  Centers for Disease Control and Prevention
Household characteristics and disability index |  2014-2018 |  Centers for Disease Control and Prevention
65 years or older |  2014-2018 |  Centers for Disease Control and Prevention
Minority status and language index |  2014-2018 |  Centers for Disease Control and Prevention
Minority |  2014-2018 |  Centers for Disease Control and Prevention
Housing type and transportation index |  2014-2018 |  Centers for Disease Control and Prevention
No health insurance coverage |  2015-2019 |  US Census Bureau
Public transportation to commute to work |  2015-2019 |  US Census Bureau
Urbanicity |  2013 |  US Department of Agriculture
Total primary care physicians |  2019 |  Health Resources & Services Administration
Life expectancy |  2016-2018 |  Robert Wood Johnson Foundation
Obesity |  2016 |  Robert Wood Johnson Foundation
Population Density |  2014-2018 |  Centers for Disease Control and Prevention
Table 5: Data sources.

(a)  (b)

Figure 7: (a) The map of Florida with colors according to the county’s overall
SVI category. (b) The daily average deaths by overall SVI in Florida from
3/25/2020 to 7/29/2020. Figure 8: The true daily average deaths by month in
the study period.

Table 6 shows the estimated coefficients as well as the 95% credible intervals
in the binary and count component. Counties with higher overall SVI levels are
more likely to report deaths (binary component). This is consistent with the
previous research (Karmakar et al., 2021). Meanwhile, the percentage of people
under age 65 without health insurance is associated with a higher risk of
reporting deaths, which is consistent with finding from other studies
(Mountantonakis et al., 2021, DuPre et al., 2021). It is also found that
obesity is associated with a lower risk of reporting deaths, although studies
have suggested that obesity is a risk factor for death from COVID-19 (Tartof
et al., 2020, Rottoli et al., 2020). This may be due to potentially
confounding effects of obesity with the overall SVI (An and Xiang, 2015).
Population density is associated with a higher risk of reporting deaths and
more deaths, which is consistent with other research (Neelon et al., 2021, Foo
et al., 2021).

Variable |  Binary component |  Count component
---|---|---
Overall SVI |  1.06 (0.011, 2.109) |  0.418 (-0.412, 1.229)
No health insurance coverage |  0.308 (0.114, 0.492) |  0.014 (-0.118, 0.148)
Urbanicity |  0.021 (-0.471, 0.484) |  -0.261 (-0.706, 0.153)
Total primary care physicians |  -0.017 (-0.038, 0.004) |  -0.005 (-0.019, 0.009)
Life expectancy |  -0.049 (-0.14, 0.051) |  -0.009 (-0.181, 0.154)
Obesity |  -0.248 (-0.399, -0.104) |  -0.062 (-0.187, 0.062)
Population density |  1.387 (0.658, 2.109) |  1.415 (0.883, 1.993)
Table 6: The parameter estimates and 95% credible intervals of the Florida
counties’ sociodemographic and health characteristics with COVID-19 deaths.

The fitted daily average death by overall SVI category is shown in panel (a)
of Figure 9. The fitted trend (solid line) is consistent with the actual
(dotted line). In the earlier time (March 2020 through the middle of July
2020), there is approximately one death each day for counties with moderate
SVI and zero deaths for counties with low or high SVI. Starting from July
2020, the death counts has increased for all counties in Florida, and counties
in the moderate category had the most significant increase than others. We
also apply a recently developed spatiotemporal model (not in a zero-inflated
setting) (Neelon et al., 2022) which used B-splines to model temporal effects
for the COVID-19 analysis. We notice that it could be difficult to set the
location and number of knots appropriately. Its fitted daily average deaths by
SVI category are shown in panel (b) of Figure 9. There is a higher predictive
uncertainty at the end of the temporal period for all three overall SVI
categories.

(a)  (b)

Figure 9: The true and fitted daily average deaths in Florida by SVI category
using the proposed model in panel (a) and the B-splines spatiotemporal model
proposed by Neelon et al. (2022) in panel (b). The dotted lines are the true
daily average, solid lines are the fitted daily average, and the grey ribbons
are the 95% credible region.

To further compare the death counts across the three SVI categories, risk
ratios (RRs) for the high and moderate SVI versus the low SVI are shown in
Figure 10, the counties in moderate SVI category have a higher risk of
reporting deaths than the counties in the low SVI category (the blue dotted
line). The counties in the high SVI category have a similar risk of reporting
death with counties in the low SVI category between March 2020 and June 2020
and lower risk in July 2020 (red dotted line).

Figure 10: The RRs for counties with high and moderate social vulnerability
versus low social vulnerability. The dotted lines are posterior mean RRs and
the grey ribbons were the 95% credible regions. RRs $>$ 1 indicate higher
risk. The black dashed line is the reference line with risk ratio $=1$.

The county level observed average death counts and posterior mean average
death counts in the overall SVI model is shown in Figure 11. The proposed
model identifies counties with high death counts such as Miami-Dade.

Figure 11: The true and fitted daily average deaths in Florida by counties.

(a)  (b)

(c)  (d)

Figure 12: The true and posterior mean daily deaths in Hernando county (a),
Broward county (b), Polk county (c), and Miami-Dade county (d) from 3/25/2020
to 7/29/2020 using our proposed method. The dotted lines are the true deaths
and the solid lines are the fitted deaths. The grey ribbons are the 95%
credible intervals.

(a)  (b)

(c)  (d)

Figure 13: The true and posterior mean daily deaths in Hernando county (a),
Broward county (b), Polk county (c), and Miami-Dade county (d) from 3/25/2020
to 7/29/2020 using the B-splines spatiotemporal model proposed by Neelon et
al. (2022). The dotted lines are the true deaths and the solid lines are the
fitted deaths. The grey ribbons are the 95% credible intervals.

In Figure 12, we select four counties from the three SVI categories to show
the individual fitted daily deaths: Hernando county (moderate SVI), Broward
county (high SVI), Polk county (low SVI), and Miami-Dade county (moderate
SVI). The fitted curves capture the overall patterns but are generally more
smooth than the observed death curves which are more wiggly. The model does
not capture the downward trend in Miami-Dade (see panel (d) of Figure 12) in
the later period of July 2020. This may be related to emergency orders of
closing certain indoor spaces and outdoor spaces where groups of people
congregated without physical distancing in July 2020, which are not included
in the analysis. For comparisons, we show the fitting results using the
spatiotemporal model proposed in Neelon et al. (2022) in Figure 13. The
posterior mean daily deaths do not capture the trends very well and the 95%
CIs fail to cover most the true values for Hernando county in panel (a), Polk
county in panel (c), and Miami-Dade county in panel (d). The result for
Broward county in panel (b) is similar to that of the proposed model.

##  5 Discussion

We have introduced a Bayesian framework for zero-inflated negative binomial
regression models with spatiotemporal effects. The proposed model is able to
simultaneously incorporate a flexible structure for the spatial effects and
temporal for the binary and count component. The Gibbs sampler is efficient
for posterior inference via the Pólya-Gamma data augmentation and latent NNGP
conditioned on these latent variables through the Bayesian hierarchical
inference. The use of the latent NNGP to approximate the covariance inverse
matrix empowers the proposed to be feasible when the spatial or temporal
dimension is large. Our simulations suggest that the proposed model is
comparable to existing methods when such comparisons are available and more
applicable in a higher dimension. It also outperforms other methods when there
is only one repetition in each sampling unit.

There are several potential improvements for future work. There is potential
for selecting important variables in this Bayesian framework. For example, by
using multivariate Bayesian models with polynomial-tailed shrinkage priors
(Wang et al., 2023) for the covariance matrix, one can derive variable
selection methods. Such methods would be useful to model ultrahigh dimensional
spatiotemporal data. Additionally, the proposed model assumes a common
temporal effect across different spatial locations and a common spatial effect
across different temporal points. Future study could be focused on dynamic
settings that capture the interaction between spatial and temporal effects.

##  6 Software

Software in the form of R code, together with a sample input data set and
complete documentation is available on request from the corresponding author
(hsin.huang@ucf.edu).

## Acknowledgments

The authors thank the reviewers for providing valuable comments. We would like
to acknowledge support for this project from the National Science Foundation
(NSF grant DMS-1924792). Conflict of Interest: None declared.

## Appendix

### The step-by-step posterior sampling procedure

The update algorithm for parameters of interests as well as random effects and
hyperparameters are shown below. In each component, we update parameters and
random effects together as the posterior distribution are both Gaussian, which
can facilitate convergence of the MCMC algorithm in simulation studies.

Step 1: Update the latent at-risk indicators
The update for $W_{j}$ depends on the value of $y_{j}$. If $y_{j}>0$, then the
$j$th subject belongs to the at-risk class and $W_{j}=1$. If $y_{j}=0$,
$W_{j}$ is updated by drawing from a Bernoulli distribution with probability

$\displaystyle Pr(W_{j}=1\mid
y_{j}=0,\text{rest})=\frac{\pi_{j}\nu_{j}^{r}}{1-\pi_{j}(1-\nu_{j}^{r})},$

where $\pi_{j}=\frac{\exp(\eta_{1j})}{1+\exp(\eta_{1j})}$ is the unconditional
probability that $W_{j}=1$ and $\nu_{j}=1-\psi_{j}$.

Step 2: Update $\alpha,\vec{a},\vec{b}$
Denote $[\alpha,\vec{a},\vec{b}]^{T}$ by $\phi_{1}$. Assume a
$\mathcal{N}(\phi_{0},\Sigma_{0})$ prior, we update $\phi_{1}$ from the
posterior distribution in two steps:
$\omega_{1j}\mid\eta_{1j}\sim\mathcal{PG}(1,\eta_{1j})$ and
$\phi_{1}\mid\Omega_{1},z_{1}\sim\mathcal{N}(\mu,\Sigma),$ where

$\displaystyle\Sigma$ $\displaystyle=$
$\displaystyle(\Sigma_{0}^{-1}+V^{T}\Omega_{1}V)^{-1},$ $\displaystyle\mu$
$\displaystyle=$
$\displaystyle\Sigma\left(\Sigma_{0}^{-1}\phi_{0}+V^{T}\Omega_{1}(z_{1}-V_{1}\vec{\epsilon}_{11}-V_{2}\vec{\epsilon}_{12})\right),$

and $z_{1j}=\frac{W_{j}-1/2}{\omega_{1j}}$ is the latent variable,
$\Omega_{1}=\text{diag}(\omega_{1j})$, and $V=[X,V_{1},V_{2}]$ is the
augmented design matrix.

Step 3: Update $\vec{\epsilon}_{11}$
Conditional on $\alpha,\vec{a},\vec{b},\vec{\epsilon}_{12}$, assuming a
$\mathcal{N}(0,\sigma^{2}_{\epsilon_{11}})$ prior, update
$\vec{\epsilon}_{11}$ from a $\mathcal{N}(\mu,\Sigma)$ distribution where

$\displaystyle\Sigma$ $\displaystyle=$
$\displaystyle\left(\frac{1}{\sigma^{2}_{\epsilon_{11}}}I+V_{1}^{T}\Omega_{1}V_{1}\right)^{-1},$
$\displaystyle\mu$ $\displaystyle=$
$\displaystyle\Sigma\left(V_{1}^{T}\Omega_{1}(z_{1}-X\alpha-
V_{1}\vec{a}-V_{2}\vec{b}-V_{2}\vec{\epsilon}_{12})\right).$

Step 4: Update $\vec{\epsilon}_{12}$
Conditional on $\alpha,\vec{a},\vec{b},\vec{\epsilon}_{11}$, assuming a
$\mathcal{N}(0,\sigma^{2}_{\epsilon_{12}})$ prior, update
$\vec{\epsilon}_{12}$ from a $\mathcal{N}(\mu,\Sigma)$ distribution, where

$\displaystyle\Sigma$ $\displaystyle=$
$\displaystyle\left(\frac{1}{\sigma^{2}_{\epsilon_{12}}}I+V_{2}^{T}\Omega_{1}V_{2}\right)^{-1},$
$\displaystyle\mu$ $\displaystyle=$
$\displaystyle\Sigma\left(V_{2}^{T}\Omega_{1}(z_{1}-X\alpha-
V_{1}\vec{a}-V_{2}\vec{b}-V_{1}\vec{\epsilon}_{11})\right).$

Step 5: Update $l_{11},\sigma_{11}$
First, we use the Metropolis-Hasting algorithm to update $l_{11}$. Using a
normal proposal distribution which is symmetric, the acceptance probability
can be written as follows:

$\displaystyle A$
$\displaystyle=\min\left(1,\frac{\pi(l^{*}_{11})}{\pi(l_{11})}\right)$
$\displaystyle=\min\left(1,\frac{p(\vec{a}\mid
K_{\sigma_{11},l^{*}_{11}}(h_{1},\ldots,h_{S}))p(l^{*}_{11}\mid
a_{l_{1}},b_{l_{1}})}{p(\vec{a}\mid
K_{\sigma_{11},l_{11}}(h_{1},\ldots,h_{S}))p(l_{11}\mid
a_{l_{1}},b_{l_{1}})}\right).$

Then, $\sigma_{11}$ can be sampled from the full conditional as

$\sigma^{2}_{11}\sim\mathcal{IG}\left(a_{\sigma_{1}}+\frac{S}{2},b_{\sigma_{1}}+\frac{\vec{a}^{T}(K_{\sigma_{11},l_{11}}(h_{1},\ldots,h_{S}))^{-1}\vec{a}}{2}\right)$

The posterior distribution of $\sigma_{11}$ requires the inverse of
$K_{\sigma_{11},l_{11}}(h_{1},\ldots,h_{S})$, which may be a computation
bottleneck when the size of spatial locations is large. Hence, we apply the
latent conjugate NNGP to approximate the inverse matrix.

Step 6: Update $l_{12},\sigma^{2}_{12}$
The update for $l_{12},\sigma^{2}_{12}$ is similar to the one for
$l_{11},\sigma^{2}_{11}$. Depending on the size of temporal points, the latent
conjugate NNGP can also be applied in the update for $\sigma^{2}_{12}$ if $T$
is large. First, we use the Metropolis-Hasting algorithm to update $l_{12}$.
Using a normal proposal distribution which is symmetric, the acceptance
probability can be written as follows:

$\displaystyle A$
$\displaystyle=\min\left(1,\frac{\pi(l^{*}_{12})}{\pi(l_{12})}\right)$
$\displaystyle=\min\left(1,\frac{p(\vec{b}\mid
K_{\sigma_{12},l^{*}_{12}}(w_{1},\ldots,w_{T}))p(l^{*}_{12}\mid
a_{l_{2}},b_{l_{2}})}{p(\vec{b}\mid
K_{\sigma_{12},l_{12}}(w_{1},\ldots,w_{T}))p(l_{12}\mid
b_{l_{2}},b_{l_{2}})}\right).$

Then, $\sigma_{12}$ can be sampled from the full conditional as

$\sigma^{2}_{12}\sim\mathcal{IG}\left(a_{\sigma_{1}}+\frac{T}{2},b_{\sigma_{2}}+\frac{\vec{b}^{T}(K_{\sigma_{12},l_{12}}(w_{1},\ldots,w_{T}))^{-1}\vec{b}}{2}\right).$

Step 7: Update $\sigma_{\epsilon_{11}}$
Assuming a $\mathcal{IG}(a_{\epsilon},b_{\epsilon})$ prior, draw
$\sigma_{\epsilon_{11}}$ from the posterior distribution:

$\displaystyle\sigma^{2}_{\epsilon_{11}}\sim\mathcal{IG}\left(a_{\epsilon}+\frac{S}{2},b_{\epsilon}+\frac{\sum_{s=1}^{S}\epsilon^{2}_{11s}}{2}\right).$

Step 8: Update $\sigma_{\epsilon_{12}}$
Assuming a $\mathcal{IG}(a_{\epsilon},b_{\epsilon})$ prior, draw
$\sigma_{\epsilon_{12}}$ from the posterior distribution:

$\displaystyle\sigma^{2}_{\epsilon_{12}}\sim\mathcal{IG}\left(a_{\epsilon}+\frac{T}{2},b_{\epsilon}+\frac{\sum_{t=1}^{T}\epsilon^{2}_{12t}}{2}\right).$

Step 9: Update $\beta,\vec{c},\vec{d}$
Denote $[\beta,\vec{c},\vec{d}]^{T}$ by $\phi_{2}$. Assume a
$\mathcal{N}(\phi_{0},\Sigma_{0})$ prior, we update $\phi_{2}$ in two steps:

$\displaystyle\omega_{2j}|\eta_{2j}$ $\displaystyle\sim$
$\displaystyle\mathcal{PG}(y_{j}+r,\eta_{2j}),\forall j\text{ s.t. }W_{j}=1$
$\displaystyle\phi_{2}|\Omega_{2},z_{2}$ $\displaystyle\sim$
$\displaystyle\mathcal{N}(\mu,\Sigma),$

where

$\displaystyle\Sigma$ $\displaystyle=$
$\displaystyle(\Sigma_{0}^{-1}+V^{*T}\Omega_{2}V^{*})^{-1}$ $\displaystyle\mu$
$\displaystyle=$
$\displaystyle\Sigma\left(\Sigma_{0}^{-1}\phi_{0}+V^{*T}\Omega_{2}(z_{2}-V^{*}_{1}\vec{\epsilon}_{21}-V^{*}_{2}\vec{\epsilon}_{22})\right),$

where $z_{2}$ is the latent variable with
$z_{2j}=\frac{Y_{j}-r}{2\omega_{2j}}$, $\Omega_{2}=\text{diag}(\omega_{2j})$,
and $V^{*}=[X^{*},V^{*}_{1},V^{*}_{2}]$ is the augmented design matrix for
observations at risk ($W_{j}=1$).

Step 10: Update $\vec{\epsilon}_{21}$
Conditional on $\beta,\vec{c},\vec{d},\vec{\epsilon}_{22}$, assuming a
$\mathcal{N}(0,\sigma^{2}_{\epsilon_{21}})$ prior, update
$\vec{\epsilon}_{21}$ from a $\mathcal{N}(\mu,\Sigma)$ distribution where

$\displaystyle\Sigma$ $\displaystyle=$
$\displaystyle\left(\frac{1}{\sigma^{2}_{\epsilon_{2S}}}I+(V^{*}_{1})^{T}\Omega_{2}V^{*}_{1}\right)^{-1},$
(8) $\displaystyle\mu$ $\displaystyle=$
$\displaystyle\Sigma\left((V^{*}_{1})^{T}\Omega_{2}(z_{2}-X\beta-V^{*}_{1}\vec{c}-V^{*}_{2}\vec{d}-V^{*}_{2}\vec{\epsilon}_{22})\right).$
(9)

Step 11: Update $\vec{\epsilon}_{22}$
Conditional on $\beta,\vec{c},\vec{d},\vec{\epsilon}_{21}$, assuming a
$\mathcal{N}(0,\sigma^{2}_{\epsilon_{22}})$ prior, update
$\vec{\epsilon}_{22}$ from a $\mathcal{N}(\mu,\Sigma)$ distribution, where

$\displaystyle\Sigma$ $\displaystyle=$
$\displaystyle\left(\frac{1}{\sigma^{2}_{\epsilon_{22}}}I+(V^{*}_{2})^{T}\Omega_{2}V^{*}_{2}\right)^{-1},$
(10) $\displaystyle\mu$ $\displaystyle=$
$\displaystyle\Sigma\left((V^{*}_{2})^{T}\Omega_{2}(z_{2}-X\beta-V^{*}_{1}\vec{c}-V^{*}_{2}\vec{d}-V^{*}_{1}\vec{\epsilon}_{21})\right).$
(11)

Step 12: Update $l_{21},\sigma_{21}$
The update of $l_{21}$ is similar to that of $l_{11}$. The Metropolis-Hasting
algorithm is used to update $l_{21}$ with acceptance probability

$\displaystyle A$
$\displaystyle=\min\left(1,\frac{\pi(l^{*}_{21})}{\pi(l_{21})}\right)$
$\displaystyle=\min\left(1,\frac{p(\vec{c}\mid
K_{\sigma_{21},l^{*}_{21}}(h_{1},\ldots,h_{S}))p(l^{*}_{21}\mid
a_{l_{1}},b_{l_{1}})}{p(\vec{c}\mid
K_{\sigma_{21},l_{21}}(h_{1},\ldots,h_{S}))p(l_{21}\mid
a_{l_{1}},b_{l_{1}})}\right).$

Hyperparameter $\sigma_{21}$ is drawn from the full conditional distribution

$\sigma^{2}_{21}\sim\mathcal{IG}\left(a_{\sigma_{1}}+\frac{S}{2},b_{\sigma_{1}}+\frac{\vec{c}^{T}(K_{\sigma_{21},l_{21}}(h_{1},\ldots,h_{S}))^{-1}\vec{c}}{2}\right),$

where $K^{-1}_{\sigma_{21},l_{21}}$ is approximated using latent NNGP.

Step 13: Update $l_{22},\sigma^{2}_{22}$
First, $l_{22}$ is updated using the Metropolis-Hasting algorithm with
acceptance probability

$\displaystyle A$
$\displaystyle=\min\left(1,\frac{\pi(l^{*}_{22})}{\pi(l_{22})}\right)$
$\displaystyle=\min\left(1,\frac{p(\vec{d}\mid
K_{\sigma_{22},l^{*}_{22}}(w_{1},\ldots,w_{T}))p(l^{*}_{22}\mid
a_{l_{2}},b_{l_{2}})}{p(\vec{d}\mid
K_{\sigma_{22},l_{22}}(w_{1},\ldots,w_{T}))p(l_{22}\mid
b_{l_{2}},b_{l_{2}})}\right).$

Then, $\sigma_{22}$ can be sampled from the full conditional as

$\sigma^{2}_{22}\sim\mathcal{IG}\left(a_{\sigma_{2}}+\frac{T}{2},b_{\sigma_{2}}+\frac{\vec{d}^{T}(K_{\sigma_{22},l_{22}}(w_{1},\ldots,w_{T}))^{-1}\vec{d}}{2}\right).$

Step 14: Update $\sigma_{\epsilon_{21}}$
Assuming a $\mathcal{IG}(a_{\epsilon},b_{\epsilon})$ prior, draw
$\sigma_{\epsilon_{21}}$ from the posterior distribution:

$\displaystyle\sigma^{2}_{\epsilon_{21}}\sim\mathcal{IG}\left(a_{\epsilon}+\frac{S}{2},b_{\epsilon}+\frac{\sum_{s=1}^{S}\epsilon^{2}_{21s}}{2}\right).$

Step 15: Update $\sigma_{\epsilon_{22}}$
Assuming a $\mathcal{IG}(a_{\epsilon},b_{\epsilon})$ prior, draw
$\sigma_{\epsilon_{22}}$ from the posterior distribution:

$\displaystyle\sigma^{2}_{\epsilon_{22}}\sim\mathcal{IG}\left(a_{\epsilon}+\frac{T}{2},b_{\epsilon}+\frac{\sum_{t=1}^{T}\epsilon^{2}_{22t}}{2}\right).$

## References

  * An and Xiang (2015) R. An and X. Xiang.  Social vulnerability and obesity among us adults.  _International Journal of Health Sciences_ , 3(3):7–21, 2015.
  * Banerjee et al. (2003) S. Banerjee, B. P. Carlin, and A. E. Gelfand.  _Hierarchical modeling and analysis for spatial data_.  Chapman and Hall/CRC, 2003.
  * Bruckhaus et al. (2022) A. A. Bruckhaus, A. Abedi, S. Salehi, T. A. Pickering, Y. Zhang, A. Martinez, M. Lai, R. Garner, and D. Duncan.  Covid-19 vaccination dynamics in the us: coverage velocity and carrying capacity based on socio-demographic vulnerability indices in california.  _Journal of immigrant and minority health_ , 24(1):18–30, 2022.
  * Casella and George (1992) G. Casella and E. I. George.  Explaining the gibbs sampler.  _The American Statistician_ , 46(3):167–174, 1992\.
  * Centers for Disease Control and Prevention (2018) Centers for Disease Control and Prevention.  CDC/ATSDR Social Vulnerability Index 2018 Database Florida.  https://www.atsdr.cdc.gov/placeandhealth/svi/data_documentation_download.html, 2018\.  Accessed: 2022-02-01.
  * Chalupka et al. (2013) K. Chalupka, C. K. Williams, and I. Murray.  A framework for evaluating approximation methods for gaussian process regression.  _Journal of Machine Learning Research_ , 14:333–350, 2013\.
  * Cheung (2002) Y. B. Cheung.  Zero-inflated models for regression analysis of count data: a study of growth and development.  _Statistics in medicine_ , 21(10):1461–1469, 2002\.
  * Cousin et al. (2016) A. Cousin, H. Maatouk, and D. Rullière.  Kriging of financial term-structures.  _European Journal of Operational Research_ , 255(2):631–648, 2016.
  * Cressie and Wikle (2015) N. Cressie and C. K. Wikle.  _Statistics for spatio-temporal data_.  John Wiley & Sons, 2015.
  * Dadaneh et al. (2018) S. Z. Dadaneh, M. Zhou, and X. Qian.  Bayesian negative binomial regression for differential expression with confounding factors.  _Bioinformatics_ , 34(19):3349–3356, 2018.
  * Dasgupta et al. (2020) S. Dasgupta, V. B. Bowen, A. Leidner, K. Fletcher, T. Musial, C. Rose, A. Cha, G. Kang, E. Dirlikov, E. Pevzner, et al.  Association between social vulnerability and a county’s risk for becoming a covid-19 hotspot—united states, june 1–july 25, 2020.  _Morbidity and Mortality Weekly Report_ , 69(42):1535, 2020.
  * Datta et al. (2016a) A. Datta, S. Banerjee, A. O. Finley, and A. E. Gelfand.  Hierarchical nearest-neighbor gaussian process models for large geostatistical datasets.  _Journal of the American Statistical Association_ , 111(514):800–812, 2016a.
  * Datta et al. (2016b) A. Datta, S. Banerjee, A. O. Finley, N. A. Hamm, and M. Schaap.  Nonseparable dynamic nearest neighbor gaussian process models for large spatio-temporal data with an application to particulate matter analysis.  _The annals of applied statistics_ , 10(3):1286, 2016b.
  * Diana et al. (2021) A. Diana, E. Dennis, E. Matechou, and B. Morgan.  Fast bayesian inference for large occupancy data sets, using the polya-gamma scheme.  _arXiv preprint arXiv:2107.14656_ , 2021.
  * DuPre et al. (2021) N. C. DuPre, S. Karimi, C. H. Zhang, L. Blair, A. Gupta, L. M. A. Alharbi, M. Alluhibi, R. Mitra, W. P. McKinney, and B. Little.  County-level demographic, social, economic, and lifestyle correlates of covid-19 infection and death trajectories during the first wave of the pandemic in the united states.  _Science of The Total Environment_ , 786:147495, 2021.
  * Foo et al. (2021) O. Foo, S. Hiu, D. Teare, A. A. Syed, and S. Razvi.  A global country-level analysis of the relationship between obesity and covid-19 cases and mortality.  _Diabetes, Obesity and Metabolism_ , 23(12):2697–2706, 2021.
  * Gelman et al. (1995) A. Gelman, J. B. Carlin, H. S. Stern, and D. B. Rubin.  _Bayesian data analysis_.  Chapman and Hall/CRC, 1995.
  * Ghosh et al. (2006) S. K. Ghosh, P. Mukhopadhyay, and J.-C. J. Lu.  Bayesian analysis of zero-inflated regression models.  _Journal of Statistical planning and Inference_ , 136(4):1360–1375, 2006.
  * Gramacy (2016) R. B. Gramacy.  lagp: large-scale spatial modeling via local approximate gaussian processes in r.  _Journal of Statistical Software_ , 72:1–46, 2016.
  * Greene (1994) W. H. Greene.  Accounting for excess zeros and sample selection in poisson and negative binomial regression models.  1994\.
  * Gu et al. (2020) X. Gu, X. Yan, L. Ma, and X. Liu.  Modeling the service-route-based crash frequency by a spatiotemporal-random-effect zero-inflated negative binomial model: An empirical analysis for bus-involved crashes.  _Accident Analysis & Prevention_, 144:105674, 2020.
  * Hughes et al. (2021) M. M. Hughes, A. Wang, M. K. Grossman, E. Pun, A. Whiteman, L. Deng, E. Hallisey, J. D. Sharpe, E. N. Ussery, S. Stokley, et al.  County-level covid-19 vaccination coverage and social vulnerability—united states, december 14, 2020–march 1, 2021.  _Morbidity and Mortality Weekly Report_ , 70(12):431, 2021.
  * Jackman et al. (2015) S. Jackman, A. Tahk, A. Zeileis, C. Maimone, J. Fearon, Z. Meers, M. S. Jackman, and M. Imports.  Package ‘pscl’.  _Political Science Computational Laboratory_ , 18(04.2017), 2015.
  * Kang and Cressie (2011) E. L. Kang and N. Cressie.  Bayesian inference for the spatial random effects model.  _Journal of the American Statistical Association_ , 106(495):972–983, 2011.
  * Karmakar et al. (2021) M. Karmakar, P. M. Lantz, and R. Tipirneni.  Association of social and demographic factors with covid-19 incidence and death rates in the us.  _JAMA network open_ , 4(1):e2036462–e2036462, 2021.
  * Khedhiri (2021) S. Khedhiri.  Statistical modeling of covid-19 deaths with excess zero counts.  _Epidemiologic Methods_ , 10(s1), 2021.
  * Kim et al. (2018) S. Kim, Z. Lu, and A. S. Cohen.  An improved estimation using polya-gamma augmentation for bayesian structural equation models with dichotomous variables.  _Measurement: Interdisciplinary Research and Perspectives_ , 16(2):81–91, 2018.
  * Kimeldorf and Wahba (1970) G. S. Kimeldorf and G. Wahba.  A correspondence between bayesian estimation on stochastic processes and smoothing by splines.  _The Annals of Mathematical Statistics_ , 41(2):495–502, 1970.
  * Lázaro-Gredilla et al. (2010) M. Lázaro-Gredilla, J. Quinonero-Candela, C. E. Rasmussen, and A. R. Figueiras-Vidal.  Sparse spectrum gaussian process regression.  _The Journal of Machine Learning Research_ , 11:1865–1881, 2010.
  * Lewsey and Thomson (2004) J. D. Lewsey and W. M. Thomson.  The utility of the zero-inflated poisson and zero-inflated negative binomial models: a case study of cross-sectional and longitudinal dmf data examining the effect of socio-economic status.  _Community dentistry and oral epidemiology_ , 32(3):183–189, 2004.
  * Liu et al. (2020) H. Liu, Y.-S. Ong, X. Shen, and J. Cai.  When gaussian process meets big data: A review of scalable gps.  _IEEE transactions on neural networks and learning systems_ , 31(11):4405–4423, 2020.
  * Lunn et al. (2013) D. Lunn, C. Jackson, N. Best, A. Thomas, and D. Spiegelhalter.  The bugs book.  _A Practical Introduction to Bayesian Analysis, Chapman Hall, London_ , 2013.
  * McCullagh and Nelder (2019) P. McCullagh and J. A. Nelder.  _Generalized linear models_.  Routledge, 2019.
  * Mountantonakis et al. (2021) S. E. Mountantonakis, L. M. Epstein, K. Coleman, J. Martinez, M. Saleh, C. Kvasnovsky, R.-M. Brown, E. McCulloch, J. Kuvin, S. Richardson, et al.  The association of structural inequities and race with out-of-hospital sudden death during the covid-19 pandemic.  _Circulation: Arrhythmia and Electrophysiology_ , 14(5):e009646, 2021.
  * Neelon (2019) B. Neelon.  Bayesian zero-inflated negative binomial regression based on pólya-gamma mixtures.  _Bayesian Analysis_ , 14(3):829–855, 2019.
  * Neelon et al. (2021) B. Neelon, F. Mutiso, N. T. Mueller, J. L. Pearce, and S. E. Benjamin-Neelon.  Spatial and temporal trends in social vulnerability and covid-19 incidence and death rates in the united states.  _Plos one_ , 16(3):e0248702, 2021.
  * Neelon et al. (2022) B. Neelon, C.-C. Wen, and S. E. Benjamin-Neelon.  A multivariate spatiotemporal model for tracking covid-19 incidence and death rates in socially vulnerable populations.  _Journal of Applied Statistics_ , pages 1–24, 2022.
  * Neelon et al. (2010) B. H. Neelon, A. J. O’Malley, and S.-L. T. Normand.  A bayesian model for repeated measures zero-inflated count data with application to outpatient psychiatric service use.  _Statistical modelling_ , 10(4):421–439, 2010\.
  * Pillow and Scott (2012) J. Pillow and J. Scott.  Fully bayesian inference for neural models with negative-binomial spiking.  _Advances in neural information processing systems_ , 25:1898–1906, 2012.
  * Polson et al. (2013) N. G. Polson, J. G. Scott, and J. Windle.  Bayesian inference for logistic models using pólya–gamma latent variables.  _Journal of the American statistical Association_ , 108(504):1339–1349, 2013.
  * Quinonero-Candela and Rasmussen (2005) J. Quinonero-Candela and C. E. Rasmussen.  A unifying view of sparse approximate gaussian process regression.  _The Journal of Machine Learning Research_ , 6:1939–1959, 2005.
  * R Core Team (2021) R Core Team.  _R: A Language and Environment for Statistical Computing_.  R Foundation for Statistical Computing, Vienna, Austria, 2021.  URL https://www.R-project.org/.
  * Rottoli et al. (2020) M. Rottoli, P. Bernante, A. Belvedere, F. Balsamo, S. Garelli, M. Giannella, A. Cascavilla, S. Tedeschi, S. Ianniruberto, E. R. Del Turco, et al.  How important is obesity as a risk factor for respiratory failure, intensive care admission and death in hospitalised covid-19 patients? results from a single italian centre.  _European Journal of Endocrinology_ , 183(4):389–397, 2020.
  * Smola and Bartlett (2001) A. J. Smola and P. L. Bartlett.  Sparse greedy gaussian process regression.  In _Advances in neural information processing systems_ , pages 619–625, 2001.
  * Snelson and Ghahramani (2007) E. Snelson and Z. Ghahramani.  Local and global sparse gaussian process approximations.  In _Artificial Intelligence and Statistics_ , pages 524–531. PMLR, 2007.
  * Stan Development Team (2021) Stan Development Team.  RStan: the R interface to Stan, 2021.  URL https://mc-stan.org/.  R package version 2.21.3.
  * Tartof et al. (2020) S. Y. Tartof, L. Qian, V. Hong, R. Wei, R. F. Nadjafi, H. Fischer, Z. Li, S. F. Shaw, S. L. Caparosa, C. L. Nau, et al.  Obesity and mortality among patients diagnosed with covid-19: results from an integrated health care organization.  _Annals of internal medicine_ , 173(10):773–781, 2020.
  * Titsias (2009) M. Titsias.  Variational learning of inducing variables in sparse gaussian processes.  In _Artificial intelligence and statistics_ , pages 567–574. PMLR, 2009.
  * Vasudevan et al. (2009) S. Vasudevan, F. Ramos, E. Nettleton, and H. Durrant-Whyte.  Gaussian process modeling of large-scale terrain.  _Journal of Field Robotics_ , 26(10):812–840, 2009.
  * Wahba (1990) G. Wahba.  _Spline models for observational data_.  SIAM, 1990.
  * Wang et al. (2023) S.-H. Wang, R. Bai, and H.-H. Huang.  Mixed-type multivariate bayesian sparse variable selection with shrinkage priors, 2023.
  * Wang et al. (2016) X. Wang, M.-H. Chen, R. C. Kuo, and D. K. Dey.  Dynamic spatial pattern recognition in count data.  In _New Developments in Statistical Modeling, Inference and Application_ , pages 185–202. Springer, 2016.
  * Yau et al. (2003) K. K. Yau, K. Wang, and A. H. Lee.  Zero-inflated negative binomial mixed regression modeling of over-dispersed count data with extra zeros.  _Biometrical Journal: journal of mathematical methods in biosciences_ , 45(4):437–452, 2003.
  * Zhang et al. (2021a) H.-Y. Zhang, A.-R. Zhang, Q.-B. Lu, X.-A. Zhang, Z.-J. Zhang, X.-G. Guan, T.-L. Che, Y. Yang, H. Li, W. Liu, et al.  Association between fatality rate of covid-19 and selenium deficiency in china.  _BMC Infectious Diseases_ , 21(1):1–8, 2021a.
  * Zhang et al. (2021b) L. Zhang, S. Banerjee, and A. O. Finley.  High-dimensional multivariate geostatistics: A bayesian matrix-normal approach.  _Environmetrics_ , page e2675, 2021b.
  * Zhang et al. (2022) Y. Zhang, H. H. Chang, A. D. Iuliano, and C. Reed.  Application of bayesian spatial-temporal models for estimating unrecognized covid-19 deaths in the united states.  _Spatial statistics_ , page 100584, 2022.
  * Zhou and Carin (2013) M. Zhou and L. Carin.  Negative binomial process count and mixture modeling.  _IEEE Transactions on Pattern Analysis and Machine Intelligence_ , 37(2):307–320, 2013.
  * Zuur et al. (2012) A. F. Zuur et al.  _Zero inflated models and generalized linear mixed models with R_.  Number 574.50182 Z8. 2012.