Transforming R-INLA code to Inlabru

[Bamidele Toba]

In INLA, I can provide a projection matrix A for a spatial effect in both estimation and prediction stacks. For example:

In Estimation Stack: A = list(1, A_matrix, diagonal(n)), where 1 represents the intercept, A_ matrix is an inla.spde.make.A from the observed data and diagonal(n) represents an IID random effect,

In Prediction Stack: A = list(1, A_pred, diagonal(n)), where 1 represents the intercept, A_pred is the SPDE projection for the prediction locations and diagonal(n)  represents the IID effect at the prediction locations.

Similarly, I can also provide a custom projection matrix for the spatial effect in both estimation and prediction stacks, where A_pred is a custom projection corresponding to the latent locations in A_matrix.

In inlabru, there is no direct way to supply the custom projection matrix A. How can I reformulate INLA codes to inlabru so that inlabru correctly maps a custom spatial latent effect to each observation for both estimation and prediction?

[Finn Lindgren]

Hi Bamidele,

with inlabru, you normally don't need to supply those matrices at all, as they are constructed automatically based on your input data, which in your case should contain the geographic locations, and the index into the prediction locations; assuming you have an sf data frame "my_data" with geometry column called "geometry" (this is usually the case, but st_geometry(data) <- "geometry" can be used to set it to whatever name you want; "geometry" is the most common, and what all the inlabru examples assume) and an index column called "ID".

Also assume the sf object has the response variable, called "y", with NA values for the rows where you have new prediction locations.

Then you can do

  matern <- inla.spde2.pcmatern(mesh, ...)

  fit_and_predict <- bru(y ~ 0 + Intercept(1) + Field(geometry, model = matern) + RandomEffect(ID, model = "iid"), data = my_data, family = ...)

The resulting linear predictor & fitted values outputs will have the same order as in "my_data", so you can easily know where the predictions are.

This is _the_ canonical way to handle such models in inlabru.

All that said, there _is_ a way (actually multiple ways) to supply custom matrices in the inlabru model construction.
Where did you get that "In inlabru, there is no direct way to supply the custom projection matrix A."?

If that's being claimed somewhere, I'd like to know, so I can tell them to correct it, as it is not true.

One way is by explicitly multiplying with a matrix in the prediction expression, and another is to use a bm_matrix() mapper for the component.

Those are, however, more advanced options, and you should only use those if you have a special reason for it.

