Information criteria for inhomogeneous spatial point processes

The density of a Poisson point process with intensity ${\rho_{n}}(\cdot;\boldsymbol{\beta}_{l})$ and observed in $W_{n}$ is given by (see e.g. Møller and Waagepetersen (2004))

$$p_{n}(\mathbf{x};\boldsymbol{\beta}_{l})=\left\{\prod_{u\in{\mathbf{x}}}{\rho_{n}}(u;\boldsymbol{\beta}_{l})\right\}\exp\left\{|W_{n}|-\int_{W_{n}}{\rho_{n}}(u;\boldsymbol{\beta}_{l})\mathrm{d}u\right\}$$

(2)

for locally finite point configurations $\mathbf{x}\subset W_{n}$. We emphasize that in the following we assume neither that $\mathbf{X}_{n}$ introduced in the previous section is a Poisson process nor that ${\rho_{n}}(\cdot;\boldsymbol{\beta}_{l})$ coincides with the intensity function $\lambda_{n}$ of $\mathbf{X}_{n}$.

Combining (1) and (2), up to a constant, the $\log$ of (2) evaluated at $\mathbf{X}_{n}$ becomes

$$\ell_{n}(\boldsymbol{\beta}_{l})=\sum_{u\in{\mathbf{X}_{n}}}\boldsymbol{\beta}_{l}^{\top}\mathbf{z}_{l}(u)-\theta_{n}\int_{W_{n}}\exp\{\boldsymbol{\beta}_{l}^{\top}\mathbf{z}_{l}(u)\}\mathrm{d}u.$$

(3)

Let $e_{n}(\boldsymbol{\beta}_{l})$ be the corresponding estimating function given by

$$e_{n}(\boldsymbol{\beta}_{l})=\frac{\mathrm{d}}{\mathrm{d}\boldsymbol{\beta}_{l}}\ell_{n}(\boldsymbol{\beta}_{l})=\sum_{u\in{\mathbf{X}_{n}}}\mathbf{z}_{l}(u)-\theta_{n}\int_{W_{n}}\mathbf{z}_{l}(u){\rho_{n}}(u;\boldsymbol{\beta}_{l})\mathrm{d}u.$$

(4)

For any $\boldsymbol{\beta}_{l}\in{\mathbb{R}^{p_{l}}}$, the sensitivity (or Fisher information) matrix is

$$\mathbf{S}_{n}(\boldsymbol{\beta}_{l})=-{\mathbb{E}}\left\{\frac{\mathrm{d}}{\mathrm{d}\boldsymbol{\beta}_{l}^{\top}}e_{n}(\boldsymbol{\beta}_{l})\right\}=\theta_{n}\int_{W_{n}}\mathbf{z}_{l}(u)\mathbf{z}_{l}(u)^{\top}{\rho_{n}}(u;\boldsymbol{\beta}_{l})\mathrm{d}u.$$

(5)

We assume $\mathbf{S}_{n}(\boldsymbol{\beta}_{l})$ is positive definite for all $\boldsymbol{\beta}_{l}$ (see also condition C5 in the next section). We can then define the estimator of $\boldsymbol{\beta}_{l}$ as

$$\hat{\boldsymbol{\beta}}_{l,n}=\operatorname*{arg\,max}_{\boldsymbol{\beta}_{l}\in{\mathbb{R}^{p_{l}}}}\,{p_{n}}(\mathbf{X}_{n};\boldsymbol{\beta}_{l})=\operatorname*{arg\,max}_{\boldsymbol{\beta}_{l}\in{\mathbb{R}^{p_{l}}}}\,\ell_{n}(\boldsymbol{\beta}_{l}).$$

(6)

If $\mathbf{X}_{n}$ is indeed a Poisson process with intensity function ${\rho_{n}}(\cdot;\boldsymbol{\beta}_{l})$, then $\hat{\boldsymbol{\beta}}_{l}$ is the maximum likelihood estimator and the sensitivity (5) equals the observed information matrix $-\mathrm{d}e_{n}(\boldsymbol{\beta}_{l})/\mathrm{d}\boldsymbol{\beta}_{l}^{\top}$.

