Nonparametric Testing of the Covariate Significance for Spatial Point Patterns under the Presence of Nuisance Covariates

Abstract

Determining the relevant spatial covariates is one of the most important problems in the analysis of point patterns. Parametric methods may lead to incorrect conclusions, especially when the model of interactions between points is wrong. Therefore, we propose a fully nonparametric approach to testing significance of a covariate, taking into account the possible effects of nuisance covariates. Our tests match the nominal significance level, and their powers are comparable with the powers of parametric tests in cases where both the model for intensity function and the model for interactions are correct. When the parametric model for the intensity function is wrong, our tests achieve higher powers. The proposed methods rely on Monte Carlo testing and take advantage of the newly introduced concepts: the covariate-weighted residual measure and nonparametric residuals. We also define a correlation coefficient between a point process and a covariate and a partial correlation coefficient quantifying the dependence between a point process and a covariate of interest while removing the influence of nuisance covariates. Supplementary materials for this article are available online.

KEYWORDS:

• Nonparametric methods
• Partial correlation coefficient
• Point process
• Random shift test
• Residual analysis

1 Introduction

1.1 Motivation and Overview

Spatial point patterns are often accompanied by spatial covariates. Determining the relevant covariates that influence the positions of points is certainly one of the most important questions of point pattern analysis. Applications include spatial epidemiology, spatial ecology, exploration geology, seismology, and many other fields.

In this article, we mainly focus on this question. Our proposed methods use nonparametric tools. The second question that we are interested in is nonparametric quantification of the spatial dependence between a point process and a covariate, both without and with the presence of nuisance covariates. Nuisance covariates are those which potentially affect the distribution of the point process but which are not of particular interest in the analysis. We define a correlation coefficient and a partial correlation coefficient (removing the influence of the nuisance covariates) between a point process and a covariate. The second problem has not been studied before, to our knowledge.

The first problem is usually solved by parametric methods (Kutoyants 1998; Schoenberg 2005; Waagepetersen and Guan 2009; Coeurjolly and Lavancier 2013), see Section 2.1 for details. However, we show in our simulation study that even when the parametric model is selected correctly, these tests of covariate significance may lead to liberality. The parametric methods have even bigger problems when: (a) the parametric model for the intensity function is incorrect, or (b) the form of interactions between points is specified incorrectly. We propose here two tests of covariate significance, a fully nonparametric one which avoids both selecting the intensity function model and the interaction model, and a semiparametric one which does not assume an interaction model but uses the log-linear intensity function model as the one predominantly used in practice. These two proposed tests do not exhibit liberality, and their powers are comparable with the powers of parametric methods in cases with correctly specified models for the intensity function and the interactions. The proposed tests also have a higher power than the parametric ones when either the intensity function model or the interaction model is misspecified.

Since the proposed nonparametric tests do not need to choose a specific model and exhibit better properties than their parametric counterparts, their use should become a standard practice in the analysis of point patterns.

For determining relevant covariates one can also use the lurking variable plots (Baddeley and Turner 2005) or appropriate information criteria (Choiruddin, Coeurjolly, and Waagepetersen 2021) but these do not provide formal tests. The only nonparametric method studying the dependence of a point process and a covariate without nuisance covariates was introduced in Dvořák et al. (2022).

We assume here that the spatial covariates are continuous. The methodology is up to a certain extent also applicable for categorical covariates, as discussed in Section 6.

1.2 Motivational Example

To illustrate the relevance of the questions posed above, we consider a part of the tropical tree dataset from the Barro Colorado Island plot (Condit 1998). We focus on the positions of 3604 trees of the Beilschmiedia pendula species in a rectangular 1000 × 500 meter sampling plot (see the top left panel of Figure 1). This part of the dataset is available in the R package spatstat (Baddeley, Rubak, and Turner 2015). Below, we call it the BCI dataset.

Fig. 1 The Barro Colorado Island dataset. From left to right, top to bottom: locations of trees, terrain elevation, terrain gradient, the soil contents of nitrogen, phosphorus, and potassium.

The intensity of point occurrence in the observation window is clearly nonconstant as the trees tend to prefer specific environmental conditions. The variation in the intensity of point occurrence may possibly be explained by the accompanying covariate information. The available covariates include the terrain elevation and gradient (available in the spatstat package) and the soil contents of mineralized nitrogen, phosphorus and potassium (Dalling et al. 2022), see Figure 1. Maybe all the covariates bring important information and should be used for inference. However, it is equally possible that some of the covariates bring redundant information (as could be expected from the nitrogen and potassium content in this dataset, see the bottom left and bottom right panel of Figure 1) or that some of the covariates, in fact, do not influence the point process. It is important to determine with high degree of confidence which covariates influence the point process and should be included in the further steps of the inference.

In certain cases, a relevant parametric model can be specified based on the available expert knowledge. However, often no such parametric model is available, or we do not want to take the risk of model misspecification. Then nonparametric methods for covariate selection need to be used. For example, the nonparametric tools proposed in this article can be used to create a backward covariate selection scheme, as described in Section 5.

1.3 Outline of the Work

In order to achieve our objectives, we propose to employ the residual analysis (Baddeley et al. 2005) with respect to the model built from the nuisance covariates. The sample (Kendall’s) correlation coefficient between the smoothed residual field and the covariate of interest then quantifies their dependence, either without or with the effect of nuisance covariates. The latter defines the partial correlation coefficient.

The testing of covariate significance is proposed to be performed via a new test statistic, the covariate-weighted residual measure, and a Monte Carlo test. The residual analysis can be performed in a parametric way, which defines our semiparametric approach, or it can be performed nonparametrically using the nonparametric estimate of the point process intensity function (Baddeley et al. 2012), which leads to our completely nonparametric approach. The nonparametric residuals are used for the first time in this work.

The replications in the Monte Carlo test are obtained through random shifts both with torus correction (Lotwick and Silverman 1982) and variance correction (Mrkvička et al. 2021). The torus correction is a standard method whereas the variance correction was recently defined, and it allows to use nonrectangular windows and better controls the significance level of the test than the torus correction.

The article is organized as follows. Section 2 recalls the concepts we need to define our procedures. Section 3 describes all new methods we are introducing in this work. That is, nonparametric residuals, spatial (partial) correlation coefficient, covariate-weighted residual measure, and tests of covariate significance with nuisance covariates. Section 4 contains a simulation study in which the exactness and power of our nonparametric methods is compared with parametric methods. Section 5 contains an example of the usage of our methods for nonparametric selection of relevant covariates. Finally, Section 6 is left for conclusions and discussion. The supplementary material contains further simulation results, discussion on the choice of the variance correction factors in the random shift tests, and a data example illustrating the use of the (partial) correlation coefficients between a point process and a covariate. The R codes providing an implementation of the proposed methods are available in the package NTSS for R.

2 Notation and Background

