CAR model does not converge

Interfaces brms
, ,
1

Hi, I’m trying to fit a CAR model, which I’m new to, using brms. The response variable is continious and on a irregular grid and they stem from model predictions. Because of that, the residual distributions are very weird and I have therefore discretized it and played around with a ordinal models.

Here are two plots to illustrate the continious and discretized data:

library(brms)
library(ggplot2)

d <- readr::read_csv("https://raw.githubusercontent.com/maxlindmark/spatial-metabolic-index/main/example_data.csv")

# Plot data
ggplot(d, aes(X, Y, fill = biom_slope_gr_gc)) +
  geom_raster()

ggplot(d, aes(X, Y, fill = biom_slope)) +
  geom_raster() + 
  scale_fill_gradient2()

The ordinal model fits really well, here’s an example:

# Fit non-spatial ordinal model
fitcod <- brm(
  formula = biom_slope_gr_gc ~ temp_slope_sc*mean_temp_ct + mean_log_biomass_ct,
  data = cod,
  chains = 3, 
  cores = 3, 
  iter = 3000,
  init = 0,
  family = cumulative(probit)
  )

# I'm following this link to get randomized quantile residuals for a discrete distribution so that I can plot them in space: https://cran.r-project.org/web/packages/tidybayes/vignettes/tidybayes-residuals.html#probability-residuals-and-quantile-residuals 

res_ns <- cod |>
  ungroup() |>
  add_predicted_draws(fitcod) |>
  mutate(.prediction = ordered(levels(biom_slope_gr_gc)[.prediction], levels = levels(biom_slope_gr_gc))) |>
  summarise(
    p_lower = mean(.prediction < biom_slope_gr_gc),
    p_upper = mean(.prediction <= biom_slope_gr_gc),
    p_residual = runif(1, p_lower, p_upper),
    z_residual = qnorm(p_residual),
    .groups = "drop_last"
  )

# QQ-plot
ggplot(res_ns, aes(sample = z_residual)) +
  geom_qq() +
  geom_abline()

# In space
ggplot(res_ns, aes(X, Y, color = z_residual)) + 
  geom_point(size = 0.6) + 
  scale_color_gradient2()

However, I wonder now if there’s an issue with spatial autocorrelation, which the underlying data and the discretized data show clear signs of. Below is code to add a CAR structure to the model. For speed I filter the data and use a simpler model just to test the CAR part, which I haven’t gotten to converge yet besides for the example in ?car.

d <- cod |> filter(X > 620 & X < 650 & Y > 6120 & Y < 6150)

# Set up distance and neighbourhood matrices
# See ?car
K <- nrow(d)

distance_grid <- d[, c("X", "Y")]

distance <- as.matrix(dist(distance_grid))
W <- array(0, c(K, K))
W[distance == 3] <- 1 # The distance between two adjacent cells is 3

d$grp <- 1:K
rownames(W) <- 1:K

fitcar <- brm(
  #formula = ordered(biom_slope_gr_gc) ~ temp_slope_sc*mean_temp_ct + mean_log_biomass_ct + car(W, gr = grp, type = 'icar'),
  formula = biom_slope ~ temp_slope_sc + car(W, gr = grp, type = 'icar'),
  data = d,
  chains = 4, 
  cores = 4, 
  iter = 4000,
  init = 0,
  data2 = list(W = W),
  #family = cumulative(probit)
  family = gaussian()
)

plot(fitcar)

Basically, doesn’t converge on any front. Any ideas on where to start adressing this?

  • Operating System: Sonoma 14.3
  • brms Version: brms_2.20.4
2

One thing that was suggested to me is to specify more informative priors, does anyone have experience with that?

3

Apologies — I put together the example too hastily and the data and code are not consistent.

Instead of editing the code I’ll provide a new example here, because we’ve made some progress. I got help from a colleague to set tighter priors on the spatial parameters, and now a small example converges nicely. However, if I make this small example a little irregular, by randomly removing points, it doesn’t converge anymore. And so I wonder if that could be a reason for it not working on my main data, which is fairly patchy (see figure in the end).

Here’s the first example (which works) on a subset of the grid that is “complete”

library(dplyr)
library(ggplot2)
library(brms)

# Read data
d <- readr::read_csv("https://raw.githubusercontent.com/maxlindmark/spatial-metabolic-index/main/data/clean/trends_car_cod.csv")

# Filter a part of this grid that is regular, no missing cells etc.
d2 <- d %>% filter(X > 600 & X < 690 & Y > 6160 & Y < 6180)

