Struggling to fit covariates with LGCP

[lvclarke]

Hi,

I'm trying to include a habitat covariate within my LGCP, but running into issues. Is there a problem with the code below?

Specify hazard-rate detection function

hn <- function(distance, lsig) {1 - exp(-(distance / (exp(lsig)))^-1)}

Define a (Matern) covariance function for the SPDE

matern <- inla.spde2.pcmatern(mesh,
prior.sigma = c(1, 0.01),
prior.range = c(10, 0.9))

Specify components of model

cmp2 <- coordinates ~
-1 +
Habitat(Covariates$Habitat, model = "factor_full") +
mySPDE(coordinates, model = matern)

Specify formula

form <- coordinates + distance ~ mySPDE + log(hn(distance, lsig)) + Intercept + log(2)

fit2 <- lgcp(
components = cmp2,
data = nests,
samplers = transects,
domain = list(
coordinates = mesh,
distance = INLA::inla.mesh.1d(seq(0, W, length.out = 30)) ),
formula = form, options = list(bru_initial = list(lsig = 2),
bru_verbose=4))

The error message I get when running this is:
Error in sort.int(x, na.last = na.last, decreasing = decreasing, ...) : 'x' must be atomic

[andyhouldcroft]

Hi,

Your components definition appears to be missing the lsig(1) and Intercept(1) parameters, not certain if this will fix the error but worth fixing to see.

Andrew

[finnlindgren]

I suspect an issue with the covariate, but need more info; the traceback() after you get the error, the sessionInfo(),
str(Covariates$Habitat), and str(mesh).

[lvclarke]

cmp2 <- coordinates ~
-1 +
Habitat(Covariates$Habitat, model = "factor_full") +
mySPDE(coordinates, model = matern) +
lsig(1) +
Intercept(1)

fit2 <- lgcp(
components = cmp2,
data = nests,
samplers = transects,
domain = list(
coordinates = mesh,
distance = INLA::inla.mesh.1d(seq(0, W, length.out = 30)) ),
formula = form,
options = list(bru_initial = list(lsig = 2),
bru_verbose=4))

traceback()
15: stop("'x' must be atomic")
14: sort.int(x, na.last = na.last, decreasing = decreasing, ...)
13: sort.default(unique(unlist(inp)), na.last = NA)
12: sort(unique(unlist(inp)), na.last = NA)
11: make_unique_inputs(inp_, allow_list = TRUE)
10: add_mapper(component$main, label = component$label, lhoods = lh, env = component$env, require_indexed = FALSE)
9: add_mappers.component(components[[k]], lhoods)
8: add_mappers(components[[k]], lhoods)
7: add_mappers.component_list(object, lhoods = lhoods)
6: add_mappers(object, lhoods = lhoods)
5: component_list.list(components, lhoods)
4: component_list(components, lhoods)
3: bru_model(components, lhoods)
2: bru(components, lik, options = options)
1: lgcp(components = cmp2, data = nests, samplers = transects, domain = list(coordinates = mesh, distance = INLA::inla.mesh.1d(seq(0, W, length.out = 30))), formula = form, options = list(bru_initial = list(lsig = 2), bru_verbose = 4))

Sorry I've had trouble c+p code so hopefully the screenshots below work for session info, str(Covariates$Habitat) and str(mesh)

[finnlindgren]

Thanks! Does it run without the covariate component? The error trace seems to indicate it's having trouble identifying the possible values of the covariate. The main issue I think is that there is no eval_spatial.RasterLayer, so inlabru doesn't realise it would have needed to do spatial evaluaton of the covariate. Fortunately, the terra style of rasters are supported, so I think you can solve the issue by simply converting it before running the model. I'm guessing that your Covariates object is a full raster object, in which case you can probably do this:

Covariates <- terra::rast(Covariates)

Otherwise you'd need to convert the individual layer.
I don't think you need to change anything else.

[finnlindgren]

In case I decide to add legacy support for the raster package to eval_spatial, could you provide the str(Covariates) output? That would help figure out what classes are involved.

[lvclarke]

Hi Finn - thanks very much for this! Your fix worked. Pasted the screenshots below with the str(Covariates) before fix was applied:

[finnlindgren]

OK, thanks! Btw, I've just updated inlabru development version in a major way that adds further sf support (the vignettes are in the process of being regenerated for the website, with some of them now using sf and terra instead of sp). Might be safest to wait a day or two until I iron out any further bugs, but then I would recommend upgrading. There are quite a few speedups, e.g. for predict() calls.

[lvclarke]

Great thank you for the heads up!