Bayesian Mapping of Mortality Clusters

In the analysis of environmental and epidemiological phenomena, researchers often gather geolocalised data in the form of multiple outcomes. For example, studies on the effects of long-term exposure to air pollution and high temperatures on public health generally collect and analyse both multiple environmental and health indicators (Künzli et al., 2000, Wong et al., 2008, Lee et al., 2009, Taylor et al., 2015). In epidemiological studies, several outcomes are used to determine the frailty characteristics of specific groups of patients (Li et al., 2021, Qian et al., 2023). Therefore, analysing the joint distribution of multiple outcomes and their interactions is a crucial step to uncover patterns in space and relationships that may not be evident when analysing each outcome on its own. This approach provides a more comprehensive understanding of the underlying environmental and epidemiological processes in a territory.

Disease mapping, in particular, aims at studying the spatial distribution of health phenomena to identify significant changes of diseases incidence or mortality. It therefore represents the first step for understanding geographical disparities in health outcomes and developing strategies for disease prevention and control. Disease mapping techniques have evolved significantly over the last decades, moving from simple representations of a disease occurrence to more advanced statistical analyses and models (Waller and Gotway, 2004, Lawson, 2013, Banerjee et al., 2015). The use of sophisticated modelling techniques yields a more detailed representation of the spatial distribution of phenomena in space and of the dependence across outcomes, resulting in more accurate explanatory and predictive performances (Aungkulanon et al., 2017, Coker et al., 2023, Griffith, 2023, Tesema et al., 2023). One of the main goals of disease mapping is discovering and locating groups of spatially proximate areas that show similar levels of the outcomes considered. We refer to these groups as spatial clusters. In public health, discovering spatial clusters of disease incidence or mortality is crucial, as it represents a step forward in the identification of possible underlying risk factors in specific areas of the territory (Kulldorff, 1997, Kulldorff et al., 2005, Robertson et al., 2010).

For instance, we consider a study of mortality levels in the Padua province, located in northeastern Italy, with data provided by the local healthcare unit ULSS6 Euganea. This institution oversees health services across the province, whose territory is divided into 106 administrative sub-areas (the city of Padua is not considered as a stand-alone municipality but is divided into six districts of size and population comparable to other municipalities in the province). The goal of ULSS6 is to optimise logistical and financial resources by identifying significant variations in mortality levels across either individual areas or clusters of areas. Their interest is particularly focused on the mortality due to three major death causes: diseases of the circulatory system, diseases of the respiratory system, and malignant neoplasms. According to Fedeli et al. (2021), these three categories collectively represent 70% of the overall mortality rate recorded in the province. Mortality data for these death causes are available from 2017 to 2019, divided by gender. The left panel of Figure 1 illustrates the province territory divided into the 106 districts managed by ULSS6. The left panel of Supplementary Figure 1 displays the geographical position of the Padua province in Italy. Bovo et al. (2023) made a first attempt to detect mortality clusters in this province by analysing each cause of death separately. Nonetheless, we recognise some important limitations in this approach since treating the variables independently can lead to identify clusters that vary substantially for each cause of death. Moreover, this method does not take full advantage of the information across variables, meaning the clustering of the territory driven by one cause of death lacks insight from the data of other causes. In contrast, gaining information from all causes of death simultaneously would lead to the identification of territorial clusters with specific mortality profiles, i.e., with estimated risk levels for each cause of death.

The case study of ULSS6 mortality data is relevant due to its depiction of a setting where areas have similar territorial characteristics, comparable population sizes and socio-economic conditions. On the other hand, in many other cases, the territory under consideration may be highly heterogeneous in one or more aspects. For instance, in the United States of America (U.S.) the administrative division of the territory into counties provides a far different scenario from the province of Padua. The distribution of population in space is heterogeneous, as well as population density, urbanisation and wealth; in addition, different U.S. states may have different health care systems, and this leads to another source of variability. The distribution of U.S. counties in space is displayed in the right panel of Figure 1. In this framework, the Centres for Disease Control and Prevention (CDC) collected data on mortality levels by counties in the U.S. from 2016 to 2019, available through the WONDER Online Databases (CDC, 2021). In such diverse contexts, it is crucial not only to quantify different risk levels across areas, but also to determine the existence of possible geographically concentrated clusters with specific epidemiological relevance.

Both the discussed examples reveal the need for a statistical method that can: i.) perform multivariate clustering of a territory based on the joint analyses of multiple outcomes, exploiting the spatial structure of data and the dependence across diseases; ii.) leverage information provided by exogenous variables that influence the outcomes; and iii.) identify which outcome variables lead to the formation of spatial clusters. This last point becomes more and more relevant as the number of outcomes increases or the size of the territory grows, making it essential to identify clusters of disease prevalence or excess mortality that are geographically concentrated in small areas while avoiding spurious signals.

In this article, we consider the problem of clustering areas of a territory based on the observed outcomes of multiple diseases. We assume that the territory is divided into polygons defined by administrative boundaries, and that the outcome variables are aggregates of measures taken within the boundaries of each area. In the ULSS6 Euganea case study presented in Section 1.1, the areas correspond to municipalities or administrative districts, while in the U.S. dataset the areas represent the U.S. counties. In both studies, the outcome variables are measurements of mortality from multiple causes of death.

A common approach for identifying groups of contiguous areas with excess of mortality, and thus detecting some spatial clusters, is through spatial scan statistics (SSS, Kulldorff, 1997). These techniques are designed to detect clusters with unexpected changes in outcome variable levels compared to the rest of the territory. This detection is carried out using suitable test statistics (Cucala, 2016, Liu et al., 2018, Ahmed and Genin, 2020). Over the past twenty years, these techniques have gained considerable attention and have been adapted for spatio-temporal data analysis in both frequentist (Kulldorff et al., 2005, Robertson et al., 2010, Amini et al., 2023) and Bayesian paradigms (Li et al., 2012). Recently, SSS have been largely extended also to the multivariate framework (Cucala et al., 2017, 2019) and to more complex data structures, such as multivariate functional data (Frévent et al., 2023). For a detailed review of applications of multivariate SSS to different contexts, readers can refer to the Introduction section of Frévent et al. (2023). Despite their extensive use, SSS are primarily designed to detect abrupt changes in outcome variables and thus they cannot be used to divide the entire map into clusters. In addition, SSS may fail to detect significant changes in areas with highly irregular shapes (Tango, 2021).

Another established approach for detecting spatial clusters is through hidden Potts models (HPM), introduced as an extension of the Gaussian mixture models (GMM) to account for spatial dependence across adjacent neighbours (Potts, 1952). In Bayesian literature, inference on HPM has been extensively explored using various approaches. Examples include Markov chain Monte Carlo (MCMC) methods (Green and Richardson, 2002), variational Bayes (McGrory et al., 2009), and approximate Bayesian computation (Moores et al., 2020a). Recent advancements also include estimation based on synthetic likelihood (Zhu and Fan, 2023). For a comprehensive review of HPM estimation methods, readers can refer to Moores et al. (2020b). Despite their popularity, the Potts model presents some relevant drawbacks, particularly related to the estimation of the inverse temperature parameter. Recently, Moores et al. (2020a) made a detailed survey of the approaches investigated so far for bypassing the estimation problems of HPM, with particular attention to Approximate Bayesian Computation algorithms, and proposed their own scalable solution. However, the current statistical literature does not adequately address the extension of HPM to the joint analysis of multiple outcomes. Additionally, to the best of our knowledge, there is no implementation of a multivariate HPM available, for instance, for the R programming language (R Core Team, 2024).

To address the need to identify mortality clusters jointly considering multiple causes of death, we propose perla (PEnalised Regression with Localities Aggregation), a multivariate Bayesian model for areal data that aggregates areas of a territory into spatially informed clusters based on different levels of outcome variables, while accounting for the presence of exogenous variables and for the marginal correlation across responses. Recently, Coker et al. (2023) applied a Dirichlet process mixture model to the analysis of areal data, without accounting for the spatial structure in the clustering phase. Instead, perla integrates the data spatial proximity directly into the clustering probabilities using the multinomial stick-breaking strategy of Linderman et al. (2015). This method consists of rewriting a multinomial model as a sequence of binomial random variables, allowing the use of the well-established Pólya-gamma data augmentation in the presence of multinomial outcomes. Not only does the multinomial stick-breaking simplifies the Bayesian inference allowing the conditional posterior distributions of the model parameters to be derived in closed form, but it also allows avoiding the use of the Potts model to account for the neighbourhood structure and consequently avoiding the known inferential problems associated with it. Alternative approaches to link a probability vector defined on a simplex to a series of Gaussian random variables include, for example, the logistic normal distribution (Lafferty and Blei, 2005, Russo et al., 2022). In addition, perla regulates the clusters intercepts through global-local shrinkage priors (Polson and Scott, 2011), allowing to affirm or inhibit the role of specific diseases in the formation of clusters. The shrinkage priors play also a regularising role in the selection of the number of clusters, as they help suppressing the effect of excess clusters. The reader can refer to Bhadra et al. (2019) and van Erp et al. (2019) for two reviews of global-local shrinkage priors.

