MICHAL et al

*Alexandra M. Schmidt,

# A Bayesian hierarchical model for disease mapping that accounts for scaling
and heavy-tailed latent effects

Victoire Michal  Laís Picinini Freitas  Alexandra M. Schmidt  Department of
Epidemiology, Biostatistics and Occupational Health, McGill University,
Montreal, Canada  Programa de Computação Científica, Fundação Oswaldo Cruz,
Rio de Janeiro, Brazil<EMAIL_ADDRESS>Michal V  Freitas L.P
Schmidt A.M

(26 April 2016; 6 June 2016; 6 June 2016)

###### Abstract

[Abstract] In disease mapping, the relative risk of a disease is commonly
estimated across different areas within a region of interest. The number of
cases in an area is often assumed to follow a Poisson distribution whose mean
is decomposed as the product between an offset and the logarithm of the
disease’s relative risk. The log risk may be written as the sum of fixed
effects and latent random effects. The commonly used BYM model further
decomposes the latent effects into a sum of independent effects and spatial
effects to account for potential overdispersion and a spatial correlation
structure among the counts. However, this model suffers from an identifiably
issue. The BYM2 model reparametrises the latter by decomposing each latent
effect into a weighted sum of independent and spatial effects. We build on the
BYM2 model to allow for heavy-tailed latent effects and accommodate
potentially outlying risks, after accounting for the fixed effects. We assume
a scale mixture structure wherein the variance of the latent process changes
across areas and allows for outlier identification. We explore two prior
specifications of this scale mixture structure in simulation studies and in
the analysis of Zika cases from the 2015-2016 epidemic in Rio de Janeiro. The
simulation studies show that, in terms of WAIC and outlier detection, the two
parametrisations always perform well compared to commonly used models. Our
analysis of Zika cases finds 19 districts of Rio as potential outliers, after
accounting for the socio-development index, which may help prioritise
interventions.

, , (2021), A Bayesian hierarchical model for disease mapping that accounts
for scaling and heavy-tailed latent effects, Statistics in Medicine,
20XX;XX:XX–XX.

###### keywords:

BYM2 model, Outliers, Scale mixture, Spatial statistics, Vector-borne disease,
Zika virus infection

††articletype: Research Article

##  1 Motivation

The first Zika cases in the Americas were identified in 2015, when it was
considered a benign disease. However, in October 2015 an unprecedented
increase in the number of microcephaly cases in neonates was reported in the
Northeast of Brazil and was later associated with the Zika virus infection
during pregnancy 1. The Zika virus is transmitted to humans by the bite of
infected Aedes mosquitoes, the same vectors that transmit dengue, chikungunya
and yellow fever. Dengue is the most prevalent Aedes-borne disease in the
world and around 3.9 billion people in 129 countries are at risk of acquiring
the disease 2. Because of climate change, the global distribution of Aedes
mosquitoes is expanding, increasing the number of people exposed to Aedes-
borne diseases. In a recent study, researchers predict that more than 8
billion people could be at risk of dengue by 2080 3.

In the city of Rio de Janeiro, Brazil, the first Zika epidemic occurred
between 2015 and 2016, with more than 35 thousand confirmed cases 4. The city
is the second-largest in Brazil, with approximately 6.3 million inhabitants,
and its main tourist destination. Rio de Janeiro has a tropical climate and a
favourable environment for the Ae. aegypti mosquitoes, which are highly
adapted to urban settings. Despite efforts to control the vector population,
the city has suffered from dengue epidemics every three to four years, in
general 5, 6, 7. The widespread presence of the mosquito also allowed the
entry and rapid dispersion of Zika and chikungunya viruses 4. This
epidemiological scenario highlights the need for novel strategies to help
design interventions that are more effective in decreasing the burden of
established Aedes-borne diseases and preventing emerging and re-emerging
arboviruses from causing new outbreaks. In this sense, we propose a model that
has the potential to help prioritise interventions by identifying areas with
outlying risks with respect to the entire region while accounting for
covariates.

Motivating the proposed model, we have available the Zika cases counts
aggregated by neighbourhood for the period of the first Zika epidemic in Rio
de Janeiro city. The data come from the Brazilian Notifiable Diseases
Information System (SINAN – Sistema de Informação de Agravos de Notificação).
In Brazil, cases attending healthcare facilities with a suspected diagnosis of
Zika are reported to this system, usually by the physician. The standardised
morbidity ratios (SMR) for the Zika counts by neighbourhoods during the study
period are presented in Figure 1. Although the epidemic affected most of the
city, some neighbourhoods seem to have been hit harder than others. The
diversity of the territory of Rio de Janeiro is possibly an important factor
influencing this. Rio’s territory is heterogeneous in terms of demographic,
socio-economic, and environmental characteristics that are involved in the
distribution of Aedes-borne diseases 8.

For this analysis, we have available the socio-development index, an index
that includes indicators related to sanitation, education and income, and for
which higher values represent better socio-economic conditions. In places with
inadequate sanitary conditions, the female Ae. aegypti can more easily find
any type of container filled with water to deposit her eggs. In Rio de
Janeiro, a very unequal city, the socio-development index ranges from 0.282
(in Grumari, a neighbourhood in the West region) to 0.819 (in Lagoa, South
region) 9. In a previous study from our group, the socio-development index
presented a strong inverse association with the spatio-temporal distribution
of chikungunya cases in Rio de Janeiro 8.

(a) Map of the SMR (b) Histogram of the SMR

Figure 1: SMR for the Zika counts across the 160 neighbourhoods of Rio de
Janeiro, 2015-2016.

###  1.1 Literature review

In the last 30 years, the area of disease mapping has experienced an enormous
growth. This is because it is an important tool for decision makers to obtain
reliable areal estimates of disease rates over a region of interest. Disease
mapping methods further help understanding the underlying associations between
covariates and the disease risk. Commonly, Bayesian hierarchical models are
used to model the disease cases observed in different areas that form a region
of interest. The number of cases in an area is assumed to follow a Poisson
distribution whose mean is decomposed as the product of an offset by the
relative risk of the disease. Further, in the log scale, the relative risk is
decomposed as the sum of covariates and latent (unobserved) areal effects. The
latent components also accommodate overdispersion as this decomposition of the
log-relative risk can be seen as a Poisson-lognormal mixture model.

Usually, these latent effects follow a spatial structure, a priori, such that
neighbouring locations will adjust in a similar fashion after accounting for
the available covariates. Indeed, it seems natural to expect that areas that
are close to each other are more correlated than areas that are further apart.
However, some areas may present outlying risks that decision makers might want
to identify after accounting for covariates. This issue of spatial
heteroscedasticity has been increasingly considered over the recent years. For
instance, regarding geostatistical data, Palacios and Steel10 proposed a log-
normal scale mixture of a Gaussian process to accommodate heavy tails. In
disease mapping, Congdon11 also considered a scaling mixture model to adjust
for locally outlying areas. In this model, the mixture parameters are assumed
to follow a Gamma distribution a priori with mean 1 such that no
overdispersion is assumed on average a priori. The scale mixture parameters
introduced by Congdon11 allow for identification of the areas that are
different from their neighbours. More recently, in social geography, Dean et
al.12 addressed these local discrepancies termed "social frontiers" by
changing the neighbouring structure according to the observed data. The
authors study the proportion of non-UK population in the neighbourhoods of
Sheffield in England. They first test for statistical differences between the
observed proportions in neighbouring areas. When two neighbouring proportions
are found to be statistically different, the spatial structure is updated such
that the two areas are not considered as neighbours. This approach differs
from Congdon’s as it is a two-step procedure that implies changing the
neighbourhood structure. Congdon11, on the other hand, introduces scaling
mixture parameters in order to deal with spatial discrepancies. In this paper,
we are interested in identifying areas with outlying risks with respect to the
entire region, after accounting for the effect of covariates. Decision makers
may be interested in detecting these outliers in order to understand the
disease distribution and prioritise interventions.

A commonly used spatial model for the latent effects that does not accommodate
outliers is the conditional auto-regressive (CAR) prior13. Let
$\bm{b}=[b_{1},\dots,b_{n}]^{\top}$ be the vector of latent effets for the $n$
areas of the region of interest. Under the CAR prior distribution, it is
assumed that

$b_{i}\mid\bm{b}_{(-i)},\sigma_{b}^{2}\sim\mathcal{N}\left(\frac{1}{d_{i}}\sum_{j=1}^{n}w_{ij}b_{j},\frac{\sigma_{b}^{2}}{d_{i}}\right),\
i=1,\dots,n,$ (1)

where $\bm{b}_{(-i)}=[b_{1},\dots,b_{i-1},b_{i+1},\dots,b_{n}]^{\top}$,
$\bm{W}=[w_{ij}]$ is the $n\times n$ matrix of weights, $w_{ij}$, that defines
the neighbourhood structure and where $d_{i}=\sum_{j=1}^{n}w_{ij}$. Note that
$\sigma_{b}^{2}$ is the variance of the conditional distribution of $b_{i}$
given its neighbours. It can be shown14 that the joint distribution of
$\bm{b}$ is proportional to
$p(\bm{b})\propto\exp\left[-(1/2\sigma^{2}_{b})\bm{b}^{\top}\bm{Q}\bm{b}\right]$,
with $\bm{Q}=\bm{D}-\bm{W}$, where $\bm{D}=diag(d_{i})$. To ease notation, let
$\bm{b}\sim\mathcal{N}(\bm{0},\sigma_{b}^{2}\bm{Q}^{-}),$ where $\bm{Q}^{-}$
is the generalised inverse of $\bm{Q}$. The CAR distribution is not a proper
multivariate normal distribution as the "precision" matrix, $\bm{Q}$, is not
positive definite. One issue with the CAR model is that it does not perform
well when there is no underlying spatial structure in the data 15.

To accommodate the presence of independent latent effects, Besag et al.16
proposed the so-called BYM model, where each areal latent effect, $b_{i}$, is
decomposed as the sum of an unstructured component, $\theta_{i}$, and a
spatially structured component, $u_{i}$:

$b_{i}=u_{i}+\theta_{i},\ i=1,\dots,n,$ (2)

with
$\bm{\theta}=[\theta_{1},\dots,\theta_{n}]^{\top}\sim\mathcal{N}(\bm{0},\sigma^{2}_{\theta}\bm{I})$,
where $\bm{I}$ is the $n$-dimensional identity matrix, independent from
$\bm{u}=[u_{1},\dots,u_{n}]^{\top}\sim\mathcal{N}(\bm{0},\sigma^{2}_{u}\bm{Q}^{-}),$
where $\sigma^{2}_{\theta}$ is the variance of the marginal distribution of
the unstructured components and $\sigma^{2}_{u}$ is the variance of the
conditional distribution of the spatially structured components. As pointed
out by MacNab17, this model presents an identifiability issue as the two
variance components cannot be distinguished.

To avoid the introduction of two random effects for each area, Leroux et al.18
proposed an alternative distribution for the latent spatial effects that
includes a spatial dependence parameter, $\lambda$. The latter is a mixing
parameter in the unit interval that allows the variance of the latent effects
to be decomposed into a weighted sum between an unstructured and a spatially
structured variance components. The conditional distribution of each latent
effect has the following form:

$b_{i}\mid\bm{b}_{(-i)},\lambda,\tau^{2}\sim\mathcal{N}\left(\frac{\lambda}{1-\lambda+\lambda
d_{i}}\sum_{j=1}^{n}w_{ij}b_{j},\frac{\tau^{2}}{1-\lambda+\lambda
d_{i}}\right),\ i=1,\dots,n,$ (3)

where $\tau^{2}$ is again the variance of the conditional distribution of
$b_{i}$ given its neighbours. This prior is a proper distribution for all
$\lambda\in[0,1)$.

Sørbye and Rue19 argued that scaling spatially structured effects is essential
to ease interpretation and prior assignment of the variance parameter of the
latent effects, independently of the neighbourhood structure. Best et al.20
discuss this sensitivity to prior assignments regarding the BYM model. Indeed,
the parameters involved in the prior distribution of the latent effects are
often incorrectly interpreted as the properties of the marginal distribution
of $\bm{b}$. These actually depend on the neighbourhood structure and do not
have the same meaning in every data application. Hence, to overcome this
issue, Riebler et al.21 proposed the BYM2 model, that decomposes the latent
effects into a weighted sum of unstructured random noises with unit variance
and scaled structured components. The vector of latent spatial effects is
scaled according to the neighbourhood structure. The decomposition of the
$i$th latent effect is as follows:

$b_{i}=\sigma\left(\sqrt{1-\lambda}\theta_{i}+\sqrt{\lambda}u_{i}^{\star}\right),\
i=1,\dots,n,$ (4)

where $\lambda\in[0,1]$ and $\bm{\theta}\sim\mathcal{N}(\bm{0},\bm{I})$ is
independent from the scaled spatially structured components,
$\bm{u}^{\star}=[u_{1}^{\star},\dots,u_{n}^{\star}]^{\top}\sim\mathcal{N}(\bm{0},\bm{Q}_{\star}^{-}).$
The matrix $\bm{Q}_{\star}^{-}$ is the generalised inverse of
$\bm{Q}_{\star}$, which is a scaled version of the CAR "precision" matrix,
$\bm{Q}$: $\bm{Q}_{\star}=h\bm{Q}$. The scaling factor, $h$, is proportional
to the generalised variance that arises from a CAR model:
$h=\exp\left[(1/n)\sum_{i=1}^{n}\ln\left(\bm{Q}_{ii}^{-}\right)\right]$. Note
that the scaling factor only depends on the neighbourhood structure. This
scaled CAR prior corresponds to
$\bm{u}^{\star}=[u_{1}/\sqrt{h},\dots,u_{n}/\sqrt{h}]^{\top},$ for
$\bm{u}\sim\mathcal{N}\left(\bm{0},\bm{Q}^{-}\right).$ As stated by Sørbye and
Rue19, this scaling process allows each structured component to have a
variance of approximately 1. It results that
$\mathbb{V}(b_{i}\mid\sigma)\simeq\sigma^{2}\left[(1-\lambda)\times
1+\lambda\times 1\right]=\sigma^{2}.$ Hence, a marginal variance,
$\sigma^{2}$, is defined for the latent effects and all the parameters can be
interpreted for all neighbourhood structures.

On the other hand, spatial heteroscedasticity is not explicitly considered in
the previous models. However, it is reasonable to imagine that some areas may
have abnormally high or low disease risks. The goal of this paper is to
propose a method to accommodate and identify outlying areas. To capture
outliers, Yan22 proposed a double CAR model, where the relative risk is
decomposed into the sum of a CAR spatial effect and an additional noise whose
variance is modelled, in the log scale, through a BYM prior. Congdon11 argued
that this model formulation may lead to identifiability issues, as the three
sets of latent effects may not be distinguishable. To allow for local
disparities, Congdon11 proposed a modification of the Leroux prior by
including scale mixture parameters. More specifically, Congdon11 assumes

$b_{i}\mid\bm{b}_{(-i)},\bm{\kappa},\lambda,\tau^{2}\sim\mathcal{N}\left(\frac{\lambda}{1-\lambda+\lambda
d_{i}}\sum_{j=1}^{n}w_{ij}\kappa_{j}b_{j},\frac{\tau^{2}}{\kappa_{i}(1-\lambda+\lambda
d_{i})}\right),\ i=1,\dots,n,$ (5)

with $\kappa_{i}\overset{i.i.d.}{\sim}\mathrm{Gamma}(\nu/2,\nu/2),\
i=1,\dots,n$ and $\nu\sim\mathrm{Exp}(1/\mu_{\nu}),$ for some value of
$\mu_{\nu}$ fixed by the analyst. These positive parameters,
$\bm{\kappa}=[\kappa_{1},\dots,\kappa_{n}]^{\top}$, allow for discrepancies in
the neighbouring estimated risks, while the usual CAR-type priors aim to
locally smooth the risk surface. The scale mixture parameters are termed
outlier indicators as $\kappa<1$ captures local outliers, namely areas that
have different rates from their neighbours. Again, $\tau^{2}$ is the variance
of the conditional distribution of $b_{i}$ given its neighbours. This implies
that the interpretation of $\tau^{2}$ differs with every spatial structure and
its prior distribution should be defined with care for every data application.
It can be shown11 that the joint distribution of the latent effects is
$\bm{b}\mid\tau^{2},\lambda,\bm{\kappa}\sim\mathcal{N}\left(\bm{0},\tau^{2}\bm{Q}_{C}^{-}\right),$
where the "precision" matrix has diagonal elements
$\bm{Q}_{C_{ii}}=\kappa_{i}(1-\lambda+\lambda d_{i})$ and off-diagonal
elements $\bm{Q}_{C_{ij}}=-\lambda w_{ij}\kappa_{i}\kappa_{j}$. The diagonal
dominance condition23 states that a sufficient condition for a symmetric
matrix $\bm{Q}_{C}$ to be symmetric positive definite is
$Q_{C_{ii}}>\sum_{j\neq i}|Q_{C_{ij}}|,\ \forall i.$ Hence, it is sufficient
that $\lambda\in(0,1)$ and
$\lambda<\min_{i}\left\\{1/\left(1-d_{i}+\sum_{j\neq
i}w_{ij}\kappa_{j}\right)\right\\},$ for $\bm{Q}_{C}$ to be a valid precision
matrix. Note that if $\bm{\kappa}=\bm{1}_{n}$, then Congdon’s prior is the
Leroux prior, which is proper for $\lambda\in[0,1)$. This mixture differs from
the commonly used normal-gamma model as the scale mixture components appear
both in the mean and in the variance of the conditional distribution.

In this paper, we propose a modification of the BYM2 prior21 that is able to
identify areas with outlying disease risks, after accounting for the effect of
covariates. A scale mixture is introduced in the variance of the BYM2 model as
in Congdon11. The proposed model keeps the appealing property of parameter
interpretation while capturing and identifying potentially outlying areas.
This paper is organised as follows: Section 2 describes the proposed model
with the proposed decomposition of the latent effects. A simulation study
showcases the performance of the proposed model in section 3 under different
scenarios. Additionally, the application of the proposed model to the data
presented in section 1 from the 2015-2016 Zika epidemic in the 160
neighbourhoods of Rio de Janeiro is shown in section 3. We conclude by a
discussion in section 4.

##  2 Proposed model

Let a region of interest be partitioned into $n$ non-intersecting areas. Let
$Y_{i}$ be the number of cases in area $i,\ i=1,\dots,n,$ and $E_{i}$, the
expected number at risk in that area. The counts are modelled through the
following Poisson model:

$Y_{i}\mid E_{i},\mu_{i}\sim\mathcal{P}(E_{i}\mu_{i}),$

where $\mu_{i}$ denotes the relative risk in area $i$ and $E_{i}$ is an
offset. Commonly, the risk is decomposed in the log scale as follows:

$\ln(\mu_{i})=\beta_{0}+\bm{x}_{i}\bm{\beta}+b_{i},$

where $\beta_{0}$ is the mean log risk, $\bm{x}_{i}$ is a $p$-dimensional
vector with the explanatory variables in area $i$, associated with the $p$
coefficients $\bm{\beta}$, and $b_{i}$ is a random effect for area $i$. The
latter is included in order to allow for overdispersion in the Poisson model
that would otherwise assume equal mean and variance for area $i$. The latent
areal effects can also accommodate an assumed underlying spatial structure in
the data. To that end, a spatial structure is defined through the matrix
$\bm{W}=[w_{ij}]$. Throughout this paper, we assume that two areas are said to
be neighbours if they share a border. This implies that $w_{ij}=1$ if areas
$i$ and $j$ are neighbours and $w_{ij}=0,$ otherwise. In this setting,
$d_{i}=\sum_{i=1}^{n}w_{ij}$ corresponds to the number of neighbours of area
$i$. To model the latent areal effects accounting for such 0-1 spatial
structure, we propose a modification of the prior specification (4) proposed
by Riebler et al.21, that is, we assume

$b_{i}=\frac{\sigma}{\sqrt{\kappa_{i}}}\left(\sqrt{1-\lambda}\theta_{i}+\sqrt{\lambda}u_{i}^{\star}\right),\quad
i=1,\dots,n,$ (6)

where $\sigma>0$ is scaled by $\kappa_{i}>0$, and where $\lambda\in[0,1]$. The
component $\theta_{i}$ is assumed independent from $u_{i}^{\star}$. In
particular,
$\bm{\theta}\equiv[\theta_{1},\dots,\theta_{n}]^{\top}\sim\mathcal{N}(\bm{0},\bm{I})$,
and
$\bm{u}^{\star}\equiv[u^{\star}_{1},\dots,u^{\star}_{n}]^{\top}\sim\mathcal{N}(\bm{0},\bm{Q}_{\star}^{-})$.
The latter two are termed the unstructured and the scaled structured
components, respectively. Like in the BYM2 model21 (4), the "precision" matrix
is such that $\bm{Q}_{\star}=h\bm{Q}$, where the scaling factor, $h$, is
computed from the neighbourhood structure (see section 1.1). It results that
$\mathbb{V}(b_{i}\mid\sigma,\kappa_{i})\simeq(\sigma^{2}/\kappa_{i})\left[(1-\lambda)\times
1+\lambda\times 1\right]=\sigma^{2}/\kappa_{i}.$ Hence,
$\sigma^{2}/\kappa_{i}$ represents the approximate marginal variance of the
latent effect of area $i$. Moreover, the variance-covariance matrix, $\bm{V}$,
of the proposed latent effects, $\bm{b}$, is given by
$\bm{V}=\sigma^{2}\bm{K}^{-1}\left[(1-\lambda)\bm{I}+\lambda\bm{Q}_{\star}^{-}\right],$
where $\bm{K}=diag(\kappa_{i})$. Thus, the parameter $\lambda$ represents the
weight of the spatial effect in the variance of the latent process. Note that
this distribution is a proper multivariate normal for small values of
$\lambda$, depending on the neighbourhood structure. Indeed, the diagonal
dominance condition23 implies that it is sufficient that $\lambda\in[0,1]$ and
$\lambda<\min_{i}\left\\{1/\left(1-\bm{Q}_{\star_{ii}}^{-}+\sum_{j\neq
i}|\bm{Q}_{\star_{ij}}^{-}|\right)\right\\}$ for the covariance matrix,
$\bm{V}$, to be valid.

In a nutshell, the proposed model uses interpretable parameters to accommodate
outlying areas while identifying them. The proposed model points at
neighbourhoods that need heavy-tailed latent effects after accounting for some
covariates. Different from Congdon11 (5), the proposed model introduces
parameters that intervene on the marginal distribution of the latent effects.
Therefore, their interpretation remains the same for any spatial structure and
the prior assignment is simplified. This concerns the marginal variance,
$\sigma^{2}$, as well as the scaling mixture parameters,
$\bm{\kappa}=[\kappa_{1},\dots,\kappa_{n}]^{\top}$.

###  2.1 Prior specification of the scale mixture component

Different priors may be considered for the scale mixture components,
$\bm{\kappa}$. A natural choice, and used by Congdon11, is to assume:

$\kappa_{i}\overset{i.i.d.}{\sim}\mathrm{Gamma}\left(\nu/2,\nu/2\right),\
i=1,\dots,n,\quad\mbox{and}\quad\nu\sim\mathrm{Exp}(1/\mu_{\nu}),$ (7)

where the hyperparameter’s mean, $\mu_{\nu}$ controls the magnitude of $\nu$.
When $\lambda=0$, and marginalising the proposed distribution (6) of the
latent effect, $b_{i}$, with respect to $\kappa_{i}$ yields a $t_{\mu_{\nu}}$.
The introduction of $\kappa_{i}$ hence allows for heavier tails than a
Gaussian distribution for the latent effects. In this case, $\mu_{\nu}$
corresponds to choosing the prior distribution of the degrees of freedom of
the resulting $t$ distribution which impact the moments of the distribution as
well as its tails. A large $\mu_{\nu}$ results in a distribution very close to
being normal, which is inadequate to capture outliers. On the other hand,
$\mu_{\nu}<3$ implies a $t$ distribution whose variance is not defined. Some
simulation studies showed that setting $\mu_{\nu}=4$ performed well, which
agrees with Gelman et al.24.

The above prior specification (7) is the same one used by Congdon11. In
Congdon’s model, the $\kappa$’s contribute with two aspects to the conditional
distribution of the latent effects: in the mean and in the variance. To
visualise the $\kappa$’s impacts in Congdon’s model, say $\kappa_{i}<1$ for
$i\sim j$, that is for neighbouring areas $i$ and $j$. First, as $\kappa_{i}$
appears in the conditional mean structure of the latent effect for area $j$,
the weight of area $i$ is decreased and area $i$ contributes less to the mean
of the $j$th latent effect than its other neighbours. Second, by being part of
the variance structure, $\kappa_{i}<1$ gives area $i$ an increased variability
and thus allows it to differ from its neighbours. On the other hand, in the
proposed model, the $\kappa$’s only intervene on the marginal variances.
Therefore, $\kappa_{i}<1$ acts as an outlier indicator for the random effect
of area $i$ with respect to the whole region of interest by increasing its
variability and hence allowing it to differ from the overall mean structure.

