GMRFLib_2order_approx: rescue NAN/INF values in logl for idx

[Andrew]

Hello all,

I am trying to run a simple LGCP model with inlabru containing one "SpatialGridDataFrame" covariate (NDVI) but when i attempt to run the model, i am faced with thousands of error codes (GMRFLib_2order_approx: rescue NAN/INF values in logl for idx) before INLA eventually crashes.

For context, the same model code runs fine with other covariates which are formatted in the same way (e.g. same resolution, origin, extent, CRS). Another note is that i have only included NDVI values ranging between 0 and 1 (no negative values). I will include the code below, does anyone have any ideas to fix this, i'm completely stuck?

Thanks!

# Set experimental inlabru mode
>
>
>
> bru_options_set(inla.mode = "experimental")
>
>
>
> # Reload data
>
>
>
> campbells <- readRDS("campbells.rds")
>
>
>
> # Assign NDVI covariate
>
>
>
> NDVI <- campbells$covariates$NDVI
>
>
>
> # Hazard rate detection function with shape parameter of 7
>
>
>
> hr <- function(distance, lsig) {
>
>   1 - exp(-(distance / exp(lsig))^-7)
>
> }
>
>
>
> # SPDE with Matern covariance structure and PC priors
>
>
>
> matern <- inla.spde2.pcmatern(campbells$mesh_fine,
>
>                               prior.sigma = c(5, 0.01),
>
>                               prior.range = c(10, 0.9))
>
>
>
> # Define and run test model
>
>
>
> components_test <- ~ SPDE(main = coordinates, model = matern) +
>
>   NDVI(NDVI, model = "linear") +
>
>   lsig(1) +
>
>   Intercept(1)
>
>
>
> formula_test <- coordinates + distance ~ SPDE + log(hr(distance, lsig)) +
>
>   NDVI +
>
>   Intercept +
>
>   log(2)
>
>
>
> model_test <- lgcp(
>
>   components = components_test,
>
>   campbells$points,
>
>   samplers = campbells$samplers,
>
>   domain = list(
>
>     coordinates = campbells$mesh_fine,
>
>     distance = INLA::inla.mesh.1d(seq(0, 0.06, length.out = 30))
>
>   ),
>
>   formula = formula_test,
>
>   options = list(bru_initial = list(lsig = -3, Intercept = 1), bru_verbose = 4, bru_method=list(rel_tol = 0.15), control.compute = list(dic = TRUE, cpo = TRUE, waic=TRUE))
>
> )
>
>
>
> summary(model_test)

[Andrew]

For more info, i allowed INLA to run through the errors and eventually reach the maximum of ten steps to converge, and this is the summary of the model it produced:

inlabru version: 2.7.0
INLA version: 22.07.23
Components:
SPDE: main = spde(coordinates)
NDVI: main = linear(covariates$NDVI)
lsig: main = linear(1)
Intercept: main = linear(1)
Likelihoods:
  Family: 'cp'
    Data class: 'SpatialPointsDataFrame'
    Predictor:

        coordinates + distance ~ SPDE + log(hr(distance, lsig)) + NDVI +
            Intercept + log(2)

Time used:
    Pre = 3.26, Running = 50.8, Post = 0.361, Total = 54.4
Fixed effects:
              mean     sd 0.025quant 0.5quant 0.975quant     mode kld
NDVI      5215.773 23.004   5170.658 5215.773   5260.888 5215.773   0
lsig         0.000 31.623    -62.017    0.000     62.017    0.000   0
Intercept -392.960 30.310   -452.403 -392.960   -333.517 -392.960   0

Random effects:
  Name   Model
    SPDE SPDE2 model

Model hyperparameters:
                 mean    sd 0.025quant 0.5quant 0.975quant   mode
Range for SPDE  14.94 0.004      14.93    14.94      14.95  14.94
Stdev for SPDE 241.21 0.143     241.08   241.20     241.42 241.15

Deviance Information Criterion (DIC) ...............: 125761.92
Deviance Information Criterion (DIC, saturated) ....: 25508.33
Effective number of parameters .....................: 124599.21

Watanabe-Akaike information criterion (WAIC) ...: -Inf
Effective number of parameters .................: 2702.78

Marginal log-Likelihood:  29822.24
CPO, PIT 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)')

[Finn Lindgren]

Hi, I just now noticed the ^-7 in your hr() definition; this is likely what causes the lsig parameter to be essentially ignored, and you’re getting the prior for that out as posterior.  Even if there might be a theoretically different posterior for losing, you’re likely running into numerical issues when evaluating it, as the log(1-exp(…)) construction is liable to numerical cancellation/truncation errors.  Sometimes one can use log1p to get more numerically stable evaluation of it, and I would strongly recommend looking into how to define a numerically stable implementation of log_hr directly, instead of relying on log(hr()), as relying solely on the numerical evaluation is likely to cause issues even with a more friendly exponent instead of -7!

Finn

On 23 Mar 2023, at 21:04, Andrew <a.houl...@hotmail.co.uk> wrote:

For more info, i allowed INLA to run through the errors and eventually reach the maximum of ten steps to converge, and this is the summary of the model it produced:

--
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/a7a54231-6aba-4014-9d38-3fdffec14aadn%40googlegroups.com.

[Andrew]

Hi Finn,

Thanks for the quick response. I have just tried running again but with the following definition for a numerically stable half-normal function:

log_hn <- function(distance, lsig) {
  (-0.5 * (distance / exp(lsig))^2)
}

However, i am unfortunately still receiving the same error messages.

Andrew