fit_lgcp: Spatial or spatiotemporal log-Gaussian Cox process (LGCP)

In cmjt/stelfi: Hawkes and Log-Gaussian Cox Point Processes Using Template Model Builder

Spatial or spatiotemporal log-Gaussian Cox process (LGCP)

Description

Fit a log-Gaussian Cox process (LGCP) using Template Model Builder (TMB) and the R_inla namespace for the SPDE-based construction of the latent field.

Usage

fit_lgcp(
  locs,
  sf,
  smesh,
  tmesh,
  parameters,
  covariates,
  tmb_silent = TRUE,
  nlminb_silent = TRUE,
  ...
)

Arguments

locs	A data.frame of x and y locations, 2 \times n. If a spatiotemporal model is to be fitted then there should be the third column (t) of the occurrence times.
sf	An sf of type POLYGON specifying the spatial region of the domain.
smesh	A Delaunay triangulation of the spatial domain returned by fmesher::fm_mesh_2d().
tmesh	Optional, a temporal mesh returned by fmesher::fm_mesh_1d().
parameters	A named list of parameter starting values: beta, a vector of fixed effects coefficients to be estimated, \beta (same length as ncol(covariates) + 1 ); log_tau, the \textrm{log}(\tau) parameter for the GMRF; log_kappa, \textrm{log}(\kappa) parameter for the GMRF; atanh_rho, optional, \textrm{arctan}(\rho) AR1 temporal parameter. Default values are used if none are provided. NOTE: these may not always be appropriate.
covariates	Optional, a matrix of covariates at each smesh and tmesh node combination.
tmb_silent	Logical, if TRUE (default) then TMB inner optimisation tracing information will be printed.
nlminb_silent	Logical, if TRUE (default) then for each iteration nlminb() output will be printed.
...	optional extra arguments to pass into stats::nlminb().

Details

A log-Gaussian Cox process (LGCP) where the Gaussian random field, Z(\boldsymbol{x}), has zero mean, variance-covariance matrix \boldsymbol{Q}^{-1}, and covariance function C_Z. The random intensity surface is \Lambda(\boldsymbol{x}) = \textrm{exp}(\boldsymbol{X}\beta + G(\boldsymbol{x}) + \epsilon), for design matrix \boldsymbol{X}, coefficients \boldsymbol{\beta}, and random error \epsilon.

Shown in Lindgren et. al., (2011) the stationary solution to the SPDE (stochastic partial differential equation) (\kappa^2 - \Delta)^{\frac{\nu + \frac{d}{2}}{2}}G(s) = W(s) is a random field with a Matérn covariance function, C_Z \propto {\kappa || x - y||}^{\nu}K_{\nu}{\kappa || x - y||}. Here \nu controls the smoothness of the field and \kappa controls the range.

A Markovian random field is obtained when \alpha = \nu + \frac{d}{2} is an integer. Following Lindgren et. al., (2011) we set \alpha = 2 in 2D and therefore fix \nu = 1. Under these conditions the solution to the SPDE is a Gaussian Markov Random Field (GMRF). This is the approximation we use.

The (approximate) spatial range = \frac{\sqrt{8 \nu}}{\kappa} = \frac{\sqrt{8}}{\kappa} and the standard deviation of the model, \sigma = \frac{1}{\sqrt{4 \pi \kappa^2 \tau^2}}. Under INLA (Lindgren and Rue, 2015) methodology the practical range is defined as the distance such that the correlation is \sim 0.139.

Value

A list containing components of the fitted model, see TMB::MakeADFun. Includes

• par, a numeric vector of estimated parameter values;

• objective, the objective function;

• gr, the TMB calculated gradient function; and

• simulate, a simulation function.

References

Lindgren, F., Rue, H., and Lindström, J. (2011) An explicit link between Gaussian fields and Gaussian Markov random fields: the stochastic partial differential equation approach. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 73: 423–498.

Lindgren, F. and Rue, H. (2015) Bayesian spatial modelling with R-INLA. Journal of Statistical Software, 63: 1–25.

See Also

fit_mlgcp and sim_lgcp

Examples

### ********************** ###
## A spatial only LGCP
### ********************** ###
if(requireNamespace("fmesher")) {
data(xyt, package = "stelfi")
domain <- sf::st_as_sf(xyt$window)
locs <- data.frame(x = xyt$x, y = xyt$y)
bnd <- fmesher::fm_as_segm(as.matrix(sf::st_coordinates(domain)[, 1:2]))
smesh <- fmesher::fm_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)))
### ********************** ###
## A spatiotemporal LGCP, AR(1)
### ********************** ###
ndays <- 2
locs <- data.frame(x = xyt$x, y = xyt$y, t = xyt$t)
w0 <- 2
tmesh <- fmesher::fm_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))
}

cmjt/stelfi documentation built on Oct. 25, 2023, 2:53 p.m.