key: cord-0812964-5mj3cv6w
authors: MacNab, Ying C.
title: Bayesian disease mapping: Past, present, and future
date: 2022-01-19
journal: Spat Stat
DOI: 10.1016/j.spasta.2022.100593
sha: a0b4608ac6ee902baf65dc4a6e3ddd93ee69a9fc
doc_id: 812964
cord_uid: 5mj3cv6w

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 analyzing and monitoring communicable disease such as the COVID-19 infection risks during an ongoing pandemic.

On the occasion of the Spatial Statistics' 10th Anniversary, I reect 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 Alfred 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, 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 eorts 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 for spatial statistics 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 hereafter) 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 Montonoli (1992) , and Clayton et al (1993) inuenced 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 rates or risks maps to establishing a methodological paradigm with powerful ideas and tools for studying spatial and spatiotemporal variations, patterns, and clusters of disease and health outcome risks in populations of geographical regions, and for identication of (risk and protective) factors that might explain the variations and patterns and clusters. The body of Bayesian disease mapping models and related methods has grown in remarkable proportions 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 stratied).

The richness of the disease mapping models and their wide range of health science applications are well documented in books (Lawson et Lee 2011 , MacNab 2011 , 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 elds 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 analyzing communicable disease such as the COVID-19 infection risks in the context of pandemic intervention and monitoring. I reect 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.

2 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 statisticians was to develop model-based procedures for stabilized cancer risk maps, from which environmental determinants of dierent 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 empirical Bayes (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 Louis 2000a, 2000b and references therein). Disease mapping models of the past and present are models of this nature.

where the death counts are compared to a set of expected counts, denoted E= (E 1 ; E 2 ; :::; E n ), which is often derived for a set of area-and age-specic at risk population and the corresponding age-specic reference rates, known as age-standardization Kaldor 1987, Manton et al 1989) . The resulting relative risks, often dened as RR i = O i =E i ; Vi, 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-specic 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: O i $ Poisson( i ); i = E i ' i ; Vi; (1) where i ; Vi, is the Poisson intensity, and ' i is the unknown relative risk.

