Joint likelihood with different meshes

[high...@highstat.com]

Hello,

I have the following example:

There are two areas, far apart from each other. Each area contains multiple sites. I have created two separate meshes — one for each area — and therefore have two different projector matrices and two different SPDEs.

My goal is to test whether the spatial correlation structure differs between the two areas. More precisely, I want to compare:

• a model where the SPDE hyperparameters (range and marginal standard deviation) are allowed to differ by area, and

• a model where the hyperparameters are constrained to be the same across both areas.

My current model looks like this:

I1 <- inla(CombinedData ~ -1 +
                           Intercept.alp +
                           f(AWTD.alp,
                             model = "rw2",
                             scale.model = TRUE,
                             hyper = PCPrior,
                             constr = TRUE) +
                           f(w.alp, model = spde.alp) +
                           Intercept.deh +
                           f(AWTD.deh,
                             model = "rw2",
                             scale.model = TRUE,
                             hyper = PCPrior,
                             constr = TRUE) +
                          f(w.deh, model = spde.deh),
           family = c("poisson", "poisson"),
           data = inla.stack.data(Stack),
           control.compute = ControlCompute,
           control.predictor = list(A = inla.stack.A(Stack),
                                    compute = TRUE,
                                    link = 1))

Given that I use two different meshes, how do do I force the hyperparameters of the two spatial terms   f(w.alp, model = spde.alp)  and f(w.deh, model = spde.deh) to be the same? I thought that this was a simple 'replicate'  option, but that does not seem to work with different meshes. I tried the 'group' option, but without success.

I know that this is straightforward to do if both areas are represented on a single common mesh — but my interest here is specifically in keeping the meshes separate (each adapted to its own spatial extent) while still sharing hyperparameters between the two spatial fields.

Alain

[Finn Lindgren]

Hi,

by the easiest option is likely to use a single mesh with disconnected components.

The only snag would be if the regions are too close to allow for a boundary effect mesh extension without overlapping the two parts.

If the regions are too close for that, you could use the barrier method to effectively disconnect the regions in terms of correlation.

A more proper solution to this is to extend the fmesher structures to non-tensor product constructions, where one has an indexing dimension for selecting values from different meshes.

The coding needed to make such a function space in fmesher isn't difficult, but someone would need to do it...

It would be similar to the existing fm_tensor() space, but with different logic for the location&index -- mesh&index fm_bary, fm_basis, and fm_fem methods.

With that done, all the inlabru methods for talking to fmesher objects should solve the rest.
(I would not want to debug code for this written in "plain INLA"; inlabru would be much easier)

But as per the above, if the mesh components don't need to overlap, a single mesh with disconnected components already works as intended.

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/2247506b-4ec4-42c5-8664-cb4f49324c6an%40googlegroups.com.

--

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

[Anderson Bussing]

Hello all, 

Inspired by this question and answer, I tried to implement this and verify if it is working correctly by using a simple simulation:

I believe it is working as intended but I am worried there may be some lurking incorrect thing in the implementation that is preventing me from getting quite correct results.

Assuming it is correct though, I was wondering if there is any way for this type of inlabru model to be extended to a situation where we actually want to give each of these dates a little bit of flexibility with their Matern hyperparameters. I don't want to go as far as letting them each learn their own Matern hyperparameters individually, but I'm worried the model in the below simulation may be a bit rigid for real-world data. 

Or do you possibly have any recommendations for a situation where we have collected covariates and response data on many different dates, usually dates relatively far from each other in time, where the data is collected from different locations of a common large geographic area on each date. I can fit any individual date just fine and learn about how the covariates and spatial dependence impact the response variable on that date, but making use of these "replicates" and having them pool information to give a more accurate picture has been hard until I came upon this idea you all gave about multiple meshes and common Matern hyperparameters. But now I feel I need to give them a BIT more wiggle room.

Thanks so much for any advice!

best regards,

Anderson 

library(mvtnorm)
library(ggplot2)
library(dplyr)
library(fmesher)
library(Matrix)
library(INLA)
library(inlabru)

# Simulate 10 different clouds of points

Nclouds <- 10

Npts_per_cloud <- 100

ourdf <- NULL

