Goodness-of-fit test for point processes first-order intensity

Computational Statistics and Data Analysis 194 (2024) 107929 Available online 1 February 2024 0167-9473/© 2024 The Author(s). Published by Elsevier B.V. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/). Contents lists available at ScienceDirect Computational Statistics and Data Analysis journal homepage: www.elsevier.com/locate/csda Goodness-of-fit test for point processes first-order intensity M.I. Borrajoa,∗, W. González-Manteiga a, M.D. Martínez-Miranda b aCITMAga, Department of Statistics, Mathematical Analysis and Optimization, Universidade de Santiago de Compostela, Spain bDepartment of Statistics and Operations Research, Universidad de Granada, Spain A R T I C L E I N F O A B S T R A C T Keywords: Point processes First-order intensity Goodness-of-fit Covariates Modelling the first-order intensity function is one of the main aims in point process theory. An appropriate model describes the first-order intensity as a nonparametric function of spatial covariates. Aformal testing procedure is presented to assess the goodness-of-fit of this model, assuming an inhomogeneous Poisson point process. The test is based on a quadratic distance between two kernel intensity estimators. The asymptotic normality of the test statistic is proved and a bootstrap procedure to approximate its distribution is suggested. The proposal is illustrated with two applications to real data sets, and an extensive simulation study to evaluate its finitesample performance. 1. Introduction The understanding of spatial point processes is crucial in many different fields: ecology (Illian et al., 2009 and Law et al., 2009); epidemiology (Lawson, 2013); seismology (Ogata and Zhuang, 2006 and Schoenberg, 2011); forestry (Stoyan and Penttinen, 2000) and geology (Foxall and Baddeley, 2002). Methodological statistics have also approached this field from different perspectives, see for example Daley and Vere-Jones (1988)focused on a measure theory approach, and, Diggle (2013)and Cressie (2015) advocating for a more exploratory perspective. A spatial point process is a stochastic process governing the location of a random number of events, 𝐗 ={𝐱1, … , 𝐱𝑁}, irregularly placed in a planar region 𝑊⊂ℝ2. Throughout this paper, point processes and patterns are denoted in bold capitals, and events are denoted in bold. The spatial distribution of events in Poisson point processes, i.e., those with independent events, is determined by the first-order intensity function, see Illian et al. (2009)and Diggle (2013). This function measures the expected number of events per unit area, and is defined as follows: 𝜆(𝑥) = lim |𝑑𝑥|→0 𝔼[𝑁(𝑑𝑥)] |𝑑𝑥|, where 𝑁is the counting measure, i.e., 𝑁is the random variable that counts the number of points lying on a given set, 𝔼denotes the expectation and | ⋅|is the Lebesgue measure of the corresponding space, area in the planar case. Poisson point processes have been widely studied. A Poisson point process verifies, by definition, that the number of events in any bounded region of the process’s underlying space follow a Poisson distribution. Moreover, considering a collection of disjoint and bounded subregions of the underlying space, the number of events of a Poisson point process in each bounded subregion is independent of all the others. This property is known as complete randomness, complete independence or independent scattering. In other words, there is a lack * Corresponding author. E-mail address: [email protected] (M.I. Borrajo). https://doi.org/10.1016/j.csda.2024.107929 Received 2 May 2023; Received in revised form 24 January 2024; Accepted 25 January 2024 Computational Statistics and Data Analysis 194 (2024) 107929 2 M.I. Borrajo, W. González-Manteiga and M.D. Martínez-Miranda of interaction among the events, which motivates the Poisson process being sometimes called a completely random process and fully characterised by its first-order intensity function. However, events in a generic spatial point process may not occur independently, and interactions between them are characterized through second-order properties such as the K-function, see Ripley (1977), and the second-order intensity function, 𝜆2(𝑥, 𝑦) = lim|𝑑𝑥|→0,|𝑑𝑦|→0 𝔼[𝑁(𝑑𝑥)𝑁(𝑑𝑦)] |𝑑𝑥||𝑑𝑦|. The conditional intensity function 𝜆𝑥(𝑥)=𝜆𝑐(𝑥|𝐲)= 𝜆2(𝑥, 𝐲) 𝜆(𝐲) determines the intensity at a point 𝑥conditional on having an event at 𝐲, see Diggle (2013), and characterizes uniquely the distribution in any spatial point process. For spatial Poisson point processes, 𝜆𝑐(𝑥, 𝐲) =𝜆(𝑥). To estimate the first-order intensity function the classical approach consists of fitting a parametric function by maximum likelihood, see Waagepetersen (2007), Moller and Waagepetersen (2003)and Diggle (2013), or pseudo-likelihood procedures, see Van Lieshout (2000). It is well-known that parametric solutions may generate unreliable estimates which do not represent the underlying true model. Hence nonparametric approaches are valuable and powerful alternatives. Diggle (1985) proposed the first kernel intensity estimator which has been widely used in exploratory analysis. The main disadvantage of this estimator is its lack of consistency. To overcome this issue (Cucala, 2006, Sect. 5.3) proposed an improvement based on the relationship between intensity and density functions: both are non-negative functions and the difference relies on the intensity function not necessarily integrating the unit. Recently, real application requirements have motivated a new scenario based on the inclusion of covariates. Considering spatially varying covariates has been a big step forward on point process theory, and it has been mostly addressed so far from a parametric perspective. See for example Waagepetersen (2007)performing inference on Neymann-Scott processes depending on a linear combination of spatial covariates; Guan and Loh (2007) estimating coefficients of covariates using the Poisson likelihood estimator, previously proposed in Schoenberg (2005); and Guan and Shen (2010) proposing a method to estimate the coefficients of the covariates improving the Poisson likelihood estimator. In the nonparametric context the inclusion of covariates has been less prolific. Guan (2008) proposes a kernel intensity estimator, assuming that the intensity function depends on some observed spatially varying covariates through an unknown continuous function. Later Baddeley et al. (2012) postulates the intensity model 𝜆(𝑥)=𝜌(𝑍(𝑥)),𝑥∈𝑊⊂ℝ2,(1) where 𝑍∶𝑊⊂ℝ2→ℝis a spatial continuous covariate that is exactly known in every point of the region of interest 𝑊, and 𝜌is an unknown real function. This proposal provides intensity estimators based on local likelihood and kernels. Borrajo et al. (2020)have detailed the theoretical properties of a related kernel intensity estimator under this intensity model, including extensions to more than one spatial covariates and appropriate data-driven bandwidth selectors, which was absolutely missing so far. Little attention has been paid to test the goodness-of-fit of intensity models with covariates. Some proposals on model selection and testing covariate significance have been made. Under a parametric approach, Yue and Loh (2015)and Thurman et al. (2015) focused on Neymann-Scott and cluster processes, respectively. A general formulation proposed by Díaz-Avalos et al. (2014) assumed that the intensity depends on a linear combination of several covariates, and they proposed a test to determine whether any of the coefficients associated to the covariates may be null. While there might exist exploratory techniques to check the adequacy of the covariate representation of the intensity (see for example Baddeley et al., 2015), to the extend of our knowledge, no formal goodness-of-fit test has been proposed in the literature. In this work we define such a test based on the comparison through a quadratic distance of the intensity model including the covariates, and the “classical” nonparametric model that only takes into account the spatial coordinates. Under the Poisson assumption, our procedure is well defined, and we are able to derive its asymptotic properties, as well as a convenient bootstrap calibration. The extension of the test to non-Poisson scenarios is also possible by specifying some more characteristics for the underlying process. This paper is organised as follows. In Section 2we introduce additional notation and briefly describe the nonparametric intensity estimators that we use to define the test statistic. Section 3presents our goodness-of-fit test, providing the asymptotic distribution of the test statistic, and a bootstrap procedure to approximate its critical values in practice. The proposed methodology is applied to real data sets in Section 4. The finite-sample performance of the test is analysed in Section 5through an extensive simulation study, where the simulated models have been defined based on the real data sets previously presented. Extensions to non-Poisson processes are discussed in Section 6, and some conclusions are drawn in Section 7. 2. First-order intensity modelling and estimation Let 𝐗be a Poisson point process defined in a region 𝑊⊂ℝ2, where 𝑊is assumed to have finite positive area. Let 𝐗1, … , 𝐗𝑁be a realisation of the process, where 𝑁is the random variable counting the number of events. Diggle (1985) proposed the first kernel intensity estimator for one-dimensional point processes, which has been easily extended to the plane as:  𝜆𝐷 𝐻(𝑥)=∑𝑁 𝑖=1 𝐊𝐻(𝑥−𝐗𝑖) 𝑝𝐻(𝑥),𝑥∈𝑊⊂ℝ2, where 𝐻is a two-dimensional bandwidth matrix, 𝐊𝐻(𝑥) =|𝐻|−1∕2𝐊 (𝐻−1∕2𝑥), with 𝐊being a bivariate kernel function, and 𝑝𝐻=´𝑊|𝐻|−1∕2𝐊(𝐻−1∕2(𝑥 −𝑦))𝑑𝑦 is the edge correction term. Here |𝐻|denotes the determinant, and 𝐻−1∕2 a power of matrix 𝐻. Computational Statistics and Data Analysis 194 (2024) 107929 3 M.I. Borrajo, W. González-Manteiga and M.D. Martínez-Miranda This intensity estimator has been widely used during decades for exploratory analysis, but its lack of consistency has limited the inference performed with it. One approach to overcome this problem involves estimating the “density of events locations” or relative density. In this scenario, the Poisson assumption is needed to develop the asymptotic theory, and guarantee the consistency of the kernel estimator of the density of event locations. Notice that we cannot distinguish between hetereogeneity and interaction in an observed point pattern unless we have some additional information, such as a parametric model, Diggle (2013). As it has been pointed out by Fuentes-Santos et al. (2015), the common practice in the analysis of spatial point patterns is assuming that the point process is Poisson, estimating the first-order intensity, and then estimating the second-order properties to test the Poisson assumption. Therefore, the Poisson assumption is not very restrictive in practice. Given that the number of events of an inhomogeneous Poisson point process, 𝑁, has distribution 𝑃 𝑜𝑖𝑠 (´𝑊𝜆(𝑥)𝑑𝑥)≡𝑃 𝑜𝑖𝑠(𝑚), Cucala (2006) defined the density of event locations, 𝜆0(𝑥) =𝜆(𝑥)∕𝑚. Using kernel estimators, this idea has been followed by Fuentes- Santos et al. (2015), and the kernel estimator of the relative density is defined:  𝜆0,𝐻 (𝑥)= 1 𝑝𝐻(𝑥)𝑁|𝐻|−1∕2 𝑁 ∑ 𝑖=1 𝐊(𝐻−1∕2(𝑥−𝑋𝑖))1{𝑁≠0}.(2) Closely related to the above estimator, Borrajo et al. (2020)have proposed a consistent kernel intensity estimator using covariates for the multivariate model 𝜆(𝑥)=𝝆(𝐙(𝑥)),𝑥∈𝑊⊂ℝ2,(3) where 𝝆is an unknown 𝑝-dimensional function and 𝐙 =(𝑍1, … , 𝑍𝑝)is a 𝑝-dimensional covariate. Each component, 𝑍𝑗∶𝑊→ℝ, is a spatial continuous covariate that is exactly known in every point of the region of interest 𝑊. Under model (3), the grounds of Borrajo et al. (2020)’s proposal to estimate the intensity function, 𝜆, rely on a result stating that, if 𝑋is a point process in 𝑊⊂ℝ2fulfilling equation (3)and 𝐙is defined as above, then the transformed point process 𝐙(𝑋) is a point process in ℝ𝑝with intensity 𝝆𝐠⋆. In this context, 𝐠⋆is defined as 𝐠⋆(⋅) =|𝑊|𝐠(⋅), where 𝐠is the first derivative of the spatial cumulative distribution function of 𝐙and 𝐆(𝐳) =1 |𝑊|´𝑊1{𝐙(𝑢)≤𝐳}𝑑𝑢. Hence, the density of events locations associated to the transformed point process is 𝐟(⋅) =𝝆(⋅)𝐠⋆(⋅) 𝑚with 𝑚 =´ℝ𝑝𝝆(𝐳)𝐠⋆(𝐳)𝑑𝐳=´𝑊𝜆(𝑥)𝑑𝑥. With these definitions a kernel estimator for the density of events locations 𝐟is defined as follows:  𝐟𝐵(𝐳)=𝐠⋆(𝐳)1 𝑁 𝑁 ∑ 𝑖=1 1 𝐠⋆(𝐙𝑖)𝐋𝐵(𝐳−𝐙𝑖)1{𝑁≠0},(4) where 𝐋is a multivariate radially symmetric kernel function, and 𝐵is a 𝑝-dimensional bandwidth matrix. From (4)the kernel estimator for the intensity function 𝜆in model (3)is given by:  𝜆𝐵(𝑥)= 𝑁 ∑ 𝑖=1 1 𝐠⋆(𝐙𝑖)𝐋𝐵(𝐙(𝑥)−𝐙𝑖).(5) 3. The proposed test The first-order intensity assumption detailed in (1)has been on the grounds of several proposals, see for example Baddeley et al. (2012), Baddeley et al. (2013)or Borrajo et al. (2020). No formal testing procedure has been proposed, and only graphical diagnostics for validating the covariate effect term in the model are available so far. In this scenario, we want to test the null hypothesis 𝐻0∶𝜆(𝑥)=𝝆(𝐙(𝑥)),𝑥∈𝑊, versus a general alternative in which the intensity function is not explained completely by the covariates, i.e., 𝐻1∶𝜆(𝑥)≠𝝆(𝐙(𝑥)),for some 𝑥∈𝑊. For simplicity in the process’ notation, and to be able to apply consistency results in Fuentes-Santos et al. (2015), our proposal is developed in ℝ2. However, the ideas and formulas can be directly extended to a generic Euclidean space, ℝ𝑑, taking into account that an extra effort on proving the consistency of kernel methods in ℝ𝑑needs to be done. As it has been introduced in Section 1, with Poisson assumption, the process under the null hypothesis is well-defined. However, in Section 6we detail some ideas on a more general version of our proposal in terms of the nature of the process. 3.1. Test statistic For convenience we express the problem in terms of the relative density of event location instead of the corresponding intensity function. Notice that the null hypothesis can be equivalently rewritten as 𝐻0∶𝜆0(𝑥)=𝝆(𝐙(𝑥))∕𝑚, Computational Statistics and Data Analysis 194 (2024) 107929 4 M.I. Borrajo, W. González-Manteiga and M.D. Martínez-Miranda where recall that 𝜆0(𝑥) =𝜆(𝑥)∕𝑚is the relative density of event locations, with 𝑚 =´𝑊𝜆(𝑥)𝑑𝑥, denoting the expected number of events in 𝑊. Let us now define the quadratic distance between the two theoretical functions of interest: ˆ 𝑊 (𝜆(𝑥)∕𝑚−𝝆(𝐙(𝑥))∕𝑚)2𝑑𝑥 =ˆ 𝑊(𝜆0(𝑥)−𝝆(𝐙(𝑥))∕𝑚)2𝑑𝑥. From this theoretical distance we construct our test statistic by plugging-in kernel estimates of the relative densities. For 𝜆0(𝑥)we use the two-dimensional kernel estimator given in (2), while for 𝝆(𝐙(𝑥)) we use the 𝑝-dimensional kernel estimator given in (5), and replace the unknown 𝑚by 𝑁. Hence, our test statistic is formally defined as follows: 𝑇𝑛=ˆ 𝑊( 𝜆0,𝐻 (𝑥)−  𝝆0,𝐵(𝐙(𝑥)))2𝑑𝑥, (6) where  𝜆0,𝐻 (𝑥)has been defined in (2)and 𝜌0,𝐵(𝑥)is defined from (5)as: 𝜌0,𝐵(𝐙(𝑥)) =  𝜆𝐵(𝑥) 𝑁1{𝑁≠0} =1 𝑁 𝑁 ∑ 𝑖=1 1 𝐠⋆(𝐙𝑖)𝐋𝐵(𝐳−𝐙𝑖)1{𝑁≠0}. Notice that 𝑇𝑛measures the discrepancy between the null and the alternative hypotheses. For big values of 𝑇𝑛we should reject the null hypothesis in favour of the alternative, i.e., the intensity of the Poisson process cannot be explained by the covariates information through model (3). On the other hand small values of 𝑇𝑛indicate that the intensity of the Poisson process may be explained by the covariates through such model. Remark 1. An important fact in the definition of this testing procedure is that it is based on a Poisson assumption. Hence, testing this condition is essential to obtain an accurate interpretation of the result of the test, i.e., when rejecting the null hypothesis, it is crucial to determine whether the rejection is due to the failure of the Poisson assumption or the falseness of 𝐻0itself. Hence, in practice, we use the tools available in the R-package spatstat to perform the test described in Loosmore and Ford (2006) (dclf.test) and Ripley (1977, 1981)(mad.test) for the K and L functions to verify the Poisson assumption before conducting the goodness-of-fit test. 3.2. Asymptotic properties For spatial and spatio-temporal point processes, three asymptotic regimes can be formulated: increasing-domain, infill asymptotics, Ripley (1988), and nearly infill sampling, as it is known the combination of both of them. As the ideas in this paper are related to density estimation, the asymptotic regime considered is infill structure asymptotic framework (Diggle and Marron, 1988), which states that the expected number of events tends to infinity for a fixed bounded observation domain. For a better understanding, take into account that the mathematical background behind this asymptotic framework is a sequence of counting processes 𝑁≡𝑁𝑚, where 𝑚 =𝔼[𝑁𝑚], and it is this expected number of events the quantity that increases to infinity. For simplicity here we develop the asymptotic theory for just one covariate (𝑝 =1). The result can be directly transferred, with similar arguments, to 𝑝 >1, but at the expense of more complexity and cumbersome notation. Some notation and regularity conditions need to be introduced. We use the notation ◦for the convolution between two functions, 𝑡𝑟(⋅)for the trace of a matrix, 𝐷2for the Hessian operator and 𝑅(ℎ) =´ℎ2(𝑥)𝑑𝑥, for any function ℎ. 𝑊=ℝ2is required for the consistency of the estimator  𝜆0,𝐻 , see Fuentes-Santos et al. (2015). Additionally we assume: (A.1) ´ℝ𝐿(𝑧)𝑑𝑧 =1; ´ℝ𝑧𝐿(𝑧)𝑑𝑧 =0and 𝜇2(𝐿) ∶= ´ℝ𝑧2𝐿(𝑧)𝑑𝑧 <∞with 𝐿a univariate kernel function. (A.2) The bandwidth value 𝑏 ≡𝑏(𝑚), denoted in the general scenario 𝑝 >1as 𝐵, verifies the following lim𝑚→∞𝑏 =0 and lim𝑚→∞ 𝐴(𝑚) 𝑏=0, where 𝐴(𝑚) ∶= 𝔼 [1 𝑁1{𝑁≠0}]. (A.3) The bandwidth matrix 𝐻is symmetric and positive-definite, such that all entries of 𝐻tend to zero and 𝑚−1|𝐻|−1∕2 →0as 𝑚 increases. (A.4) 𝑲is a continuous, symmetric, square integrable bivariate density function such that ´ℝ2𝑢𝑢𝑇𝑲(𝑢)𝑑𝑢 =𝝁𝟐(𝑲)𝐼2with 𝝁𝟐(𝑲)<∞ and 𝐼2denoting the two-dimensional identity matrix. (A.5) 𝑍(𝑥)is a continuity point of 𝜌for all 𝑥 ∈𝑊. The following theorem provides the normal limiting distribution of the statistic 𝑇𝑛. The result is based on the central limit theorem for the integrated squared error of multivariate kernel density estimators given by Hall (1984), see Appendix for the proof. In the development of this proof it is important to note that by using linearization techniques proposed in Collomb (1976), condition (A.3) implies that 𝐴(𝑚)|𝐻|−1∕2 →0as 𝑚increases: 𝔼[1 𝑁1{𝑁≠0}]=𝔼[1{𝑁≠0} 𝑁]=𝔼[𝜙 𝜉]=𝔼[𝜙] 𝔼[𝜉]+𝑐+𝑐(𝜈)+(−1)𝜈𝜎1,𝜈 𝔼[𝜉]𝜈, Computational Statistics and Data Analysis 194 (2024) 107929 5 M.I. Borrajo, W. González-Manteiga and M.D. Martínez-Miranda where the second and subsequent addends are smaller than the ratio of expectations. Hence, we can write that, in the limit, 𝔼[1 𝑁1{𝑁≠0}]≃𝔼[𝜙] 𝔼[𝜉]= 𝔼[1{𝑁≠0}] 𝔼[𝑁]=ℙ(𝑁≠0) 𝑚=1−𝑒−𝑚 𝑚≃1 𝑚. And from this, the implication of (A.3) is straightforward. Also the Poisson assumption, implying second order independence, which is essential in the proof to compute the conditional expectations. Theorem 1. Under conditions (A.1) to (A.5), and under the null hypothesis, 𝐻0∶𝜆(𝑥)∕𝑚 =𝜌(𝑍(𝑥))∕𝑚for all 𝑥 ∈𝑊, 𝑇𝑛−𝜇𝑇𝑛 𝜎𝑇𝑛 𝑚→∞ ←←←←←←←←←←←←←←←←←←←←←←←→𝑁(0,1), with 𝜇𝑇𝑛=𝐴(𝑚)|𝐻|−1∕2𝑅(𝑲)+𝝁𝟐(𝑲)ˆ𝜆0(𝑥)𝑡𝑟(𝐻𝐷2𝜆0(𝑥))𝑑𝑥 +1 4𝝁𝟐 𝟐(𝑲)ˆ𝑡𝑟2(𝐻𝐷2𝜆0(𝑥))𝑑𝑥, 𝜎2 𝑇𝑛=𝐴(𝑚)|𝐻|−1∕2 ¨𝜆2 0(𝑥)𝜆0(𝑦)(𝑲◦𝑲)(𝐻−1∕2(𝑥−𝑦))𝑑𝑥𝑑𝑦 +2𝐴(𝑚)|𝐻|−1∕2𝑅(𝜆0)𝑅(𝑲). Remark 2. To extend this result to general covariate dimension 𝑝 >1a modification of condition (A.2) is needed: (A.2)′The bandwidth matrix 𝐵is symmetric and positive-definite, such that all entries of 𝐵tend to zero and 𝐴(𝑚)|𝐵|−1∕2 →0as 𝑚 increases. This would derive into modifying the details along the proof (Appendix A) by changing the real value bandwidths of the covariatebased estimator into matrices and determinants. However, this has no direct impact on the final result of the asymptotic distribution, because the higher order terms are those of the nonparametric estimator (2) which does not include covariates. In practice, the asymptotic distribution given in Theorem 1could be approximated estimating 𝑚by the observed sample size 𝑛, and 𝐴(𝑚)by 1∕𝑛, see Cucala (2006). The slow convergence rate of the asymptotic distribution, which is detailed in the standardize factor, i.e. 𝐴(𝑚)|𝐻|−1∕2, and which is slower than parametric rate, indicates that this normal distribution may not be the best way to calibrate the test in practice, specially for small patterns. For small expected sample sizes, the real distance of the statistic to normality may not be small enough to provide with accurate results. Our proposal is to use bootstrap methods to approximate this asymptotic distribution. The bootstrap procedure, detailed in the following section, needs the result stated on Theorem 1to formally prove its consistency, see González-Manteiga and Crujeiras (2013)for a extent discussion of different goodness-of-fit tests using bootstrap methods for calibration. The idea behind the proof is to use similar arguments to those of Theorem 1using the plug-in estimators, and taking into account that the second extra step is basically averaging the resulting values from a consistent smooth bootstrap. 3.3. Bootstrap approximation to the null distribution of the test statistic Here we describe smooth bootstrap procedures inspired in Cao (1993)and Cowling et al. (1996)to resample under the null hypothesis, using the Poisson assumption. In Section 6we discuss possible extensions. We consider the general multivariate case with 𝑝 ≥1. Under the null hypothesis the intensity can be written as a function of the 𝑝covariates, 𝜆(𝑥) =𝝆(𝐙(𝑥)). So let us consider a pilot kernel intensity estimate  𝜆𝑅(𝑥)as in (5), with a pilot bandwidth matrix 𝑅. With this pilot estimate and conditional on the observed pattern 𝐗1, … , 𝐗𝑁, we suggest first a simpler bootstrap algorithm for the calibration of the test based on simulation from the null hypothesis, and second a two-stage bootstrap algorithm motivated by the insights from Baddeley et al. (2017). Algorithm 1 Simple smooth bootstrap procedure to compute the critical value of the test statistic. 1: Draw the number of events, 𝑛∗, from a random variable 𝑁∗having a Poisson distribution with intensity ´𝑊 𝜆𝑅(𝑥)𝑑𝑥. 2: Draw a bootstrap sample, ∗={𝐗∗ 1, … , 𝐗∗ 𝑛∗}, by sampling randomly 𝑛∗times from the distribution with density proportional to  𝜆𝑅in (5). 3: Compute the bootstrapped test statistic, 𝑇∗ 𝑛∗, using the expression in (6), but applied to the bootstrap sample. 4: Repeat steps 1 to 3 𝑀times leading to bootstrapped test statistics, 𝑇∗1 𝑛∗, …, 𝑇∗𝑀 𝑛∗, and approximate the critical value of the test at level 𝛼by the [(1 −𝛼)𝑀]-th order statistic of these 𝑀values. Compute the p-value as 𝑝 = 1+∑𝑀 𝑖=1 1{𝑇𝑛≤𝑇∗𝑖 𝑛∗} 𝑀+1 , where 𝑇𝑛is the value of the test statistics from the original sample defined in (6). Baddeley et al. (2017) discussed that classical Monte Carlo tests, as the derived from our algorithm above, are biased when the null hypothesis is composite. This can be overcome using a two-stage or double bootstrap. The motivation comes from the observation that the p-value derived from the algorithm is non-uniformly distributed when the null has to be estimated. Atwo-stage bootstrap will provide an adjustment of this p-value using second stage p-values (see also Davison and Hinkley, 1997). The full procedure is described below. Computational Statistics and Data Analysis 194 (2024) 107929 6 M.I. Borrajo, W. González-Manteiga and M.D. Martínez-Miranda Algorithm 2 Two-stage smooth bootstrap procedure to compute the critical value of the test statistic. 1: Run Algorithm 1to derive the first stage p-value 𝑝based on 𝑀bootstrap samples under the null hypothesis estimated with a pilot kernel intensity estimate  𝜆𝑅(𝑥). Keep the 𝑀bootstrap samples: ∗,1, … , ∗,𝑀 . 2: for 𝑖 =1, … , 𝑀with each sample ∗,𝑖 do 3: Estimate the intensity using the same kernel estimator as in the first stage, denote it by  𝜆∗,𝑖 𝑅. 4: Run Algorithm 1using  𝜆∗,𝑖 𝑅as the pilot intensity and derive the corresponding p-value 𝑝𝑖. 5: Compute the adjusted p-value as 𝑝𝑎𝑑𝑗 =1+∑𝑀 𝑖=1 1{𝑝𝑖≤𝑝} 𝑀+1 . The algorithm above admits some variations. First the bandwidth for the pilot intensity estimator could be chosen differently in the two stages. Second, the number of bootstrap samples simulated in the first and second stages could be also different, see discussion in Baddeley et al. (2017). Algorithm 2is computational much more demanding than Algorithm 1. For this reason in our simulations we restricted to the latter but computed both in the real data applications. Our simulations in Section 5will show that the bootstrap approximation of the null distribution derived from Algorithm 1performs quite well, providing a test with the correct empirical size, and reasonable power levels. In practice the two algorithms require the choice of several bandwidth parameters, including the pilot bandwidth mentioned before. Optimal bandwidth selection for tests based on kernel smoothing is still an open problem. Nevertheless in practical situations is often suitable to use good automatic bandwidth selectors designed for estimation, even though they might not be theoretically optimal for testing purposes. In point processes this has been considered for example in the former paper by Fuentes-Santos et al. (2015); in the density field for testing different characteristics, see for example Silverman (1981)and Mammen et al. (1992)for multimodality tests, Hyndman and Yao (2002)for symmetry tests and in the regression context see Eubank and Hart (1992)for goodness-of-fit. To derive our test we need two bandwidth matrices 𝐻, for the intensity estimator of Fuentes-Santos et al. (2015), and 𝐵for the covariate-based intensity estimator of Borrajo et al. (2020). Following the precedent literature we use the automatic data-driven bandwidth selectors proposed for the estimators in Fuentes-Santos et al. (2015)and Borrajo et al. (2020), respectively. In the numerical analyses showed later we will see that these choices seem to be appropriate, while for the pilot bandwidth 𝑅we have evaluated the test in a suitable range of values to validate our results. 4. Data analyses In this section we illustrate our proposal using two real data sets. The first one is the Murchison data set that has been used by A. Baddeley in his papers, books and software, see for example Baddeley et al. (2012)and Baddeley et al. (2015). The second data set consists of wildfires in Canada and it has been obtained from the Canadian Wildland Fire Information System http:// cwfis .cfs .nrcan .gc .ca /home. 4.1. Murchison gold deposits Murchison geological survey data shown in Fig. 1(left) record the spatial locations of 255 gold deposits and the surrounding geological faults. These data came from a 330 × 394 km region in the Murchison area of Western Australia and they have been obtained by Watkins and Hickman (1990). As it has been remarked in Baddeley et al. (2012), at scale 1:500000, the gold deposits spatial extension is negligible and they can be considered as points without losing generality. Note also that the real gold deposits and faults are three-dimensional, while here we use a two-dimensional projection. Moreover, some geological faults may have been missed because they are not recorded by direct observation but in magnetic field surveys or geologically inferred from discontinuities in the rock sequences. In Baddeley et al. (2012)the aim on studying these data is to “specify zones of high prospectivity to be explored for gold”, so they have already assumed that the influence of the fault information is relevant, and it may actually explain the localization of gold deposits under the one-dimensional model (1). From the faults it is possible to construct a covariate by computing the distance from every point in the observation region to the nearest fault, see Fig. 1(right). This covariate has been used to model the intensity of the process through (1). Our goal is to check whether this intensity model is indeed appropriate. As we have previously indicated, before applying our goodness-of-fit test we need to verify the inhomogeneous Poisson assumption. Hence, we have conducted the Loosmore and Ford (2006)and Ripley (1977, 1981)tests obtaining p-values of 0.12 and 0.09, respectively, and then assuming that this Murchison data set may come from an inhomogeneous Poisson point process. To apply the goodness-of-fit test described in the previous section we need to choose several bandwidth parameters. As previously discussed we propose to use automatic data-driven bandwidth selectors for the two bandwidths required to compute the test statistic in (6): the bandwidths 𝐵, which is a scalar in this case, and a two-dimensional matrix 𝐻. We use the bandwidth selectors proposed in Borrajo et al. (2020)based on a bootstrap procedure, and Fuentes-Santos et al. (2015), to estimate 𝐵and 𝐻, respectively. Additionally for the bootstrap method we need one more scalar bandwidth, this is, the pilot bandwidth 𝑅. In this case we use different values within a suitable range to thoughtfully validate our results. The range has been chosen around 0.52 (which is the value of the bandwidth selector of Borrajo et al., 2020) taking into account the scale of the data. Table 1shows the obtained p-values, considering 𝑀= 500 replications in the bootstrap algorithm. From these results we cannot reject the null hypothesis not having enough evidence against. All the bandwidths in the considered range provide high p-values. Our conclusion is therefore that, once it is assumed that the process fulfils a Poisson condition, there is statistical evidence supporting that Computational Statistics and Data Analysis 194 (2024) 107929 7 M.I. Borrajo, W. González-Manteiga and M.D. Martínez-Miranda Fig. 1. Geological faults (left) and covariate constructed from them (right) used in the Murchison data set. the geological faults are enough to explain the location of the gold deposits in the form shown in model (1). In this particular data set the Poisson assumption seems to be appropriate since the correlation in the data may be negligible, i.e., having a gold deposit at a certain location does not favour or prevent from having another gold deposit nearby. Table 1 P-values of the test for different pilot bandwidth values 𝑅in the bootstrap algorithm (Murchison data set). 𝑅0.30 0.40 0.50 0.60 0.70 p-value Algorithm 10.94 0.89 0.80 0.70 0.60 p-value Algorithm 20.80 0.87 0.72 0.49 0.38 4.2. Wildfires in Canada Forest fires are one of the most important natural disturbances since the last Ice Age and they represent a huge social and economic problem. Canada has a long tradition on recording information about their wildfires; and also studies from many different perspectives have been carried out: Walter et al. (2014), Rogers et al. (2013), Di Iorio et al. (2013), Flannigan and Harrington (1988). It is well known that fire activity in Canada mostly relies on meteorological elements such as long periods without rain and high temperatures. In this context we are interested in studying the influence of these meteorological variables on the spatial distribution of wildfires. Fig. 2. Wildfires in Canada during June 2015. The wildfire data set and also a complete meteorological information from the last decades is available at the Canadian Wildland Fire Information System website (http://cwfis .cfs .nrcan .gc .ca /home). The fire season in Canada lasts from late April until August, with a peak of activity in June and July, hence we analyse the influence of meteorological covariates on wildfires during June 2015, atotal number of 1841, see Fig. 2. In this study we focus on two covariates: temperature and precipitation, see Fig. 3left and right panels, respectively. The reason is that the above mentioned bibliography has determined them as the main necessary factor in the ignition of a wildfire. It is important Computational Statistics and Data Analysis 194 (2024) 107929 8 M.I. Borrajo, W. González-Manteiga and M.D. Martínez-Miranda to note that for inferential purposes we have removed two regions (Northwest Territories and Nunavut, mostly covered by ice layers) from the whole observation window (Canada). There are no fires registered on those regions and we cannot perform inference with such lack of information. We apply our proposed goodness-of-fit test to check the intensity dependence on the covariates in the form of model (3). We aim to conclude whether the temperature or/and the precipitation can explain the spatial distribution of the wildfires in Canada during June 2015. To this aim we check the goodness-of-fit of two intensity models. First the one-dimensional that uses the temperature as the only covariate, and then precipitation, in both cases 𝑝 =1. Then we tried a two-dimensional model using both temperature and precipitation together (𝑝 =2). The inhomogeneous Poisson assumption needs to be verified before conducting the goodness-of-fit test. We use once again the proposals in Loosmore and Ford (2006)and Ripley (1977, 1981)using in the onedimensional and the two-dimensional the corresponding intensity functions, with the effect of one or two covariates, respectively. In all the models, the tests do not reject the null hypothesis and then we can assume that the Poisson condition is fulfilled. Fig. 3. Covariates used in the intensity estimation of the wildfires in Canada during June 2015. Regarding to the first model, with just the temperature as covariate, we compute the test statistic using the same kind of bandwidth choices as for the Murchison data described before. Again we consider different bandwidth values in a suitable range for the pilot bandwidth in the bootstrap calibration. For these data it seems that a suitable range would be around 0.35, which is the value of the automatic bandwidth selector. In this case we get the same results for all the bandwidths in the range and both algorithms, rejecting the null hypothesis based on p-values that are always close to zero. Similar results are obtained when using precipitation as covariate. Our conclusion is that the temperature or precipitation cannot explain on its own the wildfires spatial distribution. We move then to the second model and include the two covariates, temperature and precipitation. In this two-dimensional context the bandwidth selectors of Borrajo et al. (2020)do not apply so we have used the plug-in rule for multivariate densities developed in Chacón and Duong (2010), which has been implemented in the function Hpi of the R-package ks. The result from our test is similar, and the null hypothesis is again rejected based on a p-value close to zero with both algorithms. This leads us to conclude that, even the covariates are likely to have an effect on the distribution of the wildfires in Canada, this is not enough to explain the first-order intensity by themselves in the way shown in (3). In this case the idea of being able to explain wildfires occurrences only with temperature and precipitation information may be too ambitious, and other meteorological covariates (or some indexes gathering several variables) should be used to analyse this process. Following this idea, we have used the Canadian Forest Fire Weather Index (FWI), which consists of six components that account for the effects of fuel moisture and weather conditions on fire behaviour. Calculation of the components is based on consecutive daily observations of temperature, relative humidity, wind speed, and 24-hour precipitation, for more details and information on this see https://cwfis .cfs .nrcan .gc .ca /background /summary /fwi. Using only this index as covariate is not neither enough to explain wildfire behaviour, possibly more complex models should be accounted for. 5. Simulation study This section is devoted to analyse the performance of our proposal through Monte Carlo simulations. Under a Poisson assumption we define two intensity models based on the two real data sets previously presented and analysed. To reduce the burden of computation we only consider a one-dimensional covariate (𝑝 =1). The first model is based on the wildfires data set and the theoretical intensity, represented in Fig. 4, has been constructed by computing the kernel intensity estimator of Borrajo et al. (2020)with the temperature as the covariate. The second model has been constructed in a similar way but using the Murchison data set, with the covariate being the distance to the nearest geological fault. The intensity function in this case is represented in Fig. 5. Fig. 4. Theoretical intensity function for the first model analysed in the simulation study, obtained by applying a kernel intensity estimator to the Canada wildfire data set. Computational Statistics and Data Analysis 194 (2024) 107929 9 M.I. Borrajo, W. González-Manteiga and M.D. Martínez-Miranda Fig. 5. Theoretical intensity function for the second model analysed in the simulation study, obtained by applying a kernel intensity estimator to the Murchison data set. To evaluate the empirical level of our test we generate samples under the null hypothesis. To evaluate the power of the test we simulate samples under the alternative. We define alternatives by including a multiplicative deviation from the null models as follows: 𝜆(𝑥)=𝜆null(𝑥)𝑡(𝑥). Here 𝜆null denotes the intensity under the null and 𝑡is a perturbation function that depends on a parameter that controls the discrepancy from the null. We have considered this function as a diagonal band that crosses the observation region nullifying the extension out of it, with a smooth change. The band is wider or thinner depending on a one-dimensional parameter: when its wide enough it covers the whole observation region, so we approach the null hypothesis, and as it becomes thinner the model gets away from it. Fig. 6. Representation of the 𝑡𝐶functions (first column) and the resulting intensity (second column) for the four values of the parameter 𝑑𝐶=6, 12, 20, 30 for the Canada model. Since the observation regions of each model are different we have defined two different functions but following the same idea. In both cases the diagonal band is based on a univariate normal density, 𝜙(⋅, 𝜇, 𝜎), depending on a parameter: 𝑑𝐶for the Canada data model and 𝑑𝑀for the Murchison data model. This gives alternatives for the first case with function 𝑡𝐶(𝑢, 𝑣) =𝜙(𝑢, 15 −𝑣 −𝑣0𝐶, 𝑑𝐶), where 𝑣0𝐶=60.40 is the middle point of the y-axis in the observation region, and 𝑑𝐶taking values 6, 12, 20 and 30. In Fig. 6we represent, for each value of the parameter 𝑑𝐶, both the 𝑡𝐶function (first column) and the final intensity function (second column). For the Murchison model, the function is defined as 𝑡𝑀(𝑢, 𝑣) =𝜙(𝑢 −𝑢0𝑀, 𝑢 −𝑣0𝑀, 𝑑𝑀), where 𝑢0𝑀= 517.69 and 𝑣0𝑀= 6900.61 are the middle point of the x and y-axis respectively, and 𝑑𝑀is the parameter taking values 10, 20, 40 and 60. We have decided to rotate the diagonal band due to the distribution of the data and the covariate information, but it could have been done in the other way. Take also into account that due to this rotation, even a thin band collects a lot of information from the covariate which may cause smaller rejection proportions. Indeed, we have performed the same simulations with a rotated band and we obtained higher rejection proportions. In Fig. 7we represent the 𝑡𝑀functions (first row) as well as the final intensities (second row). Remark that we have not included the scales in Fig. 6and 7because the values of the intensity depend on the expected sample size, so they will change in each of the situations considered in the study. Computational Statistics and Data Analysis 194 (2024) 107929 16 M.I. Borrajo, W. González-Manteiga and M.D. Martínez-Miranda Ogata, Y., Zhuang, J., 2006. Space–time ETAS models and an improved extension. Tectonophysics 413 (1), 13–23. Reitzner, M., Schulte, M., 2013. Central limit theorems for 𝑢-statistics of Poisson point processes. Ann. Probab. 41 (6), 3879–3909. Ripley, B.D., 1977. Modelling spatial patterns (with discussion). J. R. Stat. Soc., Ser. B 39, 172–212. Ripley, B.D., 1981. Spatial Statistics. John Wiley and Sons. Ripley, B.D., 1988. Statistical Inference for Spatial Processes. Cambridge University Press. Rogers, B.M., Randerson, J.T., Bonan, G.B., 2013. High-latitude cooling associated with landscape changes from North American boreal forest fires. Biogeosciences 10 (2), 699–718. Schoenberg, F.P., 2005. Consistent parametric estimation of the intensity of a spatial–temporal point process. J. Stat. Plan. Inference 128 (1), 79–93. Schoenberg, F.P., 2011. Multidimensional residual analysis of point process models for earthquake occurrences. J. Am. Stat. Assoc. 98, 789–795. Silverman, B.W., 1981. Using kernel density estimates to investigate multimodality. J. R. Stat. Soc., Ser. B, Methodol. 43 (1), 97–99. Stoyan, D., Penttinen, A., 2000. Recent applications of point process methods in forestry statistics. Stat. Sci. 15 (1), 61–78. Thurman, A.L., Fu, R., Guan, Y., Zhu, J., 2015. Regularized estimating equations for model selection of clustered spatial point processes. Stat. Sin. 25 (1), 173–188. Van Lieshout, M., 2000. Markov Point Processes and Their Applications. World Scientific. Waagepetersen, R.P., 2007. An estimating function approach to inference for inhomogeneous Neyman–Scott processes. Biometrics 63 (1), 252–258. Walter, C., Freitas, S., Kraut, I., Rieger, D., Vogel, H., Vogel, B., 2014. Influence of 2010 Canadian forest fires on cloud formation on the regional scale. In: AGU Fall Meeting Abstracts. Watkins, K.P., Hickman, A.H., 1990. Geological Evolution and Mineralization of the Murchison Province, Western Australia, vol. 1. Department of Mines, Western Australia. Yue, Y.R., Loh, J.M., 2015. Variable selection for inhomogeneous spatial point process models. Can. J. Stat. 43 (2), 288–305.