In frequentist disease mapping, the unknown relative risks in (1) are xed 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 f' i ; Vig, are the previously mentioned SMRs: ' i = O i =E i ; with associated standard errors se(' i ) = p' i =E i ; Vi: (2) It is well known that the MLEs (2) may lack stability and eciency 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 i quanties both the mean and the variance of the distribution. In 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 eld (GMRF, Besag 1974), named CAR for conditional autoregression, was proposed as a spatial prior for the log relative risks ensemble 2 = ( 1 ; :::; n ): 2 $ GMRF("; ); = 2 (s n );

where is the precision matrix, i = log(' i ); Vi, 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 prior postulates symmetric spatial dependence, accounts for spatial correlation, and facilitates spatial smoothing. 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 = 2 , 2 i = ( 1 ; :::; i 1 ; i+1 ; :::; n ), j $ i stands for the area j and i are neighbours. The area-specic conditional risk expectation of i given 2 i , Vi, 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 y. In Clayton and Kaldor (1987) , an EM algorithm (Dempster et al 1977) was used for a two-stage estimation of 2 (and 9) and its hyper-parameters.

The EB estimation is computationally straightforward for the rst example and more complex and evolved for the second one. In an inuential paper, Breslow and Clayton (1993) popularized several EB methods for inference in generalized linear mixed eects (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 9 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 rst example, where the uncertainties associated with the estimation of the hyperparameters may not be fully or adequately accounted for. This is the well known 

aive" EB estimation issue for nearly all EB approaches to estimation of 2 Gelfand 1990, 1991 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 previous mentioned gamma or CAR prior, can be postulated as an underlying process model for the underlying risks 2 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 software (Spiegelhalter et al 1995) 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 is a broad topic that concerns with spatiotemporal disease mapping (Held and Besag 1998 Here, I briey reect 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 2 = ( ij ), where j = 1; 2; :::; T for spatiotemporal disease mapping or j = 1; 2; :::; J for mapping of J diseases.

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

Some of the models are richly parameterized to hypothesize space and time interactions, in addition to decompose the 2 into (parameterized) components of space and time. Typical examples are the models of Waller et al (1997) and Held (2000) , two highly cited and inuential papers. For example, Held (2000) proposed a general model formulation it = t + ' t + i + i + it ; i = 1; 2; :::; n; t = 1; 2; :::; T; (7) and considered types of space-time interactions and associated prior options for f it ; Vi; tg. These spatiotemporal (interaction) models may be highly parameterized and, in this situation, addtional identication 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): it = t + i + it ; Vi; t; (8) in which functional analysis may be used to hierarchically model space ¢ time They illustrate modelling linear risk trends t = t and it = i t j , 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 ( Without sucient diculty, the spatiotemporal models (7) and (8) can be extended further for multi-array disease mapping. Adding age-eects is one example. Age-eects could be included in (7) or (8) Multivariate disease mapping has been an active area of research since the early 2000s. The main motivation of the eorts is to develop tools and methods for jointly modelling diseases that share common risk factors, enabling more eective 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 ij as a sum of two components: ij = i + ij ; Vi; j; (9) where i is a component shared by all diseases. This is the shared component model (SCM) proposed by Held and Best (2001) . SCMs typically hypthesize positive cross-variable local dependencies and correlations (MacNab 2010). The other more exible approach is to place a multivariate GMRF prior on vec(2) or vec(2 b ), which I discuss in Section 7.

In the context of disease mapping, including a regression part in a disease mapping model can serve serval important and practical purposes. One important purpose of a spatial regression model is to quantify disease risk variation explained and unexplained by the included risk (and/or protective) factors (Bernardinelli and Montomoli 1993 , MacNab 2003b ,c, 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 inuence 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 ( Ecological regression faces several methodological challenges. One of them is the long-standing ecological bias, refering the situation when the degree of association between risk exposure and disease at the areal-level is dierent from that of the individual level, due to within-area variability in exposures and confounders ( Another analytic challenge is spatial confounding, commonly discussed as collinearity between xed eects and random eects in a spatial GLMM, and in the presence of unmeasured confounding. Literature of ecological studies oered two main approaches to the issue. One is to deconfound the xed and random eects (Reich et Named TGMRF or TMGMRF, the basic idea of this approach is to model the relative risks 2, and its spatiotemporal or multivariate extensions, via a copula-inspired model structure that allows xed and random components to be estimated independently: The xed eects are modelled on the univariate margins, whereas the random eects 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 ( ,c, 2009 ) and, Section 7, for illustrative examples.

Of the deconfounding approach, earlier proposals are known as the restricted spatial regression (RSR) method (Reich et The key ideas of the RSR approach is to re-parameterize a n-dimension spatial prior and project the original Gaussian spatial prior onto a Gaussian spatial prior of its subspace orthogonal to the column space of the xed eects model matrix. Very recently, Zimmeraman 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 xed eects 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 eects and the included covariates are varying at dierent spatial scales and/or congurations ( This approach to deconfounding is also consistent with the Paciorek (2010) work on spatial confounding, which oered 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 tting 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 ( 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 Besag Table 1 for their full conditionals and associated precision matrices. Table 1 Three CARs commonly used in Bayesian disease mapping: = (w ij ); w ij = 1 when the ith and jth areas are neighbours (denoted i $ j) or w ij = 0 otherwise; w i+ = P j$i w ij ; hw = diag(w 1+ ; :::w n+ ). y : The covariance matrix ¦ is presented; 2 s $ iCAR( 2 s ); 2 h $ IIDN( 2 h ), yy : The covariance matrix ¦ is presented; 2 s $ iCAR( 2 ); 2 h $ IIDN( 2 ); 2 s c 2 h : 2 s and 2 h are (assumed) statistically independent.

For all models in Table 1 : The pCAR(; ) and LCAR(; ) are two-parameter and full rank GM-RFs, 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 dened on an irregular lattice of areal map, the pCAR, with its conditional expectation and precision functions

is a slight modication to the CAR (6), where fw i+ ; Vig in (10) is a set of scaling factors. As a result, the pCAR conditional expectation is an average of neighbourhood risks, rather than the sum of the neighbourhood risks as is dened in CAR (6) . And the pCAR (10) conditional precision is proportional to w i+ , 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 postulate asymmetric conditional dependency (i.e. direct inuence) of j on i , quantied by =w i+ versus i on j , dened by =w j+ , provided w i+ T = w j+ , i $ j, and 0 < < 1. That is, the direct inuence of the area j on its neighbouring area i is inversely proportional to the neighbourhood size of the area i. Put it another way, an area with a higher number of neighbours is less inuenced by its neighbour who has a lower number of neighbours. This could also be an intuitively plausible assumption, consistent with the conditional precision function Prec( i j2 i ) in (10); one might expect that an area of higher precision of risk prediction should be less inuenced 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 coecients of the conditional autoregression in (10) the coecient of inuence functions, or simply inuence functions, denoted Inuence(j; i) pCAR = =w i+ ; Vj $ i:

The LCAR conditionals also allow asymmetric spatial dependencies, provided Inuence(j; i) T = Inuence(i; j); i $ j, where Influence(j; i) LCAR = =(1 + w i+ ); Vj $ i:

That is, the coecient of inuence 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 n 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 i conditional on j ; Vj $ i (MacNab 2011). Noted in MacNab (2018), the LCAR parameterization is an ntangled" spatial and non-spatial parameterization, and as a result LCAR has limited options for multivariate generalization (also see Section 7 for a further discussion). The iCAR() is a GMRF with a singular precision matrix of rank n 1. It is the limiting distribution of the pCAR(; ) and LCAR (; ) (when 3 1).

The iCAR() is typically considered as a spatial prior for modelling spatially structure risk variation or clustered heterogeneity. It is commonly used when the n-vector of log relative risks 2 is modelled as additive components of 2 s $ iCAR( 2 s ), for modelling spatially structured heterogeneity or eects of omitted covariates that are spatially varying, and 2 h $ IIDN( 2 h ), for modelling extra-Poisson variation or eects of omitted covariates that are randomly varying. This is the well known Besag, Yorke, and Mollie (BYM) model (Besag 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 ( ( 2 ), where P (0; 1) is a weight parameter, weighting between the covariance matrix of the iCAR and the covariance matrix of the IIDN; also see Table 1 for the BYM and MBYM covariance matrices, respectively.

The LCAR has gained popularity due in part to the analytic result that when = 0, the LCAR reduces to independent and identical Gaussian priors for the log relative risks, whereas the pCAR, when = 0, reduces to n independent Gaussian priors with prior precisions proportional to w i+ ; Vi (MacNab 2006 , Lee, 2011 . However, compared to the LCAR, the pCAR has its own advantages. Its spatial and non-spatial parameters play separate and dierent 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, 2020, 2021a,b, see Sections 7 and 8).

Of the ve 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 dierent spatial dependence and correction functions (see details in MacNab 2011, 2014) and the ve 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 goodnessof-t assessment, it remains a viable and commonly taken approach to model comparison and selection.

In the broad context of multidimensional disease mapping discussed earlier, multivariate GMRFs play important roles in modelling multidimensional spa-tial dependencies and cross{covariances, enabling borrow-information over multidimensions, and facilitating spatial and cross-spatial smoothing. The ve 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. 1 ]. (gs) = hm s K ( gs), (g) = hm s K ( U g + b U g b ),(&; g) = sn & + (g), where hm = diag(m i ; :::; mn), gs = gs, gs is a symmetric matrix or diagonal matrix, g T = g b , U is the upper triangular part of ; ee b = ¦. y : 2 s $ MiCAR(¦s); 2 h $ IIDN(¦ h ); 2 s and 2 h are (assumed) statistically independent; yy : 2 s $ MiCAR(¦); 2 h $ IIDN(¦); 2 s k and 2 h k , Vk, 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 crossneighbourhood (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 (MacNab 2018). Working with a matrix of 2, one has a variety of ways to construct MGMRFs; see MacNab 2018 for a discussion and references therein.

Since the seminal works of Besag (1974) and Mardia (1988) , MGMRFs are often formulated via full conditionals (e.g. Kim 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 (g) = h m s p ( U g + b U g b ); g T = g b ; (11) provided (g) > 0, where g is a p ¢ p full matrix of spatial dependence parameters, h m = diag(m 1 ; m 2 ; :::; m n ), fm i ; Vig are area-specic scaling factors, U is the upper triangular part of , is the 

eighbourhood" connectivity or weight matrix. (g) is a spatial dependence matrix for locally dependent MGMRF constructions. Expression (11) is simplied to (g s ) = h m s p ( g s ) when g s is a symmetric matrix or a diagonal matrix of spatial dependence paramters. It also reduces to a separable model when g = cs p .

Without the ntanglement" challenges, MGMRFs of construction (11) can be readily formulated within the Besag (1974) ¦ cMGMRF vec (2 > ) (g; e) = (s n e)(g) 1 (s n e); (13) where e is a K ¢ K coregionalization coecients matrix, ee b = ¦ > 0, ¦ is a covariance matrix.

Notice that (12) and (13) exemplify one way to mix spatial and non-spatial components that characterizes and parameterizes cross-spatial interactions.

Alternatively, let the latent elds be vec( b ) $ MVN(0; s n ¦), a spatially varying (non-stationary) coregionalization (SVC) construction may be formulated via vec(2 b ) = r(g)vec( b ); (14) to have the following covariance matrix (MacNab 2018): ¦ vec (2 > ) (g; ¦) = r(g)(s N ¦)r(g) b ; (15) where r(g)r(g) b = (g) 1 . In SVC (15) the latent variables are only correlated at same locations. SVC constructions with latent (M)GMRF can also be readily formulated.

The mixing schemes of (12) and (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, 2021a,b) .

As an alternative to (12) and (14), a cMGMRF can be formulated via

to have the following precision matrix: (13) and (17).

With options of parameterization for (g), &, e, ¦, the three constructions contain all major proposals of MGMRFs to date and much more, notably 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. cMGM-RFs of constructions (13) and (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, 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, 2020, 2021a,b, Prate et al 2021) .

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

GMRFs are undirected graphical models, commonly dened 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 dened by the matrix of a given map and its neighbourhood denition. 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 exibility to allow within-and between-neighbourhood risk heterogeneities that likely happen in parts of the map ( 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 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 inuences and heterogeneities that typically manifest macro-level phenomenons such as the spatial risk discontinuity. There are important dierences between the adaptive pCARs and LCARs of this category. A common feature of the adaptive pCARs is that the adaptive parameters fc i ; Vig 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 inuence functions, as illustrated here for the pCAR(a) and pCAR(b):

Inuence(j; i) pCAR(a) = c i w i+ ; Inuence(j; i) pCAR(b) = c i c j w i+ : (18) In pCAR(a), quantify how each area is directly inuenced by its neighbours:

The higher the c i , the greater the i prediction is directly inuenced by its Table 3 Table 3 ). Similar to its non-adaptive counterpart, its coecient of inuence functions imply asymmetric pair-wise local dependencies and inuences, provided w i+ T = w j+ : An area with a higher number of neighbours is less inuenced by its neighbour who has a lower number of neighbours. Two neighbours exerts same inuence 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 duel roles of regulating locally varying spatial dependencies and heterogeneities (see Table 3 ).

For example, the Inuence(j; i) LCAR(a) = c i =(1 + c i (w i+ 1)); Vj $ i, is a nonlinear increasing function of c i , the corresponding conditional variance function is a non-linear decreasing function of c i . 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 inuence functions for interpreting locally varying spatial dependencies and inuences (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 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 = (w ij ; Vi $ j) in Table 3 , wherew ij ; Vi $ j, models locally weighted pair-wise interaction between i and j . In Brezger et al (2007) , for example, the elements fw ij ; Vi $ jg were estimated via a gamma prior. Brewer Table 3 is the Corpas-Burgos and Martinez-Beneito (2020) adaptive iCAR equivalent, with c 1=2 i being replaced by c i .

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)). 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 2 ='~ ;2 $ CAR; (19) where' = diag( 1 ; 2 ; :::; n ), and2 $ iCAR(1) or2 $ LCAR() or 2 $ pCAR(). 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 dierence among them is their inuence functions, as illustrated here for the iCAR(a) and iCAR(b):

The inuence function Inuence(j; i) iCAR(a) in (20) quanties relative inuence of area j on area i, relative to the neighbours of the area i, in terms of its direct inuence on the prediction of i . When c j is near-zero, area j exerts near-zero direct inuence on its neighbours, but is not necessarily \disconnected" from its neighbours, which was suggested by Corpas-Burgos and Martinez-Beneito (2020). This is because the area j could still be directly inuenced by its neighbours, say, when c i > 0 and i $ j, Inuence(i; j) iCAR(a) = c i = P i$j c i > 0. On the other hand, Inuence(j; i) iCAR(b) quanties direct inuence of area j on area i, which is partially controlled and partially determined by three quantities. It is an increasing function of c j but decreasing function of c i : For given c i , the higher the c j , the greater the area j exerts direct inuence on area i. The greater the c i and/or the higher the w i+ , the less the area i is inuenced by its neighbours. Area j exerts greater inuence on area i if and only if c 2 j w j+ > c 2 i w i . The LCAR(e) and LCAR(f) share a common feature that, when = 0, they both reduce to adaptive independent Gaussian distributions dened by a diagonal variance matrix diag( 2 ).

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) Table 3 illustrates, the coecients of inuence (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 oers the exibility for modelling locally varying spatial (Markovian) dependencies in the latent components2, regulated by an adaptive spatial dependence parameterization in its conditional mean functions. Linear transformation (19) introduces locally varying risk heterogeneities to 2, regulated by the adaptively parameterized coecients 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). 9 An illustrative example: Spatial and spatiotemporal modelling of COVID-19 infection risks in counties of Minnesota, USA The purpose of this illustrative example is to show (potential) utilities of some of the disease mapping models and methods for analyzing 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 inuence functions, within and between neighbourhood risks heterogeneities, and the resulting spatial risk discontinuities over the map.

Details on models and related Bayesian methods for parameter estimation, learning, and inference, and details of the data analyzed, 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 countylevel 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 rst 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. Table 3 . Overall, the adaptive CARs led to comparable goodness-of-t 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 discontinuities, 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 . 

where log(' it ) = 0t + 1t X it1 + ::: + 1t X itp + it ; p = 5;

the regression part in Equation (22) has time-varying coecients for the included (time-invariant) covariates (p = 5). The data were analyzed 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, 2 t $ adaptive CAR( t ; t ); (23) 2 t = ( 1t ; 2t ; :::; nt ) b , t = (c 1t ; c 2t ; :::; c nt ) b , Vt, 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 inuences, and (iii) spatial risk heterogeneity and discontinuity.

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 coecients presented in Fig.  3 , which illustrates a likely scenario of time-varying eects 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). Bayesian disease mapping Table 5 ). It also suggests some modest changes in the local risk inuences from week 7 to week 17, and minor or modest changes from unadjusted to adjusted risk inuences. For the second set of analysis, the spatiotemporal model dened by Equations (21) and (22) was tted for eight options of time-invariant adaptive CAR parameterization. An example of the adaptive CARs of the category I is 2 t $ adaptive CAR(; ); Vt; (24) q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q County B q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q County D q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q County E q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q County F q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q County B q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q County D q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q County E q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q Table 3 ). Red dots: Unadjusted risk inuences -estimated coecients of risk inuence based on indicated ST model without covariate; blue dots: Adjusted risk inuences -estimated coecients of risk inuence based on indicated ST model with covariates. The Minnesota county-level COVID-19 weekly data during peak period.

which is a simplication of the Equation (23) . 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 dierent from those of the timevarying parameters for week 7 (the plots on the right). They convey dierent and important information: The estimated time-invariant may serve to explain and inform locally varying neighbourhood risk inuences and spatial discontinuities over the rst 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 inuences, and spatial discontinuities. The eight adpative CAR models led to comparable and consistent posterior estimates of relative risks, although modest dierences 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 ecient, as indicated by the smaller pD scores and posterior relative risk standard deviations (not presented Bayesian disease mapping 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 inuences, 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-specic temporal CAR processes: vec(2) = (s T diag('))vec(); (25) where' = (s 1 ; :::; s n ), log(s) $ LCAR( s ; ) and the vector of latent variables vec() is modelled by a spatially adaptive temporal pCAR: vec() $ pCAR(; T ); (26) where T is a 17 by 17 rst-order temporal connectivity matrix t $ t + 1 for \TpCAR1" and t $ t+1 and/or t $ t 1 for \TpCAR2". The last model with the temporal pCAR2 led to the best DIC with notably improved statistical eciency, 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 s and in the SVC models convey important information on asymmetric spatial and temporal risks dependencies and inuences, which will be presented and discussed in a separate paper.

10 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. 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 oers unprecedented 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). Making these data accessible to statistical analysis is critically important in our eorts 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 Disease mapping models are widely known for modelling rare and noncommunicable chronic diseases such as heart disease and cancer. CARs are commonly motivated for borrowing-information and modelling smoothly varying disease risks as eects 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 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 exible spatial dependence parameterization, perhaps aided by additional covariates. The works of Condgon (2008) 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, dierences in implementation and/or uptake of public health intervention policies and preventive measures, etc, in addition to the underlying health inequality that can partially inuence 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 Bohning (2021) work presents an example.

Here, I briey illustrated another example of using adaptive spatial and spatiotemporal models to explore, explain, and monitor locally varying (cross) spatial dependencies, asymmetric local risk inuences, and spatial discontinuities. Monitoring and identication of clusters of high local risks and risk inuences could oer 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 appeal of this approach is that the resulting models have clear interpretations of locally varying spatial risk dependencies and neighbourhood risk inuences. Estimation and inference of the model parameters lead to estimation and inference of spatial dependence (local inuence) 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 dened by a given connectivity matrix, by allowing location-specic spatial dependence parameters to be zero, and thereby enables characterization of spatial risk discontinuity, as illustrated in Condgon (2008). In addition, some of the recent proposals of adaptive parameterizations of multivariate CARs, are nature extensions of the univariate adaptive CARs; the adaptive MGMRFs discussed in MacNab (2018, 2021a) are examples.

The approach to allowing the connectivity matrix to be modelled as binary variates within a CAR formulation leads to an alternative class of adaptive models. These models were initially motivated for, and were shown to be eective 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 Bohning 2021), they lack the exibility of modelling locally varying spatial dependencies that oered 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 . 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. 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) oer 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 scientic (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 n and/or p and/or T . 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 n using INLA.

Perhaps one of the most urgent need for advancing Bayesian disease mapping and its broad applications is devoted eorts for improved or new methods of Bayesian estimation, learning, and inference of (adaptive and/or multivariate) CAR. INLA oers a solution to computation for models of large n but small or modest number of (prior) model parameters (Rue et al 2017) . Recent

Mapmaking and mapthinking: cancer as a problem of place in nineteenth-century England

The Cancer Problem: Malignancy in nineteenthcentury Britain

Approximate inference for disease mapping

Neighborhood dependence in Bayesian spatial models

Mspock: Alleviating spatial confounding in multivariate disease mapping models

A Gaussian random eld mode for similarity-based smoothing in Bayesian disease mapping

Hierarchical modeling and analysis for spatial data

Empirical bayes versus fully bayesian analysis of geographical variation in disease risk

Disease mapping with errors in covariates

Spatial interaction and the statistical analysis of lattice systems (with discussions)

Bayesian image restoration, with two applications in spatial statistics

A comparison of Bayesian spatial models for disease mapping

Variation inference: a review for statisticians

Volume I -The analysis of case-control studies. Lyon: International Agency for Research on Cancer

Extra-Poisson Variation in Log-Linear Models

Approximate Inference in Generalized Linear Mixed Models

Variable smoothing in Bayesian intrinsic autoregressions

Adaptive Gaussian Markov random elds with applications in human brain mapping

Sampling strategies for fast updating of Gaussian Markov random elds

Bayesian modelling for spatially misaligned health and air pollution data through the INLA-SPDE approach

Approaches for empirical Bayes condence intervals

A sample reuse method for accurate parametric empirical Bayes condence intervals

Bayes and empirical Bayes methods for data analysis

Empirical Bayes: Past. Present and Future

Hierarchical multivariate CAR models for spatio-temporally correlated survival data (with discussion)

Empirical Bayes estimates of agestandardized relative risks for use in disease mapping

Bayesian methods for mapping disease risk

Spatial Correlation in Ecological Analysis

Bayesian statistical modelling

A spatially adaptive conditional autoregressive prior for area health data

On the use of adaptive spatial weight matrices from disease mapping multivariate analyses

Statistics for spatial data

Statistics for spatio-temporal data

Spatial disease mapping using directed acyclic graph auto-regressive (DAGAR) models (2019)

Maximum Likelihood from Incomplete Data via the EM Algorithm

Bayesian partitioning for estimationing disease risks

On predicting cancer mortality using ANOVA-type P-spline models

Skew-elliptical spatial random eect modeling for areal data with application to mapping health utilization rates

Spatial-temporal generalized additive model for modeling COVID-19 mortality risk in Toronto

Modelling spatially correlated data via mixture: A Bayesian approach

Bayesian analysis of areal data with unknown adjacencies using the stochastic edge mixed eects model

Multivariate Directed Acyclic Graph Auto-Regressive (MDAGAR) models for spatial diseases mapping

Sampling-based approaches to calculating marginal densities

Proper multivariate conditional autoregressive models for spatial data analysis

Nonstationary multivariate process modelling through spatially varying coregionalization (with discussion)

Bayesian Data Analysis

Markov chain Monte Carlo in practice

Riemann manifold Langevin and Hamiltonian Monte Carlo methods (with discussions)

Agespacetime CAR models in Bayesian disease mapping

A multivariate CAR model for improving the estimation of relative risks

Hidden Markov models and disease mapping

The geographical distribution of disease in England and Walses

The geographical distribution of cancerous disease in the British Isles. The Lancet

Modelling risk from a disease in time and space

Bayesian modelling of inseparable spacetime variation in disease risk

A shared component model for detecting joint and selective clustering of two diseases

An Adaptive MCMC Scheme for Setting Trajectory Lengths in Hamiltonian Monte Carlo

Population-weighted exposure to air pollution and COVID-19 incidence in Germany

Dimension reduction and alleviation of confounding for spatial generalized linear mixed eects models

copcar: A exible regression model for areal data

Estimating the changing nature of Scotland's health inequalities by using a multivariate spatiotemporal model

Order-free co-regionalized areal data models with application to multiple-disease mapping

A bivariate Bayes method for improving the estimates of mortality rates with a twofold conditional autoregressive model

Disease mapping and risk assessment for public health

Bayesian Disease Mapping Hierarchical Modeling in Spatial Epidemiology

A comparison of conditional autoregressive models used in Bayesian disease mapping. Spatial and spatio-temporal epidemiology

Boundary detection in disease mapping studies

CARBayes: an R package for Bayesian spatial modeling with conditional autoregressive priors

A Bayesian localized conditional autoregressive model for estimating the health eects of air pollution

Controlling for unmeasured confounding and spatial misalignment in long-term air pollution and health studies

Quantifying the small-area spatiotemporal dynamics of the Covid-19 pandemic in Scotland during a period with limited testing capacity

P-spline ANOVA-type interaction models for spatio-temporal smoothing

Graph-based multivariate conditional autoregressive models

Estimation of disease rates in small areas: a newmixedmodel for spatial dependence

Estimation of COVID-19 mortality in the United States using Spatio-temporal Conway Maxwell Poisson model

An explicit link between Gaussian elds and Gaussian Markov random elds: the stochastic partial dierential equation approach

Bayesian areal Wombling for geographical boundary analysis

Bayesian areal Wombling via adjacency modelling

Autoregressive spatial smoothing and temporal spline smoothing for mapping rates

A Bayesian hierarchical model for accident and injury surveillance

Hierarchical Bayesian spatial modelling of small-area rates of non-rare disease

Hierarchical Bayesian modelling of spatially correlated health service outcome and utilization rates

Bayesian spatial and ecological models for smallarea accident and injury analysis

Estimation in Bayesian disease mapping

An innovative application of Bayesian disease mapping methods to patient safety research: the Canadian iatrogenic injury study

Spline smoothing in Bayesian disease mapping

Regression B-spline smoothing in Bayesian disease mapping: with an application to patient safety surveillance

On empirical Bayes penalized quasi-likelihood inference in GLMMs and in disease mapping and ecological modeling

Bayesian multivariate disease mapping and ecological regression with errors in covariates

On Bayesian shared component disease mapping and ecological regression with errors in covariates

On Gaussian Markov Random Fields and Bayesian disease mapping

On identication in Bayesian disease mapping and ecological-spatial regression

Linear models of coregionalization for multivariate lattice data: a general framework for coregionalized multivariate CAR models

Linear Models of coregionalization for multivariate lattice data: Order-dependent and order-free MCARs

Some recent work on multivariate Gaussian Markov Random Fields (with discussions)

Bayesian estimation of multivariate Gaussian Markov random elds with constraint

On coregionalized multivariate Gaussian Markov random elds -Part I: Construction and parameterization

On coregionalized multivariate Gaussian Markov random elds Part II: Bayesian estimation and inference, with simulation and case studies

Empirical Bayes Procedures for Stabilizing Maps of U.S. Cancer Mortality Rates

Multi-dimensional multivariate Gaussian Markov random elds with application to image processing

Gibbs sampling on large lattice with GMRF

A general modeling framework for multivariate disease mapping

Towards a Multidimensional Approach to Bayesian Disease Mapping

Disease mapping: From foundations to multidimensional modeling

Some link between conditional and coregionalized multivariate Gaussian Markov random elds

Empirical Bayes estimation of small area prevalence of non-rare conditions

Generalized linear models (second Edition)

A geostatistical model for combined analysis of point-level and area-level data using INLA and SPDE

Scalable Bayesian modelling for smoothing disease risks in large spatial data sets using INLA

The importance of scale for spatial-confounding bias and precision of spatial regression estimators

Transformed Gaussian Markov random elds and spatial modeling of species abundance

Alleviating spatial confounding for areal data problems by displacing the geographical centroids

Non-separable spatio-temporal models via transformed Gaussian Markov random elds

Eects of residual smoothing on the posterior of the xed eects in disease-mapping models

Modeling longitudinal periodontal data: A spatially-adaptive model with tools for specifying priors and checking t

Interpreting Posterior Relative Risk Estimates in Disease-Mapping Studies

Gaussian Markov Random Fields -Theory and Applications

An adaptive spatiotemporal smoothing model for estimating trends and step changes in disease risk

A spatial analysis of multivariate lattice data

Bayesian statistics without tears: a samplingresampling perspective

Bayesian spatio-temporal joint disease mapping of Covid-19 cases and deaths in local authorities of England

A spatial analysis of multivariate output from regional climate models

Capturing spatial dependence of COVID-19 case counts with cellphone mobility data

Posterior distribution of hierarchical models using CAR(1) distributions

Use of model reparametrization to improve variational Bayes

GeoBUGS User Manual

Bayesian latent variable modelling of multivariate spatio-temporal variation in cancer mortality

Empirical bayes estimation of cancer mortality rates

Mixed model for analyzing geographic variability in disease rates

Evaluating the performance of spatio-temporal Bayesian models in disease mapping

Spatiotemporal modeling of mortality risks using penalized splines

One-dimensional, two-dimensional, and three dimensional B-splines to specify space-time interactions in Bayesian disease mapping: Model tting and model identiability

A spatial regression model for the disaggregation of areal unit based data to high-resolution grids with application to vaccination coverage mapping

Bayesian inference in multivariate spatio-temporal areal models using INLA: analysis of gender-based violence

Multivariate Bayesian spatiotemporal P-spline models to analyze crimes against women

Spatial dependence and errors-in-variables in environmental epidemiology

A statistical framework for ecological and aggregate studies

Disease mapping and spatial regression with count data

Applied spatial statistics for public health data

Hierarchical Spatio-Temporal Mapping of Disease Rates

Variational Hamiltonian Monte Carlo via score matching

On Deconfounding Spatial Confounding in Linear Models