If $\mathbf{X}_{n}$ is not Poisson, $\hat{\boldsymbol{\beta}}_{l,n}$ may be viewed as a composite likelihood estimator (Schoenberg, 2005; Waagepetersen, 2007). In the situation where the intensity function ${\rho_{n}}(\cdot;\boldsymbol{\beta}_{l})$ coincides with the true intensity function $\lambda_{n}$, asymptotic properties of maximum likelihood or composite likelihood estimators obtained as maximizers of Poisson likelihood functions have been established in various settings by Rathbun and Cressie (1994), Waagepetersen (2007), Guan and Loh (2007) and Waagepetersen and Guan (2009). In the next section we investigate the more intriguing situation where the intensity model is misspecified.

To handle the situation where ${\rho_{n}}(\cdot;\boldsymbol{\beta}_{l})$ does not coincide with the intensity function of $\mathbf{X}_{n}$, we follow Varin and Vidoni (2005) and define a (composite) Kullback-Leibler divergence between the model ${\cal M}_{l}$ with parameter $\boldsymbol{\beta}_{l}$ and the true sampling distribution. That is,

$$\mathrm{KL}_{n}(\boldsymbol{\beta}_{l})={\mathbb{E}}\left\{\ell_{n}-\ell_{n}(\boldsymbol{\beta}_{l})\right\}$$

(7)

where $\ell_{n}$ is the Poisson log-likelihood obtained with the true intensity $\lambda_{n}$. For a window $W_{n}$ and model ${\cal M}_{l}$ we let

$$\boldsymbol{\beta}_{l,n}^{*}=\operatorname*{arg\,max}_{{\boldsymbol{\beta}_{l}\in\mathbb{R}^{p_{l}}}}{\mathbb{E}}\ell_{n}(\boldsymbol{\beta}_{l})$$

denote the ‘least wrong parameter value’ under model ${\cal M}_{l}$, provided the maximum exists. It is easy to see by explicit evaluation of the right hand side that

$$-\frac{\mathrm{d}}{\mathrm{d}{\boldsymbol{\beta}_{l}}^{\top}}\left\{-\frac{\mathrm{d}}{\mathrm{d}\boldsymbol{\beta}_{l}}{\mathbb{E}}\ell_{n}(\boldsymbol{\beta}_{l})\right\}=\mathbf{S}_{n}(\boldsymbol{\beta}_{l})$$

and so condition C5 stated below implies that $\boldsymbol{\beta}_{l,n}^{*}$ is well-defined as a unique maximum when $n$ is large enough. Also it is easy to see that

$${\mathbb{E}}e_{n}(\boldsymbol{\beta}_{l,n}^{*})=0,$$

(8)

which means that $\hat{\boldsymbol{\beta}}_{l,n}$ given by (6) is a candidate to estimate $\boldsymbol{\beta}_{l,n}^{*}$.

The remainder of this section is devoted to asymptotic results for $\hat{\boldsymbol{\beta}}_{l,n}$ within the above framework of a misspecified intensity function. We thereby extend the results in the references mentioned in Section 3.1. In contrast to these references which used either increasing domain or infill asymptotics, we moreover consider a ‘double asymptotic’ framework as formalized by condition C3 presented below.

Two matrices are crucial for the asymptotic results. The sensitivity matrix $\mathbf{S}_{n}(\boldsymbol{\beta}_{l})$ is given by (5) regardless of whether the model is misspecified or not. The variance-covariance matrix $\boldsymbol{\Sigma}_{l,n}$ of (4) is, using the Campbell theorem, given by

	$\displaystyle\boldsymbol{\Sigma}_{l,n}$	$\displaystyle=\int_{W_{n}}\mathbf{z}_{l}(u)\mathbf{z}_{l}(u)^{\top}\lambda_{n}(u)\mathrm{d}u$
		$\displaystyle\qquad+\int_{W_{n}}\int_{W_{n}}\mathbf{z}_{l}(u)\mathbf{z}_{l}(u)^{\top}\lambda_{n}(u)\lambda_{n}(v)\left\{g_{n}(u,v)-1\right\}\mathrm{d}u\mathrm{d}v.$		(9)