Let X be a point process on ${\mathbb{R}}^2$ with the intensity function $\lambda (u)$. Throughout this article, we assume that the intensity function of X exists. Let $C_1,C_2,\dots ,C_m$ be the covariates in ${\mathbb{R}}^2$. These will later represent the nuisance covariates whose influence on the point process is to be removed. The covariate of interest will be denoted $C_{m+1}$ in the following. Denote by $W\subset {\mathbb{R}}^2$a compact observation window with area $\left|W\right|$ and $n(X\cap B)$ the number of points of the process X observed in the set B. We assume that the values of the covariates are available in all points of W, at least on a fine pixel grid. This can be achieved from a finite set of observations for example by kriging techniques.

2.1 Covariate Selection in Parametric Point Process Models

The dependence of the intensity function of a point process on the covariates $C_1,\dots ,C_m$ is often modeled parametrically, for example using the log-linear model (1) $$\lambda (u;\beta )=\hspace{0.17em}\mathrm{\text{exp}}\hspace{0.17em}\{{\beta }_0+{\beta }_1{C}_1(u)+\cdots +{\beta }_m{C}_m(u)\}.$$(1)

The standard approach to estimating the model parameters β_i is to maximize the Poisson likelihood (Schoenberg 2005; Waagepetersen and Guan 2009). This corresponds to the maximum likelihood approach for Poisson models, while for non-Poisson models, this constitutes a first-order composite likelihood approach. For the log-linear model (1) the estimation is implemented in the ppm function from the popular spatstat package.

For Poisson or Gibbs processes, the ppm function also provides confidence intervals for the regression parameters β_i and the p-values of the tests of the null hypothesis that ${\beta }_i=0$ for a given i, based on the asymptotic variance matrix (Kutoyants 1998; Coeurjolly and Rubak 2013). For cluster processes, the kppm function from the spatstat package provides means of model fitting. The regression parameters β_i from (1) are again estimated using the ppm function, but the asymptotic variance matrix is determined according to Waagepetersen (2008), taking into account the attractive interactions between points.

The methods discussed above provide means for formal testing of the hypothesis that ${\beta }_i=0$ for a given $i\in \{1,\dots ,m\}$, allowing one to select the set of relevant covariates to be included in the model.

2.2 Parametric Residuals for Point Processes

Residuals can be used to check whether the fitted model for the intensity function is appropriate, see Baddeley et al. (2005) or Baddeley, Rubak, and Turner (2015, sec. 11.3). In the following we employ the version of residuals based on the intensity function, as suggested by R. Waagepetersen in the discussion to the paper Baddeley et al. (2005), rather than based on the conditional intensity function as discussed in the paper itself. Let $\widehat{\beta }$ be the vector of the estimated regression parameters. The residual measure is defined as (2) $$\mathcal{R}(B)=n(X\cap B)-{\int }_B\lambda (u;\widehat{\beta })\mathrm{d}u,$$(2) where $B\subseteq W$ is a Borel set. The smoothed residual field is obtained as (3) $$s(u)=\frac{1}{e(u)}\left[\mathrm{}\sum _{x_i\in X\cap W}k(u-x_i)-{\int }_{W}k(u-v)\lambda (v;\widehat{\beta })\mathrm{d}v\mathrm{}\right]\mathrm{},$$(3) where $e(u)={\int }_{W}k(u-v)\mathrm{d}v$ is the edge-correction factor and k is a probability density function in ${\mathbb{R}}^2$. In fact, the first term in (3) gives the nonparametric kernel estimate of the intensity function, the covariates not being taken into account, while the second term gives the smoothed parametric estimate which incorporates the covariates. If the estimated model $\lambda (v;\widehat{\beta })$ describes the point process X well, the smoothed residual field s(u) is expected to fluctuate around 0 (the heuristics behind this is described in more detail in Section 3.4). Its deviations from 0 indicate a disagreement between $\lambda (v;\widehat{\beta })$ and the true intensity function in the corresponding parts of the observation window. We remark that the residuals described above are the raw residuals of Baddeley et al. (2005), where scaled versions of the residuals are also considered.

2.3 Nonparametric Estimation of the Intensity Function Depending on Covariates

As opposed to fitting a parametric model such as (1), the dependence of the intensity function on a set of covariates can be captured nonparametrically. Baddeley et al. (2012) assume that there is an unknown function $\rho :{\mathbb{R}}^m\to [0,\infty )$ such that $\lambda (u)=\rho (C_1(u),\dots ,C_m(u)).$ Assuming absolute continuity of the distribution of the vector of covariates $(C_1(u),\dots ,C_m(u))$ on ${\mathbb{R}}^m$, the function ρ can be estimated using kernel smoothing in the space of covariate values, see Baddeley et al. (2012) or Baddeley, Rubak, and Turner (2015, sec. 6.6.3). This opens up the possibility to define the nonparametric residuals in Section 3.1.

The estimation of ρ is implemented in the rhohat function from the spatstat package for m = 1 and in the rho2hat function for m = 2. We note that in these two cases, visualization of $\widehat{\rho }$ is straightforward while it is not as easy for m > 2. In our simulation experiments in Section 4 we use the spatstat implementation, while in the analysis of the real dataset with higher number of covariates we use our implementation based on the ks package (Duong 2007).

2.4 Monte Carlo Tests

When the distribution of a test statistic is too complicated to be derived analytically but there is a way of obtaining replications (simulations, permutations, …) of the data under the null hypothesis, it is possible to perform a formal test of the null hypothesis using the Monte Carlo approach (Davison and Hinkley 1997). This approach relies on the exchangeability of the vector $(T_0,T_1,\dots ,T_N)$, where T₀ is the test statistic value computed from the observed data, and $T_1,\dots ,T_N$ are obtained from the replications.

The test is performed by determining how typical or extreme the value T₀ is in the whole sample $T_0,T_1,\dots ,T_N$. For univariate test statistics, this means determining the rank of T₀, however, using functional test statistics is also possible if a suitable ranking of the functions from the most typical to the most extreme is available, as for example in Myllymäki et al. (2017). Exchangeability (invariance of the distribution with respect to permutations of the components) ensures that the Monte Carlo test matches the required significance level.

2.5 Random Shift Permutation Strategy

Random shifts provide means of nonparametric testing of independence between a pair of spatial objects, such as a pair of random fields (Upton and Fingleton 1985; Dale and Fortin 2002) or a pair of point processes (Lotwick and Silverman 1982). By randomly shifting one of the objects while keeping the other one fixed, any possible cross-correlation between them is broken while keeping the distribution of the marginals intact. Usually, at least one of the spatial objects is assumed to be stationary (see Cronie and van Lieshout (2016) for a specific example of a random shift test where both spatial objects may be nonstationary). By performing a certain amount of shifts along randomly generated vectors, one obtains replications for performing a Monte Carlo test of independence.

Assume that the spatial objects are denoted by $\Phi$ and $\Psi$ and we observe them in the window W. We denote the value of the test statistic computed directly from the observed data by $T_0=T(\Phi ,\Psi ;W)$. After producing N random shift vectors $v_1,\dots ,v_N$ we compute the value of the test statistic T_i from $\Phi$ and $\Psi$ shifted by v_i, that is $T_i=T(\Phi ,\Psi +v_i;W),i=1,\dots ,N$. Clearly, some part of $\Psi$ will be shifted outside of the observation window W and part of $\Psi +v_i$ will not overlap with $\Phi$ anymore. Hence, some form of correction is needed.

