Bayesian disease mapping: Past, present, and future.

On the occasion of the Spatial Statistics' 10th Anniversary, I reflect on the past and present of Bayesian disease mapping and look into its future. I focus on some key developments of models, and on recent evolution of multivariate and adaptive Gaussian Markov random fields and their impact and importance in disease mapping. I reflect on Bayesian disease mapping as a subject of spatial statistics that has advanced to date, and continues to grow, in scope and complexity alongside increasing needs of analytic tools for contemporary health science research, such as spatial epidemiology, population and public health, and medicine. I illustrate (potential) utility and impact of some of the disease mapping models and methods for analysing and monitoring communicable disease such as the COVID-19 infection risks during an ongoing pandemic.

Introduction

On the occasion of the Spatial Statistics’ 10th Anniversary, I reflect on the past and present of Bayesian disease mapping, now part of an established branch of spatial statistics for (areal) lattice data.

Disease mapping has many connotations. For some, typically in the contemporary contexts of geography of disease and of population and public health, it concerns with small area mapping of disease incidence and mortality rates or relative risks, commonly for a population under study and as a part of research into geographical distribution of disease. For Haviland, 1871, Haviland, 1888, for example, mapping observed disease and cancer incidence and mortality rates for the counties of England served to display geographic distribution of a disease and to illustrate possible associations between disease and the environment; see Arnold-Forster, 2020, Arnold-Forster, 2021 for recent research into Haviland’s pioneering works on disease and cancer mapping.

Maps of rare diseases (such as cancer) based on crude rates or relative risks typically lack stability in the illustrated “geographic distribution of disease”. The early efforts made to statistically ascertain localities (such as counties) of elevated or lowered disease incidence and mortality rates and risks gave rise to a new topic in spatial statistics for lattice data, known as disease mapping, a topic of my main focus here. Crude rates of a rare disease calculated for (small) geographic areas such as counties and health administrative regions are typically subject to high chance variations. Earlier statistical research of disease mapping focused on the development of empirical Bayes (EB) models and related estimation procedures for stabilizing (i.e. smoothing) maps of disease rates or relative risks (Tsutakawa et al., 1985, Clayton and Kaldor, 1987, Tsutakawa, 1988, Manton et al., 1989, Cressie, 1993, to name a few). What was notably absent during the earlier years of EB disease mapping is the statistical inference of relative risks, inference in the sense that elevated or lowered disease risks are statistically ascertained.

The seminal works of Besag et al. (1991), Clayton and Bernardinelli (1992), Bernardinelli and Montomoli (1992), and Clayton et al. (1993) influenced and shaped a paradigm shift from the frequentist to Bayesian approach to disease mapping. Known as Bayesian disease mapping, the paradigm has Bayesian hierarchical models at its core and Markov chain Monte Carlo (MCMC) simulation methods provide the computational tools for Bayesian estimation, learning, and inference of all unknown quantities and parameters. Since the 1980s, the research agenda has broadened from providing tools for producing smoothed rate or risk maps to establishing a methodological paradigm with progressive, and systemtically developed, ideas, models, and methods for studying spatial and spatiotemporal variations, patterns, and clusters of disease and health outcome risks in populations of geographical regions, and for identification of (risk and protective) factors that might explain the variations and patterns and clusters. The body of Bayesian disease mapping models and related methods for estimation, learning, and inference has grown in remarkable proportion from the early development of Bayesian hierarchical models for univariate disease mapping to multidimensional disease mapping based on data of space, time, multivariate, and multi-array (e.g. data of space, time, multiple diseases, age group and/or gender stratified).

The richness of the disease mapping models and their wide range of health science applications are well documented in books (Lawson et al., 1999, Banerjee et al., 2004, Waller and Gotway, 2004, Cressie and Wikle, 2015, Lawson, 2018, Martinez-Beneito and Botella-Rocamora, 2019) and review articles (Richardson et al., 2004, Best et al., 2005, Wakefield, 2007, Ugarte et al., 2009, Lee, 2011, MacNab, 2011, MacNab, 2018), in addition to a large quantity of research papers. Here, I review key developments of Bayesian disease mapping, with a focus on recent evolution of multivariate and adaptive Gaussian Markov random fields and their impact and importance in Bayesian disease mapping. I illustrate (potential) utility and impact of some of the disease mapping models and methods for analysing communicable disease such as the COVID-19 infection risks in the context of pandemic intervention and monitoring. I reflect on Bayesian disease mapping as a subject of spatial statistics that has advanced to date, and continues to grow, in scope and complexity alongside increasing needs of analytic tools for contemporary health science research, such as spatial epidemiology, population and public health, and medicine.

Empirical Bayes disease mapping since the 1980s

While mapping disease and cancer rates or risks began century into the past, disease mapping gained attentions of the statistics community during the 1980s when cancer registry made population-level data more readily available for small geographical (e.g. administrative) areas of a region or country. A common goal of the researchers was to develop model-based procedures for stabilized cancer risk maps, from which environmental determinants of different types of cancer and their possible etiologic factors may be hypothesized for further investigations (Tsutakawa et al., 1985, Clayton and Kaldor, 1987, Tsutakawa, 1988, Manton et al., 1989 to name a few).

The 1980s also marked the rise of popularity of EB approaches to formulate and estimate models of multi-parameters, models in which the number and dimension of unknown parameters may equal or exceed those of the data (see Carlin and Louis, 2000a, Carlin and Louis, 2000b and references therein). Disease mapping models of the past and present are models of this nature.

An example of disease mapping model: Poisson model for count data

To facilitate further discussions, let us consider a typical disease mapping setting where we study areal-level cancer mortality risks through modelling of observed area-specific death counts, denoted , where is the number of geographic areas, such as wards, counties, local health administrative areas of a region or a country under study. For brevity, and without essential loss of generality, we will consider the simplest situation where the death counts are compared to a set of expected counts, denoted , which is often derived for a set of area- and age-specific at risk population and the corresponding age-specific reference rates, known as age-standardization (Clayton and Kaldor, 1987, Manton et al., 1989). The resulting relative risks, often defined as , are the standardized mortality ratios (SMR). This is a well known method used in studies of disease epidemiology and in cancer research (Breslow and Day, 1980). In disease mapping, this approach is often taken when areal-level age-specific data are not available or extremely small.

Working with rare events such as cancer mortality, the observed death counts may be viewed as Poisson events and modelled as follows: where , is the Poisson intensity, and is the unknown relative risk.

In frequentist disease mapping, the unknown relative risks in (1) are fixed parameters to be estimated from the data. One often assumes that the mortality counts are randomly varying (i.e. statistically independent) Poisson events over geographic areas. Under this assumption, the maximum likelihood estimates (MLEs) of the relative risks, denoted , are the previously mentioned SMRs: It is well known that the MLEs (2) may lack stability and efficiency when small counts are modelled. The independence assumption is often invalid due to possible spatial correlation, e.g., the tendency that areas of close proximity (say, neighbouring areas) have similar mortality risks.

In addition, Poisson distribution has the unusual characteristic that the intensity quantifies both the mean and the variance of the distribution. In real-life applications, this variance-equals-mean assumption may be invalid. For example, disease mapping data may exhibit extra-Poisson variation, i.e. , also known as Poisson over-dispersion (Breslow, 1984, Clayton and Kaldor, 1987, Besag et al., 1991; Breslow and Clayton, 1993).

Empirical Bayes disease mapping

Model (1) is a model of multi-parameters, where the dimension of parameters equals the dimension of the data. Taking an EB approach to this problem, one would place a prior on , denoted , where is a (vector of) hyper-parameter(s) to be estimated from data. This approach can facilitate borrowing information and risk smoothing and, with certain priors, account for Poisson over-dispersion. Many EB methods for disease mapping were proposed since the 1980s (Tsutakawa et al., 1985, Clayton and Kaldor, 1987, Tsutakawa, 1988, Manton et al., 1989, Cressie, 1993, Leroux et al., 1999 to name a few). Here, the basic ideas are outlined in the following two illustrative examples, both are taken from, and presented in detail in, Clayton and Kaldor (1987).

Since relative risks are non-negative, one prior option for is the gamma distribution , where and are the respective scale and shape parameters. This is a non-spatial prior, where the Poisson–gamma mixing leads to negative binomial distribution for , with expectation and variance functions: Negative binomial distribution is well known for its use to model Poisson over-dispersion (McCullagh and Nelder, 1983), as illustrated in expression (3) where the variance of is greater than its mean. Once the hyper-parameters and are estimated from the data, say, the MLEs of the negative binomial likelihood, denoted and , the EB estimate of is readily derived from its posterior expectation, : Replacing the and in by and leads to the EB estimate