This work is strongly motivated by the two epidemiological case studies presented in Section 1.1 and can be of great usefulness in the analysis of several other social and epidemiological phenomena, including the distribution of the level of education, the incidence of diseases, and the occurrence of multiple kind of crimes. Nonetheless, perla can be applied to any kind of multivariate areal data. Two examples of the numerous possible applications outside the epidemiological context are image segmentation (Moores et al., 2015) and spatial transcriptomics (Zhao et al., 2021). For these reasons, we implement our methodological solution in the software library perla for the R programming language.

The remainder of this article is structured as follows. Section 2 introduces the perla model, highlighting the formulation of the prior distributions and their role in the detection of spatial clusters. Section 3 delineates an MCMC algorithm to perform posterior inference and explains how to operate model selection and post-processing using the posterior samples. Section 4 presents three simulation experiments. The goal of the first experiment is to demonstrate the significant improvements in interpretability of the results and inferential conclusions provided by our global-local shrinkage priors compared to standard, non-informative priors. In addition, it performs a comparison of the clustering accuracy with other competing methods. The second experiment investigates whether different a priori choices on the parameters that handle the spatial correlation of the clustering probabilities have an impact on the final clustering performance, in particular when spatial clusters present different levels of spatial association. The third experiment evaluates the performance of the proposed methodology against competing methods using data simulated from a different generating mechanism that diverges from perla in several key aspects. Section 5 presents the application of perla on the two cases studies described in Section 1.1, showing the usefulness of our methodology in responding to specific epidemiological questions. Finally, Section 6 takes some concluding remarks and outlines future perspectives.

We assume that data are collected in the form of an $n\times d$ matrix $\mathbf{Y}$, where $y_{ij}\in\mathbb{R}$ denotes the mortality level due to the $j$-th disease ($j=1,\dots,d$) measured at the $i$-th area of a map ($i=1,\dots,n$). Data are structured such that $y_{ij}>0$ denotes a mortality excess, and $y_{ij}<0$ denotes a mortality deficit. In addition, let $\mathbf{X}$ be an $n\times p$ matrix collecting $p$ exogenous variables whose values are available for the $n$ areas. We assume that the $n$ areas can be divided into $K$ territorial clusters, implying that the average mortality levels in the $i$-th spatial area due to the $d$ diseases considered vary based on the cluster to which the $i$-th area belong. We collect the information about the unknown clusters into a $n\times K$ matrix $\mathbf{Z}$, where $\mathrm{Z}_{ik}=1$ if the $i$-th area belongs to the $k$-th cluster (with $k=1,\dots,K$), and 0 otherwise. To access the rows and the columns of a generic matrix of parameters $\boldsymbol{\theta}$, we use the notation $\boldsymbol{\theta}_{i.}$ and $\boldsymbol{\theta}_{.j}$ to refer to its $i$-th row and its $j$-th column, respectively. The same notation is used also for the data matrices $\mathbf{Y}$ and $\mathbf{X}$, with the only exception that their row and column vectors are denoted with bold lowercase letters.

We assume that regions within a spatial cluster have similar mortality levels for each disease considered, whereas significant variations in mortality levels occur across clusters. In addition, we assume that exogenous variables influence the outcomes but do not regulate the formation of clusters (Coker et al., 2023). Based on these considerations, we formulate the multivariate regression model

where $\mathbf{Z}_{i.}$ and $\mathbf{x}_{i.}$ are the $i$-th row of $\mathbf{Z}$ and $\mathbf{X}$, $\boldsymbol{\mu}=\{\mu_{kj}\}_{k=1,\dots,K;j=1,\dots,d}$ is a matrix containing the cluster- and disease-specific intercepts, $\boldsymbol{\beta}=\{\beta_{\ell j}\}_{\ell=1,\dots,p;j=1,\dots,d}$ is the matrix of regression coefficients, and $\boldsymbol{\epsilon}_{i.}$ is the vector of error terms. $\boldsymbol{\Sigma}$ aims at capturing the correlation across the $d$ diseases. The formulation of the regression model in (1) is equivalent to assuming a matrix variate Gaussian distribution on the outcome matrix $\mathbf{Y}$, with mean matrix $\mathbf{Z\boldsymbol{\mu}}+\mathbf{X\boldsymbol{\beta}}$, $\mathbb{I}_{n}$ as covariance matrix of the rows, and $\boldsymbol{\Sigma}$ as covariance matrix of the columns.

By definition of $\mathbf{Z}$, we define a cluster as a set of labelled areas with the same average levels of the $d$ responses, net of the effects of the covariates $\mathbf{X}$. Since the clustering matrix $\mathbf{Z}$ is unknown, we assume that $\mathbf{Z}_{i.}$ distributes according to a multinomial distribution of size 1 and vector of probabilities $\boldsymbol{\pi}_{i}=(\pi_{i1},\dots,\pi_{iK})$. To let the clustering draw information by the spatial structure of data, we first express the multinomial model as a series of binomial probability functions using the multinomial stick-breaking representation of Linderman et al. (2015)

where $N_{i1}=1$ (there is at most one element per row that is equal to one), $\tilde{\pi}_{i1}={\pi}_{i1}$, and, for $k>1$,

By construction, $\tilde{\pi}_{iK}=1$. Then, for $k=1,\dots,K-1$, we reparametrise the probabilities $\tilde{\pi}_{ik}$ using the logistic function, thus writing $\tilde{\pi}_{ik}=\exp(\psi_{ik})/\{1+\exp(\psi_{ik})\}$. The vector $\boldsymbol{\psi}_{.k}=(\psi_{1k},\dots,\psi_{nk})$ contains $n$ realisations in space of the transformed conditional probabilities of belonging to cluster $k$, given that the observations were not assigned to the previous $k-1$ clusters $(\boldsymbol{\psi}_{.1}$ simply represents the transformed probabilities of belonging to cluster 1). To introduce spatial dependence across the transformed clustering probabilities, we assume that each vector $\boldsymbol{\psi}_{.k}$ distributes according to a conditionally autoregressive (CAR) model of the form

where $\mathbf{W}$ has null diagonal and $\mathrm{W}_{ii^{\prime}}=1$ if regions $i$ and $i^{\prime}$ are neighbours, and 0 otherwise, $\mathbf{D}$ is a diagonal matrix with $\mathrm{D}_{ii}$ equal to the number of neighbours of $i$, $\tau_{k}>0$ and $\rho_{k}\in[0,1)$.

Although the model in (3) does not force the areas within the same cluster to be connected by a path, therefore allowing even spatially distant areas to belong to the same cluster, modelling the clustering probabilities with CAR distributions as in (3) encourage the identification of spatially separated clusters. The principal role of $\rho_{k}$ is ensuring that $\mathbf{D}-\rho_{k}\mathbf{W}$ is positive semidefinite. As discussed in Section 4.3.1 of Banerjee et al. (2015), $\rho_{k}$ cannot be interpreted as a measure of spatial correlation of the data, similarly for example to what is expressed by the Moran’s I statistic. The only exception is $\rho_{k}=0$, which denotes the case of independence across spatial areas. In addition, Assunção and Krainski (2009) highlighted that $\rho_{k}$ can assume also negative values, which however can still be associated to positive correlation values. These inconsistencies were underlined also by Ver Hoef et al. (2018), who stated that, in geostatistics, it is common to set $\rho_{k}$ to be positive, as negative values of spatial correlation are generally not considered. Therefore, following also Carlin and Banerjee (2003), in this work we assume $\rho_{k}$ to be only positive. To solve these oddities, Datta et al. (2019) proposed an alternative to the CAR formulation, named DAGAR, using a directed acyclic graph (DAG) to model the precision matrix of $\boldsymbol{\psi}_{.k}$. One of the key advantages of their approach lies in the interpretation of $\rho_{k}$, which now represents the degree of spatial correlation within the data. However, for this study, we have opted to retain the CAR formulation for mainly two reasons. Firstly, the DAGAR model necessitates the ordering of map areas. While Datta et al. (2019) demonstrate that ordering does not heavily impact the final results, it is not clear how the ordering operation should be performed when the statistical units, i.e., the spatial areas, are divided into clusters. In addition, the concept of ordering spatial data may pose challenges for non-experts in the field. Given our intention for the model to be used in epidemiological studies, we prefer to maintain a standard framework that does not rely on ordering. Secondly, our decision is reinforced by the favourable performance of the CAR model compared to DAGAR, as evidenced by Datta et al. (2019). Their study indicates that both models yield comparable results from various perspectives, with the interpretability of $\rho_{k}$ emerging as the primary advantage of DAGAR over CAR. In our application, we employ the spatial model solely to incorporate spatial correlation into the parameterised clustering probabilities $\boldsymbol{\psi}_{.k}$. Thus, we are not using a spatial model for the raw data itself, but rather for a latent parameter. Therefore, while interpretable, we believe that $\rho_{k}$ would offer a few epidemiological insights.