for (j in 1:Nclouds){
 
  cntr <- runif(2, 0, 10)
 
  margsig <- runif(1, 0.5, 1.5)
 
  rhotemp <- runif(1, -0.6, 0.6)
 
  sigtemp <- rbind(c(margsig^2, rhotemp*margsig^2), c(rhotemp*margsig^2, margsig^2))
 
  xytemp <- rmvnorm(n = Npts_per_cloud, mean = cntr, sigma = sigtemp)
 
  tempdf <- data.frame(loc1 = xytemp[,1], loc2 = xytemp[,2])
  tempdf$date <- paste0("date", j)
 
  ourdf <- rbind(ourdf, tempdf)
 
}

ourdf$date <- as.factor(ourdf$date)


ggplot(ourdf, aes(x = loc1, y = loc2, color = date)) +
  geom_point(alpha = 0.7, size = 1.8) +
  coord_equal() +
  theme_minimal() +
  labs(
    title = "Locations by Date",
    x = "loc1",
    y = "loc2",
    color = "date"
  )





# Make a mesh for each date


dates10 <- sort(unique(ourdf$date))


cutoff_val <- 0.05
max_edge_val <- c(0.5, 2.0)
offset_val <- c(1, 3)
min_angle <- 21
convex_val <- -0.05

mesh_list <- vector("list", length(dates10))
names(mesh_list) <- dates10

for (k in seq_along(dates10)) {
  d <- dates10[k]
  dat_d <- ourdf %>% filter(date == d)
 
  loc_d <- unique(as.matrix(dat_d[, c("loc1", "loc2")]))
 
  bnd_d <- fm_nonconvex_hull(loc_d, convex = convex_val)
 
  mesh_list[[k]] <- fm_mesh_2d_inla(
    loc = loc_d,
    boundary = bnd_d,
    cutoff = cutoff_val,
    max.edge = max_edge_val,
    offset = offset_val,
    min.angle = min_angle
  )
}

stopifnot(all(!vapply(mesh_list, is.null, logical(1))))

xlim_all <- range(ourdf$loc1)
ylim_all <- range(ourdf$loc2)

op <- par(
  mfrow = c(5, 2),
  mar = c(0.6, 0.6, 1.6, 0.2),
  oma = c(0.2, 0.2, 1.0, 0.2)
)

for (k in seq_along(dates10)) {
  d <- dates10[k]
  m <- mesh_list[[k]]
  dat_d <- ourdf %>% filter(date == d)
 
  plot(m, asp = 1, xlim = xlim_all, ylim = ylim_all,
       main = "", axes = FALSE, xlab = "", ylab = "", bty = "n",
       draw.segments = FALSE)
 
  points(dat_d$loc1, dat_d$loc2, pch = 16, cex = 0.6, col = "blue")
 
  mtext(d, side = 3, line = 0.75, cex = 0.85)
  mtext(paste0("points=", nrow(dat_d), " | vertices=", nrow(m$loc)),
        side = 3, line = 0.05, cex = 0.65)
}

mtext("Per-date meshes (simulated loc1/loc2)", outer = TRUE, cex = 0.95, line = 0.2)
par(op)




# Simulate latent surfaces for each date

matern_cov_from_locs <- function(loc_mat, sigma2 = 1, rho = 1, nu = 1.0, nugget = 1e-8) {

  D <- as.matrix(dist(loc_mat))

  x <- D / rho
 
  Cor <- matrix(0, nrow(D), ncol(D))
  Cor[x == 0] <- 1
 
  idx <- (x > 0)
  const <- (2^(1 - nu)) / gamma(nu)
  Cor[idx] <- const * (x[idx]^nu) * besselK(x[idx], nu)

  Sigma <- sigma2 * Cor

  diag(Sigma) <- diag(Sigma) + nugget
 
  Sigma
}



for (j in 1:Nclouds){
 
  date_idcs <- which(ourdf$date == paste0("date",j))
 
  maternsig <- matern_cov_from_locs(loc_mat = as.matrix(ourdf[date_idcs, c("loc1", "loc2")]), sigma2 = 0.5, rho = 0.5, nu = 1)
 
  zsurf <- as.vector(rmvnorm(n = 1, mean = rep(0, length(date_idcs)), sigma = maternsig))
 
  zsurf <- zsurf - mean(zsurf)
 
  ourdf[date_idcs, "zsurf"] <- zsurf
 
}




zlim_all <- range(ourdf$zsurf, na.rm = TRUE)
col_fun  <- colorRampPalette(c("navy", "white", "firebrick"))
cols     <- col_fun(100)



