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.
fit_lgcp()
functionargs(fit_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)
get_fields(fit, smesh) |> show_field(smesh = smesh, sf = domain) + ggplot2::theme_void() show_lambda(fit, smesh = smesh, sf = domain) + ggplot2::theme_void()
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)
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")
sim_lgcp()
functionSimulating 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()
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.