Enhancing the SPDE modeling of spatial point processes with INLA, applied to wildfires. Choosing the best mesh for each database

Wildfires play an important role in shaping landscapes and as a source of CO 2 and particulate matter, and are a typical spatial point process studied in many papers. Modeling the spatial variability of a wildfire could be performed in different ways and an important issue is the computational facilities that the new techniques afford us. The most common approaches have been through point pattern analysis or by Markov random fields. These methods have made it possible to build risk maps, but for many forest managers it is very useful to know the size of the fire as well as its location. In this work, we use Stochastic Partial Differential Equation (SPDE) with Integrated Nested Laplace Approximation (INLA) to model the size of the forest fires observed in the Valencian Community, Spain. But the most important element in this paper is the process that needs to be carried out prior to simulating and analyzing the different point patterns, namely, the choice of the most suitable mesh for the database. We describe and take advantage of the Bayesian methodology by including INLA and SPDE in the modeling process in all the scenarios.


Introduction
Modeling the spatial 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 area burned (e.g., D ıaz-Avalos et al. 2001;D ıaz-Avalos, Juan, and Serra-Saurina 2016). Wildfires can be associated with their spatial coordinates and the corresponding covariates, and this association makes it easier to represent them as a realization of a spatial stochastic process. 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 analyze the spatial and spatio-temporal distribution of forest fires (Mandallaz and Ye 1997;Xu and Schoenberg 2011). Spatial point process theory is a rich subject that provides the mathematical fundamentals required to sustain methods that are useful to analyze phenomena with a spatial variation and make it possible to screen out the general factors that are related to such spatial variation. For a good introduction to spatial point processes, see M€ oller and Waagepetersen (2006).
Previous studies have dealt with the problem of producing risk maps or of calculating the probability of a wildfire starting at some location inside a study area D using statistical methods (e.g., Barros and Pereira 2014;. Despite their usefulness, most of the studies have not considered the burned area caused by each wildfire in the same way as in this paper. Statistical methods used to produce wildfire risk maps have included Markov Random Fields (D ıaz-Avalos et al. 2001) and spatial point processes (Møller and D ıaz-Avalos 2010;Juan, Saez, and Mateu 2012;. All these approaches have included the effect of spatial covariates. Nevertheless, the previous points in the work are the process of simulating, analyzing and modeling the three real patterns: regular, inhibition and cluster. We fitted several models, starting with models that include two simulated covariates. All the models have been compared with the use of Deviance Information Criterion (DIC) and the Watanabe-Akaike information criterion (WAIC) (Spiegelhalter et al. 2002;Watanabe 2010).
Moreover, the most important objective of this work is to show the benefits of using Stochastic Partial Differential Equation (SPDE) with Integrated Nested Laplace Approximation (INLA) for spatial data. Previous studies (Taylor and Diggle 2014) have described the difficulty of working with INLA for specific parameters, i.e., spatial parameters, but this problem could be solved with the use of SPDE in the majority of scenarios.
Previously, however, it is necessary to introduce the methodology, the simulation, and the usefulness of this methodology. This will be done step by step. Some previous papers (e.g., D ıaz-Avalos, Juan, and Serra-Saurina 2016) have identified the main problem of using this methodology as the choice of mesh, but this was only a small step. In this paper we present the procedure for the real process of analyzing the different spatial point processes more accurately, as previously defined.
Nowadays, the problem is to choose the perfect mesh formed with the SPDE of each pattern in a spatial point process, which is the main objective in this paper. Thus, it will be a valuable aid for researchers who are faced with the question: which is the correct mesh for my case?
The rest of the paper is organized as follows: Section 2, Methods, gives all the details needed to clarify the Bayesian methodology used and Spatial Point Process; Section 3, the Simulation study, is followed by the models for the simulated data in Sec. 4. Section 5 describes and discusses the results and the interpretation of applied wildfire data. We finish with Discussion and conclusions in Sec. 6. To end, the code used is presented in Appendix A and other graphical results are shown in Appendix B.

Spatial point processes
A spatial point process is a stochastic process in which we observe the locations of some events of interest within a bounded region D (Bivand, Pebesma, and Rubio 2007). Spatial point process models are useful tools to model irregularly scattered point patterns that are frequently encountered in all kinds of studies. Except in the simplest cases (Poisson processes), the points or individuals (wildfires) of a spatial pattern will be mutually stochastically dependent. The simplest of all possible models is the constant intensity Poisson process, frequently referred to as the model of 'complete spatial randomness' (CSR) (Juan, Saez, and Mateu 2012).
The next step is to decide whether the spatial point process is declared regular, cluster or inhibition, but this is not the ultimate question of this paper. This issue has been defined in many papers (e.g., Juan, Saez, and Mateu 2012) and here it is assumed that all the possibilities proposed in the spatial point process exist. Modeling and inference for spatial point processes is an issue that has been investigated broadly over the last few years. The wide range of fields of application has been the main engine driving such increased interest. Spatial Poisson processes play an important role in both statistical theory (Daley and Vere-Jones 2003) and applications (Diggle 2003). A large number of papers have discussed inference for spatial point patterns, but mostly in relation to testing the CSR hypothesis (Cressie 1993;Diggle 2003). However, in many applications, modeling issues are related to detecting patterns and interactions between points in the pattern (Baddeley and van Lieshout 1995), although this is not the most important item addressed in this paper, but instead just a previous step.

Integrated nested Laplace approximation
This work also offers the possibility of studying Spatial Point Processes in another way, by using INLA. The data can be idealized as realizations of a stochastic process indexed by: where yðÁÞ is a subset of R: The advantages of using INLA over other methods, such as basic statistical methods or more complex ones (like Markov Chain Monte Carlo (MCMC) (Tsanas and Xifara 2012), are the following: It works with reasonable computational times, thereby allowing the user to work with complex models quickly and efficiently. It allows for 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 analyzed. It does not require working with normal distributions exclusively, due to the fact that it is based on Bayesian inference.
The data can be presented by a collection of observations y ¼ fy 1 ; :::; y n g (Blangiardo et al. 2013;. A temporal correlation structure is a complicated mathematical entity and its practical estimation is very difficult if the covariates are included (Vlad, Juan, and Mateu 2015). In statistical analysis, to estimate a general model it is useful to model the mean for the unit using an additive linear predictor, defined on a suitable scale: where b 0 is a scalar, which represents the intercept, b ¼ ðb 1 ; :::; b M Þ are the coefficients of the linear effects of the covariates z ¼ ðz 1 ; :::; z M Þ on the response, and f ¼ ff 1 ð:Þ; :::f L ð:Þg is a collection of functions defined in terms of a set of covariates v ¼ vðv 1 ; :::; v L Þ: The first step in defining the structure of the data y ¼ fy 1 ; :::; y n g: A very general approach consists in specifying a distribution for y 1 characterized by a parameter u i (usually the mean Eðy i Þ) defined as a function of a structured additive predictor g i through a link function gð:Þ; such that gðu i Þ ¼ g i : The additive linear predictor g i is defined as follows (Blangiardo and Cameletti 2015): Different models can be obtained depending on the covariates, marks and pattern considered in each one. Once the battery of competing models has been obtained, the DIC and the WAIC criterium can be obtained for each of them in order to select the best one. The best models would be those with a high level of complexity and a high goodness-of-fit, that is to say, models that show the lowest WAIC and DIC should be chosen (Spiegelhalter et al. 2002;Watanabe 2010): where D h ð Þ 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 (Spiegelhalter et al. 2002). When the model is true, D h ð Þ should be approximately equal to the 'effective degrees of freedom', n À p D : For each model, a covariate X i is significant in relation to the response variable Y i if there is no change in the sign between the 0.025 quantile and the 0.975 quantile of the mean of the corresponding b 1 : In addition, the positive sign of the mean of the corresponding b 1 implies that the response variable increases when the covariate increases. The higher the mean of b 1 is, the more significant the covariate will be.

Stochastic partial differential equation
We consider that the observed data were generated by spatial point processes with covariates or marks. R software (R Core Team 2016) and the Bayesian estimation methodology, including the basic features of INLA approach and the SPDE (Lindgren, Rue, and Lindstrom 2011;Rue, Martino, and Chopin 2009), were used to fit parametric models of the observed spatial patterns of wildfires, including covariates that will play the role of risk factors showing spatial variability. This approach makes it possible to account for the effects of the risk factors on the spatial distribution of wildfire patterns (Arag o et al. 2016).
A technical difficulty to fit models using continuously indexed Gaussian random fields is the big n problem, which imposes restrictions on the size of the matrices that have to be inverted in the model fitting process. This difficulty is overcome using a discrete approximation based on SPDEs. For some Gaussian random fields with Mat ern covariance, SPDE provides an explicit link between the Gaussian random field and a discrete approximation based on a Gaussian Markov random field (GMRF) constructed on a triangulation, which can be used for computational purposes (Wist and Rue 2006;Lindgren, Rue, and Lindstrom 2011).
The SPDE approach allows a Gaussian field (GF) with the Mat ern covariance function to be presented as a discretely indexed spatial random process (GMRF), which has significant computational advantages. Gaussian Fields are defined directly by their firstand second-order moments and their implementation is highly time-consuming, which leads to the so-called "big n problem". The idea is to construct a finite representation of a Mat ern field by using a linear combination of basis functions defined in a triangulation of a given domain D. This representation leads to the SPDE approach, which is a link between the GF and the GMRF, and allows the spatio-temporal covariance function and the dense covariance matrix of a GF to be replaced with a neighborhood structure and a sparse precision matrix, respectively, both of which are typical elements that define a GMRF. This, in turn, yields substantial computational advantages (Lindgren, Rue, and Lindstrom 2011).
More particularly, the SPDE approach consists in defining the continuously indexed Mat ern GF X(s) as a discrete indexed GMRF by representing a basis function defined on a triangulation of domain D: where n is the total number of vertices in the triangulation; {u l (s)} is the set of basis functions and {x i } are the zero-mean Gaussian distributed weights. The basis functions are not random, but are instead chosen to be piecewise linear on each triangle: The key is to calculate {x i }, 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 (Simpson et al. 2016).
The problem of the choice of the best mesh for the SPDE approximation could be solved using the correct steps presented in the next sections. The steps followed in this paper are: simulation, testing the different meshes, choice of the best one for each case, modeling, and comparing the results with DIC and WAIC.

Simulation
Working with all the possibilities in spatial point processes, first of all, the data are simulated and we have Regular, Cluster and Inhibition (Figure 1), which are simulated in a regular region (type A). For each case, the numbers of points are 462, 467 and 416, respectively.
For the second type that was simulated (type B), we have an irregular region and, in this case, the numbers of points are 329, 269 and 250, respectively. The structure of the irregular region is the province of Castell on in the Valencian Community, Spain. Finally, the third row are the data simulated by rLGCP function of statstat package. In the last possibility, the data is simulated directly as Log Gaussian Cox Process as another different possibility.
The statistical description of the variables data is presented in Figure B1 in the Appendix B, with the Figures B2 and B3. These are related to the spatial distances, the most important component for the study of spatial point processes. The histogram of the distances between sites and the cumulative proportion are also shown for each case. It shows we have all types of distances in the data, which is the normal structure for spatial point processes.

SPDEs for each pattern
The next step is to choose the best mesh for the distribution of each database using SPDE. When the mesh is prepared, it is possible to use two different codes, including (or not) inla.nonconvex.hull. Indeed, the use (or not) of the boundary increases the number of cases. Next The different number of meshes is 11 possibilities for each regular and irregular region, with Regular, Cluster and Inhibition data, thus resulting in a total of 66 meshes. In the case of irregular regions, the mesh can be produced in three ways: using only data, only the boundary or both of them together. In our case, we have used the last   one, for the best model. If all the possibilities are presented, the total number of meshes is double. The numbers of vertices for each case are presented in Table 4. The dispersion in the number of vertices in each mesh is what really makes it difficult to select the best mesh.    For each case, it is necessary to decide on the best mesh. It depends on the number of points, computing time and if the data are inside the region. In the next step it will be selected with the DIC and WAIC. The next figure shows the 18 meshes selected, and the characteristics of each case can also be seen. The conditions and the differences between the best and the worst selected for each case can be observed. The differences and the importance of choosing the best mesh for modeling are quite apparent. The same is done with the data simulated with lgcp function of spatstat package. The 9 meshes selected for the regular regions are in Figure 6a-b.

Applied models
The next step is to develop the models for each case and to decide on the best one. For all cases, the models are presented with and without a spatial effect (In Table 5 the time in seconds of computation of each model is presented). For this reason, the number total of possibilities is twice that of the previous step.
For the simulation steps, two covariates are presented, both Normal distributed. Indeed, one of them is the response variable and the other, with the product of the two of them, are the two sample covariates. The different components and information for the initial case are presented in the first part. Then, all the other cases are presented together.
Case A: Regular Data with Regular region using the meshes selected in previous steps (Models M1-M3).
Case B: Regular Data with Irregular region using the meshes selected in previous steps (Models M4-M6).
The effect that variations in each simulated covariate had on the response variable was analyzed. For this purpose, INLA were obtained using the freeware statistical package R (version 3.4.3) and the R-INLA package. Models were obtained separately with or without spatial effect.
The variable response was modeled against each covariate. Models are reported in Table 6, where the structure of each model is shown in: Each cell in Table 6 reports the model for the simulated data, and for the response variables V 1 and V 1 Ã V 2 . In the same way, in this Table are presented the results for  Figure 6. (a) The 9 meshes selected for a regular region. (b) The 9 meshes selected for an irregular region. models M7: lgcp regular region and M8: lgcp irregular region data. The results are similar than the previous models (M1 to M6).
If the parameters are compared for non-spatial or spatial effect, it is clear that there are not many differences (Figures 7 and 8).
On the other hand, DIC and WAIC are presented with the same information so as to allow comparison among the models. Table 7 shows the summary results related to goodness-of-fit for the battery of models.
It could be deduced from the results obtained that in both cases including a higher number of covariates improves the model because DIC and WAIC become lower, which indicates that these models allow better prediction. This was the case, for example, when we analyzed the values of DIC and WAIC for the three models that relate V 1 with every single covariate V 1 Ã V 2 , therefore showing it is a better model. The next step is to study the spatial effect, in two formats: normal and 3 D. This is presented in Figure 9a-b.
By means of the corresponding model parameters, the spatial effect can be presented, showing the correlation with the corresponding distance ( Figure 10).
A previous step for prediction is validation. This is performed by comparing the residuals and real data, and comparing prediction and real data, as can be seen in Figure 11 below.
The residuals for lgcp data simulated are shown in next Figure 12: In the next step, the simulated data could be presented, together with the predicted model and the related standard deviation. It is shown in the next Figure 13 for a regular and irregular case for lgcp. The other cases are shown in Figure B12 in the Appendix B.
The code for all steps is presented in Appendix A.

Real data. Burned area in wildfires
A wildfire is any uncontrolled fire in combustible vegetation that occurs in the countryside or a wilderness area (Cambridge 2008). A wildfire differs from other fires by its extensive size, the speed at which it can spread out from its original source, its potential to change direction unexpectedly, and its ability to jump gaps such as roads, rivers, and fire breaks (National Interagency Fire Center 2011). Wildfires are characterized in terms of the cause of ignition, their physical properties such as speed of propagation, the combustible material present, and the effect of weather on the fire (Flannigan et al. 2006). Fire risk is very important in the Mediterranean region because it has a marked seasonality, with high temperatures and low humidity in summer. Climatic trends interact with the landscape dynamics. The process of afforestation of the different agricultural areas with the increasing abandonment of rural activities has led to a situation of extreme vulnerability to the risk of fires, especially in mountainous Mediterranean areas, where this has led to abandonment and still represents an expansion, accumulation and         greater continuity of the forest. In addition, other factors such as the development of second homes in these forest areas, the proliferation of roads and electricity networks, and an increased presence of human beings make this countryside more prone to suffer large forest fires due to the fact that they facilitate the process of starting fires (D ıaz-Delgado and Pons 2001; Moreira, Rego, and Ferrera 2001). Given that wildfires are a natural element of the Mediterranean ecosystem, the prevention and suppression of forest fires should be addressed in order to lower the risk and vulnerability levels to values that are tolerable for society.
If every 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, and controlling for the moment in time in which it was produced, one can identify all wildfires with a spatial-temporal stochastic process. Such processes, called spatial point processes, often display dependence between the spatial positions and moments in time, and interdependence between them. What is usually of more interest is the detection of trends in the intensity of fire locations, and determination of how (or whether) such trends are influenced by covariates, observable at each location of the spatial window. These covariates might include vegetation or land use, other descriptors of terrain (such as elevation, slope and orientation), and still others such as proximity to concentrations of human population or to concomitants of human activity (roads and railroads). Interaction between points may in general be of some interest in its own right. More important is the impact of the presence of interaction on statistical inference concerning trends and their dependence upon covariates. Temporal clustering of wildfires, whether deriving from multiple ignition lightning events, arson (Butry and Prestemon 2005), or other sources, combined with favorable fuel and weather conditions, can force suppression resource rationing across space. Spatial clustering can also indicate the presence of risk factors.

Data setting
We analyze the patterns produced by wildfire incidents in the Valencian Community, located on the north-east coast of the Iberian Peninsula. The region is bordered by Catalonia 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 kilometers, representing 4.6% of the Spanish national territory.
A total of 315 fires were recorded in the study area during 2015 ( Figure 14). 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.
The normal steps used for exploration of data in such cases are: 1. Data exploration 2. Outliers 3. Collinearity 4. Relationships Y vs. X Figure 15. Boxplot covariates: cause, days last rain, temp_max, Relative_H, wind speed, wind direction, combined model, danger degree, fire type.

Variance Inflation factors
6. Interactions (is the quality of the data good enough to include them?) 7. Zero inflation justification.  Figure 16. Response vs. covariates.
The covariates, related to wildfires, used for the application are: typology, cause, causative, days last rain, temp_max (maxim temperature), Relative_H (Relative Humidity), wind speed, wind direction, combined model, danger degree, fire type, total tree, total not tree, Total (tree and not tree). The relationships between the covariates are shown in Figure 15 (boxplots). The next step is the Variance Inflation factors (GVIF), seeking the optimal variables and using only values below 2 because these are informative (Table 8).
On representing the covariates in order to study the collinearity, no problematic structures or outliers are found. The next step, which is necessary for the study of real data, is the relationships. The most likely option to obtain patterns, however, is the relationship between each covariate and the response variable. In this case, the response variable is TOTAL (tree and not tree together). This is represented in Figure 16 and there is no pattern present.
In the same way, the distances are studied, and different distances are found, which gives us the chance to continue the study including all of the referenced points ( Figure 17).  Finally, before selecting the model, the decision regarding the model selected is important if the number of 0 in the data is high, since in this case the model will be Zero Inflation. There is 42.22% of zero data, and we will use a Zero-inflated model ( Figure 18). Figure 19. The 13 meshes, with respective numbers of vertices: 1612, 1430, 1612, 2021, 1732, 1356, 1592, 2071, 794, 1778, 600, 1836, 2801. Following the steps presented above, we use with simulations data is the SPDE, the same number of meshes, 9 þ 4 ( Figure 19) and Figure 20 below shows the two selected 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 this affects the computing time (Table 9). Table 10 shows the value of the parameters for each model for real data. In Figure 21, the difference between the values of the parameters is apparent, and it is important for the final results. Table 11 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 the model because DIC and WAIC become lower, which indicates that these models enable better prediction. Figure 22 shows the spatial effect of the Valencian Community in two formats.  In this case, validation is performed by comparing residuals and with the correlation between the real data and the model data (Figure 23 up). The relationship with the distances can also be seen in Figure 23. The correlation in this case, q ¼ 0.8309389, has a high value, which suggests a good result.

Conclusions and discussion
The problem of Bayesian Inference using SPDE for spatial point processes is considered in this work. The steps proposed for considering the best solution for each case are presented and modeled using simulated process models, where the latent Gaussian field is extended to a GRMF. A computationally efficient method for Bayesian inference, based on INLA and SPDE, was presented. We looked at an application for modeling a wildfire dataset in which the process is faster and clearer. The example was of interest since the point pattern shows clear signs of being affected by some covariates and spatial elements.
In a more general context of spatial point process modeling, the data structure we considered here differs from the point patterns typically analyzed in the point process literature, but it is important to show that this SPDE methodology affords researchers easier and computationally efficient results. It is very important to note the steps used in the process and the reduction in computational time.
The advantage of INLA-SPDE is that it is able to 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. Therefore, the use of triangulated meshes in INLA-SPDE may lead to the possibilities of working with simple or complex databases.
The purpose of this study was to demonstrate that the use of the INLA-SPDE package (Rue, Martino, andChopin 2009, andLindstrom 2011) to analyze all kinds of point patterns is faster, clearer, and computationally effective. The results showed that INLA-SPDE can 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 in a straightforward way with relatively few lines of code.
As a conclusion, the SPDE methodology is a relevant element for the study of the spatial data, but it could be improved with the inclusion of the temporal effect as an element not separately with the space. In the same way, the spatial effect included in the model by SPDE is difficult because there are many elements included in this methodology and it is not well known by the researchers. It is shown in the results with wildfires because the small change in a parameter could modified the results. Nevertheless, the use of INLA-SPDE may lead to two possibilities, working with simple or complex databases, and this process can be potentially applied to mapping the spatial distribution of spatial data in all kinds of point patterns, in geostatistical data and including covariates. type ¼ "l",xlab ¼ "Distance between sites",ylab ¼ "Cumulative proportion") text(0.2, 0.9, "B", cex ¼ 1.5)