# Globally intensity-reweighted estimators for K- and pair correlation functions

We introduce new estimators of the inhomogeneous K-function and the pair correlation function of a spatial point process as well as the cross K-function and the cross pair correlation function of a bivariate spatial point process under the assumption of second-order intensity-reweighted stationarity. These estimators rely on a 'global' normalization factor which depends on an aggregation of the intensity function, whilst the existing estimators depend 'locally' on the intensity function at the individual observed points. The advantages of our new global estimators over the existing local estimators are demonstrated by theoretical considerations and a simulation study.

• 1 publication
• 11 publications
• 9 publications
12/03/2020

### Second order semi-parametric inference for multivariate log Gaussian Cox processes

This paper introduces a new approach to inferring the second order prope...
12/05/2017

### Second-order and local characteristics of network intensity functions

The last decade has witnessed an increase of interest in the spatial ana...
08/04/2022

### A Hawkes model with CARMA(p,q) intensity

In this paper we introduce a new model named CARMA(p,q)-Hawkes process a...
11/18/2016

### An Empirical Study of Continuous Connectivity Degree Sequence Equivalents

In the present work we demonstrate the use of a parcellation free connec...
07/03/2020

### Prediction of Spatial Point Processes: Regularized Method with Out-of-Sample Guarantees

A spatial point process can be characterized by an intensity function wh...
01/15/2019

### Second-order variational equations for spatial point processes with a view to pair correlation function estimation

Second-order variational type equations for spatial point processes are ...
11/29/2021

### Mapping the intensity function of a non-stationary point process in unobserved areas

Seismic networks provide data that are used as basis both for public saf...

## 1 Introduction

Functional summary statistics like the nearest-neighbour-, the empty space-, and Ripley’s -function have a long history in statistics for spatial point processes (Møller & Waagepetersen, 2004; Illian et al., 2008; Chiu et al., 2013). For many years the theory of these functional summary statistics was confined to the case of stationary point processes with consequently constant intensity functions. The paper Baddeley, Møller & Waagepetersen (2000) was therefore a big step forward since it relaxed substantially the assumption of stationarity in case of the -function and the closely related pair correlation function.

Baddeley, Møller & Waagepetersen (2000) introduced the notion of second-order intensity-reweighted stationarity (soirs) for a spatial point process. When the pair correlation function exists for the point process, soirs is equivalent to that is translation invariant. However, the intensity function does not need to be constant which is a great improvement compared to assuming stationarity, see e.g. Møller & Waagepetersen (2007). When the point process is soirs, Baddeley, Møller & Waagepetersen (2000) introduced a generalization of Ripley’s -function, the so-called inhomogeneous -function which is based on the idea of intensity-reweighting the points of the spatial point process, and they discussed its estimation. The inhomogeneous -function has found applications in a very large number of applied papers and has also been generalized e.g. to the case of space-time point processes (Gabriel & Diggle, 2009) and to point processes on spheres (Lawrence et al., 2016; Møller & Rubak, 2016). Moreover, van Lieshout (2011) used the idea of intensity-reweighting to generalize the so-called -function to the case of inhomogeneous point processes.

A generic problem in spatial statistics, when just one realization of a spatial process is available, is to separate variation due to random interactions from variation due to a non-constant intensity or mean function. In general, if an informed choice of a parsimonious intensity function model is available for a point process, the intensity function can be estimated consistently. Consistent estimation of the inhomogeneous -function is then also possible when the consistent intensity function estimate is used to reweight the point process, see e.g. Waagepetersen & Guan (2009) in case of regression models for the intensity function. When a parsimonious model is not available, one may resort to non-parametric kernel estimation of the intensity function as considered initially in Baddeley, Møller & Waagepetersen (2000). However, kernel estimators are not consistent for the intensity function and they are strongly upwards biased when evaluated at the observed points. This implies strong bias of the resulting inhomogeneous -function estimators when the kernel estimators are plugged in for the true intensity.

In this paper, we introduce a new approach to non-parametric estimation of the (inhomogeneous) and -functions for a spatial point process, or of the cross -function and the cross pair correlation for a bivariate spatial point process, assuming soirs in both cases. This formalizes an approach that was used by Stone et al. (2017)