Another possible prior specification for the $\kappa$’s is to borrow ideas
from Palacios and Steel10 who proposed the inclusion of a scale mixture
component in the variance of a Gaussian process. The authors suggest the usual
gamma mixing is not always appropriate as not all positive moments exist.
Additionally, they point out that the $t$ distribution that results from
marginalising over the gamma scaling mixture parameters may still overestimate
the overall variance and struggle to detect specific outlying areas. In
particular, they assume that the scale mixture component follows a log-
Gaussian process with the same spatial structure as the one defined for the
main Gaussian process. In our proposed model, we assume a scaled log-CAR prior
distribution for the $\kappa$’s. This form of discretisation of the method
proposed by Palacios and Steel10 is applied to the latent effects, which
include both the structured and unstructured components, in order to keep the
interpretative property of the parameters. This contrasts with the method
proposed by Palacios and Steel10 as they introduced a scale mixture only for
the spatially dependent components. The scale mixture components are modelled
as follows:

$\displaystyle\begin{split}\ln(\kappa_{i})&\equiv-\frac{\nu_{\kappa}}{2}+z_{i},\
i=1,\dots,n,\\\
\mbox{where}\quad\bm{z}\equiv[z_{1},\dots,z_{n}]^{\top}&\mid\nu_{\kappa}\sim\mathcal{N}\left(\bm{0},\nu_{\kappa}\bm{Q}_{\star}^{-}\right)\quad\mbox{and}\quad\nu_{\kappa}\sim\mathrm{Exp}(1/\mu_{\nu_{\kappa}}),\end{split}$
(8)

where $\bm{Q}_{\star}=h\bm{Q}$ is again the "precision" matrix that is scaled
by $h$, which is computed according to the spatial structure. The CAR
distribution is scaled in order to approximately have that
$\mathbb{V}[\ln(\kappa_{i})\mid\nu_{\kappa}]\simeq\nu_{\kappa}\times 1$. This
distribution is not a proper distribution since $\bm{Q}_{\star}$ is not of
full rank. Similarly to Palacios and Steel10, this prior implies
$\mathbb{E}(\kappa_{i}\mid\nu_{\kappa})\simeq 1$, which corresponds to a
constant marginal variance across the areal latent effects, and
$\mathbb{V}(\kappa_{i}\mid\nu_{\kappa})\simeq\exp(\nu_{\kappa})-1,\ \forall
i$. For $\nu_{\kappa}$ close to 0, $\kappa$ is close to 1 with a small
variance. A bigger $\nu_{\kappa}$ allows the $\kappa$’s to differ greatly from
1 and to be closer to 0, when necessary. As before, smaller values of
$\kappa_{i}$ increase the latent variance for area $i$ and again,
$\kappa_{i}<1$ indicates that area $i$ is a potential outlier. Palacios and
Steel10 suggest that a reasonable prior mean for $\nu_{\kappa}$ is
$\mu_{\nu_{\kappa}}=0.2$. The many simulation studies we conducted suggest
that a sensible choice for $\mu_{\nu_{\kappa}}$ is $\mu_{\nu_{\kappa}}=0.3$ as
it yields $[0.2,2.4]$ as the 95% prior credible interval for the $\kappa$’s.
This includes $\kappa_{i}=1$ while allowing for departure from $\kappa_{i}=1$,
to accommodate the potentially outlying random effect of area $i$. This prior
specification for the $\kappa$’s allows for a detection of outlying areas with
respect to the whole region of interest while also accommodating local
discrepancies thanks to the inclusion of the spatial structure, $\bm{Q}$.

###  2.2 Inference procedure

Following the specifications discussed in the previous section the resultant
posterior distribution, regardless of the prior specification for $\kappa_{i}$
does have closed analytical form. Therefore, the resultant posterior
distributions are approximated through computational methods. In particular,
Markov Chain Monte Carlo (MCMC) methods are considered. The Hamiltonian Monte
Carlo method implemented in the R package rstan25 is used for the simulation
studies and real data application that follow. Morris et al.26 note that the
No U-Turn Sampler implemented in rstan is more efficient than other MCMC
samplers to obtain reliable estimates of the posterior distributions induced
by the complex auto-regressive type of models that are of interest in this
paper.

One way to approximate a proper posterior distribution when assigning a CAR
prior, is to add a sum-to-zero constraint on the parameters in order to
distinguish them from any added constant. This is necessary due to the
invariance of the CAR distribution to the addition of a constant23. The sum-
to-zero constraint is applied to the spatial components of the proposed model,
$\bm{u}$, that need to be distinguished from the global intercept,
$\beta_{0}$, as well as for the log-CAR scale mixture components, where
$\ln(\bm{z})$ must be different from the mean $-\nu_{\kappa}/2$. More
precisely, we add a soft sum-to-zero constraint in both cases, that is
$\sum_{i=1}^{n}u_{i}\sim\mathcal{N}(0,n/1000)$ and
$\sum_{i=1}^{n}\log(z_{i})\sim\mathcal{N}(0,n/1000)$, respectively. The rstan
implementation of the BYM2 model is discussed by Morris et al.26 and the code
for the proposed model, which is a modification of the BYM2, is available in
Appendix A.

The scaling factor, $h$, needed in the BYM2 and, hence, in the proposed model
is computed through the R package R-INLA (Integrated Laplace Approximation,
Rue et al.27, www.r-inla.org) as explained by Riebler et al.21.

##  3 Data analyses

In this section, we present the results of one of the many simulation studies
that were conducted and fit the proposed model to data obtained from the first
Zika epidemic that took place in 2015 and 2016 in Rio de Janeiro. In both
cases, we consider the two parametrisations of the proposed model, which
correspond to the two prior specifications of the scaling mixture components
described in section 2.1. First, in the simulation study, we generate data for
the 160 districts of Rio de Janeiro and contaminate 10 areas. The goal is to
check if our proposed model is able to identify these 10 generated outliers.
We compare the performance of the two parametrisations of the proposed model
to the one from Congdon11. Second, in the Zika data application, we compare
the results of our proposed model to Congdon’s as well as the BYM221 and
Leroux18 models. We identify some potentially outlying districts which might
be of interest to decision makers.

###  3.1 Simulation study

To assess the performance of the proposed model, many simulation studies are
conducted. First, in Appendix C and in Appendix D, two simulation studies show
the ability of the two versions of the proposed model to recover the true
parameters when data are generated from the model itself. Then, Appendix E
presents a simulation study without any outlying areas which results in the
proposed model performing well compared to the prior by Congdon11. Here we
present the results from a simulation study wherein some areas are
contaminated into outlying areas, to compare the performance of the proposed
model in comparison to the one proposed by Congdon. The $n=160$ districts of
Rio de Janeiro and their neighbourhood structure are used as the region of
study.

First, the latent effects are simulated according to a proper CAR model,
denoted PCAR hereafter:
$\bm{b}\sim\mathcal{N}\left(\bm{0},\sigma_{b}^{2}[\bm{D}-\alpha\bm{W}]^{-1}\right),$
with $\alpha$ and $\sigma_{b}$ arbitrarily set to $0.7$ and $\sqrt{0.7},$
respectively. Note that the parameter $\alpha$ is reasonably small as
$\alpha=1$ corresponds to a CAR prior, wherein the spatial structure is
essential, while $\alpha=0$ produces independent random effects with variances
proportional to the inverse of the number of neighbours. Then, 10 far apart
districts are contaminated; namely Estácio, Vila Kosmos, Jardim Carioca, Costa
Barros, Freguesia (Jacarepaguá), Vargem Pequena, Padre Miguel, Inhoaíba,
Rocinha and Jacarezinho. The number of cases are contaminated by adding random
uniform values to the initially generated PCAR latent effects. Finally, using
the contaminated effects, $\bm{b}$, 100 populations of size $n=160$ are
created, according to a hierarchical Poisson model. That is,
$Y_{i}\sim\mathcal{P}\left(E_{i}\exp\left[\beta_{0}+\beta
x_{i}+b_{i}\right]\right),$ with $\beta_{0}=-0.1$, $\beta=-4$ and the offsets,
$[E_{1},\dots,E_{n}]^{\top}$, and covariates, $[x_{1},\dots,x_{n}]^{\top}$,
taken from the real data application to the Zika counts that is presented in
section 3.2. The only source of randomness across the 100 replicates comes
from the repeated sampling from a Poisson distribution. Indeed, the spatial
effects are generated only once and are considered as fixed. The standardised
morbidity ratios, $Y/E$, for the 50th simulation are mapped in Figure 2.

Figure 2: Standardised morbidity ratio (SMR$=Y/E$) for the 50th simulation. *:
contaminated districts

The Congdon prior is compared to the proposed model with the two scale
mixtures described in section 2.1. The first version of the proposed model is
denoted BYM2-Gamma and the second, BYM2-logCAR. For the three models, the
intercept and regression coefficient are given a quite vague prior:
$\beta_{0},\beta\overset{i.i.d.}{\sim}\mathcal{N}(0,10^{2})$ and the mixing
parameter, $\lambda$, is assigned a uniform, $\mathcal{U}(0,1)$, prior
distribution. The same fairly informative $\mathcal{N}_{+}(0,1)$ prior is
considered for $\sigma$, which is a marginal standard deviation in the
proposed model, while it is a conditional standard deviation in Congdon’s.
Finally, in the BYM2-Gamma and Congdon models, we set a priori that
$\kappa_{i}\overset{i.i.d.}{\sim}\mathrm{Gamma}(\nu/2,\nu/2)$ and
$\nu\sim\mathrm{Exp}(1/4)$. For the BYM2-logCAR parametrisation of the
proposed model, we set $\ln(\bm{\kappa})=-(\nu/2)\bm{1}+\bm{z},$ with
$\bm{z}\mid\nu\sim\mathcal{N}(\bm{0},\nu\bm{Q}_{\star}^{-}),$ where
$\bm{Q}_{\star}$ is defined through the neighbourhood structure of Rio de
Janeiro and $\nu\sim\mathrm{Exp}(1/0.3)$.

The models are fitted through the R package rstan (Stan Development Team,
2020) and the code to fit the proposed BYM2-Gamma model is available in
Appendix A. For each dataset, the MCMC procedure consists of 2 chains of
20,000 iterations with a 10,000 burn-in period and a thinning factor of 10.
Convergence of the chains is assessed through trace plots, effective sample
sizes and the $\widehat{R}$ statistic (Gelman et al.28, Vehtari et al.29).

In terms of WAIC30, for which smaller values are preferred, the proposed
models always perform better than Congdon’s, as shown in Figure 3. On average,
the BYM2-logCAR model yields a criterion of 921.0, the BYM2-Gamma model
presents a WAIC of 925.7, versus 935.2 for Congdon’s.

Figure 3: WAIC across the 100 replicates for the proposed models and
Congdon’s. Dashed lines: mean WAIC for each model

Figure 4 and Figure 5 summarise the inference for the parameters $\beta_{0}$,
$\beta$, $\lambda$, $\nu$ and $\sigma$ for the three models across the 100
replicates. The posterior means (solid circles) and posterior 95% credible
intervals (vertical lines) are shown for every parameter for each model across
the 100 replicates. All models recover the true intercept and the true
coefficient $\beta$ as shown in Figure 4. In Figure 5, for the hyperparameter
of the scaling mixture, $\nu$, the results cannot be directly compared between
the BYM2-logCAR model and the two others as $\nu$ does not play the same role
under the different parametrisations. Both the BYM2-Gamma and Congdon’s models
lead to similar results. Note that for the latter two, there seems to be
enough information in the data to learn about this higher level parameter as
the posterior summaries differ from the prior mean of 4 with prior 95%
credible interval, $[0.1,14.7]$. The same comment applies to the
hyperparameter in the BYM2-logCAR model: the credible bounds go from
$[0.0,1.1]$ a priori to $[0.6,2.5]$ on average a posteriori. Regarding the
mixing parameter, $\lambda$, the two versions of the proposed model estimate
it to be smaller than Congdon’s. In fact, $\lambda$ is not directly comparable
between Congdon’s and the proposed models. In Congdon’s model, $\lambda$
represents the weight of the spatial structure on the means and variances of
the conditional distributions of the latent effects, while it corresponds to
the overall weight of the spatial structure on the marginal distribution of
the latent effects in the proposed model. Thus, Congdon’s model points on
average to a 75% weight of the spatial structure for the conditional
distributions of the latent effects and the two parametrisations of the
proposed model indicate a fairly small spatial structure, marginally. This
result is sensible for two reasons. First, as previously stated, the parameter
$\alpha=0.7$ is fairly small in the PCAR generating distribution of the
spatial effects. Second, by introducing contamination, we reduce yet more the
overall importance of the spatial structure in the data.

