Estimating the intensity function of spatial point processes outside the observation window

Abstract

Mapping the intensity of objects, as animal or plant species in ecological studies, is cumbersome as soon as these objects are not accessible by automated methods. The knowledge at large scale of the underlying process variability can then only be obtained through sampling and spatial prediction. Here, we aim to predict the intensity of a point process, at locations where it has not been observed, conditional to the observation using the best linear unbiased combination of the point process realization in the observation window. We show that the weight function associated to the predictor is the solution of a Fredholm equation of second kind. Both the kernel and the source term of the Fredholm equation are related to the second-order characteristics of the point process through the pair correlation function. We propose here several approximations to solve the Fredholm equation in order to obtain practical solutions and restrict the solution space to that generated by linear combinations of $\left(i\right)$ step functions, which lead to a direct solution and $\left(ii\right)$ elementary functions of a finite element basis, which provide a continuous approximation. Results are illustrated on simulations and to predict the intensity of Black Locust in a region of France.

Introduction

In many applications the study window is too large to extensively map local intensity variations of the point process of interest since observation methods may be available at a much smaller scale only. That is for instance the case when studying the spatial repartition of a bird species at a regional scale, while the observations are made in windows of few hectares; or when detecting disease at the field scale, while observations correspond to spots of a few meter squares; or when mapping the presence of plant species at the catchment scale, while the observation scale is the meter square. The intensity must then be estimated from data issued out of samples spread in the study window, and hence, from a partial realization of the point process in this window.

We thus want to predict the intensity of a stationary point process conditional to its realization within the observation window $W$ at any point $x_o⁄\in W$. In the sequel, this conditional intensity is called local intensity (Gabriel et al., 2016). It allows us, through the conditioning, to take into account the second-order structure in the prediction. As an example, let us consider the Thomas process which is a Poisson cluster process where the cluster centers (parents) are assumed to be Poisson and the offsprings are normally distributed around their parent points. This process is stationary. When the boundary of $W$ splits a cluster, the local intensity across the boundary should be larger than the global intensity, i.e. the mean number of points per unit area. Another example is a soft core process, which is a pairwise interaction process. In this case, when an observed data point is located close to the boundary of $W$, the local intensity should be smaller than the global intensity as fewer points are expected around this point due to an interpoint interaction which decreases with distance.

Few ways exist to predict the local intensity. A first way consists of using (Tscheschel and Stoyan, 2006)’s reconstruction method based on the first- and second-order characteristics of the point process. Once the empirical point pattern predicted within a given window, one can get the intensity by kernel smoothing. As it is a simulation-based method, it requires long computation times, especially when the prediction window is large and/or the point process is complex. In alternative methods, the intensity of the point process is driven by a stationary random field, and whilst based on different concepts, Bayesian Diggle and Ribeiro (2007), Diggle et al. (2013) or geostatistics Monestiez et al. (2006), Bellier et al. (2013), they are constrained within the class of Cox processes. van Lieshout and Baddeley (2001) developed for a wider class of parametric models a Bayesian approach for extrapolating and interpolating clustered point patterns. Saito et al. (2005) compared geostatistical methods for mapping object counts collected from strip transects. Although no model is then required, the direct application of kriging methods does not allow to make good prediction for repulsive processes or processes with peaky local variations. Indeed, variograms are directly computed on count data and do not take explicitly into account the structure of the point process through the pair correlation function for instance. In what follows, we propose to predict the local intensity from the first- and second-order characteristics of the point process, which we can classically estimate in practice and which allow to handle a large scope of models. A first version of such an approach has been developed in Gabriel et al. (2016), where the point process is regularized to get a count process over a grid and where the ordinary kriging is then adapted, with kriging weights defined from the structure of the point process (intensity and pair correlation function). This empirical approach thus mimics geostatistics through the link between the variogram and the pair correlation function. Whilst it can be applied in practice, it does not allow us to understand how the weights are built. The continuous approach developed in this paper generalizes the previous one and offers new (and better) approximations.

The predictor of the local intensity is defined in Section 2 as a linear combination of the point process realization. We get an exact, but not explicit, solution of the related weights, given by an integral equation (second kind Fredholm equation). Our approach clearly brings out how the pair correlation function operates and then why the proposition given in Gabriel et al. (2016) is a good approximation of our new theoretical solution. It also opens to other approximations, in particular we present here a finite element approach. Finally in Section 3, we illustrate the method on simulated and real data, including mixing processes with different interaction structures, at different scales.