Add 1D first and second order random walk x iid models (BYM2 style parameterization)

[szego]

This PR adds two new models: inla.rw1o1diid and inla.rw2o1diid. The names are formatted like:

• "rw" - random walk
• "No" - Nth order ("1o" for first order and "2o" for second order random walk)
• "1d" - one dimensional
• "iid" - iid noise

I've wanted to implement these for a while, ever since I saw this old post on the google group: https://groups.google.com/g/r-inla-discussion-group/c/girUuJhS-Q8/m/3VOkGQ2NAAAJ

The inla.rw2o1diid model is what the poster there was asking for (and what I wanted too, which lead me to that post).

Both models use BYM2-style parameterizations

tau * ( phi * rw + (1 - phi) * iid )

with PC priors.

There is an existing rw2diid model but it doesn't support 1D grids.

inla.rw1o1diid

This is just a wrapper around the existing BYM2 model, as suggested by Elias T. Krainski in the link above. It builds the chain graph internally.

inla.rw2o1diid

This is implemented as an rgeneric and uses inla.pc.dprec for the prior on tau and a new inla.pc.rw2o1diid.phi for the PC prior on phi.

The main difference between this and the RW1 model is that this uses an n-parameter latent field combining the rw and iid components, as opposed to the 2n parameters in the RW1 case. So the user can't get just the RW2 parameter estimates or just the iid parameter estimates like they can for the RW1 model.

One consequence of the n-parameter latent field is that we actually end up doing a BYM2-style mixture of the precision matrices, like Q = tau * ( phi * R + (1 - phi) * I ), where R is the (scaled) RW2 structure matrix. So we're not mixing the variances as in pure BYM2, we're mixing precisions, and tau and phi have slightly different interpretations. The marginal variance of the field is only approximately 1/tau and phi is only approximately the structured-variance fraction.

The new helper inla.pc.rw2o1diid.phi implements a PC prior for phi derived from the Gaussian KLD between the precision-mixture flexible model and the iid base model:

2*KLD(phi) = sum_i 1/(phi*lambda_i + (1-phi)) - n + sum_i log(phi*lambda_i + (1-phi))

where the lambda_i are the eigenvalues of R. The PC distance is sqrt(2*KLD(phi)), the prior on the distance is Exp(lambda), and lambda is calibrated so that P(phi < u) = alpha (default (0.5, 0.5)).

There is also a helper function inla.rw2o1diid.hyperpar to get back the user-scale posteriors for phi and tau - see example below.

Example

library(dplyr)
library(ggplot2)

set.seed(42)

n_timesteps <- 200
excluded_timesteps <- 101:120
n_obs_per_step <- 5
y_mean <- 2 * cumsum(rnorm(n_timesteps, 0, 0.2))
y <- sapply(y_mean, function(x) x + rnorm(n_obs_per_step, 0, 1))
d <-
  data.frame(
    time = rep(seq_len(n_timesteps), each = n_obs_per_step),
    y = as.numeric(y)
  ) %>%
  mutate(y = if_else(time %in% excluded_timesteps, NA_real_, y))

inla_fit <- inla(
  y ~ f(time, model = inla.rw2o1diid(n_timesteps)),
  data = d
)

d_pred <-
  inla_fit$summary.fitted.values %>%
  filter(row_number() %% n_obs_per_step == 1) %>%
  select(
    mean,
    lower = `0.025quant`,
    upper = `0.975quant`
  ) %>%
  mutate(time = row_number())

ggplot() +
  geom_point(
    data = d,
    aes(x = time, y = y)
  ) +
  geom_ribbon(
    data = d_pred,
    aes(x = time, ymin = lower, ymax = upper),
    fill = "blue",
    alpha = 0.3
  ) +
  geom_line(
    data = d_pred,
    aes(x = time, y = mean),
    inherit.aes = FALSE,
    col = "blue"
  ) +
  theme_minimal()

hp <- inla.rw2o1diid.hyperpar(inla_fit, "time")
inla.zmarginal(hp$tau)
inla.zmarginal(hp$phi)

# > inla.zmarginal(hp$tau)
# Mean            0.0646981
# Stdev           0.0148991
# Quantile  0.025 0.0402192
# Quantile  0.25  0.0540026
# Quantile  0.5   0.0630489
# Quantile  0.75  0.0735629
# Quantile  0.975 0.0984835
#
# > inla.zmarginal(hp$phi)
# Mean            0.0352956
# Stdev           0.0178236
# Quantile  0.025 0.011846
# Quantile  0.25  0.0225539
# Quantile  0.5   0.0315153
# Quantile  0.75  0.0438079
# Quantile  0.975 0.0804446

[szego]

Looks like the common practice in this repo is to squash PRs into a single commit. Happy to do that, if so.

[hrue]

Thank you, looking good!

Just one comment, I think the RW2 case can be implemented as you do for the RW1, as inla.rw(order=<>,n=<>, scale.model=T), do the right thing, and passing eigenvalues/vectors to the prior calculation (then set scale.model=FALSE, as your RW model is already scaled), and then there would be no need for 'rgeneric' (which is slow and break parallel....)

what u think?

[hrue]

aaah, now I see why this is needed for rw2... my bad.

[hrue]

Let me check the internal C-code is something can be done

[szego]

Thank you! I'll revert this to draft for now.

[hrue]

I think an alternative option is to add in the C-code the RW2+IID construction, and then the RW1+IID, will follow as a special case. The little tricky thing is the prior computations, but this is fine as we can use eigenvalues/vectors. I'll put this on the list, and try to have it done soon

…

On Sun, 2026-05-03 at 09:25 -0700, Antonio R. Vargas wrote: szego left a comment (hrue/r-inla#122) Thank you! I'll revert this to draft for now. — Reply to this email directly, view it on GitHub, or unsubscribe. You are receiving this because you commented.Message ID: ***@***.***>

-- Håvard Rue Professor of Statistics Chair of the Statistics Program CEMSE Division King Abdullah University of Science and Technology Thuwal 23955-6900 Kingdom of Saudi Arabia ***@***.*** Office: +966 (0)12 808 0640   Mobile: +966 (0)54 470 0421 Research group: bayescomp.kaust.edu.sa   R-INLA project: www.r-inla.org Zoom: kaust.zoom.us/my/haavard.rue

-- This message and its contents, including attachments are intended solely for the original recipient. If you are not the intended recipient or have received this message in error, please notify me immediately and delete this message from your computer system. Any unauthorized use or distribution is prohibited. Please consider the environment before printing this email.

[szego]

I've implemented the PC prior for this precision mixture model. The code now uses that instead of using the not-exactly-appropriate inla.pc.bym.phi() prior.