In addition to modelling Poisson over-dispersion, another feature of this method is that the EB estimate is a compromise between the SMR and the estimated mean of , . This method leads to “global smoothing” of relative risks (see Clayton and Kaldor, 1987 for additional details).

However, global smoothing may be suboptimal when data are spatially distributed and correlated. An alternative approach is to formulate models that allow for spatial correlation and enable spatial smoothing. In Clayton and Kaldor (1987), a conditionally formulated Gaussian Markov random field (GMRF, Besag, 1974), named CAR for conditional autoregression, was proposed as a spatial prior for the log relative risks ensemble : where is the precision matrix, , is a scale parameter, is a spatial dependence parameter, and is the well known adjacency matrix for a given lattice (of map) and its neighbourhood system (e.g. two areas are neighbours if they share common border(s)).

This CAR model postulates symmetric spatial dependence, accounts for spatial correlation, facilitates spatial smoothing, and enables borrowing strength for risks prediction and inference. The key features of this conditionally formulated GMRF can be understood and explained via the conditional mean and precision functions of its full conditionals: where , , stands for the area and are neighbours. The area-specific conditional risk expectation of given , , is proportional to the total risks of its neighbouring areas, regulated by a positive spatial dependence parameter . The conditional risk prediction variance is regulated by the scale parameter .

Unlike the Poisson–gamma mixing, the Poisson-log-normal mixing does not lead to an analytically known distribution for . In Clayton and Kaldor (1987), an EM algorithm (Dempster et al., 1977) was used for a two-stage estimation of (and ) and its hyper-parameters.

The EB estimation is computationally straightforward for the first example and more complex and evolved for the second one. In an influential paper, Breslow and Clayton (1993) popularized several EB methods for inference in generalized linear mixed effects (GLMM) models, a model class that contains the majority of the disease mapping models developed to date.

The main limitation of the EB estimation of is that while the point estimates are nearly unbiased, the EB standard errors typically understate the associated estimation uncertainties. This can be readily observed from the “plug-in” EB estimation of the posterior (log) relative risks in the first example, where the uncertainties associated with the estimation of the hyper-parameters may not be adequately accounted for. This is the well known “naive” EB estimation issue for nearly all EB approaches to estimation of (Carlin and Gelfand, 1990, Carlin and Gelfand, 1991, Carlin and Louis, 2000a, MacNab et al., 2004, Ainsworth and Dean, 2006, MacNab and Lin, 2009).

Into the 1990s and 2000s, some progresses were made to improve EB posterior estimation of relative risks, mainly to provide better standard error estimates by bootstrap methods (Carlin and Gelfand, 1990, MacNab et al., 2004) or via analytic adjustment (Ainsworth and Dean, 2006, MacNab and Lin, 2009). With adequately calibrated standard errors that quantify risk uncertainties, EB methods for posterior inference of relative risks may be used to statistically ascertain elevated or lowered disease risks.

BayesIan disease mapping from the 1990s: A paradigm shift

Over the broad landscape of statistical science, the 1990s saw the rise of Bayesian statistics and a new trend of Bayesian hierarchical approaches to formulating complex models, with hierarchical Bayesian methods for posterior estimation and inference of all unknown quantifies and parameters, typically implemented by Markov chain Monte Carol (MCMC) simulations (Gelfand and Smith, 1990, Smith and Gelfand, 1992, Gelman et al., 1995, Gilks et al., 1998). Amid the MCMC movement into the 1990s and the new 21st century, a shift from EB to fully Bayesian (FB) approach to disease mapping was taking place (Bernardinelli and Montomoli, 1992, Gilks et al., 1998, Lawson et al., 1999, Gelman et al., 1995, Congdon, 2001).

The move from EB to FB disease mapping was an analytically natural transition. The two are analytically comparable in the sense that the latter produces Bayesian posterior estimates for the prior parameters, but doing so by placing non-informative or weakly informative priors on the prior parameters, an approach facilitates “data-driven” posterior estimation and inference of the prior parameters (Carlin and Louis, 2000a, MacNab et al., 2004).

Analytic and conceptual transitions take place when the fully Bayesian approach constructs models hierarchically (Cressie and Wikle, 2015) and within a Bayesian hierarchical inferential framework. Within this framework, the prior distribution, such as the previously mentioned gamma or CAR prior, can be postulated as an underlying process model for the underlying risks that give rise to the observed disease incidence and mortality counts. Statistical inference utilizes the Baye’s rule to compute (approximate) posterior distributions for all unknown quantities and parameters. Conceptually, one can explore meaningful and intuitively appealing prior and hyper-prior options within a Bayesian hierarchical inferential framework for model formulation and Bayesian posterior estimation, learning, and inference of relative risks and parameters of interests and importance. This approach provides a cohesive and principled way to quantify risk uncertainties and uncertainties of all unknown parameters.

With the arrival of freely available BUGS/WinBUGS software (Spiegelhalter et al., 2003) for Bayesian hierarchical analysis and later the GeoBUGS add-on module (Thomas et al., 2004) to WinBUGS for Bayesian analysis of spatial data, in addition to various publicly available R packages, e.g. the R-INLA (Rue et al., 2017) and the CARBayes (Lee, 2013), for (approximate) Bayesian analysis of spatial areal data, fully Bayesian approach to disease mapping model formulation, estimation, and inference became mainstreamed and stayed so to this day.

Multidimensional disease mapping

Multidimensional disease mapping is a broad topic that concerns with spatiotemporal disease mapping (Held and Besag, 1998, Held, 2000, MacNab and Gustafson, 2007, Ugarte et al., 2017), multivariate disease mapping (Kim et al., 2001, Jin et al., 2007, Greco and Trivisano, 2009, Martinez-Beneito, 2013, MacNab, 2018), and multi-array disease mapping, say, when space, time, age (and gender, etc.) stratifications are all considered in a disease mapping model (Held and Besag, 1998, MacNab, 2003a, Martinez-Beneito et al., 2017, Martinez-Beneito and Botella-Rocamora, 2019). Literature on this topic is quite rich (see Lawson, 2018, Martinez-Beneito and Botella-Rocamora, 2019 and references therein). Here, I briefly reflect on common ideas and models that characterize cross-dimensional dependencies and interactions.

Without essential loss of generality, let us begin with adding another dimension to the Poisson data model (1), in which one has a matrix of log relative risks , where for spatiotemporal disease mapping or for mapping of diseases.

Spatiotemporal disease mapping models presented in the literature focused mainly on characterizing the presence of space and time risk interactions, in addition to model components of risks in space and time, respectively (Waller et al., 1997, Held and Besag, 1998, Held, 2000). These models are often motivated to facilitate borrowing information and risk smoothing within and between space and time. They are also natural extensions of the spatial risk models.

Some of the models are richly parameterized to hypothesize space and time risk interactions, in addition to decompose the into (parameterized) components of space and time. Typical examples are the models of Waller et al. (1997) and Held (2000), two highly cited and influential papers. For example, Held (2000) proposed a general model formulation and considered types of space–time interactions and associated prior options for . These spatiotemporal (interaction) models may be highly parameterized and, in this situation, additional identification constraints may be imposed (Held, 2000).

Alternatively, one can take a parsimonious way to parameterize spatiotemporal models. For example, the following model formulation was considered in the literature (MacNab and Dean, 2001): in which functional analysis may be used to hierarchically model space time interactions. For instance, illustrated in MacNab and Dean (2001), MacNab and Gustafson (2007), and MacNab (2007), space–time interactions can be modelled via spatially varying splines (e.g. B-splines, P-splines, and smoothing splines), in addition to modelling the time components by a fixed spline for non-linear risk trend. In Ugarte et al. (2010), Lee and Durbán (2011), Ugarte et al. (2017), and Feng (2021), space–time interactions are modelled by tensor-based 2d or 3d P-splines or thin plate spline. Lawson (2018) and MacNab (2018) also considered parsimonious spatiotemporal model parameterizations. They illustrate modelling linear risk trends and , which can be viewed as a special (degenerated) case of the previously mentioned spline-based approaches for modelling non-linear risk trends.

Modelling spatially varying non-linear risk trends via spatially varying random walk (RW) or auto-regressive (AR) or CAR processes are also viable options (Held, 2000, Ugarte et al., 2017, MacNab, 2018, Sahu and Böhning, 2021); see Section 9 for an example. Spatiotemporal disease mapping models were mainly developed for posterior prediction of the unobservable relative risks, although a few have been developed for short-term forecasts of disease rates (Etxeberria et al., 2015, Sahu and Böhning, 2021).

Without sufficient difficulty, the spatiotemporal models (7), (8) can be extended further for multi-array disease mapping. Adding age-effects is one example. Age-effects could be included in (7) or (8) as additional components (Held and Besag, 1998, Goicoa et al., 2016, Martinez-Beneito et al., 2017) or modelled by a spline for non-linear age-effects (MacNab, 2003a). Multivariate spatiotemporal disease mapping is another direction (Tzala and Best, 2008, Jack et al., 2019, Vicente et al., 2020, Vicente et al., 2021).

