Spatio-temporal hierarchical Bayesian analysis of wildfires with Stochastic Partial Differential Equations. A case study from Valencian Community (Spain)

ABSTRACT The spatio-temporal study of wildfires has two complex elements that are the computational efficiency and longtime processing. Modelling the spatial variability of a wildfire could be performed in different ways, and an important issue is the computational facilities that the new methodological techniques afford us. The Markov random fields methods have made possible to build risk maps, but for many forest managers, it is more advantageous to know the size of the fire and its location. In the first part of this work, Stochastic Partial Differential Equation with Integrated Nested Laplace Approximation is utilised to model the size of the forest fires observed in the Valencian Community (Spain) and so it does the inclusion of the time effect, and the study of the emergency calls. The most crucial element in this paper is the inclusion of the improved meshes for the spatial effect and the time, these are, 2d (locations) and 1d (time) respectively. The advantage of the use of spatio-temporal meshes is described with the inclusion of Bayesian methodology in all the scenarios.


Introduction
Modelling the incidence of wildfires and the fire size is necessary to understand how global warming and climate change may affect the landscape in the coming years, and to determine what factors are related to spatial incidence and the size of the burned area [10,11]. Wildfires are associated with their spatial coordinates, the time effect and the corresponding covariates. Thus, even though methods such as Markov Random Fields may also be useful to respond to some scientific questions of interest, spatial point processes are the most appealing analytical tool to investigate the spatial and spatio-temporal distribution of forest fires [16,30].
Previous studies have solved the wildfires problem by producing risk maps or by calculating the probability of a starting wildfire at some location inside a study area D using statistical methods [2,23,24]. These studies used statistical methods to produce wildfire risk maps included Markov Random Fields [11] and spatial point processes [13,18]. Despite their usefulness, most of the studies have not considered the burned area caused by each wildfire as in this paper given that it is used INLA-SPDE.
To begin with, the first part of the work is devoted to the study of the best spatial mesh for the data. The models that are introduced in this work have been compared with the use of the Deviance Information Criterion (DIC) and the Watanabe-Akaike information criterion (WAIC) [26,29].
In addition, the most essential point of this work is to show the benefits of using Stochastic Partial Differential Equation (SPDE) with Integrated Nested Laplace Approximation (INLA) for spatio-temporal wildfire data [14,15,17,21,22], including the mesh for the temporal effect, one dimension (1d -time) and for the spatial effect (2d -locations). Nowadays, the problem is to choose the perfect mesh formed with the SPDE of each pattern in spatiotemporal processes. A prerequisite for creating mesh processes is the exploratory data analysis which is fundamental in this research.
Moreover, the noteworthy elements in the paper are the different data used to this methodology (INLA -SPDE). The outline applied to wildfires in Valencian Community includes two different data set. Firstly, are the real locations and temporal wildfires, and secondly the emergency calls about the corresponding wildfires.