With the current formulation and assumptions, the stick-breaking construction offers a natural ordering of the discovered clusters. In fact, taking $\mathbb{E}(\boldsymbol{\psi}_{.k}|\tau_{k},\rho_{k})=\mathbf{0}$ is equivalent to assuming that the marginal expected value of $\tilde{\pi}_{ik}$ is $\mathbb{E}(\tilde{\pi}_{ik}|\tau_{k},\rho_{k})=0.5$, for $i=1,\dots,n$ and $k=1,\dots,K-1$. As a direct consequence, $\mathbb{E}({\pi}_{ik}|\tau_{k},\rho_{k})>\mathbb{E}({\pi}_{ik^{\prime}}|\tau_{k^{\prime}},\rho_{k^{\prime}})$ when $k>k^{\prime}$, therefore the average number of observations expected in the $k$-th cluster decreases as $k$ grows, and clusters are naturally interpreted from the most to the least predominant.

The model in Equation (1) can be rephrased also in terms of a multivariate Gaussian mixture model. In fact, it corresponds to assuming $\mathbf{y}_{i.}|\mathrm{Z}_{ik}=1\sim\mathcal{N}_{d}(\boldsymbol{\mu}_{k.}+\mathbf{x}^{T}_{i.}\boldsymbol{\beta},\boldsymbol{\Sigma})$, where $\Pr(\mathrm{Z}_{ik}=1)=\pi_{ik}$. Although, at a data level, perla and HPM are practically equivalent, they adopt distinct strategies for embedding the spatial structure of data into the clustering labels. In fact, HPM works directly with the conditional distributions $\mathrm{Z}_{ik}|\mathbf{Z}_{\mathscr{N}(i),k}$, where $\mathscr{N}(i)$ denotes the spatial neighbours of region $i$. Instead, perla assumes that the clustering probabilities are spatially correlated. In Section 3, we show how this strategy offers practical advantages in terms of Bayesian posterior inference.

In the Bayesian framework, the prior distributions of the model parameters are used either for conveying previous knowledge about the analysed phenomena or for specific modelling needs, like bounding the parameter space or performing variable selection. In this section, we discuss the a priori assumptions that we make to help the model retrieving the spatial structure of data, while inferring the underlying epidemiological phenomena.

perla is designed to disentangle spatial clusters based on the variations of disease-specific mortality in space. The elements in the Formula (1) that differentiate the clusters are the intercepts $\boldsymbol{\mu}=\{\mu_{kj}\}_{k=1,\dots,K,j=1,\dots,d}$. For instance, it is possible that not all considered diseases contribute to the formation of spatial clusters, or that certain diseases contribute to the formation of clusters only in specific areas of the map. In practice, we state that disease $j$ contributes to the formation of the $k$-th cluster if $\mu_{kj}$ deviates from zero, either positively or negatively. A positive deviation indicates a mortality excess in the cluster due to disease $j$, while a negative deviation suggests a mortality deficit. Therefore, we design an a priori setting for $\boldsymbol{\mu}$ that promotes the clustering only when it detects significant deviations from zero of the mean levels of the response variables across different areas of the map. For instance, suppose that in an epidemiological study comparing the mortality due to two diseases, the mean levels related to the first disease ($j=1$) are indicative the presence of two spatial clusters, while the mean levels related to the second disease ($j=2$) are close to zero. We aim for a prior distribution on the $\mu_{kj}$ that leverages the diversity between $\{\mu_{k1}\}_{k=1,2}$ to delineate the clusters, while shrinking the disparities between $\{\mu_{k2}\}_{k=1,2}$ toward zero, as they do not inform the clustering structure of the territory. We translate these concepts into the following global-local shrinkage prior:

where

and $\mathcal{C}_{+}$ denotes the half-Cauchy distribution. While $\phi$ represents a global shrinkage factor, shared among all the $\mu_{kj}$, $\zeta_{j}$ is a disease-specific parameter that is in common with the cluster intercepts $\mu_{1j},\dots,\mu_{Kj}$. A large value of $\zeta_{j}$ emphasises the contribution of outcome $j$ to the clustering process, and a small value of $\zeta_{j}$ diminishes its contribution. Lastly, the local factor $\gamma_{kj}$ acts specifically on the $kj$-th mean component, highlighting or shrinking the contribution of outcome $j$ to the identification of cluster $k$. We collect the shrinkage parameters appearing in (4) into the vectors $\boldsymbol{\zeta}=\{\zeta_{j}\}_{j=1,\dots,d}$ and $\boldsymbol{\gamma}=\{\gamma_{kj}\}_{k=1,\dots,K;j=1,\dots,d}$.

By introducing progressively more local factors in Formula (4), we model the influence of each disease on the overall clustering structure of data through $\boldsymbol{\zeta}$, and then we investigate the contribution of each disease to the formation of individual mortality clusters through the $\boldsymbol{\gamma}$. The local factors $\boldsymbol{\gamma}$ become particularly relevant when some diseases are informative only of a few, specific clusters. For instance, suppose that a disease $j$ is not informative of any cluster but the $k$-th. We expect the model to shrink the contribution of disease $j$ through a small value of $\zeta_{j}$, and to highlight the contribution of disease $j$ in identifying the $k$-th cluster through a large value of $\gamma_{kj}$. The total number of shrinkage parameters introduced with the prior model in (4) is $Kd+d+1$ . Without $\zeta_{j}$, our prior distribution reduces to the well-know horseshoe prior of Carvalho et al. (2010).

In alternative to (4), we consider also the prior distribution

The parameter $\sqrt{\delta_{k}}\sim\mathcal{C}_{+}(0,1)$ highlights or shrinks toward zero the differences across diseases within the whole $k$-th cluster, following a similar rationale of Denti et al. (2023). We collect all the shrinkage parameters introduced through Formula (5) into the vector $\boldsymbol{\delta}=\{\delta_{k}\}_{k=1,\dots,K}$. Not only does the prior in (5) simplify (4), but it is also more parsimonious than the horseshoe prior when $K,d\geq 2$, with at least one of them strictly greater than 2. This is because it requires only $K+d+1$ shrinkage parameters, compared to $Kd+1$ in the horseshoe prior. The assumption of separability between cluster and disease effects allows to interpret the whole contribution of disease $j$ in the formation of clusters through $\zeta_{j}$, and the amount of variability across diseases in cluster $k$ through $\delta_{k}$. In addition, despite not operating directly as a method for selecting the number of clusters, $\delta_{k}$ shrinks all the coefficients of the $k$-th cluster toward zero if the cluster finds poor evidence from the data, allowing for a clearer identification of the number of clusters in the post-processing phase. Finally, we can also obtain simpler priors taking only a subset of the shrinkage parameters in Formulas (4)-(5). For example, a prior having only the elements $\phi$ and $\boldsymbol{\zeta}$ applies a global shrinkage and a local disease-specific shrinkage.

The use of a half-Cauchy distribution could be potentially extended by allowing its scale parameter to differ from 1, thereby enabling the inferential process to incorporate some a priori knowledge about the analysed phenomena (Piironen and Vehtari, 2017). This extension can, for example, be used to specify which diseases are more likely to contribute to the clustering. However, we opted for the $\mathcal{C}_{+}(0,1)$ distribution because, a priori, there is often little or no knowledge of how many diseases are responsible of the formation of clusters. This is because, in mortality mapping, the primary question is whether clusters of mortality excess or deficit exist at all, with the number of clusters investigated subsequently (Bivand et al. 2013, Section 10.6). This differs fundamentally from other unsupervised analyses: for example, in transcriptomics, cells are clustered based on their gene expression levels and where there is typically a priori knowledge of both the number of cell types present in a tissue and of the genes known to be primarily linked to its biological processes (Zhao et al., 2021).

The parameters $\tau_{k}$ and $\rho_{k}$ calibrate the spatial distribution of $\boldsymbol{\psi}_{.k}$. In Supplementary Section 1, we discuss the roles of these two parameters on the clustering probability values. In summary, it emerges that $\tau_{k}$ has no impact on the mean of the values in $\boldsymbol{\psi}_{.k}$, but it has a large impact on their variability: the larger $\tau_{k}$, the greater the variability associated with each $\psi_{ik}$ (see for instance Supplementary Figure 4). The inverse gamma distribution is conjugate with the Gaussian model; therefore, a natural choice of a prior model for $\tau_{k}$ would be $\mathcal{IG}(\alpha_{\tau},\beta_{\tau})$. However, one could also avoid inferring an extra parameter and instead fixing $\tau_{k}$ sufficiently large to guarantee an adequate level of uncertainty in the clustering probabilities. In the simulation experiment proposed in Section 4.3, we compare different fixed values of $\tau_{k}$ with the use of its conjugate prior.

