INLA update breaks spatiotemporal model

[Eric]

I upgraded from 21.02.23 to 23.09.09 and my model spatiotemporal model gives an error when trying to draw posterior samples and the fitted parameters are off compared to the previous version which looked good.

  GitId: 5829adacd2938bbeabc52513438b1504e8438f65 - Sat Sep 9 15:30:58 2023 +0300

  Error:2 Reason: Matrix is not (numerical) positive definite

  Line:389 Function: GMRFLib_init_problem

The model formula and fit code and summary are at the end. Note that s2 is an spde in space, s is an spde on the same spatial mesh with an additional 1d mesh on time (c.f. Chapter 7.2 of Krainski et al textbook), and t is a descretized time step. The idea is to have a global temporal trend, global spatial mean, and a (separable) spatiotemporal interaction.

There are no errors in the fit but comparing the summary to the successful one on the earlier INLA version there are two differences: First, the intercept 0 +/- 10 has a huge SD (the outcome is standardized) whereas on the earlier INLA version the intercept was 0 +/- 0.03. Second, The earlier version had one quarter the range and stdev for the s2 term.

Note that there is no problem (fit agrees with earlier version and inla.posterior.samples() works) if I drop the spatiotemporal interaction f(s, ...).

Thanks for your help. I plan to try an intermediate INLA version tomorrow.

Eric

----------

mesh <- inla.mesh.2d(
  loc = coo, boundary = bnd,
  max.edge = c(1000, 1000),
  cutoff = 10
)

spde <- inla.spde2.pcmatern(
  mesh = mesh, alpha = 1.5,
  prior.range = c(10, 0.01), # P(range < 10000) = 0.01
  prior.sigma = c(5, 0.01) # P(sigma > 3) = 0.01
)

timesn <- 20
mesh.t = inla.mesh.1d(seq(min(subset_RS21$t)+0.5, max(subset_RS21$t)-0.5, length=timesn))

indexs <- inla.spde.make.index("s",
                               n.spde = spde$n.spde,
                               n.group = timesn
)

indexs2 <- inla.spde.make.index("s2",
                               n.spde = spde$n.spde
)

A <- inla.spde.make.A(mesh = mesh, loc = coo, group = subset_RS21$t, group.mesh=mesh.t)
A2 <- inla.spde.make.A(mesh = mesh, loc = coo)

stk.e <- inla.stack(
  tag = "est",
  data = list(y = subset_RS21$log_n2o %>% scale),
  A = list(1, A, A2),
  effects = list(data.frame(b0 = rep(1, nrow(subset_RS21)),
                            t=as.factor(subset_RS21$t_1h)),
                 s = indexs,
                 s2 = indexs2)
)


pred_RS21 = expand_grid(
  pixels_RS %>% sample_n(5),
  t = seq(min(floor(subset_RS21$t)), ceiling(max(subset_RS21$t)))
)

coop <- pred_RS21[, c("X","Y")] %>% as.matrix
Ap <- inla.spde.make.A(mesh = mesh, loc = coop, group = pred_RS21$t, group.mesh=mesh.t)
Ap2 <- inla.spde.make.A(mesh = mesh, loc = coop)


stk.p <- inla.stack(
  tag = "pred",
  data = list(y = NA),
  A = list(1, Ap, Ap2),
  effects = list(data.frame(b0 = rep(1, nrow(pred_RS21)),
                            t=as.factor(pred_RS21$t)),
                 s = indexs,
                 s2 = indexs2)
)

stk.full <- inla.stack(stk.e, stk.p)

formula <- y ~ 0 + b0 +
  f(s2, model=spde) +
  f(t, model="rw1") +
  f(s, model = spde, group = s.group, control.group = list(model = "rw1"))

res <- inla(formula,
            data = inla.stack.data(stk.full),
            control.predictor = list(
              compute = TRUE,
              A = inla.stack.A(stk.full)
            ),
            control.compute=list(config=TRUE)
)

summary(res)

Fixed effects:
   mean     sd 0.025quant 0.5quant 0.975quant mode kld
b0    0 10.584    -20.755        0     20.755    0   0

Random effects:
  Name    Model
    s2 SPDE2 model
   t RW1 model
   s SPDE2 model

Model hyperparameters:
                                          mean      sd 0.025quant 0.5quant
