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);
}