A second relevant aspect that emerges from the discussion in Supplementary Section 1 is that only values of $\rho_{k}$ close to 1 allow the CAR prior to handle spatially correlated data. The limit case $\rho_{k}=1$ is not admittable because $\mathbf{D}-\mathbf{W}$ is not a full-rank precision matrix, implying that the integral of the joint prior density of $\boldsymbol{\psi}_{.k}$ is not finite. Nonetheless, when the CAR model is used as a prior, Banerjee et al. (2015) encourage using the kernel of a CAR with $\rho_{k}=1$, trusting that the posterior distribution will be proper due to the contribution of the likelihood function. This solution is known in the literature as intrinsic CAR. However, in perla, the CAR serves not as a prior for a likelihood parameter, but for a parameter that defines the distribution of $\mathbf{Z}$, which is itself a random variable. Therefore, being in a second level of hierarchy in the model structure, we believe that the usage of an intrinsic CAR prior should be properly investigated and compared with a proper CAR prior that takes $\rho_{k}=0.99$ for every $k$. Alternatively, we propose also inferring $\rho_{k}$ from the data. Carlin and Banerjee (2003) suggest using $\rho_{k}\sim\mathcal{B}(18,2)$, where $\mathcal{B}$ denotes the Beta distribution, to give large probability to values of $\rho_{k}$ close to 1. It is worth noting that the density of this Beta distribution peaks at 0.944 and decreases rapidly as values approach 1, thus discouraging the limit case $\rho_{k}=1$. To construct a prior distribution for $\rho_{k}$, first it could be worth evaluating how the estimate of $\rho_{k}$ changes with respect to different levels of spatial correlation in the data. To do so, we simulated a series of spatial datasets using a DAGAR model, based on a map of $n=259$ U.S. counties on the West side of the U.S. map (see Supplementary Figure 2). The choice of the DAGAR as a data generative mechanism allows us to control the level of spatial correlation (called $\rho^{\mathrm{DAGAR}}$) in the simulations. For every level of $\rho^{\mathrm{DAGAR}}$ used as input, we produced 200 realisations from the spatial process, and we estimated the CAR model. Figure 2 reports the distribution of the maximum likelihood estimates of $\rho^{\mathrm{CAR}}$ under different $\rho^{\mathrm{DAGAR}}$. We notice that a clear correspondence between the estimate of $\rho^{\mathrm{CAR}}$ and the level of spatial correlation $\rho^{\mathrm{DAGAR}}$ does not emerge. In fact, when $\rho^{\mathrm{DAGAR}}\geq 0.4$, the distributions of the estimates of $\rho^{\mathrm{CAR}}$ concentrate above 0.9, while, when $0.1\leq\rho^{\mathrm{DAGAR}}\leq 0.3$, the estimates of $\rho^{\mathrm{CAR}}$ span over the whole $(0,1)$ interval. Lastly, when $\rho^{\mathrm{DAGAR}}=0.01$, the most likely value of $\rho^{\mathrm{CAR}}$ appears to be 0.01, although some uncertainty is also evident. Based on these considerations, a prior distribution for $\rho_{k}$ that favours values near 0 or 1 would be appropriate to capture a wide range of spatial correlation scenarios. We thus propose the following mixture prior:

We choose a symmetric structure to maintain non-informativeness about the value of $\rho_{k}$, giving equal probability to very small or very high values. The mixture component $\mathcal{B}(2,18)$ reflects the idea of the $\mathcal{B}(18,2)$ component, assigning high probability mass to small values of $\rho_{k}$, while discouraging the case of complete independence $\rho_{k}=0$, which is unusual in spatial data modelling. By placing most of the probability mass toward values close to zero or one - reflecting the almost absence or large presence of spatial correlation - the distribution in (6) retraces the idea of a continuous spike-and-slab prior (Ročková, 2018). Notice also that (6) represents a mixture model that integrates out the mixing probability $\theta_{k}$, whose prior distribution is assumed to be $\theta_{k}\sim\mathcal{B}(c_{\theta},c_{\theta})$, for any $c_{\theta}>0$. One of the goals of the simulation experiments proposed in Section 4.3 is comparing the two prior proposals for $\rho_{k}$, investigating whether they make any difference in terms of clustering accuracy and model fitting.

We complete the model specification with the prior distributions on $\boldsymbol{\beta}$ and $\boldsymbol{\Sigma}$. We generically assume the conjugate and poorly informative priors ${\beta}_{kj}\sim\mathcal{N}(0,10)$ and $\boldsymbol{\Sigma}\sim\mathcal{IW}(d,\mathbf{I}_{d})$, where $\mathcal{IW}$ denotes the inverse-Wishart distribution.

Figure 3 illustrates the hierarchical structure of perla through a DAG. Note that this scheme also includes the augmented variables $\boldsymbol{\omega}$, $\boldsymbol{\alpha}_{\phi}$, $\boldsymbol{\alpha}_{\zeta}$, $\boldsymbol{\alpha}_{\delta}$ and $\boldsymbol{\alpha}_{\gamma}$, which will be introduced in the next section to allow the posterior simulation of some model parameters to be performed via Gibbs sampling.

We carry out the inference on the model parameters simulating values from their posterior distributions. We present a Markov chain Monte Carlo scheme that allows a direct sampling from the conditional posterior distributions of most of the model parameters, requiring a Metropolis-Hastings step only for a few of them. This translates into the practical advantage of requiring fewer tuning parameters that can be hard to set when either the number of responses $d$ or the number of clusters $K$ grow. With respect to HPM, which requires particular attention for sampling the inverse temperature parameter, our method presents no critical step.

As we already mentioned, the stick-breaking representation in (2) allows to rewrite the multinomial model as a sequence of Bernoulli probability functions. To facilitate the posterior sampling, we introduce the latent variables $\omega_{ik}\sim\mathcal{PG}(N_{ik},0)$, for $k=1,\dots,K-1$ and $i=1,\dots,n$, where $\mathcal{PG}(\cdot,\cdot)$ denotes the Pólya-gamma random variable. Through this data augmentation strategy, we can update the transformed probabilities $\boldsymbol{\psi}_{.k}$ by generating from the conditional posterior distribution of $\boldsymbol{\psi}_{.k}|\omega_{k},\mathbf{Z}_{.k},\rho_{k}$, that is

where $\kappa_{ik}=\mathrm{Z}_{ik}-N_{ik}/2$. As the number of spatial areas $n$ grows, and with it the dimension of the neighbour matrix $\mathbf{W}$, it becomes more convenient to update the $\psi_{ik}$ iteratively, rather than generating the full vector $\boldsymbol{\psi}_{.k}$. This can be performed thanks to the Brook’s lemma, which allows to factorise the joint distribution in (3) into a series of univariate conditional distributions derived from the neighbourhood structure imposed by $\mathbf{W}$ (Besag, 1974). The conditional posterior of $\psi_{ik}$ is

where $\varsigma^{2}_{ik}={\tau_{k}}/\left({\tau_{k}\omega_{ik}+D_{ii}}\right)$. The augmented variables $\{\omega_{ik}\}_{i;k}$ are updated drawing from $\omega_{ik}|\mathbf{Z}_{i.},\psi_{ik}\sim\mathcal{PG}(\kappa_{ik},\psi_{ik})$ if $N_{ik}=1$, or setting $\omega_{ik}=0$ if $N_{ik}=0$.

The parameter $\rho_{k}$ under the prior model in (6) is the only case for which a closed-form updating scheme is not available. We implemented a random walk Metropolis-Hastings on the transformed parameter $\mathrm{logit}(\rho_{k})$, using a Gaussian distribution centred in the parameter value sampled at the previous iteration, and with variance $\tilde{\sigma}^{2}_{\rho}=3$, in the role of proposal.

To update the clustering labels, it is now sufficient to transform the quantities $\boldsymbol{\psi}_{.1},\dots,\boldsymbol{\psi}_{.K-1}$ into the clustering probabilities $\{{\pi}_{ik}\}_{i;k}$, then drawing a value from the multinomial distribution whose probabilities are defined as

where $\mathcal{N}_{d}(\mathbf{a}|\mathbf{b},\mathbf{C})$ denotes the density of the $d$-variate normal distribution with mean vector $\mathbf{b}$ and covariance matrix $\mathbf{C}$, evaluated in $\mathbf{a}$.

Under the prior model in (4), we update the shrinkage parameters $\boldsymbol{\zeta}$, $\boldsymbol{\gamma}$ and $\phi$ following the data augmentation strategy of Makalic and Schmidt (2016), which exploits the scale mixture representation of the half-Cauchy distribution. We thus introduce a set of augmented scale parameters for each shrinkage parameter considered: $\boldsymbol{\alpha}_{\zeta}=\{\alpha_{\zeta_{j}}\}_{j=1,\dots,d}$, $\boldsymbol{\alpha}_{\gamma}=\{\alpha_{\gamma_{kj}}\}_{k=1,\dots,K;j=1,\dots,d}$, and $\alpha_{\phi}$. Each of these augmented parameters is assumed to distribute according to an inverse gamma $\mathcal{IG}(1/2,1)$. Then, for $j=1,\dots,d$, draw from

