# Detecting Outbreaks Using a Latent Field: Part II - Scalable Estimation

Wyatt Bridgman
Sandia National Laboratories
Livermore, CA
<EMAIL_ADDRESS>
&Cosmin Safta
Sandia National Laboratories
Livermore, CA
<EMAIL_ADDRESS>
&Jaideep Ray
Sandia National Laboratories
Livermore, CA
<EMAIL_ADDRESS>

###### Abstract

In this paper, we explore whether the infection-rate of a disease can serve as
a robust monitoring variable in epidemiological surveillance algorithms. The
infection-rate is dependent on population mixing patterns that do not vary
erratically day-to-day; in contrast, daily case-counts used in contemporary
surveillance algorithms are corrupted by reporting errors. The technical
challenge lies in estimating the latent infection-rate from case-counts. Here
we devise a Bayesian method to estimate the infection-rate across multiple
adjoining areal units, and then use it, via an anomaly detector, to discern a
change in epidemiological dynamics. We extend an existing model for estimating
the infection-rate in an areal unit by incorporating a Markov random field
model, so that we may estimate infection-rates across multiple areal units,
while preserving spatial correlations observed in the epidemiological
dynamics. To carry out the high-dimensional Bayesian inverse problem, we
develop an implementation of mean-field variational inference specific to the
infection model and integrate it with the random field model to incorporate
correlations across counties. The method is tested on estimating the COVID-19
infection-rates across all 33 counties in New Mexico using data from the
summer of 2020, and then employing them to detect the arrival of the Fall 2020
COVID-19 wave. We perform the detection using a temporal algorithm that is
applied county-by-county. We also show how the infection-rate field can be
used to cluster counties with similar epidemiological dynamics.

_K_ eywords VI, variational inference; ELBO, evidence lower bound; COVID-19;
New Mexico; outbreak detection

##  1 Introduction

There have been many attempts to estimate the infection-rate of an outbreak
[1, 2, 3], especially for the COVID-19 pandemic, including our own work [4, 5,
6]. In these studies, the infection-rate is modeled as a time-varying
function, which is estimated by fitting a disease model to observed data e.g.,
case-counts of symptomatic (diagnosed) individuals per day. The infection-rate
is used to forecast the outbreak over a small time horizon, e.g., two weeks,
and then compared with the data from that period. This ability to compare
observed case-counts with model forecasts can be used to fashion an outbreak
detector for disease surveillance - if the model forecasts and reported case-
counts do not match, it could indicate a change in the epidemiological
dynamics, either due to a change in human mixing patters e.g., due to lock-
downs [6] or due to a new variant of the pathogen [4]. A shortcoming of these
studies is that they do not contain any spatial information on the variation
of the infection-rate, which impairs their usefulness in disease surveillance
- public health policy is determined globally (e.g. for a nation) and then
adapted and applied locally (e.g., in a county). Being able to estimate an
infection-rate _field_ , defined over multiple areal units (e.g., counties),
can thus be very helpful. However, the task is challenging, as the problem now
requires one to use epidemiological data e.g., case-counts, collected from
each areal unit, which could have a small population. Such data, gathered from
small populations, tends to be contaminated with high-variance noise
(reporting errors) and has to be compensated for by imposing the spatial
patterns extant in epidemiological dynamics caused by mixing of humans between
neighboring areal units. In the Part I of this paper [7, 8] (henceforth _Part
I paper_), we formulated and solved an inverse problem for estimating the
infection-rate field and demonstrated it on COVID-19 case-count data collected
between June 1, 2020 and September 15, 2020 from the counties of New Mexico,
USA. We used the infection-rate to devise an outbreak detector to detect the
start of the Fall 2020 wave of COVID-19 infections in New Mexico (NM), and
determined that it did so about a week earlier than a conventional outbreak
detector. The method employs a separate infection-rate parameterization inside
each areal unit, but links them together via a Gaussian Markov random field
(GMRF) to impose spatial correlations. The parameters are estimated using an
adaptive Markov Chain Monte Carlo (AMCMC [9]) sampler. Since AMCMC is not very
scalable in terms of the number of parameters that can be estimated, we were
limited to a study over three adjoining NM counties - Bernalillo, Santa Fe and
Valencia.

In this paper, we adopt the formulation of the inverse problem from our
previous paper [7, 8] and adapt it for high-dimensional inversion. Our
approach is based on the hypothesis that scalable mean-field variational
inference (MFVI [10]) will allow us to obtain a useful estimate of the
infection-rate field, defined over many areal units, despite the
approximations inherent in MFVI. In MFVI, one imputes a parameterized
posterior distribution for the infection-rate field; in our case, we model the
posterior distribution of (transformations of) our parameters as independent
Gaussians of unknown means and standard deviations. These means and standard
deviations are estimated by minimizing an objective function using an
iterative, gradient-based algorithm. We will test our method by reconstructing
the time-dependent infection-rate field over all 33 counties of NM, with case-
count data of diverse qualities, collected between June 1 and September 15,
2020. We will also fashion an outbreak detector using the infection-rate
fields and check its performance in detecting the Fall 2020 COVID-19 wave
against the one devised in our Part I paper [7, 8] that uses AMCMC to perform
the estimation. In doing so, we will address the following research questions:

  * •

How does the approximate posterior distribution employed in MFVI compare
against the “true” one (computed using AMCMC) from our Part I paper? By
necessity, this comparison will be limited to the three counties included on
our previous paper.

  * •

The spatial correlations embedded in the field estimation problem allow an
arial unit to “borrow information” from its neighbors and compensate for poor-
quality data. How robust is this correlation i.e., when it fails due to poor
quality data, does it do so catastrophically (i.e., the estimation process
stops) or gracefully, with a non-informative (or erroneous) estimate?

  * •

How does the approximate infection-rate field estimated by MFVI affect the
accuracy of the outbreak detector? Is MFVI sufficiently accurate to be useful,
despite its inherent approximations? Further, does the infection-rate field,
estimated over all NM counties, reveal anomalous spatial structures during the
start of the Fall 2020 wave?

The paper has three main contributions. Our first contribution is the
determination of the degree to which the predictive skill of the disease model
forecast is affected by the use of an imputed and approximate posterior
distribution; we find that it is not much vis-à-vis the AMCMC solution. Our
second contribution is the discovery that despite the approximations inherent
in MFVI, the infection-rate field so estimated retains sufficient information
on the outbreak to detect epidemiological anomalies. Note that we will not
attempt to make a proper outbreak detector in this paper; that is left to
future work. Our third contribution is of a numerical nature. Our disease
model includes a convolutional integral making it difficult (and expensive) to
compute the gradient of the objective function, analytically or by finite
differences. Our innovations lie in how the convolution and the gradient are
computed numerically via quadrature and the reparametrization of the objective
function that allows us to use an unconstrained optimization method, despite
the need to preserve non-negativity of parameters being estimated.

The paper is structured as follows. In § 2, we review relevant literature. In
§ 3 we formulate the inverse problem and the MFVI adaptation, and results are
presented in § 4. Outbreak detection is explored § 5. We conclude in § 6.

##  2 Literature Review

In this section, we review some of the literature that undergird our
spatiotemporal epidemiological model, as well as work on variation inference
algorithms that will be used to scale our model to high-dimensional inverse
problems involved in estimating infection-rate fields.

Spatial modeling in epidemiology:  Epidemiological dynamics show spatial
autocorrelation because of human mixing as well as the dependence of some
outbreaks on socioeconomic and demographic covariates which do not change
erratically in space or time. The COVID-19 pandemic was observed and recorded
with fine spatiotemporal granularity, and the data has been subjected to much
spatiotemporal modeling. Many such studies found that the spread of the
disease was mediated mostly by human mixing, rather than by socioeconomic
factors; Huang et al. [11] found it to be so in Hubei province in China, and
we found much the same result for New Mexico in our Part I paper [7, 8]. Geng
et al. [12] analysed US data and found that spatial patterns’ lengthscales
ranged from the county-level to the nation; similar results were found by
Schuler et al. [13] for Germany. McMahon et al. [14] analyzed data at the
county-level and found that the correlation lengthscales for spatial
variability changed during the course of the pandemic; further, the
correlation in epidemiological dynamics between urban centers were stronger
than elsewhere. Indika et al. [15] analyzed data from the counties of Virginia
and found that spatial autocorrelation of case-counts, as quantified by Moran
test statistics, were impacted by, and linked to, executive orders at the
state level. Thus, COVID-19 has presented us with much evidence of spatial
auto-correlation in epidemiological dynamics.

The incorporation of spatial autocorrelation in the modeling and estimation of
infection-rate fields is rare; however, it has been extensively used in
disease mapping. In disease maps, one develops a field, called _relative risk_
$r_{k}$, that is used to adjust an expected value of morbidity $e_{k}$ to the
locally observed case-counts $y^{(obs)}_{k}$ in an areal unit (e.g., county)
$k$, often via a Poisson or Negative Binomial link i.e.,
$y^{(obs)}_{k}\sim{\mathrm{Poisson}}(r_{k}e_{k})$. The expected value $e_{k}$
is usually obtained from a regional average. The relative risk $r_{k}$ is
modeled using covariates of disease activity e.g., socioeconomic conditions
$\log(r_{k})=\boldsymbol{z}_{k}\cdot\boldsymbol{\beta}+\phi_{k}$, where
$\boldsymbol{z}_{k}$ are co-variate risk factors for areal unit $k$,
$\boldsymbol{\beta}$ are regression weights and $\phi_{k}$ captures auto-
correlated random effects in space using a random field model. The simplest
random field model is iCAR (intrinsic Conditional AutoRegressive [16]), a
specific type of Gaussian Markov Random Field (GMRF). Thus

$\phi=\\{\phi_{k}\\}\sim{\mathcal{N}}\left(0,\\{\tau^{2}Q\\}^{-1}\right),\mbox{\hskip
28.45274pt}Q={\rm diag}(W\boldsymbol{1})-W,$

where $W$ is the adjacency matrix of the areal units (i.e., $w_{ij}=1$ if
areal units $i$ and $j$ share a boundary). The object of estimation from data
is $\tau^{2}$. The precision matrix $Q$ tends to be sparse. Another common
model is the Besag-York-Mollie (BYM) model [17] which decomposes $\phi_{k}$ as
$\phi=\phi^{1}+\phi^{2},\phi^{1}\sim{\mathcal{N}}(0,\\{\tau^{2}Q\\}^{-1})$ and
$\phi^{2}\sim{\mathcal{N}}(0,\sigma^{2}I)$. An adaptation of the BYM was used
in the Part I paper and is also used in this paper. More details on spatial
models is available in our Part I paper, where we develop the spatial model
employed in the current paper.

The work done by Lawson and collaborators [18, 19, 20] is the closest to our
own. Fundamentally we extend an older model of ours [6, 4], meant for a single
areal unit, to encompass multiple areal units. Each areal unit contains its
own parameterized infection-rate representation, but are “stitched” together
using a BYM model. The parameters are estimated by solving an inverse problem,
conditioned on case-count data from areal units. Our model also includes a
model for the incubation period. In contrast, Lawson and co-workers model
case-counts directly; the clearest exposition of the model in in Lawson and
Song, 2010 [18], and it has been used with COVID-19 data from South Carolina
[21] and the UK [22]. Lawson and co-workers, much like us, have used the
disagreement between calibrated model forecasts and data as signs of
epidemiological anomalies and devised metrics such as the Surveillance
Kullback-Liebler [23] (SKL) and Surveillance Conditional Predictive ordinate
[24] (SCPO) to detect them.

Part I paper:  Our current paper is an algorithmic follow-on to our Part I
paper [7, 8], which focused on developing a spatial model, with a latent
(parameterized) infection-rate field, where the parameters were estimated from
COVID-19 data from the counties of NM. As mentioned above, it consists of a
parameterized infection-rate model for each areal unit [6, 4] (i.e., a NM
county) which are linked to each other via a (modification of) BYM model, much
like the work by Lawson and collaborators [18, 19, 20]. The spatial dependence
is between nearest-neighbors only; the details of how this was arrived at are
in the Part I paper. The model was fitted to data from three adjoining NM
counties - Bernalillo, Santa Fe and Vaelncia - using AMCMC; its lack of
scalability prevented us from broadening the model to more counties. It
therefore motivated the current paper, where we use scalable, but approximate,
variational inference to address the problem of infection-rate estimation and
forecasting across all 33 NM counties. As in Lawson’s work [23, 24] we
developed metrics to detect anomalous epidemiological behavior, cast as a
disagreement between forecasts (using a calibrated model) and data.
Specifically, we used the forecating capability to detect the arrival of the
Fall 2020 wave using data collected up to September 15, 2020, and had no
diffculties in the detection within a week of its arrival. We also tested the
anomaly detector on data collected up to August 15, 2020, long before the
arrival of the COVID-19 wave, to check the detector’s susceptibility to false
positives.

Variational Inference:  Inverse problems for model calibration often require
the approximation of intractable probability densities arising from Bayesian
inference. A standard approach based on sampling is Markov Chain Monte Carlo
(MCMC) [25] but suffers from scalability issues due to slow convergence rates
for high-dimensional and/or multimodal distributions [26, 27]. Variational
inference (VI) [28] provides an alternative to sampling techniques where
approximate inference is recast as seeking a member of a family of
approximating densities which minimizes a discrepancy measure such as KL-
divergence. Originally developed for probabilistic graphical models [29] where
some degree of analytical tractability is maintained, it has more recently
been extended to many-parameter models, such as those seen in deep learning,
through gradient-based iterative schemes adapted to a probabilistic setting
[30, 31, 32]. These techniques are often termed Stochastic Variational
Inference (SVI) and can exploit the automatic differentiation available in
large ML models. They offer significantly improved scalability over MCMC while
potentially sacrificing some approximation quality depending on the set of
approximating distributions used. VI has seen successful applications to a
number high-dimensional inverse problems in areas including medical
classification [33, 34] and segmentation [35, 36], computer vision and image
processing [37], natural language processing [38, 39], and physics-based
models [40, 41].

