Advice on how to proceed when SPDE+covariate gives unrealistic values

[Rekha Mohan]

Hello,

I'm pretty new to inlabru, but I'm using it to analyse some spatial point process models using LGCP. I've been following the tutorials and vignettes online, but I've run into some issues that I'm having trouble solving. Would really appreciate any help because I'm pretty lost on how to proceed.

Some background on what I'm doing: I'm analysing lion presence in Kruger National Park. The SPDE model comes out fine with my priors( prior.sigma = c(1, 0.1), prior.range = c(10, 0.5)) , but the trouble is when I'm including covariates. The main covariates that I'm including are:

- NDVI (continuous), scaled

- Vegetation type (categorical)

When I fit the model with NDVI without SPDE, the results appear fine. But when I include SPDE, the values I get are astronomical:

> summary(fit2)
inlabru version: 2.12.0.9024
INLA version: 25.04.16
Components:
ndvi: main = linear(f.ndvi(.data.)), group = exchangeable(1L), replicate = iid(1L), NULL
mySmooth: main = spde(geometry), group = exchangeable(1L), replicate = iid(1L), NULL
Intercept: main = linear(1), group = exchangeable(1L), replicate = iid(1L), NULL
Observation models:
  Family: 'cp'
    Tag: <No tag>
    Data class: 'sf', 'tbl_df', 'tbl', 'data.frame'
    Response class: 'numeric'
    Predictor: geometry ~ .
    Additive/Linear: TRUE/TRUE
    Used components: effects[ndvi, mySmooth, Intercept], latent[]
Time used:
    Pre = 1.33, Running = 137, Post = 0.881, Total = 139
Fixed effects:
              mean    sd 0.025quant 0.5quant 0.975quant
ndvi      1393.065 1.634   1389.862 1393.065   1396.269
Intercept  186.449 1.227    184.045  186.449    188.853
              mode kld
ndvi      1393.065   0
Intercept  186.449   0

Random effects:
  Name   Model
    mySmooth SPDE2 model

Model hyperparameters:
                     mean    sd 0.025quant 0.5quant
Range for mySmooth   6.40 0.001       6.40     6.40
Stdev for mySmooth 890.50 0.150     890.21   890.50
                   0.975quant   mode
Range for mySmooth       6.40   6.40
Stdev for mySmooth     890.81 890.50

Deviance Information Criterion (DIC) ...............: 1374418.14
Deviance Information Criterion (DIC, saturated) ....: 1374416.86
Effective number of parameters .....................: 1369971.34

Marginal log-Likelihood:  310300.30
 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)')

In addition, I get a whole bunch of warnings that go on for a while:

 *** WARNING *** GMRFLib_2order_approx: rescue NAN/INF values in logl for idx=0
 *** WARNING *** GMRFLib_2order_approx: reset counter for 2 NAN/INF values in logl
 *** WARNING *** GMRFLib_2order_approx: rescue NAN/INF values in logl for idx=0

    Model        DIC    Delta.DIC
1    fit1(NDVI)  -78314.48       0.0000
2 fitnull (SPDE)  -78043.73     270.7531
3    fit3 (VEG+SPDE)  -78036.08     278.4001
4    fit5  (VEG+NDVI)-74867.66    3446.8231
5    fit2 (NDVI+SPDE)1374418.14 1452732.6209

I read a little bit about spatial confounding, and was wondering if that's the case here? Because NDVI and SPDE are trying the explain the same thing. It seems that the model with NDVI is the best model from the DIC, but that defeats the purpose of using a spatial model, right? I'm not really sure how to proceed.

Would appreciate any advice on this! Happy to share any code snippets if needed as well.

[Finn Lindgren]

Hi,

I’ve never seen as strange estimates for this type of model as those, so there’s definitely a problem of some kind.

How many points are in your data? How did you construct the mesh? What CRS are you modelling in? (It will normally manage to convert the input data to the CRS you specified to fm_mesh_2d, but if you convert it yourself you have more safety) normally, it should be some suitable local equal-area projection, converted to km units; with the latest cran fmesher version, you can convert a CRS to km with fm_crs(mycrs, units=“km”). Also, what resolution are the covariates in, compared with the mesh edges and estimated correlation range.

The first thing I would do would be to plot the posterior prediction of the SPDE component to see if it is consistently far from zero, or if it has lots of localised “spikes”.