2.5.1 Torus Correction

For a rectangular window W, one may identify its opposing edges, creating a toroidal geometry on W (Lotwick and Silverman 1982; Upton and Fingleton 1985). We denote by $[\Psi +v_i]$ the version of $\Psi$ shifted with respect to the toroidal geometry, as opposed to $\Psi +v_i$ which denotes $\Psi$ shifted with respect to the Euclidean geometry. The replications T_i are then obtained as $T_i=T(\Phi ,[\Psi +v_i];W),i=1,\dots ,N$.

As a result, all parts of the data are used for computing T_i. On the other hand, artificial cracks appear in the autocorrelation structure of the shifted marginal, as parts of the data originally far away are now “glued together”. This means that exchangeability is violated, which in turn introduces liberality of the random shift tests (Fortin and Payette 2002; Mrkvička et al. 2021). However, simulation studies show that when the spatial autocorrelations in the data are not very strong, the tests match the nominal significance level quite closely (Mrkvička et al. 2021; Dvořák et al. 2022). Traditionally, the distribution of the random shift vectors is taken to be the uniform distribution on W, but other choices are also possible.

2.5.2 Variance Correction

To remove the liberality of the torus correction, Mrkvička et al. (2021) proposed the variance correction. It uses shifts respecting the Euclidean geometry and discards those parts of the data that are shifted outside of W. No artificial cracks are introduced to the correlation structure of the data, removing the liberality of the random shift tests. Also, irregular observation windows can be considered. On the other hand, different amounts of data are dropped for different shift vectors v_i and for typical choices of the test statistic the variance of T_i varies greatly, making it impossible to perform the Monte Carlo test directly. Therefore, the variance of T_i needs to be standardized before performing the test.

Formally, we denote by W_i the smaller observation window where $\Phi$ and $\Psi +v_i$ overlap, that is $W_i=W\cap (W+v_i)$. The value T_i is computed from $\Phi$ and $\Psi +v_i$ restricted to W_i, specifically as $T_i=T(\Phi {|}_{W_i},(\Psi +v_i){|}_{W_i};W_i)$. The values $T_0,T_1,\dots ,T_N$ are then standardized to have zero mean and unit variance. This is achieved by subtracting the mean $\overline{T}=\frac{1}{N+1}{\sum }_{i=0}^N{T}_i$ and dividing by the square root of the variance: $S_i=(T_i-\overline{T})/\sqrt{\text{var}(T_i)}.$ The standardized values $(S_0,S_1,\dots ,S_N)$ are closer to exchangeability than $(T_0,T_1,\dots ,T_N)$ because their first two moments are the same. The standardized values are used to perform the Monte Carlo test. When a formula describing $\text{var}(T_i)$ as a function of the size of W_i is known, at least asymptotically, it can be directly used in the standardization. If such a formula is not available, Mrkvička et al. (2021) suggest a kernel regression approach to estimating $\text{var}(T_i)$.

Simulation studies in Mrkvička et al. (2021) and Dvořák et al. (2022) show that the random shift tests with variance correction match the nominal significance level even in the case of strong autocorrelation. In those papers, the shift vectors followed the uniform distribution on a disk with radius R centered at the origin. The choice of R is a compromise between two goals: longer shifts are more relevant for breaking the possible dependence between $\Phi$ and $\Psi$ while shorter shifts mean that a larger amount of available data is used to compute T_i. Choosing R so that $\left|W_i\right|/\left|W\right|\ge 1/4$ for all i turned out to provide satisfactory results.

2.6 Nonparametric Testing of Dependence between Point Process and a Covariate

For nonparametric testing of the null hypothesis of independence between a point process X and a covariate C₁ the paper Dvořák et al. (2022) suggests to use the random shift test with the test statistic $T=\frac{1}{n(X\cap W)}\sum _{x_i\in X\cap W}{C}_1(x_i)$, that is the mean covariate value observed at the points of the process. This test showed liberality (with torus correction) or slight conservativeness (with variance correction) in the simulation studies in Dvořák et al. (2022), with both versions having much higher power than the other tests considered there.

3 New Methods

3.1 Nonparametric Residuals for Point Processes

As discussed in Section 2.3, a nonparametric estimate of the intensity function $\widehat{\lambda }(u)=\widehat{\rho }(C_1(u),\dots ,C_m(u))$ can be used to describe its dependence on the set of covariates. Using $\widehat{\rho }$, the nonparametric version of the residual measure (2) can be defined as (4) $$\tilde{\mathcal{R}}(B)=n(X\cap B)-{\int }_B\widehat{\rho }(C_1(u),\dots ,C_m(u))\mathrm{d}u.$$(4)

The corresponding nonparametric smoothed residual field is then (5) $$\begin{array}{c}\tilde{s}(u)=\frac{1}{e(u)}[\sum _{x_i\in X\cap W}k(u-x_i)\\ \mathrm{\hspace{1em}}\mathrm{\hspace{1em}}-{\int }_W\mathrm{}k(u-v)\widehat{\rho }(C_1(u),\dots ,C_m(u))\mathrm{d}v]\mathrm{}.\end{array}$$(5)

Again, scaled versions of these residuals can be constructed as in Baddeley et al. (2005). If $\widehat{\rho }(C_1(u),\dots ,C_m(u))$ describes the intensity function of X well, meaning for example that no relevant covariate was left out, $\tilde{s}(u)$ is expected to fluctuate around 0 (the heuristics behind this is described in more detail in Section 3.4). Figure 2 illustrates that $\widehat{\rho }$ is capable of capturing the correct form of dependence even without specifying a parametric model. It shows a realization of the Poisson process on ${[0,1]}^2$ with intensity function $\lambda (x,y)=400(1-4{(x-1/2)}^2)$, the nonparametric smoothed residual field $\tilde{s}$ from (5) depending on the covariate x, the parametric smoothed residual field s from (3) with log-linear model depending on x, and the parametric smoothed residual field s from (3) with log-linear model depending on x and x².

Fig. 2 Left to right: realization of a Poisson process with a quadratic intensity function, the nonparametric smoothed residual field, the parametric smoothed residual field with a log-linear model of the intensity function, the parametric smoothed residual field with a log-quadratic model of the intensity function.

3.2 Correlation Coefficient between a Point Process and a Covariate

Assume now that no nuisance covariates are given (m = 0) and we want to investigate the strength of dependence between the intensity function of X and a given covariate C₁. Without incorporating a possible effect of C₁, the natural estimate of the intensity function is constant, $\widehat{\lambda }=n(X\cap W)/\left|W\right|$, and the smoothed residual field becomes (6) $$\tilde{s}(u)=\frac{1}{e(u)}\sum _{x_i\in X\cap W}k(u-x_i)-\widehat{\lambda }.$$(6)

If the covariate C₁ does not influence X, we expect C₁ and $\tilde{s}$ to be independent. On the other hand, if C₁ influences the intensity function of X, $\tilde{s}$ should capture the dependence structure and exhibit correlations with C₁. This motivates us to quantify the strength of dependence between X and C₁ by some measure of dependence between the two random fields C₁ and $\tilde{s}$.