In addition to this broad range of application spaces, VI has more recently
been adopted in a growing number of epidemiological modeling challenges.
Neural ODEs have been extended to a Bayesian setting where they can be
calibrated with VI and applied to state space epidemiological models [42].
Model selection and dynamic causal modeling based on evolving real-world time
series using VI have been applied to COVID-19 outbreaks [43, 44] to provide
online forecasting tools. Dynamical system inference for spatio-temporal
modeling of infectious diseases using ODE and PDE formulations has been
carried out at the state scale and combined with Bayesian neural networks
[45]. Model calibration with VI has also been explored using alternatives to
standard state-space models such as graph-coupled Hidden Markov Models [46].
In additional to predictive model calibration, generative probabilistic
modeling using VI-based approximations of distributions has also been explored
to generate mission information about disease spread [47].

##  3 Formulation

Here we propose an epidemiological model to forecast infection-rates across
adjacent geographical regions and use these forecasts to detect emergent
outbreaks. The model is an extension of previous work by Safta et al. [6] and
Blonigan et al. [4] for epidemic forecasts over a single-region wave to multi-
region outbreak detection. In this section we will first describe the single-
region model and then present statistical approaches to estimate the model
parameters over adjacent geographical regions.

###  3.1 Epidemiological Model

The epidemiological model is defined by a spatio-temporally varying infection-
rate model and an incubation model given by

$\displaystyle f_{inf}(t;t_{0}^{r},k^{r},\theta^{r})$ $\displaystyle=$
$\displaystyle\frac{(\theta^{r})^{-k}(t-t_{0}^{r})^{k-1}}{\Gamma(k^{r})}\exp\left(-\frac{t-t_{0}^{r}}{\theta^{r}}\right)$
(1) $\displaystyle F_{inc}(t;\mu,\sigma)$ $\displaystyle=$
$\displaystyle\frac{1}{2}\mathrm{erfc}\left(-\frac{\log
t-\mu}{\sigma\sqrt{2}}\right)$ (2)

where the infection rate $f_{inf}$ is a Gamma distribution with shape and
scale parameters $k^{r}$ and $\theta^{r}$, respectively.
$F_{inc}(t;\mu,\sigma)$ is the incubation period for COVID-19, taken from
Lauer et al. [48]; further details are in the Part I paper. The parameter
$t_{0}^{r}$ represents the start of the outbreak and will be inferred along
with the infection-rate parameters. Note that $1\leq r\leq R$ indexes the
spatial region. The number of people that turn symptomatic over the time
interval $[t_{i-1},t_{i}]$ is given by

$\displaystyle y_{r}(i;t_{0}^{r},N^{r},k^{r},\theta^{r})$ $\displaystyle=$
$\displaystyle N_{r}\int_{t_{0}^{r}}^{t_{i}}f_{inf}\left(\tau-
t_{0}^{r};k^{r},\theta^{r}\right)\left[F_{inc}(t_{i}-\tau;\mu,\sigma)-F_{inc}(t_{i-1}-\tau;\mu,\sigma)\right]d\,\tau$
(3)

so that $\mathbf{y}(i)=\mathbf{y}_{i}=\left[y_{1}(i)\cdots
y_{R}(i)\right]^{T}$ and $i$ represents the time-dependence of the
predictions. Here, $N^{r}$ is the fourth and final region-dependent parameter
and represents the total number of people infected during the entire epidemic
wave in spatial region $r$ normalized by the population of region $r$.

The noisy model predictions are defined as

$\mathbf{y}^{(o)}_{i}=\mathbf{y}^{(p)}_{i}+\boldsymbol{\epsilon}_{i}=\mathcal{M}(t_{i};\mathbf{m})+\boldsymbol{\epsilon}_{i},\hskip
8.53581pt\boldsymbol{\epsilon}_{i}\sim\mathcal{N}(\mathbf{0},\mathbf{\Sigma}_{i}).$
(4)

Here $\mathbf{m}=\mathrm{vec}\left(\mathbf{m}_{r}\right)$, where
$\mathbf{m}_{r}^{T}=(t_{0}^{r},N^{r},k^{r},\theta^{r})$, are the region
specific model parameters for $r=1,\ldots,R$ and $\mathcal{M}(:)$ is an
epidemiological / disease model that uses $f_{inf}(:)$. In addition,
$\mathbf{y}^{(o)}_{i}$ is the vector of observed case-counts on day $i$ for
all $R$ regions, $\mathbf{y}^{(p)}_{i}$ is the vector of $R$ case-count
predictions for the same day using the model $\mathcal{M}$ and
$\boldsymbol{\epsilon}_{i}$ is the “noise” i.e. observed case-counts on day
$i$ that cannot be explained by $\mathcal{M}$ (mostly reporting errors). To
account for spatial correlations and heteroscedastic noise seen in case-
counts, the noise is assumed to be composed of two terms
$\boldsymbol{\epsilon}_{i}=\boldsymbol{\epsilon}_{i,1}+\boldsymbol{\epsilon}_{i,2}$
where the first is given by a Gaussian Markov Random Field (GMRF) model while
the second represents temporally-varying, independent Gaussian noise. Letting
$\mathcal{D}$ and $\Theta$ represent the data and parameters, respectively,
the likelihood then takes the form

$p(\mathcal{D}\ |\
\Theta)=\prod_{i=1}^{N_{d}}\frac{1}{(2\pi)^{N_{r}/2}\mathrm{det}\mathbf{\Sigma}_{i}^{1/2}}\exp\left(-\frac{1}{2}(\mathbf{y}^{(o)}_{i}-\mathbf{y}^{(p)}_{i})\mathbf{\Sigma}_{i}^{-1}(\mathbf{y}^{(o)}_{i}-\mathbf{y}^{(p)}_{i})^{T}\right)$
(5)

where $\mathbf{\Sigma}_{i}$ is given by

$\mathbf{\Sigma}_{i}=\tau_{\Phi}P^{-1}+\mathrm{diag}\left(\sigma_{a}+\sigma_{m}\mathbf{y}^{(p)}_{i}\right)^{2}.$
(6)

The first term $P=\left[D-\lambda_{\Phi}\mathbf{W}\right]$ in Eq. (6) forms
the precision matrix of a GMRF component of the noise where the strength of
correlations induced by adjacent regions is governed by $\lambda_{\Phi}$. The
relative topology of regions is encoded by $\mathbf{W}$, the county adjacency
matrix, defined as

$w_{kk}=0\,\hskip 8.53581pt\textrm{and}\,\hskip
8.53581ptw_{kl}=\begin{cases}1&\textrm{if regions {\it k} and {\it l} are
adjacent}\\\ 0&\textrm{otherwise}\end{cases}$ (7)

Here, $D=\mathrm{diag}\\{g_{1},g_{2},\ldots,g_{R}\\}$ where $g_{i}$ is the
number of regions adjacent to region $i$.

The second term captures prediction-dependent, uncorrelated noise with
additive and multiplicative components governed by $\sigma_{a}$ and
$\sigma_{b}$, respectively. The relative contribution between the correlated
GMRF noise and the uncorrelated noise is controlled by $\tau_{\Phi}$. Full
details, including the development of $P$, can be found in the Part I paper[7,
8]. The predictive skill and flexibility of this epidemiological model, as
applied to one areal unit, are described in detail in Safta et al. [6] and
Blonigan et al. [4] and the extension to joint estimation, described in the
Part I paper, achieves much the same predictive accuracy.

###  3.2 Statistical Inference

The set of parameters $\Theta$ defining the likelihood in Eq. (5) is given by

$\Theta=\mathrm{vec}\left(\theta_{i}\right)=\mathrm{vec}\left(\begin{bmatrix}\mathbf{m}_{1}&\cdots&\mathbf{m}_{R}&\mathbf{\eta}\end{bmatrix}\right)$
(8)

where $\mathbf{\eta}=(\tau_{\Phi},\lambda_{\Phi},\sigma_{a},\sigma_{m})$ are
the global noise parameters. Inference consists of forming the posterior
distribution $p\left(\Theta\ |\ \mathcal{D}\right)=p\left(\mathcal{D}\ |\
\Theta\right)p\left(\Theta\right)/p\left(\mathcal{D}\right)$ over uncertain
parameters $\Theta$. As the posterior is intractable, we instead look to
approximate it using VI. Hence, the following sections describe how VI is
formulated carried out to approximate the posterior $p(\Theta\ |\
\mathcal{D})$ for the outbreak model as well as how the prior $p(\Theta)$ is
defined to regularize the inverse problem.

####  3.2.1 Variational Inference

We will compare the Bayesian posterior sampled with AMCMC with posterior
models obtained using MFVI which recasts approximate inference as an
optimization problem. In particular, as the exact posterior is intractable, we
consider a family of approximating densities
$\mathcal{F}=\\{q(\Theta;\boldsymbol{\phi})\ |\
\boldsymbol{\phi}\in\boldsymbol{\Phi}\subseteq\mathbb{R}^{d}\\}$ and seek to
find a density $q(\Theta;\boldsymbol{\phi}^{*})$ that minimizes the KL-
divergence with respect to the posterior

$\boldsymbol{\phi}^{*}=\mathrm{argmin}\hskip
5.69054ptD_{\mathrm{KL}}\left(q(\Theta;\boldsymbol{\phi})\ \lVert\ p(\Theta\
|\ \mathcal{D})\right)$ (9)

This can be re-expressed as minimizing the objective function
$\mathcal{L}(\boldsymbol{\phi})$ based on the evidence lower bound (ELBO) [31]

$\mathcal{L}(\boldsymbol{\phi})=-\mathbb{H}\left[q(\Theta;\boldsymbol{\phi})\right]-\mathbb{E}_{q(\Theta;\boldsymbol{\phi})}\left[\log
p(\mathcal{D}\ |\ \Theta)+\log p(\Theta)\right]$ (10)

where the first term in Eq. (10) is the entropy of the surrogate posterior and
the second, data-dependent term is an expectation with respect to the
surrogate posterior that reflects both the expected data-fit and the prior.
Here we take $\mathcal{F}$ to be the set of mean-field Gaussian distributions,
i.e.,

$q\left(\Theta;\
\boldsymbol{\phi}\right)=\prod_{i=1}^{d}q_{i}\left(\theta_{i};\
\mu_{i},\sigma_{i}\right)$ (11)

where
$q_{i}(\theta_{i};\mu_{i},\sigma_{i})=\mathcal{N}(\theta_{i};\mu_{i},\sigma_{i})$,
$\boldsymbol{\phi}=(\boldsymbol{\mu},\boldsymbol{\sigma})$ and we arrive at an
optimization problem over $2d$ parameters where $d$ is the number of
parameters defining the epidemiological model $\mathcal{M}$. To carry out the
above minimization problem, we aim to use a gradient-based iterative scheme as
the expectation in Eq. (10) cannot be evaluated explicitly due to the
nonlinearity of the forward model. Furthermore,
$\mathcal{L}(\boldsymbol{\phi})$ is potentially a non-convex objective. Note
that the gradient and expectation operators do not commute, i.e.,

$\nabla_{\boldsymbol{\phi}}\mathbb{E}_{q\left(\Theta;\boldsymbol{\phi}\right)}\left[\log
p\left(\mathcal{D}\ |\
\Theta\right)p(\Theta)\right]\neq\mathbb{E}_{q\left(\Theta;\boldsymbol{\phi}\right)}\left[\nabla_{\boldsymbol{\phi}}\log
p\left(\mathcal{D}\ |\ \Theta\right)p\left(\Theta\right)\right]$

so some care has to be taken to arrive at a Monte Carlo estimator for the
gradient $\nabla_{\boldsymbol{\phi}}\mathcal{L}(\boldsymbol{\phi})$. Two
widely used approaches are: (a) the score function estimator, described in §
A.1 (in the Appendix), which forms the basis of black-box VI and requires only
evaluations of the log-likelihood, and (b) the reparametrization approach
which requires gradients of the log-likelihood. The score function estimator
typically displays much larger variance as seen in Kucukelbir et al. [49]
where two orders of magnitude more samples were needed to arrive at the same
variance as a reparametrization estimator. A similar trend was confirmed for
the outbreak problem (Fig. 8) suggesting that the reparametrization approach
would lead to superior scalability. Reparametrization proceeds by expressing
$\Theta$ as a differentiable transformation
$\Theta=t(\boldsymbol{\epsilon},\boldsymbol{\phi})$ of a
$\boldsymbol{\phi}$-independent random variable $\boldsymbol{\epsilon}\sim
q(\boldsymbol{\epsilon})$ such that
$\Theta(\boldsymbol{\epsilon},\boldsymbol{\phi})\sim
q(\Theta,\boldsymbol{\phi})$. This allows the gradient to be expressed as

$\nabla_{\boldsymbol{\phi}}\mathcal{L}\left(\boldsymbol{\phi}\right)=-\nabla_{\boldsymbol{\phi}}\mathbb{H}\left[q\left(\Theta;\boldsymbol{\phi}\right)\right]-\mathbb{E}_{q\left(\boldsymbol{\epsilon}\right)}\left[\nabla_{\boldsymbol{\phi}}\log
p\left(\mathcal{D}\ |\
\Theta(\boldsymbol{\epsilon},\boldsymbol{\phi})\right)+\nabla_{\boldsymbol{\phi}}\log
p\left(\Theta\left(\boldsymbol{\epsilon},\boldsymbol{\phi}\right)\right)\right]$
(12)