to estimate space-time cross pair correlation functions in live-cell single molecule localization microscopy experiments with spatially varying localization probabilities. In the univariate case, our new as well as the existing estimators are given by a sum over all distinct points

and from an observed point pattern. For the new estimators, each term in the sum depends on an aggregation of the intensity function through a ‘global’ normalization factor instead of depending ‘locally’ on the intensity function at and at as for the existing estimators (a similar remark applies in the bivariate case). Intuitively one may expect this to mitigate the problem of using biased kernel estimators of the intensity function in connection to non-parametric estimation of the -function or pair correlation function. Moreover, to reduce bias when using a non-parametric kernel estimator of , we propose a ‘leave-out’ modification of our

estimator. Our simulation study shows that our new globally intensity reweighted estimators are superior to the existing local estimators in terms of bias and estimation variance regardless of whether the intensity function is estimated parametrically or non-parametrically.

The remainder of the paper is organized as follows. Some background on spatial point processes and notational details are provided in Section 2. Section 3 introduces our global estimator for the -function or the cross -function, discusses modifications to account for isotropy, and compares with the existing local estimators. Section 4 is similar but for our new global estimator of the -function or cross pair correlation function. Section 5 describes sources of bias in the local and global estimators when kernel estimators are used, and modifications to reduce bias. In Section 6, the global and local estimators of and are compared in a simulation study. Possible extensions are discussed in Section 7. Finally, Section 8 contains some concluding remarks.

## 2 Preliminaries

We consider the usual setting for a spatial point process defined on the -dimensional Euclidean space , that is, is a random locally finite subset of . This means that the number of points from falling in , denoted , is almost surely finite for any bounded subset of . For further details we refer to Møller & Waagepetersen (2004). In our examples, .

For any integer , we say that has -th order intensity function if for any disjoint bounded Borel sets ,

 E{N(A1)⋯N(An)}=∫A1⋯∫Anρ(n)(x1,…,xn)dx1⋯dxn<∞.

Since the set has Lebesgue measure 0, the specification of on this set is not so important, but for a specific model there will typically be an obvious choice (e.g. if is continuous outside this set, it may be natural to extend it to be continuous on all of ). By the so-called standard proof we obtain the -th order Campbell’s formula (see e.g. Møller & Waagepetersen, 2004): for any Borel function ,

 E≠∑x1,…,xn∈Xk(x1,…,xn)=∫⋯∫k(x1,…,xn)ρ(n)(x1,…,xn)dx1⋯dxn,

which is finite if the left or right hand side is so. Here, over the summation sign means that are pairwise distinct.

Throughout this paper, we assume that has an intensity function and a translation invariant pair correlation function . This means that for all , and , where with a symmetric Borel function. If is constant we say that is (first-order) homogeneous. In particular, if is stationary, that is, the distribution of is invariant under translations in , then is constant and is translation invariant.

Following Baddeley, Møller & Waagepetersen (2000), the translation invariance of implies that is second-order intensity reweighted stationary (soirs) and the inhomogeneous -function (or just -function) is then given by

 K(t):=∫∥h∥≤tg0(h)dh,t≥0.

This is Ripley’s -function when is stationary.

Suppose and are locally finite point processes on such that has intensity function , , and has a translation invariant cross pair correlation function for all . That is, for bounded Borel sets and denoting the cardinality of , , we have

 E{N1(A1)N2(A2)}=∫A1∫A2ρ1(x1)ρ2(x2)c(x1−x2)dx1dx2.

Then the cross -function is defined by

 K12(t):=∫∥h∥≤tc(h)dh,t≥0.

In practice are observed within a bounded window , and we use the following notation. The translate of by is denoted . For a Borel set , denotes the indicator function which is 1 if and 0 otherwise. The Lebesgue measure of (or area of when ) is denoted , and is the usual Euclidean length of .

## 3 Global and local intensity-reweighted estimators for K-functions

### 3.1 The case of one spatial point process

Considering the setting in Section 2 for the spatial point process , we define

 (1)