Again, a distinction needs to be made between the two parametrisations of the
proposed model and Congdon’s concerning the parameter $\sigma$, which appears
in the conditional variance of the latent effects in Congdon’s model and in
the marginal one in the proposed model. The BYM2-Gamma model estimates a small
$\sigma$, which corresponds to the overall standard deviation that is common
to all the districts. The true marginal standard deviation is then captured
due to the inclusion of the $\kappa$’s that increase the areal variances only
when necessary.

Figure 4: Posterior summaries of the model coefficients for the three models
across the 100 replicates.
Solid circle: posterior mean; Vertical lines: 95% posterior credible interval;
Dashed horizontal line: true value.  Figure 5: Posterior summaries of the
model parameters for the three models across the 100 replicates.
Solid circle: posterior mean; Vertical lines: 95% posterior credible interval.

Regarding the detection of outliers, which is the main focus of the proposed
model, Figure 6 shows how often each district is detected as a potential
outlier by the three models. Overall, all models are able to find the
contaminated districts, as shown by the red areas that are indicated to have
been contaminated by the white stars. However, the BYM2-logCAR model detects
Vila Kosmos as an outlier only 83% of the time, which may be due to chance
from the contamination process. The BYM2-Gamma version of the proposed model
does not find any other outliers. The BYM2-logCAR parametrisation does
identify the district of Tauá as an additional outlier in 21% of the
replicates. This is probably due to the spatial structure that is induced in
the scaling mixture components, as Tauá shares a border with Jardim Carioca,
which is one of the 10 contaminated districts. Furthermore, Congdon’s model
often points two other districts as potential outliers, namely Senador
Vasconcelos (68%) and Penha (24%). This tendency of Congdon’s prior and, at a
lower scale, the proposed BYM2-logCAR model to report more outliers than were
generated was further explored in Appendix E. In the simulation study shown in
Appendix E, data are generated from the hierarchical Poisson model with latent
effects following a PCAR distribution and no district is contaminated. Hence,
there are no outliers. This study reports that the proposed BYM2-Gamma model
does not identify any district as an outlier whereas the proposed BYM2-logCAR
model detects one district as a potential outlier in 8% of the replicates,
showing a sensitivity to the neighbourhood structure. Furthermore, Congdon’s
model detects 6 districts as potential outliers up to 8% of the times, showing
even more sensitivity to the neighbourhood structure.

(a)  (b)

(c)

Figure 6:  Maps of the outliers indicated by $\kappa_{u}<1$, where
$\kappa_{u}$ is the upper bound of the posterior 95% credible interval of
$\kappa$. a) BYM2-Gamma model; b) BYM2-logCAR model; c) Congdon’s model.

Note that all models are able to recover the observed values. Therefore, in
terms of fitted values, there is not a big difference among the models.
However, the main advantage of the proposed model is its ability to identify
areas with potentially outlying random effects while having interpretable
parameters. This interpretability helps when specifying the prior
distributions of the hyperparameters.

###  3.2 Cases of Zika during the 2015-2016 epidemic in Rio de Janeiro

The total numbers of cases of Zika were recorded across the 160 neighbourhoods
of Rio de Janeiro during the first epidemic, which took place between 2015 and
2016. Let $Y_{i}$ be the disease count in district $i=1,\dots,160.$ A
hierarchical Poisson model is fitted to these data with offsets, $E$, computed
from, $P$, the areal population sizes,
$E_{i}=P_{i}\left(\sum_{j}Y_{j}/\sum_{j}P_{j}\right)$. We consider a socio-
development index, $x$, as an explanatory variable for the number of cases.
Identifying districts with potentially outlying risks, after accounting for
covariates, may be useful for decision makers to understand how to prevent
Zika and where to start from. The distribution of Zika is described through a
map and a histogram of the standardised morbidity ratio (SMR), $Y/E$, in
Figure 1 in section 1. Some districts seem to present different SMR values
than the mean surface, such as the island Paquetá, Barra de Guaratiba and
Pedra de Guaratiba, with SMRs of 7.3, 6.5 and 5.9, respectively. In the lower
tail of the SMR distribution, three districts did not record any cases and
thus present null SMRs, namely Gericinó, Vasco da Gama and Parque Colúmbia.
However, the SMR being an exploratory tool, one cannot conclude that high or
low SMR values necessarily indicate outlying districts. Therefore, we are
interested in comparing which districts are identified as potential outliers,
after accounting for the socio-development index, by the two versions of the
proposed model and Congdon’s. The same priors are defined for the parameters
as in the simulation study presented in section 3.1 and the two versions of
the proposed model are again denoted BYM2-Gamma and BYM2-logCAR. We further
compare the performance of the three models to the BYM2 and Leroux models
which do not accommodate potential outliers.

All models are fitted in rstan (Stan Development Team, 2020) with 2 chains of
20,000 iterations thinned by 10 and of which 10,000 are burnt. As assessed by
the trace plots, the effective sample sizes and the $\widehat{R}$
statistics28, the two chains have mixed well for all five models and
convergence is attained.

The results from the fitted models are presented in Table 1 and Figure 7. In
terms of WAIC, the proposed BYM2-Gamma model performs best among the five
considered. There is an important performance gain when accommodating outliers
(BYM2-Gamma, BYM2-logCAR and Congdon: 1335, 1342 and 1337, respectively, vs
BYM2 and Leroux: 1371 and 1375, respectively). Congdon’s prior, that focuses
on local outliers, does not seem to perform significantly worse than the
BYM2-Gamma model, which concentrates on outliers with respect to the general
distribution of the spatial effects. However, Congdon’s model yields a smaller
WAIC than the BYM2-logCAR proposed prior. Interestingly, even though the
BYM2-Gamma model has 160 more parameters than Congdon’s, its effective number
of parameters is similar (80 vs 81).

Regarding the intercept, $\beta_{0}$, the proposed models and Congdon’s give
similar results whereas the Leroux and BYM2 models yield smaller posterior
means and lower credible interval bounds. This is probably due to the
difference in the spatial effects that are allowed to be more extreme in the
Congdon, BYM2-Gamma and BYM2-logCAR models. All five models indicate that
richer districts have smaller risks of Zika (posterior means for $\beta$ of -3
for the BYM2 and Leroux models and -4 for the proposed models and Congdon’s).
We cannot directly compare the parameters $\lambda$ and $\sigma$ between the
BYM2-type models and Leroux-type priors as these appear in the marginal and
conditional distributions of the latent effects, respectively. Marginally, the
BYM2-type models yield similar weights of the spatially structured components
on the latent effects (posterior means for $\lambda$ of 0.6 and 0.7). However,
a difference is noted between the Leroux-type models. The Leroux prior
estimates $\lambda$ to be very close to 0 (posterior credible interval of
$[0.00,0.01]$). This contrasts with the posterior credible interval of
$[0.5,1.0]$ that results from fitting Congdon’s model. This difference is
probably due to the presence of outliers in the data, which results in the
Leroux model finding more random noise in the latent effects. The same
observation can be made for the marginal standard deviation, $\sigma$,
regarding the BYM2-type models. The posterior credible interval for $\sigma$
is significantly higher in the BYM2 model compared to the two parametrisations
of the proposed model. Indeed, the proposed models are able to estimate a
smaller overall variance for the latent effects which is then adjusted through
the $\kappa$’s when needed. Finally, it can be noted that there is enough
information in the data to learn about the hyperparameter $\nu$. This
parameter was assigned a prior mean of 4 and prior 95% credible interval of
$[0.1,14.7]$ for the BYM2-Gamma and Congdon models and resulted in posterior
means of about 2 and posterior 95% credible intervals of about $[1,3]$. The
BYM2-logCAR model assigned an exponential distribution with mean 0.3 for
$\nu$, inducing a prior 95% credible interval of $[0.0,1.1]$, and yielded a
posterior credible interval of $[0.7,2.3]$, showing the need for some
$\kappa$’s to be different from 1, a posteriori.

| BYM2 | BYM2-logCAR | BYM2-Gamma | Congdon | Leroux
---|---|---|---|---|---
Model fit
WAIC | 1371.2 | 1342.3 | 1335.6 | 1337.5 | 1375.1
$p_{W}$ | 88.6 | 82.3 | 80.0 | 81.0 | 89.2
Parameters’ posterior summaries
| Mean (95% CI) | Mean (95% CI) | Mean (95% CI) | Mean (95% CI) | Mean (95% CI)
$\beta_{0}$ | 1.6 (0.4,2.8) | 2.5 (1.3,3.5) | 2.5 (1.7,3.4) | 2.4 (1.4,3.2) | 1.6 (0.7,2.6)
$\beta$ | -2.8 (-4.8,-0.8) | -4.2 (-5.8,-2.3) | -4.3 (-5.6,-2.9) | -4.0 (-5.4,-2.6) | -2.9 (-4.4,-1.3)
$\lambda$ | 0.7 (0.4,0.9) | 0.6 (0.2,0.9) | 0.7 (0.3,1.0) | 0.8 (0.5,1.0) | 0.0 (0.0,10-2)
$\sigma$ | 0.8 (0.7,0.9) | 0.4 (0.3,0.5) | 0.4 (0.3,0.5) | 0.6 (0.4,0.8) | 0.7 (0.6,0.8)
$\nu$ | - | 1.4 (0.7,2.3) | 2.2 (1.4,3.3) | 1.9 (1.3,2.8) | -
Table 1: Model assessment (WAIC) and parameter posterior summaries: posterior
mean and 95% credible interval (CI) for BYM2, BYM2-logCAR, BYM2-Gamma, Congdon
and Leroux.

We now focus on the outliers detected by the proposed models and Congdon’s, as
shown in Figure 7. District $i$ is again found to be a potential outlier,
after accounting for the socio-development index, if the upper bound of the
posterior 95% credible interval of $\kappa_{i}$ is below 1. The detailed
posterior summaries (posterior mean and posterior 95% credible interval) for
the scale mixture components estimated in all three models can also be found
in Table LABEL:tab:kappa_PostSum in Appendix B. In Figure 7, the blue and red
coloured districts help distinguish the detected outliers on the lower tail of
the SMR distribution from the ones on the upper tail. After accounting for the
socio-development index, some districts are pointed out by the three models,
such as Gericinó, Parque Colúmbia, Vasco da Gama and Maré, on the lower tail
of the SMR distribution, Barra de Guaratiba, Bonsucesso and Vista Alegre, on
the upper tail. However, the BYM2-Gamma and Congdon’s models both do not point
out Paquetá in the upper tail whereas the BYM2-logCAR model detects it. Note,
however, that the BYM2-Gamma model is close to identifying Paquetá as an
outlier as it results in $\kappa_{u}=1.04$ for this district, where
$\kappa_{u}$ denotes the upped bound of the 95% posterior credible interval.
Figure 8 in Appendix B shows the posterior means for the latent effects
obtained from all five models. Note that the BYM2 model and the Leroux prior,
which do not accommodate heavy-tailed random effects, both estimate a high
posterior mean for Paquetá’s latent effect. This shows that there may be no
need to allow for a higher random effect for this district. The three models
do not identify Pedra de Guaratiba which has a high SMR, as shown in Figure 1.
The BYM2-logCAR model is the closest to detecting it, with a $\kappa_{u}$ of
1.09 for that district. Interestingly, the district of São Cristóvão is
detected as an outlier by both the BYM2-logCAR and Congdon’s models whereas it
is not by the BYM2-Gamma model, with $\kappa_{u}=1.2$. The BYM2-logCAR and
Congdon’s models detect more potential outliers than the BYM2-Gamma model as
these may be more sensitive to the neighbourhood structure for the $\kappa$’s.
In the BYM2-Gamma model, the $\kappa$’s are independently and identically
distributed, a priori, similarly to Congdon’s model, and only intervene on the
variance of the latent effects. In Congdon’s model, the $\kappa$’s also have
an impact on the conditional means of the latent effects. Regarding the
BYM2-logCAR model, the $\kappa$’s intervene only on the variance of the latent
effects while incorporating the neighbourhood structure in their prior
distribution. We believe that these distinct roles and distributions for the
mixing parameters explain the differences in the outliers identified after
accounting for the socio-development index.