where gradients of the entropy term in Eq. (12) are available analytically for
the Gaussian surrogate posterior and the second term can now be approximated
with Monte Carlo given a method to compute the required gradients. For many
machine learning models, automatic differentiation can be exploited to
calculate the gradient of the log-likelihood with respect to parameters
$\Theta$. Here, the objective function involves the log of the likelihood (Eq.
(5)) where derivatives of matrix inverses and determinants with respect to
parameters are required to compute the gradient. Gradients such as these are
not available using most automatic differentiation libraries. Instead, matrix
calculus and quadrature were used to compute the derivatives of the log-
likelihood with respect to model predictions $\mathbf{y}^{(p)}_{i}$ and to
approximate the derivatives of the model predictions with respect to
parameters, respectively. For details, see § A (in the Appendix).

Note that some of the parameters comprising $\Theta$ are required to satisfy
constraints for the noise and epidemiological models to be well-defined. For
example, noise parameters $\tau_{\Phi}$, $\sigma_{a}$, and $\sigma_{m}$ as
well as model parameters $N^{r}$, $k^{r}$ and $\theta^{r}$ for $r=1,\ldots,R$
should be positive while $\lambda_{\Phi}$ should satisfy
$0\leq\lambda_{\Phi}<1$. Sampling from the mean-field Gaussian (Eq. (11))
during the Monte Carlo estimation of the gradient (Eq. (12)) may result in
violations of these constraints. To maintain the required properties without
resorting to constrained optimization, we express a constrained parameter
$\theta_{i}$ as an invertible, differentiable transformation
$\theta_{i}=f_{i}(\hat{\theta}_{i})$ of an unconstrained $\hat{\theta}_{i}$.
Hence, the distribution governing $\theta_{i}$ is the push-forward density of
$\mathcal{N}(\hat{\theta}_{i};\mu_{i},\sigma_{i})$ through $f_{i}$, i.e., the
components of the mean-field surrogate posterior (Eq. (11)) have modified
probability densities

$q_{i}(\theta_{i};\mu_{i},\sigma_{i})=\mathcal{N}(\hat{\theta}_{i};\mu_{i},\sigma_{i})|f_{i}^{\prime}(\hat{\theta}_{i})|^{-1};\hskip
5.69054pt1\leq i\leq d$ (13)

where $\hat{\theta}_{i}=f_{i}^{-1}(\theta_{i})$. This results in mean-field
approximation where some of the factors are Gaussian and others non-Gaussian.
Each factor is still defined by a $\mu_{i}$ and $\sigma_{i}$ parameter. The
transformations are listed in Table 2.

Note that MFVI has a tendency to underestimate the uncertainty / variance in
the estimated parameters[10] and in § 4, we will check if this poses a
limitation on the usefulness of the approximate solution, especially in the
predictive skill of model $\mathcal{M}$ and therefore the detection of the
start of outbreaks.

####  3.2.2 Prior distribution

The COVID-19 case count data exhibits significant noise due to inaccurate case
counting reported by hospitals. Furthermore, counties with small populations
exhibit sparse data in the sense that not many positive daily case counts were
reported. Hence, we expect the inverse problem to be ill-posed and require
regularization in the form of a prior $p(\Theta)$ over the parameters.

Because of push-forward formulation described by Eq. (13), a number of the
parameters are already constrained by transformations
$\theta_{i}=f_{i}(\hat{\theta}_{i})$. In particular, the parameters
$N^{r},k^{r},\theta^{r}$ for $r=1,\ldots,R$ and each of the noise parameters
comprising $\mathbf{\eta}$ are all constrained by transformations to take on
values in some restricted interval. For example, $\lambda_{\Phi}$ is
constrained to lie within $[0,1-\epsilon]$, for some $\epsilon>0$ so that Eq.
(6) defines a valid covariance matrix, i.e., it remains symmetric, positive
definite. The parameters $t_{0}^{r}$ are the only unconstrained variables.
Hence, we take Gaussian priors over $t_{0}^{r}$ that incorporate diffuse
assumptions about when it is reasonable for a wave to occur.

####  3.2.3 Posterior predictive tests

The mean-field variational inference (MFVI) described in § 3.2.1 results in a
multivariate Gaussian posterior distribution (Eq. (13)) that then needs to be
verified against data $\mathcal{D}$. To do so, we take samples $\Theta_{j}\sim
q\left(\Theta;\boldsymbol{\phi}\right)$, and using Eq. (4), generate
predictions
$\boldsymbol{Y}_{j}=\\{\mathbf{y}^{(p)}_{i}+\boldsymbol{\epsilon}_{i}\\}_{j},j=1\ldots
J$. For this paper, $J=100$. The time-series $\boldsymbol{Y}_{j}$ result in a
“fantail” of predictions for each areal unit $r$ which should statistically
reproduce $\mathcal{D}$ e.g., the inter-quartile range of the samples
$\boldsymbol{Y}_{j}$ should bound 50% of the observations and the $5^{th}$ and
$95^{th}$ percentiles should bound 90% of the individual data points in
$\mathcal{D}$. This test is called a posterior predictive test (PPT) which we
will use extensively in § 4. In addition, we define a score to summarize the
agreement of $\boldsymbol{Y}_{j}$ with $\mathcal{D}$. Predictive distributions
constituted out of $\\{\mathbf{y}^{(p)}_{i}\\}_{j},j=1\ldots J$ samples are
called “push-forward“ (PF) predictions.

Define $y_{r,i}^{(pred)}(\Theta)=y_{r}(i;\
\mathbf{m}_{r})+\epsilon_{r}(\boldsymbol{\eta})$ to be the model prediction
for region $r$ for day $i$ by combining Eq. (4) and Eq. (8). Let
$y_{r,i,j}=y_{r,i}^{(pred)}(\Theta_{j})$ be the prediction corresponding to a
sample $\Theta_{j}$ drawn from the posterior. Let $y^{(obs)}_{r,i}$ be the
corresponding observation (from $\mathcal{D}$). Let
$G_{r,i}(y^{(pred)};\Theta)$ be the cumulative distribution function (CDF) for
the model predictions arising from the posterior distribution
$q(\Theta;\boldsymbol{\phi})$. The Continuous Ranked Predictive Score (CRPS
[50]) is defined as

$c_{r,i}=\int\left(G_{r,i}(\upsilon)-\mathbbm{1}_{\upsilon>y^{(obs)}_{r,i}}\right)^{2}d\upsilon\hskip
8.53581pt\mbox{and}\hskip
8.53581ptC_{r}=\frac{\sum^{N_{d}}_{i}c_{r,i}}{N_{d}}.$ (14)

The CRPS has units of case-counts and will be larger for areal units with
larger total case-counts $T_{r}=\sum_{i}y^{(obs)}_{r,i}$, and we will use the
ratio $\rho_{r}=C_{r}/T_{r}$ to compare across areal units in § 4. In
practice, the empirical CDF computed using $J$ samples of $\boldsymbol{Y}_{j}$
is used to approximate $G_{r,i}(y^{(pred)};\Theta)$. This score function has
been used in judge the quality of the model to capture the spread in the data
used for calibration. [51, 6, 52] We use the implementation in the R
Statistical Software[53] (R version 4.3.2 (2023-10-31)) package
verification[54], specifically the function crpsDecomposition(), to compute
the CRPS.

##  4 Results

In this section, the calibration of the outbreak model, using the formulation
presented in § 3, is studied across several cases. The COVID-19 pandemic
arrived in NM in March 2020; Fig. 1 (left) plots the detected cases in NM over
2020. We see three clear “waves“ - the Spring wave, the Summer wave which
spanned June $1^{\text{st}}$ to September $15^{\text{th}}$, and the Fall wave
that arrived after that. Unless specified otherwise, data from the Summer wave
(i.e., June 1, 2020 to September 15, 2020) is used to perform the estimation
and the estimated infection-rate is used to forecast two weeks ahead, into the
Fall wave. The COVID-19 dataset of case-counts that we use covers the duration
from 2020-01-22 to 2022-05-13, and consists of daily (new) case-counts of
COVID-19 from each of the 33 counties of NM; the data is available online.
[55, 56]. The sudden change in epidemiological dynamics around September
$15^{\text{th}}$ implies that that the infection-rate profile from the Summer
wave will not be able to forecast the Fall wave; therefore any disagreement
between the observed data and forecasts is a sign of the arrival of the Fall
wave.

Figure 1: Left: The three COVID-19 waves in NM in 2020. We see the Spring
wave, ending around June $1^{\text{st}}$, followed by the Summer wave and the
Fall wave that started around September $15^{\text{th}}$. Our aim is to
estimate the infection-rate _field_ during the Summer wave and detect the Fall
2020 wave. The solid vertical line denotes September $15^{\text{th}}$ and the
dashed one is placed at August $15^{\text{th}}$. Right: The couties of New
Mexico. The shaded ones are where we will present results. Abbreviations: B =
Bernalillo; CI = Cibola; CU = Curry; DA = Doña Ana; RA = Rio Arriba SF = Santa
Fe and V = Valencia.

First, we consider three populated counties that are adjacent to each other
and display a large number of cases. These counties present less noise and
clearer trends resulting in a more well-posed inversion task. The aim of the
study is to compare the effect of calibrating multiple areal units jointly
using MFVI, as well as to compare against the calibration performed using
AMCMC in our Part I paper[7, 8]. Next, the coupled outbreak model is
calibrated across all 33 New Mexico counties. This represents a more
challenging task given the dimensionality of the problem as well as a
multitude of counties displaying sparse and noisy case-counts. We then
evaluate the final posterior approximation using PPT runs. The convergence of
the MVFI procedure for both the outbreak model and noise parameters is
investigated and discussed last. We next discuss anomaly detection using the
MFVI-calibrated outbreak model is carried out using COVID-19 spread-rates
across all 33 counties in New Mexico using data from the summer of 2020 to
detect the arrival of the Fall 2020 COVID-19 wave.

When estimating the infection-rate field in the 3-counties, we infer $3\times
4+4=16$ parameters – 4 parameters $(t_{0}^{r},N_{r},k^{r},\theta^{r})$ for
each county and 4 noise parameters
$(\tau_{\phi},\lambda,\sigma_{a},\sigma_{m})$; this was performed with AMCMC
in our Part I paper. For all of NM, with its 33 counties, the dimensionality
of the inverse problem is $33\times 4+4=136$ independent parameters. This is
too large for AMCMC and was the motivation for developing MFVI.

The MFVI procedure follows § 3 where the stochastic gradient descent iteration
is carried out using the Adaptive Moment Estimation (ADAM) algorithm [57].
MFVI is well-known to display mode-seeking behavior due properties of the KL-
divergence. Hence, to facilitate the convergence of MFVI, we initialize the
mean parameters for MFVI from a Maximum Likelihood Estimate (MLE) which is
readily available using gradients of the log-likelihood. During calibration
via stochastic gradient descent, $N_{s}=200$ samples were used in the Monte
Carlo estimates of the ELBO gradient Eq. (12).

###  4.1 Three-county inversion

In this section we investigate the effect of using an approximate MFVI
inversion by comparing against AMCMC solution computed in our Part I paper. We
also check the effect of estimating the infection-rate across multiple areal
units vis-à-vis independently.

Figure 2: Comparison of the predictive distribution for the MFVI inversion of
Bernalillo (left), Santa Fe (middle), and Valencia (right) done jointly using
the _GMRF_ model (top) and _independently_ for each county (bottom). Case-
count data was smoothed with a 7-day running average. The red line shows the
median predictions, the shaded green line shows the inter-quartile range, and
the dashed lines are the $5^{\text{th}}$ and $95^{\text{th}}$ percentiles. The
filled symbols are the data $\mathcal{D}$ used in the inversion. The unfilled
symbols are the observations beyond September $15^{\text{th}}$ and are used to
compare with the two-week forecasts.

Joint versus independent calibrations: First, we perform an infection-rate
estimation for three adjoining counties - Bernalillo, Santa Fe and Valencia -
using MFVI independently and jointly using the spatial (GMRF) model (Eq. (6));
see Fig. 1 (right) for their positions. In Fig. 2, the results of a PPT for
all three counties are displayed in both the joint and independently
calibrated cases. The case-count data were smoothed with a 7-day running
average. The median prediction is plotted with the red line and the dashed
lines denote the $5^{\text{th}}$ and $95^{\text{th}}$ percentiles. We see that
most of the observations (filled symbols) are within these bounds. We also see
the 2-week-ahead forecast beyond September $15{\text{th}}$ and the unfilled
symbols showing the observed data from that period. The forecasts and the data
do not match for any of the three counties, indicating a change in the
epidemiological dynamics. This change is due to the arrival of the Fall 2020
wave of COVID-19 in NM, and the figure shows that the wave arrived
approximately simultaneously in all three counties (see Fig. 1 (left) for a
clearer picture of the three COVID-19 waves encountered in 2020 in NM).
Observe that the jointly calibrated predictive distribution displays less
uncertainty in the predictions (i.e., small $\sigma_{m}$) than the
independently calibrated version. This is likely due to a combination of two
factors. The first is that incorporated spatial correlations between counties
regularizes the calibration and results in more certainty about the true
underlying model parameters. The second is that the uncorrelated MFVI
approximation is known to underestimate uncertainty for highly correlated
distributions. Note also that uncertainty is largest in the predictive
distribution for Santa Fe consistent with the markedly noisier behavior of the
case-counts for this county. The infection-rate profiles for these counties
are in Fig. 13 in the Appendix. There is not much difference between them
indicating that the infection-rate parameter estimates in the two cases might
be similar, whereas the estimate of the noise, which affects PPT results, vary
between the two formulations due to the fact that the spatial noise model has
to accommodate all the noise, across all the areal units, together. These
findings echo what we observed in our Part I paper, where the same study was
performed using AMCMC as the estimation procedure.

