Detecting multivariate interactions in spatial point patterns with Gibbs models and variable selection

TL;DR

This paper introduces a novel method combining Gibbs point process models and group lasso for detecting significant multivariate spatial interactions at multiple scales in large point pattern datasets, demonstrated on rainforest data.

Contribution

It develops a high-dimensional approach for identifying multiscale interactions in large spatial point patterns using Gibbs models and automatic variable selection.

Findings

Method effectively detects interactions at various scales.

Demonstrated on complex rainforest data with 83 species.

Shows high power and stability in identifying significant interactions.

Abstract

We propose a method for detecting significant interactions in very large multivariate spatial point patterns. This methodology develops high dimensional data understanding in the point process setting. The method is based on modelling the patterns using a flexible Gibbs point process model to directly characterise point-to-point interactions at different spatial scales. By using the Gibbs framework significant interactions can also be captured at small scales. Subsequently, the Gibbs point process is fitted using a pseudo-likelihood approximation, and we select significant interactions automatically using the group lasso penalty with this likelihood approximation. Thus we estimate the multivariate interactions stably even in this setting. We demonstrate the feasibility of the method with a simulation study and show its power by applying it to a large and complex rainforest plant population data set of 83 species.

Figures

c) we see the complex mixture of first type’s internal potential and the intra-type potential: A point of the first type would be welcome near a point of it’s own kind (attraction) but unwelcome near a point of the second kind (repulsion).

Tables

Table 1. Table 1: Experiment 1 true positive and false positive rates, mean and (sd) over 100 simulations.

	AIC0.5	CV Inverse	CV Pearson	CV Raw	SSGAM	VBSS Grouped	VBSS Ungrouped
TP	.77 (.18)	.78 (.26)	.95 (.15)	.99 (.08)	.78 (.17)	.83 (.17)	.80 (.16)
FP	.03 (.07)	.10 (.20)	.51 (.32)	.94 (.16)	.02 (.06)	.21 (.19)	.12 (.15)

Table 2. Table 2: Experiment 2 interaction detection rates, stratified by intra- and inter-type interactions and point-count per type. Mean over 50 simulations, with the standard deviation given in parenthesis. Intra-FP’s are 0 by design. Lasso penalty selected with 25-fold CV.

	$n_i=50$			$n_i=100$			$n_i=200$
Method	intra TP	inter TP	inter FP	intra TP	inter TP	inter FP	intra TP	inter TP	inter FP
raw	.51 (.20)	.97 (.04)	.04 (.06)	.65 (.21)	.99 (.02)	.03 (.05)	.89 (.10)	1.0 (.00)	.20 (.23)
inverse	.41 (.22)	.96 (.05)	.04 (.07)	.63 (.37)	.95 (.11)	.15 (.27)	.17 (.33)	.76 (.19)	.06 (.21)
Pearson	.90 (.13)	.99 (.02)	.55 (.36)	.97 (.05)	1.0 (.01)	.78 (.35)	.99 (.04)	1.0 (.00)	.90 (.27)
AIC0.5	.00 (.01)	.61 (.11)	.00 (.00)	.00 (.00)	.64 (.10)	.00 (.00)	.00 (.00)	.57 (.07)	.00 (.00)
VBSS (1,1)	.85 (.12)	.95 (.04)	.50 (.13)	.97 (.07)	.99 (.02)	.91 (.20)	.91 (.09)	1.0 (.00)	.86 (.25)
VBSS (.1,1)	.37 (.16)	.58 (.27)	.23 (.12)	.97 (.08)	.85 (.19)	.87 (.22)	.90 (.11)	1.0 (.02)	.46 (.48)
VBSS (.1,10)	.09 (.11)	.14 (.22)	.07 (.07)	.98 (.07)	.62 (.19)	.77 (.24)	.84 (.20)	.97 (.07)	.34 (.45)

Table 3. Table 3: Experiment 3, log-Gaussian Cox processes, inter-type interactions, average and sd over 50 simulations.

	corr $.6$				corr $.8$
	$n_i$=50		$n_i$=100		$n_i$=50		$n_i$=100
Method	TP	FP	TP	FP	TP	FP	TP	FP
raw	.60 (.12)	.16 (.16)	.60 (.16)	.16 (.20)	.62 (.13)	.08 (.14)	.65 (.12)	.07 (.15)
inverse	.60 (.10)	.14 (.12)	.64 (.10)	.17 (.11)	.67 (.08)	.11 (.13)	.73 (.13)	.17 (.15)
Pearson	.75 (.11)	.41 (.24)	.71 (.09)	.31 (.18)	.77 (.08)	.32 (.19)	.80 (.07)	.31 (.20)
AIC0.5	.20 (.06)	.00 (.00)	.22 (.06)	.00 (.00)	.23 (.07)	.00 (.00)	.27 (.08)	.00 (.00)

Table 4. Table 4: Experiment 4, detection rates with and without habitat effects and mostly no interaction. When type count p 𝑝 p is doubled, the inter-pair count is quadraupled ( × \times 4.22) .

	$p=10$				$p=20$
	raw	inverse	Pearson	AIC0.5	raw	inverse	Pearson	AIC0.5
	Homogeneous
FP intra	.44 (.42)	.07 (.12)	.08 (.15)	.15 (.12)	.19 (.20)	.02 (.04)	.04 (.07)	.08 (.06)
FP inter	.47 (.39)	.06 (.09)	.12 (.19)	.18 (.07)	.22 (.19)	.03 (.04)	.07 (.10)	.10 (.03)
	Homogeneous + 1 interaction
FP intra	.13 (.27)	.01 (.03)	.14 (.17)	.00 (.00)	.08 (.20)	.04 (.10)	.06 (.10)	.00 (.00)
FP inter	.07 (.13)	.00 (.01)	.13 (.16)	.00 (.00)	.06 (.17)	.02 (.06)	.06 (.07)	.00 (.00)
TP	1.0 (.00)	1.0 (.00)	1.0 (.00)	1.0 (.00)	1.0 (.00)	1.0 (.00)	1.0 (.00)	1.0 (.00)
	Inhomogeneous
FP intra	.25 (.36)	.11 (.13)	.19 (.34)	.14 (.14)	.13 (.09)	.00 (.02)	.02 (.03)	.10 (.08)
FP inter	.22 (.34)	.08 (.10)	.18 (.31)	.13 (.10)	.10 (.07)	.01 (.01)	.01 (.01)	.06 (.03)
	Inhomogeneous + 1 interaction
FP intra	.08 (.19)	.02 (.06)	.13 (.22)	.00 (.00)	.12 (.24)	.00 (.00)	.04 (.05)	.00 (.00)
FP inter	.05 (.08)	.01 (.02)	.13 (.20)	.00 (.00)	.12 (.21)	.00 (.01)	.04 (.06)	.00 (.00)
TP	1.0 (.00)	1.0 (.00)	1.0 (.00)	1.0 (.00)	1.0 (.00)	1.0 (.00)	1.0 (.00)	1.0 (.00)

Table 5. Table 5: Experiment 5, synthetic forest of independent species. Detection rates for the new approach, together with results for two Monte Carlo tests, Studentised deviance (”St’d dev.”) and rank envelope (”Rank”), with two different testing ranges. By design, all 64 intra-type interactions were non-zero, and all 2016 inter-type interactions were zero. Average and sd over 10 simulations.

		CV-raw	CV-inverse	CV-Pearson	AIC0.5
Gibbs model	TP	.76 (.34)	.91 (.06)	.96 (.08)	.35 (.06)
	FP	.01 (.01)	.05 (.10)	.20 (.12)	.00 (.00)
		St’d dev. (.5-15m)	St’d dev. (1-50m)	Rank (.5-15m)	Rank (1-50m)
Monte Carlo	TP	.74 (.05)	.89 (.03)	.80 (.05)	.96 (.01)
	FP	.02 (.00)	.00 (.00)	.03 (.00)	.02 (.00)

Table 6. Table 6: Intra- and inter-type interaction rates estimated with the proposed method using raw and inverse CV, and Monte Carlo tests using the Studentised deviance (”St’d dev.”) and rank envelopes (”rank”). Results for range-scales upto 15m and 30m separately. p = 83 𝑝 83 p=83 species, 3403 distinct pairs.

	St’d dev.	Rank	CV-raw	CV-inverse
Up to 15m
intra	0.89	0.88	0.14	0.55
inter	0.05	0.06	0.00	0.04
Up to 30m
intra	0.95	0.96	0.33	0.60
inter	0.02	0.05	0.00	0.06

Equations

$$f(\mathbf{x})\propto \exp [\sum _{i,j=1}^p{\varphi }_{ij}({\mathbf{x}}_i,{\mathbf{x}}_j)],$$

$${\int }_{\mathcal{X}}\exp [\sum _{i,j=1}^p{\varphi }_{ij}({\mathbf{x}}_i,{\mathbf{x}}_j)]{\mu }_1(d\mathbf{x}),$$

$${\varphi }_{ij}({\mathbf{x}}_i,{\mathbf{x}}_j)$$

$$f(\mathbf{x})=f_{\theta }(\mathbf{x})\propto \exp [{\theta }^T\mathbf{v}],$$

$${\mathbf{g}}_{ij}(x,\mathbf{x}\setminus x)=\sum _{y\in \mathbf{x}\setminus x}{\psi }_{ij}(||x-y||),$$

$$g(x,\mathbf{x}\setminus x)$$

$$n{e}_{ijk}(x,\mathbf{x}):=\#\{\mathbf{x}\cap [b(x,r_{ijk})\setminus b(x,r_{ij(k-1)})]\},$$

$$g_{ijk}(x,\mathbf{x})$$

$${\lambda }_{\theta }(u;\mathbf{x})$$

$$\log {\lambda }_{\theta }(u;\mathbf{x})$$

$${\omega }_{ijk}((x,i),{\mathbf{x}}_j)=g_{ijk}(x,{\mathbf{x}}_i\setminus x)+\sum _{y\in {\mathbf{x}}_j}[g_{ijk}(y,{\mathbf{x}}_i\cup x)-g_{ijk}(y,{\mathbf{x}}_i\setminus x)].$$

