spatio-temporal model with replicate and group

25 views
Skip to first unread message

Joe Lewis

unread,
Oct 2, 2025, 11:46:30 AM (yesterday) Oct 2
to R-inla discussion group
Dear all,

I'm trying to fit a multinomial model (using the Poisson trick) that varies spatio-temporally. With this model, I am trying to predict the probability of a data point being a particular 'authority' (think region/country) given its location and covariates. From this, I will e.g. predict the probability of a particular location being an authority given all other authorities (i.e. normalise across the predicted rates to get a probability).

For each data point I've replicated the rows across the number of alternative Authorities that the data point could be attributed to. A new column (obs) is created that is either 1 or 0, with 1 denoting which of the replicated rows is the observed authority. All other replicated rows within that 'strata' are given a value of 0.

Each observed data point has an associated count - the greater the count, the more confident that this location is assigned to a particular authority. So ideally, I want to weight the likelihood to take this into account. Currently, I have this count logged and included as an offset but I'm not sure if this is the correct approach?

Spatially, the data points denoted to a specific 'authority' are concentrated within a particular part of Britain. For example, one authority is concentrated in the south-west; one in the north-east, etc. Data points with different authorities might overlap, e.g. where two data points are similar spatially but have different authorities, but they will be concentrated largely separately (i.e. there will be a 'buffer' between authorities with a decay from each authority centre). 

Temporally, authorities can 'appear' and 'disappear'. For example, one authority might be present between two dates. Then it will no longer be present. Another will be present/not present at different times. I presume that when multiple authorities are present they will 'compete' with one another, i.e. they won't be independent, with a 'buffer' appearing between authorities.

I have been able to model this spatially, but not when including the temporal element (the model fitting fails):

fit_data <- readRDS("M_data_sample.rds")
mesh <- readRDS("M_mesh.rds")
spde <- readRDS("M_spde.rds")

components <- ~
  Authority(Authority_new, model = "factor_full") +
  Elevation(elevation_km, model = "linear") +
  River(rivers_1_dist2_km, model = "linear") +
  field(geometry, model = spde, replicate= Authority_new) +
  unique_id(location_id, model="iid", fixed = T)

# field, replicate is expected to fit a surface to each Authority_new, assuming that they have different spatial fields. to take into account that they are 'present' in different locations of the country.

fit_poi <- bru(
  components = components,
  formula = obs ~ -1 + unique_id + field + Authority + Elevation + River + offset(log_count),
  family = "poisson",
  data = fit_data,
  options = list(control.compute = list(dic = FALSE, waic = FALSE, cpo = FALSE)))

# include offset(log_count) to take into account the count associated with each data point. I'm unsure if this is correct, though?

For the spatio-temporal model, I do:

fit_data <- readRDS("M_data_sample.rds")
mesh <- readRDS("M_mesh.rds")
spde <- readRDS("M_spde.rds")

unique_times <- seq(-125.5, 37.5, by = 0.5) # do this so i can then predict authority extents at unobserved time stamps.
time_map <- setNames(seq_along(unique_times), unique_times)
fit_data$time_idx <- time_map[as.character(fit_data$Mid.Date)]

components <- ~
  Authority(Authority_new, model = "iid") +
  Elevation(elevation_km, model = "linear") +
  River(rivers_1_dist2_km, model = "linear") + 
  field(geometry,
        model = spde,
        replicate = Authority_new,
        group = time_idx,
        control.group = list(model = "ar1"), hyper = prec.prior) +
  unique_id(location_id, model = "iid", fixed = TRUE)

# for spde, allow it to be grouped by time_idx using an ar1 model, replicated across each Authority.
prec.prior <- list(prec = list(param = c(0.001, 0.001)))

fit_poi5 <- bru(
  components = components,
  formula = obs ~ -1 + unique_id + field + Authority + Elevation + River + offset(log_count),
  family = "poisson",
  data = fit_data,
  options = list(control.compute = list(dic = FALSE, waic = FALSE, cpo = FALSE), verbose = TRUE))

A summary of my data/full dataset (mesh, spde etc) is attached.

Any help would be greatly appreciated.

Thank you.

Joe
summary.PNG
M_data_sample.rds
M_mesh.rds
M_spde.rds

Joe Lewis