Clearly, is symmetric, that is, for all . We assume that with probability 1, for all distinct . Then, for , we can define

 ^Kglobal(t):=≠∑x,y∈X∩W1[∥y−x∥≤t]γ(y−x). (2)

If whenever , then

is an unbiased estimator of

. This follows from the second-order Campbell’s formula:

 E^Kglobal(t) =∫∫1[x∈W,y∈W,∥y−x∥≤t]γ(y−x)ρ(x)ρ(y)g0(y−x)dxdy =∫∫1[x∈W∩W−h,∥h∥≤t]γ(h)ρ(x)ρ(x+h)g0(h)dxdh =∫∥h∥≤tγ(h)γ(h)g0(h)dh=K(t).

We call the global estimator since it contrasts with one of the estimators suggested in Baddeley, Møller & Waagepetersen (2000): assuming that almost surely for distinct ,

 ^Klocal(t):=≠∑x,y∈X∩W1[∥y−x∥≤t]ρ(x)ρ(y)|W∩Wy−x|, (3)

which we refer to as the local estimator. Note that is also an unbiased estimator of provided for . In the homogeneous case,

 γ(h)=ρ2|W∩W−h|,

whereby , and in the stationary case, these estimators coincide with the Ohser & Stoyan (1981) translation estimator.

In practice and hence must be replaced by estimates. Estimators of and and the bias of these estimators are discussed in Section 5.

#### 3.1.1 Modifications to account for isotropy

In addition to soirs, it is frequently assumed that the pair correlation function is isotropic meaning that for some Borel function . We benefit from this by integrating over the sphere: for , define

 γiso(r):=∫Sd−1γ(rs)dνd−1(s)/ςd, (4)

where denotes the -dimensional unit-sphere, is the -dimensional surface measure on , and is the surface area of the unit sphere . Thus is the mean value of when

is a uniformly distributed point on the

-dimensional sphere of radius and center at the origin.

Assuming that almost surely for distinct , this naturally leads to another global estimator for when the pair correlation function is isotropic, namely

 ^Kisoglobal(t):=≠∑x,y∈X∩W1[∥y−x∥≤t]γiso(∥y−x∥). (5)

That is unbiased follows from a similar derivation as for : for any such that whenever ,

 E^Kisoglobal(t) =∫∥h∥≤tγ(h)γiso(∥h∥)g0(h)dh =∫t0g1(r)rd−1∫Sd−1γ(rs)γiso(r)dsdr (6) =∫t0g1(r)ςdrd−1dr =∫∥h∥≤tg1(∥h∥)dh=K(t), (7)

where (6) and (7) employ changes of variables to and from polar coordinates, respectively.

When is homogeneous, (5) coincides with the Ohser & Stoyan (1981) isotropic estimator. A local estimator of this form can also be defined:

 ^Kisolocal(t):=∑x,y∈X∩W1[∥y−x∥≤t]ρ(x)ρ(y)aW(∥y−x∥), (8)

where

 aW(r)=∫Sd−1|W∩W−rs|ds/ςd (9)

is an isotropized edge correction factor, and where it is assumed that almost surely for distinct . The local estimator is unbiased when for .

#### 3.1.2 Comparison of local and global estimators

The global and local estimators (2) and (3) differ in the relative weighting of distinct points . Namely, weights pairs from low-density areas more strongly than those from high-density areas, whilst for , the weight only depends on the difference . Theoretical expressions for the variances of the global and local -function estimators are very complicated, not least when the intensity function is replaced by an estimate. This makes it difficult to make a general theoretical comparison of the estimators in terms of their variances. However, under some simplifying assumptions insight can be gained as explained in the following.

Consider a quadratic observation window of sidelength . Then is a disjoint union of quadrats each of sidelength . Assume that the intensity function is constant and equal to within each , with naturally estimated by for . For fixed and large , when is replaced by its estimator , we can now approximate the local estimator:

 ^Klocal(t)= ≠∑u,v∈X∩W1[∥u−v∥≤t]^ρ(u)^ρ(v)|W∩Wu−v|≈n2∑i=1≠∑u,v∈X∩Wi1[∥u−v∥≤t]^ρ2i|W∩Wu−v| ≈ n2∑l=1≠∑u,v∈X∩Wi1[∥u−v∥≤t]^ρ2i|Wi∩(Wi)u−v|n2=1n2n2∑i=1^Ki,local(t).