Multivariate disease mapping has been an active area of research since the early 2000s. The main motivation of the efforts is to develop methods for jointly modelling diseases that share common risk factors, enabling more effective borrowing-information, and facilitating spatial and cross-spatial risk smoothing. Two main approaches to multivariate disease mapping are considered in the disease mapping literature. One is to model as a sum of two components: where is a component shared by all diseases. This is the shared component model (SCM) proposed by Held and Best (2001). SCMs typically hypothesize positive cross-variable local dependencies and correlations (MacNab, 2010). The other more flexible approach is to place a multivariate GMRF prior on vec() or vec(), which I discuss in Section 7.

Disease mapping with areal-level covariates: Ecological regression

In the context of disease mapping, including a regression part in a disease mapping model can serve several important and practical purposes. One important purpose of a spatial regression model is to quantify disease risk variations explained and unexplained by the included risk (and/or protective) factors (Clayton et al., 1993, MacNab, 2003b, MacNab, 2003c, MacNab, 2004).

Known as ecological modelling and analysis, ecological regression can serve to identify areal-level characteristics (e.g. environmental exposure, socio-economic deprivation, contextual characters) that influence the disease risks. Most of the ecological models developed in the disease mapping literature were aimed for this purpose. They are commonly motivated and applied to research into the (social and/or physical) environmental determinants of health (Clayton et al., 1993, MacNab, 2004, Lee and Sarran, 2015, Cameletti et al., 2019, Huang et al., 2021) or population-level health inequality (Lee and Lawson, 2016, Jack et al., 2019).

Ecological regression faces several methodological challenges. One of them is the long-standing ecological bias, referring the situation when the degree of association between risk exposure and disease at the areal-level is different from that of the individual level, due to within-area variability in exposures and confounders (Clayton et al., 1993, Wakefield and Salway, 2001, Wakefield, 2007, Lawson, 2018).

Another analytic challenge is spatial confounding, commonly discussed as collinearity between fixed effects and random effects in a spatial GLMM, and in the presence of unmeasured confounding. Literature of ecological studies offered two main approaches to the issue. One is to deconfound the fixed and random effects (Reich et al., 2006, Hughes and Haran, 2013, Ugarte et al., 2017, MacNab, 2018, Zimmerman and Hoef, 2021). Another is to formulate transformed GMRF or MGMRF (Prates et al., 2015, Hughes, 2015, Prates et al., 2021) so that the fixed and random effects can be modelled independently. Named TGMRF or TMGMRF, the basic idea of this approach is to model the relative risks , and its spatiotemporal or multivariate extensions, via a copula-inspired model structure that allows fixed and random components to be estimated independently: The fixed effects are modelled on the univariate margins, whereas the random effects spatial dependencies are modelled via a GMRF or MGMRF spatial dependence matrix. One limitation of this approach is that the regression part in the model is designed to explain marginal variations, not spatial dependence. In traditional spatial regression GLMMs, the included covariates may (partially) explain risk variation or spatial risk dependence or both; see MacNab, 2003b, MacNab, 2003c, MacNab, 2004, MacNab, 2009, MacNab, 2010, and Section 9, for illustrative examples.

Of the deconfounding approach, earlier proposals are known as the restricted spatial regression (RSR) method (Reich et al., 2006, Hughes and Haran, 2013, Ugarte et al., 2017, Zimmerman and Hoef, 2021). The key idea of the RSR approach is to re-parameterize a -dimensional spatial prior and project the original Gaussian spatial prior onto a Gaussian spatial prior of its sub-space orthogonal to the column space of the fixed effects model matrix (Zimmerman and Hoef, 2021). Very recently, Zimmerman and Hoef (2021) researched into the RSR approach, named it deconfounding spatial confounding, and put the method into question. In the context of spatial linear mixed model and its frequentist inference, they showed that the RSR approach can lead to underestimation of the uncertainties associated with the estimation of the fixed effects and is generally inferior to that of the original spatial regression model.

An alternative to the RSR method is to formulate spatially misaligned ecological model, misaligned in a sense that the spatial prior of the random effects and the included covariates are varying at different spatial scales and/or configurations (Lee and Sarran, 2015, Prates et al., 2019, Azevedo et al., 2021). This approach to deconfounding is also consistent with the Paciorek (2010) work on spatial confounding, which offered insight into the importance of scale for spatial-confounding bias and precision of spatial regression estimators. In Paciorek (2010), analytic and simulation results suggested evidence that bias can be reduced when fitting a spatial model in which the variation of the covariate is at a scale smaller than the scale of the unmeasured confounding.

Errors-in-covariate is another important issue in ecological regression, where areal-level covariates are commonly measured imprecisely, or because surrogates or proxies are only available and used. In the Bayesian disease mapping literature, the main approach is to formulate measurement models for the covariates and incorporate them into a hierarchically formulated Bayesian spatial regression model (Bernardinelli et al., 1998, Wakefeild and Morrise, 1999, MacNab, 2009, MacNab, 2010, Lawson, 2018).

Conditionally formulated GMRFs

Conditionally formulated GMRFs have played important roles in disease mapping and ecological regression (Clayton and Kaldor, 1987, Besag et al., 1991). They are formulated to model spatial dependencies and spatially structured correlation and covariance functions (MacNab, 2011, MacNab, 2016a, MacNab, 2016b, MacNab, 2018). Formulating GMRF via full conditionals is conceptually appealing for articulating the idea of spatial smoothing, as illustrated earlier for CAR (5) expressed by its full conditionals (the Expression (6)). Analytically, the GMRF full conditionals also facilitate coding the powerful Gibbs sampler as a computational tool for posterior estimation of the GMRFs via MCMC simulations; see [10] for an accessible explanation of Gibbs sampler.

Several CARs have been proposed since the CAR (6). Three stood the time and are commonly considered prior options for disease mapping. They are the well known intrinsic CAR (iCAR, Besag et al., 1991), proper CAR (pCAR, Cressie, 1993, Sun et al., 1999), and Leroux et al. CAR (LCAR, Leroux et al., 1999); see Table 1 for their full conditionals and associated precision matrices.

Table 1

Three CARs commonly used in Bayesian disease mapping: when the th and th areas are neighbours (denoted ) or otherwise; . : The covariance matrix is presented; , : The covariance matrix is presented; ; : and are (assumed) statistically independent. For all models in Table 1: .

Model	E(ψ_i|ψ _−i)	Var(ψ_i|ψ _−i)	Ω
iCAR(σ^2)	∑_j∼iψ_jw_i+	σ^2w_i+	σ^−2(D _w−W )
(Besag et al., 1991)
pCAR(λ,σ^2)	λ∑_j∼iψ_jw_i+	σ^2w_i+	σ^−2(D _w−λW )
(Cressie, 1993)
LCAR(λ,σ^2)	λ∑_j∼iψ_j1−λ+λw_i+	σ^21−λ+λw_i+	σ^−2(λ(D _w−W )+(1−λ)I _n)
(Leroux et al., 1999)
BYM(σ_s,σ_h)^†	ψ =ψ ^s+ψ ^h,ψ ^s⊥ψ ^h		σs2(D _w−W )^−1+σh2I _n
(Besag et al., 1991)
MBYM(c,σ)^†	ψ =λψ ^s+1−λψ ^h,ψ ^s⊥ψ ^h		σ^2(λ(D _w−W )^−1+(1−λ)I _n)
(MacNab, 2011)

The pCAR() and LCAR() are two-parameter and full rank GMRFs, where and are the respective spatial and non-spatial parameters: is a spatial dependence or weight or smoothing parameter and is a (non-spatial) Gaussian scale parameter (Cressie, 1993, Leroux et al., 1999). Typically defined on an irregular lattice of areal map, the pCAR, with its conditional expectation and precision functions is a slight modification to the CAR (6), where in (10) is a set of scaling factors. Consequently, when , the pCAR conditional expectation is proportional to an average of neighbourhood risks, rather than the sum of the neighbourhood risks as is defined in CAR (6). And the pCAR (10) conditional precision is proportional to , which implies that the greater the number of neighbours, the higher the conditional precision. That is, borrowing information from a higher number of neighbours could lead to higher precision of posterior risk prediction. This could be a plausible assumption in the context of disease mapping, where an areal map under study typically constitutes an irregular lattice with varying sizes of neighbourhoods over the map.