Data set
The patterns produced by wildfire in the Valencian Community are analysed and its location is in the north-east coast of the Iberian Peninsula. The region is bordered by Catalonia Figure 2. Boxplot covariates, in the first line are, cause, days last rain, temp_max, in the second line are Relative_H, wind speed, wind direction and in the third line are combined model, danger degree, a fire type.  1.199396 1.194232 to the north and the Iberian System range of mountains to the west. Furthermore, the region is delimited to the east by the Mediterranean Sea. It is a region with a surface area of 23,245 square kilometres, representing 4.6% of the Spanish national territory. A total of 315 fires were recorded in the studied area in 2015 ( Figure 1). In addition to the locations of the fire centroids in Cartesian coordinates (Mercator transversal projections, UTM, Datum ETRS89, zone 31-N), several covariates were also considered.
Data exploration steps used in such cases are: (1) Exploratory Data Analysis (2) Outliers (3) Collinearity (4) Relationships between the response variable (Y) and the covariates (X's) (5) Variance Inflation factors (6) Interactions (7) Zero inflation justification The covariates that directly affect wildfires are typology, cause, causative, days last rain, maximum temperature, relative humidity, wind speed, wind direction, combined model (relation between elements in a wildfire), danger degree and a fire type. The relationships between the covariates are shown in Figure 2 (boxplots). The next step is the Variance Inflation factors (GVIF) and seeks the optimal variables and using only values below 2 because these are informative (Table 1). The collinearity study of the covariates, the possible structure of outliers' problems and the possible pattern appear in Figure 3. The most likely option to obtain patterns is the relationship between each covariate and the response variable TOTAL (burned tree and not tree together), and as a result, in this illustration there is no pattern present.
In Figure 4 all the possible distances between points, extensive or short, are studied and this enables us to continue with the study.
Finally, before selecting the model, the number of zeroes determines if the variable data is high. The decision regarding the selected model depends on the number of zeroes found, and if there are many of them proximate zero or zero, the best likelihood family is Zero-inflated. In our case, there is 42.22% of zero data, and it is unavoidable the use of Zero-inflated model ( Figure 5).
On the other hand, the second data used which reinforces the analysis of the wildfires in the Valencian Community, is the use of the emergency calls in the same period for the   incidences. In Figure 6 appear the locations of the emergencies calls. The main characteristic of this data set is the position of a real wildfire which is proximate to any population, and in consequence it is necessary an early intervention and resolution of the problem.
The rest of the paper is organised as follows. Section 2, Methodology, gives all the details needed to clarify the Bayesian methodology used and Spatial Point Process. Section 3 is devoted to the Data set, and Section 4 includes the models for burned area among different wildfires scenarios. Finally, Discussions and Conclusions are in Section 5.

Integrated nested Laplace approximation (INLA)
This work offers the possibility of studying Spatial Point Processes by using integrated nested Laplace approximation (INLA) [5,6]. Rue et al. [22] develops the INLA methodology for approximate Bayesian inference as an alternative to traditional Markov chain Monte Carlo methods. INLA focuses on models that can be expressed as latent Gaussian Markov random fields (GMRF) for their computational properties, and the data applied in this case possesses such characteristics.
The data can be idealised as realisations of a stochastic process indexed by: where s i represents spatial and t i is for temporal, with both of them Y(·) is a spatio-temporal subset of R 2 × R.
The advantages of using INLA over other methods, such as basic statistical methods or more complex ones (like Markov Chain Monte Carlo (MCMC) [27]), are the following: • It works with reasonable computational times, thereby allowing the user to work with complex models quickly and efficiently. • It allows the integration of as many covariates as desired, and also the incorporation of new covariates in the model in later steps. • It allows the level of significance of covariates to be analysed.
• It does not require working with normal distributions exclusively, since its base is on Bayesian inference.
The data can be presented by a collection of observations y = {y 1 , . . . y n } [4][5][6]9,28]. In statistical analysis, to estimate a general model it is useful to shape the mean for the additive linear predictor, defined on a suitable scale: , different from the previous covariates. The first step in defining the structure of the data y = {y 1 , . . . y n }. A very general approach consists in specifying a distribution for y = {y 1 , . . . y n } characterised by a parameter φ i (usually the mean E(y i )) defined as a function of a structured additive predictor η i through a link function g(.), such as g(φ i ) = η i . The additive linear predictor η i is defined as follows [3]: where β i represents the coefficient that quantifies the effect of the covariates in the response z i . This statistical analysis can be carried out with the freeware statistical package R, version 3.4.3 [20] and the R-INLA package 2017 [21]. The priors for the formulas (1) to (3), are an important element. The fixed effects and are typically normally distributed, centred on 0 and with a large variance, while ν ij are called random effects(hyperparameters) and are typically normally distributed with an exchangeable structure, i.e. ν ij ∼ Normal(0, σ ν 2 ). With this, a prior distribution needs to be specified on the regression parameters β = {β 0 , . . . , M } including the intercept, and on the variance σ 2 of the outcome. The choice of prior is:    The aim is to perform the inferential process and to obtain the posterior distribution for β and σ 2 .
If we are interested in changing the prior for the regression parameters, for instance, reducing the variability on the prior for β 0 and β 1 , specifying β 0 ∼ Normal(0, 10,000) and β 1 ∼ Normal(0, 1), we can achieve it in R -I NLA using the option control.fixed. It is also possible to modify the specification of the prior on the outcome precision (rememberτ = 1\σ 2 ) using the option control.family of the inla command. By default, a noninformative logGamma prior is assumed on the logarithm of the precision, which is equivalent to assume a Gamma prior on the precision τ ∼ Gamma(1, 10 −5 ).
In the case of the temporal correlation is considered here the commonly used random walk (RW). This random walk structure is characterised by a variance parameter, on which we need to specify a prior distribution, in our case logGamma(1, 10 −5 ). In R-INLA the default internal representation for the SPDE parameters is log(τ ) = θ 1 and log(κ) = θ 2 , with θ 1 and θ 2 being given a joint Normal prior distribution (by default independent Normal(0, 1) priors are used).
When the battery of competing models has been obtained, the DIC and the WAIC criterium can be obtained for each one of the models to select the most suitable one, those  that occur to have a higher level of complexity and a greater goodness-of-fit. That is to say, models that show the lowest WAIC and DIC should be chosen [26,29]: where D(θ ) is the deviance evaluated at the posterior mean of the parameters and p D denotes the 'effective number of parameters', which measures the complexity of the model [26]. When the model is true, D(θ ) should be approximately equal to the 'effective degrees of freedom', n − p D .'