I would be very wary of using the DIC values, as they aren’t necessarily correct for the LGCP likelihood construction (the _is_ a version of the DIC definition that is valid for LGCP likelihoods, but it is not entirely clear whether the computation used by inla is that version as we use a trick of the interface to specify the likelihood; I hope to clear this up in the future).

Finn

On 17 Jul 2025, at 09:45, Rekha Mohan <reky...@gmail.com> wrote:

Hello,

--
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, visit https://groups.google.com/d/msgid/r-inla-discussion-group/073bd0d8-81f8-419c-9277-1ec2611893c2n%40googlegroups.com.

[Finn Lindgren]

Additional:

Just in case the problem is the same as some others have had where their domain boundary polygon ended up with a different sf geometry column name than expected, please upgrade to the latest cran version of fmesher, 0.5.0, which has added code to detect that problem (which would cause the LGCP to be defined on the entire mesh instead of the study domain, which can cause all sorts of problems…)

Finn

On 17 Jul 2025, at 10:22, Finn Lindgren <finn.l...@gmail.com> wrote:

Hi,

[Rekha Mohan]

Hi Finn,

Thank you for your reply! My fmesher version is  ‘0.5.0.9000’ so that shouldn't be the problem. I've got 5332 datapoints. Here's how I created the mesh:

> res(gcov$ndvi)

[1] 0.08864447 0.08864447

> st_crs(cleanedlion2024km)
Coordinate Reference System:
  User input: +proj=utm +zone=35 +south +units=km +datum=WGS84
  wkt:

mesh <- fm_mesh_2d(
  loc = fm_hexagon_lattice(cleanedlion2024km),
  boundary = roadskm,
  max.edge = c(10,20),
  offset=c(10,20),
  cutoff=1
)

fm_crs(mesh)<-st_crs(cleanedlion2024km)

> mesh$n
[1] 2753

I predicted the map above (fitnull model):

cmpnull<- geometry ~   mySmooth(geometry, model = matern) + 1

fitnull<- lgcp(cmpnull, cleanedlion2024km, samplers = roadlineskm, domain = list(geometry = mesh),weights=weight_scaled, options = list(
  control.compute = list(dic = TRUE, config = TRUE),
  control.predictor = list(compute = TRUE)))

intnull <- predict(fitnull, pred.df, ~ exp(mySmooth + Intercept), n.samples = 1000)
ggplot() +
  gg(intnull, geom = "tile") +
  geom_sf(data=kruger, alpha = 0, lwd = 1)

And the yellow-green map of fit2 model(NDVI+SPDE), just the SPDE component:

cmp2<- geometry ~ ndvi(f.ndvi(.data.), model = "linear") +
  mySmooth(geometry, model=matern) + 1

fit2<- lgcp(cmp2, cleanedlion2024km, samplers = roadlineskm, domain = list(geometry = mesh),weights=weight_scaled, options = list(
  control.compute = list(dic = TRUE, config = TRUE),
  control.predictor = list(compute = TRUE)))

ndvi.lp <- predict(
  fit2,
  pred.df,
  ~ list(
    smooth_ndvi= mySmooth + ndvi,
    ndvi = ndvi,
    smooth=mySmooth
  )
)

ggplot() +
  gg(ndvi.lp$smooth,geom="tile") +
  csc +
  theme(legend.position = "bottom") +
  geom_sf(data=kruger, alpha = 0) +
  ggtitle("SPDE")

Thank you again!

[Finn Lindgren]

What are you trying to do with the “weights” argument? That has a special meaning for the “cp” model. If you wanted to specify line transect weights, that needs to be done as a “weight” column in the samplers data.

