Detection function extension

[david-borchers]

Implement the flexible detection function model using SPDEs that was used in the Annals of Applied Statistics paper.

[fbachl]

Have a look at the essmod repository: essmod/tex/paper1_code/example_bwhales_semipar.R :)

[david-borchers]

Thanks Fabian.

I do not understand all the INLA commands, which leads me to ask these questions:

1. Presumably we’d want inla.spde2.matern to be replaced by inla.spde2.pcmatern, and also B.tau and B.kappa would be replaced by prior.range and prior.sigma?
2. I have never used the Cmatrix or the A.msk arguments and don’t know what "model='generic0’ “ means. Is there a simpler way of setting up the spdf in the distance dimension, using only the kinds of commands that we use in the inlabru workshop practicals?

I tried fitting the model but it fell over here:

fit = lgcp(components = cmp, data = etp$bwhale,

• samplers = etp$transect,
• domain = list(distance = c(0,6)),
• options = list(inla.call = "remote")) iinla() iteration 1 [ max: 1 ]. Error in (function (data, model, stackmaker, n = 10, result = NULL, family, : INLA returned message: The inla-program exited with an error. Unless you interupted it yourself, please rerun with verbose=TRUE and check the output carefully. If this does not help, please contact the developers at help@r-inla.org<mailto:help@r-inla.org>. In addition: Warning message: In inla.inlaprogram.has.crashed() : The inla-program exited with an error. Unless you interupted it yourself, please rerun with verbose=TRUE and check the output carefully. If this does not help, please contact the developers at help@r-inla.org<mailto:help@r-inla.org>.

[finnlindgren]

I’ll take a look at it; IIRC this is a kappa=0 model so pcmatern won’t work.

[finnlindgren]

In the INLA package, objects of type inla.model.class have elements model and f with parameters for the INLA::f() function, so that f(name, model=mod) is interpreted as f(name, model=mod$model, par1 = mod$f$par1, par2 = mod$f$par2).

We could implement the same technique in inlabru, with its own model class:

# Helper function for constructing the semiparametric model:
model.detect.semipar <- function(max.dist, n) {
  knots = seq(0, max.dist, length.out = n+1)
  mesh <- inla.mesh.1d(knots, degree = 2, boundary=c("neumann", "free"))
  spde <- inla.spde2.matern(mesh, B.tau=cbind(0), B.kappa=cbind(-10))
  Q <- inla.spde.precision(spde, theta=c())[-1,-1,drop=FALSE]
  msk = c(FALSE,rep(TRUE,length(knots)-1))
  model <- list(model = "generic0",
                Cmatrix = Q,
                A.msk = msk,
                mesh = mesh)
  class(model) <- "inlabru_model_class"
  model
}

A user would then only need to do this:

mod <- model.detect.semipar(6, 29)
cmp2 = coordinates +distance ~ spdf(map = distance, model = mod) + ...

This could also be combined with the INLA feature, and the inlabru parameters encapsulated:

model.inla <- list(model = "generic0", f = list(Cmatrix = Q))
class(model.inla) <- "inla.model.class"
model <- list(model = model.inla,
    inlabru = list(A.msk = msk, mesh = mesh))
class(model) <- "inlabru_model_class"

This version would more clearly separate basic INLA parameters from the special inlabru parameters. Care would have to be taken to ensure that the variable scoping in the formula handling of both inlabru and INLA finds the correct variables.

In the current inlabru version that doesn't have that feature, the object made with the helper function above can be use like this:

cmp2 = coordinates + distance ~ spdf(map = distance, 
                                   model = mod$model, 
                                   Cmatrix = mod$Cmatrix, A.msk = mod$msk, mesh = mod$mesh) + 
  spat(map = coordinates,
       model = inla.spde2.matern(etp$mesh)) +
  Intercept

[finnlindgren]

The "encapsulating" approach is the way to go; we need to handle hyper.default (the prior, possibly a default prior, set when constructing the model object) and hyper (a prior set in the component using the model object), and other special features, and keeping INLA::f() parameters separate from inlabru parameters makes that easier.

[finnlindgren]

Possible package code:

model.detect.semipar <- function(max.dist, n) {
  knots = seq(0, max.dist, length.out = n+1)
  mesh <- inla.mesh.1d(knots, degree = 2, boundary=c("neumann", "free"))
  spde <- inla.spde2.matern(mesh, B.tau=cbind(0), B.kappa=cbind(-10))
  Q <- inla.spde.precision(spde, theta=c())[-1,-1,drop=FALSE]
  msk = c(FALSE,rep(TRUE,length(knots)-1))
  mod.inla <- list(model = "generic0", f = list(Cmatrix = Q))
  class(mod.inla) <- "inla.model.class"
  model <- list(model = mod.inla, ## This is what INLA will see
                bru = list(A.msk = msk, ## Additional bru parameters
                           mesh = mesh))
  ## Keep track of what kind of model it is:
  class(model) <- c("bru_model_detect_semipar", "bru_model_class")
  model
}

User code:

mod <- model.detect_semipar(6, 29) # 29 intervals, to match the 30 knots in the example
cmp <- ... + spdf(map = distance, model = mod) + ...

In inlabru v2.1.4 the model object from the helper function above can be used like this:

cmp2 <- ... + spdf(map = distance, 
                   model = mod$model,
                   A.msk = mod$bru$msk,
                   mesh = mod$bru$mesh) + ...

[finnlindgren]

This might be combined with #4