(a)  (b)

(c)

Figure 7:  Maps of the outliers indicated by $\kappa_{u}<1$, where
$\kappa_{u}$ is the upper bound of the posterior 95% credible interval of
$\kappa$. The outliers on the lower tail are distinguished from the ones on
the upper tail of the SMR distribution. a) BYM2-Gamma model; b) BYM2-logCAR
model; c) Congdon’s model.

##  4 Discussion

In this paper, we propose a disease mapping model that is able to identify
areas with potentially outlying disease risks after accounting for the effects
of covariates. The proposed model is a scale mixture of the BYM2 model21. Two
different prior specifications are proposed for the scale mixture components
in order to account for potential outliers with respect to the whole region of
interest and to consider spatially structured mixture components. Our model
allows for a straightforward interpretation of the parameters, that is common
to every data application, while accommodating outliers. The parameters’
interpretation is eased by the scaling process of the latent spatially
structured components19. For one version of the proposed model, namely the
BYM2-Gamma, the focus is only on outliers with respect to the whole region of
interest, as decision makers may be interested in determining the areas where
an intervention is essential. For the other version of the proposed model,
namely the BYM2-logCAR, the potential outliers with respect to the whole
region are of interest while incorporating a spatial structure in the mixture
components. This parametrisation may help identify clusters of potentially
outlying areas.

A simulation study presents the performance of the two versions of the
proposed model in comparison with the one by Congdon11, which specifically
accommodates local outliers. That is, Congdon11 aims to identify areas that
are different from their neighbours. In this simulation study, the
neighbourhood structure of Rio de Janeiro is used and the latent effects of
some districts are contaminated to control the presence of outliers after
accounting for some covariate. The two versions of the proposed model always
perform best in terms of WAIC and are able to detect all the contaminated
districts. However, the proposed model with scale mixture inspired by Palacios
and Steel10 detects one particular contaminated district only 83% of the time.
The model by Congdon11, on the other hand, finds more outliers than are
generated, showing more sensitivity to the neighbourhood structure.

The many other simulation studies that are not presented in this paper show
that the BYM2-logCAR proposed model is fairly sensitive to the choice of prior
mean for the hyperparameter of the scaling mixture components. A smaller prior
mean tends to restrict the outliers indicators around 1 while a larger value
allows departure from 1. This departure from 1 may be towards 0, which
accommodates outliers, or may be towards infinity. Further research might help
assigning this prior mean. Finally, in all of our simulation studies, the
proposed model always performs at least as well as Congdon’s, and often
better, both in terms of WAIC and of outliers identification.

The cases of Zika that were recorded in Rio de Janeiro during the first
2015-2016 epidemic are analysed using the two parametrisations of the proposed
model as well as the model by Congdon11, the BYM221 and the Leroux prior18.
All five models find that there is a fairly strong negative association
between the socio-development index and the number of cases, meaning that
richer districts have lower disease risks. This finding is consistent with
previous studies conducted in Rio de Janeiro, one investigating the first
chikungunya epidemic in the city 8 and another also investigating Zika, but
using a different methodological approach 31. These studies, including ours,
indicate that improving sanitary conditions and reducing socio-economic
disparities are of paramount importance to fight Aedes-borne diseases.

After accounting for the effect of the socio-development index, some
neighbourhoods are detected as potential outliers by the proposed models and
Congdon’s, both in the lower and upper tails of the number of cases’
distribution across the districts. By focusing on local outliers, Congdon’s
model reports a few more potential outliers than the proposed model whose
interest lies only in the outlying areas with respect to the whole region of
interest, namely the BYM2-Gamma version. Out of the 21 neighbourhoods
identified as outliers, irrespective of the model, the proposed models
BYM2-logCAR and BYM2-Gamma identified 19 (90.5%) and 14 (66.7%), respectively.
Three outlying neighbourhoods did not report any cases (Parque Columbia, Vasco
da Gama and Gericinó) and were identified as outliers by the three models.

The detection of outliers with fewer cases than expected can potentially
indicate locations with underreporting issues. It is also possible that the
neighbourhood of residency is not accurate for some of the cases. It is not
uncommon in Rio de Janeiro for a person to report as the neighborhood of
residence a neighborhood that actually borders where they live. For example, a
person living in Vasco da Gama may report São Cristóvão as their district.
This may be a reason why some neighbourhoods were detected as outliers by the
BYM2-logCAR and the Congdon models but not by the BYM2-Gamma. If a given
district is accounting for a proportion of the cases that are in fact from the
neighbouring areas, this would increase the risk of this district while
decreasing the risk of its neighbours. Therefore, the inaccurate information
on the district of residency may artificially create outliers. However, this
issue may not bias the detection of outliers by the BYM2-Gamma model, because
it looks for outliers with respect to the whole study area, and not to their
neighbours.

The detection of outliers with more cases than expected after controlling for
the socio-development index pinpoints locations that should be further
investigated to understand why the elevated risk of Zika in these districts.
One possible explanation, discussed above, involves the inaccurate information
on the district of residency by a proportion of the cases. However, other
underlying characteristics, not included in the models, may be related to this
outlying behaviour. Also, it is important to emphasize that some socio-
environmental factors that influence the burden and distribution of Aedes-
borne diseases may be heterogeneous within the districts, our spatial unit of
analysis. For example, the same district may have areas with favelas (slums)
and areas with middle and upper class condominiums. The socio-development
index will not capture this intra-district social inequality, and a recent
study showed evidence about the presence of socio-economic inequalities in the
distribution of dengue, Zika and chikungunya in two Latin American cities 32.
Another possibility is the presence of large potential breeding sites, such as
dumps and vacant lots.

To conclude, we believe our proposed model to be useful to decision makers.
First, the parameters’ interpretation eases the use of our model regardless of
the data spatial structure. This may help decision makers to create a
systematic procedure to analyse data with our proposed model in which non
informative priors for the parameters could be defined for any spatial
structure. Then, the introduction of scaling mixture components improves the
recovering of the observed and potentially outlying disease risks, as assessed
by the model performance criterion (WAIC). Finally, these mixture components
together with high estimated risk ratios help identify all of the potential
outlying areas in which interventions may need to be prioritised.

## Acknowledgments

Schmidt and Michal are grateful for the financial support from the Natural
Sciences and Engineering Research Council (NSERC) of Canada (Discovery Grant
RGPIN-2017-04794) and the IVADO Fundamental Research Project
(PRF-2019-6839748021).

## Data Availability Statement

The Zika and the population data analysed in this study come from the
Brazilian Notifiable Diseases Information System (SINAN - Sistema de
Informação de Agravos de Notificação) and the Brazilian Institute of Geography
and Statistics (IBGE - Instituto Brasileiro de Geografia e Estatística),
respectively, and are publicly available at the Rio de Janeiro Secretariat of
Health website
(http://www.rio.rj.gov.br/dlstatic/10112/7079759/4197436/ZIKASE2015.pdf and
http://www.rio.rj.gov.br/dlstatic/10112/10617973/4260330/ZIKASE2016.pdf, for
2015 and 2016, respectively). The sociodevelopment index data come from the
Instituto Pereira Passos and can be found at www.data.rio.

## References

  * 1 Lowe R, Barcellos C, Brasil P, et al. The Zika Virus Epidemic in Brazil: From Discovery to Future Implications. International Journal of Environmental Research and Public Health 2018; 15(1): 96.  doi: 10.3390/ijerph15010096
  * 2 World Health Organization . Vector-borne diseases. https://www.who.int/news-room/fact-sheets/detail/vector-borne-diseases; 2020\.
  * 3 Colón-González FJ, Sewe MO, Tompkins AM, et al. Projecting the risk of mosquito-borne diseases in a warmer and more populated world: a multi-model, multi-scenario intercomparison modelling study. The Lancet Planetary Health 2021; 5(7): e404–e414.  doi: 10.1016/s2542-5196(21)00132-7
  * 4 Freitas LP, Cruz OG, Lowe R, Carvalho MS. Space–time dynamics of a triple epidemic: dengue, chikungunya and Zika clusters in the city of Rio de Janeiro. Proceedings of the Royal Society B: Biological Sciences 2019; 286(1912): 20191867.  doi: 10.1098/rspb.2019.1867
  * 5 Nogueira RMR, Miagostovich MP, Schatzmayr HG, et al. Dengue in the State of Rio de Janeiro, Brazil, 1986-1998. Memórias do Instituto Oswaldo Cruz 1999; 94(3): 297–304.  doi: 10.1590/s0074-02761999000300004
  * 6 Honório NA, Nogueira RMR, Codeço CT, et al. Spatial Evaluation and Modeling of Dengue Seroprevalence and Vector Density in Rio de Janeiro, Brazil. PLoS Neglected Tropical Diseases 2009; 3(11): e545.  doi: 10.1371/journal.pntd.0000545
  * 7 Santos dJPC, Honório NA, Nobre AA. Definition of persistent areas with increased dengue risk by detecting clusters in populations with differing mobility and immunity in Rio de Janeiro, Brazil. Cadernos de Saúde Pública 2019; 35(12).  doi: 10.1590/0102-311x00248118
  * 8 Freitas LP, Schmidt AM, Cossich W, Cruz OG, Carvalho MS. Spatio-temporal modelling of the first Chikungunya epidemic in an intra-urban setting: The role of socioeconomic status, environment and temperature. PLOS Neglected Tropical Diseases 2021; 15(6): e0009537.  doi: 10.1371/journal.pntd.0009537
  * 9 Prefeitura do Rio de Janeiro . Índice de Desenvolvimento Social (IDS) por Áreas de Planejamento (AP), Regiões de Planejamento (RP), Regiões Administrativas (RA), Bairros e Favelas do Município do Rio de Janeiro - 2010\. http://www.data.rio/datasets/fa85ddc76a524380ad7fc60e3006ee97; 2018\.
  * 10 Palacios MB, Steel MFJ. Non-Gaussian Bayesian geostatistical modeling. Journal of the American Statistical Association 2006; 101(474): 604–618.
  * 11 Congdon P. Representing spatial dependence and spatial discontinuity in ecological epidemiology: a scale mixture approach. Stochastic Environmental Research and Risk Assessment 2017; 31(2): 291–304.
  * 12 Dean N, Dong G, Piekut A, Pryce G. Frontiers in residential segregation: Understanding neighbourhood boundaries and their impacts. Tijdschrift voor economische en sociale geografie 2019; 110(3): 271–288.
  * 13 Besag J. Spatial interaction and the statistical analysis of lattice systems. Journal of the Royal Statistical Society: Series B (Methodological) 1974; 36(2): 192–225.
  * 14 Banerjee S, Carlin BP, Gelfand AE. Hierarchical modeling and analysis for spatial data.  CRC press .  2014\.
  * 15 Breslow N, Leroux B, Platt R. Approximate hierarchical modelling of discrete data in epidemiology. Statistical Methods in Medical Research 1998; 7(1): 49–62.
  * 16 Besag J, York J, Mollié A. Bayesian image restoration, with two applications in spatial statistics. Annals of the institute of statistical mathematics 1991; 43(1): 1–20.
  * 17 MacNab YC. On Gaussian Markov random fields and Bayesian disease mapping. Statistical Methods in Medical Research 2011; 20(1): 49–68.
  * 18 Leroux BG, Lei X, Breslow N. Estimation of disease rates in small areas: a new mixed model for spatial dependence. In: Springer. 1999 (pp. 179–191).
  * 19 Sørbye SH, Rue H. Scaling intrinsic Gaussian Markov random field priors in spatial modelling. Spatial Statistics 2014; 8: 39–51.
  * 20 Best NG, Arnold RA, Thomas A, Waller LA, Conlon EM. Bayesian models for spatially correlated disease and exposure data. 1999; 6: 131–156.
  * 21 Riebler A, Sørbye SH, Simpson D, Rue H. An intuitive Bayesian spatial model for disease mapping that accounts for scaling. Statistical methods in medical research 2016; 25(4): 1145–1165.
  * 22 Yan J. Spatial stochastic volatility for lattice data. Journal of agricultural, biological, and environmental statistics 2007; 12(1): 25–40.
  * 23 Rue H, Held L. Gaussian Markov random fields: theory and applications.  CRC press .  2005\.
  * 24 Gelman A, Carlin JB, Stern HS, Dunson DB, Vehtari A, Rubin DB. Bayesian data analysis.  CRC press .  2004\.
  * 25 Stan Development Team . RStan: the R interface to Stan. 2020.  R package version 2.21.2.
  * 26 Morris M, Wheeler-Martin K, Simpson D, Mooney SJ, Gelman A, DiMaggio C. Bayesian hierarchical spatial models: Implementing the Besag York Mollié model in stan. Spatial and spatio-temporal epidemiology 2019; 31: 100301.
  * 27 Rue H, Martino S, Chopin N. Approximate Bayesian inference for latent Gaussian models by using integrated nested Laplace approximations. Journal of the royal statistical society: Series b (statistical methodology) 2009; 71(2): 319–392.
  * 28 Gelman A, Rubin DB, others . Inference from iterative simulation using multiple sequences. Statistical science 1992; 7(4): 457–472.
  * 29 Vehtari A, Gelman A, Simpson D, Carpenter B, Bürkner PC. Rank-Normalization, Folding, and Localization: An Improved $\widehat{R}$ for Assessing Convergence of MCMC. Bayesian Analysis 2021: 1–38. .  doi: 10.1214/20-BA1221
  * 30 Watanabe S, Opper M. Asymptotic equivalence of Bayes cross validation and widely applicable information criterion in singular learning theory.. Journal of machine learning research 2010; 11(12).
  * 31 Raymundo CE, Andrade Medronho dR. Association between socio-environmental factors, coverage by family health teams, and rainfall in the spatial distribution of Zika virus infection in the city of Rio de Janeiro, Brazil, in 2015 and 2016. BMC Public Health 2021; 21(1).  doi: 10.1186/s12889-021-11249-y
  * 32 Carabali M, Harper S, Neto ASL, et al. Spatiotemporal distribution and socioeconomic disparities of dengue, chikungunya and Zika in two Latin American cities from 2007 to 2017. Tropical Medicine & International Health 2020; 26(3): 301–315.  doi: 10.1111/tmi.13530

