Spatial Statistical Models: an overview under the Bayesian Approach

Spatial documentation is exponentially increasing given the availability of Big IoT Data, enabled by the devices miniaturization and data storage capacity. Bayesian spatial statistics is a useful statistical tool to determine the dependence structure and hidden patterns over space through prior knowledge and data likelihood. Nevertheless, this modeling class is not well explored as the classification and regression machine learning models given their simplicity and often weak (data) independence supposition. In this manner, this systematic review aimed to unravel the main models presented in the literature in the past 20 years, identify gaps, and research opportunities. Elements such as random fields, spatial domains, prior specification, covariance function, and numerical approximations were discussed. This work explored the two subclasses of spatial smoothing global and local.

reshaping daily tasks, in which miniaturization devices are also placing location labels much easier than before [68].Spatial dependencies have long been identified as a component that could hinder model precision and increase bias.
Subsequent effort to account for such error created a research line in spatial statistics.
Observation oriented across space is an essential feature for the IoT data, which influences data prediction and analysis.The applications vary in complexity, and it is frequently applied in risk surfaces detection, healthcare, agriculture, urban planning, economics, and rarely applied in smart appliances that learn based on location.The complex structure is accommodated in a flexible class of models related to the observed data and the spatial dependencies.
The frequentist (classical) and the Bayesian analytical methods have been used to analyze IoT data.However, the Bayesian method is a better choice because it owes its ability to accommodate information from different sources.
In the Bayesian framework, elucidated questions are answered in the estimation procedure by combining multiples sources of information, such as previous knowledge (prior) and the acquired information in the data (likelihood) [81].
In the neuroscience field, neurorehabilitation has growing and technological improvements made given the neuronavigation, allowing personalizing the definition of transcranial accuracy [38,69].Recently, there is an increasing interest in using spatial models on the meta-analysis of brain imagery to locate hot spot regions of consistent activation on the brain for diagnostic and treatment [18,45].Another exemplification is in epidemiology, which takes advantage of patterns across geographical space to identify the areas of potentially elevated risk and create disease maps to quantify underlying risk surface [47].Moreover, the medical literature provides detailed motivation and descriptions of spatial smoothing methods by explaining the concepts, defining the technical terms, and demonstrating various visualizing spatial models.
The basic idea of the Bayesian spatial statistics is the extension of the generalized linear model, including a spatial component that accounts for spatial dependencies across a spatial domain.The component is assigned a spatial prior, usually multivariate, which accounts for spatial correlation across a region (not necessarily delimited).Afterward, the parameter estimates are smoothed across the spatial domain with a specified resolution to identify the hot spot region and provide intuition on the chain of events.Moreover, the approach is different for frequently used spatial econometric models that treat the spatial dependencies as a global correlation parameter across a spatial domain.In the heart of every spatial model have a correlation matrix that quantifies the dependencies of the spatial component and determines the complete distribution of the spatial prior.The correlation matrix, proportional to the weight matrix, has an enormous impact on spatial smoothing, and the challenge usually faced with authors is the choice of its specification.
Additionally, the model complexity can be owed to the structure of the weight matrix.A highly dense weight matrix implies high correlated spatial field and a more complex model.Besides the general specification of the diagonal element of the weight matrix set to zero, over the past 20 years, there has not been a concrete documented standard on the specification of the off-diagonal elements.The most frequently used weight matrix is the binary first neighborhood structure, which assigns 1 or 0 depending on whether the spatial locations are immediate neighbors or not.

| Objectives of this review
A Bayesian spatial statistic aims to quantify the spatial pattern and provide insight into the process generating the pattern.Given the technological advances and precisely the storage capacity, data georeferenced acquisition is commonly present in the nowadays domain sets.The Bayesian method is typically used to analyze these sets and to identify the spatial pattern.Consequently, Spatial Statistics have received considerable attention in recent years, and numerous spatial models have been proposed and applied in diverse research fields.The systematic review of these models has received little attention and specifically been conducted in epidemiology [2,21,100].This systematic review focused on the progressive development and the content analysis of the Bayesian spatial models Queries of the word "Bayesian Spatial" and "Bayesian spatial", using the Boolean operator "OR", through the year 2001-2020 was conducted.Title, abstract, and Keywords, were used in Scopus and Science Direct, topic (which entails title, abstract, and keywords) in Web of Science, and "Anywhere" in MathSciNet.
The time frame was chosen to capture the diversification of the application of Bayesian spatial models in various fields, Owing to the advancement in technology for data collection and computation.Mendeley Windows application was used to remove duplicated articles.The resulting set was further examined manually for more duplicates not identified by Mendeley's application.The titles and abstracts of articles included after removing duplicates were first screened for Bayesian spatial methodology before applying the following inclusion criteria.
• Search results that are written in English, an article published in peer-reviewed journals available online We exclude books, dissertations, thesis, conference proceedings, reviews (or any other form, not an article).
• Articles that specifically implement Bayesian spatial models excluding articles that only mentioned Bayesian spatial models.
Articles that did not meet the two inclusion criteria were excluded from the review.The search flow chart is presented in Figure 1.Using the search keywords earlier mentioned, 586 articles were retrieved from Scopus, 129 from Science direct, 492 from Web of Science, and 73 in MathSciNet.After the exclusion of duplicated articles, 590 articles were assessed for eligibility, and 38 were further excluded based on the two exclusion criteria, leaving 552 articles As a structure of the data set, 552 articles that met the eligibility criteria were classified into the following categories.

• Names of all authors,
• Publication year, • Journal title and • Response to the ten items of the conceptual classification scheme on Bayesian spatial models.This survey was divided into two parts theoretical models and empirical analyses of the published articles.These analyses scrutinized results from the last 20 years in the next section.In the first part, we discussed the different approaches under the statistical innovation and their differences, divided into eight topics; Fields of application, Spatial domains, Spatial Priors, response variables, statistical models, Prior specification, estimation methods, Simulation, and validation.We presented the various applications adopted by the existing spatial models from the systematic methodological review in part two.