Also note that your estimated range is smaller than your max.edge value, so either the mesh needs to have smaller triangles, or you need to set a more realistic priors (note that the prior needs to be set in context of what else is in the model; one cannot expect parameters to mean the same thing when adding covariates, for example; the SPDE needs to capture that which is not captured by the covariates, so the effective correlation range needed can be very different…

Finn

On 17 Jul 2025, at 15:01, Rekha Mohan <reky...@gmail.com> wrote:



To view this discussion, visit https://groups.google.com/d/msgid/r-inla-discussion-group/0c745a08-d298-4755-8bae-ddd5749e1190n%40googlegroups.com.

<Rplot01.png>

<Rplot.png>

[Rekha Mohan]

Hi,

Ah I see - the weights are referring to how much weight each sighting should have. Because it's citizen science data, I wanted to account for this bias. How should I go about incorporating that into model?

I'll tweak the mesh, and give it a go. Thank you Finn!

[Finn Lindgren]

If there is a "weight" column in the samplers argument, with the sampling weight for each line/polygon/transect, it will be incorporated in the integration automatically, i.e. treated as an "effort" multiplier to the basic sampling geometry; with no weight info, it uses "1 unit per distance/area" (depending on lines/polygons).

The weights argument you used is instead applied to the log-intensity likelihood contribution for each observed point.

It might help if you could write down what you want your likelihood model to be, precisely.

The "cp" model has likelihood

-\sum_j \int_{sampler[j]_shape} sampler[j]_weight \exp(\eta(s)) ds + \sum_i obs_weight_i \eta(y_i)

where the sampler_weight values come from the samplers weight column, and obs_weight comes from the weights argument you used.

(obs_weight is a new feature meant for things like location-uncertainty modelling, but it is possible that other interpretations/constructions lead to the same mathematical expression...)

I'm currently involved in a research project on combining other data with citizen science data; there is no fully formed "this is how one _should_ do it" answer, but several approaches, and one has to be a bit careful also about how the linear predictor is constructed when combining data to ensure that it has the same meaning for different observation types. We've seen quite a few wrong approaches in the literature, so be careful and think through what your model is! "Here be dragons..."

Finn

To view this discussion, visit https://groups.google.com/d/msgid/r-inla-discussion-group/bf5ff772-d104-4edc-b3e6-6aaaa45b1c4bn%40googlegroups.com.

--

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

[Rekha Mohan]

Thank you! I think in terms of weights, I actually might stick to the cp model, because I want to reduce the weight of the sightings themselves rather than the sampling weights (I can't quantify that). But I'll experiment with both, thanks for explaining!

For the model now:

I  decreased the resolution of my covariates to 1km by 1km and decreased the mesh size to 5km. 5km was the smallest I could manage without R crashing. Do you think that's okay? My priors are now:

matern <- inla.spde2.pcmatern(
  mesh,
  prior.sigma = c(0.3, 0.01), 
  prior.range = c(50, 0.5)
)

Which runs for just SPDE, but when I include the NDVI covariate, I get many warnings like this. Should I continue tweaking the priors so that it'll fit?

 *** WARNING *** GMRFLib_2order_approx: reset counter for 12 NAN/INF values in logl
 *** WARNING *** GMRFLib_2order_approx: rescue NAN/INF values in logl for idx=130
 *** WARNING *** GMRFLib_2order_approx: reset counter for 37 NAN/INF values in logl
 *** WARNING *** GMRFLib_2order_approx: reset counter for 37 NAN/INF values in logl

***[1] warning *** delta[0] is NAN, 'vb.correction' is aborted
*** Please (re-)consider your model, priors, confounding, etc.


***[5] warning *** delta[0] is NAN, 'vb.correction' is aborted
*** Please (re-)consider your model, priors, confounding, etc.

> summary(fit2)
inlabru version: 2.12.0.9024
INLA version: 25.04.16
Components:
ndvi: main = linear(f.ndvi(.data.)), group = exchangeable(1L), replicate = iid(1L), NULL
mySmooth: main = spde(geometry), group = exchangeable(1L), replicate = iid(1L), NULL
Intercept: main = linear(1), group = exchangeable(1L), replicate = iid(1L), NULL
Observation models:
  Family: 'cp'
    Tag: <No tag>
    Data class: 'sf', 'tbl_df', 'tbl', 'data.frame'
    Response class: 'numeric'
    Predictor: geometry ~ .
    Additive/Linear: TRUE/TRUE
    Used components: effects[ndvi, mySmooth, Intercept], latent[]
Time used:

    Pre = 1.34, Running = 1022, Post = 4.38, Total = 1027

Fixed effects:
               mean    sd 0.025quant  0.5quant 0.975quant

ndvi      -7988.618 6.655  -8001.661 -7988.618  -7975.575
Intercept -4362.284 7.402  -4376.791 -4362.284  -4347.777
               mode kld
ndvi      -7988.618   0
Intercept -4362.284   0

Random effects:
  Name   Model
    mySmooth SPDE2 model

Model hyperparameters:
                    mean    sd 0.025quant 0.5quant

Range for mySmooth  7.13 0.000       7.13     7.13
Stdev for mySmooth 94.56 0.011      94.55    94.56
                   0.975quant  mode
Range for mySmooth       7.13  7.13
Stdev for mySmooth      94.59 94.55

Deviance Information Criterion (DIC) ...............: 3040078.38
Deviance Information Criterion (DIC, saturated) ....: 3040075.53
Effective number of parameters .....................: 3033346.91

Marginal log-Likelihood:  730285.71

 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]

Did you upgrade fmesher? I still suspect there’s an issue with you data.

With covariate resolution 1km, the integration scheme needs to be finer than 5km, but the mesh itself doesn’t need to be as fine. Use geometry=fm_subdivide(mesh, 4) to get approx 1km integration scheme resolution from a 5km mesh.

To be clear, both types of weights are part of the “cp” model.

Observation weights between 0 and 1 would mean “I observed a fraction of a point” essentially. Sampler weights relate to the effective amount of search effort. Not knowing what those values should be doesn’t mean one can just ignore them though. Not setting them is the same as actively setting them to 1.

Finn

On 17 Jul 2025, at 23:11, Rekha Mohan <reky...@gmail.com> wrote:



<Rplot.png>

To view this discussion, visit https://groups.google.com/d/msgid/r-inla-discussion-group/b6643e1b-cac2-4012-89f1-22439d4af7d5n%40googlegroups.com.

<Rplot01.png>

<Rplot.png>

[Rekha Mohan]

Hi,

Thank you, that worked! I'm still getting these warnings, but perhaps I just need to have better priors?:

***[1] warning *** iterative process seems to diverge, 'vb.correction' is aborted

*** Please (re-)consider your model, priors, confounding, etc.

***[2] warning *** iterative process seems to diverge, 'vb.correction' is aborted

*** Please (re-)consider your model, priors, confounding, etc.

***[4] warning *** iterative process seems to diverge, 'vb.correction' is aborted

*** Please (re-)consider your model, priors, confounding, etc.

The summary output is more reasonable, though the predictions go up to 8000, which is still too high. I've also uninstalled fmesher and reinstalled it (the developer version this time).

> packageVersion("fmesher")
[1] ‘0.5.0.9003’

> summary(fit2)
inlabru version: 2.12.0.9024
INLA version: 25.04.16
Components:
ndvi: main = linear(f.ndvi(.data.)), group = exchangeable(1L), replicate = iid(1L), NULL
mySmooth: main = spde(geometry), group = exchangeable(1L), replicate = iid(1L), NULL
Intercept: main = linear(1), group = exchangeable(1L), replicate = iid(1L), NULL
Observation models:
  Family: 'cp'
    Tag: <No tag>
    Data class: 'sf', 'tbl_df', 'tbl', 'data.frame'
    Response class: 'numeric'
    Predictor: geometry ~ .
    Additive/Linear: TRUE/TRUE
    Used components: effects[ndvi, mySmooth, Intercept], latent[]
Time used:

    Pre = 1.53, Running = 17.4, Post = 1.45, Total = 20.4

Fixed effects:
            mean    sd 0.025quant 0.5quant 0.975quant

ndvi      -0.075 0.043     -0.159   -0.075      0.010
Intercept -4.782 0.416     -5.639   -4.768     -4.004
            mode kld
ndvi      -0.075   0
Intercept -4.768   0

Random effects:
  Name   Model
    mySmooth SPDE2 model

Model hyperparameters:
                    mean    sd 0.025quant 0.5quant

Range for mySmooth 11.96 1.059      10.04    11.91
Stdev for mySmooth  3.82 0.218       3.42     3.81
                   0.975quant  mode
Range for mySmooth      14.21 11.78
Stdev for mySmooth       4.27  3.79

Deviance Information Criterion (DIC) ...............: -56937.93
Deviance Information Criterion (DIC, saturated) ....: -56978.89
Effective number of parameters .....................: -68908.00

Marginal log-Likelihood:  -41699.71

 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)')