##  Appendix A Stan code for the proposed model

The stan code used to fit the proposed BYM2-Gamma model in the simulation
studies (section 3.1, Appendix C, Appendix E) and in the analysis of the Zika
epidemic in Rio de Janeiro (section 3.2) is presented below.

Listing 1: Stan code for the BYM2-Gamma proposed model

⬇

1data {

2 int<lower=1> N;

3 int<lower=1> N_edges;

4 int<lower=1> p; // General case where there are p covariates

5 matrix[N,p] X;

6 int<lower=1, upper=N> node1[N_edges]; // vectors of neighbourhood

7 int<lower=1, upper=N> node2[N_edges]; // structure

8 int<lower=0> y[N]; // Zika counts

9 vector<lower=0>[N] log_E; // offset

10 real<lower=0> scaling_factor; // to scale the variance of the latent
effects

11}

12

13parameters {

14 real beta0; // intercept

15 vector[p] beta;

16 real<lower=0> sigma; // marginal standard deviation

17 real<lower=0, upper=1> lambda; // mixing parameter

18 vector[N] theta; // unstructured components

19 vector[N] s; // spatially structured components

20 vector<lower=0>[N] kappa; // outlier indicator

21 real<lower=0> nu;

22}

23

24transformed parameters {

25 vector[N] convolved_re; // complete latent effect

26 for(i in 1:N){ convolved_re[i] = sqrt(1 - lambda) * theta[i] +
sqrt(lambda/scaling_factor) * s[i]; }

27}

28

29model {

30 for(i in 1:N)

31 y[i] ~ poisson_log(log_E[i] + beta0 + X[i,]*beta + convolved_re[i] *
(sigma/sqrt(kappa[i])) );

32

33 target += -0.5 * dot_self(s[node1] - s[node2]); // Prior for the spatially
structured components

34 sum(s) ~ normal(0, 0.001 * N); // Soft sum-to-zero constraint to be able to
have an intercept

35

36 for(j in 1:p){ beta[j] ~ normal(0.0, 10.0); }

37

38 beta0 ~ normal(0.0, 10.0);

39 theta ~ normal(0.0, 1.0);

40 sigma ~ normal(0.0, 1.0);

41 lambda ~ uniform(0.0, 1.0);

42 kappa ~ gamma(nu/2.0,nu/2.0);

43 nu ~ exponential(1.0/4.0);

44}

45

46generated quantities {

47 vector[N] mu_log;

48 vector[N] lik;

49 for(i in 1:N){

50 mu_log[i]=log_E[i] + beta0 + X[i,]*beta +
sigma*convolved_re[i]/sqrt(kappa[i]);

51 lik[i] = exp(poisson_log_lpmf(y[i] | mu_log[i])); // likelihood to compute the WAIC

52 }

53}

##  Appendix B Additional results from the Zika data analysis

This section shows the maps of the posterior means for the latent effects
obtained from fitting the BYM2, the proposed BYM2-Gamma, the proposed
BYM2-logCAR, Congdon’s and the Leroux models in Figure 8. Additionally, in
Table LABEL:tab:kappa_PostSum, this section provides the posterior means and
posterior 95% credible intervals for the scale mixture components,
$\bm{\kappa}$, included in the proposed BYM2-Gamma and BYM2-logCAR models as
well as in Congdon’s prior.

(a)  (b)

(c)  (d)

(e)