$$t(c):=c[1-F_a(c-1)]+a[F_a(c-1)+F_a(c-2)],$$

$$1-t(c)/2a<ϵ$$

$${\hat{\theta }}_{\gamma }={\text{argmax}}_{\theta }{L}_{\gamma }(\mathbf{x};\theta ),$$

$$L_{\gamma }(\mathbf{x};\theta ):=\log {\tilde{f}}_{\theta }(\mathbf{x})-\gamma \sum _{g=1}^m|(g){|}^{-\frac{1}{2}}\Vert {\theta }_{(g)}{\Vert }_2,$$

$${\hat{\theta }}_{\gamma ,-k}={\text{argmax}}_{\theta }{L}_{\gamma }(\mathbf{x}\setminus W_k;\theta )$$

$${\hat{R}}_{\gamma ,k}$$

$${\hat{R}}_{CV}(\gamma ):=\frac{1}{K}\sum _{k=1}^K{\hat{R}}_{\gamma ,k}^2$$

$${\hat{\gamma }}_{CV}:={\text{argmin}}_{\gamma \in \Gamma }{\hat{R}}_{CV}(\gamma ),$$

$$M=[M_{ij}]_{1}^{p},\quad\text{where}\quad M_{ij}=\left\{\begin{array}[]{rl}e_{g(ij)}&i\neq j\\ e_{g(i)}&i=j\\ \end{array}\right..$$

$$s_W(X,D;\theta ):=\sum _{u\in X\cap W}\frac{\rho (u)\mathbf{v}(u;X\setminus u)}{{\lambda }_{\theta }(u;X\setminus u)+\rho (u)}-\sum _{u\in D\cap W}\frac{\mathbf{v}(u;X){\lambda }_{\theta }(u;X)}{{\lambda }_{\theta }(u;X)+\rho (u)},$$

$$\log {\tilde{f}}_{\theta }(X\cap W)=\sum _{u\in X\cap W}\log \frac{{\lambda }_{\theta }(u;X\setminus u)}{{\lambda }_{\theta }(u;X\setminus u)+\rho (u)}+\sum _{u\in D\cap W}\log \frac{\rho (u)}{{\lambda }_{\theta }(u;X)+\rho (u)},$$

$$P(\tau (u)=1)=\frac{{\lambda }_{\theta }(u;X\setminus u)}{{\lambda }_{\theta }(u;X\setminus u)+\rho (u)}=\frac{\exp [{\theta }^T\mathbf{v}(u;X)+\log \rho (u{)}^{-1}]}{1+\exp [{\theta }^T\mathbf{v}(u;X)+\log \rho (u{)}^{-1}]}.$$

$$\log {\tilde{f}}_{\theta }(\mathbf{x})$$

$$\omega (u,\mathbf{x})$$

$$=min(c,\#[b(u,r)\cap \mathbf{x}])+\sum _{x\in \mathbf{x}}1[x\in b(u,r)]1\{\#[b(x,r)\cap \mathbf{x}]\in [0,c)\}$$

$$\mathbb{E}\omega (u,X)$$

$$=\mathbb{E}\min (c,\#[b(o,r)\cap X])+a\mathbb{E}1\{\#[b(o,r)\cap X]\in [0,c)\}$$

Full text

1 Introduction

Spatial point patterns are a common form of observation in plant ecology (Waagepetersen et al., 2016), epidemiology (Diggle et al., 2005), astrophysics (Stoica et al., 2007), seismology (Schoenberg, 2003), social science (Amburgey, 1986), medicine (Olsbo et al., 2013) and criminology (Mohler et al., 2011). While understanding single, univariate spatial point patterns and their generating point processes is important, frequently we observe labelled point processes, or more precisely, multiple types of points. The prevalence of multivariate point processes is particularly noticeable in plant ecology where there may be many tens or hundreds of types (species)  (Flügge et al., 2014; Baldeck et al., 2013b, a; Kanagaraj et al., 2011; Punchi-Manage et al., 2013). Such processes have seen much less study in the statistical literature than univariate processes, and present a number of novel challenges, as we shall explain and address in this paper.

To be able to make sense of multivariate point processes, we focus on addressing three important outstanding problems in understanding interactions: i) characterizing patterns associated with both small and large scales simultaneously, ii) characterizing multiple features spanning more than one variable, and iii) estimating such patterns stably. The effect of including these characteristics in our studies of a multivariate point process greatly increases the number of potential parameters to characterise and naturally leads to difficulties in arriving at unique and stable parameter estimates.

Early approaches to modelling and analysing more than one point processes focused on bivariate representation and analysis (Diggle and Milne, 1983; Amburgey, 1986; Brix and Møller, 2001; Gelfand et al., 2004; Shimatani, 2001; Diggle et al., 2005), and going much beyond the bivariate case has proved very challenging. To the best of our knowledge the most diverse model-based analysis of multivariate point patterns to date is of nine rainforest tree species by Waagepetersen et al. (2016). These authors used a multivariate log-Gaussian Cox process model, and fitted the model using non-parametric least squares. The reason given for confining studies to only nine species by the previous authors was to limit the computational burden of analysis. This work is inspirational, but our motivating problem in this paper is to jointly analyse an order of magnitude more species. More specifically, there are 300 species (types) in the full dataset used in  Waagepetersen et al. (2016), and we wish to extend analysis to investigate as many of these species as possible. This brings us into the realm of high dimensional statistics as the number of interactions scales exponentially in the number of species, whilst the number of points only scales linearly in the number of species. To deal with this inconvenient scaling we shall need to use shrinkage, as is commonly done in high dimensional data analysis, and has already been developed for regression problems and covariance estimation, see e.g.  van der Geer and Bühlmann (2011). In the context of point processes, because estimation is not implemented with linear methods, penalization needs to be deployed carefully. Starting from ideas of Baddeley et al. (2014) we shall use a generalized regression–based fitting approach, and so we borrow ideas from second generation penalized regression rather than from matrix shrinkage, even if we are estimating co-associations rather than a mean intensity.

A second significant problem in point process modelling is proposing models, and associated estimation methods, that yield sufficient multiscale behaviour, e.g. variability at fine scales, as well as over medium to long scales. It is all very well to posit variability for fine scales using a log-Gaussian Cox process, but as estimation is normally based on some form of averaging, unless the random intensity is very high locally, it will be impossible to estimate the log-Gaussian Cox process’ generating mechanism as we will not have enough points. For longer spatial scales the log-Gaussian Cox process is a well-suited modelling framework, but it is not a good framework for studying small-scale interactions. Instead we shall use the multivariate Gibbs point process model to discover small scale point-to-point interactions in the same Barro Colorado Island (BCI) rainforest dataset studied by Waagepetersen et al. (2016).

Key to our modelling and estimation is therefore capturing an appropriate degree of sample heterogeneity. The estimation framework we shall introduce is able to take into account variations associated with i) habitat associations, i.e. correlation of species presence in the landscape with known environmental covariates, with ii) dispersal mechanisms and competition such as seedling clustering or self-thinning, and with iii) attraction and repulsion between the small scale locations between different species. To demonstrate utility, we will fit the model to an adult plant community consisting of $83$ species, an order of magnitude more species than Waagepetersen et al. (2016).

The analysis of multivariate point patterns of more than a handful of species (usually two, e.g. Brix and Møller, 2001; Diggle and Milne, 1983; Högmander and Särkkä, 1999; Funwi-Gabga and Mateu, 2012) has mostly relied on non-parametric estimation techniques. Many ecological analyses of large and diverse rainforest datasets have been carried out using the $K$-function or similar non-parametric summaries, either directly by comparing the summary values under various null-model scenarios, or indirectly as part of minimum contrast model fitting of Cox processes (e.g. Lan et al., 2012; Flügge et al., 2014; Waagepetersen et al., 2016; Yang et al., 2016; Velázquez et al., 2016; Brown et al., 2016). In the Gibbs-model framework, the three-variate analysis of different size trees by Grabarnik and Särkkä (2009) is perhaps the most extensive modelling approach and closest in spirit to the work of this paper. We will extend their trivariate case to a full multi-scale multivariate interaction Gibbs model.

Direct likelihood inference for Gibbs models is not available due to intractable normalising constants, but several pseudo-likelihood approximations are available, such as the “Berman-Turner machine” (Baddeley and Turner, 2000) used by Grabarnik and Särkkä (2009) which casts the pseudo-likelihood estimation equation as a Poisson regression problem, and subsequently estimates the model with standard statistical software. The Poisson regression approach is surpassed in accuracy by the logistic regression approach developed by Baddeley et al. (2014), who formulates the pseudo-likelihood estimation equation as a logistic regression using auxiliary dummy point configurations. Again, very conveniently, standard statistical software can be used to fit the model.

With a model and the likelihood approximation at hand, we will next tackle the issue of high dimensional variable selection: Any reasonable model for a $p$-variate point pattern will have a high number of parameters when $p\gg 3$, scaling at least like ${\cal O}(p^{2})$. As an illustration, the model we present estimates intra-species interactions and pairwise species-to-species interactions, on 3 different spatial scales, and includes 6 covariates, giving a total of $\approx 11,000$ parameters, to be compared with the number of observations of $31,650$ points, making the number of parameters and observations on the same order. To discover significant interactions in such a high dimensional setting, we will use recent research on penalized optimization. Several techniques are available, and to show proof of concept we will be using the group lasso (Yuan and Lin, 2006; Meier et al., 2008). We will also perform a limited comparison to the Bayesian spike and slab variable selection approach (Mitchell and Beauchamp, 1988). Using penalized regression (the group lasso) rather than full posterior inference to compute a MAP yields both simplicity in interpretation and computational speed. This allows us to study several “priors”, or penalization choices, at once and thus this allows us to make fewer assumptions on the generating mechanism of the data, in fact letting us explore the properties of our modelling framework.