To this end, we consider Kendall’s correlation coefficient (Nelsen 2006, p.158) and let U₁, U₂ be independent random vectors with uniform distribution in W. Denoting $Y=C_1(U_1)-C_1(U_2)$ and $Z=\tilde{s}(U_1)-\tilde{s}(U_2)$, we define (7) $$\tau =\mathbb{P}(Y·Z>0)-\mathbb{P}(Y·Z<0).$$(7)

The empirical estimate of τ can be easily obtained if we consider a set of sampling points $\{y_1,\dots ,y_n\}$, independently and uniformly distributed in W, independent of X and C₁: (8) $$\widehat{\tau }=\frac{1}{n(n-1)}\sum _{i\ne j}\text{sgn}(C_1(y_i)-C_1(y_j))\text{ sgn}(\tilde{s}(y_i)-\tilde{s}(y_j)),$$(8) where $\text{sgn}$ is the sign function. Naturally, the values of the correlation coefficient are restricted to the interval $[-1,1]$ and allow direct comparison of the strength of dependence between different datasets.

To illustrate the use of this correlation coefficient in quantifying the strength of dependence between a point process and a covariate, we perform the following experiment. We consider the Poisson process with the intensity function proportional to $\hspace{0.17em}\mathrm{\text{exp}}\hspace{0.17em}\left\{\mathit{\text{ax}}\right\}$ in the observation window $W={[0,1]}^2$, for a given value of $a\in \mathbb{R}$, and with the expected number of points in W fixed at 200. The fixed expected number of points is achieved by adjusting the constant of proportionality in the intensity function. The covariate of interest is $C_1((x,y))=x$. The smoothed residual field $\tilde{s}$ from (6) is obtained with a large bandwidth bw = 0.5 which reflects the fact that the true intensity function of the point process is very smooth. The value of $\widehat{\tau }$ is then computed according to (8). This is repeated for 500 independent realizations of the point process for each value of a from a fine grid, and the means of $\widehat{\tau }$ are plotted as a function of a in the left panel of Figure 3. The plot shows that $\widehat{\tau }$ increases in the absolute value with increasing strength of dependence, from 0 in case of independence (a = 0) to almost 1 or –1 in case of very strong dependence. It also correctly captures the form of dependence (positive or negative association).

Fig. 3 Left: plot of the correlation coefficient $\widehat{\tau }$ as a function of the parameter a for the example in Section 3.2. Right: plot of the correlation coefficient $\widehat{\tau }$ (black curve) and the partial correlation coefficient ${\widehat{\tau }}_p$ (red curve) as functions of the parameter a for the example in Section 3.3.

3.2.1 Choice of Sampling Points

We stress that independent sampling points need to be used in this case instead of simply using the observed points of $X\cap W$. In the latter case, the preferential sampling issues could arise, resulting in biased estimates of the properties of the two random fields (Diggle, Menezes, and Su 2010; Dvořák et al. 2022). Loosely speaking, if, for example, the sampling points $\{y_1,\dots ,y_n\}$ are more likely to be chosen in locations with high values of C₁, the sample mean and sample variance of $C_1(y_1),\dots ,C_1(y_n)$ do not reflect well the true properties of C₁. This negatively affects all subsequent steps of the analysis.

Choosing the sampling points independently and uniformly distributed in the observation window is the simplest option. Under such assumption, the variance correction factors for the random shift test can be determined, see the discussion in Appendix C in the supplementary material. On the other hand, if only quantification of the dependence is desired and testing is not intended, using a regular grid of sampling points placed randomly in the observation window may result in a slightly smaller variance of $\widehat{\tau }$ compared to independent sampling points. This was confirmed by simulation experiments not reported here.

3.2.2 Choice of Measure of Dependence

Although different measures of dependence such as Pearson’s or Spearman’s correlation coefficients can be used, we suggest Kendall’s correlation coefficient. It aligns well with the nonparametric spirit of this article and has shown better performance in preliminary experiments not reported here and in previous studies on related topics (Dvořák et al. 2022).

3.2.3 Choice of Bandwidth

For the construction of the smoothed residual field $\tilde{s}(u)$ in (6) one has to select a specific kernel function k (a probability density function). The type of the kernel does not play an important role, and we use the Gaussian kernel. On the other hand, the choice of bandwidth (standard deviation of the kernel function) affects the properties of the estimates to a great extent. Traditional rules of thumb or more involved methods may be used for bandwidth selection in this case, see Baddeley, Rubak, and Turner (2015, sec. 6.5.1.2) or Cronie and Van Lieshout (2018). However, whenever available, expert knowledge about the specific problem at hand should guide the choice of bandwidth.

3.3 Partial Correlation Coefficient between a Point Process and a Covariate

When several possibly correlated covariates are available, one might be interested in assessing the strength of dependence between the point process X and the covariate of interest $C_{m+1}$ after removing the possible influence of the remaining (nuisance) covariates $C_1,\dots ,C_m$, in the spirit of the partial correlation coefficient.

The strength of dependence can be quantified by some measure of dependence between the covariate of interest $C_{m+1}$ and the smoothed residual field $\tilde{s}$ from (5) where the possible influence of the nuisance covariates $C_1,\dots ,C_m$ on X has been removed. When a parametric model for the intensity function of X is available, parametric residuals (3) may be used instead.

We suggest using Kendall’s correlation coefficient to quantify the dependence. Again, we consider a set of sampling points $\{y_1,\dots ,y_n\}$, independently and uniformly distributed in W, independent of X and $C_1,\dots ,C_{m+1}$, and define the sample version of the partial correlation coefficient as (9) $$\begin{array}{c}{\widehat{\tau }}_p=\frac{1}{n(n-1)}\sum _{i\ne j}\text{sgn}(C_{m+1}(y_i)\\ -C_{m+1}(y_j))\text{ sgn}(\tilde{s}(y_i)-\tilde{s}(y_j)).\end{array}$$(9)

The population version can be defined in a similar way as in (7). Concerning the choice of the sampling points and the choice of the measure of dependence, comments from the previous section apply here, too.

To illustrate the use of the partial correlation coefficient in quantifying the strength of dependence between a point process and a covariate of interest, after removing the influence of nuisance covariates, we performed the following experiment. The point process model is the Poisson process from Section 3.2. Its intensity function depends in a log-linear way on the covariate $C_1((x,y))=x$, now treated as a nuisance covariate. Specifically, the intensity function is proportional to $\hspace{0.17em}\mathrm{\text{exp}}\hspace{0.17em}\left\{\mathit{\text{ax}}\right\}$. The covariate of interest is $C_2((x,y))=x+y$. We consider 500 independent realizations of the point process for each value of a and compute the means of $\widehat{\tau }$ and ${\widehat{\tau }}_p$. The correlation coefficient $\widehat{\tau }$, again computed with bw = 0.5, correctly indicates that the point process depends on the covariate C₂ through the x-coordinate (black curve in the right panel of Figure 3). On the other hand, the partial correlation coefficient ${\widehat{\tau }}_p$, computed with the adaptive choice of bandwidth described below, attains approx. zero values in this case (red curve in the right panel of Figure 3), implying that the influence of the nuisance covariate C₁ was successfully removed.