where is the local estimator based on . The first approximation above follows because contributions from and , , are negligible for fixed and large, and the second approximation is justified since for , and will tend to 1 as increases. Following similar steps, we obtain for the global estimator,

 ^Kglobal(t)≈n2∑i=1^Ki,local(t)^ρ2i∑n2l=1^ρ2l.

Suppose is a Poisson process. Note that is an equally weighted average of the , but since the are independent, the optimal weighted average is obtained with weights inversely proportional to the variances of the . For large , the variance of is well approximated by (Ripley, 1988; Lang & Marcon, 2013) and the optimal weights are thus proportional to . Our global estimator is obtained from the optimal weighted average by replacing the optimal weights by natural consistent estimates. Hence one may anticipate that the global estimator has smaller variance than the local estimator. In a small-scale simulation study this was indeed the case, and the global estimator with (random) weights proportional to even had slightly smaller variance than when the optimal fixed weights were used.

### 3.2 The case of two spatial point processes

For two spatial point processes and as in Section 2, we define the following global estimator for the cross -function: for ,

 ^K12,global(t):=∑x∈X1∩W,y∈X2∩W1[∥y−x∥≤t]γ12(y−x) (10)

where

 γ12(h):=∫W∩W−hρ1(u)ρ2(u+h)du

and it assumed that almost surely for and . It is straightforwardly verified that is unbiased for any such that whenever .

The corresponding local estimator is

 ^K12,local(t):=∑x∈X1∩W,y∈X2∩W1[∥y−x∥≤t]ρ1(x)ρ2(y)|W∩Wy−x|, (11)

assuming that almost surely for and . The local estimator is unbiased when for .

Interchanging and does not affect (10): when is defined as in (10) with replaced by

 γ21(h):=∫W∩W−hρ1(u+h)ρ2(u)du.

This follows since by a change of variable, is symmetric, .

When the cross pair correlation function is also isotropic, additional unbiased estimators of are readily obtained in the same way as for the one point process case. Thus, defining

 γiso12(r):=∫Sd−1γ12(rs)dνd−1(s)/ςd,r≥0, (12)

and assuming that almost surely for and , we define an isotropic global estimator by

 ^Kiso12,global(t):=≠∑x∈X1∩W,y∈X2∩W1[∥y−x∥≤t]γiso12(∥y−x∥). (13)

This is easily seen to be unbiased when for . Finally, the isotropic local estimator is

 ^Kiso12,local(t):=≠∑x∈X1∩W,y∈X2∩W1[∥y−x∥≤t]ρ1(x)ρ2(y)aW(∥y−x∥), (14)

with as defined in Section 3.1.1, and it becomes unbiased if for .

## 4 Global and local intensity-reweighted estimators for pair correlation functions

### 4.1 The case of one spatial point process

Considering again the setting in Section 2 for the spatial point process , this section introduces global and local estimators for the translation invariant pair correlation function given by . Note that it may be easier to interpret than , but non-parametric kernel estimation of involves the choice of a bandwidth.

Let be a (normalized) kernel with bandwidth , that is, for , where

is a probability density function. We assume that

has support centered in the origin and contained in for some ; e.g.  could be a standard -dimensional normal density truncated to (this choice is convenient when is rectangular with sides parallel to the usual axes in ). Then, for ,

 E≠∑x,y∈X∩W κb(h−(y−x)) =∫W∫Wκb(h−(y−x))ρ(x)ρ(y)g0(y−x)dxdy (15) =∫W{∫W−h−xκb(−z)ρ(x)ρ(x+h+z)g0(h+z)dz}dx ≈g0(h)∫Wρ(x){∫W−h−xκb(−z)ρ(x+h+z)dz}dx (16) ≈g0(h)γ(h) (17)