op <- par(
  mfrow = c(5, 2),
  mar  = c(0.6, 0.6, 1.6, 0.2),
  oma  = c(0.2, 0.2, 1.0, 2.5)  # room on right for colorbar
)

for (k in seq_along(dates10)) {
  d <- dates10[k]
  dat_d <- ourdf %>% dplyr::filter(date == d)

  z_scaled <- (dat_d$zsurf - zlim_all[1]) / diff(zlim_all)
  z_idx    <- pmax(1, pmin(100, floor(z_scaled * 99) + 1))
 
  plot(
    dat_d$loc1, dat_d$loc2,
    col = cols[z_idx],
    pch = 16,
    cex = 0.7,
    asp = 1,
    xlim = xlim_all,
    ylim = ylim_all,
    axes = FALSE,
    xlab = "",
    ylab = "",
    main = "",
    bty = "n"
  )
 
  mtext(d, side = 3, line = 0.75, cex = 0.85)

}

mtext("Latent surface by Date", outer = TRUE, cex = 0.95, line = 0.2)



# Create covariates to go in linear predictor

ourdf$covar1 <- runif(nrow(ourdf), 1, 5)
ourdf$covar2 <- runif(nrow(ourdf), 0, 1)


# choose coefficients for these covariates

intcpt <- -0.5
coeff1 <- 0.3
coeff2 <- -0.8

# simulate bernoulli data with cloglog link

ourdf$linpred <- intcpt + coeff1 * ourdf$covar1 + coeff2 * ourdf$covar2

ourdf$berny <- rbinom(n = nrow(ourdf), size = 1, prob = 1-exp(-exp(ourdf$linpred)))



# to use inlabru with multiple meshes we need the following function rewrites

fm_mesh_indexed <- function(meshes) {
  stopifnot(is.list(meshes), length(meshes) >= 1)
  node_counts  <- vapply(meshes, function(m) nrow(m$loc), integer(1))
  node_offsets <- c(0L, cumsum(node_counts))[seq_along(meshes)]
  all_loc <- do.call(rbind, lapply(meshes, `[[`, "loc"))
  all_tv  <- do.call(rbind, lapply(seq_along(meshes), function(k) meshes[[k]]$graph$tv + node_offsets[k]))
  structure(
    list(
      meshes = meshes,
      names = names(meshes),
      node_counts = node_counts,
      node_offsets = node_offsets,
      n = nrow(all_loc),
      loc = all_loc,
      graph = list(tv = all_tv),
      manifold = "R2",
      meta = list(is.refined = TRUE)
    ),
    class = c("fm_mesh_indexed", "inla.mesh", "fm_mesh_2d")
  )
}
.fm_map_mesh_ids <- function(mesh_idx, mesh_ids) {
  if (is.character(mesh_ids)) match(mesh_ids, mesh_idx$names) else mesh_ids
}
fm_dof.fm_mesh_indexed <- function(x, ...) sum(vapply(x$meshes, function(m) nrow(m$loc), integer(1)))
fm_manifold_get.fm_mesh_indexed <- function(x, ...) "R2"



fm_bary.fm_mesh_indexed <- function(mesh, loc, crs = NULL, ..., max_batch_size = NULL) {
  if (inherits(loc, "fm_bary")) return(loc)
  loc_df <- as.data.frame(loc)
  if (!all(c("x","y","mesh_id") %in% names(loc_df)))
    stop("loc must contain columns named 'x','y','mesh_id'")
  coords <- as.matrix(loc_df[, c("x","y")])
  mesh_ids <- loc_df$mesh_id
  idx_mesh <- .fm_map_mesh_ids(mesh, mesh_ids)
  n_obs <- nrow(coords)
  tri <- rep(NA_integer_, n_obs)
  where <- matrix(NA_real_, n_obs, 3)
  unique_k <- unique(idx_mesh[!is.na(idx_mesh)])
  for (k in unique_k) {
    rows_k <- which(idx_mesh == k)
    bary_k <- fm_bary(mesh$meshes[[k]], coords[rows_k,,drop=FALSE],
                      crs = crs, ..., max_batch_size = max_batch_size)
    tri[rows_k] <- bary_k$index
    where[rows_k,] <- as.matrix(bary_k$where)
  }
  fm_bary(tibble(mesh_id = idx_mesh, index = tri, where = where),
          extra_class = "fm_bary_indexed")
}