Ah thanks for clarifying about the weights, I get it now. I'll try to include that column in the sampler then, hopefully it helps!

Thank you again, really really appreciate this!

[Finn Lindgren]

Very strange to get such large values for the SPDE; can you upgrade INLA as well, just in case?

Finn

On 18 Jul 2025, at 01:54, Rekha Mohan <reky...@gmail.com> wrote:

Hi,

To view this discussion, visit https://groups.google.com/d/msgid/r-inla-discussion-group/254d5051-f528-4737-95b0-c36f0d336e28n%40googlegroups.com.

<Rplot.png>

[Finn Lindgren]

There does seem to be some structure in the estimate though, as if it’s counteracting some feature of the NDVI covariate perhaps. But exp(5000) would lead to overflow in the calculations so there must be something in the model that prevents that. Have you checked that there aren’t errors in the NDVI data, such as missing values encoded as -999 or similar? That’s the only real explanation I can think of that could cause this behaviour.

Finn

On 18 Jul 2025, at 07:44, Finn Lindgren <finn.l...@gmail.com> wrote:

Very strange to get such large values for the SPDE; can you upgrade INLA as well, just in case?

[Rekha Mohan]

Hi Finn,

Thank you for all your help so far! So I redownloaded my covariates and reran the model. I kept getting the warnings regarding the wrong choice of priors so I iteratively changed the priors and managed to get to a point where there were no warnings (is that the right way to go about choosing priors?). I settled on the following. I thought they were pretty tight, but it seemed to work on the model.