C O N C E P T U A L S C H E M E F O R S PAT I A L M O D E L S
This research focused on content analysis in Bayesian Spatial model to systematically assess the content of a large volume of recorded information in this field.It aims to provide a deep insight into the contributions and identify the key research directions and opportunities in the Bayesian spatial methodology.To accomplish this, we applied a conventional approach to content analysis [40] by scrutinizing through samples of the articles to clearly define the characteristics that better explain the scope and richness of the literature, identify the key concepts and patterns.This stands as the initial characteristics clustering.As we inspect the entire articles, subsequent updates on the clustering were included when new data that did not fit into the defined characteristics were encountered.This approach gives room for the literature to be classified without presumption.
Every Bayesian spatial analysis aims to estimate the spatial pattern over an extended geographical region to identify regions with extreme realization.In Bayesian hierarchical models, the spatial pattern is represented with a component that uses the same set of smoothing parameters across the entire study region.This type of smoothing is referred to as global smoothing.In some geographical setting, a global smoothing may be inappropriate, owing to the complexity of the geographical setting, and the spatial pattern is likely to exhibit a localized behavior.Thus, localized regions are smoothed using different parameters, and it is referred to as local smoothing [56].
TA B L E 1 Summary of the Spatial Models and its variation.Mixture Model Green and Richardson [34] Spatial Partition Model Leonhard and Raßer [50] Asymmetric Laplace Kuzobowski and Pogorski [51] Student-t Fonseca [25] Log-Gamma Bradley et al .[10] Dirichlet Gelfand et al., 2005 [30] For instance, let's suppose as a special case of the spatial models, a general framework of the parametric spatial model given response variable Y , for instance, let's consider θ as a generalization of the response variable, and a set of a linearly related covariate is formulated as where θ is the transformed quantity of interest, g (.) is a link function, X is the fixed effect design matrix, β is the unknown fixed effects vector, Z is the random effect design matrix (particularly, the spatial effect), and γ is the latent spatial variable to be estimated.The probability distribution of Y, Y ∼ f y (. |θ), determines the function g (.).In a frequentist estimation procedure, it involves maximizing the joint log-likelihood l (y , γ |β ).In a Bayesian framework, all the parameters are considered random, which can be guided by an informative structure (prior).
The spatial effect is assigned a spatial prior γ ∼ f sp at (. |φ, Σ(α)), where φ is the vector of hyperparameters, and Σ(α) is a covariance matrix that determines the spatial dependencies of γ across a spatial domain with smoothing parameter α.The f sp at (.), usually multivariate, assumes a parametric distribution (not necessarily an exponential family) such as Gaussian, Asymmetric Laplace, Student-t, Log-Gamma, and more.The α is a global smoothing parameter if the same set of α smooths the entire spatial region, whereas it is a local smoothing if different sets, represented by the vector α, smooth the spatial region.Additionally, the spatial process can be modeled in a semi-parametric framework, such as the spatial mixture model, and non-parametric, such as the Dirichlet process.
In the review process, as a research methodology, we first classify each article into one of the two disjoint classes, "Theoretical" and "Applied".The theoretical methods involve investigating fundamental principles and reasons for the occurrence of an event, random phenomenon, or processes.On the contrary, applied research involves solving a particular problem with known or accepted theories and principles.

| Spatial Statistics Fields of Application
Bayesian spatial statistics is a useful tool to determine the dependence structure and hidden patterns over space, through the prior knowledge and data likelihood.In some cases, the hypotheses of interest of a random phenomenon do not directly relate to the effect of spatial dependencies.However, it is crucial to adjust for spatial variation [82].
Adjustment for spatial patterns in modeling random occurrence has since been practiced across various fields such as Agriculture, Medicine, Biology, Epidemiology, Geography, Geology, Economics, Climatology, Ecology, among others [49].Moreover, spatial dependence in the Agriculture experiment has long received consideration.RA.Fisher identified spatial variations and used blocking to mitigate the effect of spatial dependencies in a randomized experimental design [4,24].
In many Biological and Medical experiments, such as gene classification, brain mapping, the randomized blocking technique may not be a viable alternative.Moreover, in demography, disease mapping, image analysis, remote sensing, fabrication engineering, and species detection, the variation due to spatial proximity cannot be neglected.It may result in bias and inconsistent estimates.Responses at close range tend to have similar behavior and variation.The homogeneity of the variation depreciates with the increased distance apart.An efficient procedure to tackle the effect of spatial proximity is to consider random field Statistical models.Random field Statistical models, known as spatial models, describe the distribution of a random phenomenon over a spatial domain.
Spatial models have long be applied in various fields.In 1949, Isard described the general theory of the spatial formation of economic activities focusing on the geographic distribution of costs, prices, and location of industries [43].
Spatial statistics applied to economics, often referred to as spatial econometrics, have gained more attention in recent years to analyze economic data over a wide range of spatial domain [92].Similarly, In 1950, DA.Krige took advantage of nearby variations to pursue the spatial prediction of gold distribution in South Africa, basing predictions practically on lognormal-de Wijsian spatial models [52].In Epidemiology and Public Health, spatial statistics have gained increasing importance in predicting disease outbreak [33,65,72,70,89].The problems that arose from these fields lead to the motivation of several intuitions that gave birth to the consistent improvements in the spatial model literature.Some of the most advances are identifying spatial risk factors, disease surveillance, and spatial predictive models, [101].
In this research, the fields of the application were classified into five major parts.

| Spatial Domains
Geographically reference data, also known as spatial data, is a collection of a stochastic process indexed by space.In other words, Suppose Z (s) is a random process observed at location s, the set Z (s) ≡ {z (s), s ∈ D } is a spatial data, where D , a subset of d , usually ( but not necessarily ) fixed and represents a spatial domain.According to Blangiardo [8], the spatial domain are distinguished as follows This process may sometimes be referred to as a discretization of D .
• Geostatistical or Point-reference data: z (s) is a realization at a specific location s in a continuous spatial domain D .
The location s is considered a coordinate made up of longitudes and latitudes, and sometimes includes altitudes.
The location s could also be represented in Cartesian coordinates.
• Spatial point pattern: The realization z (s) represents the occurrence or non-occurrence of an event at location s.In this case, the location itself is considered to be random.The random realization is a location indicator of the presence or absence of a phenomenon of interest in the domain D .In Agriculture, for example, the interest may be the distribution of a specific tree species, where each realization is the presence or absence of the tree specie in domain D .In epidemiology, the realization may be the house address of a patient having a particular disease [3,14].

