Fast, Approximate Maximum Likelihood Estimation of Log-Gaussian Cox Processes

Abstract

Log-Gaussian Cox processes (LGCPs) offer a framework for regression-style modeling of point patterns that can accommodate spatial latent effects. These latent effects can be used to account for missing predictors or other sources of clustering that could not be explained by a Poisson process. Fitting LGCP models can be challenging because the marginal likelihood does not have a closed form and it involves a high dimensional integral to account for the latent Gaussian field. We propose a novel methodology for fitting LGCP models that addresses these challenges using a combination of variational approximation and reduced rank interpolation. Additionally, we implement automatic differentiation to obtain exact gradient information, for computationally efficient optimization and to consider integral approximation using the Laplace method. We demonstrate the method’s performance through both simulations and a real data application, with promising results in terms of computational speed and accuracy compared to that of existing approaches. Supplementary material for this article is available online.

KEYWORD:

• Automatic differentiation
• Fixed rank kriging
• Integrated nested Laplace approximation
• Point patterns
• Spatial point process
• Variational approximation

1 Introduction

Researchers often wish to model a spatial point pattern—a collection of event locations in some domain of interest—as a function of a set of predictors measured over the spatial domain. A motivating example in this article is a survey of gorilla nesting sites in Kagwene Gorilla Sanctuary, Cameroon (Funwi-Gabga and Mateu 2012, Figure 1(A); further detail can be found in Section 5). We would like to model this point pattern as a function of elevation and other spatial variables in order to better understand where gorillas nest and why. This can be done using a spatial point process model. However, point patterns often exhibit clustering that cannot be explained by measured predictors alone, rendering the widely-used Inhomogeneous Poisson Process (IPP) model inappropriate, as we see for the gorilla data (Figure 1(B), top panel).

Fig. 1 A: An observed point pattern of gorilla nesting locations in Kagwene Gorilla Sanctuary in Cameroon. B: An inhomogeneous K function for the point pattern (Baddeley and Turner 2000), with simulation envelope (computed as in Myllymäki et al. 2017), to check goodness of fit under an IPP model (top panel) and LGCP model (bottom panel). For model details see Section 5. Note there is evidence of additional spatial clustering to that accounted for by the IPP model, but not so for the LGCP.

The log-Gaussian Cox Process (LGCP; Møller, Syversveen, and Waagepetersen 1998) offers a way to incorporate additional spatial clustering into point process models that is not captured by measured covariates. This is achieved by including a Gaussian Random Field (GRF) to induce additional spatial correlation between observations—effectively acting as a spatially correlated error term in the model. LGCP models are particularly appropriate in instances where clustering arises from missing or unmeasured environmental processes/phenomena, as opposed to those in which clustering/dispersal is due to interactions between the point events themselves (Diggle et al. 2013). Here we are focusing on a LGCP framework in the context of modeling a realized point process (via its point pattern) driven by a collection of measured and unmeasured predictors throughout some spatial domain, rather than trying to model, say, potential parent and offspring processes that sometimes underlie interactions between points (and see also Adams, Murray, and MacKay 2009, Gonçalves and Gamerman 2018, who fit different types of Cox process outside of the log-Gaussian framework).

Fitting LGCP models is computationally challenging. The largest burden in fitting LGCPs arises from the variance-covariance structure of a finite-dimensional realization of the latent Gaussian process. Taylor et al. (2013) provide an R package that fits a LGCP via Markov-chain Monte Carlo (MCMC) methods, initially using Metropolis-Hastings algorithms via Gibbs samplers, although a Metropolis-adjusted Langevin algorithm (MALA) is now more common (Shirota and Gelfand 2017). These algorithms can be very costly computationally, and do not scale well as the size of point patterns increases, because of the need to sample the high-dimensional Gaussian random field at each iteration. Faster algorithms have been proposed using a few general strategies: a reduced rank approximation to the latent variance-covariance matrix (Chakraborty et al. 2011, Simpson et al. 2016); inducing sparsity using, for example, nearest-neighbor Gaussian processes (Datta et al. 2016, Shirota and Banerjee 2019); aggregating point events to counts (Johnson, Diggle, and Giorgi 2019); or approximating the posterior distribution using variational Bayes mean-field (Teng, Nathoo, and Johnson 2017) or innovative applications of the Laplace method (Rue, Martino, and Chopin 2009). Of these options, Integrated Nested Laplace Approximation (INLA; Rue, Martino, and Chopin 2009, Simpson et al. 2016) is currently the most widely used in practice, exploiting work by Lindgren, Rue, and Lindström (2011) to approximate a Gaussian-Markov Random Field (GMRF) using a reduced rank approach that is motivated as a solution to a stochastic partial differential equation. In simulation studies, Taylor and Diggle (2014) found INLA to be generally less accurate but faster than MCMC-MALA assuming covariance parameters are known a priori, and Teng, Nathoo, and Johnson (2017) found INLA to have similar accuracy but substantially faster computation times than variational Bayes mean-field and Hamiltonian Monte Carlo approaches. Bachl et al. (2019) wrote a software wrapper (inlabru) to facilitate use of INLA to fit LGCP models by non-specialists.

In the machine learning literature, sparse Gaussian processes are a popular approach (Quinonero-Candela and Rasmussen 2005) closely related to other reduced-rank Gaussian processes used to induce correlation in hierarchical models (Banerjee et al. 2008, Datta et al. 2016). Variational Bayesian approximations have also been explored in this context (Titsias 2009, Hensman, Matthews, and Ghahramani 2015), including for point patterns (Lloyd et al. 2015). In this article we also propose a variational approach and exploit sparsity to reduce computation times, but we do so in the context of log-Gaussian Cox Processes, fitted by (approximate) maximum likelihood.

Within the frequentist paradigm, Waagepetersen and Guan (2009) proposed a two-step fitting procedure for this problem, which involves estimating first-order parameters of the mean intensity function using a Poisson likelihood score function, then estimating higher-order parameters of the latent field to minimize contrasts between the theoretical and parameterized K functions. This idea is related to method-of-moments estimation, and while it can achieve fits in seconds for quite complex models, it does so by emphasizing pragmatism rather than statistical efficiency.

In this article we propose a fast, novel maximum-likelihood approach to fitting LGCP to spatial point patterns, involving three innovations outlined in Section 3. First, variational approximation (VA) permits a closed form approximation to the marginalized log-likelihood. Second, we use a reduced-rank approximation to the large spatial variance-covariance matrices that arise. Finally, automatic differentiation is used to quickly obtain gradient information for efficient optimization and inference. Performance in fitting LGCP is also examined in comparison to INLA in simulation (Section 4) and application (Section 5). We conclude the article with a brief discussion (Section 6).