where is defined in (1). Here, (15) follows from the second-order Campbell’s formula, and the approximations (16) and (17) are obtained for a small bandwidth . When , the bounded support of shrinks to whereby the approximations (16) and (17) become exact; the former approximation is expected to be more accurate, but the latter is simpler to compute.

From (17) we conclude that can be estimated by the following global estimator,

 ^gglobal(h):=≠∑x,y∈X∩Wκb(h−(y−x))/γ(h),

provided . This contrasts with the local estimator

 ^glocal(h):=≠∑x,y∈X∩Wκb(h−(y−x))/{ρ(x)ρ(y) |W∩Wx−y|},

which is analogous to the estimator suggested in Baddeley, Møller & Waagepetersen (2000) for an isotropic pair correlation function, see also Section 4.1.1.

#### 4.1.1 Modifications to account for isotropy

For isotropic point processes as defined in Section 3.1.1, the global pair correlation function estimator may be modified to estimate the isotropic pair correlation function given by : for such that , define

 ^gisoglobal(r):=1ςdrd−1≠∑x,y∈X∩W~κb(r−∥x−y∥)/γiso(r), (18)

where for , , , for a probability density with support centered at 0 and contained in the interval for some constant , and where is defined in (4). This definition is motivated by the following derivation:

 E≠∑x,y∈X∩W~κb(r−∥y−x∥) =∫W∫W~κb(r−∥y−x∥)ρ(x)ρ(y)g1(∥y−x∥)dydx (19) =∫W{∫∞0~κb(r−ξ)g1(ξ)ξd−1∫Sd−1ρ(x)ρ(x+ξs)1[x+ξs∈W]dνd−1(s)dξ}dx (20) ≈g1(r)ςdγiso(r)rd−1∫∞0~κb(r−ξ)dξ (21) ≈g1(r)ςdγiso(r)rd−1, (22)

using the second-order Cambell formula in (19), a ‘shift to polar coordinates’ in (20), the assumption that is small in (21), and that the kernel is a probability density function in (22). Note regarding (22) that

 ∫∞0~κb(r−ξ)dξ=∫r−∞~κb(ξ)dξ

which is not in general. Since for , the integral is 1 if . From (22) we obtain (18).

In the isotropic case the most commonly used local estimators (Baddeley, Møller & Waagepetersen, 2000) are

 ^gisolocal(r)=1ςdrd−1≠∑x,y∈X∩W~κb(r−∥y−x∥)ρ(x)ρ(y)|W∩Wx−y|

and

 ~gisolocal(r)=1ςd≠∑x,y∈X∩W~κb(r−∥y−x∥)ρ(x)ρ(y)|W∩Wx−y|∥y−x∥d−1,

assuming that almost surely for distinct . These estimators suffer from strong positive respectively negative bias for values of close to 0.

### 4.2 Two point processes

A similar derivation is possible for the cross pair correlation function of a bivariate point process , yielding similar global and local estimators of : for ,

 ^cglobal(h):=∑x∈X1∩W,y∈X2∩Wκb(h−(y−x))/γ12(h);

for ,

 ^cisoglobal(r):=1ςdrd−1∑x∈X1∩W,y∈X2∩W~κb(r−∥y−x∥)/γiso12(r);

and for almost surely when and ,

 ^clocal(h)=∑x∈X1∩W,y∈X2∩Wκb(h−(y−x))/{ρ1(x)ρ2(y) |W∩Wx−y|}

and

 ^cisolocal(r)=1ςdrd−1∑x∈X1∩W,y∈X2∩W~κb(r−∥y−x∥)ρ1(x)ρ2(y)|W∩Wx−y|.

Also an intermediate estimator is possible, with the intensity weighting for one of the processes applied locally, and the other applied globally: with , , and as above, we have

 E∑x∈X1∩W,y∈X2∩W κb(h−(y−x))ρ2(y) =∫W∫Wκb(h−(y−x))c(y−x)ρ1(x)dxdy =∫W∫W−x−hκb(−z)c(h+z)ρ1(x)dzdx ≈c(h)∫W∩W−hρ1(x)dx

