Chapter 4 Log-Gaussian Cox process

A Log-Gaussian Cox Process (LGCP) is a doubly stochastic spatial point process. In the simplest case, the intensity of the point process over space is given by: \[\Lambda(s) = \text{exp}(\beta_0 + G(s) + \epsilon)\] where \(\beta_0\) is a constant, known as the intercept, \(G(s)\) is a Gaussian Markov Random Field (GMRF) and \(\epsilon\) an error term.

It is conventional to use Matérn covariance function to define the covariance of the random field. This takes two parameters \(\tau\) and \(\kappa\), commonly reported as \(r=\frac{\sqrt{8}}{\kappa}\) and \(\sigma=\frac{1}{\sqrt{4\pi\kappa^2\tau^2}}\), where \(r\) is the range and \(\sigma\) is the standard deviation.

4.1 The fit_lgcp() function

args(fit_lgcp)

## function (locs, sf, smesh, tmesh, parameters, covariates, tmb_silent = TRUE, 
##     nlminb_silent = TRUE, ...) 
## NULL

4.1.1 Fitting a spatial only LGCP

data(xyt, package = "stelfi")
domain <- sf::st_as_sf(xyt$window)
locs <- data.frame(x = xyt$x, y = xyt$y)
bnd <- INLA::inla.mesh.segment(as.matrix(sf::st_coordinates(domain)[, 1:2]))
smesh <- INLA::inla.mesh.2d(boundary = bnd,
                                max.edge = 0.75, cutoff = 0.3)
fit <- fit_lgcp(locs = locs, sf = domain, smesh = smesh,
                                parameters = c(beta = 0, log_tau = log(1),
                                               log_kappa = log(1)))
get_coefs(fit)

##             Estimate Std. Error
## beta       2.4481848 0.09034657
## log_tau   -1.3156570 0.33288053
## log_kappa  0.9547072 0.29440167
## range      1.0887319 0.32052449
## stdev      0.4047190 0.05027126

get_fields(fit, smesh) |>
    show_field(smesh = smesh, sf = domain) + ggplot2::theme_void()

show_lambda(fit, smesh = smesh, sf = domain) + ggplot2::theme_void()

4.1.2 Spatiotemporal LGCP

The LGCP model can also be used for spatiotemporal modelling where there is autoregressive temporal dependence.

To achieve this, we choose an arbitrary number of time knots. The equation for an AR(1) process is as follows: \[\Lambda_i(s) = \text{exp}(\beta_0 + G_i(s) + \epsilon)\] where \(i\) indexes the time knot. \(\Lambda_i(s)\) is the field intensity at time knot \(i\), and \(G_i(s)\) the GMRF at the same time knot. Each \(G_i(s)\) shares common values for \(\tau\) and \(\kappa\).

Successive random fields are correlated through the formula \[G_i(s)=\rho G_{i-1}(s) + \epsilon_i\] where \(\rho\) is a constant between -1 and +1, and \(\epsilon_i\) is normally distributed with mean 0.

ndays <- 2
locs <- data.frame(x = xyt$x, y = xyt$y, t = xyt$t)
bnd <- INLA::inla.mesh.segment(as.matrix(sf::st_coordinates(domain)[, 1:2]))
w0 <- 2
smesh <- INLA::inla.mesh.2d(boundary = bnd,
                                max.edge = 0.75, cutoff = 0.3)
tmesh <- INLA::inla.mesh.1d(seq(0, ndays, by = w0))
fit <- fit_lgcp(locs = locs, sf = domain, smesh = smesh, tmesh = tmesh,
                    parameters = c(beta = 0, log_tau = log(1),
                                   log_kappa = log(1), atanh_rho = 0.2))
get_coefs(fit)

##             Estimate Std. Error
## rho        0.3100841 0.20738457
## beta       1.6852264 0.13685037
## log_tau   -1.0683016 0.22402076
## log_kappa  0.4785904 0.22707189
## range      1.7526525 0.39797812
## stdev      0.5087488 0.05765947

show_lambda(fit, smesh = smesh, sf = domain, tmesh = tmesh, timestamp = 1) +
    ggplot2::theme_void() + ggplot2::ggtitle("t = 1")

show_lambda(fit, smesh = smesh, sf = domain, tmesh = tmesh, timestamp = 2) +
    ggplot2::theme_void() + ggplot2::ggtitle("t = 2")

4.2 The sim_lgcp() function

Simulating a spatiotemporal LGCP

Option 1 using sim_lgcp

parameters <- c(beta = 1, log_tau = log(1), log_kappa = log(1), atanh_rho = 0.2)
simdata <- sim_lgcp(parameters = parameters, sf = domain, smesh = smesh, tmesh = tmesh)
simdata$x[,1] |>
    show_field(smesh = smesh) +
    ggplot2::theme_void()

Option 2 directly from the fitted model

simdata <- fit$simulate()
simdata$x[,1] |>
    show_field(smesh = smesh) +
    ggplot2::theme_void()