We start in Section 2 by introducing the Gibbs model, then recall the chosen likelihood approximation and discuss the chosen variable selection techniques in detail. In Section 3 we recall for comparison the non-parametric Monte Carlo technique often used for analysing $p$-variate patterns when $p\gg 1$. In Section 4 we test the method on several increasingly complex simulation scenarios, to get a better understanding of the performance of the method. In Section 5 we apply the method to the Barro Colorado Island rainforest data (Condit, 1998), that inspired these developments. We conclude in Section 6, and discuss outstanding problems and future avenues of investigation.

2 The Gibbs model and Gibbs model fitting with variable selection

Let the observed multivariate point pattern be a set of labelled point locations $\mathbf{x}=\{(x,t)\}$ where $x\in W$ are the observed locations inside a known, bounded observation window $W\subset\mathbb{R}^{d}$, $d=1,2,3,...$ (for the problems we shall study, $d=2$), and $t\in\{1,...,p\}$, $p\geq 1$, are categorical labels (types, or in our case, species) attached to each location. Denote the type $i$ sub-pattern by $\mathbf{x}_{i}=\{x:(x,i)\in\mathbf{x}\}$. Write $n_{i}(A)=\#(\mathbf{x}_{i}\cap A)$ for the count of points of type $i$ in a set $A\subset\mathbb{R}^{d}$, and write $n_{i}=n_{i}(W)$. As is standard let $b(u,r)$ denote a ball of radius $r$ centred at $u\in\mathbb{R}^{d}$.

2.1 Gibbs models for point patterns

Our model for multivariate spatial point patterns is part of the Gibbs point process model family; we refer to van Lieshout (2000) who details general properties of Markov Point Processes (where Gibbs point processes are a special case), and Illian et al. (2008) and Chiu et al. (2013) who discuss general point processes. We will use a Gibbs model, defined from a set of potential functions $\phi_{ij}$, to have a probability density of the form

[TABLE]

with respect to the unit rate Poisson process $\mu_{1}$ on $W$. The normalising constant for the density,

[TABLE]

with $\mathcal{X}$ the space of all locally finite multitype point patterns, is in practice intractable for all but the homogeneous Poisson models where $\phi_{ij}\equiv\alpha_{i}\in\mathbb{R}$. The functions $\phi_{ij}$ are used to specify the model class member. We will use the special form of

[TABLE]

In (4) the parameters $\alpha_{i}\in\mathbb{R}^{K_{i}+1}$ regulate intensity and covariate effects (so-called first order effects), and $\beta_{ij}\in\mathbb{R}^{K_{ij}}$ are parameters for interactions between the locations (second order effects). The $\beta_{ij}$’s give the magnitudes of the interactions, whereas the vector-valued functions $\mathbf{g}_{ij}$, which we specify later in Section 2.2, determine the form, spatial scales and orders of the interactions. The covariate effects $\mathbf{z}_{i}(u)=[1\ z_{i1}(u)\ z_{i2}(u)\ldots z_{iK_{i}}(u)]^{T}$ represent the baseline effect and any covariate and trend effect values we have at locations $u\in W$. In this formulation we assume that the covariates are available everywhere in the window $W$. This is usually achieved by interpolation from prior data collection efforts.

To get a heuristic understanding of the model, assume first that all parameters except $\alpha_{i}\in\mathbb{R}$ are zero and $\mathbf{z}_{i}\equiv 1$. Then the density (1) becomes $\exp[\sum_{i=1}^{p}n_{i}\alpha_{i}]=\prod_{i=1}^{p}e^{n_{i}\alpha_{i}}$, which is the likelihood of a collection of $p$ independent homogeneous Poisson processes. Now add some non-constant covariates $z_{i2}(u),z_{i3}(u),\ldots$ and set $\alpha_{i2},\alpha_{i3},\ldots\neq 0$: The model becomes a collection of independent and inhomogeneous Poisson processes. We subsequently add intra-type interaction terms by letting $\beta_{ii}\neq 0$: The independent components are no longer Poisson, but exhibit internal, within type point-to-point interactions (attraction or repulsion depending on the sign of $\beta_{ii}$). Finally we can add inter-type interaction terms by setting $\beta_{ij}\neq 0$: The locations of different types are no longer independent.

The model of (1) is log-linear in parameters $\alpha$ and $\beta$. We collect the covariate parameters in the vector $\theta_{0}=[\alpha_{1}^{T}\ldots\alpha_{p}^{T}]^{T}$, the intra-type interaction parameters to $\theta_{1}=[\beta_{11}^{T}\ldots\beta_{pp}^{T}]^{T}$ and inter-type interaction parameters, assuming for now that $i$ and $j$ interact symmetrically, to $\theta_{2}=[\beta_{12}^{T}\ldots\beta_{(p-1)p}^{T}]^{T}$. Subsequently we collect them in the vector $\theta=[\theta_{0}^{T}\ \theta_{1}^{T}\ \theta_{2}^{T}]^{T}$. Then the density in (1) can be written in the form

[TABLE]

where the vector $\mathbf{v}=\mathbf{v}(\mathbf{x})=[\mathbf{s}_{0}^{T}\ \mathbf{s}_{1}^{T}\ \mathbf{s}_{2}^{T}]^{T}$ has components

• •

$\mathbf{s}_{0}=[\mathbf{s}_{1}^{T}\ldots\mathbf{s}_{p}^{T}]^{T}$ with $\mathbf{s}_{i}=\sum_{x\in\mathbf{x}_{i}}\mathbf{z}_{i}(x)$

• •

$\mathbf{s}_{1}=[\mathbf{s}_{11}^{T}\ldots\mathbf{s}_{pp}^{T}]^{T}$ with $\mathbf{s}_{ii}=\sum_{x\in\mathbf{x}_{i}}\mathbf{g}_{ii}(x,\mathbf{x}_{i}\setminus x)$

• •

$\mathbf{s}_{2}=[\mathbf{s}_{12}^{T}\ldots\mathbf{s}_{(p-1)p}^{T}]^{T}$ with $\mathbf{s}_{ij}=\sum_{x\in\mathbf{x}_{i}}\mathbf{g}_{ij}(x,\mathbf{x}_{j})$, $i<j$.

The vector $\mathbf{v}=\mathbf{v}(\mathbf{x})$ is written as a matrix $[\mathbf{v}^{T}]$, and takes the role of a design matrix in the standard regression setting. The matrix has only one row when the point process is observed only once. Independent replicates would result in additional rows in the matrix.

Given data $\mathbf{x}$, covariates $\mathbf{z}$ and the interaction functions $\mathbf{g}_{ij}$, the vector $\mathbf{v}$ is fixed and the non-normalised model (5) is log-linear in the unknown coefficients $\theta$. The model therefore belongs to the family of exponential Gibbs models, and we can apply inference techniques designed for the exponential Gibbs family.

2.2 The interaction functions

We now define the exact model we will use in our examples. Several definitions are available for the interaction functions $\mathbf{g}_{ij}$ in (4). We assume $\mathbf{g}_{ij}$ are non-negative functions to remove sign ambiguity when estimating $\beta_{ij}$. The most popular class of models is the pairwise interacting models with

[TABLE]

for some functions $\psi_{ij}:\mathbb{R}_{+}\mapsto\mathbb{R}_{+}^{K_{ij}}$. The most common choice in the case of $K_{ij}=1$, is the Strauss model

[TABLE]

effectively counting the number of $r$-close pairs in the pattern. In the univariate case the Strauss model is only valid when the corresponding interaction coefficient $\beta<0$ so that fewer point pairs leads to a higher likelihood. Trying to model positive interactions, or clustering, leads to an unstable model that produces patterns of singular mega-clusters, so the case $\beta>0$ is excluded for the simple Strauss model (Gates and Westcott, 1986).

To circumvent this limitation Geyer (1999) introduced a model he called the saturation model that still defines a locally stable process even with a positive interaction parameter $\beta>0$. In the simple univariate case the saturation model is defined via

[TABLE]

where the range $r$ is the reach of the Euclidean neighbourhood, and $c$ is a saturation level. In this model, each point contributes to the likelihood a factor relative to the number of $r$-neighbours or $c$, which ever is smaller (hence the saturation). In ecological terms, the Geyer model is able to capture the fact that individuals may cluster at some distances, but are likely to segregate at shorter distances due to intense competition, and the saturation parameter reproduces the feature that the neighbourhood must eventually saturate with individuals as resources are finite. The model belongs to a class of interacting neighbour models (Grabarnik and Särkkä, 2001), so-named because the conditional intensity (10) for this model depends not just on the local neighbourhood of a point in $u$ (which it does for pairwise models), but also on the neighbourhoods of the neighbours of $u$.

To model several types of points and more than one spatial scale, we generalise the models in two ways by adding (i) multiple ranges, and (ii) cross-type interactions when $p\geq 2$. Let $\mathbf{r}_{ij}=\{r_{ijk}:0=r_{ij0}<r_{ij1}<...<r_{ijK_{ij}}\}$ be a fixed increasing vector of ranges for $i\leq j$. Let $c_{ij}=\{c_{ijk}\in\mathbb{N}\}$ be the saturation parameters. For a given range vector, write

[TABLE]

for the number of neighbours of $\mathbf{x}$ a point $x$ has in the annulus between ranges $r_{ij(k-1)}$ and $r_{ijk}$. Then the multi-step multi-type Strauss model is defined using the interaction function made of components

[TABLE]

and the multi-step multi-type saturation model of components

[TABLE]

in definitions of $\mathbf{g}_{ij}$. Note that the Strauss model is a special case of the saturation model with $c\rightarrow\infty$.

Several other forms of $g_{ijk}$ can be used, and multiple forms can be combined as described in Baddeley et al. (2013). We will not pursue them here as either a) they can be approximated by the Strauss or saturation model as the interaction functions $\beta_{ij}^{T}\mathbf{g}_{ij}$ are step-functions over spatial scales; b) there is not sufficient data available to estimate very fine details over many spatial scales; or c) they are computationally costly (e.g. morphological functions). However, if forms such as the area interaction model (Baddeley and Lieshout, 1995) seem more appropriate for a specific application, the proposed framework is still valid and can be adapted to be used in this setting.