Figure 3: Comparison of the predictive distribution for the joint inversion of
Bernalillo (left), Santa Fe (middle), and Valencia (right) done jointly using
_MFVI_ (top) and _AMCMC_ (bottom). Case-count data was smoothed with a 7-day
running average. The red line shows the median predictions, the shaded green
line shows the inter-quartile range, and the dashed lines are the
$5^{\text{th}}$ and $95^{\text{th}}$ percentiles. The filled symbols are the
data $\mathcal{D}$ used in the inversion. The unfilled symbols are the
observations beyond September $15^{\text{th}}$ and are used to compare with
the two-week forecasts. The AMCMC results are taken from our Part I paper[7,
8].

AMCMC versus MFVI estimation: Next, we study the effect of using our
approximate MFVI method versus the estimates computed using AMCMC in the Part
I paper. In Fig. 3 we plot the PPT results from the joint MFVI (top) and AMCMC
(bottom) for the same three counties. Here we see that the MFVI estimates a
larger $\sigma_{m}$ \- the $5^{\text{th}}$ and $95^{\text{th}}$ percentile
bounds (dashed lines) are far wider for MFVI results vis-à-vis AMCMC results
below. By dint of having wider bounds, the MFVI estimate is also better at
bounding the data used to compute the infection-rate _field_ for the three
counties. Again the arrival of the Fall 2020 wave is clearly discerned in the
figure. The infection-rate profiles for these counties are in Fig. 14 in the
Appendix and there is not much difference between them.

Figure 4: Comparison of the posterior over model parameters $t_{0}$, $N$, $k$,
and $\theta$ for 3-county inference of Bernalillo (left), Santa Fe (middle),
and Valencia (right) using both MCMC (orange) and MFVI (blue). Both joint
(solid lines) and independent calibration (dashed lines) are results are
displayed. $t_{0}$ values are negative as it is measured from June
$10^{\text{th}}$, 2020, and the PDFs imply that infections for the Summer wave
started in late May. ’VI indep.’ in the legend implies an estimate obtained
for each county independently using MFVI.

In Fig. 4, we summarize the marginalized posteriors for the infection-rate
parameters for Bernalillo, Santa Fe, and Valencia computed using MFVI and
AMCMC, jointly and independently. $t_{0}^{r}$ assumes negative values as it is
measured from June $10^{\text{th}}$, 2020; the PDFs peak around -20 (for AMCMC
results) and imply that the infections for the Summer wave started in late
May, about 20 days before June $10^{\text{th}}$. The MFVI results also agree
approximately with this estimation. We see that apart from $N$, the MFVI and
AMCMC posteriors do not match. The MFVI posteriors are extremely narrow
providing a spurious degree of certainty in the estimates; this arises from
the form of the posterior distribution - independent Gaussians - that we
postulate in Eq. 11. In addition, as is clear from the AMCMC results, the
“true” posteriors are not Gaussian. The marginalized posteriors computed via
AMCMC, for joint and independent estimation, do match (sometimes very well, as
in the case of Valencia), but they are wide apart for the MFVI, showing the
effect of the Gaussian approximation. Yet the infection-rate profiles in Fig.
13 (in the Appendix) do not show much of a difference, nor do the PPT results
in Fig. 2, leading us to conjecture that the influence of some of these
parameter on the infection-rate may be muted. This can also be surmised from
the marginal distributions computed from AMCMC - they are quite wide.

Finally, we summarize the predictive skill of the joint estimates computed
using MFVI and AMCMC in Table 1 using CRPS i.e. $C_{r}$. The CRPS is computed
over data between June $1^{\text{st}}$ and September $15^{\text{th}}$, 2020
i.e., the CRPS $C_{r}$ summarizes the agreement of the PPTs with observations
for each of the counties. We see that the CRPS (a measure of the error between
predictions and data) is about 10 cases per day, for Bernalillo (where case-
counts peaked at about 100 cases/day; see Fig. 2) and 2.5 cases a day for
Valencia and Santa Fe (which peaked at about 12 cases a day). What is
remarkable is the difference in the CRPS as computed using MFVI and AMCMC - it
is less than a case-count per day. This small change in the predictive skill
leads us to believe that the scalable, but approximate, MFVI approach might be
sufficiently accurate to allow us to detect the arrival of the Fall 2020 wave
in an automated manner.

_To summarize_ , joint estimation may lead to more certain forecasts (narrower
PPTs) vis-à-vis independent estimation of infection-rates in areal units
(counties). MFVI provides spuriously low levels of uncertainty in estimated
parameters. The marginalized posteriors of infection-rate parameters obtained
from MFVI and AMCMC do not agree, but the infection-rate does not seem to be
very sensitive to the disagreement. Consequently, the PPTs from the AMCMC and
MFVI inversion are very similar in terms of their predictive skill.

Table 1: Predictive skill of PPTs generated using infection-rate estimates
computed using different procedures. All estimations are performed using the
joint formulation, using the spatial model. The PPTs are scored using CRPS.
They have units of “case-counts”.

Procedure | Bernalillo | Santa Fe | Valencia
---|---|---|---
AMCMC, 3-county${}^{{\textsuperscript{\bf a}}}$ | 11.3 | 2.65 | 1.76
MFVI, 3-county${}^{{\textsuperscript{\bf b}}}$ | 10.75 | 2.54 | 1.70
MFVI, 33-county | 10.7 | 2.35 | 1.47

  * a

Taken from Table 2 of our Part I paper[7, 8].

  * b

Computed from Fig. 2 (top).

###  4.2 Joint inversion of all NM counties

Figure 5: $C_{r}/T_{r}$ plotted as a function of $T_{r}$. Both the axes are
log-axes. We see, as expected, that the predictive skill of the disease model,
post calibration, is better for counties with larger case-counts $T_{r}$ where
noise variance is low. The horizontal dashed lines are the first and third
quartiles of $\rho=C/T$.

Next, calibration of the full 33-county model using MFVI was carried out,
followed by PPT runs which were then summarized using CRPS i.e., $C_{r}$ was
computed for each NM county. In Fig. 5, we plot the CRPS normalized by the
total number of cases i.e., $C_{r}/T_{r}$ as a function of $T_{r}$, for
$r=1,\ldots,R$. Data between June $1^{\text{st}}$ and September
$15^{\text{th}}$ was used in the infection-rate estimation as well as the
computation of CRPS i.e., $C_{r}$ is a measure of the "goodness of model fit"
to data. We see that the CRPS, as a ratio of the total cases, decreases with
increasing number of total cases, as the disease model fits larger outbreaks
in counties like Bernalillo. Others, like Cibola, do not agree with model
predictions, due to flaws in the data, as we will see later in Fig. 11. The
straight line fit to data has the form

$\log(\rho)=\log\left(\frac{C}{T}\right)\approx 1.1-0.28\log(T),$

where $C$ and $T$ are the CRPS and total number of cases. This scaling shows
that $\rho$ scales (approximately) as the fourth root of the total number of
cases, a slow reduction indeed.

The horizontal lines in Fig. 5 show the first and third quartiles of $\rho$.
The three counties studied in our Part I paper, Bernalillo, Santa Fe and
Valencia, are marked and fall in the lowest quartile i.e., their data is good
and the calibrated disease model is predictive. We will use the county of Doña
Ana as another member of the “good” class, while Curry and Rio Arriba, which
fall in the inter-quartile range of $\rho$, will serve as exemplars of the
“middling” class of calibration. Cibola, which falls in the last quartile,
will be an exemplar of the “bad” class of calibration. We now examine the
quality of the inversion in each of these counties.

Figure 6: Means (dark) and standard deviation (gray) of various infection-rate
parameters, for select counties marked in Fig. 5. Top left: The start time of
the Summer wave $t^{r}_{0}$. Top right: The total size of the Summer outbreak
$N_{r}$. Bottom left: $k^{r}$, the shape parameter of the Gamma profile of the
infection-rate. Bottom right: $\theta^{r}$, the scale parameter of the Gamma
profile.

In Fig. 6 we plot the infection-rate parameters, as estimated from the Summer
wave case-count data, for 7 counties marked in Fig. 5. We see that the
standard deviations are spuriously tiny, in line with what was observed in
Fig. 4; thus MFVI consistently underestimates the uncertainty in the parameter
estimates. Further, this spuriously low uncertainty is pervasive - counties
with high CRPS such as Cibola and Rio Arriba show much the same estimation
uncertainties as counties such as Bernalillo and Santa Fe with CRPSs a factor
of three smaller. Thus the uncertainty in the parameters’ estimates do not
seem credible and we will omit them from further discussion. However, despite
the quality of the case-count data, the inversion completes _stably_ and
provides plausible results. However, as Fig. 5 shows, the infection-rate may
not be estimated very accurately in some counties and this might hamper the
task of detecting the Fall 2020 wave reliably.

Figure 7: Some infection-rate and noise parameters, estimated jointly among
three counties, compared with their counterparts from a 33-county inversion.
The noise parameter $\sigma_{m}$ is markedly larger in the 33-county
inversion.

In Fig. 7 we compare some infection-rate parameters
$(t^{r}_{0},N_{r},k^{r},\sigma_{m})$ for Bernalillo, Santa Fe and Valencia, in
the 3-county (as plotted in Fig. 4) and 33-county inversions. MFVI was used
for both the computations. We see that the parameters are not the same and
there does not seem to be a clear trend, except that $\sigma_{m}$ increases
when we include counties (some with poor quality data) in the estimation. This
implies that the PPTs for the “good” counties will likely be wider than what
could be achieved with AMCMC and this might have repercussions regarding
detecting the Fall 2020 wave. The use of the AMCMC results to detect the Fall
2020 wave is described in our Part I paper - the wave was detected within a
week of its arrival.

_To summarize_ , the uncertainty in the estimated parameters, computed using
MFVI, are not very credible. However, the MFVI inversion is stable, though the
uncertainty in the PPTs will be larger because of an inflated $\sigma_{m}$,
required to accommodate counties with very noisy data.

###  4.3 Algorithmic results

Finally, we investigate the numerical aspects of the reparametrized algorithm
described in § 3.2.1. The convergence of MFVI is depicted in Fig. 8 where the
ELBO and the norm of its gradient are shown as a function of gradient descent
iterations.

Figure 8: (Top) Convergence of the ELBO for the 33-county inversion along with
the norm of the ELBO gradient as a function of gradient descent iterations for
the reparametrized gradient formulation of MFVI. (Bottom) Convergence of the
ELBO for a 1-county inverse problem along with the norm of the ELBO gradient
for the black box formulation of MFVI. In both cases, $N_{s}=300$ samples were
used for the MC estimators of the gradient.

The top of Fig. 8 shows the ELBO and gradient for the 33-county calibration
using the reparametrization formulation while the bottom provides a comparison
to using black box VI for 1-county calibration of Bernalillo. In both cases,
$N_{s}=300$ samples were used for the reparametrization and score function MC
estimators of the gradient. Note that even for a single county, black box VI
shows significantly higher variance in the gradient leading to poor
convergence in comparison to the much larger 33-county problem calibrated with
reparametrization. The convergence of the ELBO in the top row is also quite
smooth suggesting that significantly less samples could be used to obtain good
estimates of the gradient with the reparametrization approach. Hence, it is
clear that reparametrization is necessary to scale the calibration to the
33-county inversion despite the added complexity of obtaining gradients of the
log-likelihood.

The top two rows of Fig. 9 display convergence information from two counties,
Bernalillo and Rio Arriba, taken from the 33-county inversion. Bernallilo and
Rio Arriba were chosen as they have larger and smaller populations,
respectively. The mean of the initial condition for MFVI is given by a MLE
solution shown in red. Intermediate solutions are shown in blue along with the
final solution in green. Observe that while similar to the MLE, MFVI subtly
expands the shape of the wave to better cover the tail of the outbreak. This
is potentially an effect of the tendency of the KL-divergence to increase the
overlap of the surrogate and true posterior distributions at some expense of
the mean prediction fitting the data less accurately. Comparing Fig. 2 to Fig.
9, we can see that the coupled, 33-county inversion introduces a bias in the
parameter estimates for Bernalillo. This is due to the multitude of NM
counties that are sparsely populated and display significant noise in their
daily case counts. Despite this effect, the outbreak detection performs well
on the 33-county inversion data suggesting that the predictive uncertainties
provided by the calibration remain informative.

Figure 9: (Top) Visualization of the convergence of 33-county inversion using
MFVI in terms of model predictions for two counties: Bernalillo (left) and Rio
Arriba (right). The mean of the initial condition for MFVI is given by a MLE
solution shown in red. Intermediate solutions are shown in blue along with the
final mean solution in green. (Middle) Corresponding convergence of the model
parameters for Bernalillo and Rio Arriba. (Bottom) The ELBO objective (left)
and its gradient (right) as a function of iteration for the full 33 county
inversion. Here, day zero of the inversion is defined as June $1^{\text{st}}$,
2020.

_To summarize_ , reparametrization provided us with the scalability needed to
solve the high-dimensional inverse problem conditioned on data from all 33 NM
counties.

##  5 Discussion