| Spatial Priors
It is necessary to determine or specify a prior distribution for the posterior distribution's complete estimation in an empirical or full Bayesian approach.For example, in hierarchical models, the prior distribution assumed for a random field ( spatial component ) z (s) is termed, spatial model.We encountered several types of spatial models (priors) in the literature, and most were sub-class of the Gaussian Markov Random Field (GMRF) defined as a Gaussian random field with Markov property [14,86].
In the literature, due to the large class of priors, we collapsed the encountered priors based on the most frequently The non-GMRF is a large class that consists of non-trivial prior models that utilize a spatial correlation function to determine the covariance matrix of the spatial process in a continuous space, including the Asymmetric Laplace, Log-Gamma, skewed normal, Student-t process, and Dirichlet process.We created a class O t her s to accommodate unspecified models and those that do not belong to the above classes.
A response variable is a quantity used to describe a random process to relate it to a deterministic process mathematically.In statistical modeling, the most frequently used response variables are the discrete (categorical), ordinal, and continuous variables.The type of variables used in modeling a random phenomenon is intuitive from the process under study.The statistical models used to describe a random phenomenon vary depending on the quantity and parameters of interest.
The Bernoulli distribution is often used for modeling the random phenomenon of two possible outcomes.The Binomial, Negative Binomial, Hypergeometric, and Poisson distribution are frequently used for modeling count cases such as disease occurrence, wildlife, signal, and more.The Poisson distribution has been used to approximate Binomial distribution for a large sample size [98,11].The equality of mean to variance restriction imposed by the Poisson distribution considers the Negative Binomial a better choice to model a random variable that exhibits over-dispersion.
The Multinomial distribution is often used to model a phenomenon of more than two categories usually encountered in Biological experiments.It is a generalization of the Binomial distribution.
In the continuous case, a large class of distribution of the exponential family is used, such as Gaussian, Exponential, Student t, Weibull, Gamma, and more.However, according to the Central limit theorem, the Gaussian distribution is used to approximate both discrete and continuous distribution for large sample size [53].
An analyst's interest is to quantify the association of a random phenomenon and explanatory processes, describing a random phenomenon according to a set of explanatory variables.In the literature, the statistical models encountered are the Generalized Linear Mixed Model and the Hierarchical model, Survival model, and Spatial Econometrics models.
The details of each model are presented in the appendix.Additionally, we created a Proposed U nspeci f i ed and O t her s classes to accommodate proposed and validated models and unspecified models.The class of O t her s accommodates statistical models outside the above listed classes.
An appropriate prior distribution specification in a Bayesian inference remains a challenge in various fields of application.A prior distribution is associated with the representation of uncertainty of the interest parameters before data are observed.The elicitation of an appropriate prior distribution is a non-trivial task [78], and such challenges are accumulated in spatial models due to the large number of associated parameters involved.In our review, we came across four main approaches.
One way to set a prior distribution is to assume ignorance about the appropriate model, that is, (θ ∝ 1), and allows the data model to carry all the information.Such an approach is not always advantageous because inference on the parameters can be improved by performing prior elicitation based on identified characteristics or expert opinion.The elicitation procedure is termed elicited prior.
In elicited priors, convenient prior distributions are sometimes a choice and have spread across literature and have been set as default priors in most simulation packages.As a result, subsequent authors use such prior distribution verbatim.However, several authors did not explicitly state the type of prior used and were classified as not available.

| Estimation Method
In Bayesian inference, the prior information expressed through the prior distribution (θ) and the data likelihood (y |θ) are used for inferences.In simple or convenient cases, the posterior distribution given in (1) is easy to evaluate.
However, in practice, evaluation of the posterior distribution to make inferences is a non-trivial task because it usually contains compound integrands with complicated support that is not analytically integrable [27].In this sense, authors explore different approaches to make inferences.In the literature, we encountered several estimation methods and Also, the unspecified class was added to accommodate articles that neither discuss nor state the approach used in the estimation procedure.The Others class comprises of estimation methods that do not fit into the defined classes. (1)

| Simulation Study and Validation
A simulation study is a systematic and scientific computer procedure that involves fixing model parameters to generate data by pseudo-random sampling [73].It comprises two main steps; the data generation and the estimation.In the first step, a set of parameters is fixed and used to generate pseudo-random data.In the second step, the generated data is fed back to the model to estimate the "unknown" parameters and check for bias and model error.
A simulation study is usually carried out for proposed models and methods.The articles reviewed were classified into two; "Yes, " if the article contains statistical simulation studies, and "No" if it does not.
In addition to the simulation studies, we also investigated how Bayesian spatial models were validated using real data.It is a procedure to test for overfitting or underfitting.A model overfits if it performs well in the training set and badly in the test set, while it underfits if it performs poorly in the training set.A classical approach to cross-validation is to form a disjoint subset of a whole data into training and testing sets.The model is fitted on the former and tested on the later set.Doing this process k times until all observations in the data set participate in training and testing, once, it is called K-fold cross-validation.The whole data is split into k disjoint subsets, where the combined k − 1 sets serve as the training set and the remaining set of size n k = n/k serve as the test set n is the data size.A particular case to the k-fold is the Leave-One-Out cross-validation, where one observation serves as the test set, and the remainder n − 1 serves as the training set.After going through all subsets, the validation measures are statistically combined to make a valid conclusion.
Since the spatial models are often modeled in a Bayesian framework, we included the Posterior predictive check [32] class.In a predictive posterior check, a test statistic is chosen and computed for the observed data process.The same statistic is computed for replicated posterior predictions of the process.The model is said to present a good fit if the posterior prediction average is close to the test statistic for the observed data [27].We added None or not applicable class to accommodate articles that did not conduct cross-validation and Others to accommodate the validation method not mentioned above.