2 Log-Gaussian Cox Processes

We denote the domain of interest by $\mathcal{D}$ which is continually indexed by spatial parameter s so that we write the presence-only data as $S_n=\left\{{\mathit{s}}_i\right\}$ for $i=1,\dots ,n$, that is, a collection of locations at which a species of interest is recorded. In this article we consider $\mathcal{D}\subseteq {\mathbb{R}}^2$, however, the methods could be applied more generally.

The target of inference for point process models is their spatially varying intensity function $\lambda (\mathit{s})$, which describes the limiting number of point events per unit area in some infinitesimally small region around the point s. We assume the log-intensity is a linear combination of predictor variables, $\mathit{X}(\mathit{s})$ with coefficients $\mathit{\beta }$. For conciseness, any intercept term is absorbed into $\mathit{\beta }$, with a corresponding indicator variable in $\mathit{X}(\mathit{s})$. For LGCPs, the intensity inherits stochasticity from a zero-mean, latent Gaussian process, $\xi (\mathit{s})$, characterized by some valid covariance function ($f_{\xi }$) governed by parameters ${\mathit{\theta }}_{\xi }$. Note this model assumes that there are two stochastic outcomes to the observed point pattern, $S_{N=n}$, that is, the realization of the Gaussian field and the realization of the Poisson process. The LGCP intensity is given by (1) $$\hspace{0.17em}\ln \hspace{0.17em}\lambda \left(\mathit{s}\right)=\mathit{X}\left(\mathit{s}\right)\mathit{\beta }+\xi \left(\mathit{s}\right).$$(1)

We further assume that conditional on ξ, S_n is the realization of an IPP with intensity $\lambda (\mathit{s})$.

We will use $\pi (·)$ throughout to denote the probability density function of a random variable (or stochastic process) in general, with different probability density functions distinguished by the arguments of $\pi (·)$. It can be shown that the point pattern has (conditional) probability density: (2) $$\pi \left(S_n|\xi \right)=\left\{\prod _{i=1}^n\hspace{0.17em}\exp \hspace{0.17em}\{\mathit{X}\left({\mathit{s}}_i\right)\mathit{\beta }+\xi \left({\mathit{s}}_i\right)\}\mathrm{}\right\}\hspace{0.17em}\exp \hspace{0.17em}\left\{\mathrm{}\left|\mathcal{D}\right|-{\int }_{\mathcal{D}}\lambda \left(\mathit{t}\right)\mathrm{d}\mathit{t}\mathrm{}\right\}$$(2) with respect to a unit rate Poisson process, where $\left|\mathcal{D}\right|$ denotes the size of the domain of interest, see Daley and Vere-Jones (2007) for details. The proportionality constant $\hspace{0.17em}\exp \hspace{0.17em}\left\{\right|\mathcal{D}\left|\right\}$ is dropped in subsequent use for convenience.

The spatial integral that describes the mean is often approximated using numerical quadrature (e.g., Berman and Turner 1992, Illian, Sørbye, and Rue 2012): (3) $${\int }_D\lambda \left(t\right)\mathrm{d}t\approx \sum _{i=n+1}^{n+q}{w}_i\hspace{0.17em}\exp \hspace{0.17em}\left\{\mathit{X}\left({\mathit{s}}_i\right)\mathit{\beta }+\xi \left({\mathit{s}}_i\right)\right\}$$(3) where $\mathcal{D}$ has effectively been discretized into q regions represented by quadrature points ${\left\{{\mathit{s}}_i\right\}}_{i=n+1}^{n+q}$, in addition to the n point events, so that each $i=n+1,\dots ,n+q$ occupies a neighborhood of area $w_i\equiv \left|{\mathit{s}}_i\right|$ such that $\left|\mathcal{D}\right|={\sum }_{i=n+1}^{n+q}{w}_i$. The quadrat sizes, $\left|{\mathit{s}}_i\right|,i=n+1,\dots ,n+q$, need not be constant (i.e., forming a regular grid) and could, for example, be estimated using a Dirichlet tessellation (Simpson et al. 2016). In this article, for simplicity, we will choose quadrature points on a regular grid such that they have equal quadrat size $\left|{\mathit{s}}_i\right|=\frac{\left|\mathcal{D}\right|}{q}$ for $i=n+1,\dots ,n+q$, and we set $\left|{\mathit{s}}_i\right|=0$ for presence points, $i=1,\dots ,n$. Quadrature weights set in this way are not a function of the observed point pattern S_n, which is convenient for simulation, or analysis of multiple point patterns in a common domain.

Applying the approximation in (3) involves an n + q dimensional realization to the Gaussian random field $\xi (\mathbf{s})$, at both the observed point events S_n and the quadrature points, which we write as $\mathit{\xi }={\left[\xi ({\mathit{s}}_1),\dots ,\xi ({\mathit{s}}_{n+q})\right]}^{\mathrm{T}}$. This quickly becomes problematic when computing the marginal log-likelihood of the point pattern, given by (4) $$\mathcal{l}\left(\mathit{\beta }\right)\approx \hspace{0.17em}\log \hspace{0.17em}\int \mathrm{ }\sim \pi \left(S_n|\mathit{\xi }\right)\pi \left(\mathit{\xi }\right)\mathrm{d}\mathit{\xi }$$(4) where $\pi (\mathit{\xi })\sim \mathcal{N}(\mathit{0},\mathit{C})$; $\sim \pi (S_n|\mathit{\xi })$ is (2) approximated with (3); and C is the variance-covariance matrix induced by the covariance function, $f_{\mathit{\xi }}$. Note that C has dimension $(n+q)\times (n+q)$ which becomes a computational bottleneck when either n or q is large. The number of quadrature points q is necessarily large in settings where the predictors, hence, estimated intensity, are not spatially smooth. A further difficulty is that (4) does not have a closed form.

3 Fast Approximate Estimation of Log-Gaussian Cox Processes

Our proposed strategy is to maximize an approximation to the marginal likelihood that addresses key computational bottlenecks in a pragmatic fashion.

3.1 Approximate Marginalization

We consider addressing the intractability of (4) using two types of approximation—a Gaussian variational approximation (Ormerod and Wand 2010, VA), or a Laplace approximation.