matern <- inla.spde2.pcmatern(
  mesh,
  prior.sigma = c(0.03, 0.01), 
  prior.range = c(40, 0.5) 
)

Using DIC, it seemed that the best models were the ones without SPDE, but I understand that I should be wary about using DIC, so I think I'll include the models with SPDE and do posterior predictive checks. How do we choose between models with different covariates without looking at DIC though?

> summary(fit2)
inlabru version: 2.12.0.9024
INLA version: 25.04.16
Components:
ndvi: main = linear(f.ndvi(.data.)), group = exchangeable(1L), replicate = iid(1L), NULL
mySmooth: main = spde(geometry), group = exchangeable(1L), replicate = iid(1L), NULL
Intercept: main = linear(1), group = exchangeable(1L), replicate = iid(1L), NULL
Observation models:
  Family: 'cp'
    Tag: <No tag>
    Data class: 'sf', 'tbl_df', 'tbl', 'data.frame'
    Response class: 'numeric'
    Predictor: geometry ~ .
    Additive/Linear: TRUE/TRUE
    Used components: effects[ndvi, mySmooth, Intercept], latent[]
Time used:

    Pre = 2.31, Running = 43.2, Post = 3.01, Total = 48.5

Fixed effects:
            mean    sd 0.025quant 0.5quant 0.975quant

ndvi      -0.020 0.043     -0.104   -0.020      0.064
Intercept -3.999 0.206     -4.423   -3.994     -3.610
            mode kld
ndvi      -0.020   0
Intercept -3.994   0

Random effects:
  Name   Model
    mySmooth SPDE2 model

Model hyperparameters:
                   mean    sd 0.025quant 0.5quant

Range for mySmooth 8.93 0.437       8.11     8.92
Stdev for mySmooth 2.03 0.061       1.92     2.03
                   0.975quant mode
Range for mySmooth       9.83 8.90
Stdev for mySmooth       2.15 2.03

Deviance Information Criterion (DIC) ...............: -52267.90
Deviance Information Criterion (DIC, saturated) ....: -52308.32
Effective number of parameters .....................: -65383.14

Marginal log-Likelihood:  -40697.53

 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)')

The predictions seem reasonable now:

> range(ndvi.lp$ndvi_smooth$mean)
[1]  0.005736814 42.034006825

Just wanted to say thank you again for all your help!

[Finn Lindgren]

I’m confused.

Are you still getting the nonsensical large values for the SPDE model?

Did you check the actual values of the covariate?

Without checking those things, I would not rely on the results _at all_, and definitely not trust any summary metrics such as DIC.

Finn

On 21 Jul 2025, at 15:38, Rekha Mohan <reky...@gmail.com> wrote:

Hi Finn,

To view this discussion, visit https://groups.google.com/d/msgid/r-inla-discussion-group/a03e0946-8c7e-4519-b10a-5f75ab181469n%40googlegroups.com.

[Rekha Mohan]

Hello,

Nope, my values for the SPDE are very normal now. What I changed was instead of 

exp(mySmooth + 1), I put exp(mySmooth + Intercept) and it seemed to work! Here are my predicted values now: