Bayesian Inference for Multivariate Spatial Models with INLA

Bayesian methods and software for spatial data analysis are well-established now in the broader scientific community generally and in the spatial data analysis community specifically. Despite the wide application of spatial models, the analysis of multivariate spatial data using the integrated nested Laplace approximation through its R package (R-INLA) has not been widely described in the existing literature. Therefore, the main objective of this article is to demonstrate that R-INLA is a convenient toolbox to analyse different types of multivariate spatial datasets. This will be illustrated by analysing three datasets which are publicly available. Furthermore, the details and the R code of these analyses are provided to exemplify how to fit models to multivariate spatial datasets with R-INLA.


Introduction
Multivariate spatial models have been studied by several authors.For example, Van Lieshout and Baddeley (1999) describe dependence between multivariate point patterns by proposing novel summary statistics.Additionally, multivariate log-Gaussian Cox processes have been used to analyse multivariate point patterns (Diggle et al., 2013;Waagepetersen et al., 2016;Gómez-Rubio et al., 2015).Furthermore, several studies have been published for analysing multivariate lattice data.For example, in MacNab (2018) an insight of the generalization of the univariate models to multivariate models is extensively discussed and in Martínez-Beneito et al. (2017) a framework to analyse multiple variables is proposed in the context of disease mapping.A review of the multivariate spatial models in disease mapping is performed in Martínez-Beneito and Botella-Rocamora (2019).These studies use Markov Chain Monte Carlo (MCMC, Gilks et al., 1995) methods to exemplify their proposals.In the context of spatial modeling, these methods can lead to a high computation cost.
In recent times, the integrated nested Laplace approximation (INLA, Rue et al., 2009) has provided an alternative way of fitting Bayesian hierarchical models.Furthermore, when estimating continuous processes, the stochastic partial differential equation approach (SPDE, Lindgren et al., 2011) combined with INLA can be employed.The different spatial models which can be fit with R-INLA (Rue et al., 2020) have been compiled by several authors (Lindgren et al., 2015;Blangiardo and Cameletti, 2015;Bakka et al., 2018).For a recent review, the reader is referred to Chapter 7 of Gómez-Rubio (2020).Krainski et al. (2019) provide an exhaustive tutorial about how to adjust advanced spatial models in R-INLA.Additionally, INLA can be combined with other algorithms such as MCMC techniques in order to adjust models which can not be adjusted solely with INLA (Gómez-Rubio and Palmí-Perales, 2019).
The aim of this paper is to provide a summary of how to analyze the different types of multivariate spatial data with R-INLA.The increasing interest (and application) of multivariate models have been the main motivation of this work.Therefore, the main objective of this article is to describe how to analyse any multivariate spatial dataset using R-INLA, which appears as an interesting toolbox to analyse this models' class in a Bayesian framework.
The remainder of the manuscript is structured as follows.First, a brief description of the different multivariate models is detailed in the case of areal, geostatistical and point pattern data.The next sections discuss the prior distribution choice and the structure the data.Then, three examples about how to analyse multivariate spatial models with R-INLA are detailed illustrating the main objectives.Finally, a brief summary of the conclusions of this work appear in the last part of this manuscript.

Multivariate Spatial Models
The R programming language offers a wide range of standalone packages for analysing spatial datasets.Several of them focus on a particular type of spatial data.For instance, point pattern analysis can be performed with spatstat (Baddeley et al., 2015) and spatialkernel (Gómez-Rubio et al., 2017).Geostatistical data can be modeled using gstat (Pebesma and Wesseling, 1998;Pebesma, 2004) or spBayes (Finley et al., 2007(Finley et al., , 2015)).Other R packages such as CARBayes (Lee, 2013)

Areal data
When analysing lattice data, the domain is divided into non-overlapping areas in which the data are collected.It is usually considered that two areas are neighbours if they share a common boundary.This adjacency structure is often included to account for spatial autocorrelation (Banerjee et al., 2014).
When the value of several variables are recorded in each area the resulting data become a multivariate lattice dataset.The joint analysis of the spatial distribution of several variables allows to detect similar (spatial) patterns between some of these variables (Banerjee et al., 2014, Chapter 10) while estimating the spatial effects.We will illustrate the analysis of this type of data using a Poisson regression model commonly employed in spatial epidemiology to analyse count data but other similar models can be proposed for binary or continuous outcomes.
Given the d-th variable of interest (with d = 1, . . ., D) and area i (with i = 1, . . ., n), the response of interest Y d,i can be modeled using a Poisson distribution with mean µ d,i : The mean is usually modeled as a sum of different terms through a link function.The selection of these terms depends on the available data and the proposed model structure.For instance, one option can be described as Here, ψ(•) is a link function (which in this case is the logarithm function), α d is a specific intercept, u i a shared (between all or a group of variables) spatial term and v d,i a variable-specific random effect.Note that restrictions may be required on v d,i to make all effects identifiable (Rue and Held, 2005).
In the context of disease mapping, the usual variables of interest are the counts of mortality or incidence of different diseases over the study region.Now d represents the specific disease, therefore, following the above structure, the observed number of cases of the dth disease in the i area, Y d,i , can be modeled as where, E d,i and θ d,i are the expected number of cases and the relative risk of disease d in area i, respectively.As before, α d is a disease specific intercept (to account for differences in the total number of observed cases), u i a shared spatial term (which does not depend on the disease) and v d,i a diseasespecific spatial random effect.
Several authors (see, for example, Martínez-Beneito, 2013, and the references therein) have proposed different approaches to model multiple diseases in space and time.Gómez-Rubio et al. (2019) propose a separable spatio-temporal model with weighted shared components that can be used to detect diseases with similar patterns.In Palmí-Perales et al. (2021), the authors have developed a R package (INLAMSM) which builds on top on R-INLA to adjust some of the most common multivariate lattice approaches.

Geostatistics
Geostatistical datasets are built with values of variables which vary continuously over a spatial domain.These datasets contain observations which are geographically referenced, i. e. both the value and where it is collected (the coordinates) appear in the dataset.Then, the estimate of the variation surface of these variables is computed using geostatistical models.
The R Journal Vol.XX/YY, AAAA 20ZZ ISSN 2073-4859 Similarly to lattice data, geostatistical multivariate models can be fit with INLA by sharing common terms.Here, the first variable can be modeled so that the mean includes a shared spatial term assumed to be a Gaussian process with a covariance defined using a Matérn function, and all the other variables can depend on this shared spatial term plus specific spatial effects.Hence, for example, K variables of interest (Y k for k = 0, . . ., K − 1), with a general likelihood function P(•), measured at n different locations can be written as: Then, the mean of the baseline variable (µ i,0 ) will be modelled as a sum of an intercept (α 0 ) and a shared spatial effect (u i,0 ).Furthermore, the mean of observation i and variable j (µ i,j ) will be modeled through an intercept for each variable (α j ), the shared spatial effect (u i,0 ) and a specific spatial effect (u i,j ) as follows: Here, u i,0 represents the shared (between variables) spatial term, while u i,j (j ≥ 1) represents specific terms that can be used to assess departures from the shared spatial term.These random effect terms are assumed to be Gaussian processes with covariance defined using a Matérn covariance function, which for a generic spatial random effect u(s) is defined as: where s p − s l represents the Euclidean distance between points s p and s l and σ is the marginal variance of the latent Gaussian process.κ is the scaling parameter, which is related to the range, and K λ represents the modified Bessel function of the second kind and order λ, which measure the smoothness of the process.Furthermore, this latent Gaussian field u(s) can be approximated using the SPDE approach defined in Lindgren et al. (2011) that relies on the finite element method through (an appropriate choose of) deterministic basis functions defined in a triangulation of the domain where φ k (s) are the basis functions (pairwise linear functions), m is the total number of nodes (triangle vertices) and w k are zero-mean Gaussian distributed weights.For more details, the reader is referred to Lindgren et al. (2011).
In this example, it is considered that the different variables are measured at the same n spatial points, but the measurements of each variable can be located in different locations of the study region leading in a misalignment framework.
A multivariate geostatistical analysis can be performed using the R package gstat.However, the models implemented in this R package are based on a classic and frequentist statistical approach.

Point patterns
A point pattern is defined as a group of points (geographically located) which are a single realization of a stochastic process called point process.A multivariate point pattern can be defined as a group of several point patterns where each point pattern has a different origin, i.e. each point pattern is caused by different processes.These are also referred to a specific case of marked point pattern (Baddeley et al., 2015) where each point pattern is labelled with a categorical mark.
For a completely random point process, points appear independently of each other and uniformly over the study region.This is also known as a homogeneous Poisson process with (constant) intensity λ, which measures the average number of points per unit area.It is also possible to consider a spatially varying intensity, λ(x), so that the process becomes an inhomogeneous Poisson process.Other types of point processes can be more complex (sse, for example, Baddeley et al., 2015).
Several methods have been used to model the intensity function, λ(s).Here, it will be considered as a continuous process over all the study region and it will be analysed using log-Gaussian Cox processes (Møller et al., 1998;Diggle et al., 2013).Log-Gaussian Cox processes can be fit by including The R Journal Vol.XX/YY, AAAA 20ZZ ISSN 2073-4859 spatial terms using the SPDE approach implemented in INLA (Simpson et al., 2016).The analysis of the intensity as a continuous function over the study region is similar to the case of multivariate geostatistics.
Given K point patterns in a region D, an example of how to structure a multivariate point patterns model is: where λ 0 (s) is the intensity of the baseline point pattern and λ j (s) the intensity function of the jth point pattern.Moreover, α terms represent the intercepts and u(s)s are the spatial effects.Specifically, u 0 (s) is the shared (between the different point patterns) spatial term and u j (s) are the specific spatial terms which will catch the differences between each point pattern and the baseline one.
In spatial epidemiology, researchers pursue to analyse whether a distribution of cases follows the spatial distribution of a set of controls, or it depends on exposure to pollution sources or other risk factors (Palmí-Perales et al., 2021).This is an application of the model described here, the log-intensity of the controls can be modeled using a shared spatial term and the log-intensity of the cases can include this shared spatial term plus a disease-specific spatial term.Furthermore, the linear predictor can include other terms to account for risk factors and exposure to pollution sources.

Prior distribution selection
Prior choice is an essential step in a Bayesian inference process.In both, geostatistical and point patterns analysis, the SPDE approximation requires to set a prior distribution to the nominal range, r, and the nominal variance, σ.Penalized complexity prior distributions (PC-priors, Simpson et al., 2017) can be chosen for both parameters.In few words, PC priors are based on the idea of penalising the complexity from a baseline model, i.e., the prior density is related to the distance from a baseline model.A remarkable benefit of the PC-priors are their high intuitiveness in their definition.
In the case of point patters and geostatistical data, we establish a high probability of the range, r, been lower than the half of the maximum distance (d m ) of the study region, that is, P(r < d m /2) is almost 1.In the case of the variance, following Sørbye et al. (2019), an upper limit for the variation of the intensity (U α ) is considered.Therefore, the probability of the standard deviation been greater than this upper limit (P(σ > U α )) is almost 0.
In the case of multivariate lattice data, a prior distribution should be assigned to the variance/standard deviation/precision of each spatial effect.Some authors have discussed the most appropriate vague prior distributions in these cases.For instance, Gelman (2006) proposes to avoid inverse Gamma distributions on the precision and propose some alternatives.In this article, flat uniform prior distributions are assigned to the standard deviation parameters.

The structure of the multivariate spatial data in R-INLA
The spatial data sets must be formatted properly in order to be analysed with R-INLA.Their structure depends on the spatial data type.In this section, a brief explanation of how to tackle each dataset will be given.The following code will exemplify how to structure the spatial data in practice.

Areal data
Areal datasets usually are structured as a matrix where each row corresponds to a single area and each column corresponds to a single variable such as the number of cases of a specific disease.However, this is not the appropriate structure to analyse lattice data in R-INLA.
Consider a dataset with D variables measured in a study region divided in n areas.In order to analyse this dataset with R-INLA, a matrix with D columns and Dxn rows has to be built.Specifically, this matrix will store the n values of the first variable (D = 1) in the first column between the first and the nth row, then the following values of this column will be N A. In R, N A represent an empty space.
The R Journal Vol.XX/YY, AAAA 20ZZ ISSN 2073-4859 Then, the data of the second variable (D = 2) would be placed in the second column and between the n + 1 and the 2n rows.The rest of the values of the second column will be filled with N As.Therefore, the filling of the rest of the matrix will be done following this procedure.
As a toy example, let's consider n = 2 and D = 3, the original dataset is structured as a 2x3 matrix: 1 4 3 2 6 5 .
Following, the above procedure, the matrix that should be passed to R-INLA would be the following:

Geostatistics
Geostatistical datasets contain the values of the different variables and the locations where they are measured.In this case, the main objective is to estimate the value of the variables as a continuous surface over the study region.The SPDE (Lindgren et al., 2011) will be used in order to analyse the multivariate geostatistical data using a discretization (a mesh) of the surface.
Following the SPDE approach, data must be stacked using the appropriate format.The helper function inla.stack() can be used in order to built a stack object which contains the data, the effects considered in the model and the A matrix which define the location of the data on the mesh, that is, A shows which triangle of the mesh contains each data point.Specifically, a stack object will be created for each variable, then, the stack function will be used once more to combine all the stacks in a definitive stack object.
The measured values of each variable will be included in its stack using a vector which will follow the same procedure of the above subsection.As a general example, let's consider n measures and K variables.In this case, all the vectors will have nxK length and, for instance, the vector of the stack of the first variable will store the n values of the first component and the rest will be filled with N As.
As a toy example, let's consider n = 2 and K = 3, and the values of the three variables which are stored in this matrix: 1.2 4.8 3.7 2.1 6.5 5.4 The three vectors (each will be placed in the corresponding stack) are:

Point patterns
Multivariate point patterns datasets contain the locations of different point patterns.Here, the objective is to estimate the intensity surface of each point pattern.Hence, the SPDE is once more required.Therefore, the same procedure using the inla.stack()R function will be followed in this case.However, there is a special detail that has to be carefully tackled.
In the case of point patterns, the data included in the stack function have a specific structure.In this case, a list with two elements will be included for storing the data of each point pattern.
The first element of the list will be a matrix with N v + N i rows and K columns where N i are the total number of points of the ith point pattern, K is the number of different point patterns and N v is The R Journal Vol.XX/YY, AAAA 20ZZ ISSN 2073-4859 the number of points (vertices) of the mesh (triangulation) needed in the SPDE approach.Then, the matrix of the stack of the ith point pattern will be filled with N As except for the ith column.This ith column will contain firstly N v zeros corresponding with the point of the mesh.After these zeros, N i ones related with the N i points of the ith point pattern will be placed.
The second element of the list will be a vector containing an offset.Specifically, the length of this vector is also N v + N i where the first N v elements will contain the "weights" of the mesh points.The reader is referred to Simpson et al. (2016) for more details about this approximation.Then, the rest of point will be zeros.This list has to been built for each stack of the each point pattern.Then, as before, all the single stacks are combined in a single stack object.As a toy example, consider a mesh with N v = 3, two point patterns K = 2 with three and four elements, respectively (N 1 = 3 and N 2 = 4).Then the data passed to the second stack is a list with these two elements:

Multivariate lattice data
Gómez-Rubio and Palmí-Perales (2019) study the spatial risk variation of three types of cancer in peninsular Spain.In particular, they consider oral cavity, esophagus and stomach cancer at the province level.In order to assess similar spatial variation, that may lead to the identification of shared risk factors, models with shared and disease specific spatial patterns can be proposed.
Data are available in a RData file available from GitHub (see Gómez-Rubio and Palmí-Perales, 2019, for details).The following code creates the response variable (by staking the three vectors of observed cases) as well as the expected counts (other variables are needed but not shown here).Note that due to the different structure of the variables involved, data are stored together in a list object instead of a data.frame.# Load data load("dismap_sim_data.RData") # Set shorter names names(OyE.sim)[1:6]<-c("Obs.Cb", "Esp.Cb", "Obs.Eso", "Esp.Eso", "Obs.Est", "Esp.Est") As an example, we will consider a model in which the log-relative risk of oral cavity cancer (θ o,i ) is modeled using an intercept (α o ) plus a shared spatial term (u i ) using an ICAR specification following the model of the areal or lattice data analysis section.Furthermore, log-relative risks of esophagus (θ e,i ) and stomach cancer (θ s,i ) are modeled using disease-specific intercepts (α e and α s ) plus the shared spatial term and cancer-specific ICAR spatial terms (v e,i and v s,i ).Specifically, the relative risks are measured as follows: The R Journal Vol.XX/YY, AAAA 20ZZ ISSN 2073-4859 These spatial specific terms can be used to assess departures from the shared spatial term.The chosen prior distributions for the standard deviation of all the effects are flat uniform distributions.
The model formula is defined in the model below.The rf object captures the value of the intercepts (the value of the above αs ).Latent effects of type "copy" are used to define shared terms in the model.Furthermore, some latent effects have the uniform prior for the standard deviation defined in object prior.prec.
res <-inla(formula = form, data = d, family = rep("poisson", 3), E = d$EXP) Figure 1 shows the different spatial terms in the model.The shared specific term represents the spatial variation of the risk of oral cavity cancer and also serves as a baseline for the other types of cancer.The esophagus specific spatial term is quite mild, which indicates that these two types of cancer have a very similar spatial pattern.The specific spatial pattern for stomach cancer shows some provinces in the center of the country with a higher mortality than oral cavity (and esophagus).