Stochastic Partial Differential Equation (SPDE)
The Stochastic Partial Differential Equation (SPDE) approach is used for the study of spatial effect with the Matérn covariance function. The triangulation presented allows the spatio-temporal covariance function and the dense covariance matrix of a Gaussian Field (GF) to be replaced with a neighbourhood structure and a sparse precision matrix. This yields substantial computational advantages [15,24], and this approach makes possible to detect the risk factors effects in the spatial distribution of wildfire patterns [1]. In this case, the covariance structure of the Matérn type for the dispersion matrix , that is, if h ij = ||x i -x j || denotes the distance between two arbitrary wildfires within W the covariance of their fire sizes is given by where K τ denotes a Bessel function of second kind and order τ , which controls the smoothness of the process. The parameters τ and κ relate empirically to the nominal range of the Figure 11. Spatial effect.  spatial covariance (r = √ 8τ /κ). The Matérn covariance is a general model that encompasses covariance models such as the Exponential and the Gaussian, commonly used in geostatistical analyses [14].
The structure and the basis functions used are defined on a triangulation of domain D: where n is the total number of vertices in the triangulation; {ϕ l (s)} is the set of base functions and {ω l } are the zero-mean Gaussian distributed weights. These basis are presented as:  The key is to calculate {ω l }, which reports the value of the spatial field at each vertex of the triangle. The values inside the triangle will be determined by linear interpolation [24,25]. The choice problem of the best mesh for the SPDE approximation could be solved using the correct elements of the mesh shown in Figure 7, where it is different for 1d mesh used for the temporal effect because it is easier for a one-dimension mesh.