For the purpose of fitting the model (Section 2.3), we need to define the conditional (or Papangelou) intensity of the model: At any point $u=(x,i)\in\mathbb{R}^{d}\times\{1,...,p\}$ let

[TABLE]

with $0/0:=0$. Heuristically, $\lambda_{\theta}(u;\mathbf{x})du$ can be understood as the conditional probability of observing a point $u$, given the rest of the pattern $\mathbf{x}$. For the exponential family Gibbs models the conditional intensity has the rather simple form of

[TABLE]

where we use the notation $\mathbf{v}(u;\mathbf{x})=\mathbf{v}(\mathbf{x}\cup u)-\mathbf{v}(\mathbf{x}\setminus u)$. Notice that the intractable normalising constant cancels out in this expression.

The conditional intensity (10) with the step-wise components (8 and 9) is at any marked point $u=(x,i)$

[TABLE]

with

[TABLE]

The functions $\omega_{ijk}\in[-c_{ijk},n_{i}c_{ijk}]$ track the changes in the neighbourhood inclusion counts within and between different types. An illustration of the model is given in Appendix C.

For both models the non-canonical parameters $r_{ijk}$ and $c_{ijk}$ cannot be directly inferred by the pseudo-likelihood methods. Thus the parameters need to be fixed as part of the model definition. In data analysis one usually has a priori information on relevant ranges $\{r_{ijk}\}$ (e.g. Uriarte et al. 2004), but the saturation level $\{c_{ijk}\}$ is harder to set. To reduce this complexity we propose to choose the saturation levels $\{c_{ijk}\}$ automatically depending on the abundances $n_{j}$. The idea is that if the abundance $n_{j}$ is high the neighbourhoods $ne_{ijk}$ and therefore $\omega_{ijk}$ saturate often even under the independence assumption. Under full independent assumption, the expectation of $\omega_{ijk}$ is a function of $c=c_{ijk}$

[TABLE]

where $a=a_{ijk}=|b(o,r_{ijk})\setminus b(o,r_{ij(k-1)})|n_{j}/|W|$ and $F_{a}$ is the cumulative distribution function of a $Poisson(a)$ random variable (cf. Appendix B). The function $t$ is non-decreasing and $t\rightarrow 2a$ when $c\rightarrow\infty$, so that $1-t(c)/2a$ is the cdf of the event that saturation occurs under independence. In order to avoid saturation due to high abundances alone, a sensible choice of $c$ is a value for which saturation under independence is unlikely. Therefore, we set a small $0<\epsilon<1$ such that

[TABLE]

and, by further using the approximation $t(c)/2a\approx F_{a}(c-1)$, we choose $c=c_{ijk}$ to be the $(1-\epsilon)$-quantile of $F_{a}$. In the examples we use $\epsilon=.01$. Note that with the modification the interaction functions $\mathbf{g}_{ij}$’s become asymmetric, but we will treat the $\beta_{ij}$’s symmetrically in our examples for simplicity.

2.3 Inference: Approximating the likelihood at its mode

The likelihood in Equation (1 equiv. 5) is not computationally tractable due to an unspecified normalising constant. In order to carry out standard likelihood inference the constant can be approximated by Monte Carlo techniques, but these tend to be computationally costly even for the univariate case. The more commonly used approach is to use pseudo-likelihood techniques, which replace the function to maximise with something that approximates the likelihood at its mode. We will use the recent developments proposed by Baddeley et al. (2014), which in practice conveniently leads to a logistic regression formulation.

We have summarised the details of the method in Appendix A. For the purposes of the discussion, it suffices to know that the pseudo-likelihood function $\tilde{f}_{\theta}$ is formally a likelihood of a logistic regression function. It involves additional sets of random dummy points per type, of which construction, particularly their intensity $\rho_{i}$, is an additional user decision. Baddeley et al. (2014) discuss several potential options to be used for the dummy distributions and intensities. We will use the recommended homogeneous stratified uniform distributions and, if not otherwise stated, intensities that are four times the intensity of data.

To address boundary effects due to censoring near the edges of observation window $W$, we will exclude components of $\tilde{f}_{\theta}$, say $\mathbf{b}(u)$, for which the distance from $u$ to the border of $W$ is less than a range $r_{bor}>0$, i.e. in the sum in (Appendix A eq. A) the $W$ is replaced by $W\ominus b(o,r_{bor})$. The range $r_{bor}$ is taken to be the maximal interaction range in the model, and is determined by the $\mathbf{g}_{ij}$’s.

2.4 Penalised inference for grouped coefficients

Our main goal is the detection and estimation of significant within-type and between-type interactions in multivariate point pattern data. To this end the previous section allows us to estimate groups of coefficients, such as the inter-type interaction vectors $\beta_{ij}\in\mathbb{R}^{K_{ij}}$ for type pairs $i\neq j$. Subsequent to that, we need to determine if $\beta_{ij}=0$ or if this is not the case.

In a general setting, let the length of an unknown coefficient vector $\theta$ be $M$, and split $\theta$ into $m$ smaller groups using a partition given by $\{\pi_{1},...,\pi_{m}\}$ of $\{1,...,M\}$. For shorthand we will write $(g)=\pi_{g}$, where $g=1,...,m$. Write $e_{g}=0$ when $\theta_{(g)}=0$ and $e_{g}=1$ when at least one of the coefficients $\theta_{(g)}$ is non-zero. In our case the partition is given by the types and type-pairs, and the $e_{g}$’s are connected to the events $\beta_{ij}\neq 0$. Our task is to determine which $e_{g}$’s are non-zero in a given data $\mathbf{x}$, a task known as grouped variable selection (Yuan and Lin, 2006).

Variable selection and shrinkage in high dimension, usually implemented using penalised optimisation, has become an important computational statistical technique due to an increase in high data-throughput applications such as genomics and text analytics. Variable selection has already been applied to point pattern analysis, see for example the work by Yue and Loh (2014). The latter authors applied variable selection in the univariate case for covariate selection. Several methods exist for the particular problem of group-level selection (or penalisation) that we are addressing in this context (see e.g. Breheny and Huang, 2009). From the available selection of optimization criteria, we picked the group lasso for logistic regression (Meier et al., 2008), as a suitable penalized version of Baddeley et al. (2014). The group lasso is an extension of the original $\ell_{1}$-penalisation by Tibshirani (1996) for individual coefficients, and has been further extended to a mixed level penalisation by Simon et al. (2013). To be concrete, the group lasso estimator is defined by

[TABLE]

with the group-penalised likelihood function

[TABLE]

where $\gamma>0$ is a penalisation parameter and $\|\cdot\|_{2}$ is the Euclidean distance. Some groups $\theta_{(g)}$ are shrunk to exactly 0, depending on the strength of the penalisation, so the group lasso does provide us with binary indicators $e_{g}=1(\hat{\theta}_{(g)}\neq 0)$.

The penalisation parameter $\gamma$ is directly linked to the amount of non-zero $e_{g}$’s that the algorithm outputs and it needs to be chosen by the user. With some extra computational cost, we compute what is known as the lasso-path, a set of estimates for each value of $\gamma\in\Gamma=[0,\gamma_{max}]$. The maximal penalty $\gamma_{max}$ is the penalty level below which the first penalised group is let in to the model. In practice, we will use by default a 100 step log-linearly decreasing equidistant grid from $\gamma_{max}$ to $0.001\gamma_{max}$.

Typically when using a lasso-based variable selection analysis one would use Cross Validation (CV) to choose the level of penalization that gives the best balance between model fit and prediction quality, or minimises the expected risk (Giraud, 2014, Sec. 5.2), also called the extra-sample error (Hastie et al., 2001, Sec. 7.1.). In the point pattern context, Yue and Loh (2014) also chose their penalty using CV, but unfortunately no details were given how the data was partitioned or what error metric was used for the prediction. Conducting CV by splitting the constructed data frame that was used in the logistic regression part of the inference is not justified as the rows are dependent due to spatial correlation. This flaw leads us to over-estimate the complexity of the model. Note also that it has been observed that as cross-validation is aimed at prediction, rather than model selection, and so too many variables are often retained (van der Geer and Bühlmann, 2011, Section 2.5.1).

Let us define a CV procedure for spatial point pattern models that coincides with the CV for classical regression. Two elements are needed; a way to split the data into meaningful subsets of independent data; and a way to measure prediction error, or risk, with the already explained caveats. Let the observation window $W$ be partitioned into disjoint quadrats $W_{1},...,W_{K}$, and write for each $k=1,...,K$

[TABLE]

for the group lasso estimate using all the data except points inside $W_{k}$. For the CV risk we need to define prediction quality, for example the residual sum of squares used in linear regression. Baddeley et al. (2004) introduced the $h$-residual measure for point patterns which we can use to define a residual sum of squares. Define the CV $h$-residual as

[TABLE]

where $h$ is a non-negative function called the test function. Baddeley et al. (2004) and Coeurjolly and Lavancier (2013) list four options for the test function, of which we will look at three. The raw residuals with $h=1$ assess only the trend part of the model, and the inverse residuals with $h(u;\mathbf{x},\theta)=\lambda_{\theta}(u;\mathbf{x})^{-1}$ assess the interactions. The Pearson residuals with $h(u;\mathbf{x},\theta)=\lambda_{\theta}(u;\mathbf{x})^{-1/2}$ is similar to the inverse but with variance that is in theory independent of $\lambda$ for the Poisson process.

We define the $K$-fold CV estimator of the prediction risk as the mean of the squared residuals

[TABLE]

and we choose the minimizer of the estimated risk,

[TABLE]

as the cross-validated penalty level. Note that one can easily weight the residuals by relative importance of each type or the quadrat size.

Figure 8 in the appendix illustrates this approach for an example pattern generated by Experiment 1. The individual quadrats’ residuals vary—this is not unexpected as the number of data points per quadrat is small—but the minimum average risk still leads to a reasonable penalisation.