For $k=1,\dots,K$ and $j=1,\dots,d$, draw from

Lastly, draw from

Let $\mathcal{C}_{k}=\{i=1,\dots,n:\mathrm{Z}_{ik}=1\}$ be the set denoting which observations belong to cluster $k$, and $n_{k}=|\mathcal{C}_{k}|$. The cluster intercepts are updated drawing from

where $\boldsymbol{\Sigma}^{*}_{k}=\left[\left\{\phi\;\mathrm{diag}(\zeta_{1}\gamma_{k1},\dots,\zeta_{d}\gamma_{kd})\right\}^{-1}+n_{k}\boldsymbol{\Sigma}^{-1}\right]^{-1}$.

If instead the prior model in (5) is considered, we introduce the set of augmented scale parameters $\boldsymbol{\alpha}_{\delta}=\{\alpha_{\delta_{k}}\}_{k=1,\dots,K}$, each of which distributes according to an inverse gamma $\mathcal{IG}(1/2,1)$. Then, for $k=1,\dots,K$, draw from

The rest of the updates can be performed as in Formulas (7), (9) and (10), just replacing $\gamma_{kj}$ with $\delta_{k}$.

Finally, the updates of $\boldsymbol{\Sigma}$ and $\boldsymbol{\beta}$ are straightforward to perform, thanks to the conjugancy between the likelihood and the selected priors. Details can be found for example in Marin and Robert (2007).

In our software implementation, we initialise the MCMC algorithm as follows. The $n$ observations are randomly assigned to the $K$ clusters; each element in $\boldsymbol{\psi}$ and $\boldsymbol{\beta}$ is generated from $\mathcal{N}(0,1)$; the Pólya-gamma parameters in $\boldsymbol{\omega}$ are drawn from $\mathcal{G}(1,1)$; and all shrinkage parameters are drawn from $\mathcal{C}_{+}(0,1)$. The cluster centroids $\boldsymbol{\mu}$ and the covariance matrix $\boldsymbol{\Sigma}$ are the only two parameters not randomly initialised. The former are taken equal to the ordinary least squares estimate $(\mathbf{Z}^{T}\mathbf{Z})^{-1}\mathbf{Z}^{T}\mathbf{Y}$, while the latter are initialised as $\sum_{i=1}^{n}(\mathbf{y}_{i.}-\overline{\mathbf{y}})^{T}(\mathbf{y}_{i.}-\overline{\mathbf{y}})/(n-1)$, where $\overline{\mathbf{y}}$ denotes the vector of sample means for the $d$ outcome variables.

To provide point estimates of the model parameters, we consider the posterior mode. With the notation $\mathcal{M}(\theta|\mathrm{data})$, we denote the posterior mode of a generic model parameter $\theta$, given the data. A relevant quantity in the analysis of mortality data is also the posterior probability of observing an excess mortality in a specific region, compared to the whole territory. Following the considerations of Richardson et al. (2004), we state that cluster $k$ shows evidence of mortality excess due to disease $j$ if $\Pr(\mu_{kj}>0|\mathrm{data})\geq 0.95$, and a deficit of mortality if $\Pr(\mu_{kj}>0|\mathrm{data})\leq 0.05$. These quantities can be computed straightforwardly from the output of the MCMC algorithm.

The presented MCMC scheme is susceptible to the common label-switching phenomenon, which arises due to the unidentifiability of the mixture components and leads the posterior distributions of $\boldsymbol{\mu}$ to be multimodal. In practice, this makes inference on model parameters challenging. To disentangle the mixture components, we explore the ECR algorithm of Papastamoulis and Iliopoulos (2013), considering the partition associated to the largest log-likelihood value as the pivotal one. In addition to this method, the R package label.switching (Papastamoulis, 2016) implements numerous alternative algorithms. Not only are these tools useful to relabel the posterior draws, but they are also designed to return an optimal partition of data based on the MCMC output. This process can also be accomplished using other algorithms, such as the SALSO greedy search method proposed by Dahl et al. (2022).

perla can be fitted with different setups, varying the number of clusters $K$ and the prior distributions of the model parameters. Therefore, a criterion that allows users to evaluate different model setups and determine which is preferable for the analysed data is necessary. This approach aligns with the philosophy of parsimonious clustering (Celeux and Govaert, 1995, Scrucca et al., 2016), which involves comparing the same type of model with different degrees of flexibility and determining the setup that best fits the data.

We conduct the model selection using the $\mathrm{DIC}_{3}$ information criterion of Celeux et al. (2006), which is recommended for mixture models, and whose properties have been studied by Kim (2021). A detailed description on how to compute $\mathrm{DIC}_{3}$ for perla is given in Supplementary Section 2. Models that return smaller $\mathrm{DIC}_{3}$ values are preferable to those producing larger values.

We already discussed in Section 2.2 that different setups for $\boldsymbol{\mu}$ imply different numbers of model parameters. Assuming that the global shrinkage parameter $\phi$ is always inserted into the model, we report in Table 1 the possible combinations of shrinkage factors that we consider for perla, along with the corresponding code notation used in our R package perla. Given the model dimension, users can evaluate all the five combinations of shrinkage factors shown in the figure, starting with simpler models and increasing the complexity.

Lastly, it is crucial to underline that model selection cannot rely solely on an information criterion, but must be performed considering also other factors, such as the quality of model fitting, evaluated through the MCMC output, and the interpretability of the clusters.

The data generating mechanism that we adopt follows the scheme of perla displayed in Figure 3, but using the DAGAR prior in place of the CAR to generate the clustering probabilities. This choice allows us to control the amount of spatial correlation simulated in each cluster. The use of the DAGAR requires the ordering of the areas. We thus order the counties from the most southern to the most northern. If $i=1,\dots,259$ is the index denoting the counties, then $i=1$ corresponds to the most southern county, and $i=259$ to the most northern. Being $K^{\text{true}}$ the true number of clusters in the map, the simulation model generates the transformed clustering probabilities as

where $k=1,\dots,K^{\text{true}}-1$, $\mathscr{N}(i)$ denotes the set of spatial areas that are neighbours of $i$ and precede $i$ according to the imposed ordering, and $n_{<i}=|\mathscr{N}(i)|$. The quantity $\rho^{\text{true}}_{k}\in(0,1)$ denotes the amount of spatial correlation of the probabilities of belonging to cluster $k$. For further details on the DAGAR model, the reader can refer to Datta et al. (2019).

By transforming the quantities $\{\psi^{\text{true}}_{ik}\}_{i=1,\dots,259;k=1,\dots,K^{\text{true}}-1}$ into the clustering probabilities $\{\pi^{\text{true}}_{ik}\}_{i=1,\dots,259;k=1,\dots,K^{\text{true}}}$ through the inverse transformation of the stick-breaking weights, we would obtain highly unbalanced clustering probabilities. To simulate approximately balanced clusters, the inverse transformation of the stick-breaking weights must be applied to $\{\psi^{\text{true}}_{ik}-\log(K-k)\}_{i=1,\dots,259;k=1,\dots,K^{\text{true}}-1}$. Then, the model allocates the spatial areas into clusters, storing the labels into the matrix $\mathbf{Z}^{\text{true}}$. Lastly, given the matrix of centroids $\boldsymbol{\mu}^{\text{true}}$ and the covariance matrix of diseases $\boldsymbol{\Sigma}^{\text{true}}$, the model generates the observed data in the form of the logarithm of mortality rates of $d$ diseases (named also death causes) as

In all the proposed simulation experiments, we take our conclusions using 20 replicates of the data generated under the same experimental conditions.

The first experiment serves a dual purpose. On one hand, it compares perla with some competing clustering models, such as mixture models and SSS. On the other hand, it examines the gain brought by the use of shrinkage priors, such as those appearing in Formulas (4) and (5), in posterior inference compared to using non-informative priors when the number of diseases considered in the analysis is large.

To accomplish this, we propose a scenario where the formation of mortality clusters is driven only by a subset of the diseases considered. We use $d=10$ diseases and $K^{\text{true}}=3$ clusters. For each of the 20 simulated datasets, we induce large spatial correlation across clustering probabilities, taking ${\rho}^{\text{true}}_{k}\sim\mathcal{U}(0.8,1)$, for $k=1,2$. Next, we randomly determine the diseases that do not contribute to cluster formation: for $j=1,\dots,10$, we introduce a random variable $I_{j}\sim\mathcal{B}e(0.5)$, which takes value 1 if disease $j$ informs the clustering, and 0 if it doesn’t. Then, we draw $\mu^{\text{true}}_{kj}|I_{j}=1\sim\mathcal{U}(-0.5,0.5)$, while we set $\mu_{kj}=0$ if $I_{j}=0$, for all $k$. To generate the correlation matrix of the diseases $d$, we draw values from the uniform distribution $\mathcal{U}(-0.1,0.1)$ to generate the correlations across diseases, which are used to fill the off-diagonal elements of an identity matrix of order $d$. Then, we check that the resulting matrix is semi-positive definite. Finally, we multiply all the elements of the correlation matrix by a factor of 0.05 to obtain the covariance matrix $\boldsymbol{\Sigma}^{\text{true}}$. The 20 synthetic datasets are finally generated following the simulation model described in Section 4.1. The left panel of Figure 4 illustrates the distribution of the outcome variables in one of the 20 datasets, where $I_{j}=0$ for $j=1,\dots,5$ and $j=7$.