The results in § 4 show that MFVI can estimate an infection-rate field and the
posterior distribution can be used to produce PPT runs. The MFVI parameter
estimates do not quite agree with the AMCMC estimates from the Part I paper[7,
8], but their effect on the PPT runs is muted, as seen in the comparison in
Fig. 3 and the CRPS summaries in Table 1. Given that the inversion is a
smoothing operation, i.e., we learn the infection-rate from historical data,
any forecast produced with the estimated infection-rate will be predictive
only if the epidemiological dynamics do not change. Consequently, if the
forecast and data disagree, it could indicate the arrival of a new wave of
infection.

Temporal detection: The argument above was used to fashion an outbreak
detector using PPT runs in the Part I paper. The detector works as follows. We
sample the posterior in the same manner as for PPT runs, and use
$\boldsymbol{Y}$ to compute an “outlier boundary”; in line with our Part I
paper, we define it as the ${\rm 99^{th}}$ percentile of the forecasts.
$\boldsymbol{Y}$ is computed using data from June $1^{\text{st}}$ to an end
date (usually August $15^{\text{th}}$ or September $15^{\text{th}}$) where we
test for a change in the epidemiological dynamics. The actual test consists of
comparing a two-week-ahead forecast with the data that was observed during
that period. Any day with case-counts above the “outlier boundary” is deemed
an outlier. Three consecutive outlier days cause an “alarm”, corresponding to
an anomalous change in the disease dynamics. Using this detector, in the Part
I paper, we found that we could detect the arrival of the Fall 2020 wave
correctly, when tested using data up to September $15^{\text{th}}$. When using
data up to August $15^{\text{th}}$, a month _before_ the arrival of the Fall
2020 wave, we did not encounter a false positive. The method was also compared
against a conventional detector[58], which performed poorly - the reason for
this performance is discussed in the Part I paper. The infection-rate
estimation was performed using AMCMC, for Bernalillo, Santa Fe and Valencia
jointly. In Fig. 10 we repeat the same test, but the infection-rate is
estimated using MFVI. The top row depicts the outbreak detector being applied
beyond September $15^{\text{th}}$, 2020. We see outliers and alarms for all
three counties with a week of September $15^{\text{th}}$ i.e., the Fall 2020
wave was easily detected within a week’s worth of data. The bottom row repeats
the test, but applied to August $15^{\text{th}}$. We see that Santa Fe and
Valencia incur false positives, whereas Bernalillo does not. This implies that
the approximations in MFVI may lead to erroneous detections for borderline
cases (i.e., counties with small-count data and high variance noise) but areal
units with large case-counts might be unaffected.

Figure 10: The infection-rate detector implemented using the infection-rate
from the joint MFVI estimation of three counties (corresponding to Fig. 3
(top)). The symbols are the case-counts, the solid blue line the median
forecast with the infection-rate, post-calibration and the solid red line is
the alarm boundary. Outliers are circled and an alarm, corresponding to three
consecutive outliers, is encased in a square. The vertical line is the point
beyond which we forecast and thus test for the arrival of the Fall 2020 wave.
Top row: Testing for arrival after September $15^{\text{th}}$, 2020. Bottom
row: Testing for arrival after August $15^{\text{th}}$, 2020.

The reason for the false positives is simple - the data is very noisy and the
outbreak detector makes no attempt to reduce the variance in the noise.
Comparing the correct detection on the top row of Fig. 10 with the false
positives in the bottom row, we see that outliers and alarms are plentiful
after September $15^{\text{th}}$ and sporadic after August $15^{\text{th}}$. A
slightly more sophisticated detector that performed temporal averaging could
eliminate such false positives. We will address this below, using a simple
spatio-temporal method.

Figure 11: Performance of the outbreak detector for a “good” (Doña Ana), two
“middling” (Curry and Rio Arriba) and one “bad” county (Cibola), as determined
using CRPS in Fig. 5. Left column: Testing for Fall 2020 wave arrival on
September $15^{\text{th}}$, 2020. Right column: The same, but detection
performed on August $15^{\text{th}}$. We do not see any false positives in the
right column.

In Fig. 11, we demonstrate the outbreak detector on a “good” (Doña Ana), two
“middling” (Curry and Rio Arriba) and one “bad” county (Cibola), as determined
using CRPS in Fig. 5. The plotting conventions are the same as in Fig. 10. We
see that we correctly detect the Fall 2020 wave when testing on September
$15^{\text{th}}$ and do not incur any false positive on August
$15^{\text{th}}$. This performance is rather fortuitous for Rio Arriba and
Cibola, which have unexplained case-count peaks in mid-July (Rio Arriba) and
in late July - early August (Cibola). It would be difficult to estimate an
infection-rate profile for either of these counties independently (this is
especially true for Cibola, where the _reported_ case-counts show no wave-like
structure); however, our spatial GMRF model regularizes the inversion and
constructs an infection-rate stably. Forecasts using this infection-rate do
not match the case-count data at all, as is clear in the plots for Cibola, but
provides us with stable PPT results and a perhaps meaningless detection.
However, numerical and algorithmic stability in the face of poor quality data
encourage us to believe that the MFVI algorithm and outbreak detection can be
automated.

Figure 12: Top: Plot of mean exceedance $\overline{\gamma_{r}}$ for NM, where
$\overline{\gamma_{r}}$ is averaged over August $15^{\text{th}}$ and August
$31^{\text{st}}$ (left) before the arrival of the Fall 2020 wave and September
$16^{\text{th}}$ and October $1^{\text{st}}$ (right), when the Fall 2020 wave
had arrived. Bottom left: Dendrogram from hierarchical clustering, with the
cut at a height of 0.6, resulting in clustering of counties. Bottom right:
Disease clusters from the dendrogram.

Detecting spatial patterns:  The results above show that the MFVI estimates of
the infection-rate have the ability to detect the Fall 2020 wave, especially
for the counties in the first quartile (plotted in Fig. 5), though our crude
detector, which does not smooth / de-noise the data, might suffer from false
positives due to high variance noise in observed case-count. However, the MFVI
infection-rate _field_ allows prediction in space and time, which allows
spatio-temporal assimilation of data. We address this next.

Define “exceedance” $\gamma_{r,i}=y^{obs}_{r,i}/y^{(99)}_{r,i}$ for region $r$
and day $i$. Here $y^{obs}_{r,i}$ is the observed case-counts and
$y^{(99)}_{r,i}$ is the alarm boundary (the ${\rm 99^{th}}$ percentile
computed from $\boldsymbol{Y}$). $\gamma_{r,i}>1$ would denote an outlier. Let

$\overline{\gamma_{r}}=\frac{1}{N_{smooth}}\sum_{i=1}^{N_{smooth}}\gamma_{r,i}$

where $N_{smooth}$ is a time-period over which we will average out the
temporal noise and detect spatial patterns. Fig. 12 (top row) shows a map of
$\overline{\gamma_{r}}$ computed over a two-week period after August ${\rm
15^{th}}$ (on the left) and a two-week period after September ${\rm 15^{th}}$
(on the right) when the Fall 2020 wave had arrived. We see from the colormap
that the $\overline{\gamma_{r}}\leq 1$ before the arrival in all counties
except one; specifically, the false positives seen in Fig. 10 for Santa Fe and
Valencia are no longer visible. After the arrival, we see clear spatial
structures on the right. Thus simply averaging the data in time provides a
more robust detection, though with a loss of timeliness. However, we see two
adjacent counties (in shades of yellow) with very high
$\overline{\gamma_{r}}$.

Epidemiological activity spreads due to population mixing, and as we found in
our Part I paper, this is largely between adjoining NM counties (because of
their large spatial expanse). Therefore adjacent counties are expected to have
similar epidemiological behavior. If the epidemiological activity is
represented using a derived / inferred / estimated quantity, and adjoining
areal units differ greatly, it is either a consequence of erroneous estimation
or an artifact of noise in data i.e., spatial autocorrelation might be a
simple mechanism for removing erratic behavior in data. To this end, we
subject the $\overline{\gamma_{r}}$ map in Fig. 12 (top right) to clustering.
The centroids of the counties serve as the spatial feature whereas
$\overline{\gamma_{r}}$ serves as a measure of how strong the Fall 2020 wave
is, in each county. The data was Z-scored and subjected to hierarchical
clustering using the R Statistical Software[53] (R version 4.3.2 (2023-10-31))
package stats, specifically the function hclust(). The resulting dendrogram is
in Fig. 12 (bottom left), which is cut at a height of 0.6 (corresponding to a
quantile of 0.15) to reveal 3 clusters of counties. These are plotted in Fig.
12 (bottom right) and reveal the same clustering as was seen in Fig. 12 (top
right). However, the counties in yellow were eliminated - their level of
$\overline{\gamma_{r}}$, though high, were too different and thus violated the
need for spatial autocorrelation. Two of the clusters surround the county of
Bernalillo, the population center of NM, where COVID-19 was very active in the
city of Albuquerque.

##  6 Conclusion

In this paper, we have developed a scalable but approximate method to estimate
the infection-rate field of an outbreak, defined over a collection of areal
units. The method was demonstrated using COVID-19 data from the counties of
NM. The purpose of estimating the infection-rate field was to detect a sudden
change in the epidemiological dynamics, corresponding to the arrival of a new
wave of infections. Contemporary methods use case-counts to perform this
detection and are plagued by stochasticity and reporting errors in the data,
especially when case-counts are low (as may be expected from sparsely-
populated regions); in contrast, the infection-rate is governed by human
mixing patterns that do not vary erratically day-to-day. Our method is based
on mean-field variational inference (MFVI), but required some innovations for
computational efficiency and for enforcing non-negativity constraints in the
variables being estimated. The MFVI builds on a disease modeling formalism
that was developed in the Part I paper[7, 8], which also addressed the issue
of spatial autocorrelation in the COVID-19 dataset with a Gaussian Markov
Random field. This previous paper used Adaptive Markov chain Monte Carlo
(AMCMC) as the estimation algorithm and could not scale to all the counties of
NM.

The MFVI method obtains its scalability (in terms of the number of variables
being estimated) by imputing a parametrized form for the posterior density, a
set of independent Gaussians in our case; thereafter we estimate their means
and standard deviations using scalable gradient-based algorithms. We compared
the predictive skill of disease models that used the infection-rate field
estimated using AMCMC versus MFVI to answer the research questions posed in §
1. Our findings are listed below.

  * •

We find that the predictive skills of AMCMC- and MFVI-calibrated
epidemiological models are similar, even though the parameters’ posterior
distributions estimated by the two method are somewhat different. The
uncertainties in the parameters estimated from MFVI are too small to be
credible, in line with previous findings[10]. Both the approaches also
estimate the noise in the data (i.e., the component of data variability that
cannot be represented by the disease model), and the MFVI estimate is far
larger than the AMCMC counterpart. In this manner MFVI compensates for its
spuriously low uncertainty estimates in the parameters of the disease model,
and achieves similar levels of predictive skill as the AMCMC method (which
makes no approximations).

  * •

The GMRF spatial model, developed in the Part I paper[7, 8], plays an
important role in stabilizing the inversion. It regularizes the estimation
problem with a Gaussian Markov Random Field, and is sufficient to stabilize
the estimation when the observed case-counts from certain areal units
(counties) bear no resemblance to the waxing and waning of the COVID-19
pandemic that is clearly observed in case-count data aggregated to the state
level. Counties like Cibola and Rio Arriba, where errors / shortcomings in the
data do not allow a reliable estimation of the infection-rate are regularised
by it, and in the case of Rio Arriba also provides a credible forecast of the
epidemiological dynamics. This robustness to occasional low quality data,
which is to be expected for large inversions, allows the automation of the
estimation process without any manual “cleaning” of the data and imputations
of “cleaned” data values.

  * •

The infection-rate field, estimated using MFVI, is used to detect the arrival
of the Fall 2020 wave of COVID-19 in NM. To do so, we design a crude temporal
anomaly detector which contrasts model forecasts with observed data; large
discrepancies imply a change in epidemiological dynamics from the past. We had
no difficulty in detecting the wave when it was present, but suffered from
false positives when tests were conducted using data from before the arrival.
These false positives were caused by high-variance noise in low case-count
data and could be removed by temporal averaging. However, by doing so, our
detections were no longer very timely. In addition, the availability of an
infection-rate _field_ allowed us to exploit the spatial autocorrelation to
remove counties with spuriously high levels of epidemiological activity when
their neighbors were quiescent. It is clear that MFVI yields useful
information about the outbreak signature and does so in a scalable fashion.

Note that we made no attempt to design a proper anomaly detector with tunable
parameters to trade-off specificity versus sensitivity and plot Receiver
Operation Characteristic curves; our aim was to merely test of the existence
of certain information in the infection-rate estimates. However, the study
revealed that a proper anomaly detector would either have to smooth out high-
variance noise in low case-count data, or be formulated using Negative
Binomial or Poisson assumptions; these will be necessary to suppress false
positives when the data has high-variance noise. False negatives, on the other
hand, may not require any special considerations; this is certainly true for
the NM COVID-19 dataset.

## Author contributions

Wyatt Bridgman formulated the problem, wrote the software to solve it,
generated the figures and wrote the paper. Cosmin Safta assisted with
variational inference, interpretation of results, and contributed to writing
the paper. Jaideep Ray posed the problem, assisted with the epidemiological
interpretation, and suggested the calibration approach and metrics.

## Acknowledgments

