Function rLGCP fails when trying to generate an LGCP with covariate data using spatstat

Question

I'm trying to simulate a log gaussian cox process based on an image array for the covariate field, and a matern covariance function for the GMRF field, based pretty much exactly on the example given in "Advanced Spatial Modeling with Stochastic Partial Differential Equations Using R and INLA" (https://becarioprecario.bitbucket.io/spde-gitbook/ch-lcox.html). However the code as its written in the online textbook is throwing a warning when it generates the LGCP.

The code used is here:

library(spatstat)

win <- owin(c(0, 3), c(0, 3))
npix <- 300
spatstat.options(npixel = npix)

beta0 <- 3
sigma2x <- 0.2
range <- 1.2
nu <- 1

set.seed(7)

x0 <- seq(0, 3, length=npix)
y0 <- seq(0, 3, length=npix)
gridcov <- outer(x0, y0, function(x,y) cos(x) - sin(y - 2))

beta1 <- 0
sum(exp(beta0 + beta1 * gridcov) * diff(x0[1:2]) * diff(y0[1:2]))

lg.s.c <- rLGCP('matern', im(beta0 + beta1 * gridcov, xcol = x0,
  yrow = y0), var = sigma2x, scale = range / sqrt(8), 
  nu = 1, win = win)
xy.c <- cbind(lg.s.c$x, lg.s.c$y)[, 2:1]
n.c <- nrow(xy.c)

Lam <- attr(lg.s.c, 'Lambda')
rf.s.c <- log(Lam$v)

The returned warning is the images ‘e1’ and ‘e2’ were not compatible and the resulting intensity field rf.s.c has dimensions 562 x 562 whereas the the input x and y arrays are 1x100 (the intensity field should be 100 x 100), and is entirely NA values.

Could someone help me understand this warning and why the example is failing?

Answer 1

This may very well be a bug in spatstat. Recently, new code was added for simulation of LGCPs due to the package RandomFields no longer being on CRAN. We will investigate further and you can track the status of the issue on GitHub: https://github.com/spatstat/spatstat.random/issues/2

Answer 2

Part of the problem here is that the input arguments win and mu are incompatible in this example.

The spatstat code has not automatically repaired the incompatibility, and I guess we could call that a bug.

In your example, win is the square [0,3] x [0,3], while mu was constructed by something like

x0 <- seq(0, 3, length=npix)
y0 <- seq(0, 3, length=npix)
MU <- im(......., xcol=x0, yrow=y0)

The example seems to expect that MU will be an image inside [0,3] x [0,3] but this is not true:

> MU
real-valued pixel image
300 x 300 pixel array (ny, nx)
enclosing rectangle: [-0.0050167, 3.005] x [-0.0050167, 3.005] units

In an image object in spatstat, the coordinates xcol and yrow are the coordinates of the centres of the pixels. For an image that is defined inside [0,3] x [0,3], the coordinates xcol will lie strictly inside the interval (0,3).

> as.im(win)$xcol[1:4]
[1] 0.005 0.015 0.025 0.035

I have now fixed the code in rLGCP (in spatstat.random version 3.2-2.002 available from the GitHub repository) so that it now automatically fixes the incompatibility.

Also, some advice: to make a pixel image where you want the pixel values to be defined by a function f(x,y) of the spatial coordinates, it would be safer and neater to use as.im.function:

meanfun <- function(x,y, beta0, beta1) {
       beta0 + beta1 * (cos(x) - sin(y - 2))
}
MU <- as.im(meanfun, dimyx=300, W=win, beta0=3, beta1=0)

You can specify the raster dimensions using arguments like dimyx, eps or specify the raster by using an existing pixel image as the argument W.