A N A LY S E S
As described in the search procedure section, a total of 552 articles were selected after filtering with the exclusion criteria (duplication and context).After a careful watch, the articles were categorized into the first class, labeled as being only an applied paper, theoretical, or both, where 4 (1%) of the papers showed no application (only theoretical with synthetic data), 188(34%) showed an improvement in the field with real-world application and 360(65%) only applied the existing methodologies.
The papers were sub-divided into five classes of application field: Agricultural and Environmental Science, Economics and Humanities, Medical Science, Physical Science and Engineering and Other.Mainly, three fields hold the majority of the publications, which are Agricultural and Environmental Science (30.1%),Economics and Humanities (30.6%), and Medical Science, which includes epidemiology, (33.7%).Moreover, the spatial domain used was also taken into account.The Area/lattice, Geostatistical, and Spatial Point Patterns.The Area/lattice occurred 65.6% and Geostatistical 31.2% of the F I G U R E 2 Distribution of the response variable type.reviewed articles, and in combination, they hold 95.8% of the publications.It is important to note that more than one spatial domain could be used in an article, such as the 1% observed in this review.This procedure is common when a continuous spatial domain is discretized to lower computational burden.
The core part of this review is the spatial priors (models) used in the literature.While the literature contains numerous spatial priors applied in various problems across various research fields, this systematic review only presented models encountered during the review.Gelfand et al. [29] explicit the first incorporation of a spatial prior, and it is gain.
Table 1 summarized the spatial models encountered and were classified into distinguishable groups.Among the spatial models, the CAR family, usually used for area/lattice spatial domain, Stochastic Partial Differential Equations (SPDE), Gaussian Markov Random Field (GMRF) except for the CAR family, models with author-specific defined covariance structure often for a continuous spatial domain, and Non-Parametric methods.The CAR family appeared in 44.2% of the published articles reviewed.Of this percentage, the CAR and the BYM model appeared 96.3%.The SPDE appeared about 3.9%, and the non-parametric procedure appeared 1.3% of the articles reviewed.The GMRF, exempting the CAR family and the SPDE, appeared in the literature 4.0%.Consequently, to the diverse application of the spatial model, authors specified the type of covariance structure based on the prior knowledge of the interactions between the phenomenon of interest that appeared 31.3% of the articles reviewed.Only a single documentation of the application of spatial models on robotic technology was encountered [95].Other models that could not fit into any of these groups appeared 7.2%, and 9.3% of the articles did not described or state the type of model adopted.
The observation or response variable's nature dictates the statistical model class to be adopted to make inferences.
In our search, as shown in Figure 2, the discrete (Countable) was the most used, then the continuous and binary response variable types.Additionally, the most adopted statistical model to spatial analysis, withing a Bayesian framework, is the GLMM.Owing to the integral complexity of the posterior marginal distribution, the MCMC estimation method is the most frequently adopted numerical integration, Figure (3).
Maybe the most critical question regarding Spatial Models is related to specifying Spatial Prior.Elucidating events in an area, some times even under a few frequencies, is desirable through direct probabilistic statements that may unravel hidden patterns [47].Whereas the inter-dependence may conduct spatial correlation with Bayesian Models F I G U R E 3 Class of models often used the spatial modelling and its numerical estimation method distribution.Given the class of GLMM/Hierarchical Models, Markov Chain Monte Carlo (MCMC) is the most used computation intensive technique.
across the spatial fields, and the Bayesian approach enables the information expertise to be allocated with the acquired data.The results obtained in this systematic review, as Table 2 shows that the knowledge of the expert is used (30.43%),although it can be better explored.

Prior specified
Freq.

Elicitated from experts or from the problem 168
No explicit use or reference/not applicable 101 Used verbatim from the literature 166 Vague prior (Non-informative) 119 The Authors who were more recurrent in the literature in this past 20 years were James Law, ACA Clements, and Hei Huang.F I G U R E 5 Publication Growth, regarding the Bayesian Spatial Models, across time.

C O N C L U D I N G R E M A R K S
Spatial statistics has gained tremendous attention in recent years owing to the efficiency in collecting spatial dependence data.Neglecting such dependencies may result in bias and consequently leads to wrong inference.The Bayesian approach to analyzing spatial data often outperforms the frequentist approach, given that the prior information is taking into account.In a Bayesian framework, spatial priors play a significant roll in accounting for space dependencies.The consistency in the improvement of data collection and computational tools in analyzing spatial data, Bayesian spatial statistics will further penetrate numerous fields and becomes one of the leading tools for analyzing data.
Many authors account for spatial dependencies assuming a Gaussian random field.In many real data applications, the Gaussian random field may be inappropriate, especially in extreme data, skewed data, data with spikes and heavy tails.
Examining a different random field such as the Laplace, Student-t, Pareto, Nakagami-m, and more may improve inference.
A significant complication of assuming these distributions is the non-trivial method of fixing a prior distribution for the model parameters.For instance, the prior distribution for the degrees of freedom of the student-t distribution is not a trivial task [76].Regardless of the prior distribution, eliciting priors for the parameters is critical, and when wrongly assumed, it could lead to misleading results and inference.To circumvent these, which is also not a trivial task, it is essential to consider objective priors for the random field parameters and hyper-parameters to improve inference.
The Bayesian spatial literature lacks sufficient information on the objective priors, such as Jeffery's prior, reference priors, matching priors, and more.These priors stand out to elicit experts ideas that could improve the inferences of the spatial dependencies.To derive an objective prior distribution for a spatial random field's smoothing parameter is currently an open problem that needs urgent attention.
The Markov property is a useful tool to lower the computational cost in Bayesian inferential statistics by subjecting the immediate neighbors' spatial dependencies.However, the realization of some random phenomenon exhibits strong spatial dependencies beyond the immediate neighbors, and truncating such dependencies structure will result in bias and incorrect inferences.Moreover, there is an insufficient standard approach to determine the covariance matrix structure of the spatial effects.Thus, spatial smoothing is not trivial to compare across different models.
In Neuroscience, the application of Bayesian spatial statistics to brain experiments is gaining interest [19,46,91,97].
The complexity of the brain structure has prevented the application of classical spatial models.In other words, the primary assumption of spatial contiguity in the analysis may result in incorrect inference owing to brain complexity.
Moreover, a response received at one location on the skull may be due to brain activity in the opposite location.Beyond complexity, the dynamics of the body system of the subject influences the experiments.Thus, accounting for such complex structures is an open problem that requires further studies.
Despite the increasing availability of the Big IoT data, the spatial model has not received adequate attention as a classification algorithm and regression machine learning models.Application of machine learning to spatial data could unravel hidden patterns, and applied to efficiently navigates a robotic technology [95] through space, which creates new research lines.

C O N FL I C T O F I N T E R E S T
There is no potential conflict involved in this report

R E F E R E N C E S
[1] Charles E. Antoniak.Mixtures of Dirichlet processes with applications to Bayesian nonparametric problems.The Annals of Statistics, 2:1152-1174, 1974.

A P P E N D I X A
The coding of the characteristics was framed according to the conceptual classification scheme developed by Hachicha and Ghorbel [37] and applied by Tiago et al .[27].In this framework, bias is limited in survey data and clarifies reporting results and findings to draw concrete conclusions.Such classification is useful to researchers as it provides an overview of the research methodology's application and can reveal research gaps in the literature.
This review employed a conceptual classification scheme of 10 items presented in Table 3.For lucidity, each item in the classification scheme is discussed in detail.
List of questions for conceptual classification scheme 1.0 Is it only an application?