unread,
Oct 2, 2025, 12:33:57 PM (yesterday) Oct 2
to R-inla discussion group
Dear all,

The returned error from:

fit_data <- readRDS("M_data_sample.rds")
mesh <- readRDS("M_mesh.rds")
spde <- readRDS("M_spde.rds")

unique_times <- seq(-125.5, 37.5, by = 0.5)
time_map <- setNames(seq_along(unique_times), unique_times)
fit_data$time_idx <- time_map[as.character(fit_data$Mid.Date)]

components <- ~
  Authority(Authority_new, model = "iid") +
  Elevation(elevation_km, model = "linear") +
  River(rivers_1_dist2_km, model = "linear") +
  field(geometry,
        model = spde,
        replicate = Authority_new,
        group = time_idx,
        control.group = list(model = "ar1"), hyper = prec.prior) +
  unique_id(location_id, model = "iid", fixed = TRUE)

prec.prior <- list(prec = list(param = c(0.001, 0.001)))

fit_poi5 <- bru(
  components = components,
  formula = obs ~ -1 + unique_id + field + Authority + Elevation + River + offset(log_count),
  family = "poisson",
  data = fit_data,
  options = list(control.compute = list(dic = FALSE, waic = FALSE, cpo = FALSE), verbose = TRUE))

Compute initial values...
        Iter[0] RMS(err) = 1.000, update with step-size = 0.301
        Iter[1] RMS(err) = 0.452, update with step-size = 0.453
        Iter[2] RMS(err) = 0.510, update with step-size = 0.325
        Iter[3] RMS(err) = 0.577, update with step-size = 0.487
        Iter[4] RMS(err) = 0.598, update with step-size = 0.333
        Iter[5] RMS(err) = 0.617, update with step-size = 0.494
        Initial values computed in 0.1304 seconds
                x[0] = 0.1932
                x[1] = 0.2002
                x[2] = 0.2007
                x[3] = 0.2004
                x[4] = 0.1998
                x[1382359] = 0.2991
                x[1382360] = 0.1902
                x[1382361] = -0.7424
                x[1382362] = 0.1750
                x[1382363] = 0.0077

Optimise using DEFAULT METHOD
Segmentation fault (core dumped)
Warning in bru_log_warn(paste0("iinla: Problem in inla:\n", result)) :
  iinla: Problem in inla:
Error in inla.core.safe(formula = formula, family = family, contrasts = contrasts,  :
  The inla-program exited with an error. Unless you interupted it yourself, please rerun with verbose=TRUE and check the output carefully.
  If this does not help, please contact the developers at <he...@r-inla.org>.
The inla program failed and the maximum number of tries has been reached.
iinla: Problem in inla:
1: bru(components = components, formula = obs ~ -1 + unique_id + field + Auth [...]
2: iinla(model = info[["model"]], lhoods = info[["lhoods"]], inputs = info[[" [...]
3: fm_try_callstack(...)
4: do.call(INLA::inla, inla.options.merged, envir = environment(model$effects))
5: (function (formula = NULL, family = "gaussian", contrasts = NULL, data = N [...]
6: inla.core.safe(formula = formula, family = family, contrasts = contrasts,  [...]
7: stop(paste0(r$message, "\n", "The inla program failed and the maximum numb [...]
iinla: Giving up and returning last successfully obtained result for diagnostic purposes.

hopefully this might help/give some guidance to the issue?

Thank you.

Kind regards,
Joe

INLA help

unread,
Oct 2, 2025, 2:03:21 PM (yesterday) Oct 2
to Joe Lewis, R-inla discussion group
Please retry with the new testing version of today 

Haavard Rue
--
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/e26ee3d2-0e3e-449b-ac73-e526ca542c40n%40googlegroups.com.

Joe Lewis

unread,
Oct 2, 2025, 6:35:28 PM (23 hours ago) Oct 2
to R-inla discussion group
Dear Haavard,

Thank you for the reply.

I updated to the test version of INLA using 'inla.upgrade(testing=TRUE)'

The new error message:
 *** inla.core.safe:  The inla program failed, but will rerun in case better initial values may help. try=1/1

Warning in bru_log_warn(paste0("iinla: Problem in inla:\n", result)) :
  iinla: Problem in inla:
Error in inla.core.safe(formula = formula, family = family, contrasts = contrasts,  :
  Unknown keyword in `hyper' ` prec '. Must be one of  theta theta1 theta2 theta3 theta4 theta5 theta6 theta7 theta8 theta9 theta10 theta11 theta12 theta13 theta14 theta15 theta16 theta17 theta18 theta19 theta20 theta21 theta22 theta23 theta24 theta25 theta26 theta27 theta28 theta29 theta30 theta31 theta32 theta33 theta34 theta35 theta36 theta37 theta38 theta39 theta40 theta41 theta42 theta43 theta44 theta45 theta46 theta47 theta48 theta49 theta50 theta51 theta52 theta53 theta54 theta55 theta56 theta57 theta58 theta59 theta60 theta61 theta62 theta63 theta64 theta65 theta66 theta67 theta68 theta69 theta70 theta71 theta72 theta73 theta74 theta75 theta76 theta77 theta78 theta79 theta80 theta81 theta82 theta83 theta84 theta85 theta86 theta87 theta88 theta89 theta90 theta91 theta92 theta93 theta94 theta95 theta96 theta97 theta98 theta99 theta100 t1 t2 t3 t4 t5 t6 t7 t8 t9 t10 t11 t [... truncated]

iinla: Giving up and returning last successfully obtained result for diagnostic purposes.


I assume this is related to:
components <- ~
  Authority(Authority_new, model = "iid") +
  Elevation(elevation_km, model = "linear") +
  River(rivers_1_dist2_km, model = "linear") +
  field(geometry,
        model = spde,
        replicate = Authority_new,
        group = time_idx,
        control.group = list(model = "ar1"), hyper = prec.prior) +
  unique_id(location_id, model = "iid", fixed = TRUE)

prec.prior <- list(prec = list(param = c(0.001, 0.001)))


So I ran it without the prec.prior

-library/4.4/INLA/bin/linux/64bit/libreadline.so.8)
/home/jl2094/R/x86_64-pc-linux-gnu-library/4.4/INLA/bin/linux/64bit/inla.mkl: /lib/x86_64-linux-gnu/libc.so.6: version `GLIBC_2.38' not found (required by /home/jl2094/R/x86_64-pc-linux-gnu-library/4.4/INLA/bin/linux/64bit/libtirpc.so.3)
/home/jl2094/R/x86_64-pc-linux-gnu-library/4.4/INLA/bin/linux/64bit/inla.mkl: /lib/x86_64-linux-gnu/libc.so.6: version `GLIBC_2.38' not found (required by /home/jl2094/R/x86_64-pc-linux-gnu-library/4.4/INLA/bin/linux/64bit/libgomp.so.1)
/home/jl2094/R/x86_64-pc-linux-gnu-library/4.4/INLA/bin/linux/64bit/inla.mkl: /lib/x86_64-linux-gnu/libc.so.6: version `GLIBC_2.38' not found (required by /home/jl2094/R/x86_64-pc-linux-gnu-library/4.4/INLA/bin/linux/64bit/libgssapi_krb5.so.2)
/home/jl2094/R/x86_64-pc-linux-gnu-library/4.4/INLA/bin/linux/64bit/inla.mkl: /lib/x86_64-linux-gnu/libc.so.6: version `GLIBC_2.38' not found (required by /home/jl2094/R/x86_64-pc-linux-gnu-library/4.4/INLA/bin/linux/64bit/libkrb5.so.3)
/home/jl2094/R/x86_64-pc-linux-gnu-library/4.4/INLA/bin/linux/64bit/inla.mkl: /lib/x86_64-linux-gnu/libc.so.6: version `GLIBC_2.38' not found (required by /home/jl2094/R/x86_64-pc-linux-gnu-library/4.4/INLA/bin/linux/64bit/libk5crypto.so.3)
/home/jl2094/R/x86_64-pc-linux-gnu-library/4.4/INLA/bin/linux/64bit/inla.mkl: /lib/x86_64-linux-gnu/libc.so.6: version `GLIBC_2.38' not found (required by /home/jl2094/R/x86_64-pc-linux-gnu-library/4.4/INLA/bin/linux/64bit/libcom_err.so.2)
/home/jl2094/R/x86_64-pc-linux-gnu-library/4.4/INLA/bin/linux/64bit/inla.mkl: /lib/x86_64-linux-gnu/libc.so.6: version `GLIBC_2.38' not found (required by /home/jl2094/R/x86_64-pc-linux-gnu-library/4.4/INLA/bin/linux/64bit/libkrb5support.so.0)
/home/jl2094/R/x86_64-pc-linux-gnu-library/4.4/INLA/bin/linux/64bit/inla.mkl: /lib/x86_64-linux-gnu/libc.so.6: version `GLIBC_2.36' not found (required by /home/jl2094/R/x86_64-pc-linux-gnu-library/4.4/INLA/bin/linux/64bit/libresolv.so.2)
/home/jl2094/R/x86_64-pc-linux-gnu-library/4.4/INLA/bin/linux/64bit/inla.mkl: /lib/x86_64-linux-gnu/libc.so.6: version `GLIBC_ABI_DT_RELR' not found (required by /home/jl2094/R/x86_64-pc-linux-gnu-library/4.4/INLA/bin/linux/64bit/libresolv.so.2)

Warning in bru_log_warn(paste0("iinla: Problem in inla:\n", result)) :
  iinla: Problem in inla:
Error in inla.core.safe(formula = formula, family = family, contrasts = contrasts,  :
  The inla-program exited with an error. Unless you interupted it yourself, please rerun with verbose=TRUE and check the output carefully.
  If this does not help, please contact the developers at <he...@r-inla.org>.
The inla program failed and the maximum number of tries has been reached.
iinla: Problem in inla:
1: bru(components = components, formula = obs ~ -1 + unique_id + field + Auth [...]
2: iinla(model = info[["model"]], lhoods = info[["lhoods"]], inputs = info[[" [...]
3: fm_try_callstack(...)
4: do.call(INLA::inla, inla.options.merged, envir = environment(model$effects))
5: (function (formula = NULL, family = "gaussian", contrasts = NULL, data = N [...]
6: inla.core.safe(formula = formula, family = family, contrasts = contrasts,  [...]
7: stop(paste0(r$message, "\n", "The inla program failed and the maximum numb [...]
iinla: Giving up and returning last successfully obtained result for diagnostic purposes.

I'm running it on a Linux HPC if that helps? do I need to e.g. do: 

> inla.binary.install()
* Looking for Version_25.10.01 and os='<choose interactively>'
  Available alternatives:
         Alternative 1  is  ./CentOS Linux-7 (Core)/Version_25.10.01/64bit.tgz
         Alternative 2  is  ./Fedora Linux-43 (Workstation Edition Prerelease)/Version_25.10.01/64bit.tgz
         Alternative 3  is  ./Rocky Linux-10.0 (Red Quartz)/Version_25.10.01/64bit.tgz
         Alternative 4  is  ./Rocky Linux-8.10 (Green Obsidian)/Version_25.10.01/64bit.tgz
         Alternative 5  is  ./Rocky Linux-9.6 (Blue Onyx)/Version_25.10.01/64bit.tgz
         Alternative 6  is  ./Ubuntu-22.04.5 LTS (Jammy Jellyfish)/Version_25.10.01/64bit.tgz
         Alternative 7  is  ./Ubuntu-25.04 (Plucky Puffin)/Version_25.10.01/64bit.tgz
  Chose alternative [1:7]
        1: 1.

I'm  running the model currently after installing Alternative 1. I will let you know how I get on.

Thank you .

Kind regards,
Joe

INLA help

unread,
Oct 2, 2025, 6:47:40 PM (23 hours ago) Oct 2
to Joe Lewis, R-inla discussion group
I see you’re using  R-4.4.  Please upgrade to 4.5 

Haavard Rue

Finn Lindgren

unread,
Oct 2, 2025, 6:53:49 PM (22 hours ago) Oct 2
to R-inla discussion group
For the issue with setting the prior (unrelated to the R version issue etc), the issue is most likely that the group ar1 model doesn't have a separate precision parameter;
The precision for the joint main+group+replicate model is set for the main model, i.e. the parameters of the "spde" model, from inla.spde2.pcmatern()

Finn



--
Finn Lindgren
email: finn.l...@gmail.com
Reply all
Reply to author
Forward
0 new messages