In a CV procedure the model needs to be fitted $K$ times so computational cost and CV risk estimate stability are to be balanced. In the spatial setting an additional problem arises due to border correction, viz. we use $W_{k}\ominus b(o,r_{bor})$ instead of $W_{k}$ in (2.4) with some $r_{bor}>0$. Since each sub-window needs to be reduced by the same border correction range as the original window to get truly independent subsets of data, the loss of data for estimating $\hat{R}_{\gamma,k}$’s limits the number of splits that can be done on $W$, as some data is simply lost in the process. For example, in a square window, a border correction range that is 5% of the window’s side together with a $3\times 3$ partitioning will effectively lead to a 50% loss of data when estimating the risk (cf. Appendix Figure 9). Furthermore, the varying abundances of the types need weighting to account for importances of different types, and heterogeneity of data should be considered as well. It is therefore very hard to give general guidelines for partitioning the window that would work in all scenarios.

For comparison with the CV penalty selection, we also include a rule-of-thumb penalty selection based on the AIC. We computed the AIC for the group lasso as described by Breheny and Huang (2009). In our experiments, penalties with the lowest AIC, say $\gamma_{AIC}$, consistently led to too dense solutions. This suggests that either the pseudo-likelihood approximation or the spatial dependency, or both, lead to underestimation of the effective number of parameters. As a plug-in rule of thumb, initial trials indicate that a penalty around $(\gamma_{AIC}+\gamma_{max})/2$ leads to reasonable penalisation; we will report these results with the label AIC0.5.

We also constructed a Bayesian algorithm for computing approximate posteriors $Pr(e_{g}=1|\mathbf{x})$. We call this tailored algorithm VBSS (Variational Bayes Spike and Slab). It uses a quadratic approximation to the logistic function (Jaakkola and Jordan, 2000), selects variables via the spike and slab priors (Mitchell and Beauchamp, 1988) and uses the mean-field variational approximation to infer the maximum-a-posteriors (Ormerod and Wand, 2010). In Experiment 1 we also checked a non-grouped version of VBSS and compared it to the MCMC based SSGAM (Scheipl, 2011) that has a slightly different implementation of the spike and slab priors. All of these alternatives solve the same optimisation problem as group lasso (Equation 13), but with different choices of penalisation.

3 Alternative approach for multi-type interaction detection: Non-parametric Monte Carlo

The current state-of-the-art methodology for analysing high dimensional multivariate point pattern data is based on non-parametric summary statistics and Monte Carlo (MC) testing. This approach splits the analysis into a set of bivariate tests for no co-associations (second order interactions). As a reference to our proposed method we apply a MC testing scheme in the BCI rain forest -like, large inhomogeneous simulation experiment of Section 4.5 and the actual BCI rain forest data analysis Section 5. We give here a short description of MC testing but the interested reader should seek more in-depth texts on the topic see e.g. Illian et al. (2008) Sec. 7.5; Baddeley et al. (2014); Velázquez et al. (2016); Brown et al. (2016).

Monte Carlo -based testing procedures in general, and tests of spatial type-to-type interactions in particular, are built on three components: 1) A null hypothesis that describes a counter factual pattern without the interactions of interest. 2) A statistical summary that measures the interaction of interest. 3) A statistical test based on the summary for measuring data’s departure from the distribution of the summary when the null hypothesis holds.

In the context of potentially inhomogeneous patterns, a popular null model, especially in ecology, is some form of the inhomogeneous Poisson process (IPP) (Wiegand et al., 2012). A common, asymmetric version can be simulated by keeping type $i$ locations fixed while sampling type $j$ locations from IPP independently from all other data. The idea is to nullify short range interactions but keep longer range, environment-related associations via the spatially-varying intensity. The intensity surface needed for the IPP simulations of $\mathbf{x}_{j}$ is estimated with a fixed smoothness that reflects understanding of the division between ”short” and ”long” ranges. We implement the IPP simulations following Wiegand et al. (2012): For each type $i=1,...,p$, an intensity field $\eta_{i}(u),u\in W=[0,1000]\times[0,500]m$ is kernel-estimated on a $2\times 2m$ grid using the border-corrected Epanechnikov kernel with bandwidth fixed to $30m$ (this is of course not optimal for all types; we simply emulate Wiegand et al. 2012). Then $n_{i}$ points are distributed on $W$ with density relative to $\eta_{i}(u)$, i.e. we use a conditional IPP as is commonly done. We simulate 999 patterns this way for each type $i=1,...,p$.

For the statistical summary of pattern interaction, a popular choice is the cross-Ripley’s $K_{ij}$-function. It describes the number of type $j$ neighbours for an average point of type $i$, over different neighbourhood ranges. So, assuming isotropy, for each pair of bi-variate patterns the $K_{ij}(r)$-function estimate is a sampled curve over spatial scales $r>0$. In our analyses, we estimate the cross-type $K_{ij}(r)$ curves on a range grid $r=r_{1},...,r_{max}$ with translation edge correction. For each species $i=1,...,p$ and $j=1,...,p$, we estimate the curve $K_{ij}^{0}=\{K_{ij}^{0}(r_{1}),...,K_{ij}^{0}(r_{max})\}$ from the bivariate data pattern $\mathbf{x}_{i}\cup\mathbf{x}_{j}$. When $i=j$ we estimate the univariate Ripley’s $K_{ii}=K$. Subsequently, we estimate the curves $K_{ij}^{b}$ from the bivariate synthetic patterns $\mathbf{x}_{i}\cup\mathbf{x}_{j}^{b}$ where $\mathbf{x}_{j}^{b}$ is a simulation of the null model for $j$ as described above, and $b=1,...,999$. With the set of curves $(K_{ij}^{0},...,K_{ij}^{999})$ for all $p^{2}$ combinations of $i$ and $j$, we do the variance stabilising $\sqrt{K(r)/\pi}$ transform to increase statistical power.

The third component of the MC testing framework is then needed for a proper combined test for the $K_{ij}$ curves, i.e. to determine if $K^{0}_{ij}$ is different from the null model curves. A family of such tests are called deviation tests or envelope tests (Myllymäki et al., 2017; Baddeley et al., 2014), and several options are available. We will use the Studentised deviation test, which measures the $L^{2}$-distance of the Studentised curves (scaled with respect to the null model), and the rank envelope test, which is a multidimensional analogue of a rank-test (cf. Myllymäki et al. 2017). Each test leads to a $p$-value per test, say $p_{ij}$. We then report the values $e_{ij}:=1(p_{ij}<.05)$ as indicators of interaction.

This concludes our description of our methodological framework, and explains our automated approach to the selection of “active” interactions that are important to explain the observed spatial pattern.

4 Simulation trials

To check the model fitting procedure before data-analysis we illustrate its characteristics via multiple simulation trials. We will focus on the estimation of the cross-interaction terms, which we can write as a square matrix

[TABLE]

We call $M$ the interaction matrix. Note that the coefficients related to each interaction function $\beta_{ij}^{T}\mathbf{g}_{ij}$ are in groups, meaning that either the whole step-function is estimated to be $\hat{\beta}_{ij}=0$, or one or more of the components $\beta_{ijk}$ is estimated to be non-zero. This is an example of the group lasso (Yuan and Lin, 2006).

The key quality metrics we will look at to assess group level performance are the true positive rate (TP) and the false positive rate (FP), stratified to intra- and inter-type interactions when relevant.

We simulated patterns in either $W=[0,10]^{2}$ or $W=[0,1000]\times[0,500]m$ window with different settings for interactions. Experiments 1,2,4 study the method under the scenario of a correctly specified model, i.e. we simulate and estimate the multi-range multivariate model. In Experiments 3 and 5 we simulate Cox models to see how general and flexible the step-function (shrinkage) approach can be for interaction discovery. Simulations of our model were carried out using the Birth-and-Death algorithm for Figure 7, but for the trials we will use a Metropolis-Hastings algorithm with fixed point counts to keep intensities at desired levels (for more detail we refer to Illian et al., 2008, p. 147-154).

4.1 Experiment 1: Interactions in a small pattern

We first simulated a small ($p=4$) example to familiarise the reader with our methodology and check that the method works as intended. We simulated 100 realisations of the multi-range multivariate saturation model with per type point counts $(100,100,50,150)$ in a $[0,10]^{2}$ window. The saturation levels were all set to $c=1$, producing low levels of interaction. Intra-type ranges were set to $\mathbf{r}_{i}=(.1,.2,.3)$, and inter-type ranges were set to $\mathbf{r}_{l}=(.1,.4)$. Types 1 and 2 were set to exhibit internally a mixture of short range repulsion and medium range clustering ($\beta_{ii}=(-1,1,0),i=1,2$), type 3 had some medium range clustering ($\beta_{33}=(0,1,0)$), and type 4 had no internal correlation. A positive inter-type correlation ($\beta_{ij}=(.6,.3)$) was set between types 1 and 2 and types 3 and 4. The true range vectors were used for fitting. Strongly penalising hyper-prior $Pr(e_{g}=1)\sim Beta(.1,10)$ was chosen for the Bayesian methods. A $4\times 4$ partitioning of $W$ was used for CV to keep the data loss around 50%.

Figure 1 depicts the average detection rates per coefficient $\beta_{ijk}$. Apart from Pearson residuals, the rates of the grouped methods are very similar, indicating the estimated effects do not depend on the chosen algorithm (group lasso or spike and slab). The non-zero structure in types 1 and 2 are detected well by the method. The un-grouped SSGAM and VBSS capture well the medium range clustering for types 1 and 2 but not the short range regularity. This indicates that the medium range clustering is the main signal at the group level. The type 3 sub-patterns had interaction only at medium range, and it seems to be difficult to uncover by any of the methods used. We posit that this difficulty is due to a smaller point count and mixed interaction with type 4. The grouping helps discover the inter-type interactions between the 1-2 and 3-4 pairs, evident from the lower detection rates of the un-grouped methods.

A summary of the true and false positive rates at group level for all methods are given in Table 1. The un-grouped outputs are considered non-zero per group if any of the group members were estimated non-zero. We used a .5 threshold for the Bayesian posterior probabilities for classification. The raw and Pearson residuals clearly are too prone to false positives in this example to be useful.

