LGCPs - An example in two dimensions
David Borchers and Finn Lindgren
Generated on 2023-12-03Source:
For this vignette we are going to be working with a dataset obtained
spatstat. We will set up a
two-dimensional LGCP to estimate Gorilla abundance.
For the next few practicals we are going to be working with a dataset
obtained from the
contains the locations of 647 gorilla nests. We load the dataset like
data(gorillas, package = "inlabru")
This dataset is a list containing a number of
including the locations of the nests, the boundary of the survey area
INLA mesh - see
details. Extract the the objects we need from the list, into other
objects, so that we don’t have to keep typing
nests <- gorillas$nests mesh <- gorillas$mesh boundary <- gorillas$boundary
Plot the points (the nests).
Fit an LGCP model to the locations of the gorilla nests, predict on the survey region, and produce a plot of the estimated density - which should look like the plot shown below.
Recall that the steps to specifying, fitting and predicting are:
Specify a model, comprising (for 2D models)
coordinateson the left of
~and an SPDE
+ Intercept(1)on the right. Please use the SPDE prior specification stated below.
lgcp( ), passing it (with 2D models) the model components, the
SpatialPointsDataFramecontaining the observed points and the
SpatialPolygonsDataFramedefining the survey boundary using the
predict( ), passing it the fitted model from 2., locations at which to predict and an appropriate predictor specification. The locations at which to predict can be a
SpatialPixelsDataFramecovering the mesh, obtained by calling
fm_pixels(mesh, format = "sp").
You should get a plot like that below (the command below assumes that
the prediction is in an object called
pred <- predict( fit, fm_pixels(mesh, mask = boundary, format = "sp"), ~ data.frame( lambda = exp(mySmooth + Intercept), loglambda = mySmooth + Intercept ) ) pl1 <- ggplot() + gg(pred$lambda) + gg(boundary) + ggtitle("LGCP fit to Points", subtitle = "(Response Scale)") pl2 <- ggplot() + gg(pred$loglambda) + gg(boundary, alpha = 0) + ggtitle("LGCP fit to Points", subtitle = "(Linear Predictor Scale)") multiplot(pl1, pl2, cols = 2)
You can plot the median, lower 95% and upper 95% density surfaces as
follows (assuming that the predicted intensity is in object
ggplot() + gg(cbind(pred$lambda, data.frame(property = "q0.500")), aes(fill = median)) + gg(cbind(pred$lambda, data.frame(property = "q0.025")), aes(fill = q0.025)) + gg(cbind(pred$lambda, data.frame(property = "q0.975")), aes(fill = q0.975)) + coord_equal() + facet_wrap(~property)
Plot the SPDE parameter and fixed effect parameter posteriors.
int.plot <- plot(fit, "Intercept") spde.range <- spde.posterior(fit, "mySmooth", what = "range") spde.logvar <- spde.posterior(fit, "mySmooth", what = "log.variance") range.plot <- plot(spde.range) var.plot <- plot(spde.logvar) multiplot(range.plot, var.plot, int.plot)
Look at the correlation function if you want to:
corplot <- plot(spde.posterior(fit, "mySmooth", what = "matern.correlation")) covplot <- plot(spde.posterior(fit, "mySmooth", what = "matern.covariance")) multiplot(covplot, corplot)
Finally, estimate abundance using the
As a first step we need an estimate for the integrated lambda. The
weight values are contained in the
Given some generous interval boundaries (500, 800) for lambda we can estimate the posterior abundance distribution via
Get its quantiles via
inla.qmarginal(c(0.025, 0.5, 0.975), marginal = list(x = Nest$N, y = Nest$mean)) #>  599.8821 673.8480 753.1823
… the mean via
and plot posteriors:
The true number of nests in 647; the mean and median of the posterior distribution of abundance should be close to this if you have not done anything wrong!