A P P E N D I X B
In this Appendix, we made an overview of the primary spatial statistical models explored in this systematic review.

| Class of Gaussian Markov Random Field (GMRF)
In standard form, the covariance matrix of a Gaussian random field is positive definite and often dense.This restriction makes it difficult to construct an appropriate covariance matrix.Furthermore, the computational cost is high due to the cost of O (n 3 ) to factorize the covariance matrix [63].To overcome the above challenges, it is necessary to construct a sparse covariance matrix of less computational cost and closely approximate the dense covariance matrix.
A possible choice is to adopt a Markovian property on the covariance matrix to have a GMRF approximation.It is a discreetly indexed finite random vector that has a joint Gaussian distribution with a sparse precision matrix.Following the definition given by Rue [85], the set Z (s) ∈ d over a spatial domain s ∈ D is said to be a GMRF with respect to graph (V , D ), centered at µ and precision matrix Q > 0 and density of the form where Q i ,j 0 if i and j are considered neighbors, V is a set of vertices, and ξ is the edges of graph G representing D .Usually, the covariance function is constructed such that it is a function of the relative positions of location i and j in D .Q is considered stationary when it is only a function of the distance between the locations in D .Let δ i = {z : z is a neighbor of i }, the full conditionals is given as for every z i ∈ Z, and Z −i is the vector Z excluding z i .Some of the subclasses of GMRF identified in the literature for spatial smoothing are the proper and Intrinsic Conditional Autoregressive (ICAR), Beseg York Mollie (BYM and BYM2), Lourex CAR (Leroux CAR), and Dean's CAR.

-Conditional Autoregressive (CAR)
Let the spatial domain D be partitioned into n disjointed areas E such that D = ∪ n i =1 E i .E could be regular or not, as in Lattice or area respectively.Let z i be a random variable observed at E i and the vector Z = (z i , ..., z n ) T with µ = (µ 1 , ..., µ n ).
Let Z −i be an (n − 1) dimensional vector such that Z −i = (z 1 , ..., z i −1 , z i +1 , ..., z n ) T , the full conditionals of a proper CAR model is given as where .., n and c i j is a function of the adjacency between E i and E j .To make π(•) a proper distribution, C = (c i j ) and κ are carefully chosen.A frequent choice for C is a function of contiguity between the areas.Let A = (a i j ) be an n × n matrix, such that a i i = 0, a i j = 1 if z j ∈ δ i for i j , and 0 otherwise.Noticed that A is symmetric, since if E j is a neighbor of E i , then E i is a neighbor of E j .We define c i j = ρ σ 2 and d i = j a i j [5,28,85].The full conditionals and the joint distribution of ( 4) is given by where ). and |ρ | < 1 is chosen such that ( − C) −1 > 0. Beseg, York and Mollie [6], proposed the intrinsic CAR (ICAR), by setting ρ = 1, then c i j = a i j d i . The full conditionals of the ICAR The model inferred that the conditional expectation of z i equals the weighted deviations of its neighbors in addition to its mean.

-Beseg York Mollie (BYM)
The BYM, proposed by Beseg [6], is a variant of the CAR model that incorporates an additional term to control for overdispersion in spatial data.Suppose z is partitioned such that z i = u i + v i , the unstructured term v is modeled as and the structured component is modeled as The BYM poses identifiable problem such that each observation is represented by u i and v i , thus, only the sum u i + v i is identifiable [55].Setting appropriate hyperparameters is challenging.However, constraining φ to sum to zero allows the confounding problem to be avoided and both components to be fitted.

-Simpson Conditional Autoregressive
For the purpose of scaling and interpretablility of the hyper prior, Simpson et al [90] proposed a modification of the BYM which avoid the problems possed by BYM model.The combined random effect and the covariance matrix are given by, where σ 0, φ ∈ [0, 1] and Q is the ICAR covariance matrix such that V ar (u i ) ≈ 1.

-Leroux Conditional Autoregressive (Lourex CAR)
Leroux et al .[60] proposed a variant of the CAR model which, unlike the BYM, only requires a single set of spatial component.The full conditionals is given by which is a limiting case of the I C AR model.The model avoids the identifiability problem in the BYM model by the specification of a single hyper parameter for random effect [13].The specification (11) smooths the neighboring random effect weighted by ρ and the global mean weighted by 1 − ρ.

-Conditional Autoregressive dissimilarity
The variations of the CAR model earlier discussed exhibits a global degree of spatial smoothing.In many instances, a global smoothing may be inappropriate, especially areas that exhibit locally constrained spatial structure.Lee and Mitchell [57] proposed a CAR dissimilarity to smooth area elevated risks which, depends on local spatial parameters.
It is one among many proposed local spatial smoothings, such as locally adaptive model [58], localized conditional autoregressive [59], to list a few.

| Class of Gaussian Non-Markov Random Field Models
Markov property is known to lighten the burden of factorizing a dense covariance matrix.However, a random realization with a strong correlation structure outside its neighborhoods can be poorly accounted for.Also, the discretization imposed by the GMRF sometimes does not account for the presence of discontinuities in spatial domains.These could be addressed by relaxing the Markov assumption and possibly allowing a different distribution other than Gaussian.

-Spatial Weight Matrix
In spatial Econometrics, a weight matrix and a correlation parameter are often used to specify a dependency structure between economic variables of interest observed at different spatial locations in a Spatial lag model (SLM) and Spatial Error Model (SER) [36].The weight matrix is essential in the covariance matrix specification for parametric models, such as the class of GMRF prior.The spatial weight configuration specification varies with the relationships between geographical locations, such as spatial distance, interactions, nearest neighbors, and contiguity.The specification is categorized into four: adjacency-based weights, weight-based on geographical distance, the distance between covariate values, and hybrid of geographical distance and covariates [21].In most cases, nearby neighbors to a spatial reference location receive higher weights compared to farther neighbors.Though the diagonal elements of the spatial weight matrix are universally accepted to be set to zero [35], in the literature, there is no standard to define the weight matrix the specification is dominated by choice of computational convenience [36].