For each of the 20 scenarios, we fit the following models:

• •

  perla, with shrinkage factors $\phi\zeta$, $\phi\delta$, $\phi\gamma$, $\phi\zeta\delta$, and $\phi\zeta\gamma$. To recall them, we use the notation reported in the right column of Table 1;

• •

  perla, with the non-informative priors $\mu_{kj}\sim\mathcal{N}(0,10)$, for all $k$ and $j$;

• •

  a Gaussian mixture model (GMM), selected with the BIC criterion, using the R library mclust (Scrucca et al., 2016);

• •

  the k-means algorithm, considered also by Aungkulanon et al. (2017);

• •

  the SSS of Gómez-Rubio et al. (2019), implemented in the R library Dclusterm. This algorithm fits a series of generalised linear models using the observed death counts as outcomes, and the expected death counts as offsets. However, our simulation model produces continuous data. Therefore, we applied Dclusterm within a linear regression model, rather than a Poisson regression. In addition, this method does not handle the analysis of multiple causes of death simultaneously, so it needs to be applied to every disease separately.

To the best of our knowledge, a library software that implements the multivariate HPM is not available in the R computing environment. In addition, the current version of bayesImageS (Moores et al., 2020a) does not fully support non-regular lattice data. This does not allow us to compare perla with HPM.

We apply perla, GMM and k-means using $K=4$ to test how effectively each model can retrieve the clustering structure of data when the number of clusters is incorrectly specified. The right panel of Figure 4 displays the adjusted Rand index (ARI, Hubert and Arabie 1985) distribution on the 20 datasets given by each model. We notice that all the configurations of perla lead to comparable results in terms of classification accuracy. This result can be explained by the presence of some $\mu^{\text{true}}_{kj}$ whose value is close to 0 for some $k$ even when $I_{j}=1$. For example, in the left plot of Figure 4 we see that the distribution of $\{y_{i6}:\mathrm{Z}^{\text{true}}_{i}=3\}$, the death cause 6 in cluster 3, is much closer to zero than the distributions in the two other clusters. Therefore, in this simulation context, it is reasonable that different combinations of the shrinkage parameters can be equally efficient in retrieving the clustering structure. As expected, all the configurations of perla globally perform better than mixture models not informed by the spatial data structure. We do not report the results obtained with the Dclusterm algorithm, as it failed to retrieve the clustering structure in all of the 20 datasets.

Despite the ARI values obtained with different perla setups are comparable, we now show that the contribute brought by different shrinkage priors extends beyond the overall clustering accuracy and consists of retrieving which diseases are responsible of the formation of mortality clusters. We display in Figure 5 the posterior distributions of the $\mu_{kj}$ parameters across the 20 simulated maps for which $I_{j}=0$. We compare the results obtained from the five penalised versions of perla, and perla with non-informative priors (denoted ‘Perla NI’ in the figures). Ideally, since the corresponding $\mu^{\text{true}}_{kj}$ are zero, the posterior distributions should be concentrated around zero. Figure 5a displays the posterior modes and the logarithm of the 95% HPD interval lengths. From top-left panel, it emerges that ‘Perla NI’ provides highly variable posterior estimates of the $\mu_{kj}$, and large HPD intervals. This is also evident from the marginal distributions of the posterior modes and interval lengths, displayed in Figures 5b and 5c. We notice that, among the estimates provided by ‘Perla NI’, some parameters show very large HPD interval sizes. These parameters are associated to the fourth mixture component, which does not find a correspondence with a real cluster. This effect is no longer present when using perla with only the cluster-specific shrinkage factors (setup (c)), as seen in the top-central panel of Figure 5a, although the posterior modes remain generally variable. With the setups (cd), (d), and (c,d), many posterior modes are shrunk toward zero, while a few of them remain insufficiently shrunk, and the HPD interval lengths decrease. In addition, note from Figure 5c that the HPD intervals of perla (cd), which uses the horseshoe prior, remain wider that the ones provided by the two other setups. Finally, our proposed prior (4), implemented in (d,cd) effectively shrinks all posterior modes towards zero, with narrow HPD intervals.

Figure 6 displays the posterior distributions of the $\mu_{kj}$ parameters for cases with $I_{j}=1$, corresponding to diseases that effectively contribute to the formation of mortality clusters. Since $\mu^{\text{true}}_{kj}\mid I_{j}=1\sim\mathcal{U}(-0.5,0.5)$, some coefficients may be generated very close to zero, and we expect our shrinkage priors to estimate them as null. This is observed in Figure 6a, where several posterior distributions have mode on zero, especially under the priors (cd) and (d,cd). The reason is that these setups employ the shrinkage factor $\gamma_{kj}$ to act directly on $\mu_{kj}$, without affecting the other elements in the same row, $\boldsymbol{\mu}_{k\cdot}$, and column, $\boldsymbol{\mu}_{\cdot j}$. In contrast, the setup (d) applies the same amount of shrinkage to all the coefficients of the same disease, while (c,d) applies to $\mu_{kj}$ a combination of the shrinkage factors $\delta_{k}$ and $\zeta_{j}$, which account for all the other elements in $\boldsymbol{\mu}_{k.}$ and $\boldsymbol{\mu}_{.j}$. Finally, the version (c) effectively shrinks to zero all coefficients that do not match a real cluster, which is something that ‘Perla NI’ fails to achieve. Overall, the marginal distributions of the posterior modes and interval lengths from the five penalised versions of perla appear similar (see Figures 6b and 6c).

We conclude that clustering accuracy is not the only metric to consider for evaluating the model performance. Although perla provides stable results in terms of ARI using different prior setups, our results demonstrate that our shrinkage priors improve interpretability by indicating which diseases do, and which do not, contribute to cluster formation. If disease $j$ does not affect the overall clustering, a prior setup that uses just its specific shrinkage parameter $\zeta_{j}$ is sufficient to shrink the means $\boldsymbol{\mu}_{.j}$ towards zero while recovering the clustering using information from informative diseases. Differently, if disease $j$ is informative for only a subset of the $K$ clusters, then the model requires the use of additional cluster-disease shrinkage parameters $\gamma_{kj}$. Our results also indicate that combining disease-specific shrinkage factors $\zeta_{j}$ with cluster-disease shrinkage factors $\gamma_{kj}$ yields better performance in distinguishing between informative and uninformative diseases, compared to using cluster-disease shrinkage factors alone. Finally, the cluster-specific factors $\delta_{k}$ efficiently shrink additional mean parameters, a feature that is particularly useful when the number of mixture components is taken larger than the number of clusters in the data.

Model selection can be guided by the $\mathrm{DIC}_{3}$ criterion, as discussed in Section 3.2. In our practical applications presented in Section 5, we select the model based on both the $\mathrm{DIC}_{3}$ and an examination of the posterior distributions of the $\{\mu_{kj}\}_{k;j}$, evaluating which models produce well-separated mixture components that can be used to describe the underlying epidemiological phenomena.

The second simulation experiment aims to study the performance of perla when the clustering probabilities have different levels of spatial correlation. This setting reflects a scenario where some clusters exhibit strong spatial connection, resulting in geographically concentrated and interconnected areas, while others may show weak spatial connection, with areas scattered across the whole territory. Since perla accounts for the spatial structure of the data through the CAR priors in (3), an additional scope of this experiment is investigating the effects of different prior configurations for the parameters $\tau_{k}$ and $\rho_{k}$.

To accomplish this, we generate 20 simulated maps of mortality rates, considering $d=3$ death causes and $K^{\text{true}}=4$ clusters. To induce different levels of spatial correlation among clustering probabilities, we consider three equispaced values of spatial correlation $\boldsymbol{\rho}^{\text{true}}=(0.01,0.455,0.9)$, that remain fixed across the 20 replicates. To induce marginal correlation across causes of death, we generate $\boldsymbol{\Sigma}^{\text{true}}$ as described in Section 4.2: correlation values are drawn from a uniform distribution $\mathcal{U}(-0.2,0.2)$. After ensuring the semi-positive definiteness of the correlation matrix, we scale all its elements by a factor of 0.07 to obtain the covariance matrix. We generate the cluster- and disease-specific intercepts $\mu_{kj}$, for $k=1,\dots,4$ and $j=1,2,3$, from $\mathcal{U}(-1,1)$. Lastly, each dataset is drawn from the simulation model detailed in Section 4.1. Supplementary Figure 5 displays, on the left panel, the spatial clusters in one of the 20 simulated scenarios and, in the right panel, the distribution of the three outcomes within each cluster appearing in the left panel. The map shows that counties in cluster 1 are spread throughout the territory, while those in clusters 2 and 3 are more geographically concentrated due to the larger spatial correlation values employed. From the spatial shape of the clusters, we highlight the importance of using a model that leverages spatial information without strictly enforcing clusters to be made solely of contiguous areas. For example, cluster 3 (shown in green in the left panel) consists of separate agglomerates of contiguous areas that are not directly connected. It is essential that the clustering model recognises these separate agglomerates as part of a unique cluster, since they share the same epidemiological characteristics.