This paper (SAND2024-08888O) describes objective technical results and
analysis. Any subjective views or opinions that might be expressed in the
paper do not necessarily represent the views of the U.S. Department of Energy
or the United States Government. This article has been authored by an employee
of National Technology & Engineering Solutions of Sandia, LLC under Contract
No. DE-NA0003525 with the U.S. Department of Energy (DOE). The employee owns
all right, title and interest in and to the article and is solely responsible
for its contents. The United States Government retains and the publisher, by
accepting the article for publication, acknowledges that the United States
Government retains a non-exclusive, paid-up, irrevocable, world-wide license
to publish or reproduce the published form of this article or allow others to
do so, for United States Government purposes. The DOE will provide public
access to these results of federally sponsored research in accordance with the
DOE Public Access Plan https://www.energy.gov/downloads/doe-public-access-
plan.

## Financial disclosure

None reported.

## Conflict of interest

The authors declare no potential conflict of interests.

## References

  * [1] M. L. Daza-Torres, M. A. Capistrán, A. Capella, and J. A. Christen, “Bayesian sequential data assimilation for covid-19 forecasting,” Epidemics, vol. 39, p. 100564, 2022.
  * [2] Z. Wang, X. Zhang, G. Teichert, M. Carrasco-Teja, and K. Garikipati, “System inference for the spatio-temporal evolution of infectious diseases: Michigan in the time of covid-19,” Computational Mechanics, vol. 66, pp. 1153–1176, 2020.
  * [3] P. Chen, K. Wu, and O. Ghattas, “Bayesian inference of heterogeneous epidemic models: Application to covid-19 spread accounting for long-term care facilities,” Computer Methods in Applied Mechanics and Engineering, vol. 385, p. 114020, 2021.
  * [4] P. Blonigan, J. Ray, and C. Safta, “Forecasting multi-wave epidemics through bayesian inference,” Archives of Computational Methods in Engineering, vol. 28, no. 6, pp. 4169–4183, 2021.
  * [5] Y. Lin, J. Neumann, E. Miller, R. Posner, A. Mallela, C. Safta, J. Ray, G. Thakur, S. Chinthavali, and W. Hlavacek, “Daily forecasting of regional epidemics of coronavirus disease with bayesian uncertainty quantification, united states.,” Emerging Infectious Diseases, vol. 27, no. 3, pp. 767–778, 2021.
  * [6] C. Safta, J. Ray, and K. Sargsyan, “Characterization of partially observed epidemics through bayesian inference: Application to covid-19,” Computational Mechanics, vol. 66, no. 5, pp. 1109–1129, 2020.
  * [7] C. Safta, J. Ray, and W. Bridgman, “Detecting outbreaks using a latent field: Part I - Spatial modeling,” Statistics in Medicine, 2024.  Submitted.
  * [8] C. Safta, W. Bridgman, and J. Ray, “Detecting outbreaks using a latent field: Part i – spatial modeling.” https://arxiv.org/abs/2406.12810, 2024.
  * [9] H. Haario, E. Saksman, and J. Tamminen, “An adaptive metropolis algorithm,” Bernoulli, pp. 223–242, 2001.
  * [10] W. Han and Y. Yang, “Statistical inference in mean-field variational bayes,” arXiv preprint arXiv:1911.01525, 2019.
  * [11] X. Huang, H. Zhou, X. Yang, W. Zhou, J. Huang, and Y. Yuan, “Spatial characteristics of coronavirus disease 2019 and their possible relationship with environmental and meteorological factors in hubei province, china,” GeoHealth, vol. 5, no. 6, p. e2020GH000358, 2021.
  * [12] X. Geng, G. G. Katul, F. Gerges, E. Bou-Zeid, H. Nassif, and M. C. Boufadel, “A kernel-modulated sir model for covid-19 contagiousspread from county to continent,” Proceedings of the National Academy of Sciences, vol. 118, no. 21, p. e2023321118, 2021.
  * [13] L. Schuler, J. M. Calabrese, and S. Attinger, “Data driven high resolution modeling and spatial analyses of the covid-19 pandemic in germany,” PLOS ONE, vol. 16, pp. 1–14, 08 2021.
  * [14] T. McMahon, A. Chan, S. Havlin, and L. K. Gallos, “Spatial correlations in geographical spreading of COVID-19 in the United States,” Scientific Reports, vol. 12, no. 1, p. 699, 2022.
  * [15] S. H. S. Indika, N. Diawara, H. A. Jeng, B. D. Giles, and D. S. K. Gamage, “Modeling the spread of covid-19 in spatio-temporal context,” Mathematical Biosciences and Engineering, vol. 20, no. 6, pp. 10552–10569, 2023\.
  * [16] A. Lawson and D. Lee, “Chapter 16 - bayesian disease mapping for public health,” in Disease Modelling and Public Health, Part A (A. S. Srinivasa Rao, S. Pyne, and C. Rao, eds.), vol. 36 of Handbook of Statistics, pp. 443–481, Elsevier, 2017.
  * [17] J. Besag, J. York, and A. Mollié, “Bayesian image restoration, with two applications in spatial statistics,” Annals of the institute of statistical mathematics, vol. 43, pp. 1–20, 1991.
  * [18] A. B. Lawson and H.-R. Song, “Bayesian hierarchical modeling of the dynamics of spatio-temporal influenza season outbreaks,” Spatial and spatio-temporal epidemiology, vol. 1, no. 2-3, pp. 187–195, 2010.
  * [19] A. B. Lawson, “Evaluation of predictive capability of bayesian spatio-temporal models for covid-19 spread,” BMC Medical Research Methodology, vol. 23, no. 1, p. 182, 2023.
  * [20] J. Kim, A. B. Lawson, B. Neelon, J. E. Korte, J. M. Eberth, and G. Chowell, “Evaluation of bayesian spatiotemporal infectious disease models for prospective surveillance analysis,” BMC Medical Research Methodology, vol. 23, no. 1, p. 171, 2023.
  * [21] A. B. Lawson and J. Kim, “Space-time covid-19 bayesian sir modeling in south carolina,” PLoS One, vol. 16, no. 3, p. e0242777, 2021.
  * [22] B. Sartorius, A. Lawson, and R. Pullan, “Modelling and predicting the spatio-temporal spread of covid-19, associated deaths and impact of key risk factors in england,” Scientific reports, vol. 11, no. 1, p. 5378, 2021\.
  * [23] C. Rotejanaprasert, A. Lawson, S. Bolick-Aldrich, and D. Hurley, “Spatial bayesian surveillance for small area case event data,” Statistical methods in medical research, vol. 25, no. 4, pp. 1101–1117, 2016.
  * [24] A. Corberán-Vallet and A. B. Lawson, “Conditional predictive inference for online surveillance of spatial disease incidence,” Statistics in medicine, vol. 30, no. 26, pp. 3095–3116, 2011.
  * [25] S. Brooks, “Markov chain monte carlo method and its application,” Journal of the royal statistical society: series D (the Statistician), vol. 47, no. 1, pp. 69–100, 1998.
  * [26] D. Van Ravenzwaaij, P. Cassey, and S. D. Brown, “A simple introduction to markov chain monte–carlo sampling,” Psychonomic bulletin & review, vol. 25, no. 1, pp. 143–154, 2018.
  * [27] G. O. Roberts and R. L. Tweedie, “Exponential convergence of langevin distributions and their discrete approximations,” Bernoulli, pp. 341–363, 1996.
  * [28] D. M. Blei, A. Kucukelbir, and J. D. McAuliffe, “Variational inference: A review for statisticians,” Journal of the American Statistical Association, vol. 112, no. 518, pp. 859–877, 2017.
  * [29] M. I. Jordan, Z. Ghahramani, T. S. Jaakkola, and L. K. Saul, “An introduction to variational methods for graphical models,” Machine Learning, vol. 37, no. 2, pp. 183–233, 1999.
  * [30] C. Blundell, J. Cornebise, K. Kavukcuoglu, and D. Wierstra, “Weight uncertainty in neural networks,” in International conference on machine learning, pp. 1613–1622, PMLR, 2015.
  * [31] D. P. Kingma, M. Welling, et al., “An introduction to variational autoencoders,” Foundations and Trends® in Machine Learning, vol. 12, no. 4, pp. 307–392, 2019.
  * [32] J. M. Hernández-Lobato and R. Adams, “Probabilistic backpropagation for scalable learning of bayesian neural networks,” in International conference on machine learning, pp. 1861–1869, PMLR, 2015.
  * [33] R. Tanno, D. E. Worrall, A. Ghosh, E. Kaden, S. N. Sotiropoulos, A. Criminisi, and D. C. Alexander, “Bayesian image quality transfer with cnns: exploring uncertainty in dmri super-resolution,” in Medical Image Computing and Computer Assisted Intervention- MICCAI 2017: 20th International Conference, Quebec City, QC, Canada, September 11-13, 2017, Proceedings, Part I 20, pp. 611–619, Springer, 2017.
  * [34] Ł. Rączkowski, M. Możejko, J. Zambonelli, and E. Szczurek, “Ara: accurate, reliable and active histopathological image classification framework with bayesian deep learning,” Scientific reports, vol. 9, no. 1, p. 14347, 2019.
  * [35] O. Ozdemir, B. Woodward, and A. A. Berlin, “Propagating uncertainty in multi-stage bayesian convolutional neural networks with application to pulmonary nodule detection,” arXiv preprint arXiv:1712.00497, 2017.
  * [36] G. Luo, S. Dong, W. Wang, K. Wang, S. Cao, C. Tam, H. Zhang, J. Howey, P. Ohorodnyk, and S. Li, “Commensal correlation network between segmentation and direct area estimation for bi-ventricle quantification,” Medical image analysis, vol. 59, p. 101591, 2020.
  * [37] E. D. Carvalho, R. Clark, A. Nicastro, and P. H. Kelly, “Scalable uncertainty for computer vision with functional variational inference,” in Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pp. 12003–12013, 2020.
  * [38] P. Liang, S. Petrov, M. I. Jordan, and D. Klein, “The infinite pcfg using hierarchical dirichlet processes,” in Proceedings of the 2007 joint conference on empirical methods in natural language processing and computational natural language learning (EMNLP-CoNLL), pp. 688–697, 2007.
  * [39] S. Hu, X. Xie, S. Liu, J. Yu, Z. Ye, M. Geng, X. Liu, and H. Meng, “Bayesian learning of lf-mmi trained time delay neural networks for speech recognition,” IEEE/ACM Transactions on Audio, Speech, and Language Processing, vol. 29, pp. 1514–1529, 2021.
  * [40] X. Meng, H. Babaee, and G. E. Karniadakis, “Multi-fidelity bayesian neural networks: Algorithms and applications,” Journal of Computational Physics, vol. 438, p. 110361, 2021.
  * [41] L. Yang, X. Meng, and G. E. Karniadakis, “B-pinns: Bayesian physics-informed neural networks for forward and inverse pde problems with noisy data,” Journal of Computational Physics, vol. 425, p. 109913, 2021.
  * [42] R. Dandekar, K. Chung, V. Dixit, M. Tarek, A. Garcia-Valadez, K. V. Vemula, and C. Rackauckas, “Bayesian neural ordinary differential equations,” arXiv preprint arXiv:2012.07244, 2020.
  * [43] K. Friston, A. Costello, and D. Pillay, “‘dark matter’, second waves and epidemiological modelling,” BMJ Global Health, vol. 5, no. 12, p. e003978, 2020.
  * [44] K. J. Friston, G. Flandin, and A. Razi, “Dynamic causal modelling of covid-19 and its mitigations,” Scientific reports, vol. 12, no. 1, p. 12419, 2022\.
  * [45] Z. Wang, X. Zhang, G. Teichert, M. Carrasco-Teja, and K. Garikipati, “System inference for the spatio-temporal evolution of infectious diseases: Michigan in the time of covid-19,” Computational Mechanics, vol. 66, pp. 1153–1176, 2020.
  * [46] K. Fan, C. Li, and K. Heller, “A unifying variational inference framework for hierarchical graph-coupled hmm with an application to influenza infection,” in Proceedings of the AAAI Conference on Artificial Intelligence, vol. 30, 2016.
  * [47] I. Biazzo, A. Braunstein, L. Dall’Asta, and F. Mazza, “A bayesian generative neural network framework for epidemic inference problems,” Scientific Reports, vol. 12, no. 1, p. 19673, 2022.
  * [48] S. A. Lauer, K. H. Grantz, Q. Bi, F. K. Jones, Q. Zheng, H. R. Meredith, A. S. Azman, N. G. Reich, and J. Lessler, “The Incubation Period of Coronavirus Disease 2019 (COVID-19) From Publicly Reported Confirmed Cases: Estimation and Application,” Annals of Internal Medicine, 2020.
  * [49] A. Kucukelbir, D. Tran, R. Ranganath, A. Gelman, and D. M. Blei, “Automatic differentiation variational inference,” Journal of machine learning research, 2017.
  * [50] T. Gneiting and A. E. Raftery, “Strictly proper scoring rules, prediction, and estimation,” Journal of the American statistical Association, vol. 102, no. 477, pp. 359–378, 2007.
  * [51] C. Safta, D. M. Ricciuto, K. Sargsyan, B. Debusschere, H. N. Najm, M. Williams, and P. E. Thornton, “Global sensitivity analysis, probabilistic calibration, and predictive assessment for the data assimilation linked ecosystem carbon model,” Geoscientific Model Development, vol. 8, no. 7, pp. 1899–1918, 2015.
  * [52] A. Hegde, E. Weiss, W. Windl, H. N. Najm, and C. Safta, “A bayesian calibration framework with embedded model error for model diagnostics,” International Journal for Uncertainty Quantification, vol. 14, no. 6, pp. 37–70, 2024.
  * [53] R Core Team, R: A Language and Environment for Statistical Computing.  R Foundation for Statistical Computing, Vienna, Austria, 2023.
  * [54] M. Salmon, D. Schumacher, and M. Höhle, “Monitoring count time series in R: Aberration detection in public health surveillance,” Journal of Statistical Software, vol. 70, no. 10, pp. 1–35, 2016.
  * [55] “Coronavirus (Covid-19) Data in the United States.” https://github.com/nytimes/covid-19-data.  Accessed: 2023-01-01.
  * [56] “COVID-19 Data Repository by the Center for Systems Science and Engineering (CSSE) at Johns Hopkins University.” https://github.com/CSSEGISandData/COVID-19.  Accessed: 2023-01-01.
  * [57] D. P. Kingma and J. Ba, “Adam: A method for stochastic optimization,” arXiv preprint arXiv:1412.6980, 2014.
  * [58] M. Höhle and M. Paul, “Count data regression charts for the monitoring of surveillance time series,” Computational Statistics and Data Analysis, vol. 52, no. 9, pp. 4357–4368, 2008.
  * [59] K. B. Petersen, M. S. Pedersen, et al., “The matrix cookbook,” Technical University of Denmark, vol. 7, no. 15, p. 510, 2008.