-Spatial Covariance Function
In this section, we briefly discussed the frequently used stationary and isotopic correlation functions to construct a valid correlation matrix for a random field in a Geostatistical models.A covariance function C (s i , s j ) is said to be stationary when C (s i , s j ) = C (s i + l , s j + l ) for any lag l , and isotopic when it only depends on the distance s i and s j .A covariance function must be positive definite and symmetric.
Consider a Gaussian random field Z (s), such that the realizations {z (s i ), s i ∈ D, i = 1, ..., n }, D ⊆ 2 are of interest.
Note that the information we seek can be described by the mean ¾(Z (s)) and the covariance matrix C ov (s i , s j ).Thus, we could describe the mean by a linear function of covariates, and the covariance matrix as C ov where ρ( | |s i − |s j |) is the correlation function.
Let d (s i , s j ) be Euclidean distance from site i to j .The four most commonly used correlation function [14,86] are given by, where K v is a modified Bassel function of the second kind of order v , S v is the scaling factor, and r > 0 defines a significant range of correlation for exponential, Gaussian, and Matern correlation function.While in spherical correlation function, r is defined as the correlation length [86], and v is the smoothness parameter, which measures the differentiability of the Gaussian random field [94].
-Skewed Gaussian random field Skewed form of the Gaussian random field has been used to model a skewed air pollution data.A skewed Gaussian random field cannot be defined in the same way as a Gaussian random field due to its marginals' dependence on its component parameter.For condense information, see [48,49,83].

-Geostatistical model
The geostatistical model proposed by Clement et al .incorporates an exponential distance decay on the average of a Gaussian process in a geostatistical model given by where d i j is the distance between points i and j , φ controls the decay rate and k is the smoothing parameter.The decay rate is modeled as uni f or m(0.1, 6).However, the parameter could be intuitively determined from the random process of interest [13].Moreover, other valid spatial covariance function could replace the exponential function.

-Stochastic Partial differential Equation (SPDE)
SPDE, proposed by Lindgren [64], allows fitting a GRF with a continuously and smoothly decaying covariance function while gaining computation benefits from a GMRF representation.It represents a continuous random field using a discretely indexed spatial process.According to Blangiardo and Cemeletti [8], the SPDE model is defined as follows.Let z(s) be as defined previously for a continuous spatial points s ∈ D ⊂ R d .The SPDE model is given by is the Laplacian, α controls the smoothness, τ controls the variance.(s) assumes a Gaussian spatial white noise process.However, (s) can also assume a Laplace noise, especially, when the data exhibit spikes and heavy tail.The solution to the above differential equation is a stationary Gaussian random field z(s) with Matern covariance structure σp(s i , s j ), where p(s i , s j ) is a Matern correlation fuction between site i and j as defined previously, where k = Sv r , S v , r and v are as previously defined.The solution to the SPDE model is approximated by a basis function representation defined on a triangulation of the domain D .That is z (s) = G g =1 ϕ g (s) xg , G is the total number of vertices of the triangulation, {ϕ } is the set of basis function, and { xg } are zero-mean Gaussian distributed weights.

| Class of Non-Gaussian random Fields Models
A non-Gaussian process has increasingly been useful for modeling extreme ill-behaved random processes, such as predicting an earthquake and atmospheric temperature.This section discussed the asymmetric Laplace process, inverse Wishart distribution, Pareto process, and Dirichlet process.

-Asymmetric Laplace Process
Suppose Y (s) be an Asymmetric Laplace random vector with parameter p and τ, AL(p, 0, τ).According to Kuzobowski and Pogorski [51] and Fontanella et al [26] Y (s) can be written as a sum of normal and exponential process, given by, where ρ z (s, s * ; θ) is a valid spatial covariance function.Z (s) is a Gaussian random field to accounts for spatial errors and it exist as a standard normal in its marginal form.ξ(s) is marginally exponential distribution with rate τ, and it is conditionally independent of Z (s).Thus, the conditional distribution of Y p (s) given ξ(s) is given by, Kristian and Alan [66] discussed approaches to model ξ(s) in a generalized quantile regression.In the spatial case, they defined ξ(s) through CDF or copula transformation by letting ξ(s , where V ξ (s) is again a Gaussian process with a valid spatial covariance, Φ is a standard normal CDF, and F τ is an exponential distribution CDF.
According to Bradley et al .[10], a n dimensional vector q = µ + V γ is a multivariate log-Gamma random variable and its distribution, denoted by M LG (µ, V, κ, α), mean and variance are given by, where |A | is a determinant of a matrix A.
Consider a spatial random process Z (s) distributed as a multivariate log-Gamma, Yang et al .[41], and Hu and Bradley ) is a function of location s and s * , θ is a vector of some parameters of interest, and is an identity matrix of appropriate dimension.An appropriate prior distribution are assigned to α and κ.
-Student-t Process Some random processes exhibit heavy tail property, which may be inappropriate for using a Gaussian distribution, hence applying Student-t distribution in modeling spatial random process.Suppose Z (s), as defined previously, be a random process that exhibits heavy tail property, observed at location s ∈ D. This can be represented by setting Z (s) to be distributed as a multivariate Student-t distribution Z(s) ∼ t n (0, Σ, v), where Σ is a covariance matrix which is determined by a valid covariance function.The i , j t h element of Σ is given by C ov (Z (s i ), Z (s j )), and v is the degree of freedom.The Σ is used to account for spatial dependencies.Moreover, the covariance matrix can further account for spatially structured and unstructured (nugget) effects [88].The main challenge of the Student-t process is the difficulty of assigning appropriate prior specification on v to make inferences in a Bayesian analysis.Fonseca [25] derived a Jeffreys-rule before v, which lead to a proper posterior distribution.Moreover, Ordoñez [77] extended it to a spatial domain and derived a reference prior to the joint spatial hyper-parameter and the degrees of freedom v .

