INLA poisson model with many covaraites

[Alokesh Manna]

Hi INLA Group,

I was trying to model as follows:

For the i-th replicate and  time t, the count y_{i,t} can be modeled as 

y_{i,t} ~ Pois ( exp (f(time, rw1) + f(i, iid) + X'{i,s} \beta_{s} ) ), with \beta_{s} has a neighborhood structure generated by inla graph. Here we want to put a besag kernel for smoothing over \beta_{s}.

Generally, for besag we need to construct an ID for the locations s and implement besag. Here important to note that the covariates only depend on spatially, not temporally. The number of covariates is huge, for example, 700.

p.s. if we make id for each location s, and make it a dataframe in long format and try to implement, we say that y_{i,s,t} but actually the count only depends on i and t.

Can we implement such a model?

[Finn Lindgren]

Hi Alokesh,

your pseudo-math/code model definition isn't sufficiently clear to tell what model you actually want.

Can you spell out the model more carefully (preferably in a way that doesn't involve pseudo-code, but instead proper maths?

(There are lots of possible models and most syntactically valid code corresponds to _some_ model, but to know if it is the _intended_ model one needs a clear model definition that doesn't involve the code).

Your initial statement doesn't have "s" as an index on "y" but your RHS expression does, in a way that leaves it "dangling".
I can't tell if you want some weighted sum over all possible s, or if each "i" has a specific associated value of "s", is beta a single spatial field, or a separate field for each covariate, or something entirely different, etc.

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/2c469127-0101-473c-b9ff-696382fce496n%40googlegroups.com.

--

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

[Finn Lindgren]

To be more specific, it's the
  X'{i,s} \beta_{s}

formulation that is undefined. I have the feeling that you either intend some kind of sum over "s", _or_

that each i has an associated value that we could write as s_{i}, and that each covariate is spatial.

In that latter case, a more complete statement would be

where X_j(.) is the spatial function defining covariate nr j.

But since that expression could be replaced by
   \beta_{s_i} \sum_{j=1}^J X_j(s_i)

and the value of the sum pre-computed for each "i", I'm guessing that's _not_ what you intended.

Since your original statement is ambiguous, you need to spell out the details of it.

Finn

[Finn Lindgren]

Sorry, missed one line in the middle. It should have said this:

"

In that latter case, a more complete statement would be

  \sum_{j=1}^J X_j(s_i) \beta_{s_i}

where X_j(.) is the spatial function defining covariate nr j.
"

[Alokesh Manna]

Dear Professor Finn,

Let me know the following formulation works to resolve the particular problem?

P.S. I also conversed with Dr. Rue. But it seems the current INLA version does not support me A.local. I was looking an alternate way to resolve this problem.

"##

set.seed(42)

# ---- Parameters ----
N_img      <- 80       # number of distinct images
repeats    <- 10        # repeats per image
N_trials   <- N_img * repeats
T          <- 100      # time bins per trial
grid_size  <- 3        # 3x3 grid

# ---- Stimuli X_{n,i,j} ----
# grayscale intensities in [0,1]
X <- array(runif(N_trials * grid_size * grid_size),
           dim = c(N_trials, grid_size, grid_size))

# ---- Beta_{i,j,t} ----
# Diagonals: strongly positive constant over time
beta_ijt <- array(0, dim = c(grid_size, grid_size, T))
for (i in 1:grid_size) {
  for (j in 1:grid_size) {
    if (i == j) {
      beta_ijt[i, j, ] <- 2.5 # + rnorm(T, mean = 0, sd = 0.05)  # diagonal excitatory
    } else {
      beta_ijt[i, j, ] <- 0#-1.5 + rnorm(T, mean = 0, sd = 0.05) # off-diagonal inhibitory
    }
  }
}

# ---- Per-trial random effect b_n ----
b_n <- rnorm(N_trials, mean = 0, sd = 0.1)

# ---- Baseline and simulate y_{n,t} ~ Poisson(exp(linear)) ----
beta0 <- -2
y <- matrix(0, nrow = N_trials, ncol = T)

for (n in 1:N_trials) {
  for (t in 1:T) {
    lin <- beta0 + b_n[n]
    # add spatial contributions
    for (i in 1:grid_size) {
      for (j in 1:grid_size) {
        lin <- lin + beta_ijt[i, j, t] * X[n, i, j]
      }
    }
    lambda_nt <- exp(lin)
    y[n, t] <- rpois(1, lambda = lambda_nt)
  }
}

# ---- Build dataframe ----
library(dplyr)

df <- expand.grid(
  trial = 1:N_trials,
  time  = 1:T
)
df$y <- as.vector(y)

# ---- Add image_id and rep_id ----
df$image_id <- (df$trial - 1) %/% repeats + 1
df$rep_id   <- (df$trial - 1) %% repeats + 1

# ---- Add flattened stimulus features per trial ----
X_mat <- matrix(NA, nrow = N_trials, ncol = grid_size * grid_size)
colnames(X_mat) <- paste0("x_", rep(1:grid_size, each = grid_size), "_", rep(1:grid_size, times = grid_size))
for (n in 1:N_trials) {
  X_mat[n, ] <- as.vector(t(X[n, , ]))
}

df <- cbind(df, df$trial %>% sapply(function(tr) X_mat[tr, ]) %>% t())
colnames(df)[(ncol(df) - ncol(X_mat) + 1):ncol(df)] <- colnames(X_mat)

# ---- Quick check ----
mean_beta <- apply(beta_ijt, c(1, 2), mean)
print("Mean beta (3x3): diagonals should be large positive:")
print(round(mean_beta, 3))

diag_means <- diag(mean_beta)
offdiag_means <- mean(mean_beta[upper.tri(mean_beta) | lower.tri(mean_beta)])
cat("Diagonal means:", round(diag_means, 3), "\n")
cat("Average off-diagonal mean:", round(offdiag_means, 3), "\n")

cat("Overall mean spike count:", mean(df$y), " (sd:", sd(df$y), ")\n")

# ---- Optional: visualize beta ----
image(1:3, 1:3, mean_beta, main = "Mean β across pixels", xlab = "i", ylab = "j")
contour(1:3, 1:3, mean_beta, add = TRUE)

# ---- Check new columns ----
head(df)


### getting estimated by spde , model recovery


df_real <- df %>%
  rename(
    time_bin_num = time,
    firing_count = y
  )


df <- df_real


library(INLA)
library(Matrix)

# Identify pixel columns
pixel_cols <- grep("^x_\\d+_\\d+$", colnames(df), value = TRUE)

#  Extract pixel coordinates
pixel_coords <- do.call(rbind, lapply(pixel_cols, function(name) {
  x <- as.numeric(sub("^x_(\\d+)_\\d+$", "\\1", name))
  y <- as.numeric(sub("^x_\\d+_(\\d+)$", "\\1", name))
  c(x, y)
}))
colnames(pixel_coords) <- c("x", "y")

# Build SPDE mesh
mesh <- inla.mesh.2d(loc = pixel_coords, max.edge = c(1, 1), cutoff = 0.01)
spde <- inla.spde2.matern(mesh = mesh, alpha = 2)

#  Projector matrix for pixels
A_pix <- inla.spde.make.A(mesh = mesh, loc = pixel_coords)

#  Observation × SPDE design
X_spatial <- as.matrix(df[, pixel_cols]) %*% A_pix

#  INLA stack with rep_id
stack <- inla.stack(
  data = list(y = df$firing_count),
  A = list(1, X_spatial, 1),   # intercept + spatial field + rep_id
  effects = list(
    intercept = rep(1, nrow(df)),
    spatial_field = 1:spde$n.spde,
    rep_id = df$rep_id          # include rep_id here
  ),
  tag = "spatial_stack"
)
# Assuming 'time_col' is the column in df representing time (t)

#  Formula including rep_id as IID random effect
formula <- y ~ 1 +
  f(rep_id, model = "iid") +
  f(spatial_field, model = spde)

#  Fit the model
result <- inla(
  formula,
  data = inla.stack.data(stack),
  family = "poisson",
  control.predictor = list(A = inla.stack.A(stack), compute = TRUE),
  control.compute = list(dic = TRUE, waic = TRUE, config = TRUE),
  control.inla = list(strategy = "simplified.laplace")
)


#  Posterior mean of spatial field ## beta retrival?
field_mean <- result$summary.random$spatial_field$mean

[Helpdesk (Haavard Rue)]

This option has been there quite a while (more than a year), so its there in
recent versions but might require R-4.5. If you're running version previous to
R-4.5, this is the moment to get up to date

On Mon, 2025-09-29 at 07:22 -0700, Alokesh Manna wrote:
> P.S. I also conversed with Dr. Rue. But it seems the current INLA version does
> not support me A.local. I was looking an alternate way to resolve this
> problem.

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