Another feature of the pCAR, which is not often discussed in the literature, is that it postulates asymmetric conditional dependency (i.e. direct influence) of on , quantified by versus on , defined by , provided , , and . That is, the direct influence of the area on its neighbouring area is inversely proportional to the neighbourhood size of the area . Put it another way, an area with a higher number of neighbours is less influenced by its neighbour who has a lower number of neighbours. This could also be an intuitively plausible assumption, consistent with the conditional precision function in (10); one might expect that an area of higher precision of risk prediction should be less influenced by an area with a lower precision of (predicted) risk. The spatial dependence parameter in pCAR is often called a spatial smoothing parameter; it regulates smoothly varying relative risks over the map.

To facilitate discussion, I name the coefficients of the conditional autoregression in (10) the coefficient of influence functions, or simply influence functions, denoted

The LCAR conditionals also allow asymmetric spatial dependencies, provided , where That is, the LCAR influence function is a non-linear but increasing function of . The spatial dependence parameter has triple roles and is also a spatial weight or smoothing parameter; it weights a precision matrix of the iCAR() and a precision matrix of IIDN() (normal) variates, a mixing of purely local smoothing and global smoothing (Congdon, 2008). The LCAR conditional variances are non-linear decreasing functions of (see Table 1); that is, also partially controls the risk prediction variance of conditional on (MacNab, 2011). Noted in MacNab (2018), the LCAR parameterization is an “entangled” spatial and non-spatial parameterization, and, as a result, LCAR has limited options for multivariate generalization (see Section 7).

The iCAR() is a GMRF with a singular precision matrix of rank . It is the limiting distribution of the pCAR() and LCAR () (when ). The iCAR() is typically considered as a spatial prior for modelling spatially structured risk variation or clustered heterogeneity. It is commonly used when the -vector of log relative risks is modelled as additive components of , for modelling spatially structured heterogeneity or effects of omitted covariates that are spatially varying, and , for modelling extra-Poisson variation or effects of omitted covariates that are randomly varying. This is the well known Besag, Yorke, and Mollie (BYM) model (Besag et al., 1991), widely used in Bayesian disease mapping and ecological regression. This model is highly parameterized. The BYM is typically used as a prior in Bayesian disease mapping, where estimation and inference of the relative risks may be implemented with (weakly) informative priors for (functions of) the scale parameters and/or (Thomas et al., 2004, MacNab, 2011).

To improve computational efficiency and parameter identification, MacNab (2011) proposed a modification to the BYM, named MBYM: , where is a weight parameter, weighting between the covariance matrix of the iCAR and the covariance matrix of the IIDN; see Table 1 for the BYM and MBYM covariance matrices.

The LCAR has gained popularity due in part to the analytic result that when , the LCAR reduces to independent and identical Gaussian priors for the log relative risks, whereas the pCAR, when , reduces to independent Gaussian priors with prior precisions proportional to (MacNab, 2003c, MacNab, 2011, Lee, 2011).

However, compared to LCAR, pCAR has its own advantages. Its spatial and non-spatial parameters play separate and different roles, one regulates spatial dependency, the other controls non-spatial variance. With its rich options of multivariate generalization, the pCAR, and some of its adaptive and multivariate extensions, have theoretical and practical appeals for modelling and interpreting spatial dependencies (MacNab, 2018, MacNab, 2020, MacNab, 2021a, MacNab, 2021b see Sections 7, 8).

Of the five CAR models, none has shown to outperform the others in all disease mapping situations. This is consistent with the analytic results that the iCAR, pCAR, and LCAR have different spatial dependence and correction functions (see details in MacNab, 2011, MacNab, 2014) and the five CAR models may play nuanced roles in Bayesian disease mapping applications. Bayesian sensitivity analysis under various prior and hyper-prior options for Bayesian estimation, learning, and inference is an important part of Bayesian disease mapping application (MacNab, 2011); and, with goodness-of-fit and complexity assessment, say, based on deviance informion criterion (the devance, pD, DIC scores, [125], also see Section 9), it remains a viable and commonly taken approach to model evaluation, comparison, and selection.

Multivariate GMRFs for multidimensional disease mapping

In the broad context of multidimensional disease mapping discussed earlier, multivariate GMRFs play important roles in modelling multidimensional spatial dependencies and cross-covariances, enabling borrow-information over multi-dimensions, and facilitating spatial and cross-spatial smoothing. The five commonly used disease mapping models in Table 1 have their multivariate extensions presented in the literature. A list of illustrative models are presented in Table 2, where their joint covariance matrices are given.

Table 2

A selected multivariate CARs. is a matrix, for -variate CARs.

Model	ψ	Σ _M
MiCAR(Σ )	vec(ψ ^⊤)	(D _w−W )^−1⊗Σ
(G+V 2003)1
MpCAR(c,Σ )	vec(ψ ^⊤)	(D _w−cW )^−1⊗Σ
(G+V 2003)1
MLCAR(c,Σ )	vec(ψ ^⊤)	(c(D _w−W )+(1−c)I _n)^−1⊗Σ
(M+G 2007)2
MLCAR(c ,Σ )	vec(ψ ^⊤)	(I _n⊗A )Bdiag[(c_k(D _w−W )+(1−c_k)I _n)^−1](I _n⊗A )
(MacNab, 2018)
MpCAR(C _s,Σ )	vec(ψ )	(A ⊗I _n)S (C _s)^−1(A ^⊤⊗I _n)
(Jin et al., 2007)
MpCAR(C ,Σ )	vec(ψ )	(A ⊗I _n)S (C )^−1(A ^⊤⊗I _n)
(G+T 2009)3
MpCAR(C ,ρ ,σ )	vec(ψ )	(σ ^2⊗I _n)S (ρ ,C )^−1(σ ^2⊗I _n)
(Sain et al., 2010)
MultBYM(Σ _s,Σ _h)^†	vec(ψ )=vec(ψ ^s)+vec(ψ ^h)	(D _w−W )^−1⊗Σ _s+I _n⊗Σ _h
(Thomas et al., 2004)

1: Gelfand and Vounatsou (2003). 2: MacNab and Gustafson (2007). 3: Greco and Trivisano (2009). is a covariance matrix. . , , , where , , is a symmetric matrix or diagonal matrix, , is the upper triangular part of ; . : ; and are (assumed) statistically independent; : ; and , , are (assumed) statistically independent.

For mapping two or more diseases, for example, conditionally formulated multivariate GMRFs, such as the 2-fold CAR of Kim et al. (2001), the proper multivariate CARs of Gelfand and Vounatsou (2003), and the coregionalized MpCARs of Jin et al. (2007), were proposed for modelling multivariate disease risks that are correlated at co-locations as well as within- and cross-neighbourhood (i.e. cross-correlation between diseases at neighbouring locations). It is worth mentioning that the three groups of authors took three approaches to arrive their MCARs, all are multivariate generalizations of the pCAR. Working with a matrix of , one has a variety of ways to construct MGMRFs; see [89] for a discussion and references therein.

Since the influential works of Besag (1974) and Mardia (1988), MGMRFs are often formulated via full conditionals (e.g. Kim et al., 2001, Gelfand and Vounatsou, 2003, Sain and Cressie, 2007, Sain et al., 2011). However, as shown in MacNab (2018), formulation and implementation of MGMRFs face several challenges, most notably the entanglement of spatial and non-spatial parameterizations, compliance with symmetry requirement, and enforcement for positivity. With the exception of separable MCARs, neither the Besag (1974) nor the Mardia (1988) framework readily enables multivariate generalizations of all CARs.

One example is the LCAR, for which non-separable multivariate generalization of full conditionals is not readily derived, due in part to its entangled spatial and non-spatial parameterizations. The MLCAR() presented in Table 2 was derived via linear transformation of independent LCARs (MacNab, 2016a), a method known as coregionalization (Jin et al., 2007).

Linear transformation is a well known and powerful idea for multivariate generalizations of univariate distributions and models. In a series of recent papers (MacNab, 2016a, MacNab, 2016b, MacNab, 2018, MacNab, 2021a, MacNab, 2021b), and built on the works of Gelfand et al., 2004, Jin et al., 2007, Greco and Trivisano, 2009, MacNab develops a coregionalization framework that connects and unifies the various lines of MGMRF development. The framework enables a systematic development, categorization, and implementation of a broad range of multivariate models via linear or spatially varying coregionalization, with extensions to locally adaptive models. The coregionalization methods readily facilitate multivariate generalizations of (any) GMRF and contain the major MCARs proposed in the literature.

To briefly explain the key ideas, consider a -variate latent MGMRF construction with the following precision matrix (Greco and Trivisano, 2009, MacNab, 2018) provided , where is a full matrix of spatial dependence parameters, , are area-specific scaling factors, is the upper triangular part of , is the “neighbourhood” connectivity or weight matrix. is a spatial dependence matrix for locally dependent MGMRF constructions. Expression (11) is simplified to when is a symmetric matrix or a diagonal matrix of spatial dependence parameters. It also reduces to a separable model when .