-Spatial Mixture model
Green and Richardson [34] broaden the application of hidden Markov models for the random component z i by extending it to a spatial domain.It utilizes the model benefits of the Hidden Pott model.The allocation of the mixture components to clusters utilizes a spatial dependence structure, and the numbers of clusters are considered to be random.According to Best and Thomson [7] the model specification is given in (17).Areas estimated to belong in the same clusters need not be contiguous, where ψ > 0, c max is the upper bound of the number of clusters, U (w) = i ∼i [z i = z i ] is the number of same labels pairs to i in a neighbouring area i i , parameter η j is associated with each component, and δ w (ψ) = l og ( w ∈{1,2,...,k } n exp ψU (γ) ) is the normalizing constant, where the sum is the total possible ways for the allocation for the n areas.The estimation of the numbers of clusters or components are obtained through reversible jump Markov chain Monte Carlo.Notice that p(w |ψ, k ) represents the probability function of the Pott model [62].
be the realizations with replicates of a spatial random field at distinct locations s.Gelfand et al .[31] described the spatial Dirichlet process as follows.Let Z t = (Z t (s 1 ), ..., Z t (s n )) , where t represents replicates at site s i .A Dirichlet process is a random probability measure defined on a measurable space (Ω, »), denoted by D P (v G 0 ), where v > 0 is the precision parameter and G o is a specific base probability distribution over (Ω, »).Let k i , i = 1, 2, ... be independent and identically distributed Bet a(1, v ).The resulting random process for Z (s) from Dirichlet process D P (v G 0 ) defined on (Ω, ») can be written as ∞ i =1 λ i δ θ i ,D , where δ k represent a point mass at k and ), i = 2, 3, ... , and {θ i (s) : s ∈ D } are realizations from base probability G 0 which can possibly be a stationary Gaussian process.Notice that the DP process are independent, however, MacEachern [67] relaxed this condition by deriving the dependent D P .Moreover, Gelfand et al .[31] extended the dependent D P to account for spatial dependence.The property that DP process is almost surely discrete restricts its applications to a wide class of continuous problems.Hence, Antoniak [1] derived a mixtures of DP to circumvent the problems, and thus extended it to handle continuous cases.

| Statistical Models -Generalized Linear Mixed Model
A classical linear model formulation is given by and the matrix form is given by, where the latent field β = (β 0 , β 1 , ..., β p ) defines the relationship between response variable Y and covariate X. [20].Each outcome y i is assumed to be generated according to a Gaussian distribution.The mean depends on related covariates through ¾(y |µ).A generalized form of (20) relaxes the assumption that the errors are Gaussian distributed and each outcome is generated from a non-Gaussian distribution such as Binomial, Poisson, Beta, among others [22].It allows these random processes to be modeled through a link function and allows the magnitude of the measurement error to be a function of the predicted estimates.That is, where g is an appropriate link function such as the l og i t for a Binomial model and l og e for the Poisson model.The mixed form of Equation ( 21) incorporates a function f ( ) to relax the linearity assumption on covariates or to introduce a random effect usually modeled as a random walk, auto-regressive, or penalized spline models [22].The mixed form is given in (22) y i ∼ (. |θ) where z k is the k t h random effect and v i is a spatial effect assigned spatial prior distribution.The latent field of interest is given as Θ = {β , f , v }.Equation 22 is termed the Generalized Linear Mixed Model (GLMM).In a Bayesian setting, all the parameters in the model are considered random, and appropriate prior distributions are assigned.In the absence of prior information, a non-informative prior is assigned to β ∼ N or mal (0, 10 6 ) and l og (1/σ 2 ) ∼ l og g amma(1, 10 −5 ) [8].
A Bayesian hierarchical model contains the data, prior, and hyperprior stages, where θ i ∈ Θ |y is of main interest.These models are usually referred to as latent Gaussian models, which are flexible and can accommodate a wide range of models.

-Survival Model
A survival data analysis models the time to the event occurrence.This can be time to death of subject under study, process failure time, or time to radioactive emission.Let T i be a random variable of survival times, then S (t |τ) = (T > t |τ) is the survival function that, for example, determines the probability that a patient survives over time t .It is assumed that all subject are alive, S (0|τ) = 1, and all subject will eventually die lim t →+∞ S (t |τ) = 0.The survival function is expressed as a distribution function F (t ) = (T < t |τ) = 1 − S (t |τ) with probability distribution f (t |τ).The hazard function h(t ) measures the probability that an event will occur at a small instance ∆t after the subject has survived through time t .That is, where H (t ) is the cumulative harzard function.
Let (t i , δ i , x i ) represents a survival data, where T i represents survival times of the i t h subject of interest, x i is a predictor variables, and δ i = (T i ≤ C ), where C is a censoring threshold set a priori.Considering the covariates, the likelihood of the model is expressed by, where θ = {β 1 , ..., β j , τ, Φ }, β i , is a fixed effect, z k is a random effect, v i with parameter Φ, is a spatial random effect assigned a spatial prior in a Bayesian framework [15,75].f (t i |τ) can assume Exponential, Weibull, Logistics, or lognormal distribution to list a few.The censoring status δ controls the contribution to the likelihood of subjects that experienced the event and those that survived through the entire study period.Kaplan-Meier provides non-parametric estimates of the survival curves.The proportional hazards model and accelerated failure time model are the most frequently used frequentist parametric methods to analyze survival data [53].

-Bayesian Spatial Econometrics
Spatial analytical tools are widely used in Economics to quantify the spatial dependencies and heterogeneity of economic variables.Referred to as spatial econometrics, it extends the traditional econometrics to a spatial domain.The most frequently used models are the Spatial Autoregressive (SAR), Spatial Error Model (SEM), and Spatial Durbin Model (SDM).The models are briefly described as follows.
Let Y be a response variable assuming a linear relationship with explanatory variables X, the SAR model is represented by, W be a spatial weight matrix, ρ is a spatial autocorrelation parameter, and β is a vector of the regression parameters.
The SAR model assumes that the dependent variable is spatially autocorrelated [44].In contrast, the SEM model assumes the error is a spatial correlation.The SEM formulation is given by, The SDM is an extension of the SAR model which assumes that the dependent variable and covariates are spatially correlated.The formulation is given by, According to LeSage and Pace [61], SDM is appropriate when the included covariates are correlated with spatially correlated variables not included in the model.In the next sections, we present each class's main idea of the spatial literature's estimation methods.