3.3.1 Choice of Bandwidth

Construction of the smoothed residual field requires choosing a bandwidth for the smoothing kernel. Again, standard recommendations may be employed, or the available expert knowledge may be used. However, in our pilot experiments with a single nuisance covariate C₁, we observed that the influence of C₁ was usually not completely removed from X during the construction of the smoothed residual field $\tilde{s}(u)$, in the sense that the empirical Kendall’s correlation coefficient of $\{(\tilde{s}(y_j),C_1(y_j)),j=1,\dots ,n\}$ was nonzero. Its value was strongly influenced by the value of bandwidth.

To remove this effect, we suggest selecting the bandwidth value (from a given finite set of candidate values) that minimizes the absolute value of the empirical Kendall’s coefficient of $\{(\tilde{s}(y_j),C_1(y_j)),j=1,\dots ,n\}$, denoted $\widehat{\tau }(\tilde{s},C_1)$ in the following. In this way, we select the bandwidth value that removes the influence of C₁ on X the most successfully and it can be seen as a conservative version of the correlation coefficient. This is important mostly in cases where the nuisance covariate is correlated with the covariate of interest. For independent covariates, this procedure has very little effect on the performance of the random shift tests. We apply this approach to bandwidth selection in our simulation experiments below. When more than one nuisance covariate is available, this adaptive bandwidth procedure can be generalized by minimizing ${\sum }_{i=1}^m\widehat{\tau }{(\tilde{s},C_i)}^2$.

3.4 Covariate-Weighted Residual Measure

While ${\widehat{\tau }}_p$ is useful for quantifying the strength of dependence between X and the covariate of interest $C_{m+1}$ after removing the influence of nuisance covariates $C_1,\dots ,C_m$, the random shift test using ${\widehat{\tau }}_p$ as the test statistic turned out to have a rather low power in our simulation studies. The reason lies in the applied smoothing and the deliberate removal of the preferential sampling effects—the association between the points of X and the covariate $C_{m+1}$ brings important information.

To overcome these issues, we define the following characteristic that we call the covariate-weighted residual measure of W: (10) $$\begin{array}{c}\mathit{\text{CWR}}={\int }_W{C}_{m+1}(u)\tilde{\mathcal{R}}(\mathrm{d}u)=\sum _{x_i\in X\cap W}{C}_{m+1}(x_i)\\ -{\int }_W{C}_{m+1}(u)\widehat{\rho }(C_1(u),\dots ,C_m(u))\mathrm{d}u.\end{array}$$(10)

This can be viewed as a generalization of the test statistic T from Section 2.6 which also includes the sum of covariate values, but does not take into account possible nuisance covariates. By sampling the values of $C_{m+1}$ at the points of X we take advantage of any possible preferential sampling effects, and no smoothing is performed when computing the value of CWR, hence, we avoid the problem of bandwidth selection.