Without the “entanglement” challenges, MGMRFs of construction (11) can be readily formulated within the Besag (1974) or Mardia (1988) framework (MacNab, 2018).

Linear transformation of , denoted leads to a coregionalization MGMRF (cMGMRF) with covariance matrix (Jin et al., 2007, Greco and Trivisano, 2009, MacNab, 2016a, MacNab, 2016b, MacNab, 2018): where is a coregionalization coefficients matrix, , is a covariance matrix.

Notice that (12), (13) exemplify one way to decomponse and mix spatial and non-spatial components that characterizes and parameterizes cross-spatial interactions. Alternatively, let the latent fields be , a spatially varying (non-stationary) coregionalization (SVC) construction may be formulated via to have the following covariance matrix (MacNab, 2018): where . In SVC (15) the latent variables are only correlated at same locations. SVC constructions with latent (M)GMRF can also be readily formulated.

The decomposition and mixing schemes of (12), (14) are typically used to produce valid covariance and cross-covariance functions for cMGMRFs. The precision matrices of the cMGMRFs loss their appeal of interpreting cross-spatial dependencies (MacNab, 2018, MacNab, 2021a, MacNab, 2021b).

As an alternative to (12), (14), a cMGMRF can be formulated via to have the following precision matrix: where , is a by non-spatial (with-in location) partial correlation matrix, is a diagonal matrix of scale parameters (MacNab, 2018, MacNab, 2020, MacNab, 2021a). Notice that is a precision (dependence) matrix of a latent MGMRF; its full conditionals are readily available (Sain et al., 2011, MacNab, 2018, MacNab, 2020, MacNab, 2021a). It models spatial dependencies by and non-spatial dependencies by . A main feature of this approach to separating and mixing the cMGMRF spatial and non-spatial components, with its spatial and non-spatial parameterizations, is that the cMGMRF preserves the appeal of interpreting cross-spatial dependencies in the latent MGMRF with precision matrix ; also see MacNab, 2018, MacNab, 2021a for connections between the cMGMRFs (13), (17).

With options of parameterization for , , , , the three constructions contain all major proposals of MGMRFs to date and much more; they are rich options of multivariate generalizations of pCAR. These constructions may be parameterized to build multivariate spatial and spatiotemporal smoothers, as well as multivariate spatial and spatiotemporal covariance models. cMGMRFs of constructions (13), (15) could be parameterized for latent spatial components and factor analysis, include, but not limited to, principal spatial components analysis and dimensionality reduction (see MacNab, 2018, MacNab, 2021a for illustrative examples).

Of the cMGMRF (17) construction, options of parameterizations for the spatial and non-spatial components may be considered for their intuitively plausible characterizations and interpretations of conditional spatial dependencies within and between variables, as well as local interactions in space and time (Sain et al., 2011, MacNab, 2018, MacNab, 2020, MacNab, 2021a, MacNab, 2021b, Prates et al., 2021).

A -variate GMRF that allows for asymmetric cross-covariance functions can be formulated by linear coregionalization of latent MGMRF with a full matrix of spatial and cross-spatial dependence parameters or by spatially varying coregionalization of independence GMRFs; the latter is a non-stationary adaptive MGMRF with adaptive non-spatial parameterization (Martinez-Beneito, 2020, MacNab, 2021a). In addition, as I show in Section 8, non-stationary GMRFs that allow for asymmetric pair-wise spatial dependencies can be formulated via adaptive parameterization.

Adaptive CARs

GMRFs are undirected graphical models, commonly defined for a lattice system of nodes and edges (Rue and Held, 2005). Of the GMRFs commonly used in disease mapping (say, those in Table 1), the edges are defined by the matrix of a given map and its neighbourhood definition. These GMRFs often have single spatial and/or scale parameter. In the context of disease mapping, they postulate pair-wise interactions of neighbouring risks. They typically lead to smoothly varying posterior relative risks prediction and lack the flexibility to allow within- and between-neighbourhood risk heterogeneities that likely happen in parts of the map (MacNab et al., 2006, Brewer and Nolan, 2007).

To mitigate the mentioned limitations, adaptive CARs have been proposed with applications to analysis of imaging data (Brezger et al., 2007), analysis of correlated data (Reich and Hodges, 2008), and disease mapping (Lu and Carlin, 2005, MacNab et al., 2006, MacNab, 2018, Lu et al., 2007, Brewer and Nolan, 2007, Congdon, 2008, Lee and Mitchell, 2012, Rushworth et al., 2017, Gao and Bradley, 2019, Corpas-Burgos and Martinez-Beneito, 2020, among others).

Adaptive CARs proposed to date can be broadly grouped under two motivating directions. One is to hierarchically formulate adaptive CAR(GMRF) with unknown adjacency matrix to be modelled and estimated from data. This is the direction taking by Lu and Carlin (2005), Lu et al. (2007), Lee and Mitchell (2012), Rushworth et al. (2017), and Gao and Bradley (2019), among others. With the exception of Gao and Bradley (2019), these authors proposed adaptive iCAR, CAR (5), and LCAR in which the non-zero elements of were modelled as Bernoulli random variates or random variates of unit interval (i.e. , Rushworth et al., 2017). In Gao and Bradley (2019), an adaptive iCAR was formulated for unknown adjacency matrix , where , were modelled as Bernoulli random variates. This line of research primarily concerned with adjacency modelling for boundary analysis and detection (Lu et al., 2007, Lee and Mitchell, 2012, Gao and Bradley, 2019) and/or analysis of spatial discontinuity in disease risk maps (Lee et al., 2014, Rushworth et al., 2017, Sahu and Böhning, 2021).

An adaptive LCAR was also proposed in Baptista et al. (2016), in which the matrix is replaced by a similarity matrix , where the elements , are functions of covariates that characterize risk similarities over a map.

Formulating adaptively parameterized CARs over a given is the second approach. Under this approach, adaptive CARs have been proposed for each of the CARs presented in Table 1; the key proposals are presented in Table 3. This group of adaptive CARs may be used to model locally varying (asymmetric) spatial dependencies, interactions, and heterogeneities. These adaptive CARs might capture micro-level phenomenons such as locally varying (asymmetric) risk dependencies and influences and heterogeneities that typically manifest macro-level phenomenons such as the spatial risk heterogeneities and discontinuity. They can be broadly sub-grouped into three categories. One contains CARs with adaptive spatial (dependence) parameters; they were typically motivated to model locally varying spatial dependencies (MacNab et al., 2006, Congdon, 2008, MacNab, 2018, MacNab, 2021a). Another includes CARs of adaptive scale or precision parameters that serve as spatial weights for modelling weighted pair-wise risk interactions and locally varying risk heterogeneities (Corpas-Burgos and Martinez-Beneito, 2020). CARs with adaptive spatial and scale parameters are among the third category. Here, I highlight the key features of the three categories of adaptive CARs.

Table 3

Key options of adaptively parameterized CAR in the literature.

Model	Influence(j,i)	Var(ψ_i|ψ _−i)	Σ
pCAR(a)	c_i/w_i+	σ^2/w_i+	–
pCAR(b) (MacNab, 2018)	λc_ic_j/w_i+	σ^2/w_i+	σ^2(D _w−λW ^c)^−1
LCAR(a) (MacNab et al., 2006)	1/(ci−1+(w_i+−1))	σ^2/w_i+(c_i)	–
LCAR(b) (Congdon, 2008)	c_j/(ci−1+(w_i+−1))	σ^2/w_i+(c_i)	σ^2(c (D _w−W )c +(I _n−c ))^−1
LCAR(c) (MacNab, 2018)	c_j/(ci−1/2+(w_i+−1))	σ^2/w_i+(c_i)	σ^2(c ^1/2(D _w−W )c ^1/2+(I _n−c ))^−1
iCAR(a) (C-B 2020)2	c_j/∑_j∼ic_j	σ^2/wi+c	σ^2(D wc−W ^c)^−1
pCAR(c)	λc_j/∑_j∼ic_j	σ^2/wi+(c)	σ^2(D wc−λW ^c)^−1
LCAR(d) (C-B 2020)2	λc_ic_j/wi+(c,λ)	σ^2/wi+(c,λ)	σ^2(λ(D wc−W ^c)+(1−λ)I _n)^−1
LCAR(e) (C-B 2020)2	λc_j/wi+(c~,λ)	σ^2/(c_iwi+(c~,λ))	σ^2(λ(D wc−W ^c)+(1−λ)c )^−1
iCAR(b)^† (MacNab, 2018)	c_j/(c_iw_i+)	(σs_i)^2/w_i+	σ ~(D _w−W )^−1σ ~
pCAR(d)^† (MacNab, 2018)	λc_j/(c_iw_i+)	(σs_i)^2/w_i+	σ ~(D _w−λW )^−1σ ~
LCAR(f) ^† (MacNab, 2018)	λc_j/(c_iwi+λ)	(σs_i)^2/wi+λ	σ ~(λ(D _w−W )+(1−λ)I _n)^−1σ ~
pCAR(e) MacNab (2018)	λc_ic_j(σ_iσ_j)^−1/w_i+	σi2/w_i+	σ (D _w−λW ^c)^−1σ
Brezger et al. (2007) iCAR	w~_ij/w~_i+	σ^2/w~_i+	σ^2(D _w~−W ~)^−1,w~_ij=c_ijw_ij
Brewer (2007)1 iCAR	τ_ij/wi+τ	1/wi+τ	(D wτ−W ^τ)^−1