First, we study the effects of different a priori assumptions on $\rho_{k}$, keeping $\tau_{k}=1$ for all $k$. For this experiment, we consider three prior setups: i.) $\rho_{k}=0.99$ $\forall k$, that is a proper CAR model with large spatial correlation; ii.) $\rho_{k}=1$ $\forall k$, that is an improper CAR; and iii.) the prior given in Formula (6) (denoted as ‘S&S’ in the figures). We evaluate whether allowing $\rho_{k}$ to be learnt from the data offers advantages in clustering accuracy and model fitting compared to assuming a constant value, especially in a framework where clusters are characterised by different levels of spatial correlation. In addition, we check whether assuming a proper or an improper CAR model yields meaningful differences in performance. We fit only perla (d,cd), as the evaluation of different prior models for $\mu_{kj}$ is already treated in Section 4.2.

The left panel of Figure 7 displays the ARI index distribution based on 20 simulated datasets. We observe that, when the number of clusters is correctly specified, perla accurately recovers the underlying clustering structure (median ARI $>0.9$). As expected, the performance is generally worse when the number of clusters is misspecified (median ARI index $\in(0.85,0.9)$). All the prior setups perform similarly in terms of clustering accuracy for $K=3$ and $K=4$. Since clusters are characterised by different levels of spatial correlation, we also measure the accuracy in recovering each cluster separately. To do so, let $\mathbf{Z}^{\mathrm{post}}$ be the posterior estimate of the matrix $\mathbf{Z}$, and $h(k)=\arg\max_{\ell=1,\dots,K}\mathrm{ARI}(\mathrm{Z}^{\text{true}}_{.k},\mathrm{Z}^{\mathrm{post}}_{.\ell})$ be the index of the estimated cluster that is most associated with the $k$-th real cluster. The top line of Figure 8 displays the ARI distribution between $\mathrm{Z}^{\text{true}}_{.k}$ and $\mathrm{Z}^{\mathrm{post}}_{.h(k)}$, for each true cluster $k=1,\dots,K^{\text{true}}$. These graphs help to visualise which clusters are more challenging to be recovered. When $K=3$, the third cluster remains the most difficult to retrieve, and the model employing the spike-and-slab-type prior ‘S&S’ performs worse than the models with assume a fixed $\rho_{k}$. When $K=4$, we do not observe significant differences across clusters or models.

Last, we test if using different prior setups for $\rho_{k}$ has an impact or not on model fitting. To do so, we display in the second row of Figure 8 the mean squared error (MSE) distribution of each real cluster, that is computed as $\sum_{i:\mathrm{Z}^{\text{true}}_{ik}=1}(y_{ij}-\mathbf{Z}_{i.}^{\mathrm{post}}\boldsymbol{\mu}^{\mathrm{post}}_{.j})^{2}$, where $\boldsymbol{\mu}^{\mathrm{post}}_{.j}$ denotes the posterior estimate of $\boldsymbol{\mu}_{.j}$. When the number of clusters is misspecified, all models struggle to fit the data in cluster 3. However, the model denoted ‘S&S’ consistently performs poorly in all clusters and causes of death compared to the CAR priors with fixed $\rho_{k}$ values, which perform similarly. This trend is also observed when $K=4$.

The second step involves studying the effects of different prior models on $\tau_{k}$, keeping $\rho_{k}=1$ $\forall k$. For this experiment, we assume $\tau_{k}=\tau$ $\forall k$, and we investigate the following options: i.) fixing $\tau$ respectively to 0.1, 1, and 5, ii.) allowing $\tau$ to be learned from the data, assuming the conjugate prior $\mathcal{IG}(2.1,3.1)$, which implies $\mathcal{M}(\tau)=1$ and a 95% HPD interval that ranges approximately between 0.3 and 7.85.

The right panel of Figure 7 displays the distribution of the ARI across the 20 simulated datasets. We observe that any choice of $\tau$ leads essentially to equivalent results in terms of the overall classification accuracy. We display also in the third row of Figure 8 the distribution of the ARI between $\mathrm{Z}^{\text{true}}_{.k}$ and $\mathrm{Z}^{\mathrm{post}}_{.h(k)}$, for each true cluster $k=1,\dots,K^{\text{true}}$. When $K=3$, the third cluster remains the most difficult to recover, and the model assuming $\tau\sim\mathcal{IG}$ performs slightly worse in recovering it than the models with fixed $\tau$. When $K=4$, we do not observe significant differences across clusters or across models.

Finally, the last row of Figure 8 displays the distribution of the mean squared error (MSE) within each true cluster. The third cluster remains the most difficult to recover when $K=3$, while the quality of the model fitting remains stable across every cluster when $K=4$.

We conclude that, when the data-generating mechanism presents several clusters with different levels of spatial correlation, perla performs well in recovering the clustering structure when the number of clusters is correctly specified. In addition, it appears not to be influenced by the varying levels of spatial correlation. On the contrary, when the assumed number of clusters is lower than the true number, both the clustering accuracy and the model fitting are more affected in areas with high spatial correlation. This experiment has served also to frame the impact of the hyperparameters values of the CAR prior on the model performance. Our simulations show that specifying a prior distribution on $\rho_{k}$ or fixing it at a constant value yield equivalent results in terms of clustering accuracy, while the model fitting is poorer when $\rho_{k}$ is not taken as constant. In addition, we have observed that assuming an improper CAR ($\rho_{k}=1$) is practically equivalent to assuming a proper CAR with a sufficiently large $\rho_{k}$. Finally, our model appears to be robust to different choices of $\tau_{k}$, producing practically equivalent results whether $\tau_{k}$ is assumed variable or fixed in advance.

In epidemiological studies, data are typically collected as counts. For instance, mortality studies rely on the number of deaths attributed to various diseases within each area of the analysed territory. However, these raw counts are not meaningful per se unless compared with an appropriate reference measure. For example, the observed number of deaths can be assessed against the expected number in the territory, which accounts for population stratification by age and allows to identify areas with mortality excess or deficit.

From a statistical modelling perspective, Poisson regression is a natural choice for count data, with the chosen reference measure incorporated as an offset. In practice, however, while raw counts are often available, suitable comparison measures are frequently inaccessible due to privacy restrictions. This limitation motivated us to develop a model that works with continuous outcome variables, such as rates, which are both readily available and unaffected by privacy constraints. Nonetheless, we believe it essential to assess our model’s performance in scenarios where both observed counts and a valid comparison measure are available.

In this section, we propose an experiment that stresses perla by simulating the data using a generating mechanism that substantially deviates from our model in three aspects: first, it uses Formula (11) to simulate the log-risks, incorporating also three covariates into the generative mechanism; then, it simulates the observed deaths from a Poisson model; and last, it assumes that different diseases yield different spatial clustering patterns. We consider $d=5$ causes of death. The first two diseases yield two clusters, obtained by dividing the map along a line at longitude 117°8’2.4” W. This line lies at the midpoint of the latitude range, so one cluster occupies the northern part of the map and the other the southern part. A representation of these two clusters is provided in the left panel of Supplementary Figure 6. We store the clustering information into the $n\times 2$ matrix ${\mathbf{Z}}^{(1:2)^{\text{true}}}$, and the mean levels of the two diseases observed in the two clusters into the $2\times 2$ matrix ${\boldsymbol{\mu}^{(1:2)}}^{\text{true}}$. The next two diseases yield two clusters, obtained by dividing the map along a line at latitude 40°10’33.888” N. As this line lies at the midpoint of the longitude range, one cluster occupies the western part of the territory and the other the eastern part. A representation of these two clusters is shown in the central panel of Supplementary Figure 6. We store the clustering information into the $n\times 2$ matrix ${\mathbf{Z}}^{(3:4)^{\text{true}}}$, and the mean levels into the $2\times 2$ matrix ${\boldsymbol{\mu}^{(3:4)}}^{\text{true}}$. The last disease does not show any clustering pattern.

The proposed data-generating scheme simulates the log-risks as

The independent variables are generated as $x_{i\ell}\sim\mathcal{N}(-0.2,\,0.5^{2})$ and the regression coefficients as $\beta_{\ell j}\sim\mathcal{N}(0.02,\,0.05^{2})$, for $\ell=1,2,3$. These values were manually tuned to avoid unrealistically high mortality risks. The covariance matrix $\boldsymbol{\Sigma}^{\text{true}}$ is assumed to be block diagonal, with the blocks $\boldsymbol{\Sigma}^{\text{true}}_{1:2,1:2}$ and $\boldsymbol{\Sigma}^{\text{true}}_{3:4,3:4}$ generated as in Simulation 1, and with $\Sigma^{\text{true}}_{55}=0.05$.

