Data integration help (joint likelihood)

[simonrolph]

I have put this in here as a discussion rather than issue - hope that's appropriate.

I want to fit a model with joint likelihood. I am testing with simulated data. Here's a rendered markdown to help illustrate (it was html but I can't seem to upload html to I've saved it as a pdf)
MAMBO sim framework tests.pdf

I've simulated species with NLMR (for the environment) and virtualspecies (for the species). Then simulating sampling via 3 methods: presence only (eg. citizen science), presence/absence (eg. a survey), presence/absence via something similar to an automated monitoring system eg. acoustic recorder which has very few locations but multiple sampling events.

So this is what I have:

Key:

• Back solid circles - citizen science PO data
• Black ‘X’ - presence timed counts
• Black ‘A’ - presence automated trap
• Red ‘X’ - absence from timed count
• Red ‘a’ - absence from automated trap

When I fit separate models these are the result (each row is a different model):

They vary in success but that's expected, especially as the automated trap only has 4 locations and so can sometimes predict quite poorly.

However I thought this would be a good scenario where data integration would work well however the integrated models don't look great (each of these is a different integrated combining different data sets, then all datasets):

This is my fully integrated model, the other are simply with dataset removed from the pipe

mod_all <- distribution(background) %>%
  add_predictors(env = env_variables,
                 transform = "scale", 
                 derivates = "none") %>%
  add_biodiversity_poipo(sp1_po,
                       name = "Citizen sci po data",
                       field_occurrence = "detected") %>%
  add_biodiversity_poipa(sp1_pa_ami,
                         name = "AMI trap pa data",
                         field_occurrence = "detected") %>%
  add_biodiversity_poipa(sp1_pa_timed,
                         name = "Timed count data",
                         field_occurrence = "detected") %>%
  engine_inlabru() %>%
  train(runname =  "sp1_po_pa_ami_timed",
        verbose = T,
        aggregate_observations = F)

Why are the predicted values so large, up to 10000? I found when I log the predicted values the map actually looks okay so do I have a scaling issue?

The data integration docs page seems to imply that I shouldn't need to be generating pseudoabsences as it's modelling with a PPM. Is that correct?

Do I need to do some more additional set-up, parameterisation, priors etc.?

Any help would be appreciated, thanks in advance!

[Martin-Jung]

Hi Simon,

Cool that you testing out the package for your usecase. Here are some quick reactions.

• Generally fitting integrated SDMs with Inlabru is hard, especially as INLA uses a mesh to link observations to predictors and critically such predictions vary depending on the density of the mesh, see also here and there also general questions whether joint likelihood models as fitted through INLA even help. To me the issue is also that INLA as an modelling paradigm does not support regularizing Bayesian priors (e.g. 'Horseshoe', 'Laplace' or 'Spike-and-Slab' priors), such as supported for example by engine_breg() or engine_stan(). Without such regularization, the SDMs fitted with INLA(bru) are more less similar to a GLM in structure (although they combine different likelihoods together). There is some plan to have a INLA alternative fully programmed for stan in engine_stan() in ibis.iSDM, but owing to time constraints I haven't yet implemented this as of now. Might happen when we have a project where we need it.

Why are the predicted values so large, up to 10000? I found when I log the predicted values the map actually looks okay so do I have a scaling issue?

• This is a common issue when fitting PPMs and predicting at the response scale (not only for INLA but sometimes also other engines), which includes measurement uncertainty. Usually most people in papers using INLA don't predict for the response but rather at the predictor scale (you can enable that by setting engine_inlabru(type = "predictor"). In addition I would also check for your coefficients via coefficient(mod_all) or summary(mod_all). If there are some values that seem rather large, you will usually find that the prediction will improve considerably if those variables are removed (mod |> rm_predictor(variable)). Lastly, I also note that you are currently fitting INLA without any spatial latent effect, which is another thing you could try. Here adding for example a SPDE to the model could help (see here the help pages from the underlying package). In ibis.iSDM this simply works by adding latent_spatial(method = "spde") to your modelling pipe, although as always parameter settings matter and model convergence is not assured.

The data integration docs page seems to imply that I shouldn't need to be generating pseudoabsences as it's modelling with a PPM. Is that correct?

• Correct, all PPMs regardless of engine automatically set background points (not pseudo-absence points) during estimation. See here for an overview of the how and why. By default in ibis.iSDM the settings specified in ibis_options()$ibis.pseudoabsence are used for background sampling. If you want to use alternative settings, simply define them via pseudoabs_settings() and add them to a biodiversity call (e.g., add_biodiversity_poipo(pseudoabsence_settings = abs)).

Do I need to do some more additional set-up, parameterisation, priors etc.?

• It is technically working, but as mentioned before you could try a few additional tweaks and see if that helps with convergence. Defining better mesh parameters can help. If you want to set a (gaussian) prior, you can do so for example via add_priors(priors(INLAPrior("variable",type = "normal", hyper = c(2, 0.001)))). I also often observed that joint-likelihood models have rather skewed posterior distributions, so make sure to also visualize the other moments via plot(mod_all, "q50"), plot(mod_all, "q05") or plot(mod_all, "mode")

Good luck, Fitting models with INLA is a pain :)

[simonrolph]

Many thanks for your reply Martin. Lots of things to try out!

[simonrolph]

