ipoints error when trying distance sampling model

[ASeatonSpatial]

I ran into a problem trying to refit a distance sampling model using the newest version of inlabru. I've created a reprex using the gorillas data.

When I call lgcp I get the following error:

Error in ipoints(NULL, domain[[nm]], name = nm, int.args = int.args) : 
  The number of domain integration points, or a full mesh specification, must be provided
In addition: Warning message:
In proj4string(x) : CRS object has comment, which is lost in output

Reprex:

library(inlabru)
library(INLA)
library(rgeos)
init.tutorial()

data(gorillas)

pp = gorillas$nests
mesh = gorillas$mesh
boundary = gorillas$boundary

# create some point transects
inner_boundary = gBuffer(boundary, width = -0.5)   # don't put them near the boundary
samplers = spsample(inner_boundary, 16, type = "regular")
W = 0.25   # transect radius
samplers_buffered = gBuffer(samplers, width = W, byid = TRUE)

ggplot() + 
  gg(boundary) +
  gg(pp) +
  gg(samplers_buffered) +
  coord_equal()

### SPDE model### 
matern <- inla.spde2.pcmatern(mesh, 
                              prior.sigma = c(0.01, 0.01),    
                              prior.range = c(5, 0.01))

# Try basic plot sampling:
# -----------------------
cmp <- coordinates ~ grf(map = coordinates, model = matern) + Intercept

obs = pp[samplers_buffered,]

fit <-  lgcp(components = cmp, 
             data = obs,
             samplers = samplers_buffered,
             domain = list(coordinates = mesh))    # this works

# Try distance sampling:
# ----------------------

# Log half-normal
log_hn = function(distance, lsig){ 
  -0.5*(distance/exp(lsig))^2
}

# Half-normal
hn <- function(distance, lsig) exp(log_hn(distance, lsig))

# Include r for switching to polar coords, point transect sampling
# - this is what is called in the model components
dsamp = function(distance, lsig){ 
  log(distance) + log_hn(distance, lsig)
}

# calculate distances
samplers$weight <- rep(1, nrow(samplers@coords))  # add weight = 1 column
samplers$plotID <- 1:nrow(samplers@coords)
samplers_buffered$plotID = 1:nrow(samplers)
# Associate each nest with a sampler
IDs <- over(obs, samplers_buffered)
obs$plotID <- IDs$plotID

distances <- c()
for (i in 1:nrow(obs)) {

  pt <- obs[i,]
  transect <- subset(samplers, plotID == pt$plotID)
  distances[i] <- gDistance(pt, transect)
}

# calculate probability of detection for these distances
lsig_true = -2
pdet <- hn(distances, lsig_true)

obs$distance = distances
obs$pdet = pdet

# subsample detectable to create detected
detected = rbinom(nrow(obs), 1, pdet)
sum(detected)

obs_detected = subset(obs, as.logical(detected))

# does this look sensible?
ggplot() +
  gg(boundary) +
  gg(samplers_buffered) +
  gg(obs, colour = "blue") +
  gg(obs_detected, colour = "green") +
  coord_equal() +
  ggtitle("blue not detected, green detected")

# looks okay to me, more detected near centre of transects

# Model components
cmp <- ~ grf(map = coordinates, model = matern) + 
  lsig + Intercept

fml <- coordinates + distance ~ grf +
  dsamp(distance, lsig) + 
  log(2*pi) +   # 2*pi for not knowing angle theta
  Intercept

starting_values <- data.frame(grf = 0, lsig = -1, Intercept = 0)

distance_domain = seq(.Machine$double.eps, W, length.out = 30)

fit = lgcp(components = cmp,
           formula = fml, 
           data = obs_detected,
           domain = list(coordinates = mesh,
                         distance = distance_domain),
           options = list(result = starting_values))    # this throws ipoints error

I actually have an old Rmd file with this something very similar to this code that used to work. The main difference is that now I must say domain = list(coordinates = mesh,...) explicitly, in the old version of inlabru I only needed to specify the distance domain and the coordinates domain was assumed. It seems like this is causing a problem, ipoints doesn't seem to think I have supplied the mesh.

Traceback is

