Issues with WAIC and DIC calculations in integrated SDM

[Peter Stewart]

Hi folks,

Myself and a few others have been encountering issues when using the PointedSDMs package (based on inlabru - https://github.com/PhilipMostert/PointedSDMs) to fit integrated species distribution models. Our models have been producing WAIC values which are "-Inf" or "NaN", or which are extremely large (e.g. 2.56e+182). The DIC values for these models also tend to be NA. The issue does not appear to be confined to a particular dataset or specific covariates. 

We have been discussing the issue over on GitHub (https://github.com/PhilipMostert/PointedSDMs/issues/5) but haven't yet managed to pinpoint the cause, so I thought I would also post here in case anybody is able to help.

An example with very large WAIC and DIC = NA can be reproduced using the Solitary Tinamou example dataset which comes with PointedSDMs:

library(PointedSDMs)
library(terra)

data("SolitaryTinamou")
projection <- "+proj=longlat +ellps=WGS84"
species <- SolitaryTinamou$datasets
NPP <- scale(terra::rast(system.file('extdata/SolitaryTinamouCovariates.tif',
                                     package = "PointedSDMs"))$NPP)
mesh <- SolitaryTinamou$mesh

model <- intModel(species, spatialCovariates = NPP, Coordinates = c('X', 'Y'),
                  Projection = projection, Mesh = mesh, responsePA = 'Present')
modelRun <- fitISDM(model, options = list(control.inla = list(int.strategy = 'eb'),
                                          safe = TRUE))

The model output from this example is:

Summary of 'bruSDM' object: inlabru version: 2.10.1 INLA version: 23.09.09 Types of data modelled: eBird Present only Parks Present absence Gbif Present only Time used: Pre = 5.28, Running = 277, Post = 0.38, Total = 283 Fixed effects: mean sd 0.025quant 0.5quant 0.975quant mode kld NPP 13.595 0.455 12.704 13.595 14.485 13.595 0 eBird_intercept -16.881 1.910 -20.624 -16.881 -13.138 -16.881 0 Parks_intercept -22.089 2.001 -26.010 -22.089 -18.168 -22.089 0 Gbif_intercept -18.490 1.913 -22.240 -18.490 -14.741 -18.490 0 Random effects: Name Model shared_spatial SPDE2 model Model hyperparameters: mean sd 0.025quant 0.5quant 0.975quant mode Theta1 for shared_spatial -2.36 0.008 -2.37 -2.36 -2.34 -2.37 Theta2 for shared_spatial -5.82 0.013 -5.84 -5.82 -5.79 -5.84 Deviance Information Criterion (DIC) ...............: NA Deviance Information Criterion (DIC, saturated) ....: NA Effective number of parameters .....................: NA Watanabe-Akaike information criterion (WAIC) ...: 4.64e+39 Effective number of parameters .................: 2.32e+39 Marginal log-Likelihood: -790.12 is computed Posterior summaries for the linear predictor and the fitted values are computed (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')

Many thanks,

Peter

[Finn Lindgren]

Hi Peter,

since WAIC is related to leave-one-out cross-validation (see e.g.
https://arxiv.org/pdf/1507.04544.pdf which I think has since been
published), it is not a well-defined concept for point process models,
where the continuous _space_ (together with the individual points) is
the "observation". The likelihood construction in inlabru is not based
on spatial gridding and aggregating counts (which would yield a
poisson count model as the approximate model, for which WAIC would be
well defined) but rather a higher order approximation method. Work is
underway for a more well-defined cross-validation quantity, but as it
is, WAIC is not appropriate to use for general point process models,
whether in inla/inlabru or other software. Only basic aggregation
approximations have well defined WAIC quantities.
(Note: earlier materials for inlabru did compute WAIC, since we hadn't
realised the problem with how WAIC is defined, or rather not defined
at all, in this context!)

For DIC, the situation _should_ be better, but I suspect that one
would need to keep more careful track of normalisation constants.
Until someone has time to do the necessary maths&coding, I prefer
other model assessment methods based on more explicit spatially
resolved prediction evaluation. See for example the point process
residuals vignette at
https://inlabru-org.github.io/inlabru/articles/2d_lgcp_residuals.html

There have also been other recent work on point process forecast
assessments; usually these are based on prediction of aggregated
counts on a specific spatial scale.
The vignette on prediction score computation with inlabru,
https://inlabru-org.github.io/inlabru/articles/prediction_scores.html,
doesn't cover point processes explicitly, but it does cover count
models, which can be applied to point process predictions.

When working with point processes, one has to be mindful that the
observations are of a fundamentally different nature than for
point-referenced observations at given locations, and the usual
turnkey solutions don't always apply!

Finn

> --
> You received this message because you are subscribed to the Google Groups "R-inla discussion group" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to r-inla-discussion...@googlegroups.com.
> To view this discussion on the web, visit https://groups.google.com/d/msgid/r-inla-discussion-group/1d17aa60-658f-4e69-9798-a67758fc8dd6n%40googlegroups.com.



--
Finn Lindgren
email: finn.l...@gmail.com

[Peter Stewart]

Hi Finn,

Thank you for your reply, this is very useful to know! 

I'll look into applying the methods from the residuals vignette instead.

Best wishes,

Peter