With the mapper approach (I'm using fm_basis() here instead of inla.spde.make.A(), as the latter is effectively deprecated for inlabru models):

  fit_and_predict <- bru(y ~ 0 + Intercept(1) +

    Field(fm_basis(mesh, geometry), model = matern, mapper = bm_matrix(fm_dof(mesh))) +

    RandomEffect(ID, model = "iid"), data = my_data,family = ...)

This is essentially what the default mapper (bm_fmesher()) for the spde models does; it calls the fm_basis() function that returns the A-matrix.
By computing the matrix based on the data supplied to bru() it will work also when calling predict(), whereas pre-computing the matrix would prevent normal operation of predict().

Explicit matrix multiplication:

  A_precomputed <- fm_basis(mesh, my_data$geometry)

  fit_and_predict <- bru(~ 0 + Intercept(1) +

    Field_basis_weights(seq_len(fm_dof(mesh)), model = matern, mapper = bm_index(fm_dof(mesh))) +

    RandomEffect(ID, model = "iid"),
    bru_obs(y ~ Intercept + A_precomputed %*% Field_basis_weights + RandomEffect, family = ..., data = my_data, allow_combine = TRUE)
)

Note: In the current inlabru version, this last version will be slower, as it (currently) forces inlabru to consider that the general R expression given as predictor expression in bru_obs() might be non-linear.

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, visit https://groups.google.com/d/msgid/r-inla-discussion-group/8908d7a9-49c2-469b-809e-9d9b902c9c96n%40googlegroups.com.

--

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

[Bamidele Toba]

Thank you for your feedback Finn,

 

In SPDE models with inlabru/INLA, the precision matrix is built automatically from the mesh and Matérn parameters, and the projection matrix maps the mesh nodes to the observations. I want to replicate this manually in inlabru using a generic1 latent effect: how can I explicitly set up a generic1 component where I provide (i) a custom Cmatrix to define the dependence structure, (ii) a diagonal identity mapping so each latent effect links to its own observation, and (iii) a projection matrix linking latent effects to the observed data?

This is the code

#Precision Matrix

Cmatx <- sparseMatrix(i=tri$i, j=tri$j, x=tri$x, dims=c(n, n))

hyp <- list(
theta1 = list(prior = "pc.prec", param = c(1, 0.01)),
theta2 = list(prior = "gaussian", param = c(0, 1))
)

data1[, id := 1:.N]

components <- ~
intercept(1) +
s(id,
model = "generic1",
Cmatrix = Cmatx, # (i) a custom Cmatrix to define the dependence structure
constr = FALSE,
hyper = hyp) +
eps(1, model = "iid")

fit_bru <- bru(
components = components,
family = "binomial",
data = nodes,
formula = y ~ intercept + covariate + s + eps,
Ntrials = data1$n,
options = list(
control.predictor = list(
A = list(
s = Diagonal(n) # (ii) a diagonal identity mapping so each latent effect links to its own observation
),
compute = TRUE
),
control.compute = list(dic = TRUE, waic = TRUE, cpo=TRUE, config=TRUE),
verbose = FALSE
)
)

But I get the error message below, and there are no missing values in Cmatx
“Error in Matrix::sparseMatrix(i = triplets$i, j = triplets$j, x = triplets$x, : 'i' and 'j' must not contain NA”

Hence, I could not proceed to predict and
(iii) How can I integrate a projection matrix linking latent effects to the observed data?

#The projection Matrix
k <- 5
knn <- get.knnx(data =location, query = data2_locs, k = k)
distance <- knn$nn.dist
a <- 1 / distance
A <- a/ rowSums(a)
# Build sparse S of size N x n
N <- nrow(data2)
i <- rep(1:N, each = k)
j <- as.vector(knn$nn.index)
x <- as.vector(a)
S <- sparseMatrix(i = i, j = j, x = x, dims = c(N, n))

pred_bru <- predict(fit_bru,
newdata = data2,
formula = ~ plogis(intercept + covariate + s + eps),
n.samples = 1000)

[Bamidele Toba]

Thank you for your reply. 

Meanwhile, in SPDE models with inlabru/INLA, the precision matrix is automatically built from the mesh and Matérn parameters, and the projection matrix maps the mesh nodes to the observations. I want to replicate this manually using a generic1 latent effect: how can I explicitly set up a generic1 component where I provide (i) a custom Cmatrix to define the dependence structure, (ii) a diagonal identity mapping so each latent effect links to its own observation, and (iii) a projection matrix linking latent effects to the observed data?

On Sunday, September 7, 2025 at 8:28:38 PM UTC+1 Finn Lindgren wrote:

[Finn Lindgren]

Hi,

Your component definition is lacking a "covariates" item; perhaps that's the source of the NA error.

Most features in R-INLA work the same in inlabru (e.g. the Cmatrix argument). It's mostly the mapping between latent variables and observations that is easier to specify in inlabru, as we usually do _not_ need to work explicitly with the matrices. Please see the documentation for inlabru at http://inlabru-org.github.io/inlabru/
That includes https://inlabru-org.github.io/inlabru/articles/bru_mapper.html that includes an example of even defining a custom mapper.

If you don't want to define a custom mapper (it looks like you really _should_ do that though, as it makes the code more general), you can use precomputed matrices, using mapper=bm_matrix(size of the model)

s(A_your_custom_estimation_data_mapping_matrix,

mapper=bm_matrix(nrow(Cmatx)),

model = "generic1",
Cmatrix = Cmatx, # (i) a custom Cmatrix to define the dependence structure
constr = FALSE,
hyper = hyp)

and you should _not_ supply a control.prediction A element.

In the prediction call, you then need to explicitly do the matrix multiplication operating on the latent variables directly (this wouldn't be needed if you define your own custom mapper for the s-component):

pred_bru <- predict(fit_bru,
newdata = data2,

formula = ~ plogis(intercept + covariate + A_your_prediction_matrix %*% s_latent + eps),
n.samples = 1000)

You'll need to do the same for other components that need a special matrix.
Since you're defining your own matrix weighting for the prediction, it seems to me that you should consider implementing a proper bru_mapper (see example in the link I gave above for how to properly define a custom mapper class).
You only need to implement one or two methods:
  ibm_n.your_mapper_class() should return the model size; if you store the size as an element "n" in the mapper object, you don't need to implement this method.
  ibm_jacobian.your_mapper_class() should return the mapping matrix.

Finn

To view this discussion, visit https://groups.google.com/d/msgid/r-inla-discussion-group/f4543f2d-8450-4b51-be9f-b2b591b25bb6n%40googlegroups.com.