Multivariate geostatistics
The meuse dataset in the gstat package gives the locations and measurements of topsoil heavy metals collected in a flood plain by the Meuse river, close to the village of Stein (The Netherlands).After loading the gstat package, the data can be loaded using the following code: # Load the data data(meuse) # Create the spatial object coordinates(meuse) = ~x+y proj4string(meuse) <-CRS("+init=epsg:28992") These measurements are highly correlated and we will explore in this example how to fit geostatistical models with the SPDE approach.Note that now observations do not need to be in a regular grid.
The R Journal Vol.XX/YY, AAAA 20ZZ ISSN 2073-4859 Instead of a grid, a mesh is defined to apply the SPDE approach.The boundary of the study region is stored in object meuse.bdy(see accompanying code).The definition of the mesh is done using the coordinates of the boundary of the area following the code below: # Create the mesh mesh <-inla.mesh.2d(boundary= meuse.bdy,loc = coordinates(meuse), max.edge = c(250, 500), offset = c(250, 500), n=c(32, 32)) The left plot of Figure 2 shows the mesh built with the above code for this example.In particular, the model will consider concentrations of lead and zinc.This concentrations have been measured in diferent locations which are displayed in the right part of Figure 2.However, more heavy metals can be analysed following the same structure.Values are considered in the log-scale.These log-transformed concentrations (y l and y z ) are assumed to be normally distributed.The mean of the concentration of lead is modeled using an intercept α l plus a Gaussian process with a Matérn correlation, u i,s , while the mean of the concentration of zinc is modeled using an intercept, α z (different from the previous one), the shared spatial effect, u i,s , plus another spatial Gaussian process with Matérn covariance, u i,z .This will allow us to assess differences in the spatial distribution of the concentration of both heavy metals.
Following the geostatistics subsection of the multivariate spatial models section, the model can be written as: where µ i,l and µ i,z represent the mean of the concentration of lead and zinc, respectively and which are modelled as follows: In this example, the prior choices are: P(r < 2394.16)= 0.95 P(σ > 1000) = 0.05 where a nominal range higher than the half of the maximum distance (i.e., 2394.16meters) of the domain is unlikely.Similarly for the nominal standard deviation, which is really unlikely to be higher than 1000 mg/kg of soil (ppm).These prior distributions are specified using the inla.spde2.pcmatern()function: The R Journal Vol.XX/YY, AAAA 20ZZ ISSN 2073-4859 Estimates of the posterior mean of the shared spatial effect (bottom-left) and the zinc-specific spatial effect (bottom-right).Specifically, in this example, the intensity of all of the types of forest fires is estimated by considering lightning fires as the baseline pattern.Furthermore, INLA is used to assess the similarities or differences between their spatial patterns.First of all, the mesh is built using the coordinates of the boundary bdy.SP of the dataset by following the below code.mesh <-inla.mesh.2d(boundary = list(bdy.SP,NULL), cutoff = 2, max.edge = c(20, 50), min.angle= 27, offset = c(1, 50), n=c (16,16)) Following the model structure detailed in the point patterns subsection of the multivariate models section, the log-intensity of the lightning fires λ l (s) will be modeled using an intercept, α l , and a spatial Gaussian effect with Matérn covariance, u 0 (s) as follows log λ l (s) = α + u 0 (s); s ∈ D Spatial effect u 0 (s) will also be shared in the linear predictor of the other types of fire.Similarly, the log-intensity of the accidental, intentional and the other fires (with unknown cause) will be modeled using specific intercepts, plus the shared spatial effect plus a specific spatial effect as follows: where these disease specific spatial effects can be used to assess any deviation from the shared pattern.Following the prior specification section, the chosen PC-priors in this example are: The R Journal Vol.XX/YY, AAAA 20ZZ ISSN 2073-4859 where a nominal range higher than half the maximum distance (e.g., 200 kilometres) of the domain is unlikely.Similarly for the nominal standard deviation, which is really unlikely to be higher than 100 fires per unit area in this context.These prior distributions have been specified using the function inla.spde2.pcmatern(): spde <-inla.spde2.pcmatern(mesh= mesh, prior.range= c(200, 0.95), prior.sigma= c(10, 0.05)) where the argument prior.rangesets the prior for the range and the argument prior.sigmasets the prior for the standard deviation of the spatial effect.Then, the stack objects are built.
As stated above, when analysing point patterns it is necessary to assign some weights to the observed points.This is done by creating a Voronoi tesselation using these points, so that the area of the associated polygon becomes the associated weight.The following code illustrates how to obtain the Voronoi tesselation and associated weights: require(deldir) dd <-deldir(mesh$loc[, 1],mesh$loc[, 2]) # Create a list of tiles in a tessellation mytiles <-tile.list(dd)if (!require(gpclib)) install.packages("gpclib",type = "source") pl.study <-as(bdy, "gpc.poly")# Class for polygons area.poly(pl.study)# Computing the area of the whole polygon # Compute weight as area of the polygon given as an # interaction between Voronoi tiles and domain polygon w <-unlist(lapply(mytiles, The R Journal Vol.XX/YY, AAAA 20ZZ ISSN 2073-4859 function(p) area.poly(intersect(as(cbind(p$x,p$y), "gpc.poly"),pl.study) ) ) ) These computed weights are introduced as the expected weights on the mesh points.Furthermore, the data should be structured as follows Another element of the SPDE approach is the projector matrix.When analysing multivariate point patterns this matrix has two parts: one for the mesh points (named imat) and the other for the points of the point pattern (named lmat).Then, the projector matrix is the combination of these two matrices: #imat: define imat imat <-Diagonal(nv, rep(1, nv)) #lmat: define lmat lmat.lig<-inla.spde.make.A(mesh, pts.lig) lmat.acc<-inla.spde.make.A(mesh, pts.acc) lmat.int<-inla.spde.make.A(mesh, pts.int) lmat.oth<-inla.spde.make.A(mesh, pts.oth) #Projector matrix: Put together imat and lmat A.lig <-rbind(imat, lmat.lig)A.acc <-rbind(imat, lmat.acc)A.int <-rbind(imat, lmat.int)A.oth <-rbind(imat, lmat.oth)Once all the elements of the stack object are set, the inla.stack()function is used to built the different stack objects.For instance, the stack objects of two (out of four) types of forest fires are shown  examples/tree/master/dismap_example.The idea of illustrating our objective with straightforward examples was the main cause of the choose of these models.However, more complex spatial and spatio-temporal structures can be proposed and, obviously, R-INLA is able to adjust more complex spatial and spatio-temporal models.
The main goal of this work has been to illustrate how to perform multivariate spatial Bayesian inference using R-INLA.Other alternatives based on MCMC algorithms may be highly computational demanding in a multivariate spatial context, therefore, becoming unwise in analysing larger spatial datasets.Hence, it has been shown that R-INLA is a interesting and a worthwhile toolbox to apply multivariate spatial models.Additionally, the necessary R scripts to reproduce the examples are published in https://github.com/FranciscoPalmiPerales/Mult-Sp-INLA.
In conclusion, the present work explains how to perform Bayesian multivariate spatial inference using R-INLA.We have exemplified how to perform these analysis using public datasets paying special attention to the main steps of the code.

Figure 1 :
Figure 1: Posterior mean of the shared spatial pattern (left), the esophagus-specific (middle) and stomach-specific (right) spatial patterns.

Figure 2 :
Figure 2: Mesh used in the estimation of the concentration of heavy metals around river Meuse, only with the boundary of the study region (left) and jointly with the survey locations (right).

Figure 3 :
Figure 3: Estimate of the posterior mean of the log-concentration of lead (top-left) and zinc (top-right).Estimates of the posterior mean of the shared spatial effect (bottom-left) and the zinc-specific spatial effect (bottom-right).

Figure 4 :
Figure 4: The mesh used is shown alone (top-left) and with all the fire types (bottom-left).Furthermore, the mesh is also displayed for each fire type separately: lightning fires (top-middle), accidental fires (bottom-middle), intentional (top-right) and other fires (bottom-right).

Figure 5 :
Figure 5: Posterior mean of the log-intensity of the lightning (top-left), accidental (top-right), intentional (bottom-left) and other (bottom-right) fires.

Figure 6 :
Figure 6: Posterior mean of the shared spatial effect (top-left), the accidental-specific (top-right), the intentional-specific (bottom-left) and the other-specific (bottom-right) spatial effects.
are designed to analyse lattice data.