Figure 8:  Maps of the latent effects’ posterior means. a) BYM2 model; b) BYM2-Gamma model; c) BYM2-logCAR model; d) Congdon’s prior; e) Leroux prior. Table 2: Posterior summaries (mean and 95% credible interval) of the scale mixture components, the $\kappa$’s, in each district for the proposed BYM2-Gamma and BYM2-logCAR models and Congdon’s District |  BYM2-Gamma |  BYM2-logCAR |  Congdon
---|---|---|---
Saude |  1.4 (0.11,4.37) |  4.9 (0.12,30) |  1.1 (0.07,3.82)
Gamboa |  1.5 (0.15,4.58) |  2.4 (0.12,11.37) |  1.2 (0.09,3.84)
Santo Cristo |  0.5 (0.05,1.34) |  0.4 (0.08,1.22) |  0.5 (0.03,1.48)
Caju |  0.9 (0.06,2.59) |  0.2 (0.03,0.65) |  0.7 (0.04,2.69)
Centro |  1.3 (0.14,3.48) |  1.0 (0.17,3.16) |  1.5 (0.12,3.77)
Catumbi |  0.8 (0.08,2.14) |  0.7 (0.09,2.19) |  1.1 (0.09,2.99)
Rio Comprido |  0.9 (0.09,2.63) |  0.8 (0.11,2.33) |  1.8 (0.16,4.22)
Cidade Nova |  1.5 (0.19,4.61) |  1.0 (0.13,3.45) |  1.6 (0.10,4.69)
Estacio |  0.9 (0.09,2.51) |  0.8 (0.09,2.57) |  1.1 (0.08,3.06)
Sao Cristovao |  0.4 (0.03,1.19) |  0.2 (0.04,0.69) |  0.3 (0.01,0.97)
Mangueira |  0.8 (0.05,2.44) |  0.4 (0.04,1.47) |  0.8 (0.04,3.08)
Benfica |  1.1 (0.08,3.50) |  0.2 (0.04,0.72) |  1.0 (0.05,3.84)
Paqueta |  0.3 (0.03,1.04) |  0.3 (0.01,0.94) |  0.4 (0.02,1.16)
Santa Teresa |  0.6 (0.06,1.77) |  0.8 (0.20,1.94) |  0.6 (0.03,1.64)
Flamengo |  1.1 (0.07,3.35) |  1.7 (0.14,6.80) |  0.9 (0.05,3.16)
Gloria |  1.3 (0.10,3.84) |  1.7 (0.15,6.39) |  1.2 (0.08,3.67)
Laranjeiras |  1.2 (0.09,3.73) |  1.6 (0.16,6.11) |  1.1 (0.07,3.68)
Catete |  1.4 (0.12,4.29) |  2.2 (0.17,8.46) |  1.3 (0.09,4.20)
Cosme Velho |  1.1 (0.06,3.66) |  1.7 (0.10,7.78) |  1.0 (0.05,3.45)
Botafogo |  0.9 (0.07,2.81) |  1.2 (0.25,3.33) |  0.6 (0.03,2.05)
Humaita |  1.2 (0.10,3.68) |  1.3 (0.21,4.10) |  1.2 (0.07,3.82)
Urca |  1.1 (0.09,3.28) |  1.8 (0.13,8.16) |  1.1 (0.08,3.36)
Leme |  1.2 (0.11,3.66) |  2.1 (0.18,8.96) |  1.4 (0.12,4.28)
Copacabana |  1.5 (0.12,4.40) |  2.1 (0.25,7.69) |  1.5 (0.11,4.48)
Ipanema |  1.4 (0.11,4.66) |  2.4 (0.13,11.67) |  1.5 (0.08,4.89)
Leblon |  0.5 (0.03,1.63) |  0.7 (0.07,2.25) |  0.4 (0.03,1.60)
Lagoa |  1.2 (0.08,3.89) |  1.4 (0.19,4.80) |  1.2 (0.06,4.07)
Jardim Botanico |  1.4 (0.08,4.28) |  1.7 (0.15,7.09) |  1.3 (0.08,4.39)
Gavea |  1.0 (0.07,3.25) |  0.9 (0.12,2.94) |  0.9 (0.05,3.36)
Vidigal |  0.9 (0.07,3.02) |  0.8 (0.07,2.81) |  0.7 (0.03,2.39)
Sao Conrado |  1.1 (0.07,3.79) |  0.7 (0.09,2.14) |  1.0 (0.04,3.65)
Praca da Bandeira |  1.2 (0.11,3.40) |  0.7 (0.11,2.02) |  1.4 (0.11,4.14)
Tijuca |  1.0 (0.09,3.02) |  0.8 (0.15,2.34) |  0.8 (0.06,2.71)
Alto da Boa Vista |  1.2 (0.10,3.80) |  0.8 (0.16,2.10) |  1.0 (0.04,3.25)
Maracana |  1.3 (0.10,4.18) |  0.7 (0.08,2.37) |  1.2 (0.06,4.13)
Vila Isabel |  1.3 (0.10,3.99) |  0.7 (0.15,2.04) |  0.8 (0.04,2.68)
Andarai |  1.0 (0.07,3.29) |  1.1 (0.08,4.11) |  1.0 (0.05,3.71)
Grajau |  1.5 (0.12,4.54) |  0.9 (0.14,2.87) |  1.3 (0.08,4.19)
Manguinhos |  0.3 (0.01,0.88) |  0.2 (0.04,0.60) |  0.2 (0.01,0.82)
Bonsucesso |  0.3 (0.01,0.90) |  0.2 (0.03,0.65) |  0.2 (0.01,0.58)
Ramos |  1.1 (0.10,3.39) |  0.5 (0.04,1.82) |  0.8 (0.04,2.95)
Olaria |  0.8 (0.06,2.50) |  0.5 (0.07,1.48) |  0.7 (0.04,2.37)
Penha |  1.4 (0.12,4.38) |  0.7 (0.08,2.61) |  1.2 (0.07,4.25)
Penha Circular |  0.6 (0.06,1.88) |  0.5 (0.09,1.41) |  0.5 (0.03,1.77)
Bras de Pina |  1.0 (0.08,3.10) |  0.6 (0.08,1.99) |  0.8 (0.04,2.58)
Cordovil |  1.5 (0.12,4.36) |  0.7 (0.08,2.60) |  1.4 (0.08,4.70)
Parada de Lucas |  1.0 (0.08,3.23) |  0.6 (0.06,1.96) |  1.0 (0.06,3.32)
Vigario Geral |  0.6 (0.04,1.70) |  0.4 (0.04,1.42) |  0.6 (0.04,2.13)
Jardim America |  1.5 (0.10,4.77) |  0.5 (0.05,2.26) |  1.5 (0.09,4.94)
Higienopolis |  1.2 (0.08,3.60) |  0.6 (0.09,2.12) |  1.2 (0.07,3.99)
Jacare |  0.9 (0.07,2.83) |  0.6 (0.10,1.74) |  1.4 (0.10,3.85)
Maria da Graca |  1.2 (0.09,3.95) |  0.6 (0.07,2.05) |  1.5 (0.11,4.71)
Del Castilho |  1.4 (0.11,4.55) |  1.2 (0.11,5.04) |  1.2 (0.07,4.20)
Inhauma |  1.5 (0.12,4.59) |  1.1 (0.15,3.56) |  1.4 (0.08,4.21)
Engenho da Rainha |  1.4 (0.12,4.44) |  0.8 (0.14,2.43) |  1.1 (0.06,3.70)
Tomas Coelho |  0.9 (0.07,2.91) |  0.9 (0.14,2.64) |  1.1 (0.07,3.45)
Sao Francisco Xavier |  0.6 (0.05,1.78) |  0.4 (0.05,1.28) |  0.3 (0.01,1.12)
Rocha |  1.0 (0.09,3.04) |  0.5 (0.06,1.52) |  0.9 (0.06,3.10)
Riachuelo |  1.3 (0.11,3.74) |  0.9 (0.09,3.27) |  1.6 (0.11,4.62)
Sampaio |  0.8 (0.07,2.40) |  0.7 (0.09,2.48) |  1.2 (0.11,3.08)
Engenho Novo |  0.9 (0.08,2.80) |  0.7 (0.14,2.03) |  1.1 (0.07,2.92)
Lins de Vasconcelos |  1.5 (0.12,4.61) |  1.1 (0.14,3.81) |  1.5 (0.11,5.06)
Meier |  0.7 (0.05,2.22) |  0.7 (0.11,2.05) |  0.5 (0.03,1.71)
Todos os Santos |  1.0 (0.07,3.23) |  1.0 (0.09,3.49) |  0.8 (0.04,2.80)
Cachambi |  1.4 (0.09,4.39) |  0.8 (0.15,2.17) |  1.2 (0.07,3.80)
Engenho de Dentro |  1.4 (0.11,4.29) |  1.0 (0.23,3.06) |  1.3 (0.07,3.92)
Agua Santa |  1.3 (0.09,4.13) |  1.1 (0.17,3.84) |  1.2 (0.06,3.88)
Encantado |  1.4 (0.09,4.64) |  2.2 (0.19,9.33) |  1.3 (0.07,4.24)
Piedade |  1.4 (0.12,4.20) |  1.3 (0.23,3.79) |  1.3 (0.08,3.98)
Abolicao |  1.2 (0.08,3.77) |  1.7 (0.18,6.52) |  1.3 (0.08,3.95)
Pilares |  1.2 (0.12,3.78) |  1.2 (0.18,4.06) |  1.3 (0.10,3.78)
Vila Kosmos |  0.5 (0.04,1.53) |  0.5 (0.08,1.53) |  0.6 (0.03,2.22)
Vicente de Carvalho |  0.8 (0.05,2.53) |  0.7 (0.11,2.08) |  1.3 (0.08,4.16)
Vila da Penha |  1.2 (0.08,3.99) |  0.7 (0.09,2.56) |  0.8 (0.03,2.93)
Vista Alegre |  0.2 (0.01,0.69) |  0.3 (0.04,0.76) |  0.2 (0.01,0.55)
Iraja |  0.9 (0.05,2.83) |  0.5 (0.12,1.13) |  0.4 (0.02,1.38)
Colegio |  1.4 (0.13,4.31) |  1.3 (0.08,6.12) |  1.7 (0.10,5.31)
Campinho |  1.4 (0.10,4.38) |  1.3 (0.17,4.72) |  1.2 (0.08,3.98)
Quintino Bocaiuva |  0.8 (0.06,2.65) |  0.9 (0.16,2.49) |  0.8 (0.04,2.31)
Cavalcanti |  1.1 (0.10,3.44) |  1.0 (0.18,2.98) |  1.5 (0.11,4.43)
Engenheiro Leal |  0.7 (0.05,2.46) |  0.8 (0.08,3.12) |  1.0 (0.05,3.00)
Cascadura |  1.5 (0.11,4.49) |  1.2 (0.17,4.20) |  1.3 (0.08,4.31)
Madureira |  1.2 (0.11,3.67) |  1.0 (0.17,3.16) |  1.6 (0.12,4.53)
Vaz Lobo |  0.9 (0.07,2.68) |  0.8 (0.11,2.55) |  0.9 (0.06,3.25)
Turiacu |  1.5 (0.13,4.53) |  1.1 (0.13,4.18) |  1.3 (0.08,4.17)
Rocha Miranda |  0.9 (0.06,2.85) |  0.8 (0.16,2.17) |  0.6 (0.04,1.91)
Honorio Gurgel |  1.4 (0.11,4.06) |  1.1 (0.16,3.70) |  1.4 (0.09,4.26)
Osvaldo Cruz |  0.7 (0.05,2.09) |  0.8 (0.12,2.28) |  0.6 (0.03,1.82)
Bento Ribeiro |  1.4 (0.12,4.18) |  1.3 (0.19,4.25) |  1.3 (0.10,4.13)
Marechal Hermes |  1.4 (0.10,3.97) |  1.1 (0.17,3.59) |  1.4 (0.10,4.19)
Ribeira |  1.5 (0.11,4.47) |  3.3 (0.02,22.45) |  1.6 (0.09,5.80)
Zumbi |  0.2 (0.01,0.91) |  0.3 (0.04,1.04) |  0.2 (0.01,0.59)
Cacuia |  1.2 (0.08,3.71) |  1.1 (0.12,3.78) |  1.2 (0.07,3.81)
Pitangueiras |  0.9 (0.07,2.82) |  0.7 (0.05,2.74) |  1.2 (0.06,3.57)
Praia da Bandeira |  0.9 (0.06,2.92) |  0.9 (0.06,3.85) |  1.0 (0.05,3.14)
Cocota |  1.1 (0.08,3.74) |  1.4 (0.12,5.38) |  0.6 (0.03,2.61)
Bancarios |  1.4 (0.12,4.22) |  2.3 (0.10,10.38) |  1.5 (0.09,4.32)
Freguesia (Ilha) |  1.0 (0.08,3.18) |  1.6 (0.04,6.78) |  1.0 (0.06,2.88)
Jardim Guanabara |  1.4 (0.10,4.20) |  1.9 (0.16,8.57) |  1.3 (0.06,4.12)
Jardim Carioca |  1.4 (0.10,4.27) |  1.8 (0.20,7.86) |  1.4 (0.09,4.40)
Taua |  0.7 (0.04,2.16) |  1.0 (0.11,3.42) |  0.5 (0.03,1.60)
Monero |  1.4 (0.11,4.27) |  2.8 (0.15,11.91) |  1.4 (0.07,4.30)
Portuguesa |  1.1 (0.09,3.59) |  1.8 (0.16,7.16) |  1.0 (0.05,3.23)
Galeao |  1.4 (0.10,4.37) |  0.7 (0.10,2.50) |  1.2 (0.07,3.90)
Cidade Universitaria |  1.2 (0.07,3.82) |  0.5 (0.03,2.14) |  1.2 (0.06,4.55)
Guadalupe |  1.4 (0.13,4.18) |  0.8 (0.12,2.67) |  1.4 (0.09,4.29)
Anchieta |  1.2 (0.10,3.70) |  0.6 (0.07,2.01) |  1.7 (0.10,4.85)
Parque Anchieta |  0.3 (0.03,1.08) |  0.3 (0.04,0.82) |  0.4 (0.02,1.36)
Ricardo de Albuquerque |  0.9 (0.08,2.70) |  0.6 (0.07,1.87) |  1.1 (0.08,3.40)
Coelho Neto |  1.4 (0.12,4.78) |  0.6 (0.10,2.01) |  1.4 (0.10,4.65)
Acari |  1.4 (0.11,3.97) |  0.6 (0.05,2.39) |  1.4 (0.07,4.53)
Barros Filho |  0.8 (0.06,2.50) |  0.6 (0.09,1.88) |  0.7 (0.04,2.33)
Costa Barros |  0.3 (0.02,1.02) |  0.3 (0.04,0.92) |  0.2 (0.01,0.58)
Pavuna |  1.3 (0.12,4.02) |  0.5 (0.08,1.47) |  1.3 (0.09,4.26)
Jacarepagua |  0.2 (0.01,0.69) |  0.4 (0.10,0.79) |  0.1 (0.00,0.20)
Anil |  0.5 (0.05,1.66) |  0.4 (0.04,1.33) |  0.5 (0.02,2.16)
Gardenia Azul |  0.6 (0.04,1.85) |  0.4 (0.06,1.35) |  0.5 (0.02,2.15)
Cidade de Deus |  1.0 (0.06,3.20) |  0.7 (0.09,2.42) |  0.7 (0.03,2.84)
Curicica |  0.5 (0.03,1.63) |  0.4 (0.03,1.34) |  1.0 (0.05,3.87)
Freguesia (Jacarepagua) |  1.3 (0.08,4.03) |  0.8 (0.15,2.38) |  0.7 (0.04,2.77)
Pechincha |  1.2 (0.09,3.48) |  0.9 (0.11,3.47) |  1.1 (0.06,3.61)
Taquara |  0.6 (0.04,1.81) |  0.5 (0.09,1.39) |  1.0 (0.05,3.22)
Tanque |  1.4 (0.12,4.27) |  1.0 (0.16,3.03) |  2.1 (0.15,6.09)
Praca Seca |  0.7 (0.06,2.30) |  0.8 (0.14,2.46) |  0.8 (0.04,2.50)
Vila Valqueire |  1.4 (0.11,4.52) |  1.1 (0.16,3.65) |  1.4 (0.08,4.55)
Joa |  1.0 (0.04,3.65) |  1.3 (0.04,8.09) |  1.0 (0.04,3.56)
Itanhanga |  0.2 (0.01,0.66) |  0.3 (0.04,0.72) |  0.2 (0.01,0.72)
Barra da Tijuca |  1.4 (0.11,4.34) |  0.7 (0.09,2.21) |  1.4 (0.09,5.16)
Camorim |  0.5 (0.03,1.59) |  0.4 (0.06,1.30) |  0.3 (0.02,1.15)
Vargem Pequena |  1.4 (0.10,4.28) |  1.1 (0.06,5.08) |  1.3 (0.06,4.51)
Vargem Grande |  1.2 (0.10,3.89) |  0.5 (0.10,1.66) |  1.2 (0.08,3.87)
Recreio dos Bandeirantes |  1.4 (0.11,4.36) |  0.7 (0.07,2.49) |  1.4 (0.08,4.64)
Grumari |  0.9 (0.03,3.58) |  0.7 (0.02,4.39) |  0.9 (0.03,3.53)
Deodoro |  1.4 (0.11,4.38) |  1.0 (0.10,3.59) |  1.4 (0.10,4.29)
Vila Militar |  0.2 (0.01,0.58) |  0.2 (0.05,0.67) |  0.1 (0.00,0.37)
Campo dos Afonsos |  1.2 (0.09,3.62) |  0.9 (0.11,3.60) |  1.1 (0.05,3.72)
Jardim Sulacap |  0.8 (0.06,2.58) |  0.6 (0.10,1.72) |  0.9 (0.05,2.85)
Magalhaes Bastos |  1.3 (0.09,3.93) |  0.9 (0.04,4.45) |  1.2 (0.07,3.97)
Realengo |  1.5 (0.12,4.63) |  0.5 (0.08,1.34) |  1.3 (0.06,4.66)
Padre Miguel |  0.9 (0.08,2.70) |  0.5 (0.03,2.15) |  1.7 (0.07,5.62)
Bangu |  0.5 (0.04,1.41) |  0.3 (0.04,0.85) |  0.8 (0.05,2.62)
Senador Camara |  0.3 (0.02,0.91) |  0.3 (0.06,0.89) |  0.2 (0.01,0.74)
Santissimo |  1.2 (0.09,3.86) |  0.7 (0.07,2.60) |  1.5 (0.09,4.78)
Campo Grande |  1.4 (0.12,4.19) |  0.7 (0.13,1.90) |  1.4 (0.08,4.40)
Senador Vasconcelos |  1.2 (0.10,3.62) |  0.8 (0.06,3.66) |  0.8 (0.05,2.89)
Inhoaiba |  1.4 (0.11,4.48) |  1.8 (0.11,9.18) |  1.3 (0.08,4.54)
Cosmos |  0.9 (0.05,3) |  0.9 (0.11,2.88) |  0.6 (0.03,1.90)
Paciencia |  1.2 (0.10,3.55) |  1.1 (0.09,5.07) |  0.8 (0.05,2.49)
Santa Cruz |  1.5 (0.12,4.59) |  1.4 (0.10,6.01) |  1.2 (0.07,4.17)
Sepetiba |  0.7 (0.06,2.06) |  0.6 (0.05,2.03) |  0.7 (0.04,2.35)
Guaratiba |  1.2 (0.08,3.58) |  0.6 (0.11,1.97) |  1.7 (0.13,4.56)
Barra de Guaratiba |  0.2 (0.01,0.64) |  0.2 (0.03,0.65) |  0.2 (0.01,0.57)
Pedra de Guaratiba |  0.4 (0.03,1.26) |  0.3 (0.02,1.09) |  0.6 (0.03,2.06)
Rocinha |  0.9 (0.07,2.76) |  0.7 (0.06,2.78) |  0.7 (0.03,2.53)
Jacarezinho |  0.1 (0.01,0.35) |  0.2 (0.03,0.46) |  0.1 (0.00,0.26)
Complexo do Alemao |  1.1 (0.08,3.38) |  0.6 (0.09,1.97) |  1.0 (0.05,3.20)
Mare |  0.1 (0.00,0.16) |  0.1 (0.03,0.30) |  0.0 (0.00,0.09)
Parque Columbia |  0.0 (0.00,0.08) |  0.1 (0.01,0.23) |  0.0 (0.00,0.05)
Vasco da Gama |  0.0 (0.00,0.06) |  0.0 (0.00,0.11) |  0.0 (0.00,0.05)
Gericino |  0.0 (0.00,0.21) |  0.0 (0.00,0.15) |  0.1 (0.00,0.78)