Brewer 2007: Brewer and Nolan (2007) and Reich and Hodges (2008); C-B 20202: Corpas-Burgos and Martinez-Beneito (2020), is replaced by . ; when or otherwise; ; ; ; , , ; in Brewer and Nolan (2007), in Reich and Hodges (2008). , ; , ; , . : , . E.

Adaptive CARs of the category I

There are important differences between the adaptive pCARs and LCARs of this category. A common feature of the adaptive pCARs is that the adaptive parameters control spatial dependencies in the conditional mean functions, whereas its scale parameter regulates (log relative) risk prediction variances in the conditional variance functions (see Table 3). The notable differences between the pCARs are their influence functions, as illustrated here for the pCAR(a) and pCAR(b): In pCAR(a), quantify how each area is directly influenced by its neighbours: The higher the , the greater the prediction is directly influenced by its neighbouring . When is near-zero, prediction of is not directly influenced by its neighbouring .

The pCAR(b) is a variant of the adaptive pCAR proposed in MacNab (2018). Here, the pCAR(b) has an additional spatial dependence parameter that characterizes neighbourhood risk dependence under the new weight matrix (see Table 3). Similar to its non-adaptive counterpart, its coefficient of influence functions imply asymmetric pair-wise local dependencies and influences, provided : An area with a higher number of neighbours is less influenced by its neighbour who has a lower number of neighbours. Two neighbours exert same influence on each other if and only if they have the same number of neighbours.

Similar to the LCAR(a), the adaptive pCAR(a) full conditionals do not lead to a GMRF due to noncompliance with the GMRF symmetry requirement. Nevertheless, the adaptive pCAR(a) or LCAR(a) can be used as a prior in a Bayesian hierarchical model for characterizing locally varying spatial dependencies and heterogeneities (which I illustrate in Section 9).

In the LCARs, the adaptive spatial parameters serve dual roles of regulating locally varying spatial dependencies and heterogeneities (see Table 3). For example, the , is a non-linear increasing function of , the corresponding conditional variance function is a non-linear decreasing function of .

Illustrated in Table 3, the LCAR(b) and LCAR(c) are comparable GMRFs. While the LCAR(c) reduces to its non-adaptive LCAR, the LCAR(b) does not. Compared to the LCAR(a), the two have more complex influence functions for interpreting locally varying spatial dependencies and influences (given in Table 3). An alternative interpretation of the two is that, as adaptive GMRFs, their adaptive parameters are mixing parameters that allows adaptive mixing of spatial and non-spatial smoothing (Congdon, 2008).

Adaptive CARs of the category II

Adaptive CARs of this category began with proposals of adaptive iCAR, primarily motivated for modelling spatial discontinuity and heterogeneity. The basic idea is to replace the connectivity matrix in iCAR by adaptive parameterization of a weight matrix of non-negative elements, denoted in Table 3, where , models locally weighted pair-wise interaction between and . In Brezger et al. (2007), for example, the elements were estimated via a gamma prior. Brewer and Nolan (2007), Reich and Hodges (2008), and Corpas-Burgos and Martinez-Beneito (2020) proposed more parsimonious parameterizations to : in Brewer and Nolan (2007), in Reich and Hodges (2008), and in Corpas-Burgos and Martinez-Beneito (2020). The Corpas-Burgos and Martinez-Beneito (2020) adaptive iCAR can be viewed as a re-parameterization of the Reich and Hodges (2008) adaptive iCAR, by letting , where is an additional scale parameter. The iCAR(a) in Table 3 is the Corpas-Burgos and Martinez-Beneito (2020) adaptive iCAR equivalent, with being replaced by .

With an adaptive iCAR, one can derive its extensions of adaptive LCAR or pCAR. For example, Corpas-Burgos and Martinez-Beneito (2020) extended their adaptive iCAR to two adaptive LCARs; the LCAR(d) and LCAR(e) in Table 3 (as extensions of the iCAR(a)) are the equivalents. The pCAR(c) in Table 3 is a pCAR extension of the iCAR(a).

Similar adaptive iCAR, pCAR and LCAR of this category can be formulated via linear transformation of non-adaptive CARs. For example, let where , and or or . Linear transformation (19), named linear (model of) coregionalization in MacNab (2018), leads to adaptive iCAR or pCAR or LCAR. The adaptive iCAR(b), pCAR(d), and LCAR(f) in Table 3 are examples.

A common feature of this category of CARs is that the adaptive parameters serve dual roles: They regulate both the spatial risk dependencies in the conditional mean functions and risk heterogeneities in the conditional variance functions (see Table 3). A noteworthy difference among them is their influence functions, as illustrated here for the iCAR(a) and iCAR(b): The influence function in (20) quantifies relative influence of area on area , relative to the neighbours of the area , in terms of its direct influence on the prediction of . When is near-zero, area exerts near-zero direct influence on its neighbours, but is not necessarily “disconnected” from its neighbours, which was suggested by Corpas-Burgos and Martinez-Beneito (2020). This is because area could still be directly influenced by its neighbours, say, when and , .

On the other hand, quantifies direct influence of area on area , which is partially controlled and partially determined by three quantities. It is an increasing function of but decreasing function of : For given , the higher the , the greater the area exerts direct influence on area . The greater the and/or the higher the , the less the area is influenced by its neighbours. Area exerts greater influence on area if and only if .

The LCAR(e) and LCAR(f) share a common feature that, when , they both reduce to adaptive independent Gaussian distributions defined by a diagonal variance matrix diag().

It is worth mentioning that the adaptive LCARs in the Category I are similar to the adaptive CARs of this category in the sense that their adaptive parameters regulate both spatial dependencies and heterogeneities.

In addition, in the context of Bayesian hierarchical modelling, the Rushworth et al. (2017) adaptive LCAR may be viewed as an adaptive CAR of this category. For example, a LCAR (extension) of the adaptive iCAR with is a Rushworth et al. (2017) adaptive LCAR equivalent. It is readily seen this LCAR (and the Rushwoth et al. model) may be excessively parameterized. To gain identifiability, the LCAR can be simplified by letting , which is the Corpas-Burgos and Martinez-Beneito (2020) adaptive LCAR(d) in Table 3 (with being replaced by ).

Adaptive CARs of the category III

Adaptive CARs of this category can be derived via linear transformation of adaptive CARs of the Category I. The adaptive pCAR and LCAR proposed in MacNab (2018) are examples; they can be readily derived via linear transformation (19), where or , pCAR() and LCAR() are adaptive CARs without scale parameter.

As the pCAR(e) in Table 3 illustrates, the coefficients of influence (in CARs of this category) are functions of adaptive spatial and scale (or precision) parameters; the adaptive scale parameters regulate risk heterogeneities in the conditional variance functions.

The coregionalized adaptive pCAR offers the flexibility for modelling locally varying spatial (Markovian) dependencies in the latent components , regulated by an adaptive spatial dependence parameterization in its conditional mean functions. Linear transformation (19) introduces locally varying risk heterogeneities to , regulated by the adaptively parameterized coefficients of coregionalization .

In the context of disease mapping, CARs of this category may be excessively parameterized. However, they may be useful when additional information, such as covariates, are available to infer these parameters (MacNab, 2018).

An illustrative example: Spatial and spatiotemporal modelling of COVID-19 infection risks

The purpose of this illustrative example is to show (potential) utilities of some of the disease mapping models and methods for analysing and monitoring communicable disease such as the COVID-19 infection risks during an ongoing epidemic or pandemic. I focus on potentials of the adaptive CAR models presented herein, in the context of Bayesian mapping of spatial and spatiotemporal COVID-19 infection risks, without or with covariates.

The overall analytic approach aimed for small area pandemic minoring by describing and understanding locally varying spatial risk dependencies, asymmetric local risk influence functions, within and between neighbourhood risks heterogeneities, and the resulting spatial risk discontinuities.