The expectation of CWR is close to 0 if the covariates $C_1,\dots ,C_m$ capture all variation in $\lambda (u)$. The heuristic reasoning behind this is the following. Assume that the covariate $C_{m+1}$ is deterministic, then the expectation of the sum in (10) is ${\int }_W{C}_{m+1}(u)\lambda (u)\mathrm{d}u$ according to Campbell’s theorem, where $\lambda (u)$ is the true intensity function of X. The integral ${\int }_W{C}_{m+1}(u)\widehat{\rho }(C_1(u),\dots ,C_m(u))\mathrm{d}u$ is close to ${\int }_W{C}_{m+1}(u)\lambda (u)\mathrm{d}u$ if $\widehat{\rho }(C_1(u),\dots ,C_m(u))$ is close to $\lambda (u)$. Thus, if $\widehat{\rho }(C_1(u),\dots ,C_m(u)$ is close to $\lambda (u)$, CWR is expected to be close to 0, and to differ from 0 otherwise. This enables testing the significance of $C_{m+1}$ after removing the influence of $C_1,\dots ,C_m$.

3.5 Testing the Covariate Significance Under the Presence of Nuisance Covariates

Now we focus on the null hypothesis that X and $C_{m+1}$ are independent, conditionally on $C_1,\dots ,C_m$. We employ the random shift test described in Section 2.5, either with torus or variance correction. The test statistic can be ${\widehat{\tau }}_p$ in which case the two spatial objects to be shifted against each other are the two random fields $\Phi =\tilde{s}$ and $\Psi =C_{m+1}$. Alternatively, one can use the covariate-weighted residual measure of W as the test statistic. In this case $\Phi =\tilde{R}$ is a measure and $\Psi =C_{m+1}$ is a random field. If v_i is a shift vector, the shift of the random field $\Psi$ should be interpreted in both cases as $(\Psi +v_i)(u)=\Psi (u-v_i)$. The choice of the correction factors for the variance correction is discussed in Appendix C in the supplementary material, including Proposition 1 which studies the variance of CWR for Poisson processes and an empirical study for log-Gaussian Cox processes. The assumption of stationarity of one of the spatial objects is discussed in detail in Section 6.

4 Simulation Study

To assess the performance of the proposed tests, we present below a set of simulation experiments, both under the null hypothesis and under various alternatives. The models range from clustering through complete spatial randomness to regularity, even combining clustering and inhibition on different scales. The null hypothesis states that X and $C_{m+1}$ are independent, conditionally on the nuisance covariates. For simplicity, we focus on the situation with a single nuisance covariate. The proposed nonparametric tests are compared with the parametric methods available in standard software represented by the spatstat package.

4.1 Simulation Study Design

The following notation and choices are used in all simulation experiments. $Z_1,Z_2,\dots$ are independent, identically distributed Gaussian random fields, centered, unit variance, with exponential covariance function with scale 0.1. The observation window is $W={[0,1]}^2$. The expected number of points in W is equal to $\hspace{0.17em}\mathrm{\text{exp}}\hspace{0.17em}\left\{5\right\}\doteq 148.4$ for Poisson and clustered models and approximately equal to $\hspace{0.17em}\mathrm{\text{exp}}\hspace{0.17em}\left\{5\right\}$ for models exhibiting regularity.

For each model, we simulate 5000 independent realizations, and for each realization, we perform a set of tests on the 5% nominal significance level. In the tables of results we report the fractions of rejections for the individual tests, rounded to three decimals. To assess the liberality or conservativeness of the tests, one can compare the reported rejection rates (in experiments performed under the null hypothesis) with the interval based on the 2.5% and 97.5% quantiles of the binomial distribution with parameters n = 5000 and p = 0.05, that is, with the interval $[0.0440,0.0562]$.

We investigate the performance of the random shift tests with either ${\widehat{\tau }}_p$ or CWR as the test statistic, with either parametric or nonparametric version of residuals (denoted by the symbol “p” or “n” in the tables of results) and with either torus or variance correction (denoted by “tor” or “var” in the tables of results). The tests are based on 999 independent random shift vectors with uniform distribution on a disk with radius 0.5. The values of ${\widehat{\tau }}_p$ are computed with the bandwidth selected by the adaptive procedure from Section 3.3.1 and with 100 sampling points chosen uniformly and independently in W. The random shift tests are compared with the parametric tests provided by the functions ppm (for Poisson, Strauss, and hardcore Strauss processes) and kppm (for log-Gaussian Cox processes, denoted LGCP in the following) from the spatstat package, see Section 2.1. The test assuming Poisson interactions is denoted “ppm, Poisson” in the tables of results below. Similarly, “ppm, Strauss(R)” denotes the test assuming the Strauss interactions with interaction distance R, “ppm, HC Strauss(R)” denotes the test assuming the hardcore Strauss model with interaction distance R and the true hardcore distance, and “kppm, LGCP” denotes the test assuming the LGCP model.

To mimic the practical issues with model specification, we consider these parametric tests both with the correct interaction model and with an incorrect interaction model of a similar type. Specifically, in addition to fitting the correct model to the LGCPs we also fit a Matérn cluster process (the test is denoted “kppm, MC” in the tables of results); for the Strauss and hardcore Strauss processes we fit the models with the interaction distance fixed to either the correct or incorrect value, specified in the tables of results in the column “Variant”. For the hardcore Strauss processes the correct hardcore distance is always assumed in the model fitting. We also fit an inhomogeneous Poisson process to all datasets to investigate the effect of ignoring the interaction structure. On the other hand, we do not try fitting clustered models to clearly regular datasets and vice versa.

All the parametric tests assume the log-linear model for the intensity function (1), even though for some point process models we consider below this does not hold. This also illustrates possible issues with model misspecification in practice.

4.2 Significance Level Under Independent Covariates

In the following, we let the nuisance covariate influencing the intensity function of the point process be $C_1(u)=Z_1(u)$ and the covariate of interest be $C_2(u)=Z_3(u)$, which means that the covariate of interest C₂ is independent of the nuisance covariate C₁ and the point process X (we study the significance level under dependent covariates in the Appendix A in the supplementary material). For the construction of the LGCP models we also use the random field Z₂, which is responsible for interactions in the point process rather than variation in its intensity function. For the Poisson and LGCP models, the covariate C₁ influences the intensity function directly. For the Strauss and hardcore Strauss models, it directly influences the trend function $\beta (u)$. This influence is transformed to the intensity function in a nontrivial way. We consider the following models:

• (P₁)  Poisson process with intensity function $\lambda (u)=\hspace{0.17em}\mathrm{\text{exp}}\hspace{0.17em}\{4.5+Z_1(u)\}$

• (P₂)  Poisson process with intensity function $\lambda (u)=\hspace{0.17em}\mathrm{\text{exp}}\hspace{0.17em}\left\{5\right\}·Z_1{(u)}^2$

• (L₁)  LGCP with driving intensity function $\Lambda (u)=\hspace{0.17em}\mathrm{\text{exp}}\hspace{0.17em}\{4.0+Z_1(u)+Z_2(u)\}$

• (L₂)  LGCP with driving intensity function $\Lambda (u)=\hspace{0.17em}\mathrm{\text{exp}}\hspace{0.17em}\{4.5+Z_2(u)\}·Z_1{(u)}^2$

• (S₁)  Strauss process with interaction parameter $\gamma =0.5$, interaction range R = 0.05 and trend $\beta (u)=220·\hspace{0.17em}\mathrm{\text{exp}}\hspace{0.17em}\{Z_1(u)\}$

• (S₂)  Strauss process with $\gamma =0.5,R=0.05$ and $\beta (u)=350·Z_1{(u)}^2$

• (H₁)  Strauss process with hardcore distance hc = 0.01, interaction parameter γ = 4, interaction distance R = 0.02 and trend $\beta (u)=180·\tilde{Z}(u)$, where $\tilde{Z}(u)=c·\hspace{0.17em}\mathrm{\text{exp}}\hspace{0.17em}\{Z_1(u)/5\}$ and c is chosen for each realization so that the maximum of the given realization of $\tilde{Z}(u)$ over W is 1.

• (H₂)  Strauss process with hardcore distance hc = 0.01, γ = 4, R = 0.02 and $\beta (u)=120·\tilde{Z}(u)$, where $\tilde{Z}(u)=c·\mathrm{\text{max}}(1-Z_1{(u)}^2/5,0)$, again scaled for each realization to attain the maximum value of 1 over W.

Note that in the notation for the models the letter represents the type of interaction in the point process while the subscript specifies whether the covariate C₁ influences the intensity function in a log-linear way (denoted by 1) or in a quadratic way (denoted by 2).

Since the covariate of interest C₂ is independent of X, the tests should reject in 5% of cases. Table 1 shows the fractions of rejection. We make the following observations:

Table 1 Estimated sizes of the tests for experiments with independent covariates.

Test	Variant	P1	P2	L1	L2	S1	S2	H1	H2
${\widehat{\tau }}_p$	p, tor	0.053	0.058	0.059	0.065	0.056	0.062	0.063	0.061
${\widehat{\tau }}_p$	n, tor	0.049	0.057	0.059	0.063	0.058	0.060	0.061	0.062
${\widehat{\tau }}_p$	p, var	0.040	0.043	0.031	0.047	0.043	0.049	0.047	0.045
${\widehat{\tau }}_p$	n, var	0.043	0.046	0.035	0.047	0.046	0.047	0.047	0.050
CWR	p, tor	0.049	0.053	0.055	0.051	0.048	0.051	0.057	0.048
CWR	n, tor	0.049	0.050	0.054	0.050	0.045	0.045	0.056	0.051
CWR	p, var	0.050	0.049	0.050	0.054	0.047	0.050	0.056	0.048
CWR	n, var	0.048	0.051	0.049	0.048	0.048	0.045	0.055	0.049
ppm	Poisson	0.048	0.166	0.268	0.331	0.021	0.052	0.198	0.161
kppm	LGCP	0.020	0.023	0.080	0.027	–	–	–	–
kppm	MC	0.041	0.025	0.086	0.030	–	–	–	–
ppm	Str(0.02)	0.039	0.111	–	–	0.019	0.040	–	–
ppm	Str(0.05)	0.044	0.070	–	–	0.038	0.060	–	–
ppm	Str(0.10)	0.041	0.063	–	–	0.028	0.034	–	–
ppm	HC Str(0.02)	–	–	–	–	–	–	0.134	0.137
ppm	HC Str(0.05)	–	–	–	–	–	–	0.095	0.079
ppm	HC Str(0.10)	–	–	–	–	–	–	0.093	0.092

NOTE: Different rows correspond to different tests, either the nonparametric tests proposed here (${\widehat{\tau }}_p$, CWR) or the parametric tests based on different assumptions, denoted by the name of the spatstat function used to perform them (ppm, kppm). Different columns correspond to different models.

• The nonparametric tests match the nominal significance level correctly for all models, the tests based on CWR match it slightly more precisely than those based on ${\widehat{\tau }}_p$. Both the torus correction and the variance correction perform well, with only a slight tendency toward liberality observed for the torus correction and the tests based on ${\widehat{\tau }}_p$.

• Parametric tests assuming correct interaction structure and correct model for the intensity function (denoted by 1) match the nominal significance level correctly for the Poisson model P₁ (the estimated size of the test 0.048) while being liberal for the LGCP model L₁ (0.080) or the hardcore Strauss model H₁ (0.095). The parametric test is conservative for the Strauss model S₁ (0.038).

• Parametric tests assuming correct interaction structure and incorrect model for the intensity function (denoted by 2) may exhibit very strong liberality (P₂ – 0.166, H₂ – 0.079) or conservativeness (L₂ – 0.027).

• Parametric tests assuming incorrect interaction structure, even with the correct model for the intensity function, may exhibit very strong liberality, for example assuming Poisson interactions for an LGCP (L₁ – 0.268) or for a hardcore Strauss process (H₁ – 0.198). Also, they may exhibit conservatives, for example assuming attractive interactions for a Poisson process (P₁ – 0.020, 0.041).

These observations illustrate that parametric tests are prone to perform poorly under model misspecification either in terms of the interaction structure or the intensity function. However, even when both of these model components are specified correctly, there is a risk of strong liberality of the parametric tests with the sample sizes considered here. From this point of view, the nonparametric tests are preferable as they match the nominal significance level correctly for all models in this study.

4.3 Power under Dependent Covariates

In this section, we study the power of the tests in situations where the covariate of interest C₂ influences the intensity function of X even after removing the effect of the nuisance covariate C₁. We consider models similar to those in the previous sections and let $C_1(u)=Z_1(u)$ and $C_2(u)=Z_1(u)+2Z_3(u)$. We focus on the case with dependent covariates, which is more challenging for our proposed tests as they show conservativeness in this case, see Appendix A in the supplementary material. The models depend on a parameter a > 0 that controls the strength of dependence between X and the covariate of interest C₂. The value of a is chosen so that all the tests exhibit nontrivial powers, that is not close to 0.05 and not close to 1.00. The models are given as follows:

• ($P_1^p$)  Poisson process with intensity function $\lambda (u)=\hspace{0.17em}\mathrm{\text{exp}}\hspace{0.17em}\{4.5+Z_1(u)+a{Z}_3(u)-a^2/2\}$ with a = 1/4.

• ($P_2^p$)  Poisson process with intensity function $\lambda (u)=\hspace{0.17em}\mathrm{\text{exp}}\hspace{0.17em}\{5.0+a{Z}_3(u)-a^2/2\}·Z_1{(u)}^2$ with a = 1/4.

• ($L_1^p$)  LGCP with driving intensity function $\Lambda (u)=\hspace{0.17em}\mathrm{\text{exp}}\hspace{0.17em}\{4.0+Z_1(u)+Z_2(u)+a{Z}_3(u)-a^2/2\}$ with a = 1/2.

• ($L_2^p$)  LGCP with driving intensity function $\Lambda (u)=\hspace{0.17em}\mathrm{\text{exp}}\hspace{0.17em}\{4.5+Z_2(u)+a{Z}_3(u)-a^2/2\}·Z_1{(u)}^2$ with a = 1/2.

• ($S_1^p$)  Strauss process with interaction parameter $\gamma =0.5$, interaction range R = 0.05 and trend $\beta (u)=210·\hspace{0.17em}\mathrm{\text{exp}}\hspace{0.17em}\{Z_1(u)+a{Z}_3(u)\}$ with a = 1/4.

• ($S_2^p$)  Strauss process with $\gamma =0.5,R=0.05$ and $\beta (u)=350·\hspace{0.17em}\mathrm{\text{exp}}\hspace{0.17em}\{a{Z}_3(u)\}·Z_1{(u)}^2$ with a = 1/4.

• ($H_1^p$)  Strauss process with hardcore distance hc = 0.01, interaction parameter γ = 4, interaction distance R = 0.02 and trend $\beta (u)=190·\tilde{Z}(u)$, where $\tilde{Z}(u)=c·\hspace{0.17em}\mathrm{\text{exp}}\hspace{0.17em}\{Z_1(u)/5+a{Z}_3(u)-a^2/2\}$ and c is chosen for each realization so that the maximum of the given realization of $\tilde{Z}(u)$ over W is 1.

• ($H_2^p$)  Strauss process with hardcore distance hc = 0.01, γ = 4, R = 0.02 and $\beta (u)=170·\tilde{Z}(u)$, where $\tilde{Z}(u)=c·\hspace{0.17em}\mathrm{\text{exp}}\hspace{0.17em}\{a{Z}_3(u)-a^2/2\}·\mathrm{\text{max}}(1-Z_1{(u)}^2/5,0)$, again scaled for each realization to attain the maximum value of 1 over W.

Table 2 shows the fractions of rejections for the individual tests for the eight models specified above. We make the following observations:

Table 2 Estimated powers of the tests for experiments with correlated covariates.

Test	Variant	$P_1^p$	$P_2^p$	$L_1^p$	$L_2^p$	$S_1^p$	$S_2^p$	$H_1^p$	$H_2^p$
${\widehat{\tau }}_p$	p, tor	0.198	0.126	0.153	0.155	0.127	0.114	0.101	0.128
${\widehat{\tau }}_p$	n, tor	0.215	0.214	0.164	0.199	0.137	0.152	0.104	0.134
${\widehat{\tau }}_p$	p, var	0.160	0.105	0.113	0.121	0.100	0.096	0.076	0.100
${\widehat{\tau }}_p$	n, var	0.177	0.178	0.125	0.164	0.111	0.122	0.080	0.109
CWR	p, tor	0.758	0.444	0.820	0.702	0.608	0.396	0.310	0.405
CWR	n, tor	0.793	0.677	0.846	0.791	0.690	0.514	0.331	0.418
CWR	p, var	0.753	0.441	0.808	0.698	0.604	0.392	0.304	0.402
CWR	n, var	0.789	0.675	0.835	0.783	0.688	0.511	0.324	0.416
ppm	Poisson	0.774	0.720	0.960	0.943	0.507	0.483	0.618	0.714
kppm	LGCP	0.644	0.254	0.841	0.477	–	–	–	–
kppm	MC	0.745	0.264	0.848	0.499	–	–	–	–
ppm	Str(0.02)	0.682	0.579	–	–	0.432	0.381	–	–
ppm	Str(0.05)	0.614	0.456	–	–	0.467	0.399	–	–
ppm	Str(0.10)	0.515	0.387	–	–	0.344	0.263	–	–
ppm	HC Str(0.02)	–	–	–	–	–	–	0.346	0.375
ppm	HC Str(0.05)	–	–	–	–	–	–	0.302	0.340
ppm	HC Str(0.10)	–	–	–	–	–	–	0.307	0.355

NOTE: Different rows correspond to different tests, either the nonparametric tests proposed here (${\widehat{\tau }}_p$, CWR) or the parametric tests based on different assumptions, denoted by the name of the spatstat function used to perform them (ppm, kppm). Different columns correspond to different models.

• For both ${\widehat{\tau }}_p$ and CWR, the versions of the test based on nonparametric residuals exhibit higher power than those based on parametric residuals.

• The tests based on ${\widehat{\tau }}_p$ have very low power due to the smoothing and removal of the preferential sampling effects.

• The tests based on CWR exhibit very high power comparable to the parametric tests with correct interaction model and correct model for the intensity function (for P, L, and H models) or even higher power (S).

• When the parametric tests are used with the correct interaction model and incorrect model for the intensity function, the nonparametric tests based on CWR have much higher power (L, S), slightly higher power (H), or direct comparison is not possible due to severe liberality of the parametric test (P).

• The torus correction and the variance correction perform nearly equivalently for tests based on CWR, while for tests based on ${\widehat{\tau }}_p$ the torus correction shows slightly higher power, which can be explained by the small liberality of these tests observed in Section 4.2.

These observations indicate that the random shift tests based on CWR with nonparametric residuals and either torus or variance correction can be preferred in practice to parametric tests since the possible issues with model misspecification are avoided without compromising the power of the test.

5 Nonparametric Covariate Selection for the BCI Dataset

To illustrate the possibility to use the proposed random shift tests for covariate selection, we consider now the BCI dataset described in Section 1.2. Five covariates are available that possibly influence the intensity function of the point process. A possible way to select the set of covariates that have a significant effect on the intensity function is the backward selection procedure described in the following. The numerical results are given in Table 3.

Table 3 Backward selection of covariates for the BCI dataset.

	Stage 1	Stage 2	Stage 3	Stage 4
Elevation	0.168	0.186	0.018	0.016
Gradient	0.228	0.240	0.232	–
Nitrogen	0.426	0.388	–	–
Phosphorus	0.066	0.128	0.220	0.048
Potassium	0.598	–	–	–

Individual cells show the p-values of the random shift tests. The row indicates the covariate of interest, while all other covariates considered in the given stage (column) are considered to be the nuisance covariates.

We start in stage 1 with all five covariates, and for each of those we perform the random shift test where the given covariate is the covariate of interest and the remaining four covariates are considered to be the nuisance covariates. We use the test based on CWR with nonparametric residuals and torus correction, with 999 random shifts where the shift vectors have uniform distribution on a disk with radius 250 meters. The covariate with the highest p-value (potassium in this case, printed in italics in Table 3) is removed, and the procedure is repeated in stage 2 with the four covariates. In this stage, the nitrogen covariate is removed, then the gradient covariate, and finally in stage 4 where only two covariates are considered (elevation, phosphorus), both covariates are found significant on the 5% significance level, see Table 3. We conclude that these two covariates significantly affect the intensity function of the point process and should be included in the further steps of the inference. Other covariates can be disregarded without losing important information. Figure 4 shows the locations of the trees together with the nonparametric estimate of the intensity function depending on all five available covariates (see Section 2.3 for details about the estimator) and the nonparametric estimate depending on the two covariates selected above (elevation, phosphorus). We observe that the two-covariate model describes the data quite similarly to the full, five-covariate model, and we are not losing important information by disregarding the three covariates identified as redundant by the selection procedure above.

Fig. 4 The Barro Colorado Island dataset. Left to right: locations of trees, nonparametric estimate of the intensity function in the full model, nonparametric estimate in the final model with two covariates.

Concerning computational times, performing the five tests in stage 1 took approx. 66 sec on a regular laptop. For the consequent stages, the total computational time for performing all the required tests in the given stage was 13.7, 9.50, and 6.6 sec, respectively. Up to three nuisance covariates, computing the test statistic values for the randomly shifted data takes the most time. With higher number of nuisance covariates, computing the nonparametric estimate of the intensity function depending on the nuisance covariates (see Section 2.3) becomes the bottleneck. For example, computing the intensity function estimate in the full model with five covariates in Figure 4 took approx. 73 sec.

For comparison, we have also fitted the log-linear model (1) with the five covariates considered here, using the kppm function from the spatstat package as described in Section 2.1. We assume the Thomas type of interactions as suggested in Baddeley, Rubak, and Turner (2015, sec. 12.4.4). With this approach, three covariates are found significant on the 5% significance level: elevation, gradient and phosphorus, with p-values 0.019, 0.044 and ${10}^{-4}$, respectively. Two of these covariates were also found significant by the nonparametric procedure described above (elevation and phosphorus, see Table 3). On the other hand, the gradient covariate was found borderline significant by the parametric approach and not significant by the nonparametric procedure.

6 Conclusions and Discussion

Parametric methods for testing the covariate significance may be severely liberal in cases where the interaction structure is misspecified. This is confirmed by the simulation experiments presented here. The same holds in cases where the assumed model for the intensity function is incorrect. Furthermore, the simulation study indicated that even when both the interaction model and the intensity function model are correctly specified, the parametric methods may be liberal due to the limited data size and asymptotic nature of the parametric methods.

The methods proposed in this article allow quantification and testing of the significance of the correlations between a point process and a covariate of interest, possibly after removing the influence of nuisance covariates. We stress that the proposed methods can be used without specifying any model for the data. The simulation experiments reported in Section 4 show that the random shift tests based on ${\widehat{\tau }}_p$ or CWR match the nominal significance level correctly in all situations.

Concerning power, the nonparametric tests based on CWR exhibit comparable or even higher power than parametric tests under the correct model, while showing higher power than parametric tests under incorrect models (where either the interaction or the intensity function is misspecified). This indicates the superiority of the CWR tests over parametric tests in practical applications where the true model is not known. Hence, using the proposed nonparametric CWR tests for covariate selection, for example as discussed in Section 5, provides more reliable results than the available parametric tests are able to provide, and the selected covariates can be used in the further steps of inference with greater confidence.

The only assumption of the proposed random shift tests is that at least one of the objects is stationary under the null hypothesis so that its distribution is not affected by the shifts. Either the covariate of interest can be assumed to be stationary or the covariate-weighted residual measure or the smoothed residual field can be assumed to be close to stationarity if all the relevant covariates are used in the construction of the residuals.

In practice, the true process generating the data is not known. Parametric methods rely on model assumptions which may or may not be appropriate. Even when an appropriate model is determined, the corresponding parametric test may not be available or even may not be developed yet. The nonparametric tools proposed in this article do not suffer from this issue and can be used for any type of process.

A natural question is whether the proposed methods are applicable to categorical covariates. If one of the nuisance covariates is categorical, nonparametric estimation of the intensity function may be performed separately on the individual subregions of W determined by the categorical covariate, allowing all the proposed methods to be used as described above. If the categorical covariate is the covariate of interest, computing $\widehat{\tau }$ or ${\widehat{\tau }}_p$ is not relevant due to the ties in the data. However, the observation window W can be separated into subregions $W_1,\dots ,W_k$ determined by the values of the covariate of interest. The values $V_i=n(X\cap W_i)-{\int }_{W_i}\widehat{\rho }(C_1(u),\dots ,C_m(u))\mathrm{d}u,i=1,\dots ,k,$ can be used to form a vector test statistic $(V_1,\dots ,V_k)$ and the random shift test can be performed for example by means of the global envelope test (Myllymäki et al. 2017). This approach corresponds to the determination of differences between point process intensities in subregions $W_1,\dots ,W_k$.

Supplementary Materials

Supplementary text: This supplementary material contains several appendices to the above manuscript. (pdf file)

R-code for the example:R-code to compute the BCI data example. The data are available in the R-package spatstat. (r file)

Disclosure Statement

The authors declare there are no competing interests, financial or otherwise.

Additional information

Funding

No projects supported this work.