With stan, it's currently unclear how to fit a stan model from the docs.

On https://iiasa.github.io/ibis.iSDM/reference/engine_stan.html it says that 'Stan code has to be written separately and this function acts as compiler to build the stan-model' but not sure where I define this.

There's this:
https://iiasa.github.io/ibis.iSDM/reference/run_stan.html but how do I put the model into the engine?

If I just substitute in engine_stan so my model looks like this (which as described above I expect to fail):

Background extent: 
     xmin: -1.5, xmax: 53,
     ymin: 0, ymax: 54
   projection: +proj=longlat +datum=WGS84 +no_defs
 --------- 
Biodiversity data:
   Point - Presence only <73 records>
 --------- 
  predictors:     env1, env2, env3, env4 (4 predictors)
  priors:         <Default>
  latent:         None
  log:            <Console>
  engine:         <STAN>

and try train(mod) then I get this rather unclear error which should probably say something about model not being defined?:
Error: object 'x' not found

[Martin-Jung]

Thanks, the error with object x is very likely a bug. There is no way to add a custom stan model to engine_stan and joint likelihood modelling with engine_stan() as written above not yet implemented as of now.

[simonrolph]

Not really a conclusion to this thread as lots of things that I am yet to try out but just some results from trying type=predictor.

Using type = predictor does seem to give more useful predictions (in terms of eyeballing the real species distribution) however the cv seems to be less useful. Also noticing that the cv of the type = response seems to be the same as the sd of the type=predictor fit.

Before:

After:

And here's summary(mod_all):

summary(mod_all)
# A tibble: 8 x 8
  variable                       mean      sd     q05    q50      q95   mode   kld
  <chr>                         <dbl>   <dbl>   <dbl>  <dbl>    <dbl>  <dbl> <dbl>
1 Intercept                     1.20  27.4    -43.8    1.20  46.2      1.20      0
2 Intercept_sp1_po_poipo        1.20  27.4    -43.8    1.20  46.2      1.20      0
3 Intercept_sp1_pa_timed_poipa  1.20  27.4    -43.8    1.20  46.2      1.20      0
4 Intercept_sp1_pa_ami_poipa    1.20  27.4    -43.8    1.20  46.2      1.20      0
5 env1                         -1.53   0.0847  -1.67  -1.53  -1.39    -1.53      0
6 env2                         -0.469  0.0769  -0.595 -0.469 -0.342   -0.469     0
7 env3                          0.884  0.0615   0.783  0.884  0.985    0.884     0
8 env4                         -0.102  0.0661  -0.211 -0.102  0.00651 -0.102     0

Attempting to add SPDE resulted in model failure:

Error in inla.set.hyper(model, section, hyper = hyper.default, debug = debug) : 
  Unknown keyword in `hyper' ` theta2 '. Must be one of  theta theta1 b beta .

 *** inla.core.safe:  inla.program has crashed: rerun to get better initial values. try=1/2 
Error in inla.set.hyper(model, section, hyper = hyper.default, debug = debug) : 
  Unknown keyword in `hyper' ` theta2 '. Must be one of  theta theta1 b beta .

 *** inla.core.safe:  inla.program has crashed: rerun to get better initial values. try=2/2 
Error in inla.set.hyper(model, section, hyper = hyper.default, debug = debug) : 
  Unknown keyword in `hyper' ` theta2 '. Must be one of  theta theta1 b beta .
Warning in iinla(model = info[["model"]], lhoods = info[["lhoods"]], options = info[["options"]]) :
  iinla: Problem in inla: Error in inla.core.safe(formula = formula, family = family, contrasts = contrasts,  : 
  *** Fail to get good enough initial values. Maybe it is due to something else.
iinla: Problem in inla: 1: fit_and_train_model(env_variables, po_datasets = list(sp1_po = sp1_po), 
       pa_datasets = list(sp1_pa_timed = sp1_pa_timed, sp1_pa_ami = sp1_pa_ami))
2: train(mod, runname = runname, verbose = F, aggregate_observations = F)
3: train(mod, runname = runname, verbose = F, aggregate_observations = F)
4: x$engine$train(model, settings)
5: res(x, ...)
6: try({
       fit_bru <- inlabru::bru(components = comp, likelihoods, options = options)
   }, silent = FALSE)
7: tryCatch(expr, error = function(e) {
       call <- conditionCall(e)
       if (!is.null(call)) {
           if (identical(call[[1L]], quote(doTryCatch))) 
               call <- sys.call(-4L)
           dcall <- deparse(call, nlines = 1L)
           prefix <- paste("Error in", dcall, ": ")
           LONG <- 75 [... truncated]
iinla: Giving up and returning last successfully obtained result for diagnostic purposes.
Error in res(x, ...) : 
  Model did not converge. Try to simplify structure and check priors!

[Martin-Jung]

Brief response on top of my head. I think the similar CV/SD makes sense as it is calculated based on the ratio of the mean and sd from the posterior, and these relative differences should differ if a predictor or response unit is used. Will reassess when I have a moment.
For the spde this is likely related to the mesh specification and/or the prior setting. SPDE priors can be set like this priors(INLAPrior('spde','prior.range', c(0.01,0.05)), INLAPrior('spde','prior.sigma', c(2,0.05))), but I would first try altering the mesh parameters in engine_inlabru.