##  Appendix A Variational Inference

###  A.1 Score gradients of the ELBO

We briefly review the score estimator, or black-box, approach to estimating
the gradient of the ELBO. Recall that the ELBO is given by 10 which, for the
sake of clarity, can be written in a more generic form

$\mathcal{L}(\boldsymbol{\phi})=\mathbb{E}_{q(\boldsymbol{\theta};\boldsymbol{\phi})}\left[f(\boldsymbol{\theta})\right]$
(15)

where $f(\boldsymbol{\theta})$ encapsulates the dependence on the random
vector $\boldsymbol{\theta}$. To derive an estimator for the gradient, we can
carry out the following manipulations

$\displaystyle\nabla_{\boldsymbol{\phi}}\mathcal{L}(\boldsymbol{\phi})$
$\displaystyle=\nabla_{\boldsymbol{\phi}}\mathbb{E}_{q(\boldsymbol{\theta};\boldsymbol{\phi})}\left[f(\boldsymbol{\theta})\right]$
(16) $\displaystyle=\nabla_{\boldsymbol{\phi}}\int
q(\boldsymbol{\theta};\boldsymbol{\phi})f(\boldsymbol{\theta})d\,\boldsymbol{\theta}$
(17)
$\displaystyle=\int\nabla_{\boldsymbol{\phi}}q(\boldsymbol{\theta};\boldsymbol{\phi})f(\boldsymbol{\theta})d\,\boldsymbol{\theta}$
(18) $\displaystyle=\int
q(\boldsymbol{\theta};\boldsymbol{\phi})\frac{\nabla_{\boldsymbol{\phi}}q(\boldsymbol{\theta};\boldsymbol{\phi})}{q(\boldsymbol{\theta};\boldsymbol{\phi})}f(\boldsymbol{\theta})d\,\boldsymbol{\theta}$
(19)
$\displaystyle=\mathbb{E}_{q(\boldsymbol{\theta};\boldsymbol{\phi})}\left[f(\boldsymbol{\theta})\nabla_{\boldsymbol{\phi}}\log
q(\boldsymbol{\theta};\boldsymbol{\phi})\right]$ (20)

Hence, the gradient can be expressed as an expectation with respect to
$q(\boldsymbol{\theta};\boldsymbol{\phi})$ where only the log of the surrogate
posterior needs to be differentiated with respect to the variational
parameters $\boldsymbol{\phi}$.

###  A.2 Reparametrization gradients of the ELBO

The likelihood and log likelihood are given by

$\displaystyle p(\mathcal{D}|\boldsymbol{\theta})$
$\displaystyle=\prod_{i=1}^{N_{d}}2\pi^{-N_{r}/2}\det(\mathbf{\Sigma}_{i})^{-1/2}\exp\left(-\frac{1}{2}(\mathbf{y}_{i}^{(o)}-\mathbf{y}_{i})^{T}\mathbf{\Sigma}_{i}^{-1}(\mathbf{y}_{i}^{(o)}-\mathbf{y}_{i})\right)$
(21) $\displaystyle l(\boldsymbol{\theta})$
$\displaystyle=-\frac{N_{d}N_{r}2\pi}{2}-\frac{1}{2}\sum_{i=1}^{N_{d}}\log\det(\mathbf{\Sigma}_{i})+(\mathbf{y}_{i}^{(o)}-\mathbf{y}_{i})^{T}\mathbf{\Sigma}_{i}^{-1}(\mathbf{y}_{i}^{(o)}-\mathbf{y}_{i})$
(22)

Using the reparametrization trick, we can write the ELBO (10) and its gradient
in the form

$\displaystyle\mathcal{L}(\boldsymbol{\phi})$
$\displaystyle=-\mathbb{H}[q(\boldsymbol{\theta};\boldsymbol{\phi})]-\mathbb{E}_{q(\boldsymbol{\epsilon})}[\log
p(\mathcal{D}|\boldsymbol{\theta}(\boldsymbol{\epsilon},\boldsymbol{\phi}))+\log
p(\boldsymbol{\theta}(\boldsymbol{\epsilon},\boldsymbol{\phi}))]$ (23)
$\displaystyle\nabla_{\boldsymbol{\phi}}\mathcal{L}(\boldsymbol{\phi})$
$\displaystyle=-\nabla_{\boldsymbol{\phi}}\mathbb{H}[q(\boldsymbol{\theta};\boldsymbol{\phi})]-\mathbb{E}_{q(\boldsymbol{\epsilon})}[\nabla_{\boldsymbol{\phi}}\log
p(\mathcal{D}|\boldsymbol{\theta}(\boldsymbol{\epsilon},\boldsymbol{\phi}))+\nabla_{\boldsymbol{\phi}}\log
p(\boldsymbol{\theta}(\boldsymbol{\epsilon},\boldsymbol{\phi}))]$ (24)

where $\boldsymbol{\phi}=(\boldsymbol{\mu},\boldsymbol{\rho})$,
$\boldsymbol{\theta}=\boldsymbol{\mu}+\boldsymbol{\sigma}(\boldsymbol{\rho})\odot\boldsymbol{\epsilon}$
with $\boldsymbol{\epsilon}\sim\mathcal{N}(\mathbf{0},\mathbf{I})$. Here,
$\sigma$ is a positive transformation of the unconstrained variable
$\boldsymbol{\rho}$ to ensure the variance is constrained to be positive. A
Monte Carlo estimator of the gradient can then be written as

$\displaystyle\nabla_{\boldsymbol{\phi}}\mathcal{L}(\boldsymbol{\phi})$
$\displaystyle\approx-\nabla_{\boldsymbol{\phi}}\mathbb{H}[q(\boldsymbol{\theta};\boldsymbol{\phi})]-\frac{1}{N_{s}}\sum_{i=1}^{N_{s}}\nabla_{\boldsymbol{\phi}}\log
p(\mathcal{D}|\boldsymbol{\theta}(\boldsymbol{\epsilon}_{i},\boldsymbol{\phi}))+\nabla_{\boldsymbol{\phi}}\log
p(\boldsymbol{\theta}(\boldsymbol{\epsilon}_{i},\boldsymbol{\phi}))$ (25)
$\displaystyle=-\nabla_{\boldsymbol{\phi}}\mathbb{H}[q(\boldsymbol{\theta};\boldsymbol{\phi})]-\frac{1}{N_{s}}\sum_{i=1}^{N_{s}}\left(\nabla_{\boldsymbol{\theta}}\log
p(\mathcal{D}|\boldsymbol{\theta}(\boldsymbol{\epsilon}_{i},\boldsymbol{\phi}))+\nabla_{\boldsymbol{\theta}}\log
p(\boldsymbol{\theta}(\boldsymbol{\epsilon}_{i},\boldsymbol{\phi}))\right)\odot\nabla_{\boldsymbol{\phi}}\boldsymbol{\theta}(\boldsymbol{\epsilon}_{i},\boldsymbol{\phi})$
(26)

where the last line is given by the chain rule and the fact that
$\boldsymbol{\theta}$ is defined by an element-wise transformation of
$\boldsymbol{\phi}$. Observe that

$\displaystyle\nabla_{\boldsymbol{\mu}}\boldsymbol{\theta}$
$\displaystyle=\mathbf{1}$ (27)
$\displaystyle\nabla_{\boldsymbol{\rho}}\boldsymbol{\theta}$
$\displaystyle=\nabla_{\boldsymbol{\rho}}\boldsymbol{\sigma}(\boldsymbol{\rho})\odot\boldsymbol{\epsilon}_{i}$
(28)
$\displaystyle\nabla_{\boldsymbol{\mu}}\mathbb{H}[q(\boldsymbol{\theta};\boldsymbol{\phi})]$
$\displaystyle=\mathbf{0}$ (29)
$\displaystyle\nabla_{\boldsymbol{\rho}}\mathbb{H}[q(\boldsymbol{\theta};\boldsymbol{\phi})]$
$\displaystyle=\frac{1}{\boldsymbol{\sigma}(\boldsymbol{\rho})}\odot\nabla_{\boldsymbol{\rho}}\boldsymbol{\sigma}(\boldsymbol{\rho})$
(30)

so that it remains to compute the gradients of the log-likelihood and prior.

###  A.3 Gradients of the Log Likelihood

As the log likelihood factors independently across the data
$i=1,\ldots,N_{d}$, it suffices to compute the gradients
$\nabla_{\boldsymbol{\theta}}\log\det\mathbf{\Sigma}_{i}$ and
$\nabla_{\boldsymbol{\theta}}(\mathbf{y}_{i}^{(o)}-\mathbf{y}_{i})^{T}\mathbf{\Sigma}_{i}^{-1}(\mathbf{y}_{i}^{(o)}-\mathbf{y}_{i})$
for a particular day $i$. The differentials of the these two terms are
computed using matrix calculus [59] as

$\displaystyle\partial\log\det\mathbf{\Sigma}_{i}$
$\displaystyle=\text{Tr}(\mathbf{\Sigma}_{i}^{-1}\partial\mathbf{\Sigma}_{i})$
(31)
$\displaystyle\partial(\mathbf{y}_{i}^{(o)}-\mathbf{y}_{i})^{T}\mathbf{\Sigma}_{i}^{-1}(\mathbf{y}_{i}^{(o)}-\mathbf{y}_{i})$
$\displaystyle=-(\mathbf{y}_{i}^{(o)}-\mathbf{y}_{i})^{T}\mathbf{\Sigma}_{i}^{-1}(\partial\mathbf{\Sigma}_{i})\mathbf{\Sigma}_{i}^{-1}(\mathbf{y}_{i}^{(o)}-\mathbf{y}_{i})-2(\mathbf{y}_{i}^{(o)}-\mathbf{y}_{i})^{T}\mathbf{\Sigma}_{i}^{-1}\partial\mathbf{y}_{i}$
(32)

where the differential of $\mathbf{\Sigma}_{i}$ is

$\displaystyle\partial\mathbf{\Sigma}_{i}$
$\displaystyle=(\partial\tau_{\Phi})[\mathbf{D}-\lambda_{\Phi}\mathbf{W}]^{-1}-(\partial\lambda_{\Phi})\tau_{\Phi}[\mathbf{D}-\lambda_{\Phi}\mathbf{W}]^{-1}\mathbf{W}[\mathbf{D}-\lambda_{\Phi}\mathbf{W}]^{-1}$
(33)
$\displaystyle+2\text{diag}(\sigma_{a}+\sigma_{m}\mathbf{y}_{i})[(\partial\sigma_{a})\mathbf{I}+(\partial\sigma_{m})\text{diag}(\mathbf{y}_{i})+\sigma_{m}\text{diag}(\partial\mathbf{y}_{i})]$
(34)

We have the following derivatives

  * •

Case:
$\hat{\mathbf{m}}_{r}^{T}=(\hat{t}_{0}^{r},\hat{N}^{r},\hat{k}^{r},\hat{\theta}^{r})$
are the unconstrained model variables for region $r$ and $\hat{\theta}_{r}$ is
a particular variable.

$\displaystyle\partial_{\hat{\theta}_{r}}\mathbf{y}_{i}$
$\displaystyle=(0,\ldots,\partial_{\theta_{r}}y_{r}(i;t_{0}^{r},N_{r},k^{r},\theta^{r})\theta^{\prime}_{i}(\hat{\theta}_{i}),\ldots,0)$
(35) $\displaystyle\nabla_{\hat{\mathbf{m}}_{r}}\log\det\mathbf{\Sigma}_{i}$
$\displaystyle=2\sigma_{m}[\mathbf{\Sigma}_{i}^{-1}]_{ii}\nabla_{\hat{\mathbf{m}}_{r}}y_{r}(i;t_{0}^{r},N^{r},k^{r},\theta^{r})$
(36)
$\displaystyle\partial_{\hat{\theta}_{r}}(\mathbf{y}_{i}^{(o)}-\mathbf{y}_{i})^{T}\mathbf{\Sigma}_{i}^{-1}(\mathbf{y}_{i}^{(o)}-\mathbf{y}_{i})$
$\displaystyle=-2\sigma_{m}(\mathbf{y}_{i}^{(o)}-\mathbf{y}_{i})^{T}\mathbf{\Sigma}_{i}^{-1}\text{diag}(\sigma_{a}+\sigma_{m}\mathbf{y}_{i})\text{diag}(\partial_{\hat{\theta}_{r}}\mathbf{y}_{i})\mathbf{\Sigma}_{i}^{-1}(\mathbf{y}_{i}^{(o)}-\mathbf{y}_{i})$
(37)
$\displaystyle-2(\mathbf{y}_{i}^{(o)}-\mathbf{y}_{i})^{T}\mathbf{\Sigma}_{i}^{-1}\partial_{\hat{\theta}_{r}}\mathbf{y}_{i}$
(38)

  * •