Precision for the Gaussian observations  22.51   0.811     20.945   22.502
Range for s2                             88.35 124.444      7.895   51.580
Stdev for s2                              1.46   1.743      0.114    0.926
Precision for t                          82.81  16.183     55.508   81.279
Range for s                             154.44  17.037    123.881  153.425
Stdev for s                               2.34   0.190      1.991    2.331
                                        0.975quant    mode
Precision for the Gaussian observations      24.14  22.491
Range for s2                                399.57  20.046
Stdev for s2                                  6.01   0.312
Precision for t                             118.97  78.303
Range for s                                 190.88 151.239
Stdev for s                                   2.74   2.313

[Eric]

PS could you clarify what model this specifies: f(s, model = spde, group = s.group, control.group = list(model = "rw1")) ? I would expect there to be a parameter for the precision of the random walk increments but don't see it in the summary.

[Eric]

One more reply to myself: I upgraded to testing (24.01.18) and the error persists as in 23.09.09.

[INLA help]

This can happen, that numerically it tip over and gives an error.

I would like to rerun the model, can you share the complete R code and data?

—

Haavard Rue

HelpDesk 

help@r-inla. org

--
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/7c543ad5-78dc-4eac-8e4a-bddfce22a8d1n%40googlegroups.com.

[Finn Lindgren]

Hi Eric,

When rw1 is used as a group model, one can interpret it as the main model defining the increment distribution; the model as a whole can only have one common precision scaling parameter, so the precision controlling parameter of the spde model (I.e. the stddev parameter) is the only one in the combined model.

Finn

On 19 Jan 2024, at 05:04, Eric <eric....@gmail.com> wrote:

PS could you clarify what model this specifies: f(s, model = spde, group = s.group, control.group = list(model = "rw1")) ? I would expect there to be a parameter for the precision of the random walk increments but don't see it in the summary.

--

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/81ef963b-3190-4dd9-9cd7-ad08949a4c60n%40googlegroups.com.

[Eric]

Thanks Haarvard for helping me improve the model.

I now need some guidance on constraints for grouped terms. I could not find this in the documentation. A few questions:

1) How is the constraint for a grouped term defined? It seems for example that (with `t` time and `location` a discrete (numeric) location) the following are different:

 f(t, model = "rw1", constr=TRUE,  group = location, control.group = list(model = "iid"))

 f(location, model = "iid", constr=TRUE,  group = t, control.group = list(model = "rw1"))

For my model the first one does the trick.

2) But I want an spde in space, so I am working with this term:

f(s, model = spde, group = s.group, control.group = list(model = "rw1"))

But how can I implement a constraint as if the order were switched?

3) More generally how can I add constraints to the full spde x time? it seems that the extraconstr argument is looking for a constraint with input dimension equal to the number of mesh locations so I am not sure how to add constraints in the time domain?

Thanks,

Eric

[Håvard Rue]

On Mon, 2024-01-22 at 12:41 -0800, Eric wrote:
> Thanks Haarvard for helping me improve the model.
>
> I now need some guidance on constraints for grouped terms. I could not
> find this in the documentation. A few questions:
>
> 1) How is the constraint for a grouped term defined? It seems for
> example that (with `t` time and `location` a discrete (numeric)
> location) the following are different:
>
>  f(t, model = "rw1", constr=TRUE,  group = location, control.group =
> list(model = "iid"))
>  f(location, model = "iid", constr=TRUE,  group = t, control.group =
> list(model = "rw1"))

constraints are replicated over replications and groups. so if

f(s, constr=T)

there is a sum to zero constraint for s, like \sum_i s_i = 0

for

f(s, constr=T, replicate=r)

there are sum to zero constraints for each replications,

\sum_i s_ir = 0, for r=1, ..., nrep

same with group.

>
> For my model the first one does the trick.
>
> 2) But I want an spde in space, so I am working with this term:
>
> f(s, model = spde, group = s.group, control.group = list(model =
> "rw1"))
>
> But how can I implement a constraint as if the order were switched?

you cannot. not all models can be in control.group, and spde is not one
of them.

>
> 3) More generally how can I add constraints to the full spde x time?
> it seems that the extraconstr argument is looking for a constraint
> with input dimension equal to the number of mesh locations so I am not
> sure how to add constraints in the time domain?

you cannot.

the reasoning is that the normalizing constant is now known if you
impose general constraints, whereas if you do that in the replicate
and/or group, setting, it is.

one may argue, that is not that relevant as the overall computational
cost is larger, so this is less relevant, which is true.

its all about implementation. how much do we want to code.

I think if there is a new implementation of inla, then this should be
added, yes.

Håvard Rue
hr...@r-inla.org