Our variational approximation involves re-expressing the marginal likelihood (4) as (5) $$\begin{array}{c}\mathcal{l}\left(\mathit{\beta }\right)=\hspace{0.17em}\ln \hspace{0.17em}\pi \left(S_n\right)=\hspace{0.17em}\ln \hspace{0.17em}\pi \left(S_n\right)·\int \pi \left(\mathit{\xi }|S_n\right)\mathrm{d}\mathit{\xi }\\ =\int \hspace{0.17em}\ln \hspace{0.17em}\left[\frac{\pi \left(S_n|\mathit{\xi }\right)\pi \left(\mathit{\xi }\right)}{\pi \left(\mathit{\xi }|S_n\right)}\right]\pi \left(\mathit{\xi }|S_n\right)\mathrm{d}\mathit{\xi }\end{array}$$(5) then replacing the conditional or “posterior” density of the latent field, $\pi (\mathit{\xi }|S_n)$ with some multivariate Gaussian density function, ${\pi }_{\text{VA}}(\mathit{\xi })\sim \mathcal{N}({\mathit{m}}_{\text{VA}},{\mathit{C}}_{\text{VA}})$. This leads to a closed form approximation to (4), ${\underset{¯}{\mathcal{l}}}_{\text{VA}}(\mathit{\beta })$, which can be written as (6) $$\begin{array}{c}{\underset{¯}{\mathcal{l}}}_{\text{VA}}\left(\mathit{\beta }\right)=\int \hspace{0.17em}\ln \hspace{0.17em}\left(\frac{\pi \left(S_n|\mathit{\xi }\right)\pi \left(\mathit{\xi }\right)}{{\pi }_{\text{VA}}\left(\mathit{\xi }\right)}\right){\pi }_{\text{VA}}\left(\mathit{\xi }\right)\mathrm{d}\mathit{\xi }\\ \approx {\mathbb{E}}_{\text{VA}}\left[\sum _{i=1}^n\mathit{X}\left({\mathit{s}}_i\right)\mathit{\beta }+\xi \left({\mathit{s}}_i\right)\right]\\ \hspace{1em}-{\mathbb{E}}_{\text{VA}}\left[\sum _{i=n+1}^{n+q}{w}_i\hspace{0.17em}\exp \hspace{0.17em}\left\{\mathit{X}\left({\mathit{s}}_i\right)\mathit{\beta }+\xi \left({\mathit{s}}_i\right)\right\}\right]\\ \hspace{1em}+{\mathbb{E}}_{\text{VA}}\left[\ln \hspace{0.17em}\left(\frac{\pi \left(\mathit{\xi }\right)}{{\pi }_{\text{VA}}\left(\mathit{\xi }\right)}\right)\right]\\ =\sum _{i=1}^n\mathit{X}\left({\mathit{s}}_i\right)\mathit{\beta }+{\left[{\mathit{m}}_{\text{VA}}\right]}_i-\sum _{i=n+1}^{n+q}{w}_i\\ \hspace{1em}\hspace{0.17em}\exp \hspace{0.17em}\left\{\mathit{X}\left({\mathit{s}}_i\right)\mathit{\beta }+{\left[{\mathit{m}}_{\text{VA}}\right]}_i+\frac{1}{2}{\left[{\mathit{C}}_{\text{VA}}\right]}_{i,i}\right\}\\ \hspace{1em}+\frac{1}{2}[\ln \hspace{0.17em}\left|{\mathit{C}}^{-1}{\mathit{C}}_{\text{VA}}\right|-\text{Tr}\left({\mathit{C}}^{-1}{\mathit{C}}_{\text{VA}}\right)\\ \hspace{1em}-{\mathit{m}}_{\text{VA}}^{\mathrm{T}}{\mathit{C}}^{-1}{\mathit{m}}_{\text{VA}}+\left(n+q\right)]\end{array}$$(6) where ${\mathbb{E}}_{\text{VA}}$ denotes expectation with respect to the variational density; ${\left[{\mathit{m}}_{\text{VA}}\right]}_i$ denotes the ith element of vector ${\mathit{m}}_{\text{VA}}$; and ${\left[{\mathit{C}}_{\text{VA}}\right]}_{i,i}$ denotes the ith diagonal element of matrix ${\mathit{C}}_{\text{VA}}$. We estimate parameters to satisfy $$\left\{\widehat{\mathit{\beta }},{\widehat{\mathit{m}}}_{\text{VA}},{\widehat{\mathit{C}}}_{\text{VA}},\widehat{\mathit{C}}\right\}=\underset{\left\{\mathit{\beta },{\mathit{m}}_{\text{VA}},{\mathit{C}}_{\text{VA}},\mathit{C}\right\}}{\arg \max }{\underset{¯}{\mathcal{l}}}_{\text{VA}}.$$

That is, we find estimates that simultaneously maximize ${\underset{¯}{\mathcal{l}}}_{\text{VA}}$ for the fixed effects, $\mathit{\beta }$, as well as the variational parameters, ${\mathit{m}}_{\text{VA}}$ and ${\mathit{C}}_{\text{VA}}$, and the prior variance-covariance matrix C. Maximizing with respect to the variational parameters can be motivated as minimizing the Kullback-Leibler distance between ${\pi }_{\text{VA}}(\mathit{\xi })$ and $\pi (\mathit{\xi }|S_n)$ (Ormerod and Wand 2010). A more detailed derivation of these results is available in the supplementary material: Section S1.

An alternative to variational approximation is to use a Laplace approximation for the marginal likelihood. If we again re-express the marginal likelihood in (4) as $$\begin{array}{cc}\mathcal{l}\left(\mathit{\beta }\right)& =\hspace{0.17em}\log \hspace{0.17em}\left[\text{const}.·\int \hspace{0.17em}\exp \hspace{0.17em}\left\{g\left(\mathit{\xi }\right)\right\}\mathrm{d}\mathit{\xi }\right],\end{array}$$then its Laplace approximation is given by (7) $${\underset{¯}{\mathcal{l}}}_{\text{Laplace}}\left(\mathit{\beta }\right)\approx \hspace{0.17em}\log \hspace{0.17em}\left[\text{const}.·\frac{{\left(2\pi \right)}^{\frac{k}{2}}}{{\left|-H_g\left({\mathit{\xi }}_0\right)\right|}^{\frac{1}{2}}}\exp \hspace{0.17em}\left\{g\left({\mathit{\xi }}_0\right)\right\}\right]$$(7) where H_g is the Hessian matrix of $g(\mathit{\xi })=\hspace{0.17em}\ln \hspace{0.17em}\pi (S_n|\mathit{\xi })+\hspace{0.17em}\ln \hspace{0.17em}\pi (\mathit{\xi })$, twice differentiable with respect to $\mathit{\xi }$. The approximation is centered at ${\mathit{\xi }}_0$ which maximizes g and so in this case is the mode of the joint distribution of the latent field and point pattern, $\pi (S_n,\mathit{\xi })$. As the approximation involves only derivatives with respect to $\mathit{\xi }$, we can implement this approach easily using automatic differentiation—see Section 3.3.

