Car models and non-centred parameterisation

I have been using the following code, which is basically the same as in
https://mc-stan.org/users/documentation/case-studies/mbjoseph-CARStan.html

However, I have been having problems when running it. Specifically:

check_all_diagnostics(fit_stan_car_n)
[1] “n_eff / iter looks reasonable for all parameters”
[1] “Rhat looks reasonable for all parameters”
[1] “0 of 16000 iterations ended with a divergence (0%)”
[1] “4 of 16000 iterations saturated the maximum tree depth of 10 (0.025%)”
[1] " Run again with max_depth set to a larger value to avoid saturation"
[1] “Chain 1: E-FMI = 0.0470903724171197”
[1] “Chain 2: E-FMI = 0.0797831522453168”
[1] “Chain 3: E-FMI = 0.0742560409244954”
[1] “Chain 4: E-FMI = 0.0735366147434497”
[1] " E-FMI below 0.2 indicates you may need to reparameterize your model"

I wonder if these problems could be solved by using a non-centered parameterization.

Has anyone implemented a non-centered parameterization for this model?

Thank you in advance for your help



functions {

/**

* Return the log probability of a proper conditional autoregressive (CAR) prior

* with a sparse representation for the adjacency matrix

*

* @param phi Vector containing the parameters with a CAR prior

* @param tau Precision parameter for the CAR prior (real)

* @param alpha Dependence (usually spatial) parameter for the CAR prior (real)

* @param W_sparse Sparse representation of adjacency matrix (int array)

* @param n Length of phi (int)

* @param W_n Number of adjacent pairs (int)

* @param D_sparse Number of neighbors for each location (vector)

* @param eigen Eigenvalues of D^{-1/2}*W*D^{-1/2} (vector)

*

* @return Log probability density of CAR prior up to additive constant

*/

real sparse_car_lpdf(vector phi, real tau, real alpha, array[,] int W_sparse, vector D_sparse, vector eigen, int n, int W_n) {

row_vector[n] phit_D; // phi' * D

row_vector[n] phit_W; // phi' * W

vector[n] ldet_terms;

phit_D = (phi .* D_sparse)';

phit_W = rep_row_vector(0, n);

for (i in 1 : W_n) {

phit_W[W_sparse[i, 1]] = phit_W[W_sparse[i, 1]] + phi[W_sparse[i, 2]];

phit_W[W_sparse[i, 2]] = phit_W[W_sparse[i, 2]] + phi[W_sparse[i, 1]];

}

for (i in 1 : n) {

ldet_terms[i] = log1m(alpha * eigen[i]);

}

return 0.5 * (n * log(tau) + sum(ldet_terms)

- tau * (phit_D * phi - alpha * (phit_W * phi)));

}

}

data {

real<lower=0> prior_a;

real<lower=0> prior_b;

int<lower=1> n;

array[n] int<lower=0> y;

matrix<lower=0, upper=1>[n, n] W; // adjacency matrix

int W_n; // number of adjacent region pairs

}

transformed data {

array[W_n, 2] int W_sparse; // adjacency pairs

vector[n] D_sparse; // diagonal of D (number of neigbors for each site)

vector[n] eigen; // eigenvalues of invsqrtD * W * invsqrtD

{

// generate sparse representation for W

int counter;

counter = 1;

// loop over upper triangular part of W to identify neighbor pairs

for (i in 1 : (n - 1)) {

for (j in (i + 1) : n) {

if (W[i, j] == 1) {

W_sparse[counter, 1] = i;

W_sparse[counter, 2] = j;

counter = counter + 1;

}

}

}

}

for (i in 1 : n) {

D_sparse[i] = sum(W[i]);

}

{

vector[n] invsqrtD;

for (i in 1 : n) {

invsqrtD[i] = 1 / sqrt(D_sparse[i]);

}

eigen = eigenvalues_sym(quad_form(W, diag_matrix(invsqrtD)));

}

}

parameters {

real<lower=0> tau;

real<lower=0, upper=1> alpha;

vector[n] phi;

}

model {

phi ~ sparse_car(tau, alpha, W_sparse, D_sparse, eigen, n, W_n);

tau ~ gamma(prior_a, prior_b);

alpha ~ uniform(0,1);

y ~ poisson_log(phi);

}

generated quantities{

vector<lower=0>[n] lambda; // To use as hotspot threshold

lambda = exp(phi);

}