Problem when including spatial component + 1D SPDE

[Silas C Principe]

Dear all, hi. 

I'm fitting LGCP models through inlabru/INLA and everything was working relatively fine. However, I have recently updated my dataset, which nows include 3x more points than the previous one (now it is >1000), and the model stopped working when I include the spatial component (SPDE) [it crashes or, if changing the pc priors, it works but produces absurd results].

It seems to be an identifiability issue, as Finn previous highlighted to me in another case, but including an "extraconstr.int" argument is not helping. Of course, it may be the case I'm not setting this argument correctly, as I never had to use that before.

An important note is that if I exclude two covariates that are considered to better explain the spatial distribution of the species, then the model runs fine and I'm able to recover the lambda estimates. One of those variables have a non-linear relationship with the occurence of points, and is being modelled used a 1D SPDE (instead of a rw2 model).

I would like to know if there is any sensible way of solving this problem so that both the spatial component and the covariates could be included. I appreciate any tips you can give.

I have attached both the code and data needed for it.

A further note: I'm receiving a lot of "GMRFLib_2order_approx: rescue NAN/INF values in logl for idx=" errors when including the spatial component.

All the best,

Silas

[Silas C Principe]

Dear all, hi!  Just to add some updates to my question:

I tried to make some changes to the mesh to have less points and a simpler inner boundary. Doing that I was abble to at least reach a model which run with just 1 try and that also returned less absurd results. However it is still not working as expected.

Something that I realized is that the range for the 2D SPDE model (spatial) and the standard deviation for the 1D spde model (sst, a variable with non-linear relation with the response) are giving almost the same values (see summary at the end). Is this a sign of non identifiability or of a mispecification in the model?

Thanks again for all the help!

--- Summary of the model:

inlabru version: 2.5.0
INLA version: 22.03.28
Components:
  Intercept: Model types main='linear', group='exchangeable', replicate='iid'
  sal: Model types main='linear', group='exchangeable', replicate='iid'
  sst: Model types main='spde', group='exchangeable', replicate='iid'
  dist: Model types main='linear', group='exchangeable', replicate='iid'
  spatial: Model types main='spde', group='exchangeable', replicate='iid'
Likelihoods:
  Family: 'cp'
    Data class: 'SpatialPointsDataFrame'
    Predictor: coordinates ~ .
Time used:
    Pre = 16.2, Running = 568, Post = 1.52, Total = 585
Fixed effects:
            mean    sd 0.025quant 0.5quant 0.975quant mode   kld
Intercept -2.939 0.058     -3.053   -2.939     -2.824   NA 0.130
sal        0.099 0.042      0.016    0.099      0.182   NA 0.001
dist      -0.198 0.031     -0.259   -0.198     -0.137   NA 0.006

Random effects:
  Name      Model
    sst SPDE2 model
   spatial SPDE2 model

Model hyperparameters:
                    mean    sd 0.025quant 0.5quant 0.975quant mode
Stdev for sst     229.90 0.365     229.20   229.89     230.63   NA
Range for spatial 229.27 0.344     228.60   229.26     229.95   NA
Stdev for spatial   6.78 0.011       6.76     6.78       6.80   NA

Deviance Information Criterion (DIC) ...............: -39096.56
Deviance Information Criterion (DIC, saturated) ....: -63646.33
Effective number of parameters .....................: -39096.56

Watanabe-Akaike information criterion (WAIC) ...: 14608.38
Effective number of parameters .................: 4677.48

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

[Finn Lindgren]

That looks a bit strange. Certainly having a std.dev. of >200 for an
spde component (sst) as part of a log-Gaussian inhomogeneous Poisson
point process is unreasonable; it would indicate an extraordinarily
wide intensity range:

> exp(c(-2,2)*230)
[1] 1.67702e-200 5.962957e+199

But your results don't show a range parameter for sst. Perhaps there's
a problem with the model definition for that component?

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/76063299-8961-4721-afdc-85c4937ac050n%40googlegroups.com.



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

[Silas Candido Principe de Souza]

Hi Finn! Thanks for helping with this.

In this particular case I fixed the range parameter to a certain large value, following advice from another INLA-group thread about 1D SPDE (https://groups.google.com/g/r-inla-discussion-group/c/YPf4CZKexB4/m/09m9Dj1nCQAJ). The problem in this thread was pretty much the same as I was having... The idea of using a 1D spde instead of a rw2 model was to avoid the binning of the covariate.

Anyway, I tried to run a model with non-fixed range and the problems persist (i.e. awkward range and SD estimates for both the 1D and 2D spde).

The model without the 2D SPDE (spatial comp.) runs fine (i.e., the 1D SPDE component works). In the opposite direction it is the same (i.e., I can run a model with the spatial component, but leaving the 1D SPDE out).

Any light on what can be happening?

Thanks again,

Silas

[Finn Lindgren]

Yes, but fixing the range parameter shouldn't cause the stdev parameter to become that large.

Finn

[Silas Candido Principe de Souza]

Hi Finn!

After the last message I tried some different things and I was able to reach a model that is working more properly. The spatial component effect and SD is still a little too high [see figure attached] compared to the models I was running before getting a dataset with more points, but I'm not sure if this is or not a problem.

Basically I did 2 things to get better results:

1. Increased the number of mesh nodes (decreasing cutoff and max. edge len.) from ~1700 to ~7300. After that, the model was already returning more acceptable estimates for the hyperpriors.

2. Reduced the value of sigma in the pcmatern from 0.5 to 0.1 (and 0.01 prob.).

Thus, I would like to know if (1) there is any general guidance on how many nodes we need on a mesh on a LGCP setting (when the nodes are integration points). The area is quite big in my case (from ~45°S up to 45°N, almost the whole of America - see attached figure) and I'm using a projection in km. Maybe I need to include even more nodes?

And (2) how should I use the 'extraconstr' or 'extraconstr.int' when using inlabru. I think I still need to constrain it to be able to make the model with the 1D-SPDE component work more properly... Should the 'Aconstr' argument contain the data only for the integration points and the data points or for the whole mesh nodes?

Thanks again for all the help! I hope this thread will also help others working with coastal areas in the future!

Silas