We also studied how varying the range vectors and dummy intensities affect the outcome. The maximum ranges in the simulations were $R=.3$ for intra- and $R=.4$ for inter-interactions. To simulate miss-specification of the ranges, we altered these by multiplying by a factor .5 or 1.5, and created range vectors for fitting with 1,2,3,4 or 5 equidistant steps. Then we fitted the model using the group lasso with the AIC0.5 rule for penalty with dummy intensity factor $\rho$ either 2, 4 or 8.

Two main features emerged from these studies (cf. Appendix Figure 10). First, choosing ranges that are too short reduces the detection quality quite markedly. Second, the number of grid steps in the interaction functions can be misspecified without large variation in the results, giving a type of numerical robustness. Many small steps in the interactions functions is not recommended, due to the small number of pairs in the data that hit every annuli. For dummy intensity factor $2$ the joint point configuration did not have enough point pairs to fill every bin in the range grid of 5 steps, even after re-sampling the dummies (10 repeated attempts). As we also saw in the results for type $i=3$ above, using many steps may result in lower true positive rates at group level because the group lasso penalises the whole group equally over all its members. The false positive rates are higher for dummy intensity factor $\rho=8$ than for $\rho=4$, a side effect of the AIC0.5 rule-of-thumb. A check of the AIC curves showed that the minimum was often achieved with lower penalty when $\rho=8$, likely due to there being more observations in the logistic regression design matrix, an unsurprising indication that the standard model selection tools are not to be trusted when using pseudo-likelihood. Since we are focusing on the grouped analysis in this study we continue with group Lasso and VBSS.

4.2 Experiment 2: Finding interactions in blocks

Next, we increased the type count to $p=10$, and added interactions in two blocks with in-block pairwise interactions but no interaction between the blocks. The main task was to discover the two blocks with minimal amount of false positives in their cross-section.

The ranges for simulating from the multi-range saturation model were all set to a two step vector $\mathbf{r}=(.25,.50)$, and saturations fixed to $c=1$. The first block of 5 types had short and medium range clustering, $\beta_{ii}=(1,.5)$ for $i=1,...,5$, the second block of 5 types had short range repulsion with mild medium range clustering $\beta_{ii}=(-1,.5)$ for $i=6,...,10$. The inter-type interaction for each pair in both blocks was positive correlation with $\beta_{ij}=(.5,.25)$, with no correlation between blocks. We simulated three intensity scenarios, having point count per type $n_{i}$ either $50,100$ or $200$, so that total point counts per simulated pattern were either $500$, $1000$ or $2000$. The window was again $[0,10]^{2}$.

Fitting was done with the misspecified range vectors $\mathbf{r}=(.15,0.3)$, to increase the challenge (cf. Experiment 1). CV was conducted with a 5$\times$5 partitioning to keep expected data loss around 50%. VBSS was fitted with three hyper-priors $\pi_{l}\sim Beta(.,.)$ ranging from flat (1,1) to medium (.1, 1) and strong (.1, 10) preference for no interactions. The choices correspond to increasing the penalisation in group lasso and facilitate comparisons between the methods.

The rates of interaction detection are shown in Table 2. VBSS produces many false positives, even with the strong prior that should penalise towards sparsity. The raw CV method works well, as does the inverse CV method, but the Pearson CV has a high false positive rate. AIC0.5 over-penalises, thereby missing all of the intra-type interactions, indicating that the rule-of-thumb is not generally useful.

Figure 2 depicts the interaction matrix estimate for a single realisation and the mean interaction over the simulations when $n_{i}=100$, as given by VBSS with the strong prior and inverse CV method. We plot the main diagonal of the interaction matrix from the southwest corner to the northeast corner, situating the $(1,1)$ pixel at the southwest corner, as per Figure 1 in Flügge et al. (2014). Note that white indicates that the detection rate was unity.

We note that the first block is clear in VBSS matrices but the second block, with intra-type regularity mixing with the inter-type clustering, is not so clear and many false positives have been detected. The group lasso is not so efficient in detecting the intra-type interactions but the inter-type interactions are much better detected, with only a few false positives on average.

Further tailoring of ranges will improve quality, as per Experiment 1. Adjusting the priors helps VBSS but provides little practical gain over the automatic penalty rules for the group lasso which work reasonably well on all cases. As the priors lead to varying results, analysis with VBSS needs to include a sensitivity analysis, similar to the lasso penalisation selection. The major issue is that instead of having one tuning parameter there are several, and cross validating the space becomes infeasible as (at least with our implementation) the VBSS algorithm is not much faster for a single prior choice than the group lasso algorithm is for the whole penalty path. The usual Bayesian benefits such as posterior distributions and variances are not worth pursuing either due to the likelihood approximation. For these practical reasons we will continue with only the group lasso to perform computational feasible variable selection.

4.3 Experiment 3: Detecting blocks of interacting Log-Gaussian Cox processes

Next we check if interactions can be detected in correlated log-Gaussian Cox process data. We do not expect a large power here due to the local Poisson distribution of the patterns, but it is important to check how flexible the detection is under model misspecification as the rainforest data potentially exhibits a variety of spatial mechanisms.

We simulated homogeneous, stationary multivariate log-Gaussian Cox processes of $p=24$ types, structured into 3 correlated blocks of 8 types each. Inside a block the 8 latent Gaussian fields are correlated LMC fields (Gelfand et al., 2002), for which we set the cross-field correlation levels to either $.6$ or $.8$. All fields were driven by a Matérn covariance function with smoothness $\nu=10$, marginal variance $\sigma^{2}=2$ and correlation range $3$ (so that $cor(r>3)<.1$). For illustration, four sub-patterns from two different blocks are plotted on the top row of Figure 3, overlaid on their generating fields. We ran the experiment with constant $n_{i}=50$ and $n_{i}=100$ points per type. The window was $[0,10]^{2}$.

Estimation range vectors for intra- and inter-type interactions were set to $\mathbf{r}=(.25,.5)$ after examining one realisation’s cross pair correlation function. The CV’s were carried out using a $3\times 3$ partitioning to keep the expected data loss at 50%.

Table 3 lists the detection rates. The block structure is detected to some extent, with somewhat elevated false positive rates. The raw residual CV performs the best, Pearson CV resulting in high false positive rates, and inverse CV landing overall somewhere in between. The AIC0.5 rule penalises too much to detect more than approximately 25% of the block structure. As expected, an increase in inter-field correlation improves the detection, as does doubling the point count. The bottom row of Figure 3 depicts the estimated interaction matrices for the raw CV when point counts were $n_{i}=100$. Even though the inter-type interactions are not always detected correctly, the intra-type clustering are detected quite well.

We noticed that improvements are easily achieved by adjusting the range vectors used in the fitting; but eye-balling a single realisation’s pair correlation function is obviously not optimal for a repeated experiment. A practical solution is to choose the ranges per realisation based on a proper exploration of the data, in which case we expect quite a good performance considering how misspecified the model is for the generating mechanism. During our development of the methods, further improvements in this experiment were also achieved by increasing the level of clustering in the patterns, either by increasing the variance or by decreasing the range. For example, setting $\sigma^{2}=3$ (the original level was $\sigma^{2}=2$) reduces false positive rates and increased true positive rates up to 10% across the table. The presented results were chosen as a generic illustration, and, by eye at least, the simulation parameters produced clustering similar to that observed in the rainforest data.

4.4 Experiment 4: Rare interactions

To focus more specifically on the false positive rates of the candidate methods, we simulated patterns where, at most, there is only one pairwise inter-type interaction that is non-zero. We considered four scenarios: (1) homogeneous multivariate Poisson patterns; (2) homogeneous multivariate Poisson patterns with an extra intra-correlated pattern; (3) inhomogeneous multivariate Poisson patterns, and (4) inhomogeneous multivariate Poisson patterns with an extra intra-correlated pattern.

For this and the following Experiment 5 we used real covariate maps from the BCI dataset. Thirteen different soil variables (such as soil Ph and magnesium concentration), together with elevation and elevation gradient maps, are available on a $20\times 20m$ grid covering the observation window $W=[0,1000]\times[0,500]m$. For this experiment (and the data analysis later on) we decomposed the 15 covariate maps using a singular value decomposition of the point-wise measurement matrix, and kept the six component maps corresponding to the six largest singular values, capturing approximately 70% of the features. These six pointwise independent PCA maps were then used as covariates.

To generate trends for the inhomogeneous multivariate Poisson processes, we combined the PCA covariates linearly with random coefficients taking values $-1$, 0 or $+1$, with probabilities $.25,.5,$ and $.25$, respectively. For each simulation of the multivariate process, two trend maps were generated and exponentiated to work as intensity surfaces. Each of the two intensity surfaces was assigned to half of the $p$ sub-processes. This way each simulated multivariate pattern has two blocks of types, $p/2$ each, with very strong intra-block correlation that is fully explained by the covariates. For each simulation the point counts per type were set to range from $50$ to $300$, increasing log-linearly so that lower counts were more common. We repeated the experiment for $p=10,20$.

The extra intra-correlated patterns were generated with the multi-type multi-scale model with homogeneous intensity of 100 points, and exhibited very short range repulsion and medium range clustering ($\mathbf{r}=(1,20),\beta=(-10,1),c=(3,3)$). For fitting we set all ranges to $\mathbf{r}=(7,15)m$.

The results for the four scenarios are given in Table 4. With the Poisson data, the methods all produce some false positives. The inverse CV is best, with around 6% false positive rates. With the exception of Pearson CV, including the extra interacting type sharpens the results by reducing the false positive rates. Every method discovers the extra interaction perfectly. Including covariate dependency, and producing correlation, does not increase the false positive rates overall. Again, the inclusion of a single interacting type reduces the false positive rates while itself being clearly detected. Doubling the type count roughly quadruples the parameter count, but it does not seem to affect the quality as false positive rates go down, albeit not always by a factor of 4.

In the estimation procedure we did not penalise the covariate coefficients, because we found that penalising them using the group lasso lead to under-prediction of the trend due to the shrinkage effects of lasso, consequently leading to under-penalisation by the cross validation. When the covariate effects are to be considered more closely we suggest a bias correction step to the lasso or the use of a less strongly penalizing added term such as the SCAD or MCP (Breheny and Huang, 2009).