Observe that $\boldsymbol{\Sigma}_{l,n}$ does not depend on $\boldsymbol{\beta}_{l}$ (whence its notation). Our results will be based on the following assumptions where for a square matrix $\mathbf{M}$, $\nu_{\min}(\mathbf{M})$ (resp. $\nu_{\max}(\mathbf{M})$) stands for the smallest (resp. largest) eigenvalue. We use $a_{n}\asymp b_{n}$ to denote that $a_{n}={\mathcal{O}}(b_{n})$ and $b_{n}={\mathcal{O}}(a_{n})$.

1. [C1]

   As $n\to\infty$, $\sup_{n\geq 1}\|\boldsymbol{\beta}_{l,n}^{*}\|=\mathcal{O}(1)$.

2. [C2]

   $\mathbf{z}_{l}:\mathbb{R}^{d}\to\mathbb{R}$ is continuous, $0<\inf_{u\in\mathbb{R}^{d}}\|\mathbf{z}_{l}(u)\|<\sup_{u\in\mathbb{R}^{d}}\|\mathbf{z}_{l}(u)\|<\infty$.

3. [C3]

   The sequence $\left(\tau_{n}:=\theta_{n}|W_{n}|\right)_{n\geq 1}$ is an increasing sequence, such that$\liminf_{n\to\infty}\theta_{n}>0$ and $\lim_{n\to\infty}\tau_{n}=\infty$. The sets $W_{n}$ are convex and compact.

4. [C4]

   As $n\to\infty$, $\tau_{n}^{-1}e_{n}(\boldsymbol{\beta}_{l,n}^{*})\to 0$ almost surely.

5. [C5]

   For any $\boldsymbol{\beta}_{l}\in\mathbb{R}^{p_{l}}$ and $n\geq 1$, $\mathbf{S}_{n}(\boldsymbol{\beta}_{l})$ is positive definite. In addition, for any $n\geq 1$, there exists a set $B_{n}\subseteq W_{n}$ such that $|B_{n}|\asymp|W_{n}|$ and a $c>0$ such that $\inf_{n}\inf_{\boldsymbol{\phi}\in\mathbb{R}^{p_{l}},\|\boldsymbol{\phi}\|=1}\inf_{u\in B_{n}}|\boldsymbol{\phi}^{\top}\mathbf{z}_{l}(u)|\geq c$. Finally, we assume that $\liminf_{n\to\infty}\nu_{\min}\left(\tau_{n}^{-1}\boldsymbol{\Sigma}_{l,n}\right)>0$.

6. [C6]

   As $n\to\infty$, $\|\boldsymbol{\Sigma}_{l,n}\|=\mathcal{O}(\tau_{n})$.

7. [C7]

   As $n\to\infty$, $\boldsymbol{\Sigma}_{l,n}^{-1/2}e_{n}(\boldsymbol{\beta}_{l,n}^{*})\to N(0,\mathbf{I}_{p_{l}})$ in distribution.

We can then state the following asymptotic result which is verified in Appendix A.

We stress that Theorem 1 (ii)-(iii) do not require the strong consistency of $\tau_{n}^{-1}e_{n}(\boldsymbol{\beta}_{l,n}^{*})$ to 0, i.e. condition C4. We conclude by some remarks regarding the assumptions C1-C7. Condition C1 ensures that the sequence of ‘least wrong parameter value‘ does not diverge with $n$.

Condition C3 is different from existing conditions as it embraces both of the standard asymptotic frameworks considered in the literature:

• •

  infill asymptotics: $W_{n}=W$ with $W$ a bounded set of $\mathbb{R}^{d}$ and $\theta_{n}\to\infty$ as $n\to\infty$.

• •

  increasing domain asymptotics: $\theta_{n}=\theta>0$ and $(W_{n})_{n\geq 1}$ is a sequence of bounded domains of $\mathbb{R}^{d}$ such that $|W_{n}|\to\infty$.