Details on models and related Bayesian methods for parameter estimation, learning, and inference, and details of the data analysed, and additional results and discussions, will be presented elsewhere. Here, for brevity, is a cursory report that outlines two sets of analysis, one is spatial modelling of county-level aggregates of daily infection cases for eighty-seven counties of Minnesota, USA, over the period of 2020/01/22–2021/02/14, named the cumulative period (data) hereafter. The second is spatiotemporal modelling of weekly aggregates of infection cases for the period of 2020/09/30–2021/01/20, a period the State-level 7-day averaged number of new cases exceeded 1000, named the peak period (data) hereafter (the State’s first major infection wave). I illustrate spatiotemporal modelling of county-level weekly aggregates of 17 weeks. The covariates presented here are tentative and only used for illustrative purpose; they are census-based aggregates and will be explained and further discussed in a separate paper.

Spatial analysis of the cumulative period data

The Poisson data model (1) with a regression part was fitted, in which the vector of log relative risks was modelled by each of the adaptive CARs in Table 3. Overall, the adaptive CARs led to comparable goodness-of-fit among the spatial models. Fig. 1 presents posterior estimates of the adaptive parameters and relative risks for two illustrative models; the results from the LCAR(a) illustrate spatial discontinuity, and those from the pCAR(b) illustrate locally varying neighbourhood risk dependencies. Fig. 1 and Table 4 suggest that the included (time-invariant) covariates explained modest (but noteworthy) amount of neighbourhood risk dependencies, modelled by the spatial dependence parameter ; they also suggest that the included covariates explained modest amount of risk variability, modelled by the scale parameter .

Fig. 1

Posterior estimates of the county-specific (87 counties) adaptive parameters and relative risks for indicated models. Red line or dot: without covariate, blue line: with five covariates. The Minnesota county-level COVID-19 cumulative period data. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Table 4

Posterior estimates, median and standard deviation (sd), of the model parameters for the adaptive pCAR(b). The Minnesota bcounty-level COVID-19 cumulative period data.

Parameter	Without covariate	With covariates
	Median sd	Median sd
Intercept	0.02 0.01	0.02 0.02
Private transportation to work		0.88 0.68
Age 55–64		−4.32 0.96
Education less than high school		3.94 0.70
College education		0.95 0.60
Unemployment		−4.56 1.35
λ	0.91 0.12	0.67 0.28
σ	0.40 0.03	0.32 0.03

Spatiotemporal analysis of the peak period data

The Poisson data model (1) was extended to a spatiotemporal data model where the regression part in Eq. (22) has time-varying coefficients for the included (time-invariant) covariates ().

The data were analysed via three sets of spatiotemporal models. I began with modelling the weekly data as time independent, and with weekly (i.e. time-varying) adaptive CARs, say, , , , where CAR (23) is an example of adaptive CARs of the Category I. The main purpose of the analysis was to describe and understand weekly changes in (i) locally varying spatial risk dependencies, (ii) asymmetric neighbourhood risk influences, and (iii) spatial risk heterogeneities and discontinuities.

As anticipated, and illustrated in Fig. 2 for indicated adaptive CARs, the spatial patterns of the adaptive parameters and relative risks varied over the 17 weeks. The included covariates explained modest amount of variations in the adaptive parameters and the relative risks during the week 7, which was among the peak weeks of this period. The observed patterns were also consistent with the posterior estimates of time-varying regression coefficients presented in Fig. 3, which illustrates a likely scenario of time-varying effects of covariates, such as the small area and neighbourhood characteristics included in the model. Fig. 3 also indicates that the included covariates led to modest reductions in the time-varying scale parameters over the peak weeks (the weeks 7 to 10).

Fig. 2

Illustrations of posterior estimates of county-specific (87 counties) adaptive parameters and relative risks for the indicated models. For the adaptive parameter estimates: Red lines—without covariate, blue lines—with covariates. For the relative risk estimates: Red dots: without covariate, blue lines: with covariates. The Minnesota county-level COVID-19 weekly data during peak period. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 3

Posterior estimates of the week-specific (17 weeks) regression coefficients and scale parameters for the indicated models. For the regression coefficients: dashed line and dots: posterior median, solid lines: lower and upper limits of 95% credible intervals. For the scale parameters: Red line: without covariate, blue line: with covariates. The Minnesota county-level COVID-19 weekly data during peak period. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 4 suggests modest variations in the estimated coefficients of influence for a selection of counties, using the results from the adaptive pCAR(d) (a model that had one of the lowest DIC, see Table 5). It also suggests some modest changes in the local risk influences from week 7 to week 17, and minor or modest changes from unadjusted to adjusted risk influences.

Fig. 4

Posterior estimates of the coefficients of influence for indicated counties, calculated using the posterior medians of the adaptive parameters in pCAR(d) (see Table 3). Red dots: Unadjusted risk influences—estimated coefficients of risk influence based on indicated ST model without covariate; blue dots: Adjusted risk influences—estimated coefficients of risk influence based on indicated ST model with covariates. The Minnesota county-level COVID-19 weekly data during peak period.

Table 5

DIC results for indicated models. The Minnesota COVID-19 peak period data.

Model	Without covariate	With covariates
	Dbar pD DIC	Dbar pD DIC
(1) pCAR(a)	10693 1237 11930	10687 1242 11929
(2) pCAR(b)	10687 1243 11930	10682 1247 11929
(3) pCAR(c)	10688 1223 11911	10690 1225 11915
(4) pCAR(d)	10692 1212 11904	10689 1214 11903
(5) LCAR(a)	10704 1227 11931	10692 1230 11922
(6) LCAR(b)	10694 1227 11921	10682 1230 11912
(7) LCAR(e)	10677 1228 11905	10676 1227 11903
(8) LCAR (f)	10697 1208 11905	10693 1212 11905
(9) SVC LCAR TpCAR1	10706 1222 11928	10706 1215 11921
(10) SVC LCAR TpCAR2	10671 1121 11792	10675 1120 11795

For the second set of analysis, the spatiotemporal model defined by Eqs. (21), (22) was fitted for eight options of time-invariant adaptive CAR parameterization. An example of the adaptive CARs of the category I is which is a simplification of Eq. (23). Table 5 presents the DIC results for the indicated adaptive CARs (1)–(8). Time- or variable-invariant adaptive CARs were previously used in Bayesian spatiotemporal models (Reich and Hodges, 2008, Rushworth et al., 2017, Sahu and Böhning, 2021) or multivariate spatial models (Corpas-Burgos and Martinez-Beneito, 2020). In this study, an intention of the analysis was to observe and compare the estimated time-invariant adaptive parameters with their counterparts of time-varying parameters. Fig. 5 illustrates the results for the pCAR(a). The patterns displayed from the estimated time-invariant spatial dependence parameters (the plots on the left) are notably different from those of the time-varying parameters for week 7 (the plots on the right). They convey different and important information: The estimated time-invariant may serve to explain and inform locally varying neighbourhood risk influences and spatial discontinuity over the first wave of the pandemic, while the weekly estimates of the adaptive parameters may serve to explain retrospectively, or monitor on going, temporal dynamics of locally varying risk dependencies, varying neighbourhood risk influences, and spatial discontinuities.

Fig. 5

Posterior estimates (median and standard deviation) of county-specific (87 counties) adaptive parameters for the indicated models. Red dots: Unadjusted risk influences—estimated coefficients of risk influence based on indicated ST model without covariate; blue dots: Adjusted risk influences—estimated coefficients of risk influence based on indicated ST model with covariates. The Minnesota COVID-19 weekly data during peak period.

The eight adaptive CAR models led to comparable and consistent posterior estimates of relative risks, although modest differences were observed. Overall, the CARs of the Category II performed slightly better than the CARs of the Category I. The LCARs also had slightly smaller DIC scores, compared with their pCAR counterparts. Of note is that the LCAR(e) (Corpas-Burgos and Martinez-Beneito, 2020) and the LCAR(d) and pCAR(f) of linear models of coregionalization (MacNab, 2018) led to comparable DIC results, although the LCAR(d) and pCAR(f) were shown to be more efficient, as indicated by the smaller pD scores and posterior relative risk standard deviations (not presented here). The overall DIC comparison, and the estimated adaptive parameters and relative risks, consistently suggested notable degrees of locally varying spatial dependencies, asymmetric local risk influences, and spatial discontinuities.

The third set of analysis concerned with CARs of the Category III, they are illustrative examples of exploring and explaining adaptive spatiotemporal risk interactions via a SVC model with spatially adaptive areal-specific temporal CAR processes: where , log and the vector of latent variables is modelled by a spatially adaptive temporal pCAR: where is a 17 by 17 first-order temporal connectivity matrix for “TpCAR1” and and/or for “TpCAR2”. The last model with the temporal pCAR2 led to the best DIC with notably improved statistical efficiency, as indicated by a notably reduced pD (in Table 5) and smaller posterior relative risk standard deviations (not presented here). The two sets of adaptive parameters and in the SVC models convey important information on asymmetric spatial and temporal risks dependencies and influences, which will be presented and discussed in a separate paper.