4.5 Experiment 5: Independent patterns in a rainforest landscape

As a final synthetic example we simulated $p=64=4\times 16$ independent, inhomogeneous patterns in $W=[0,1000]\times[0,500]m$ window with first order interaction depending on the covariates in the BCI data. We selected half ($2\times 16$) of the patterns to come from a Thomas cluster process (Illian et al., 2008, Sec. 6.3.2) with two different dispersal ranges and points-per-cluster rates, to arrive at an appropriate degree of heterogeneity. We set the first block’s model to generate patterns with a few large clusters (”Thomas 1”), and the second block’s model to generate patterns with many small clusters (”Thomas 2”). The second half of the patterns was generated using the multi-type multi-range model with either repulsion followed by clustering (”Geyer 1”) or just repulsion (”Geyer 2”). The ranges of the multi-range models were inversely dependent on target point count: This way small amount of points could spread out more realistically, and if target point count was high then maximal packing density would not be violated (i.e. only a finite amount of ”discs” with certain radius fit inside $W$).

To simulate habitat effects we connected the patterns to the maps of four covariates $Mn,P,pH$ and $grad$ which have been found to be relevant covariates for the rainforest population (Schreeg et al., 2010). The covariate values were first standardized. Then for each simulated species we sampled uniformly $t\in\{0,..,4\}$ of them, and summed them point-wise with weights $(t,...,1)/(t+1)$ to produce a surface. The surface was then used as the unscaled intensity field for Thomas processes generator points, and as first order field for the multi-step multivariate Geyer models. The covariate choices and the coefficients were kept fixed over the replicates of the multivariate simulations.

In a single realisation of the process, each of the 4 models were simulated 16 times, with target intensities within each 16 type block ranging from 50 to 1000 points, median 225. The realised point counts ranged from 34 to 1024 due to edge effects while simulating the Thomas process (points outside $W$ were dropped). The point count in a single multivariate realisation was around 20800.

Figure 4 depicts 8 sub-patterns in a typical realisation, each overlaid on its inhomogeneity generating habitat field. ”Thomas 1” patterns have about half the number of clusters of the ”Thomas 2” patterns. Note that it is hard to see from the plots the very short scale features operating at ranges 1-10m as the area is so large.

For fitting we used ranges $\mathbf{r}=(10,20)$ and 7$\times$4 partitioning to keep the estimated data loss at 50%. Table 5 shows the detection rates for the proposed methods together with the MC based Studentised deviance and rank envelope test, both based on the cross-$K$ functions estimated over ranges .5-15m and 1-50m (corresponding to short and medium-long range testing scenarios) with 999 simulations of the IPP null-model, as described in Section 3. Apart from the Pearson CV, and the AIC0.5 rule-of-thumb which penalises too much, the detection rates are good. Around or below 5% of the false inter-type interactions are detected while 74-96% of the intra-type interactions are detected. The model based and the Monte Carlo based estimates are very close to each other.

Somewhat surprisingly, the Thomas patterns, for which the model is miss-specified, are not the harder of the two families to discover. Figure 5 shows the averages and examples of interaction matrix estimates for the Monte Carlo tests and raw and inverse CV methods. We placed the $(1,1)$ interaction at the southwest corner of each image. The matrices are ordered by model (Thomas 1, Thomas 2, Geyer 1, Geyer 2) and increasing point count from bottom left to top right within the model blocks. For intra-type interactions (on the diagonal), the rarer species’ interaction within the Geyer 2 model seem to be harder to discover with the CV methods. This could be due to the fixed $\mathbf{r}$ that was used for estimation when the processes had varying ranges. The MC null-hypothesis design, where the types corresponding to columns in the interaction matrix were kept fixed and the types corresponding to rows were randomised, leads to non-symmetric estimates. This is an important feature of the MC method, and as we return to in the data analysis below, can lead to interpretation issues. The Geyer 1 block, where there is short scale repulsion, and medium scale attraction, is the major source of false positives for the MC tests, and types with higher abundance show generally higher false positive rates. This is potentially important, because such a mixture of repulsion and attraction is likely to be common in many plant communities.

4.6 Summary of the simulation trials

The results from the simulation studies are very encouraging despite the sometimes large dimensionality of the underlying problems, especially demonstrated with Experiment 5. In Experiment 1 we established feasibility of our model, and the agreement between penalisation methods validate the approach. Experiment 2 demonstrated the model’s ability to discover block community structures in a moderately complicated 10 type pattern. Experiment 3 showed that the model can detect significant interactions in not just realisations of itself, but also in realisations of Cox processes. Experiment 4 showed that false positive rate does not need to be high for us to be able to discover rare events. Finally in Experiment 5 we saw that the method handles rainforest-like patterns with covariate effects, varying point counts and various interaction mechanisms and models well.

The experiments suggest that CV with Pearson residuals is not reliable in practice. The raw residuals work well in some cases (Experiments 2 and 3), but fails completely in others. The inverse residuals produced most consistently good detection. When the raw and inverse residuals both worked well, the variance of the inverse residuals were higher, so computing both is advisable.

As noted in the Experiments, we did not specifically tailor the range parameters for each fitting task per se. Adjusting the model’s range-scales according to specific analysis is highly recommended as they are expected to improve the results in practice.

5 Data example: Rainforest interactions

The BCI rainforest dataset consists of multivariate point patterns corresponding to censuses of rainforest plants (shrubs and trees) living in a $W=[0,1000]\times[0,500]m$ area of Barro Colorado Island, Panama. Censuses have been taken regularly since 1981 (Hubbell et al., 2005; Condit, 1998; Condit et al., 1999). In each census, woody plants (shrubs, trees) with diameter at breast height (dbh) over 1cm were cataloged, noting their dbh, location, species, condition and some other details not relevant to our studies in this article.

The total number of species is about 300, with slight variation over the years due to immigration and extinctions in $W$, which is physically an open area allowing for migration . For the example data analysis we chose the 2005 census, and selected some specific species. Recent studies (Kanagaraj et al., 2011; Yang et al., 2016) suggest that spatial features of the plants vary with life stage. Normally the distribution of young plants is more clustered due to seed dispersal and adult plants’ distribution more regular due to competition based self thinning, and that resource requirements also change with maturation. We therefore selected a subset of species for which an estimate of reproducible size, a surrogate for juvenile/adult life stage, was available (unpublished data by R. Foster, available as supplements for Flügge et al., 2014). We included only adult plants and excluded species with less than 50 adults in the region, leading to a multivariate pattern of $p=83$ species. Point counts vary from 50 to 8784, with a median of 118 and a total of 31650.

5.1 Interaction detection using the introduced methods

We carried out the interaction estimation with a multi-range multivariate saturation model using group lasso penalisation as in Experiment 5. Habitat effects were accounted for by including the 6 PCA maps as covariates as described in Experiment 4 (using PCA maps avoids problems with co-linearity of covariates). We set the range vector equal across all 3486 intra- and inter-type interactions, and fitted the model twice, with range vector $\mathbf{r}=(7,15)m$ and with range vector $\mathbf{r}=(7,15,30)m$, the ranges being concordant with previous results of neighbourhood dependent growth models within the BCI forest (Table 4, Uriarte et al. 2004). We used the dummy intensity factor 4, with a minimum of 500 dummies per type to avoid singularities. The saturation parameters are set according to Section 2.3. To correct for edge effects we implemented a 30m buffer zone, and, to keep expected data loss below 50%, the CV partitioning was fixed to 6$\times$3. For comparison, we also implemented the MC tests as described in Section 3.

5.2 Results

The resulting interaction matrices for the two MC tests are depicted on the first two columns of Figure 6. The group lasso estimates with raw and inverse residual CV penalty selection are depicted on the two rightmost columns. The top row results involve two ranges up to 15m, while the bottom row involves three ranges up to 30m. The MC test interaction matrices are not symmetric by design. For example, according to the rank envelope test at ranges .5-30m, species $i=10$ interacts significantly with 1 species when the other species in a pair is randomized in the test, but it interacts significantly with 22 species when itself is the one randomised in the test. Although, biological interactions (competition) can be asymmetric, with one species being a superior competitor to a second species, the spatial correlations that emerge from these interactions should be symmetric. The non-symmetric spatial interaction matrices produced by the MC method are therefore problematic when it comes to interpretation. In comparison, the inverse CV model-based method estimates the number of interactions for species $i=10$ to be 5, of which 4 are amongst the 22 species indicated by the MC test. Another difference is that the MC test outcome suggests that more abundant species interact more. The model based results do not indicate such a trend (cf. Appendix Figure 11); so when the number of inter-type interactions are of similar order, the distributions are different. Finally, nearly all species are deemed internally structured by MC, but around or less than half by the model-based approach (Table 6).

We have designed a complete analysis pipeline that can be used to answer a number of specific ecological questions, and in particular uncover important processes affecting the forest. For example, we could now proceed by analysing the covariate effects (Appendix Figure 12) to see how important environment (soil, elevation) effects are in structuring the forest. To then consider structure caused by second order effects, we could further study the interaction matrices to discover sub-communities and interacting groups as was done by Flügge et al. (2014) and Morueta-Holme et al. (2016) based on different MC approaches.

6 Discussion

We have developed a flexible model-based method for detecting small scale interactions in multivariate spatial point patterns, demonstrated its potential on several synthetic examples, and applied it to a large rainforest data set. This is in contrast to using log-Gaussian Cox processes for multivariate modeling, as in for example Waagepetersen et al. (2016), where small scale interactions are harder to capture. In so-doing we have greatly extended the potential applicability of model-based inference for multitype point patterns, allowing both for the analysis of more interactions, and more complex multiscale behaviour.

In the Thomas cluster process, and in Cox processes in general, the two-stage generation of ”offspring” given ”mothers” is a natural model for plant regeneration in natural plant populations. However, as the offspring are assumed mutually independent, the model is suitable for only a particular type of pattern, such as just germinated individuals. The individual-to-individual competition and survival over a plant’s lifetime prior adulthood is better captured with explicit point-wise interaction models, such as the Gibbs models. Alternative models that emulate natural thinning in addition to natural clustering can be constructed (Stoyan, 1979; Andersen and Hahn, 2016; Lavancier and Møller, 2016), but fitting such models rely heavily on non-parametric methods for which simultaneous multivariate analysis is currently not well understood. Our approach takes a more data-analytical approach by including all effects simultaneously and addressing the high dimensionality with penalised inference.