fm_basis.fm_mesh_indexed <- function(x, loc, weights = NULL, ..., full = FALSE) {
  loc_df <- as.data.frame(loc)
  if (!all(c("x","y","mesh_id") %in% names(loc_df)))
    stop("loc must contain 'x','y','mesh_id'")
  coords <- as.matrix(loc_df[, c("x","y")])
  mesh_ids <- loc_df$mesh_id
  idx_mesh <- .fm_map_mesh_ids(x, mesh_ids)
  n_obs <- nrow(coords)
  n_tot <- fm_dof(x)
  i <- integer(0); j <- integer(0); xv <- numeric(0); ok <- logical(n_obs)
  for (r in 1:n_obs) {
    mk <- idx_mesh[r]
    if (is.na(mk)) next
    submesh <- x$meshes[[mk]]
    bary <- fm_bary(submesh, coords[r,,drop=FALSE], ...)
    if (!is.na(bary$index[1])) {
      tri <- submesh$graph$tv[bary$index[1], ]
      v_off <- x$node_offsets[mk]
      gverts <- v_off + tri
      i <- c(i, rep(r, 3))
      j <- c(j, gverts)
      xv <- c(xv, as.numeric(bary$where[1, ]))
      ok[r] <- TRUE
    }
  }
  A <- sparseMatrix(i=i, j=j, x=xv, dims=c(n_obs, n_tot))
  if (!is.null(weights)) {
    if (length(weights) == 1) A <- A * weights
    else if (length(weights) == n_obs) A <- Diagonal(n=n_obs, x=weights) %*% A
    else stop("weights must be 1 or nrow(loc)")
  }
  if (full) { list(A=A, ok=ok) } else A
}



fm_fem.fm_mesh_indexed <- function(mesh, order = 2, ...) {
  stopifnot(inherits(mesh, "fm_mesh_indexed"))
  K <- length(mesh$meshes); if (K == 0) stop("no submeshes")
  fem_list <- lapply(mesh$meshes, fm_fem, order=order, ...)
  piece_names <- names(fem_list[[1]])
  for (k in seq_len(K)) if (!identical(names(fem_list[[k]]), piece_names))
    stop("FEM piece names differ across submeshes")
  .is_sparse <- function(M) inherits(M, c("dgCMatrix","dgTMatrix","dsCMatrix","dtCMatrix"))
  out <- vector("list", length(piece_names)); names(out) <- piece_names
  for (nm in piece_names) {
    blocks <- lapply(fem_list, `[[`, nm)
    if (.is_sparse(blocks[[1]])) {
      out[[nm]] <- as(Matrix::bdiag(blocks), "CsparseMatrix")
    } else if (is.numeric(blocks[[1]]) || (is.matrix(blocks[[1]]) && is.numeric(blocks[[1]]))) {
      vecs <- lapply(blocks, function(v) if (is.matrix(v)) as.numeric(v[,1]) else as.numeric(v))
      out[[nm]] <- unlist(vecs, use.names = FALSE)
    } else stop(paste0("Unsupported FEM block: ", nm))
  }
  if (!("cc" %in% names(out)) && ("c0" %in% names(out))) out$cc <- out$c0
  out
}


make_sumzero_constraints <- function(mesh_idx, weighted = TRUE) {
  K <- length(mesh_idx$meshes)
  nv_k <- mesh_idx$node_counts
  off <- mesh_idx$node_offsets
  Ntot <- sum(nv_k)
  if (weighted) {
    fem <- fm_fem(mesh_idx)
    w_all <- as.numeric(colSums(fem$c0))
  } else {
    w_all <- rep(1, Ntot)
  }
  i <- integer(0); j <- integer(0); x <- numeric(0)
  for (k in seq_len(K)) {
    cols_k <- (off[k] + 1L):(off[k] + nv_k[k])
    i <- c(i, rep(k, nv_k[k])); j <- c(j, cols_k); x <- c(x, w_all[cols_k])
  }
  list(A = sparseMatrix(i=i,j=j,x=x,dims=c(K, Ntot)), e = rep(0, K))
}





# now we can set up the inlabru

mesh_idx <- fm_mesh_indexed(mesh_list)


zc <- make_sumzero_constraints(mesh_idx, weighted = TRUE)

spde_multi <- inla.spde2.pcmatern(
  mesh = mesh_idx,
  alpha = 2,
  prior.range = c(1, 0.9),
  prior.sigma = c(1, 0.1)
)


