CAR model in STAN

Hello everyone,

I am working on a Multinomial spatial classification problem, where we have considered a CAR prior on the regression coefficients \beta s. I have attached the stan code below. While the code is syntactically alright and it runs fine, the chains are not mixing properly. The data has 3 classes. M=336 and N=30K+.

Now I know, there might be several reasons for the chains not mixing, but I feel like I am also not specifying the model correctly. I was going through the STAN manual, and I saw a point about the K vs K-1 specification of the \beta s. However, I feel like the K-1 specification might not be appropriate given the CAR prior. Is there a workaround on this? Any other suggestion to improve the MCMC sampling is also welcome. Thank you in advance.

The CAR prior has been set up as per suggestion in Banerjee,Carlin, and Gelfand.

data {
  int<lower=1> N; // Number of data
  int<lower=1> M; // Number of covariates
  matrix[N, M] X;
  int<lower=1,upper=3> y[N];
  matrix [M,M] DW;     // Diagonal matrix, entries equal the number of neighbors
  matrix [M,M] W;     // Neighborhood matrix
  real lmin;         // minimum eigenvalue of  DW^{-1/2}W DW^{-1/2}
  real lmax;        // maximum eigenvalue of  DW^{-1/2}W DW^{-1/2}
}
transformed data{
 vector[M] zeros=rep_vector(0, M);
}

parameters {
  vector[M] beta;
  real alpha;             
}
transformed parameters{
  matrix[M,M] prec= inverse(DW-alpha*W);
}
model {
  beta ~  multi_normal(zeros,prec);
  alpha ~ uniform(1/lmin,1/lmax);
  y ~ categorical_logit(X*beta);
}

for the ICAR model, we found that the “soft-sum-to-zero” constraint was preferable - the hard sum-to-zero constraint imposes a geometry that’s extremely difficult for the sampler to fit.

have you looked at the Stan case studies on the CAR and ICAR models?

https://mc-stan.org/users/documentation/case-studies/mbjoseph-CARStan.html

https://mc-stan.org/users/documentation/case-studies/icar_stan.html#about-conditional-autoregressive-models

2 Likes

Hi. Thank you for your reply. I have looked into the links you have preferred before we set out to do this. The current model is actually an extension of a previous model, which was a Binomial classification and for that we used a CAR prior exactly of the form and it worked perfectly. So we wanted to try and see if a simple extension into the multinomial case would work or not.