| Estimation Method -Maximum Likelihood Estimation
The maximum likelihood estimation approach involves an optimization problem to determine the best sets of distribution parameters representing data.In spatial statistics, the MLE method is usually used to estimate global spatial dependencies in a spatial econometric model.Authors often compare its estimates with one obtained in a Bayesian inferential framework and highlight the MLE approach's inadequacies to spatial statistical inferences [74,79,80,103].
Let y i ∼ f (. |θ i ), where θ i ∈ Θ is a parameter of interest and Θ is the parameter space.y i need not be independent and identical.The estimation procedure involves computing the likelihood (log-likelihood) of the data distribution and determine the sets of parameter θ * = {θ * i : θ * i ∈ Θ, i = 1, 2, ..., n } where the likelihood is maximum.For independent y i , i = 1, 2, ..., n, the likelihood is given by, An equivalent approach to the classical MLE when the dimension of the model parameters is large or complex for estimation is the quasi-Maximum likelihood estimation (QMLE).Unlike the MLE, the QMLE maximizes a function l og L * (θ |y) that is related to the logarithm of the likelihood function L(θ |y).Su and Yang [93] proposed a QMLE of dynamic panel models with spatial errors and was further broadened by Yu, De Jong, and Lee et .al[103].An equivalent approach to MLE in the Bayesian setting is the Maximum a Posteriori (MAP).It involves maximizing the posterior conditional probability, where p(θ) is the prior distribution and p(y|θ) is the data likelihood defined as L(y|θ) in ( 28).

-Expectation Maximization (EM Algorithm)
The maximum likelihood estimation limitation is the assumption that the variables that generate the process are all observable.In practice, this assumption rarely holds.One possible choice to overcome the limitation is the estimation through the EM algorithm.
The EM algorithm, proposed by Dempster et al .[17], fits a model to a latent representation of a data rather than merely fitting distribution models.It can work well in data that contains unobserved (latent) variables.The algorithm, an iterative method, has two major stages: estimating the latent variables (E-step) and maximizing the model parameters given the data and the estimated variables (M-step).
Let θ be initialized model parameter, the E-step is used to update the latent space variables z , usually discrete or cluster in particular, through p(z |y, θ), where y is the observed data.To update θ, the expectation ¾ z|y,θ * l og (p(y|z, θ)) The iteration continues until the difference between the current and the previous expectation is lesser than > 0 set at the initial stage.

-Markov Chain Monte Carlo (MCMC)
Given a posterior distribution p(θ, ψ |y) ∝ p(y |θ, ψ) × p(θ, ψ), (32) and assuming that the posterior p(θ |y , Ψ) is of known form, such as a standard probability distribution, we can resort to Monte Carlo approach to approximate posterior quantities h(θ), which could be the mean, median or higher moments.The procedure consist of simulating an m random samples from p(θ |y ), say {θ 1 , θ 2 , ..., θ m } and evaluate the unknown quantity h(θ) using the empirical average ¾( ĥ(θ) |y) Under the Law of Large Numbers, the empirical distribution will converge to the true distribution.In the case of a joint posterior distribution, an approximation of the posterior marginals are achieved by first sampling from a marginal distribution of a subset of parameters given its complements and then use it to evaluate the full joint distribution.
In practice, the posterior distribution's functional form is unknown or complex, and independent samples are not feasible.An alternative approach comprises generating independent samples from an importance distribution (θ |y ) which is a close distribution to (θ |y ).Empirically, the quantity h(θ) is obtained as This approach is not trivial for a large number of dimensions of θ [39].A more widely used alternative approach comprises generating correlated samples by running a Markov chain whose stationary distribution converges to the posterior density.The posterior summaries are computed from these samples using the empirical method, as described above.Suppose χ is the state space of the posterior distribution.As stated by Blangiardo [8], the convergence of the Markov chains stationary distribution to the posterior distribution requires that the Markov chains are irreducible (the chain has a positive probability of reaching all region of χ regardless of the starting point), recurrency (the limit of the probability of the chain visiting a subset χ infinitely many times is 1), and aperiodic (the chain do not circles when exploring χ).The highlighted procedure is referred to as MCMC.Gibbs sampler and Metropolis-Hastings algorithm are the most frequently used standard MCMC algorithm in Bayesian inference literature.For a description of these algorithms, see [8] (pp.[91][92][93][94][95][96][97][98][99][100][101][102][103].

-Integrated Nexted Laplace Approximation (INLA)
INLA, proposed by Rue et al [87], is an alternative approach to the estimation of posterior marginals.It has gained considerable attention and has been proven to outperform the MCMC approach in computational speed [87] .The availability of the R − I N LA simplifies the implementation of the approach which authors from the diverse field have The objective is to obtain the posterior marginals p(θ i |y ) for each parameter in the vector and the estimates of the hyperparameters given by, p(ψ k |y) The INLA approach utilizes the model assumptions to approximate the marginal posterior distribution and its moments based on Laplace approximation [99].According to [8,9] INLA approximation follows the following steps.
Firstly, the posterior marginals of the hyperparameters are approximated, that is, is then used to obtain p(θ i |ψ, y ).
used and seldom-used classes, which are the Conditional Autoregressive (CAR), Beseg York Mollie (BYM), Lorex CAR, Stochastic Partial Differential Equation (SPDE), GMRF (none of the above), Non-GMRF and Others.See table (1) for summary.Details of each model are presented in the Appendix.The Lorex CAR class comprises priors with similar specifications, such as Dean's and Simpson CAR models.The GMRF class consists of GMRF priors except for the CAR family and the SPDE earlier stated.
classify them into the Markov Chain Monte Carlo (MCMC), Integrated Nested Laplace Approximation (INLA), Expectation-Maximization (EM) and Maximum (Penalized quasi) Likelihood Method classes.Details of these estimation methods are presented in the appendix.The MCMC class comprises of all numerical approximation that utilizes Monte Carlo method.
Figures 4 displays the visualization towards the relevant authors in the Bayesian Spatial Models publication field.The top 5 journals which most contained articles, from the obtained data set, which combined theoretical methodology with real-world application publishing on Spatial and Spatio-temporal Epidemiology (#15), Accident Analysis and Prevention (#14), PLoS ONE (#14), Spatial Statistics (#11), and Environmentrics (#10).Whereas, across time, the spatial modeling publication rate using the Bayesian approach proliferates, as shown in Figure 5.The year 2020 refers only to the first half of the year.F I G U R E 4 Most frequent authors on the Bayesian Spatial Models.Left-hand graph is a Tag Cloud for the 50 most frequent authors in the past 20 years.Right-hand graph is a Barplot displaying the Top 10 authors and its relative frequencies.

•
Area or Lattice data: it is a simple way to represent spatial data in the domain D .In this type of spatial domain, z (s) is a random aggregated realization across an area s of distinct boundaries.For area data, the boundaries are irregular, such as administrative divisions, while for lattice, the boundaries are a regular division of D .For simplicity purposes, it may be necessary to aggregate other types of Spatial domain realizations to form area or lattice data.