> traceback(max.lines = 10)
9: stop("The number of domain integration points, or a full mesh specification, must be provided")
8: ipoints(NULL, domain[[nm]], name = nm, int.args = int.args)
7: FUN(X[[i]], ...)
6: lapply(nosamp.dim, function(nm) ipoints(NULL, domain[[nm]], name = nm, 
       int.args = int.args))
5: ipmaker(samplers, domain = domain, dnames = response, data = NULL, 
       model = NULL, int.args = options$int.args)
4: (function (family, formula = . ~ ., data = NULL, mesh = NULL, 
       E = NULL, Ntrials = NULL, samplers = NULL, ips = NULL, domain = NULL, 
       include = NULL, exclude = NULL, options = list()) 
   {
       options <- do.call(bru.options, options)
       inla.family <- family
       linear <- as.character(formula)[length(as.character(formula))] == 
           "."
       if (!linear) {
           expr <- parse(text = as.character(formula)[length(as.character(formula))])
    ...
3: do.call(like, c(lhoods, list(options = options)))
2: bru(components, family = "cp", formula = formula, data = data, 
       samplers = samplers, E = E, ips = ips, domain = domain, options = options)
1: lgcp(components = cmp, formula = fml, data = obs_detected, domain = list(coordinates = mesh, 
       distance = distance_domain), options = list(result = starting_values))

> sessionInfo()
R version 4.0.3 (2020-10-10)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 18.04.5 LTS

Matrix products: default
BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.7.1
LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.7.1

locale:
 [1] LC_CTYPE=en_GB.UTF-8       LC_NUMERIC=C              
 [3] LC_TIME=en_GB.UTF-8        LC_COLLATE=en_GB.UTF-8    
 [5] LC_MONETARY=en_GB.UTF-8    LC_MESSAGES=en_GB.UTF-8   
 [7] LC_PAPER=en_GB.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=en_GB.UTF-8 LC_IDENTIFICATION=C       

attached base packages:
[1] parallel  stats     graphics  grDevices utils     datasets  methods  
[8] base     

other attached packages:
[1] rgeos_0.5-5        INLA_20.10.12-1    foreach_1.5.1     
[4] Matrix_1.2-18      inlabru_2.1.14.900 sp_1.4-4          
[7] ggplot2_3.3.2     

loaded via a namespace (and not attached):
 [1] pillar_1.4.6       compiler_4.0.3     iterators_1.0.13  
 [4] tools_4.0.3        digest_0.6.27      lifecycle_0.2.0   
 [7] tibble_3.0.4       gtable_0.3.0       lattice_0.20-41   
[10] pkgconfig_2.0.3    rlang_0.4.8        rstudioapi_0.11   
[13] yaml_2.2.1         rgdal_1.5-18       withr_2.3.0       
[16] dplyr_1.0.2        MatrixModels_0.4-1 generics_0.1.0    
[19] vctrs_0.3.4        grid_4.0.3         tidyselect_1.1.0  
[22] glue_1.4.2         R6_2.5.0           purrr_0.3.4       
[25] farver_2.0.3       magrittr_1.5       scales_1.1.1      
[28] codetools_0.2-16   ellipsis_0.3.1     splines_4.0.3     
[31] colorspace_1.4-1   labeling_0.4.2     munsell_0.5.0     
[34] crayon_1.3.4    

[finnlindgren]

Yes, the new code has tightened the requirements, so that it's not enough to provide a number sequence; "The number of domain integration points, or a full mesh specification, must be provided" At the moment, that just means you can wrap the sequence in a call to inla.mesh.1d; distance_domain = inla.mesh.1d(seq(.Machine$double.eps, W, length.out = 30))

[finnlindgren]

The reason is that in some cases it's impossible to distinguish between a "number of integration points" and "a domain that only has one point" (although I think there is a bug in inla.mesh.1d that might prevent it from creating such a mesh at the moment).

[ASeatonSpatial]

Ah okay I was assuming it was coordinates but it was the distance domain. The model fitted with distance_domain = inla.mesh.1d(seq(.Machine$double.eps, W, length.out = 30)), might be wroth expanding on the documentation for the domain argument in the lgcp docs

[finnlindgren]

I think I can make the error message include the name of the problem domain to make it easier. The integration code is still in some flux.

[finnlindgren]

The relevant code has drastically changed, and the whole integration code is being rewritten, so this particular issue is no longer relevant.