Note that the Laplace approximation assumes that the integrand is Gaussian and centeres the approximation about the mode of the joint probability distribution, while a Gaussian variational approximation assumes the posterior distribution of the latent field (given the point pattern) is Gaussian, with the approximation centred about the mean. The idea of assuming the posterior is Gaussian is particularly appealing since its prior distribution is Gaussian, so this remains a plausible approximation in data-poor settings, which need not be true of the Gaussian approximation to the integrand.

We also note that Laplace approximations can perform poorly for high dimensional integrals (Shun and McCullagh 1995), and the integral in (4) would often be high-dimensional. However, in the next section we propose a rank reduction approach where we effectively reduce the dimension to $k\ll n+q$, which addresses this issue.

3.2 Rank Reduction

As mentioned previously, a key computational bottleneck for both the variational and Laplace approximation is that $\mathit{\xi }$ has dimension n + q, and n or q is often large. We address this problem using a reduced rank approximation. That is, we assume (8) $$\xi \left(\mathit{s}\right)\approx \mathit{Z}\left(\mathit{s}\right)\mathit{u}$$(8) where $\mathit{Z}(\mathit{s})=\left[Z_1(\mathit{s}),\dots ,Z_k(\mathit{s})\right]$ is a set of k basis functions and $\mathit{u}={\left[u_1,\dots ,u_k\right]}^{\mathrm{T}}\stackrel{\mathit{iid}}{\sim }\mathcal{N}(\mathbf{0},{\sigma }_{\text{prior}}^2\mathbf{I})$ are random coefficients.

This approach is commonly referred to as fixed rank kriging (FRK) after Cressie and Johannesson (2008). The software package FRK (Zammit-Mangion and Cressie 2021) can assist in implementation, providing options for generating basis functions under a range of different configurations, including approximations for spatio-temporal GPs by indexing the random coefficients over time. Closely related methods are predictive processes (Banerjee et al. 2008), which construct basis functions around knots using kriging, and the GMRF mesh of Lindgren, Rue, and Lindström (2011), which represents a Gaussian process as a collection of piecewise linear basis functions. A key difference here is that the basis functions in (8) are not chosen as a function of covariance parameters, so they do not need to be updated in optimization, which has advantages in computational efficiency and stability.