##  Appendix C Simulation study: generating data from the proposed BYM2-Gamma
model

To assess the proposed BYM2-Gamma model’s ability to recover the truth, a
simulation study is conducted wherein data are generated from the proposed
BYM2-Gamma model. Again, the $n=160$ districts of Rio de Janeiro and their
neighbourhood structure are used. The latent effects’ unstructured and scaled
spatially structured components are generated once:

$\bm{\theta}\sim\mathcal{N}(\bm{0},\bm{I}),\quad\mbox{and}\quad\bm{u}^{\star}\sim\mathcal{N}(\bm{0},\bm{Q}_{\star}^{-}),$

where $\bm{Q}_{\star}=h(\bm{D}-\bm{W})$, with $h$, the scaling factor,
entirely defined by the spatial structure of Rio de Janeiro. An algorithm to
generate from the CAR prior is presented in Chapter 2 of Rue and Held23. The
outlier indicators, $\bm{\kappa}$, are independently generated once from a
$\mathrm{Gamma}(\nu/2,\nu/2)$, with an arbitrary $\nu=4$ to allow for fairly
heavy tails. The latent effects are then computed as

$b_{i}=\left[\sqrt{1-\lambda}\theta_{i}+\sqrt{\lambda}u_{i}^{\star}\right]\times\sigma/\sqrt{\kappa_{i}},\
i=1,\dots,n,$

where $\lambda=0.8$ and $\sigma=0.3$. Finally, 100 replicates of populations
of size $n=160$ are generated from the Poisson model

$Y_{i}\sim\mathcal{P}\left(E_{i}\exp\left[\beta_{0}+b_{i}\right]\right),$

with $\beta_{0}=-0.1$ and the offsets, $[E_{1},\dots,E_{n}]^{\top}$, taken
from the analysis of the Zika counts. The proposed BYM2-Gamma model and
Congdon’s are both fitted on the 100 replicated datasets using the same
inference procedure as in section 3.1.

Figure 9 shows that the WAIC is able to always choose the model that generated
the data, namely the BYM2-Gamma model. This simulation study is specifically
conducted to assess whether the proposed BYM2-Gamma model is able to estimate
the true values of the generated parameters. Figure 10 presents the posterior
summaries obtained from the BYM2-Gamma model across the 100 replicates for the
intercept, $\beta_{0}$, the mixing parameter, $\lambda$, the hyperparameter,
$\nu$, and the overall standard deviation, $\sigma$. For $\beta_{0}$ and
$\nu$, the BYM2-Gamma model is able to recover the true values. Regarding
$\lambda$ and $\sigma$, the BYM2-Gamma model tends to underestimate these
parameters on average but still recovers the truth in the posterior 95%
credible intervals. The interest lies particularly on the main parameters of
the model, such as the outlier indicators, $\bm{\kappa}$. Figure 11 shows the
estimated densities from the true $\kappa$’s (bold line) along with the ones
from the first 500 MCMC iterations (gray lines) for three different replicated
datasets, namely the 1st, 50th and 100th simulations. The proposed BYM2-Gamma
model seems to recover the true distribution of the $\kappa$’s. Regarding the
complete latent effects, $\bm{b}$, the estimated densities are again displayed
in Figure 12. The BYM2-Gamma model estimates accurately the distribution of
the latent effects.

Figure 9: WAIC across the 100 replicates for the proposed BYM2-Gamma model and
Congdon’s regarding the simulated data from the proposed BYM2-Gamma model.
Dashed lines: mean WAIC for each model Figure 10: Posterior summaries of the
parameters for the proposed BYM2-Gamma model across the 100 replicates
regarding the simulated data from the proposed BYM2-Gamma model.
Solid circle: posterior mean; Vertical lines: 95% posterior credible interval;
Solid horizontal line: true value.  Figure 11: Comparison between the
estimated density of the generated $\kappa$’s (solid black line) and the
estimated ones from the proposed BYM2-Gamma model for the 1st, 50th and 100th
replicates across the first 500 MCMC iterations (solid gray lines) Figure 12:
Comparison between the estimated density of the generated $b$’s (solid black
line) and the estimated ones from the proposed BYM2-Gamma model for the 1st,
50th and 100th replicates across the first 500 MCMC iterations (solid gray
lines)

##  Appendix D Simulation study: generating data from the proposed BYM2-logCAR
model

We now assess the proposed BYM2-logCAR model’s ability to recover the truth.
Similar to Appendix C, a simulation study is conducted wherein data are
generated from the proposed BYM2-logCAR model using the $n=160$ districts of
Rio de Janeiro. The unstructured and spatially structured components,
$\bm{\theta}$ and $\bm{u}^{\star}$ respectively, are independently generated
once, like in Appendix C. The scaling mixture components, $\bm{\kappa}$, are
generated once using the spatial structure as follows:

$\bm{z}\sim\mathcal{N}(\bm{0},\bm{Q}_{\star}^{-})\quad\mbox{and}\quad\kappa_{i}=\exp\left(-\frac{\nu_{\kappa}}{2}+\sqrt{\nu_{\kappa}}z_{i}\right),\
i=1,\dots,n,$

with an arbitrary $\nu_{\kappa}=0.3$ to allow the $\kappa$’s to depart from 1.
Like in Appendix C, the latent effects are then computed as
$b_{i}=\left[\sqrt{1-\lambda}\theta_{i}+\sqrt{\lambda}u_{i}^{\star}\right]\times\sigma/\sqrt{\kappa_{i}},\
i=1,\dots,n,$ where $\lambda=0.8$ and $\sigma=0.3$. Finally, 100 replicates of
populations of size $n=160$ are generated from the Poisson model,
$Y_{i}\sim\mathcal{P}\left(E_{i}\exp\left[\beta_{0}+b_{i}\right]\right),$ with
$\beta_{0}=-0.1$ and the offsets, $[E_{1},\dots,E_{n}]^{\top}$, taken from the
analysis of the Zika counts. The proposed BYM2-logCAR model and Congdon’s are
both fitted on the 100 replicated datasets using the same inference procedure
as in section 3.1.

Figure 13 shows that the WAIC always favours the proposed BYM2-logCAR model,
which generated the data. Figure 14 shows how well the proposed BYM2-logCAR
model is able to recover the true values of the model parameters through the
posterior summaries across the 100 replicates for the intercept, $\beta_{0}$,
the mixing parameter, $\lambda$, the hyperparameter, $\nu_{\kappa}$, and the
overall standard deviation, $\sigma$. Across the 100 replicates, the proposed
BYM2-logCAR model always captures the truth as the posterior 95% credible
intervals (vertical lines) always cover the true values of the parameters
(solid horizontal lines). Regarding the scaling mixture components,
$\bm{\kappa}$, Figure 15 shows the estimated densities from the true
$\kappa$’s (bold line) along with the ones from the first 500 MCMC iterations
(gray lines) for three different replicated datasets, namely the 1st, 50th and
100th simulations. Note that the choice of the prior mean
$\mu_{\nu_{\kappa}}=0.3$ for the hyperparameter $\nu_{\kappa}$ allows the
$\kappa$’s to depart from 1. This departure may be towards 0, to accommodate
potential outliers, as well as towards infinity. To ease visualisation, Figure
15 shows the estimated densities on a limited scale, cropping the $\kappa$’s
bigger than 8. The proposed model seems to recover the true distribution of
the $\kappa$’s. The proposed BYM2-logCAR model again estimates accurately the
distribution of the latent effects as shown in Figure 16.

Figure 13: WAIC across the 100 replicates for the proposed BYM2-logCAR model
and Congdon’s regarding the simulated data from the proposed BYM2-logCAR
model. Dashed lines: mean WAIC for each model Figure 14: Posterior summaries
of the parameters for the proposed BYM2-logCAR model across the 100 replicates
regarding the simulated data from the proposed BYM2-logCAR model.
Solid circle: posterior mean; Vertical lines: 95% posterior credible interval;
Solid horizontal line: true value.  Figure 15: Comparison between the
estimated density of the generated $\kappa$’s (solid black line) and the
estimated ones from the proposed BYM2-logCAR model for the 1st, 50th and 100th
replicates across the first 500 MCMC iterations (solid gray lines) Figure 16:
Comparison between the estimated density of the generated $b$’s (solid black
line) and the estimated ones from the proposed BYM2-logCAR model for the 1st,
50th and 100th replicates across the first 500 MCMC iterations (solid gray
lines)

##  Appendix E Simulation study: no outlying areas

To confirm that the proposed model does not detect outliers when unnecessary,
a simulation study is again conducted on the map of Rio de Janeiro without
contaminating any district. Data are generated 100 times as follows:

$Y_{i}\sim\mathcal{P}\left(E_{i}\exp[\beta_{0}+b_{i}]\right),\ i=1,\dots,n,$

with $n=160$, $\beta_{0}=-0.1$, $\bm{E}=[E_{1},\dots,E_{n}]^{\top}$ taken from
the Zika data analysis. The latent effects,
$\bm{b}=[b_{1},\dots,b_{n}]^{\top}$, are simulated once from a PCAR
distribution:

$\bm{b}\sim\mathcal{N}\left(\bm{0},\sigma_{b}^{2}\left[\bm{D}-\alpha\bm{W}\right]^{-1}\right),$

with $\sigma_{b}=\sqrt{0.2}$ and $\alpha=0.7$. Figure 17 shows the map of the
50th replicate of the simulated dataset, where no district seems to be an
outlier with respect to the whole city. Again, the two parametrisations of the
proposed model are compared to Congdon’s, using the same prior distributions
as described in section 3.1.

In terms of WAIC, the proposed models seem to perfom best, as shown in Figure
18. For this simulation study, the interest lies particularly in comparing the
outliers detections from the two versions of the proposed model and Congdon’s.
Figure 19 presents the districts that are found to be outliers by the
BYM2-Gamma proposed model (a), the BYM2-logCAR proposed model (b) and
Congdon’s (c). The BYM2-Gamma model only identifies one district, Freguesia on
the Governor’s island, to be a potential outlier in 2% of the replicates. The
BYM2-logCAR and Congdon’s models on the other hand detect Freguesia up to 8%
of the times, showing more sensitivity to the neighbourhood structure.
Congdon’s model further identifies 5 districts as potential outliers although
no district was contaminated.

Figure 17: Standardised mordibity ratio for the 50th simulation without
outliers. Figure 18: WAIC across the 100 replicates for the proposed models
and Congdon’s for the simulation without outliers. Dashed lines: mean WAIC for
each model

(a)  (b)

(c)

Figure 19:  Maps of the percentages of outliers as indicated by
$\kappa_{ur}<1$ across the $r=1,\dots,100$ replicates, where $\kappa_{ur}$ is
the upper bound of the posterior 95% credible interval of $\kappa$ in the
$r$th replicate of the simulated dataset without outliers. a) BYM2-Gamma
model; b) BYM2-logCAR model; c) Congdon’s model.