for a small bandwidth , which suggests the partially-reweighted estimator

 ^cpartial(h):=∑x∈X1∩W,y∈X2∩Wκb(h−(y−x))ρ2(y)∫W∩W−hρ1(x)dx,

provided . This estimator may be useful when is much easier to estimate than , e.g. when is homogeneous.

## 5 Sources of bias when ρ is estimated

All of the estimators of , , , and discussed above are unbiased (at least when are sufficiently small) when the true intensity function is used to compute the weight functions in the local estimators or , , , or in the global estimators. However, in most applications is not known, and must be replaced by an estimate. When the source of inhomogeneity is well understood, it is recommended to fit a model with an appropriate parametric intensity function and use it as the estimate, cf. Baddeley, Møller & Waagepetersen (2000) and Waagepetersen & Guan (2009).

In the absence of such a model, the most common alternative is a kernel estimator

 ^ρ(x):=∑y∈X∩Wκσ(y−x)wW(x;y) (23)

where is a symmetric kernel on with bandwidth , and where is an appropriate edge correction weight. We take the standard choice from Diggle (1985),

 wW(x;y)=∫Wκσ(u−x)du,

see also Van Lieshout (2012) (other types of edge corrections may depend on both and which is why we write although the weight here only depends on .)

In the following we discuss estimators for and with particular focus on the implications of estimation bias when kernel estimators are used to replace the true or in the global and local estimators.

### 5.1 Bias of local estimators with estimated ρ

We start by considering a single spatial point process . For each point pair (), the corresponding term in the local - and pair correlation function estimators is normalized by the product . While an exact expression for the bias of the estimators with estimated is not analytically tractable, we can understand major sources of bias by considering the expression , which appears in each of the local estimators.

First, following Baddeley, Møller & Waagepetersen (2000), we note that as defined in (23) is subject to bias when evaluated at the points of , and that a ‘leave-one-out’ kernel estimator given by

 ¯ρ(x):=∑y∈(X∩W)∖{x}κσ(y−x)wW(x;y),x∈W, (24)

is a better choice, with reduced bias in most cases.

Second, we note that

 E(1/¯ρ(x))>1/E(¯ρ(x))

(if exists; in some cases it may be infinite). This follows from Jensen’s inequality, since is strictly convex for . In addition, note that the leading contribution to is proportional to (Liao & Berg, 2019). This discrepancy leads to a strong positive bias of the local - and pair correlation function estimators, especially at large , where and are almost independent. This effect becomes more pronounced for smaller , since typically increases as decreases.

Third, we note that for distinct points that are close compared to the bandwidth , the covariance of and leads to bias. For the local (and global) estimators, we consider sums over distinct , which leads us to condition on as follows (for details, see Coeurjolly, Møller & Waagepetersen, 2017). By conditioned on distinct points with , we mean that is equal to in distribution, where follows the second-order reduced Palm distribution of at :

 P(X∈F∣x,y∈X)=P(Xxy∪{x,y}∈F).

Assuming has -th order joint intensity functions for , has intensity function and second order joint intensity function . Now, for distinct with , neglecting the edge correction in (24) for simplicity, we obtain the following by the first and second-order Campbell’s formulas for and using that is symmetric:

 E [¯ρ(x)¯ρ(y)∣∣x,y∈X]=% E⎧⎨⎩∑u∈(Xxy∩W)∪{y}κσ(x−u)∑v∈(Xxy∩W)∪{x}κσ(y−v)⎫⎬⎭ =E≠∑u,v∈Xxy∩Wκσ(x−u)κσ(y−u)+E∑u∈Xxy∩Wκσ(x−u)κσ(y−u) +κσ(x−y)κσ(y−x) (25) +κσ(x−y)E∑v∈Xxy∩Wκσ(y−v)+κσ(y−x)E∑u∈Xxy∩Wκσ(x−u) =∫W∫Wκσ(x−u)κσ(y−v)ρ(4)(x,y,u,v)ρ(2)(x,y)dudv (26) +∫Wκσ(x−u)κσ(y−u)ρ(3)(x,y,u)ρ(2)(x,y