Spatial mesh with SPDE
A wildfire is any uncontrolled fire in combustible vegetation that occurs in the countryside or a wilderness area [8,13]. A wildfire differs from other fires by its extensive size, and they are characterised in terms of the cause of ignition, their physical properties or the weather effect on the fire [12,19].
Wildfires are a natural element of the Mediterranean ecosystem, and their prevention and suppression is in the benefit of lowering the levels of risk and its vulnerability to values that are tolerable for society [13].
A wildfire is associated with its spatial coordinates, longitude and latitude of the centroid of the burned area or the place where it was detected, along with other variables such as size or cause of the forest fire. The spatial-temporal stochastic is the process by which controlling the moment in time when it was produced, all wildfires can be identified. Temporal clustering of wildfires, whether deriving from multiple ignition lightning events, arson [7], or other sources, combined with favourable fuel and weather conditions, can force suppression resource rationing across space. Spatial clustering can also indicate the presence of risk factors. The temporal effect is included as a covariate mesh (1d) in the model.
The steps for modelling the application that includes the spatial effect created with the mesh, are firstly creating the spatial locations, by Matérn covariance, and then, is created the spatial mesh structure [15]. In the next steps, the covariates are included, and the SPDE spatial model is done and introduced as a function in the final model [22].
Following the possibilities studied in previous works, the meshes are shown in Figure 8, and in Figure 9, following the natural instructions [14], are selected the best meshes. The models applied for the real data, that is, wildfires in the Valencian Community in 2015, have both non-spatial and spatial effects and these affects the computing time (Table 2). Table 3 shows the value of real data parameters for each model. In Figure 10, the difference between the values of the parameters appears, and it is crucial for the final results. Table 4 below shows the summary results related to goodness-of-fit for the battery of models. DIC and WAIC.
With these results, it can be deduced that including a higher number of covariates improves any statistical model because DIC and WAIC become lower and provides a better prediction.
In this case, validation is performed by comparing residuals and the correlation between the real data and the model data. The relationship with distances can also be seen in Figure 11, above. Figure 12 bellow shows the spatial effect of the Valencian Community in two formats. The correlation in this case, ρ = 0.8309389, has a high value which suggests a good result.

Spatio-temporal meshes
The second part for the applied data is the analysis of emergency calls about the wildfires in 2015 in the Valencian Community. In this case, for the inclusion of the time effect, the temporal 1d mesh is developed in 4 parts (blue lines in Figure 13). The black lines are the mouthparts of the data. The next step is developing the spatial 2d mesh (left) and regionalised mesh (right) when the contour of the region is not used (Figure 14). The leftward mesh is used for the inclusion of the spatial effect in the models and this reduces the computing time.
With these data and with k = 4 (four temporal parts), the models tested are shown in Table 5 and the structure of their parameters in Figure 15. The inference to the parameters is another advantage of the use of Bayesian methodology. In Figure 16 the region models are presented without the use of region contour.  Then, the introduction of the contour in the model defines the studied region. The 2d mesh is different, and the outside part of the region is not necessary. In Figure 17 the new mesh and regionalised zones are shown, and the red dots are the regions without points.
With these data and with k = 4 (four temporal parts), the models tested are (Table 6 and Figure 18): Finally, in Figure 19 the models with the contour included is presented.

Discussions and conclusions
In this work, the study of wildfires with Bayesian methodology, including SPDE for spatial and temporal effect, is done. The phases proposed for considering the best solution for each case are presented and modelled using latent Gaussian field which is extended to Gaussian Random Markov Fields. A computationally efficient method for Bayesian inference, based on INLA and SPDE, was presented and the newest element is the inclusion of meshes in time (one dimension -1d) and space (two dimensions -2d). We looked at an application for modelling a wildfire data set in which the process is faster and more precise. This method is faster than the other ones proposed in previous works, and it a basic advantage for new research.
The advantage of INLA-SPDE is that it can predict the subsequent marginal distributions of the model parameters as well as the model responses without carrying out extensive simulations. This methodology can be potentially applied to mapping the spatial distribution of environmental variables or in all kinds of spatial point patterns in geostatistical issues, including covariates. The use of triangulated meshes in INLA-SPDE may lead to two possibilities, working with simple or complex databases.
As a conclusion, the wildfire data and emergency calls have a similar behaviour in space and time as it is shown in this study. Hence, it is relevant for following research since data set about real locations or emergency calls of wildfires depend on the accessibility of getting both. The obtained models in this investigation show that the time covariate is an essential element for the behaviour of the wildfires, and not only the other covariates as the elevation, human effects or climatological effect.
The results show that INLA-SPDE could also be a complementary tool in the wildlife biologist's analytical toolkit, where models are specified using a syntax that should be familiar to users of R, and where data are formatted straightforwardly with relatively few lines of codes, and that implies a leading advance in spatial statistics.