Bias in inlabru predictions compared to INLA

[lionel68]

Inspired by the code to fit LGCP in INLA from: https://becarioprecario.bitbucket.io/inla-gitbook/ch-spatial.html#sec:pointpat, I simulated some data and tried to compare the predictions form INLA and inlabru.

It seems that there is a slight negative bias for the inlabru predictions as compared to INLA, do you have an idea where this could come from?

The R code used:

library(INLA)
library(inlabru)
library(rgeos)
library(spatstat)
library(SDraw)
library(deldir)
library(ggplot2)
library(viridis)

# set a region
coords_pol <- matrix(c(0,0,1,1,0,0,1,1,0,0), ncol = 2, byrow = FALSE)
bndry <- SpatialPolygons(list(Polygons(list(Polygon(coords = coords_pol)), ID = "a")))

# simulate the point process
pp <- rLGCP("matern", nu = 1, var = 2, scale = 0.1, mu = 5, win = owin(c(0,1),c(0,1)))
# put coords in a SpatialPoints format
pp_pts <- SpatialPoints(coords = data.frame(x = pp$x, y = pp$y))
# and in a matrix format
pp_mat <- cbind(pp$x, pp$y)

# create the mesh
mesh <- inla.mesh.2d(loc.domain = bndry@polygons[[1]]@Polygons[[1]]@coords,
                     max.edge = c(0.1, 0.5), offset = c(0.1, 0.5))

# the integration point on the mesh
mesh_mat <- as.matrix(mesh$loc[,1:2])
all_mat <- rbind(pp_mat, mesh_mat)

# dimensions of points and integration points
n_pp <- nrow(pp_mat)
n_int <- nrow(mesh_mat)

# define the spde
spde <- inla.spde2.pcmatern(mesh, 
                            prior.range = c(0.05, 0.1),
                            prior.sigma = c(5, 0.1))

# get weights of the integration points
mytiles <- voronoi.polygons(SpatialPoints(mesh_mat))

w <- sapply(1:length(mytiles), function(p) {
  aux <- mytiles[p, ]  

  if(gIntersects(aux, bndry) ) {
    return(gArea(gIntersection(aux, bndry)))
  } else {
    return(0)
  }
})

# put together
y.pp <- rep(0:1, c(n_int, n_pp))
e.pp <- c(w, rep(0, n_pp))

# projection matrix
lmat <- inla.spde.make.A(mesh, pp_mat)
imat <- Diagonal(n_int, rep(1, n_int))
A.pp <- rbind(imat, lmat)

# generate index for integration points
spde.index <- inla.spde.make.index("spatial", n.spde = spde$n.spde)

# create stack
pp_stack <- inla.stack(data = list(y = y.pp, e = e.pp),
                       A = list(A.pp),
                       effects = list(spde.index),
                       tag = "lgcp")
# predictions
pts_pred <- expand.grid(x = seq(0, 1, length.out = 50),
                        y = seq(0, 1, length.out = 50))

A.pred <- inla.spde.make.A(mesh, loc = as.matrix(pts_pred))
pred_stack <- inla.stack(data = list(y = NA),
                         A = list(A.pred),
                         effects = list(spde.index),
                         tag = "pred")

# put together
joint_stack <- inla.stack(pp_stack, pred_stack)

# fit the model
m_spat <- inla(y ~ 1 + f(spatial, model = spde),
               family = "poisson", data = inla.stack.data(joint_stack),
               E = inla.stack.data(joint_stack)$e,
               control.predictor = list(A = inla.stack.A(joint_stack), compute = TRUE, link = 1))

# plot the spatial predicted fields
idx <- inla.stack.index(joint_stack, 'pred')$data
# MOdel with no covariates
pts_pred$pred <- m_spat$summary.fitted.values[idx, "mean"]
pts_pred$sd_pred <- m_spat$summary.fitted.values[idx, "sd"]

ggplot(pts_pred, aes(x=x,y=y)) +
  geom_raster(aes(fill=log(pred)), interpolate = TRUE) + # seems a bit weird to have to log
  scale_fill_viridis(direction = -1) +
  geom_point(data = as.data.frame(pp_pts))

# try now with inlabru
m_bru <- lgcp(coordinates ~ Intercept + spatial(map = coordinates, model = spde),
              data = pp_pts, samplers = bndry)

lambda <- predict(m_bru, SpatialPoints(pts_pred), ~ spatial + Intercept)
lbd_df <- as.data.frame(lambda)

ggplot(lbd_df,aes(x=x,y=y)) +
  geom_raster(aes(fill=mean), interpolate=TRUE) +
  scale_fill_viridis(direction = -1) +
  geom_point(data = as.data.frame(pp_pts))

plot(log(pts_pred$pred), lbd_df$mean,
     xlab = "INLA prediction",
     ylab = "inlabru prediction")
abline(0,1)

[finnlindgren]

You’re using the log of the posterior mean for the INLA code, but the posterior mean of the log in the inlabru code. These should not be expected to be the same. For a proper comparison I would suggest comparing the posterior medians instead, since those are consistent under transformation.

[lionel68]

Wow, thanks for your super fast answer. You are right, when comparing the prediction on the link scale (exp), then all falls into line