# Plot in space. We now have a binomial response variable as an example
ggplot(d2, aes(X, Y, color = biom_slope_cat)) +
  geom_point(size = 4)
  
# We can fit a model now with covariates and a CAR component
K <- nrow(d2)
distance_grid <- d2[, c("X", "Y")]
distance <- as.matrix(dist(distance_grid))
W <- array(0, c(K, K))
W[distance == 3] <- 1 # The distance between two adjacent cells is 3

d2$grp <- 1:K
rownames(W) <- 1:K

fitcar <- brm(
  biom_slope_cat ~
    temp_slope_sc*mean_temp_ct + mean_log_biomass_ct +
    car(W, gr = grp, type = 'icar'),
  data = d2,
  chains = 4,
  cores = 4,
  iter = 2000,
  data2 = list(W = W),
  family = bernoulli(),
  prior = c(
    prior(normal(0, 1), class = "sdcar"),
    prior(normal(0, 1), class = "b")
  ),
  control = list(adapt_delta = 0.9, max_treedepth = 12)
)

# Looks great!
plot(fitcar)

Now if we try again with the same model fitted to the same data but with some cells removed to make it irregular, it doesn’t converge.

set.seed(1)

d3 <- d2[sample(nrow(d2), 130), ]

ggplot(d3, aes(X, Y, color = biom_slope_cat)) +
  geom_point(size = 4)

K <- nrow(d3)
distance_grid <- d3[, c("X", "Y")]
distance <- as.matrix(dist(distance_grid))
W <- array(0, c(K, K))
W[distance == 3] <- 1 # The distance between two adjacent cells is 3

d3$grp <- 1:K
rownames(W) <- 1:K

fitcar3 <- brm(
  biom_slope_cat ~
    temp_slope_sc*mean_temp_ct + mean_log_biomass_ct +
    car(W, gr = grp, type = 'icar'),
  data = d3,
  chains = 4,
  cores = 4,
  iter = 2000,
  data2 = list(W = W),
  family = bernoulli(),
  prior = c(
    prior(normal(0, 1), class = "sdcar"),
    prior(normal(0, 1), class = "b")
  ),
  control = list(adapt_delta = 0.9, max_treedepth = 12)
)

# Looks not great!
plot(fitcar3)

I wonder if the grid I have to work is too irregular, even for the ‘icar’ method. Here’s an example of the full data:

grid
grid1616×1465 209 KB

4

hey @max_lindmark

It might help to hear more about your raw data and how the spatial grid was created (or how observations were made or data collected with respect to spatial order). It looks like the raw outcome variable is no longer in the github file you’re sharing

The response variable is continious and on a irregular grid and they stem from model predictions. Because of that, the residual distributions are very weird and I have therefore discretized it and played around with a ordinal models.

I’m wondering if discretizing the outcome is really necessary, or if there is another way to address whatever issues you faced there; and whether you’d consider something different (like a SAR model or proper CAR model, or maybe a geostatistical model). The icar model just might not be the best choice for your data.

One issue with the icar model is that it can’t properly handle disconnected graph structures (except with some modification GitHub - ConnorDonegan/Stan-IAR: Spatial intrinsic autoregressive models in Stan ) and that’s an issue with your example model

This R package might be helpful if you haven’t seen it already, especially if you don’t require the flexibility of brms: Bayesian Spatial Analysis • geostan

1 Like
5

Hi! Sorry, the data are on the repository again!

The grid is a 3*3 km UTM grid, created for predicting biomass density of a fish from a spatiotemporal GLMM. In a secondary analysis, I wanted to explore how grid-level trends in biomass density over time were related to grid-level environmental conditions and the trends of those. The irregularity in the attached figure stems from filtering the grid to only have grid cells that cumulatively make up most of the total biomass (95%), to focus on the key habitats.

Initially we thought about staying in the same spatiotemporal GLMM framework, but the spatial random effects had issues capturing the very fine scale spatial autocorrelation that was present when using the raw trends as a response. Because it was also extremely difficult to find a distribution that would work (best I can describe the residuals for this species is a skewed pyramid… and other species had other unique shapes), we started exploring discretizing the response to “positive or negative biomass trends” or ordinal models for more categories.

Models converge nicely without the spatial structure, so I’m thinking that’s where the issue lies!

Alright, so as long as there are only few missing “cells” here and there it’s fine as long as they are not complete islands? In my case that might be an issue then, because for other species the grid is much more made up of isloated bands along the coast.

I think I’m going to read up and explore my options again. Thank you so much for taking the time to think about this, and for pointing me to further reading! I very much appreciate it.

2 Likes