To generate the observed deaths, we consider the expected ones as offset. These are obtained from $E_{ij}\sim\mathcal{G}(1,0.03)$. This distribution was chosen because it assumes that the average number of expected deaths is 30 and approximately the 95% of the areas in the territory have at most 100 deaths, while the remaining 5% can have much higher values, representing a few areas with large population sizes (an consequently more deaths). The observed counts in the $i$-th area of the map and for the $j$-th cause of death are simulated from $O_{ij}\sim{\mathcal{P}oi}(\theta^{\text{true}}_{ij}E_{ij})$. Finally, to analyse these data using perla, we compute the log-SMR as $y_{ij}=\log(O_{ij}/E_{ij})$ or, in case $O_{ij}=0$, $y_{ij}=\log\{(O_{ij}+1/2)/E_{ij}\}$, as suggested by Bivand et al. (2013).

We generate 20 replicates of the experiment, and we apply perla (d,cd) to ensure flexibility, incorporating the model covariates. For comparison, we apply the SSS of Gómez-Rubio et al. (2019), already considered in Section 4.2, that must be applied to each outcome separately and works directly with counts data $O_{ij}$. In order to obtain a final clustering, we combine the labels obtained on the five outcomes into a unique clustering configuration. In addition, we consider also the Gaussian mixture model that assumes $\mathbf{y}_{i\cdot}\mid\mathrm{Z}_{ik}=1\sim\mathcal{N}_{d}\Bigl{(}\boldsymbol{\mu}_{k\cdot}+\mathbf{x}_{i\cdot}^{T}\boldsymbol{\beta},\,\boldsymbol{\Sigma}_{k}\Bigr{)}$, with $\Pr(\mathrm{Z}_{ik}=1)=\pi_{k}$ and the covariance structure $\boldsymbol{\Sigma}_{k}$ selected using the BIC criterion. Similarly to perla, this model assumes that the covariates are not directly related to the clustering. All the model parameters are estimated via maximum likelihood.

Due to the complex clustering structure used to generate these data, evaluating the clustering accuracy of the models is not immediate. By combining ${\mathbf{Z}^{(1:2)}}^{\text{true}}$ and ${\mathbf{Z}^{(3:4)}}^{\text{true}}$, we obtain a partition of the territory into four clusters (see the right panel of Supplementary Figure 6) that we use to evaluate the clustering models. We fit perla and GMM using $K=2,3,4$. The ARI distributions in the left panel of Figure 9 show that the highest clustering quality is achieved by perla with $K=4$. This result is particularly remarkable considering that the SSS is works directly with raw counts data. We notice also that both perla and GMM perform better when the number of clusters is enlarged. This is due to the nested clustering structure that characterises this simulation; the gain in the ARI level is particularly visible moving from 2 to 3 clusters. We conclude that, when different the clusters vary with the outcome variables, an effective strategy is to increase the number of clusters assumed by the model. In addition, we observe that transforming the outcomes into continuous variables rather that using directly the counts data does not lead to any loss in clustering accuracy.

In addition, we report in the right panel of Figure 9 the distribution of the computational cost of perla as a function of the number of clusters. Each model estimation involves running four Markov chains of length 10,000 sequentially using the algorithm detailed in Section 3.1 and combining the results with the post-processing algorithm described in Section 3.2. Consequently, the runtime for a single Markov chain of length 10,000 can be roughly deduced by dividing the values in the graph by four.

We investigate the presence of mortality clusters in Padua province using age-adjusted standardised mortality ratios $SMR_{ij}=O_{ij}/E_{ij}$, derived as the ratio between observed death counts ($O_{ij}$) and expected death counts ($E_{ij}$) in the $i=1,\dots,106$ areas and for the $j=1,2,3$ diseases considered. The expected number of deaths represents the number of deaths per age class that would occur if the mortality rates per age class of the entire province were applied. Thus, assuming $T$ age classes, the expected number of deaths from disease $j$ in the $i$-th area is $E_{ij}=\sum_{t=1}^{T}n_{it}r_{jt}$, where $r_{jt}$ is the provincial mortality rate of the $j$-th disease of interest within the $t$-th age class, measured in the $i$-th area, and $n_{it}$ is the population size within the $i$-th area and the $t$-th age class. For our analysis, we consider

the logarithm is applied so that $y_{ij}=0$ represents the case of concordance between observed and expected deaths, $y_{ij}>0$ denotes a mortality excess compared to the province level, and $y_{ij}<0$ denotes a mortality deficit compared to the province level. In cases where $O_{ij}=0$ for some $i$ and $j$, one can apply the alternative formula $SMR_{ij}=(O_{ij}+1/2)/E_{ij}$ to enable the computation of the logarithm of the standardised mortality ratios. This approach is described, for example, in Section 5 of Bivand et al. (2013). However, it is important to clarify that, although we know how the standardised mortality ratios were calculated, the raw data $O_{ij}$ and $E_{ij}$ were not provided by ULSS6 Euganea due to privacy considerations; only the $SMR_{ij}$ values are available for the analysis.

We also consider for the analysis two variables that are representative of the socioeconomic status of individuals in an area. Specifically, we consider:

• •

  Deprivation index, quantified through the Caranci’s index (Caranci and Costa, 2009). This index quantifies the deprivation, which represents socio-economic disadvantages, by considering various factors affecting residents in the territory. It is calculated as the standardised sum of five indicators, with higher scores indicating higher deprivation (Rosano et al., 2020). These five indicators include the level of education among individuals aged 15-60, the unemployment status, the condition of single-parent family with underage children, house renting and housing density per 100 m². The index typically ranges between $-5$ and 40, with most cases between $-$2 and 2. The data used to compute the index, sourced from the Italian 2011 census, are available at census section level and were aggregated to obtain data at municipality level.

• •

  Sale price of real estate per squared meter, which serves as an indicator of the general economic status of municipalities. These price data were scraped from www.immobiliare.it, a web portal that collects property sales listings in Italy, using the average of the prices in euro/m² on January 1${}^{\makebox{st}}$ of 2017, 2018, and 2019 (immobiliare.it, 2022). The specific time frame was chosen as it is characterised by stable price trends, without notable increases or decreases.

On male data, we fit perla considering $K=2$ and $K=3$, mainly motivated by the relatively small number of areas in the map. Results of the model selection are shown in Figure 10. The $\mathrm{DIC}_{3}$ criterion clearly supports $K=3$. In fact, the plot of the posterior distribution of $\boldsymbol{\mu}$ under $K=2$ (top panel of Supplementary Figure 7) displays multiple modes, which is a clear indication that a model with an additional component should be considered. We thus apply perla with $K=3$ and shrinkage factor (c,d), which is more parsimonious than (d,cd) while providing in practice the same value of $\mathrm{DIC}_{3}$. The final classification of the map is reported in the top-left plot of Figure 11. While cluster 1 extends over the whole territory, particularly on the northern area of the province, clusters 2 and 3 are prevalently located in the south-west. The posterior densities of $\{\exp(\mu_{kj})\}_{k=1,2,3;j=1,2,3}$ are displayed in the top-right panel of Figure 11, showing that cluster 1 is characterised by an excess mortality due to respiratory diseases ($\mathcal{M}(\mu_{12}|\mathrm{data})=1.216$, $\Pr\{\exp(\mu_{12})>1|\mathrm{data}\}=0.997$), and cluster 2 by a mortality excess due to circulatory diseases ($\mathcal{M}(\mu_{21}|\mathrm{data})=1.166$, $\Pr\{\exp(\mu_{21})>1|\mathrm{data}\}=0.966$) and by a mortality deficit due to respiratory diseases ($\mathcal{M}(\mu_{22}|\mathrm{data})=0.318$, $\Pr\{\exp(\mu_{22})>1|\mathrm{data}\}=0$). Lastly, cluster 3 does not show evidence of neither excess nor deficit of mortality due to the three causes of death considered. This is visible from the posterior distributions of the elements in $\{\exp(\mu_{3j})\}_{j=1,2,3}$, which are shrank toward one thanks to the cluster-specific shrinkage factor $\delta_{3}$. In Table 2, we report all the posterior SMR estimates and the probabilities of observing an excess mortality. Notice that, according to our model, deaths for tumours appear not to impact on the formation of spatial mortality clusters. This is highlighted by the disease-specific shrinkage factor $\zeta_{3}$, which collapses the posterior estimates of $\{\exp(\mu_{k3})\}_{k=1,2,3}$ toward one. We report also the trace plots of the cluster centroids $\boldsymbol{\mu}$ in the top panel of Supplementary Figure 8. Chains appear stable, demonstrating consistency across different parallel runs.