It is also valid if both asymptotics are considered at the same time. The assumed convexity in C3 enables the use of the mean value theorem in the proofs of our theoretical results.

In Condition C2, assuming an upper bound for $\|\mathbf{z}_{l}(u)\|$ is quite standard. The upper bound further implies that assuming a lower bound is not really restrictive since for each covariate $z_{j}$ we can always find some $k$ so that $|z_{j}(u)+k|\neq 0$ for any $u\in\mathbb{R}^{d}$. Replacing $z_{j}$ by $z_{j}+k$ while changing the intercept from $\beta_{0}$ to $\beta_{0}-\beta_{j}k$ leaves the model unchanged. The continuity assumption is used to prove the strong consistency of $\hat{\boldsymbol{\beta}}_{l,n}$.

Condition C4 is also used to ensure the strong consistency of $\hat{\boldsymbol{\beta}}_{l,n}$. This condition can be seen as a law of large numbers. In Example 1, we present a class of models where such an assumption is valid under the generalized asymptotic condition C3. We can observe that if there exists $p>0$ such that ${\mathbb{E}}\{\|\tau_{n}^{-1/2}e_{n}(\boldsymbol{\beta}_{l,n}^{*}\|^{p}\}=\mathcal{O}(1)$ and $\sum_{n}\tau_{n}^{-p/2}<\infty$, then C4 ensues from an application of Borel-Cantelli’s lemma. At least in the increasing domain framework, assuming such a moment assumption is quite standard to derive a central limit theorem.

Condition C5 is very similar to the assumption required by Rathbun and Cressie (1994) under the Poisson case or by Waagepetersen and Guan (2009) for more general point processes, within the increasing domain asymptotic framework. Note in particular that C5 combined with C1-C2 ensures that $\liminf_{n\to\infty}\nu_{\min}\{\tau_{n}^{-1}\mathbf{S}_{n}(\boldsymbol{\beta}_{l,n}^{*})\}>0$.

Under the Poisson case, $g_{n}=1$ and so C6 is obviously satisfied if $\sup_{u\in W_{n}}\lambda_{n}(u)=\mathcal{O}(\theta_{n})$. For more general point processes, assume further that $g_{n}=g$ does not depend on $n$ and is invariant under translations. Then, thanks to C2, the assumption

$$\int_{\mathbb{R}^{d}}\{g(o,w)-1\}\mathrm{d}w<\infty.$$

will imply C6. This assumption is quite standard and satisfied by a large class of models, see e.g. Waagepetersen and Guan (2009).

Continuing within the increasing domain framework, C7 was established by Rathbun and Cressie (1994) in the Poisson case, by Guan and Loh (2007) and Waagepetersen and Guan (2009) for $\alpha$-mixing point processes, and by Lavancier et al. (2020) for determinantal point processes. The infill asymptotic framework was used to establish C7 in case of Poisson cluster processes in Waagepetersen (2007). However, the ‘double’ asymptotic framework has never been considered. Below, we provide an example of a model which satisfies C4, C6 and C7 under the new general asymptotic setting C3.

The BIC criterion (see e.g. Schwarz (1978); Lebarbier and Mary-Huard (2006)) defines the best model $\mathcal{M}_{\mathrm{BIC}}$ as

$$\mathcal{M}_{\mathrm{BIC}}=\operatorname*{arg\,max}_{\mathcal{M}_{l},l=1,\ldots,2^{p}}\;\mathrm{P}_{n}\left(\mathcal{M}_{l}\mid\mathbf{X}_{n}\right).$$

(15)

From Bayes formula,

$$\mathrm{P}_{n}(\mathcal{M}_{l}\mid\mathbf{X}_{n})=\frac{p_{n}(\mathbf{X}_{n}|{\mathcal{M}}_{l})\mathrm{P}(\mathcal{M}_{l})}{p_{n}(\mathbf{X}_{n})}.$$

Letting $p(\boldsymbol{\beta}_{l}|{\mathcal{M}}_{l})$ denote the prior density of $\boldsymbol{\beta}_{l}$ given ${\mathcal{M}}_{l}$,

$$p_{n}(\mathbf{X}_{n}|{\mathcal{M}}_{l})=\int_{\mathbb{R}^{p_{l}}}p_{n}(\mathbf{X}_{n};\beta_{l})p(\boldsymbol{\beta}_{l}|{\mathcal{M}}_{l})\mathrm{d}\boldsymbol{\beta}_{l}$$

since $p_{n}(\mathbf{X}_{n};\beta_{l})$ is the conditional density of $\mathbf{X}_{n}$ given $\mathcal{M}_{l}$ and $\boldsymbol{\beta}_{l}$. We assume that the prior distribution over models is non-informative, so that the BIC criterion defines the best model as

$$\mathcal{M}_{\mathrm{BIC}}=\operatorname*{arg\,max}_{{\mathcal{M}_{l}},l=1,\dots,2^{p}}\;p_{n}\left(\mathbf{X}_{n}\mid\mathcal{M}_{l}\right).$$

(16)

In principle, one could evaluate the $p_{n}\left(\mathbf{X}_{n}\mid\mathcal{M}_{l}\right)$ using numerical quadrature and then determine $\mathcal{M}_{\mathrm{BIC}}$. However, this is computationally costly and also the need to elicit a specific prior $p(\boldsymbol{\beta}_{l}|{\mathcal{M}}_{l})$ for each model may be a nuisance. Our next result therefore proposes an asymptotic expansion of $\log p_{n}\left(\mathbf{X}_{n}\mid\mathcal{M}_{l}\right)$. The methodology is standard (basically a Laplace approximation) and well-known in the literature, see Tierney and Kadane (1986) or e.g. Lebarbier and Mary-Huard (2006) and the references therein. However, due to our spatial framework and the double asymptotic point of view considered in this paper, the standard results do not apply straightforwardly. Due to the use of asymptotic results we again need to rely on the notions of a true model and the ‘least false parameter value’ $\boldsymbol{\beta}_{l,n}^{*}$.

We impose the following conditions on the prior for $\boldsymbol{\beta}_{l}$ and the mean $\mu_{n}={\mathbb{E}}\mathbf{X}_{n}$.

1. [C8]

   The prior density $p(\boldsymbol{\beta}_{l}\mid\mathcal{M}_{l})$ of $\boldsymbol{\beta}_{l}$ given ${\mathcal{M}}_{l}$ is continuously differentiable on $\mathbb{R}^{p_{l}}$.

2. [C9]

   $\mu_{n}\asymp\tau_{n}$.

Note that $\mu_{n}$ is the marginal mean of $\mathbf{X}_{n}$ under the true intensity model $\lambda_{n}$ as discussed in Section 2. Condition C9 seems reasonable since it would hold if $\lambda_{n}$ coincided with any of the specified parametric intensity models $\rho_{n}(\cdot;\boldsymbol{\beta}_{l})$. We also need to slightly strengthen assumption C1 and replace it by

1. [C1^′]

   As $n\to\infty$, there exists $\boldsymbol{\beta}_{l}^{*}\in\mathbb{R}^{p_{l}}$ such that $\lim_{n\to\infty}\boldsymbol{\beta}_{l,n}^{*}=\boldsymbol{\beta}_{l}^{*}$.

The following result is verified in Appendix D using Łapiński (2019, Theorem 2), which is a rigorous statement of a multivariate Laplace approximation.

The criterion (16) is defined entirely within the specified Bayesian framework. Hence no reference to a true model and no need for asymptotic results. However, this is changed when we derive the expansion (17) for (16). The expansion is around $\hat{\boldsymbol{\beta}}_{l,n}$ and for technical reasons, when applying the Laplace approximation, convergence of $\hat{\boldsymbol{\beta}}_{l,n}$ is needed. Therefore we need the assumptions that ensure strong consistency of $\hat{\boldsymbol{\beta}}_{l,n}$.

Since $\mu_{n}\asymp\tau_{n}$ and from the strong consistency of $\hat{\boldsymbol{\beta}}_{l,n}$, we have that $\log\mathrm{det}\{\mu_{n}^{-1}\mathbf{S}_{n}(\hat{\boldsymbol{\beta}}_{l,n})\}=\mathcal{O}_{\mathrm{P}}(1)$ while C8 ensures that $\log p(\hat{\boldsymbol{\beta}}_{l,n}\mid\mathcal{M}_{l})$ is $\mathcal{O}_{\mathrm{P}}(1)$. So, if we neglect terms which are $\mathcal{O}_{\mathrm{P}}(1)$ in (17), we follow the standard heuristic and suggest to define a first version of the BIC criterion as

$$-2\ell_{n}(\hat{\boldsymbol{\beta}}_{l,n})+p_{l}\log(\mu_{n})$$

(18)

where we remind that $p_{l}$ is the length of $\boldsymbol{\beta}_{l}$.

In practice $\mu_{n}$ is not known. However, since ${\mathbb{V}\mathrm{ar}}N(W_{n})=\mu_{n}$ it follows that $N(W_{n})/\mu_{n}-1=\mathcal{O}_{\mathrm{P}}(\mu_{n}^{-1/2})=o_{\mathrm{P}}(1)$. This justifies to define the $\mathrm{BIC}$ criterion in the following natural way

$$\mathrm{BIC}_{l}=-2\ell_{n}(\hat{\boldsymbol{\beta}}_{l,n})+p_{l}\,\log\left\{N(W_{n})\right\}.$$

(19)

Suppose $\mathbf{X}_{n}$ has an intensity function of the form (1) but is not a Poisson process. Then as mentioned in Section 3.1, (3) may be viewed as a composite likelihood score for estimating $\boldsymbol{\beta}_{l}$. In this case, following Gao and Song (2010), (16) may be viewed as a composite likelihood BIC. Again we obtain (19) from (16) by Laplace approximation. However, Gao and Song (2010) suggest to replace $p_{l}$ by the ‘effective degrees of freedom’ $p_{l}^{*}$ considered in Section 4. Thus our proposed composite likelihood BIC is

$$\mathrm{CBIC}_{l}=-2\ell_{n}(\hat{\boldsymbol{\beta}}_{l,n})+p_{l,\text{approx}}^{*}\,\log\big{\{}N(W_{n})\big{\}}$$

which becomes equal to the ordinary BIC in (19) for a Poisson process. From a practical point of view, we simply estimate $p_{l,\text{approx}}^{*}$ as in (14).

We consider two different scenarios to illustrate both types of asymptotics.

• •

  Scenario 1 (infill asymptotics). $W=[0,1]\times[0,0.5]$. We adjust $\omega=\theta\exp(\beta_{0})$ such that $\mu$ equals either $50$ or $200$.

• •

  Scenario 2 (increasing domain asymptotics). $W=[0,500]\times[0,250]$ or $[0,1000]\times[0,500]$ and $\omega$ is chosen so that $\mu=200$ or $\mu=800$.

For each scenario and the choices of $W$ and $\omega$, 500 simulations are generated from inhomogeneous Poisson point processes with intensity function given by (20) using the function $\mathtt{rpoispp}$ from the $\mathtt{spatstat}$ $\mathtt{R}$ package. We set $\beta_{1}=0.5$ and $\beta_{2}=-0.25$ to represent moderate effects of elevation and slope, and set $\beta_{3}=\dots=\beta_{6}=0$. For each simulation, parameters are estimated by maximizing the Berman-Turner approximation (see e.g.  Baddeley and Turner, 2000) of the Poisson log-likelihood (3) using a number $m$ of quadrature dummy points to approximate the integral in (3). The estimation is done using the $\mathtt{ppm}$ function. Then, the model is selected according to the AIC and BIC-type criteria. We consider different variants of the BIC criterion, namely

$$\mathrm{BIC}_{l}(\pi)=-2\ell(\boldsymbol{\hat{\beta}}_{l})+p_{l,\text{approx}}^{*}\;\log(\pi)\quad l=1,\cdots,64,$$

where $\pi$ represents a penalty. Note that (omitting dependence on $l$) $\mathrm{AIC}=\mathrm{BIC}\{\exp(2)\}$ and that $\mathrm{BIC}(N)$ ($N=N(W)$) corresponds to the criterion used by Thurman et al. (2015) and is also the criterion suggested by the present paper. We also consider $\mathrm{BIC}(|W|)$ used by Choiruddin et al. (2018) and $\mathrm{BIC}(N+m)$ considered by Jeffrey et al. (2018).

For both scenarios, we perform estimation with $m=4\mu$ (the rule of thumb suggested by spatstat) and model selection with the criteria $\mathrm{AIC}$, $\mathrm{BIC}(N)$, $\mathrm{BIC}(N+4\mu)$ and $\mathrm{BIC}(|W|)$. Similar results are obtained with $\mathrm{BIC}(N)$ and $\mathrm{BIC}(N+4\mu)$ and we omit the results for the latter. We also perform estimation and selection with $m=400\mu$ and only report results for $\mathrm{BIC}(N+400\mu)$ since the results with the criteria $\mathrm{AIC}$, $\mathrm{BIC}(N)$ and $\mathrm{BIC}(|W|)$ are similar to those obtained when the estimation is performed with $m=4\mu$.

When $|W|$ is small, especially when $\log|W|<0$, the criterion $\mathrm{BIC}(|W|)$ obviously fails as it selects the most complex model regardless of the value of $\mu$. Thereby the FP rate becomes 100%. In addition, as indicated in the second and third rows of Table 1 where the point patterns have the same average number of points, it is worth noticing that $\mathrm{BIC}(|W|)$ selects pretty different models. In particular this criterion has an undesirable strong dependence on the choice of length unit.

The criterion $\mathrm{BIC}(N+m)$ with a large $m$ also fails since the TPR with this criterion and $m=400\mu$ is very small compared to the other criteria, especially when the expected number of points is small or moderate. The AIC criterion achieves a high TPR in all situations but fails since it suffers from a high FPR, even in the scenario 2 where $\mu=800$.

In all cases, $\mathrm{BIC}(N)$ provides the best trade-off between TPR and FPR and the results improve when $\omega$ or $|W|$ is increased. The minimal values of MKL and MISE are further always obtained with $\mathrm{BIC}(N)$. In case of a Poisson process, we therefore recommend $\mathrm{BIC}(N)$ (simply denoted BIC in the following).

To generate a simulation from a Thomas point process with intensity (20), we first generate a parent point pattern from a stationary Poisson point process $\mathbf{C}$ with intensity $\kappa>0$. Given $\mathbf{C}$, clusters $\mathbf{X}_{c}$, $c\in\mathbf{C}$, are generated from inhomogeneous Poisson point processes with intensity functions

$\displaystyle\rho_{c}(u;\boldsymbol{\beta})=\omega\exp\{\beta_{1}z_{1}(u)+\cdots+\beta_{6}z_{6}(u)\}k(u-c;\gamma)/\kappa,$

where $k(u-c;\gamma)=(2\pi\gamma^{2})^{-1}\exp(-\|u-c\|^{2}/(2\gamma^{2}))$ is the density for $\mathcal{N}(0,\gamma^{2}\mathbf{I}_{2})$. Finally, $\mathbf{X}=\cup_{c\in\mathbf{C}}\mathbf{X}_{c}$ is an inhomogeneous Thomas point process with intensity (20). The regression parameters are set as follows: $\beta_{1}=2$, $\beta_{2}=-1$, $\beta_{3}=\dots=\beta_{6}=0$. We consider $\kappa=4\times 10^{-4}$ and two scale parameters $\gamma=5$ and $\gamma=15$. A lower value for $\gamma$ tends to produce more clustered patterns. We consider the observation domains $W=[0,500]\times[0,250]$ and $W=[0,1000]\times[0,500]$ with $\omega$ adjusted to give expected numbers of points $\mu=400$ and $\mu=1600$ for the two windows. The chosen value of $\kappa$ implies on average 50 parent points on $W=[0,500]\times[0,250]$ and 200 parent points on $W=[0,1000]\times[0,500]$.