Case: $\hat{\theta}_{r}=\hat{\tau}_{\Phi}$

$\displaystyle\partial_{\hat{\theta}_{r}}\mathbf{y}_{i}$
$\displaystyle=\mathbf{0}$ (39)
$\displaystyle\partial_{\hat{\theta}_{r}}\log\det\mathbf{\Sigma}_{i}$
$\displaystyle=\tau^{\prime}_{\Phi}(\hat{\tau}_{\Phi})\text{Tr}(\mathbf{\Sigma}_{i}^{-1}[\mathbf{I}-\lambda_{\Phi}\mathbf{W}]^{-1})$
(40)
$\displaystyle\partial_{\hat{\theta}_{r}}(\mathbf{y}_{i}^{(o)}-\mathbf{y}_{i})^{T}\mathbf{\Sigma}_{i}^{-1}(\mathbf{y}_{i}^{(o)}-\mathbf{y}_{i})$
$\displaystyle=\tau^{\prime}_{\Phi}(\hat{\tau}_{\Phi})(\mathbf{y}_{i}^{(o)}-\mathbf{y}_{i})^{T}\mathbf{\Sigma}_{i}^{-1}[\mathbf{I}-\lambda_{\Phi}\mathbf{W}]^{-1}\mathbf{\Sigma}_{i}^{-1}(\mathbf{y}_{i}^{(o)}-\mathbf{y}_{i})$
(41)

  * •

Case: $\hat{\theta}_{r}=\hat{\lambda}_{\Phi}$

$\displaystyle\partial_{\hat{\theta}_{r}}\mathbf{y}_{i}$
$\displaystyle=\mathbf{0}$ (42)
$\displaystyle\partial_{\hat{\theta}_{r}}\log\det\mathbf{\Sigma}_{i}$
$\displaystyle=\lambda^{\prime}_{\Phi}(\hat{\lambda}_{\Phi})\tau_{\Phi}\text{Tr}(\mathbf{\Sigma}_{i}^{-1}[\mathbf{I}-\lambda_{\Phi}\mathbf{W}]^{-1}\mathbf{W}[\mathbf{I}-\lambda_{\Phi}\mathbf{W}]^{-1})$
(43)
$\displaystyle\partial_{\hat{\theta}_{r}}(\mathbf{y}_{i}^{(o)}-\mathbf{y}_{i})^{T}\mathbf{\Sigma}_{i}^{-1}(\mathbf{y}_{i}^{(o)}-\mathbf{y}_{i})$
$\displaystyle=\lambda^{\prime}_{\Phi}(\hat{\lambda}_{\Phi})\tau_{\Phi}(\mathbf{y}_{i}^{(o)}-\mathbf{y}_{i})^{T}\mathbf{\Sigma}_{i}^{-1}[\mathbf{I}-\lambda_{\Phi}\mathbf{W}]^{-1}\mathbf{W}[\mathbf{I}-\lambda_{\Phi}\mathbf{W}]^{-1}\mathbf{\Sigma}_{i}^{-1}(\mathbf{y}_{i}^{(o)}-\mathbf{y}_{i})$
(44)

  * •

Case: $\hat{\theta}_{r}=\hat{\sigma}_{a}$

$\displaystyle\partial_{\hat{\theta}_{r}}\mathbf{y}_{i}$
$\displaystyle=\mathbf{0}$ (45)
$\displaystyle\partial_{\hat{\theta}_{r}}\log\det\mathbf{\Sigma}_{i}$
$\displaystyle=2\sigma^{\prime}_{a}(\hat{\sigma}_{a})\text{Tr}(\mathbf{\Sigma}_{i}^{-1}\text{diag}(\sigma_{a}+\sigma_{m}\mathbf{y}_{i}))$
(46)
$\displaystyle\partial_{\hat{\theta}_{r}}(\mathbf{y}_{i}^{(o)}-\mathbf{y}_{i})^{T}\mathbf{\Sigma}_{i}^{-1}(\mathbf{y}_{i}^{(o)}-\mathbf{y}_{i})$
$\displaystyle=-2\sigma^{\prime}_{a}(\hat{\sigma}_{a})(\mathbf{y}_{i}^{(o)}-\mathbf{y}_{i})^{T}\mathbf{\Sigma}_{i}^{-1}\text{diag}(\sigma_{a}+\sigma_{m}\mathbf{y}_{i})\mathbf{\Sigma}_{i}^{-1}(\mathbf{y}_{i}^{(o)}-\mathbf{y}_{i})$
(47)

  * •

Case:$\hat{\theta}_{r}=\hat{\sigma}_{m}$

$\displaystyle\partial_{\hat{\theta}_{r}}\mathbf{y}_{i}$
$\displaystyle=\mathbf{0}$ (48)
$\displaystyle\partial_{\hat{\theta}_{r}}\log\det\mathbf{\Sigma}_{i}$
$\displaystyle=2\sigma^{\prime}_{m}(\hat{\sigma}_{m})\text{Tr}(\mathbf{\Sigma}_{i}^{-1}\text{diag}((\sigma_{a}+\sigma_{m}\mathbf{y}_{i})\odot\mathbf{y}_{i}))$
(49)
$\displaystyle\partial_{\hat{\theta}_{r}}(\mathbf{y}_{i}^{(o)}-\mathbf{y}_{i})^{T}\mathbf{\Sigma}_{i}^{-1}(\mathbf{y}_{i}^{(o)}-\mathbf{y}_{i})$
$\displaystyle=-2\sigma^{\prime}_{m}(\hat{\sigma}_{m})(\mathbf{y}_{i}^{(o)}-\mathbf{y}_{i})^{T}\mathbf{\Sigma}_{i}^{-1}\text{diag}((\sigma_{a}+\sigma_{m}\mathbf{y}_{i})\odot\mathbf{y}_{i}))\mathbf{\Sigma}_{i}^{-1}(\mathbf{y}_{i}^{(o)}-\mathbf{y}_{i})$
(50)

The variable transformations are listed in Table 2.

Variable | Constraint | Transformation
---|---|---
$t_{0}^{r}$ | None | $id(\hat{t}_{0}^{r})$
$N^{r}$ | $N^{r}\geq 0$ | $\exp(\hat{N}^{r})$
$k^{r}$ | $k^{r}\geq 2$ | $s(\hat{k}^{r})+2$
$\theta^{r}$ | $k^{r}\geq\epsilon$ | $s(\hat{\theta}^{r})+\epsilon$
$\tau_{\Phi}$ | $\tau_{\Phi}\geq 0$ | $\exp(\hat{\tau}_{\Phi})$
$\lambda_{\Phi}$ | $0\leq\lambda_{\Phi}\leq 1$ | $(1+\exp(-\hat{\lambda}_{\Phi}))^{-1}$
$\sigma_{a}$ | $\sigma_{a}\geq 0$ | $\exp(\hat{\sigma}_{a})$
$\sigma_{m}$ | $\sigma_{m}\geq 0$ | $\exp(\hat{\sigma}_{m})$
Table 2: Variable transformations. $s(x):\log(\exp(x)+1)$ is the softplus
function and satisfies $s(x)\in[0,1]$.

###  A.4 Approximation of model predictions and gradients via quadrature

The model predictions $\mathbf{y}_{i}$, given by (3), involve a convolution
integral that cannot be expressed in closed form. Hence, we approximate the
predictions by integral quadrature

$\displaystyle y_{r}(i;t_{0}^{r},N^{r},k^{r},\theta^{r})$
$\displaystyle=\int_{t_{0}^{r}}^{t_{i}}N_{r}f_{inf}(\tau-
t_{0}^{r};k^{r},\theta^{r})(F_{inc}(t_{i}-\tau;\mu,\sigma)-(F_{inc}(t_{i-1}-\tau;\mu,\sigma))d\,\tau$
(51)
$\displaystyle\approx\sum_{j=1}^{n}N_{r}w_{j}f_{inf}(\tau_{j}-t_{0}^{r};k^{r},\theta^{r})(F_{inc}(t_{i}-\tau_{j};\mu,\sigma)-(F_{inc}(t_{i-1}-\tau_{j};\mu,\sigma))$
(52)

where $w_{j}$, $\tau_{j}$ are quadrature weights and points given by a method
such as Gaussian quadrature. As the function
$(F_{inc}(t_{i}-\tau;\mu,\sigma)-(F_{inc}(t_{i-1}-\tau;\mu,\sigma))$ does not
depend on parameters $\boldsymbol{\theta}$, we can write it as
$\tilde{F}_{inc}(\tau)$ for simplicity of notation. By the Leibniz integral
rule, we can write the derivatives of the model predictions as

$\displaystyle\partial_{\hat{\theta}_{r}}\mathbf{y}_{i}$
$\displaystyle=\int_{t_{0}^{r}}^{t_{i}}[\partial_{\hat{\theta}_{r}}N_{r}f_{inf}(\tau-
t_{0}^{r};k^{r},\theta^{r})]\tilde{F}_{inc}(\tau)d\,\tau$ if
$\hat{\theta}_{r}\neq t_{0}^{r}$ (53)
$\displaystyle\partial_{\hat{\theta}_{r}}\mathbf{y}_{i}$
$\displaystyle=\int_{t_{0}^{r}}^{t_{i}}[\partial_{\hat{\theta}_{r}}N_{r}f_{inf}(\tau-
t_{0}^{r};k^{r},\theta^{r})]\tilde{F}_{inc}(\tau)d\,\tau-
f_{inf}(0;k^{r},\theta^{r})$ if $\hat{\theta}_{r}=t_{0}^{r}$ (54)

This requires that the functions $N_{r}f_{inf}(\tau-
t_{0}^{r};k^{r},\theta^{r})]\tilde{F}_{inc}(\tau)$ and
$\partial_{\hat{\theta}_{r}}N_{r}f_{inf}(\tau-
t_{0}^{r};k^{r},\theta^{r})]\tilde{F}_{inc}(\tau)$ are continuous in $\tau$
and $\hat{\theta}_{r}$ in a region of the $\tau$-$\hat{\theta}_{r}$ plane
including $t_{0}\leq\tau\leq t_{i}$, a condition that’s easily met by ensuring
certain constraints are satisfied by the parameters
$(t_{0}^{r},N^{r},k^{r},\theta^{r})$ for $r=1,\ldots,N_{r}$. These constraints
are enforced via the variable transformations in (13). Hence, we can
approximate (53) and (54) via quadrature in the form

$\displaystyle\partial_{\hat{\theta}_{r}}\mathbf{y}_{i}$
$\displaystyle\approx\sum_{j=1}^{n}w_{j}[\partial_{\hat{\theta}_{r}}N_{r}f_{inf}(\tau_{j}-t_{0}^{r};k^{r},\theta^{r})]\tilde{F}_{inc}(\tau_{j})$
if $\hat{\theta}_{r}\neq t_{0}^{r}$ (55)
$\displaystyle\partial_{\hat{\theta}_{r}}\mathbf{y}_{i}$
$\displaystyle\approx\sum_{j=1}^{n}w_{j}[\partial_{\hat{\theta}_{r}}N_{r}f_{inf}(\tau_{j}-t_{0}^{r};k^{r},\theta^{r})]\tilde{F}_{inc}(\tau_{j})-f_{inf}(0;k^{r},\theta^{r})$
if $\hat{\theta}_{r}=t_{0}^{r}$ (56)

Hence, in this implementation of VI, gradients of the ELBO have a pseudo-
analytic form where “outer" gradients of the log likelihood with respect to
model predictions are exact and “inner" gradients of the model predictions
with respect to parameters are approximated via quadrature. This allows for
accurate gradient approximations that can be caculated efficiently leading to
a scalable VI algorithm that can be applied to the high-dimensional inverse
problem for the outbreak model.

##  Appendix B Infection-rate estimates

§ 4.1 describes the estimation of the infection-rate field over Bernalillo,
Santa Fe and Valencia using MFVI jointly (where the GMRF spatial model is
used) and individually. Fig. 2 plots the PPT runs, driven by a distribution of
infection-rate fields. The corresponding fields are in Fig. 13.

Figure 13: Comparison of the infection rate curves that determine predictions
in Fig. 2. These are for a 3-county inference of Bernalillo (left), Santa Fe
(middle), and Valencia (right) done jointly using the GMRF model (top) and
independently for each county (bottom). The red line shows the median
predictions, the shaded green line shows the inter-quartile range, and the
dashed lines are the $5^{\text{th}}$ and $95^{\text{th}}$ percentiles. Not
much of a difference can be discerned visually.

The same section describes the estimation of the infection-rate field over
Bernalillo, Santa Fe and Valencia using AMCMC and the MFVI. Fig. 3 plots the
PPT runs, driven by a distribution of infection-rate fields. The corresponding
fields are in Fig. 14.

Figure 14: Comparison of the infection rate curves that determine predictions
in Fig. 3. These are for a 3-county inference of Bernalillo (left), Santa Fe
(middle), and Valencia (right) done jointly using the MFVI (top) and AMCMC
(bottom), taken from our Part I paper[7, 8] The red line shows the median
predictions, the shaded green line shows the inter-quartile range, and the
dashed lines are the $5^{\text{th}}$ and $95^{\text{th}}$ percentiles. Not
much of a difference can be discerned visually.