Choice of number and type of basis functions $Z_r(\mathit{s})$ is a key step, these govern the spatial dependence induced in $\mathit{\xi }$. For $r=1,\dots ,k$, we define a regular grid of “knots” (${\mathit{s}}_r^{*}\in \mathcal{D}$) and compute local bi-square functions of the form (9) $$Z_r\left(\mathbf{s}\right)=\{\mathrm{}\begin{array}{cc}{\left[1-{\left(\frac{d\left(\mathit{s},{\mathit{s}}_r^{*}\right)}{\phi }\right)}^2\right]}^2\hfill & \left|d\left(\mathit{s},{\mathit{s}}_r^{*}\right)\right|\le \phi \hfill \\ 0\hfill & \left|d\left(\mathit{s},{\mathit{s}}_r^{*}\right)\right|>\phi \hfill \end{array}$$(9) where $\phi$ is the function radius, which enforces sparsity. We have found that use of sparse basis functions results in considerable computational savings. If the knots have been organized in a regular grid with spacing e then we choose $\phi =1.5e$ as in Zammit-Mangion and Cressie (2021). By choosing $\phi$ as a function of knot configuration, the only remaining tuning parameter in our basis function configuration is the number of basis functions k. An alternative approach could be to consider the knots as variational parameters as in Titsias (2009).

Alternative basis function configurations are obviously possible. In our explorations we have found negligible effect of the functional form of the basis functions, but the number of basis functions and their radius is important. We suggest these choices be informed by the data, as detailed later. One alternative configuration worthy of noting is to specify basis functions at multiple resolutions (usually two or three), with k_l basis functions of radius ${\phi }_l$ at level l. When k_l varies across resolutions this configuration can capture variation at a range of scales of dependence (Cressie and Johannesson 2008).

Applying the assumption in (8), together with quadrature approximation as in (3), our model simplifies to a type of Generalized Linear Mixed Model (GLMM) with weighted observations. Specifically, the marginal likelihood we wish to maximize has the form: (10) $$\begin{array}{c}{\mathcal{l}}_{\mathit{u}}\left(\mathit{\beta }\right)\approx \int \left\{\sum _{i=1}^n\hspace{0.17em}\ln \hspace{0.17em}{\lambda }_{\mathit{u}}\left({\mathit{s}}_i\right)-\sum _{i=n+1}^{n+q}{w}_i{\lambda }_{\mathit{u}}\left({\mathit{s}}_i\right)\right\}\pi \left(\mathit{u}|S_n\right)\mathrm{d}\mathit{u}\\ +\int \hspace{0.17em}\ln \hspace{0.17em}\left[\frac{\pi \left(\mathit{u}\right)}{\pi \left(\mathit{u}|S_n\right)}\right]\pi \left(\mathit{u}|S_n\right)\mathrm{d}\mathit{u}\end{array}$$(10) where $\hspace{0.17em}\ln \hspace{0.17em}{\lambda }_{\mathit{u}}(\mathit{s})=\mathit{X}(\mathit{s})\mathit{\beta }+\mathit{Z}(\mathit{s})\mathit{u}$. Our reduced rank assumption (8) makes our likelihood approximations much more tractable, as we now marginalize over k-dimensional u rather than n + q-dimensional $\mathit{\xi }$.

For our variational approximation to (10), we replace $\pi (\mathit{u}\mid S_n)$ with $\mathcal{N}(\mathit{\mu },\Sigma )$. These type of variational approximations to GLMMs have previously been studied by Ormerod and Wand (2012) and, given that our model takes such a form, we have the following result:

Theorem 1. A

mong all variational approximations to (10) taking the form $\mathit{u}\mid S_n\stackrel{\text{VA}}{\sim }\mathcal{N}(\mathit{\mu },\Sigma )$, the optimal lower bound is achieved when $\Sigma =\text{diag}\left\{{\sigma }_1^2,\dots ,{\sigma }_k^2\right\}$.

This result follows directly from Theorem 1 of Ormerod and Wand (2012), since we have assumed that the $u_1,\dots ,u_k$ are independent. We can then set the variational variance-covariance matrix of u to be diagonal and our variational approximation reduces to (11) $$\begin{array}{c}{\underset{¯}{\mathcal{l}}}_{\mathit{u},\text{VA}}\left(\mathit{\beta },\mathit{\mu },{\mathit{\sigma }}^2,{\sigma }_{\text{prior}}^{-2}\right)\\ \hspace{1em}=\sum _{i=1}^n\mathit{X}\left({\mathit{s}}_i\right)\mathit{\beta }+\mathit{Z}\left({\mathit{s}}_i\right)\mathit{\mu }\\ \mathrm{\hspace{1em}}\mathrm{\hspace{1em}}-\sum _{i=n+1}^{n+q}{w}_i\hspace{0.17em}\exp \hspace{0.17em}\left\{\mathit{X}\left({\mathit{s}}_i\right)\mathit{\beta }+\mathit{Z}\left({\mathit{s}}_i\right)\mathit{\mu }+\frac{1}{2}\sum _{r=1}^k{\sigma }_r^2{Z}_r{\left({\mathit{s}}_i\right)}^2\right\}\\ \mathrm{\hspace{1em}}\mathrm{\hspace{1em}}+\frac{1}{2}[k\hspace{0.17em}\ln \hspace{0.17em}\left({\sigma }_{\text{prior}}^{-2}\right)+\left(\sum _{r=1}^k\hspace{0.17em}\ln \hspace{0.17em}\left({\sigma }_r^2\right)\right)\\ \mathrm{\hspace{1em}}\mathrm{\hspace{1em}}-{\sigma }_{\text{prior}}^{-2}\left(\sum _{r=1}^k{\mu }_r^2+{\sigma }_r^2\right)+k].\end{array}$$(11)

Derivation of this can be found in the supplementary material: Section S1. The assumed independence in the random coefficients is key to our use of (11). By Theorem 1, this is the optimal variational approximation if the coefficients are truly independent. Our simulation study (Section 4) aims to demonstrate the quality of the approximation empirically in cases where the true data-generation process is an LGCP with a general GRF.

3.3 Automatic Differentiation

The final component of our methodology for rapid fitting of LGCP models is the use of automatic differentiation (AD; see Griewank 1989) via the Template Model Builder (TMB) package in R (Kristensen et al. 2016). This computes exact derivatives of a given function quickly, which has a few benefits in the maximum likelihood context: gradient information can speed up optimization algorithms (e.g., Shanno 1970); it enables automated approximation of intractable integrals using the Laplace approximation, given that these are a function only of first and second derivatives as in (7); the Hessian matrix can be obtained automatically, for likelihood-based inference about parameters.

TMB software requires the objective function of interest to be coded in C++. We implemented both approximations previously discussed—the variational approximation in (11), and also a Laplace approximation, computed automatically from the integrand using AD techniques.

4 Simulation Study

We undertook a simulation study to compare the computational and statistical efficiency of our proposed methodology against other approximate methods for fitting LGCP within the statistical computing language R. Computations were executed on a high performance computing cluster running a 64 bit version of the CentOS distribution of Linux (Katana, UNSW Sydney 2010) with allowances for 20-60GB of RAM. Comparisons with exact fitting routines were not computationally feasible for the simulation settings here, for example, using a MALA routine (via the R package lgcp: Taylor et al. 2013) was found to require many millions of samples in each simulation to achieve adequate mixing properties, as suggested in Taylor et al. (2015), which would take at least several hours to fit to a moderate-sized dataset, and would add months to computation time in our simulations.

Data were simulated on the unit square as a Poisson point process. The intensity function was a log-linear function of a single, deterministic covariate X, and a zero-mean GRF (ξ): $$\hspace{0.17em}\ln \hspace{0.17em}\lambda \left(\mathit{s}\right)={\beta }_0+{\beta }_{1}X\left(\mathit{s}\right)+\xi \left(\mathit{s}\right).$$

ξ was assumed to have an isotropic Gaussian covariance function, $f_{\xi }=\hspace{0.17em}\exp \hspace{0.17em}\left\{-d{(s,t)}^2/{\rho }^2\right\}$, where ρ is a parameter that controls the correlation range. We additionally explored simulations under an exponential covariance function, $f_{\xi }=\hspace{0.17em}\exp \hspace{0.17em}\left\{-d(s,t)/\rho \right\}$ to reflect a broader range of latent field surfaces—see supplementary material: Section S2.1 for further details. We treated ξ as an unmeasured/unobserved covariate, and hence this formed a LGCP.

In our simulation design we examined the interplay between the spatial scale of the covariate X and the latent random field ξ using a 2 × 2 simulation design, where each of X and ξ was either chosen to be wiggly (W), with a correlation range of $\approx 0.05$, or smooth (S), with a correlation range of $\approx 0.3$ (Figure 2). We expect models to more accurately estimate the true data simulation process when the spatial scales of X and ξ do not coincide.

Fig. 2 The 2 × 2 simulation design showing the scenarios examined. Deterministic functions are used for the covariate ($X(\mathit{s})$) while particular examples of the latent field, $\xi (\mathit{s})$, are shown here. “Smooth” means a range of effect $\approx 0.3$ while “wiggly” means a range of effect $\approx 0.05$, with the spatial domain being the unit square.

We used the spatstat package in R (Turner and Baddeley 2005) to simulate 1000 point patterns from this LGCP within each scenario. We controlled the size of point patterns simulated through β₀ to examine our simulation scenarios with expected number of points $\mathbb{E}\left[N(\mathcal{D})\right]=1000,10,000,100,000$. We standardized the covariate and set the marginal variance of the latent field to 1 so that the magnitudes of model components were roughly equal, to assist convergence in optimization. We fixed the single covariate effect (${\beta }_1=1.25$) as preliminary investigation revealed no change in relative performance of the competing models for a varied fixed effect size. To ensure an adequately fine quadrature approximation, (3), we used a regular square 101 × 101 grid of quadrature points.

We fitted point process regression models to the data, assuming intensity was a log-linear function of X and unobserved ξ, using four different procedures:

INLA via the package INLA, a common approach to fitting a LGCP regression model (Lindgren, Rue, and Lindström 2011, Simpson et al. 2016). The GRF was approximated using piecewise linear basis functions generated via INLA::inla.mesh.create(), applied under default settings to our regular grid of quadrature point locations. We used penalized complexity priors on covariance parameters of the latent field following Fuglstad et al. (2019), setting $\mathrm{Pr}(\sigma >5)=0.01$ and $\mathrm{Pr}(\rho <0.025)=0.01$.

VA the proposed methodology of Section 3, using a variational approximation to the likelihood as in (11) programmed in C++ and fit using TMB and quasi-Newton optimization. FRK was used with either a sparse grid of basis functions (7 × 7) or a dense grid (14 × 14) as in (9).

Lp the proposed methodology of Section 3, using a Laplace approximation to the likelihood as in (7) with the integrand programmed in C++ and fit using TMB and quasi-Newton optimization. FRK was used with either a sparse grid of basis functions (7 × 7) or a dense grid (14 × 14) as in (9).spNNGP (Finley, Datta, and Banerjee 2020) uses MCMC-MALA to fit spatial regression models. Poisson models cannot be fitted in spNNGP, so we assumed binomial data with a logistic link, which approximates an LGCP when the number of quadrature points is large (Warton and Shepherd 2010). The latent field is approximated using a nearest-neighbor Gaussian process (Datta et al. 2016). Priors on the latent field’s covariance parameters are somewhat informative, for details see supplementary material: Section S2.3.

IPP an inhomogeneous Poisson process, that is, with ξ omitted, to study the implications of failing to account for the missing covariate when fitting the model. This was fitted using (2) and (3), each excluding ξ, programmed in C++ and fit using TMB and quasi-Newton optimization.

We expect Lp and VA to perform better when the spatial scale of the basis functions matches the scale of ξ, that is, using a set of basis functions set along a regular 7 × 7 grid should perform better when ξ is smooth (W,S or S,S in Figure 2), whereas using basis functions set along a 14 × 14 grid should perform better when ξ is wiggly (S,W or W,W). We don’t expect spNNGP to perform as well because it only approximates a LGCP, and we only sampled a short chain (taking $10,000$ samples after $20,000$ burn-in) to keep computation time manageable in a simulation setting. As mentioned previously, we also expect difficulties teasing apart effects of X and ξ when they operate at similar spatial scales (S,S or W,W).

We defined how well a method fits by: Root Mean Square Error (RMSE) of $\widehat{{\beta }_1}$; coverage probability of Wald confidence intervals for $\widehat{{\beta }_1}$; and Kullback-Leibler divergence between the fitted and true $\lambda (\mathcal{D})$, which has the form: (12) $$\begin{array}{c}{\int }_{\mathcal{D}}\lambda \left(\mathit{s}\right)\hspace{0.17em}\ln \hspace{0.17em}\left[\frac{\lambda \left(\mathit{s}\right)}{\widehat{\lambda }\left(\mathit{s}\right)}\right]\mathrm{d}\mathit{s}\approx \sum _{i=n+1}^{n+q}{w}_i\lambda \left({\mathit{s}}_i\right)\hspace{0.17em}\ln \hspace{0.17em}\left[\frac{\lambda \left({\mathit{s}}_i\right)}{\widehat{\lambda }\left({\mathit{s}}_i\right)}\right]\\ \hspace{1em}-\sum _{i=n+1}^{n+q}{w}_i\left[\lambda \left({\mathit{s}}_i\right)-\widehat{\lambda }\left({\mathit{s}}_i\right)\right].\end{array}$$(12)

The spNNGP method, fitting a different model, is biased for the intercept term and $\lambda (s)$. We manually corrected for this by adjusting $\widehat{\lambda }(s)$ for spNNGP such that its integral over $\mathcal{D}$ was equal to n, the size of the simulated point pattern. Results summarizing the various scenarios under investigation can be found in Figures 3 and 4. The full and detailed results tables can be found in supplementary material: Section S2.4.

Fig. 3 Average computation times over variously sized point patterns ($\mathbb{E}\left[N(\mathcal{D})\right]=1 000,10,000,100,000$)—see supplementary material: Section S2.2 for details. Models include IPP, INLA, spNNGP and our proposed method using variational (VA) and Laplace (Lp) approximations. The latter two models were fitted both using a coarse regular grid of basis functions, 7 × 7 and a fine regular grid, 14 × 14. For small point patterns our methods were at least 50 times faster than INLA on average. For large point patterns they were as much as 1000 times faster than INLA on average.

Fig. 4 Performance of point and interval estimators for the various models, across simulation scenarios using either a wiggly (W) or a smooth (S) covariate and latent field (labeled with the covariate first, for example, “W,S” means the covariate is wiggly and the latent field smooth). A: Root mean squared error estimating the slope coefficient β₁. B: Kullback-Leibler divergence from the true intensity field $\lambda (\mathcal{D})$ to that fitted by each model. C: Coverage probabilities of 95% Wald intervals for the slope estimate, $\widehat{{\beta }_1}$ (dotted vertical line indicates nominal coverage). We found that, provided we choose the most appropriate basis function configuration for the scenario, our proposed methods performed comparably to INLA in point estimation of β₁, and for Lp, also in estimation of intensity and interval coverage.

Simulations were very computationally intensive, with the generation and fitting of 1000 simulated datasets taking up to 48 thousand hours in some simulation settings.

4.1 Simulation Results

In simulations, our methods of fitting a LGCP model were orders of magnitude faster than INLA, converting computation times in some cases from half an hour to a few seconds (Figure 3). IPP was faster again, but note it does not incorporate spatially correlated errors, which comes at considerable costs to performance. The number of basis functions used in our proposed methods (VA and Lp) appeared to affect computation time, more so than the size of the point pattern being fitted, but computation times remained substantially faster than INLA for all basis function configurations we considered. Despite fitting spNNGP with a relatively short chain length, it was still the slowest method considered, taking over an hour per dataset. In the remaining results we restrict our focus to point patterns with $\mathbb{E}\left[N(\mathcal{D})\right]=1000$ however, trends in results are broadly consistent with the more computationally demanding settings—see supplementary material: Section S2.2.

Figure 4(A)–(B) compares point estimates across fitting techniques. It was typically easier to estimate the parameters when the covariate and latent fields were more distinct, that is, acting at different spatial scales, as we found RMSE to be larger in W,W compared to W,S and S,S compared to S,W. Likewise with confidence intervals (Figure 4(C)), all models had poor coverage when the spatial scales of the covariate and latent fields were similar.

It does appear important that we select an appropriate number of basis functions for the type of latent field we are approximating—we needed many basis functions (14 × 14) when the latent field was wiggly (S,W or W,W) and fewer (7 × 7) when the latent field was smooth (S,S or W,S). Further, the cost of using too few basis functions (7 × 7 for S,W or W,W), is more than the cost of using more than are needed (14 × 14 for S,S or W,S).

When a basis function configuration was used that was appropriate for the simulation setting, the proposed Laplace approximation methods produced generally comparable results to INLA (Figure 4). Variational approximation was competitive in RMSE (Figure 4(A)), but noticeably less accurate as an estimator of intensity when many basis functions were used (Figure 4(B)), and its confidence intervals tended to be too short (Figure 4(C), with supplementary material: Section S2.4 for interval widths).

While choice of the number of basis functions (k) is key, the data could be used to guide this decision by selecting k to maximize the likelihood. We looked at how reliable a data-driven strategy was at choosing the more appropriate basis function configuration for each simulation scenario, and found that the more appropriate configuration could be chosen approximately 99% of the time using a Laplace approximation (see supplementary material: Table S1 for further details). The variational approximation had difficulty in the W,W simulations. We also note that the Laplace approach was sensitive to non-convergence, we found approximately a 1.7% fit failure rate for point patterns with $\mathbb{E}\left[N(\mathcal{D})\right]=1000$—we discuss this further in Section 6.

Table 1 Four-fold cross-validated predicted conditional log-likelihoods of the various models and corresponding computation times.

	INLA	LP	VA	IPP
Predicted $\mathcal{l}(\mathit{\beta }|\mathit{\xi })$	2307.7	2299.8	2301.4	1727.5
Comp. Time	46.79 min	31.88 sec	13.17 sec	0.67 sec

NOTE: We see that our proposed methods (LP and VA) achieved comparable predicted conditional log-likelihoods, in under a minute, compared to INLA that took over 45 min.

5 Real Data Example: Gorilla Nesting Locations

We illustrate our proposed method by analyzing the gorilla nesting dataset (Figure 1(A)) provided in the R package inlabru (Bachl et al. 2019). We wish to model intensity of gorilla nest locations as a log-linear function of three predictors—elevation above sea level; distance to nearest water source; and average temperature, as a three level ordinal factor (Figure 5(A)–(C)).

Fig. 5 A: Spatial covariate describing elevation in meters above sea level. B: Spatial covariate describing distance in meters to nearest water source. C: Factor describing heat category: 1 = Coolest, 2 = Moderate, 3 = Warmest. D: Latent field (conditional on the data) estimated by VA method. E: Latent field (conditional on the data) estimated by INLA. F: Spatially blocked four-fold cross-validation.

We modeled the gorilla nesting point pattern using four methods for comparison—an IPP, as well as an LGCP using the INLA package on R, and both the variational and Laplace versions of our proposed methodology. As in Section 3.2 we used a regular grid of k local bi-square basis functions, but here we chose k to maximize our variational approximation to the marginal likelihood, arriving at k = 63 (Figure 6). Note that models with different values of k are not nested, because basis functions radius is linked to k (Section 3.2), so when choosing k to maximize the likelihood we are in effect using the data to estimate the spatial range of the latent field ξ.

Fig. 6 The variational approximation to the marginal log-likelihood of the proposed LGCP model as a function of increasingly dense regular grids of local bi-square basis functions. The vertical dashed line denotes the basis functions at which the maximum is found. Secondary plot and axis (gray) indicates the cumulative computational time, in seconds, taken to fit models using each of the basis function configurations.

We compared the fitting methods in terms of predictive performance, via a hold-one-out, four-fold cross-validation (CV) that predicts the likelihood in held out test areas of the domain. The folds used are a collection of spatial blocked regions, aligned diagonally across the domain (Figure 5(F)). We compared model fits using predicted conditional log-likelihood: (13) $$\begin{array}{c}\sum _{h=1}^4[\sum _{i=1}^{n_h}\mathit{X}\left({\mathit{s}}_i^{\left(h\right)}\right){\widehat{\mathit{\beta }}}^{\setminus h}+{\widehat{\xi }}^{\setminus h}\left({\mathit{s}}_i^{\left(h\right)}\right)\hspace{1em}\\ \hspace{1em}-\sum _{i=n_h+1}^{n_h+q_h}{w}_i^{(h)}\hspace{0.17em}\exp \hspace{0.17em}\left\{\mathit{X}\left({\mathit{s}}_i^{\left(h\right)}\right){\widehat{\mathit{\beta }}}^{\setminus h}+{\widehat{\xi }}^{\setminus h}\left({\mathit{s}}_i^{\left(h\right)}\right)\right\}]\hspace{1em}\end{array}$$(13) where for fold h we denote as ${\left\{s_i^{(h)}\right\}}_{i=1}^{n_h}$ the point events and ${\left\{s_i^{(h)}\right\}}_{i=n_h+1}^{n_h+q_h}$ the quadrature points in this fold, and $w_i^{(h)}$ are their corresponding weights. The ${\widehat{\mathit{\beta }}}^{\setminus h}$ are fixed effect parameters estimated from the data excluding fold h, and ${\widehat{\xi }}^{\setminus h}$ is likewise the estimated latent field from the training data without fold h.

We found that our proposed method reduced the computation time required for a single model fit from minutes to seconds, with negligible loss of predictive performance (Table 1). Specifically, our four-fold cross-validation procedure took over 45 min to complete on INLA using code from Krainski et al. (2018), whereas the same procedure, using a variational approximation with the appropriate number of basis functions, took just over 10 sec for a comparable predicted conditional likelihood—when run on Katana, UNSW Sydney (2010) with allowance for 20GB of RAM. Note that this computation time does not take into account the basis function search of Figure 6, which took an additional two and a half minutes. INLA used a mesh structure to approximate the latent field, provided by Bachl et al. (2019), which act like basis functions (k = 1479) as described earlier. This larger value for k is one of the reasons for the difference in computation time, explored further in Table 2. In light of Figure 6 this seems excessive, but we found that INLA performed poorly if given a comparable number of basis functions. Conversely, our proposed method performed poorly when using $k\approx 1479$ due to overfitting (Table 2). Omitting a latent field altogether (IPP) was much faster to fit than other methods, but with very low predictive conditional likelihood.

Table 2 The computation time and likelihood for a single model fit for our proposed model for which we toggled on/off key components—a variational approximation for the marginalization (VA, vs. Laplace approximation, Lp); automatic differentiation (AD); sparse basis approximation for the latent field (Sparse/Dense Z); and large number of basis functions (Large k).

Method	Sparse/dense Z	AD	Large k	Comp. Time (sec)	$\mathcal{l}(\mathit{\beta })$
VA	Sparse	✓		4.30	2323.78
Lp	Sparse	✓		16.10	2333.62
VA	Dense	✓		41.80	2271.96
VA	Sparse			93.60	2323.78
VA	Sparse	✓	✓	96.19	2187.01
INLA			✓	718.87	2336.70

NOTE: INLA is also included for comparison.

We also constructed point and interval estimates of the fixed effects parameters in all models (Figure 7). The most conspicuous difference was that the effect of elevation (β₁) was estimated to be larger for IPP, with an unrealistically small standard error. Amongst the LGCP fits, parameter estimates were generally similar, although confidence limits were longer for INLA than Lp or VA, especially for elevation. Note that elevation was the covariate that was smoothest spatially (Figure 5(A)), varying at a comparable spatial scale to our estimated latent fields (Figure 5(D)–(E)), so spatial confounding may be an issue here (analogous to S,S in simulations). This was our worst-case scenario from simulations (Figure 4) and we should not put a lot of weight in inferences about the elevation effect using any of the methods considered.

Fig. 7 Estimated fixed effects 95% Wald confidence intervals for the various models being compared. (a) Elevation effect. (b) Distance to water effect. (c) Heat category—contrast effect of Moderate to Coolest. (d) Heat category—contrast effect of Warmest to Coolest. We see little difference in the estimated fixed effects for the LGCP models (INLA, VA, Lp), so much of the differences in predictive performance came down to the way the latent effect was modeled (or was not modeled for the IPP).

Finally, we want to get a sense of which components of our novel methodology contributed most to the substantial improvements in the speed with which models can be fit. We have proposed three main innovations: using a variational or Laplace approximation; a fixed rank kriging approximation to the latent field, using a relatively small number of sparse basis functions; and automatic differentiation. We looked at the contributions of each of these innovations to computation time for the gorilla nesting LGCP model, by comparing fits for: variational versus Laplace approximations; sparse basis functions (9) versus dense basis functions using thin plate splines as in (as in Tzeng and Huang 2018); optimizing ${\underset{¯}{l}}_{\text{VA}}$ using automatic differentiation via TMB versus a quasi-Newton optimization, using function values only, of an expression for ${\underset{¯}{l}}_{\text{VA}}$ written in R; and using a small number of basis functions (k = 63) rather than a number closer to that used by Bachl et al. (2019) (k = 1470). When we examine the various component models we see that each contributes substantial relative gains in computation speed. The largest speed gains appear to be due to use of AD and use of a small number of basis functions, relative to INLA. The variational approximation yielded the most modest speed gain here, but it still offered an almost four-fold improvement in computation time compared to using a Laplace approximation.

6 Discussion

We have proposed methodology for fitting LGCP models to point patterns that is orders of magnitude faster than INLA, a common method for fitting these models in practice. This speed-up came with negligible loss in statistical efficiency, if a Laplace approximation were used together with an appropriate basis function configuration. We implemented three key techniques—likelihood approximation (via a variational or Laplace approach), rank reduction, and automatic differentiation, each of which contributed to the computational speed-up (Table 2). The largest speed advantages seem to come from rank reduction and from using automatic differentiation (and associated coding of the objective function in C++). We also found considerable computational savings when using sparse basis functions (Cressie and Johannesson 2008) in rank reduction.

The two likelihood approximations considered here had different strengths—the Laplace approach provided more reliable inferences (Figure 4(C)), whereas the variational approach was faster and more stable (supplementary material: Table S1), especially for large k, with some suggestion that it also gave more precise point estimates (see supplementary material). The variational approximation tended to underestimate error in point estimates (Section 4), a result also noticed in Hui et al. (2019), when studying variational approximations to generalized additive models. Our approach has similarities to the Poisson family of these models, if our approximation to the latent field were viewed as a 2D penalized smoother. Hui et al. (2019) suggest using the method of Louis (1982) to obtain a more robust estimator of the observed information matrix. The Louis (1982) method was designed for missing data analysis, so by using this method we are arguing that the variational parameters represent missing data, which is somewhat defensible given that the variational parameters are used to estimate an unobserved random effect. Until such an approach is pursued we recommend that the Laplace approximation be used for inference while the variational approach primarily be used for exploratory work, for example, when doing a basis function search as in Figure 6 (also supplementary material Figure S6). A variational approach could also be used to provide warm starts for parameters in Laplace fits to reduce nonconvergence and speed up optimization for the Laplace approach, which we have found to have some success.

A key decision when implementing our method is choice of number of basis functions, and we recommend using a data-driven approach to decide this as in our application (Figure 6). In simulations, performance of our proposed approach depended strongly on which basis function configuration was used (Figure 4), but we found that the data could be used to quite reliably identify the relevant configuration to use for a given simulation scenario (supplementary material, Table S1). Fitting repeated models at different basis function configurations arguably erodes some of the computational gains of our method. But note that choice of basis function configuration might also be important to performance of INLA, but exploring this issue is computationally prohibitive in many applications—recall that a single INLA fit to our example dataset took more than 10 min (Table 2) compared to the two and a half minutes taken to determine an appropriate basis configuration for our approach (Figure 6).

Decisions were also required concerning the functional form of our basis functions and their knot locations. Cressie and Johannesson (2008) suggest that the precise choice of form of basis function is not critical, nor is its location, as seemed to be the case in our exploratory analyses. The important consideration however is the number of basis functions and their spatial range, as flagged previously, which we determined using the data (Section 5). Cressie and Johannesson (2008) also recommended using multiple resolutions of basis functions, which could readily be implemented here. Including multiple resolutions in a single fit could be particularly useful in settings where there might be multiple missing covariates, operating at different spatial scales, or in massive data situations where it is not practical to fit multiple models to search for an appropriate resolution for basis functions. In our application, we instead searched for an appropriate resolution for basis functions, and in simulations (Section 4), our latent field was stationary and Gaussian hence it operated at a single spatial scale. Note that in our simulations we found point estimation and inference to be difficult when the scale of the covariate and latent field coincided (scenarios S,S and W,W), essentially, basis functions became confounded with the covariate. Including multiple resolutions of basis functions may only increase the chance of such spatial confounding (Hodges and Reich 2010). Perhaps the merits of including multiple resolutions depends whether the main objective is to accurately predict the intensity of point events at certain locations, or to understand their relationship to particular measured predictors.

Our proposed approach to fitting LGCPs offers speed gains so large that we can consider using LGCPs in different ways to how they have been used previously. We have already exploited this speed gain to perform model selection to choose the number of basis functions and their range. This would not be computationally feasible using INLA, where a fit at a single basis function configuration took over 10 min. Other opportunities are in extending the model to handle more complex data structures, such as flexible models for spatio-temporal point patterns, and integrating different data types, as is often desired when modeling species distribution in ecology (Dorazio 2014, Fithian et al. 2015).

Supplementary Materials

All supplemental files are contained in a single archive.

Fast LGCP Modeling—Supplementary.pdf: containing S1: Proofs and working for the variational approximation to the marginal log-likelihood for a proposed rank-reduced LGCP Model. S2: Complete and detailed results for the simulation study (Section 4). S3: Figure of complete basis function search for applied data analysis (Section 5).

Fast LGCP Code and Data.zip: Gorilla nesting dataset, R coding scripts, and saved objects to conduct the applied data analyses presented in the manuscript.

readme.txt: Description of the Supplementary Material archive, including directions for data access and use of code.

Acknowledgments

The authors would like to thank the anonymous reviewers for their contribution to the improved manuscript.

Disclosure Statement

The authors report there are no competing interests to declare.

Additional information

Funding

This work was supported by the Australian Government Research Training Program (RTP) Scholarship, as well as the Australian Research Council Discovery Project scheme (DP180103543).