components <- ~
  intcpt(1) +
  beta1(covar1) +
  beta2(covar2) +
  spatial_smooth(
    cbind(x = loc1, y = loc2, mesh_id = date),
    model = spde_multi,
    extraconstr = zc
  )



joint_model <- bru(
  components = components,
  bru_obs(
    formula = berny ~ intcpt + beta1 + beta2 + spatial_smooth,
    family = "binomial",
    data = ourdf,
    control.family = list(link = "cloglog")
  ),
  options = list(
    control.compute = list(config = TRUE, dic = TRUE, waic = TRUE),
    control.inla = list(strategy = "gaussian", int.strategy = "ccd"),
    verbose = FALSE
  )
)



summary(joint_model)

[Finn Lindgren]

I have limited time to investigate this right now as I have exam marking plus need to get a new inlabru version onto CRAN, but could you please take a look at the fm_collect() feature? https://inlabru-org.github.io/fmesher/reference/fm_collect.html

It was implemented for this specific application, and should be constructing the right fem matrices, and inlabru mapper, etc.

The mapper takes a tibble or list with two elements, loc, and index.

The only thing that your code does that might be missing from it is an automated way to get define per-submesh constraints.

Finn

On 14 Dec 2025, at 21:39, Anderson Bussing <andyb...@gmail.com> wrote:

Hello all, 

To view this discussion, visit https://groups.google.com/d/msgid/r-inla-discussion-group/c28b42c3-76d6-4ded-aceb-26f3f62068e8n%40googlegroups.com.

[Anderson Bussing]

Thank you so much! I am looking at fm_collect() now.

About the second part of my question-- the part about the different dates (meshes) maybe needing a little wiggle room Matern parameters-- this is not urgent for you to answer, but is there capability in inlabru for such a thing? I was thinking of writing an rgeneric to put a prior on the range parameters of the different date's Materns, but that would violate the INLA assumption of the latent layer being Gaussian.  Or maybe there is another type of inlabru model for dealing with these types of "replicates" where the locations of the observations on each date are different but they're coming from a common geographic domain and the interest is not in the spatio but not temporal aspect.

Please ignore this question for as long as you need (or forever of course)... it is not urgent!

Thanks again,

Anderson

[Finn Lindgren]

The technique in the classic inla.spde2.matern code for allowing covariates in the range and variance parameters would work also for fm_collect; the B.tau and B.kappa matrices just need to be set to give separate parameter effects for each submesh.

No special inlabru code is needed, and r/cgeneric also isn’t needed. One would only need r/cgeneric code to implement e.g. pcpriors instead of lognormals.

Finn

On 14 Dec 2025, at 22:39, Anderson Bussing <andyb...@gmail.com> wrote:

Thank you so much! I am looking at fm_collect() now.

To view this discussion, visit https://groups.google.com/d/msgid/r-inla-discussion-group/476b4177-bb42-4c8f-a3e1-f8700f080623n%40googlegroups.com.

[Finn Lindgren]

The case of a fixed spatial domain and temporal evolution is handled either by separambe covariance models via the “group” feature, or with nonseparable covariance models, from the INLAspacetime package.

Finn

On 14 Dec 2025, at 23:21, Finn Lindgren <finn.l...@gmail.com> wrote:

The technique in the classic inla.spde2.matern code for allowing covariates in the range and variance parameters would work also for fm_collect; the B.tau and B.kappa matrices just need to be set to give separate parameter effects for each submesh.

[Anderson Bussing]

Thank you so much for your responses.

I am worried with the "separate parameter effects for each submesh" idea that dates with not much information in the observations will struggle if we don't have any pooling for the matern hyperparameters. That is why the common Matern hyperparameters among all dates was appealing. But then I was wondering about giving each date a little deviation around the pooled common Matern hyperparameter values but not so much deviation to the point that uninformative dates would have to learn Matern hyperparams completely on their own-- uninformative dates should just get their hyperparams sucked to to the common values.

My apologies for the continued bother; this is my last question I promise.

best regards

Anderson

[Finn Lindgren]

This is doable with the B.tau/B.kappa construction, as you can have joint parameters controlling all of them, as well as extra parameters for the perturbations, with stricter priors.

Finn

On 15 Dec 2025, at 00:03, Anderson Bussing <andyb...@gmail.com> wrote:

Thank you so much for your responses.

To view this discussion, visit https://groups.google.com/d/msgid/r-inla-discussion-group/f2276cf8-0b56-4755-af82-af401ec3e49cn%40googlegroups.com.