The model can be tailored for specific applications with relatively low effort since apart from the point-to-point indicator functions the analysis pipeline is fixed. For example, adding Strauss components might detect clustering more efficiently as the saturation level of the Geyer components caps the clustering tendency. Such a model would not work as a generative model and simulations would not be stable (singular clusters would form, Gates and Westcott 1986), but it could be effective as a statistical interaction discovery tool.

The pipeline is also computationally practical. Group lasso algorithms are very fast and work well with computational underlying calculations for sparse matrices, which the approximation design matrices for Gibbs models are by construction. On a regular laptop it took us days to compute all MC estimates to compare to our model, but it took only hours to estimate the full model. Of course, one can parallelise the MC estimation since the type-to-type interaction tests are performed independently.

An advantage of our model-based approach is the explicit treatment of covariates that may affect the intensity but not the interactions directly. This means it is possible to infer which covariates are important for the distributions of the different species. In contrast the MC approach that predominates the ecological literature relies upon inhomogeneous Poisson processes to capture the first order (environmental) processes, and the user needs to determine/estimate the scale over which these processes are affecting the intensities. The current trend has been to use one scale for all species, but this is unlikely to be optimal. However, we note that such covariates are not always available, and the MC method might be preferable when covariates have not been measured but are thought to be important. The within-type spatial structure that is commonly found in the data is also broken down when conducting the null model test in the MC approach, and it is not clear what effect this might have on results. Nonetheless the user also has to make important decisions in our model-based approach, such as the range of scales over which to carry out the analysis, and as we discuss below the model validation and penalisation methods also are an area of future exploration.

The computational burden of model fitting was increased tenfold because we determined the best penalty level using cross validation. To simplify the penalty selection, we tried including a few extra Poisson ”noise” sub-patterns into the data to select the penalisation so that the noise stayed independent. The result of this approach was inconclusive, and remains an area of future investigation. We also tried using cross validation with the constructed design matrix of the logistic regression likelihood, as is common for the the lasso applications, but much like with the AIC, too many false positives resulted for this approach. Standard model selection approaches do not do well because the likelihood is an approximation, rather than the correct form of the likelihood. Correcting the score and Hessian of the pseudo-likelihood while doing penalisation would be useful in this regard, and some work has already been done in the unpenalised case by Coeurjolly and Rubak (2013). Using a corrected Hessian would also allow us to do better inference on the covariates, as (approximate) confidence intervals could then be constructed.

Replacing the group lasso with some other penalisation could have several benefits. Some penalisations shrink non-zero coefficients less which is ideal both for the predictions needed in the cross validation step, and downstream analysis of the estimated effects. For example, we suspect the shrinkage issue prevented the method from working properly when covariates were penalised and consequently the inhomogeneous trends became too flat. A further limitation of the group lasso is that all group members were penalised equally, this in turn leading to low group detection with many range-steps and low amount of data per range annuli. A more refined group penalisation method such as sparse-group lasso or minimax concave penalty is more sensitive to individual group members being non-zero (Breheny and Huang, 2009) and would allow for more detailed interaction functions with more steps. Adding steps does however require using more dummy points so that each annuli registers something and numerical problems are avoided. We could also replace the step functions with overlapping components, such as radial basis functions, and gain not only numerical stability but also smoother interaction function estimates. Further refinement to the CV optimisation penalty grid is also needed in actual applications. In our high-$p$ examples the 100 step grid was often too coarse to include/exclude individual groups.

The model is constructed via small point-to-point interactions, and any large scale unobserved variability is not equally well captured. It might be possible to add spatial random effects, such as Gaussian processes, to the first-order interaction and still use the same pipeline. This would bridge the interface between Gibbs models and log-Gaussian Cox processes, an exciting area of future investigation. This paper therefore stands as a further step to understand general heterogeneous and high dimensional point process observations.

We would also like to point out that the modular construction of the point-to-point interaction functions provides a potentially useful connection to dynamic modeling of ecological communities defined through generations of dispersal related birth and competition based death events (Law and Dieckmann, 2000; Law et al., 2003). If the transition probabilities are modelled log-linearly, the stationary distribution (if it exists) of the corresponding birth-and-death process following the detailed balance condition would be a Gibbs process of the kind we have discussed here in the more general, potentially inhomogeneous and highly multivariate setting.

7 Acknowledgements

The BCI forest dynamics research project was founded by S.P. Hubbell and R.B. Foster and is now managed by R. Condit, S. Lao, and R. Perez under the Center for Tropical Forest Science and the Smithsonian Tropical Research in Panama. Numerous organizations have provided funding, principally the U.S. National Science Foundation, and hundreds of field workers have contributed. We thank the EPSRC for support via EP/N007336/1 and EP/L001519/1, and S. C. Olhede also acknowledges support from the 7th European Community Framework Programme (a European Research Council Fellowship via Grant CoG 2015- 682172NETS).

We thank the two anonymous reviewers for their very helpful comments.

8 Software and data

All computations were implemented in R (v.3.3.1). The random fields for the log-Gaussian Cox models were simulated using the R-package ”RandomFields” (v.3.1.3), and uniform simulation and some utilities were used from the R-package ”spatstat” (v1.46.1). For the group lasso we adapted the local coordinate descent algorithm in the R-package ”grpreg” (v.3.0-2) with the inclusion of the offset terms. The VBSS was implemented by hand. For the SSGAM we used the R-package ”spikeSlabGAM” (v.1.1-11), with hyper-prior Beta-distribution settings as for the VBSS and otherwise using default parameters. An R-package implementing the method pipeline is available from the first author.

The BCI data is available for research purposes from

http://ctfs.si.edu/webatlas/datasets/bci/.

Appendix A Details about the likelihood approximation

Let $X$ be the point process generating the pattern $\mathbf{x}$ and let $D$ be a marked dummy point process in the window $W$ which is independent of $X$ and has known intensity functions $\rho_{i},i=1,...p$. Baddeley et al. (2014) proposed to solve likelihood optimization problem by using the estimating function

[TABLE]

where $\rho(u)=\rho_{i}(x)$ for $u=(x,i)$. With the help of the Campbell formula and the Georgii-Nguyen-Zessin theorem (formulas 1.5.10 and 6.6.2 Illian et al., 2008), it can be shown that $s_{W}$ is an unbiased estimation function and that finding the root of $s_{W}$ gives an unbiased estimate of the maximum of the likelihood (5).

The $s_{W}$ as a function of $\theta$ is the derivative of

[TABLE]

which formally is the likelihood of a logistic regression with variables $\tau(u)=1(u\in X)$ for $u\in X\cup D$ and

[TABLE]

In practice the method works as follows: We sample a set of dummy points $\psi=\{(x,t)\}=\cup\mathbf{\psi}_{i}$ with sub-patterns $\mathbf{\psi}_{i}$ having a known distribution with constant intensity $\rho_{i}$ in $W$, for each type $i=1,...,p$. We then calculate the vectors $\mathbf{b}(u):=\mathbf{v}(u;\mathbf{x})$ for each $u\in\mathbf{x}\cup\psi$. Let $N=\#(\mathbf{x}\cup\psi)$. The log-likelihood of the logistic regression can then be written in a compact vector form corresponding to

[TABLE]

where $t(u)=1(u\in\mathbf{x})$ indicate the data points, $B$ is a row matrix of $\mathbf{b}(u)$’s for each $u\in\mathbf{x}\cup\psi$, and $o(u)=-\log\rho_{i}(u)$ are offset terms.

Appendix B Derivation of saturation approximation of the Papangelou terms under independence

We show how the probability of saturation of the components in the log-Papangelou conditional intensity can be derived using independence assumption. First we note that for a homogeneous Poisson process neighbour counts in any set depend only the volume of the set. Thus to simplify from the annuli notation it suffices to consider only $r=r_{k}>r_{k-1}=0$. Then

[TABLE]

where $1()$ denotes the indicator function. Assume now that $\mathbf{x}$ comes from a Poisson point process with intensity $\lambda$. Write $a=\lambda|b(o,r)|$. The Georgii-Nguyen-Zessin formula (Illian et al., 2008, p.399) gives

[TABLE]

Now, the random variable $y:=\#[b(o,r)\cap X]$ is $Poisson(a)$ distributed. Denote its CDF with $F_{a}$. Then

[TABLE]

giving the function $t(c)$ in the text. Note that we also use the fact $\lambda\approx n/|W|$.

Appendix C Illustration of the interaction and potential functions

Figure 7 (top row) depicts two particular potential shapes with different $\mathbf{r}$’s and $\beta$’s, to illustrate their role. The first function corresponds to a decreasing attraction in range, and the second function has both a repulsion and an attraction component. How the attraction and repulsion impact the point pattern depends on the choice of the individual $g_{ijk}$’s: In this example we use the saturation model $\mathbf{g}_{ij}$’s with $c_{ijk}\equiv 1$. The first two subfigures in the middle row of Figure 7 show univariate simulations from each of the two interaction functions. The third subplot in the middle row shows a bivariate simulation with intra-type interaction given by the first potential and inter-type interaction given by the second potential. The bottom row shows the conditional intensity (12) at the window locations $u\in W$ given the simulated patterns.

In the case of the bivariate pattern the conditional intensity is for points of the first type. The conditional intensities exhibit various features: For the first potential function, high potential locations are near data-points, but only if the data-points do not already have neighbours (those data-points’ neighbourhoods are already saturated). With the second potential any location too near the data-points has a low potential (repulsion), but being too far from the data-points is not encouraged by the potential either (attraction). In subfigure 3c) we see the complex mixture of first type’s internal potential and the intra-type potential: A point of the first type would be welcome near a point of it’s own kind (attraction) but unwelcome near a point of the second kind (repulsion).

Appendix D Additional figures