Discussion and looking into the future

The importance and popularity of Bayesian disease mapping in spatial epidemiology, population and public health, and medicine are seen in the large and still growing literature on spatial and spatiotemporal small area mapping of risks of COVID-19 infection and related health outcomes, and on applications of CAR models in COVID-19 related research. Identifications of high risk areas and patterns of disease incidence and mortality can provide valuable information for public health planning such as priority setting for allocating funds and localized disease prevention or intervention. In addition to mapping disease and illness, novel disease mapping models have been developed for spatial modelling of child growth (Gelfand and Vounatsou, 2003), geographical mapping of gene frequencies (Gelfand and Vounatsou, 2003), mapping health care utilization rates (Nathoo and Ghosh, 2012), and vaccination coverage mapping (Utazi et al., 2018), among others.

The (still on-going) COVID-19 pandemic, being an unprecedented global crisis, also offers unparallelled opportunities for the statistics community to put the established and newly developed disease mapping methods and tools, particularly the commonly used CAR, MCAR, and adaptive (M)CAR models, to rigorous assessment and test, to improve the existing ideas, methods, and tools, and to broaden the subject further. The colossal amount of COVID-19 related data available around the world are intrinsically spatial, spatiotemporal, multivariate (e.g. infection, hospitalization, mortality, uptake of vaccination, etc.), and multi-array (e.g. concerns spatial, temporal, and multivariate outcomes, in addition to other dimensions such as age, gender, risk/protective factors, and countries, for the varying public health intervention policies, updates of preventative measures, and availabilities of vaccines, among others). Making these data accessible to statistical analysis is critically important in our efforts to better understand the current (and future) pandemic and its impact on population and public health and to inform on pandemic intervention policies and strategies.

We see in issues of the Spatial Statistics, and (bio)statistical journals elsewhere, that innovative (multivariate) spatial and spatiotemporal disease mapping models have been developed to estimate, explain, and map the risks or spread of COVID-19 infection or related mortality (to name some most recent ones, Li and Dey, 2021, Slater et al., 2021, Huang et al., 2021; Feng, 2021, Lee et al., 2021, Sahu and Böhning, 2021) and for short-term forecasts of infection or infection-related mortality rates (Sahu and Böhning, 2021).

Disease mapping models are widely known for modelling rare and non-communicable diseases such as heart disease and cancer. CARs are commonly motivated for borrowing information and modelling smoothly varying disease risks as effects of omitted covariates. However, typical Bayesian disease mapping models are hierarchically formulated Bayesian spatial GLMMs that allows a variety of data models to be considered and applied. For example, in addition to the Poisson model (1) for rare events, we also have the options of binomial models for non-rare events Martuzzi and Elliott, 1996, MacNab, 2003b, frailty models for time-to-event data (Carlin and Banerjee, 2003, Lawson, 2018), and Gaussian models for continuous data (MacNab, 2021b).

When non-rare events are studied, (M)CARs and adaptive (M)CARs may be considered to allow for nuanced (intelligent) smoothing and for characterization and learning of (cross) spatial risks dependencies and/or (cross) spatial risk correlations, as well as spatial risk discontinuity. For instance, areal data of adequate numbers of events may contain useful information to inform on more flexible spatial dependence parameterization, perhaps aided by additional covariates. The works of [29] and Rushworth et al. (2017), and the present case study of spatiotemporal modelling of peak period COVID-19 infection risks, are examples.

Risks of communicable diseases such as the COVID-19 infection are quite likely to be spatially varying and to exhibit spatial discontinuity. This is due to the nature of communicable disease outbreak and transmission in proximities, differences in implementation and/or uptake of public health intervention policies and preventive measures, etc., in addition to the underlying health inequality that can partially influence the risks of infection or infection-related hospitalization or mortality. For these reasons, adaptive CARs might be more plausible options for modelling small area COVID-19 related data, without or with associated (protective) risk factors. The Sahu and Böhning (2021) work presents an example.

Here, I briefly illustrated another example of using adaptive spatial and spatiotemporal models to explore, explain, and monitor locally varying (cross) spatial dependencies, asymmetric local risk influences, and spatial discontinuities. Monitoring and identification of clusters of high local risks and risk influences could offer timely and valuable information that might assist and guide targeted public health interventions such as local restrictions and lockdowns, perhaps avoiding unnecessary state- or nation-wide restrictions and lockdowns.

In this article I gave a more detailed review on formulating adaptive CARs via options of adaptive spatial dependence or weight parameterization for a known connectivity matrix . One attraction of this approach is that the resulting models have clear interpretations of locally varying spatial risk dependencies and neighbourhood risk influences. Estimation and inference of the model parameters lead to estimation and inference of spatial dependence (local influence) functions, with associated estimation uncertainties. The adaptive spatial dependence parameterization discussed herein may mitigate the impact of a potentially unrealistic assumption of local risk dependencies defined by a given connectivity matrix, by allowing location-specific spatial dependence parameters to be zero, and thereby enables characterization of spatial risk discontinuity, as illustrated in [29]. In addition, some of the recent proposals of adaptive parameterizations of multivariate CARs, are natural extensions of the univariate adaptive CARs; the adaptive MGMRFs discussed in MacNab, 2018, MacNab, 2021a are examples.

The approach to allowing the connectivity matrix to be modelled as (binary) random variates within a CAR formulation leads to an alternative class of adaptive models. These models were initially motivated for, and were shown to be effective in, adjacency modelling and boundary detection (Lu and Carlin, 2005, Lu et al., 2007, Lee and Mitchell, 2012). While adaptive CARs of this class were also used for modelling spatial risk discontinuity (Lee et al., 2014, Sahu and Böhning, 2021), they lack the flexibility of modelling locally varying spatial dependencies that offered by adaptive spatial dependence parameterization. Nevertheless, the two classes of adaptive CARs are useful tools for their intended purposes and may be considered as complementary approaches in real-life applications. Adaptive CARs that accommodate both adaptive spatial dependence (or weight) parameterization and adjacency modelling are possible options to be further explored and developed.

The disease mapping community is poised to embrace new challenges and opportunities brought by the rise of data science and a (big) rich data era. Earlier disease mapping literature presented alternative ideas of modelling structured spatial discontinuities. Examples are hidden Markov models (Green and Richardson, 2002) and mixture models (Denison and Holmes, 2001, Fernandez and Green, 2002). These earlier ideas and methods have shown limited advantage in mapping rare disease. However, they may regain their recognition and appeal in the context of modelling (big) data of rich content, information, and structure. Together with (adaptive) CARs and their multivariate extensions, these models and related methods for estimation and inference may be further explored and better understood for their roles in disease mapping applications.

Recent advances in multivariate CAR/GRMFs, a topic of my focus here and a (new) comprehensive book has written for (Martinez-Beneito and Botella-Rocamora, 2019), provide ideas, models, methods, and tools for jointly modelling multidimensional spatial data in general, and the COVID-19 infection and related data in particular. Recent proposals of the INLA-SPDE approach to Bayesian analysis of spatial data over high resolution grids mark a new frontier of Bayesian disease mapping and its applications (Lindgren et al., 2011, Moraga et al., 2017, Utazi et al., 2018)

In addition, recent proposals of directed acyclic graph auto-regressive models in the context of (multivariate) disease mapping (Datta et al., 2019, Gao et al., 2021) offer alternative ways to model spatial dependencies. In the context of multivariate disease mapping, Liang (2019) proposed MCARs in which the non-spatial dependencies were characterized by a directed acyclic graph representation. The fusion of the Datta et al. (2019) and Liang (2019) ideas might broaden the scope and complexity of Bayesian disease mapping and widen its applicability in applied scientific (e.g. health science) research. The data science and big data movements are also inspiring new ideas of Bayesian disease mapping methods for handling data of high dimensions, e.g. areal data of large and/or and/or . A recent example is the idea of “divide and conquer” (Orozco-Acosta et al., 2021) for scalable approximate Bayesian inference of disease mapping of large using INLA.

Perhaps one of the most urgent need for advancing Bayesian disease mapping and its broad applications is devoted efforts for improved or new methods of Bayesian estimation, learning, and inference of (adaptive and/or multivariate) CAR. INLA offers a solution to computation for models of large but small or modest number of (prior) model parameters (Rue et al., 2017). Recent proposals of Gibbs sampling of GMRFs on large lattice (Marcotte and Allard, 2018, Brown et al., 2021) offer ideas and options to be explored, adapted, and improved. Alternative computational options that might be fruitfully explored and furthered are the hybrid or Hamiltonian Monte Carlo sampling methods (Girolami and Calderhead, 2011, Zhang et al., 2018, Hoffman et al., 2021) and the variational inference (Blei et al., 2017, Tan, 2020).