Hello Stan Community! I am new to stan and I am trying to replicate the work of Mitzi Morris on BYM2 from https://mc-stan.org/users/documentation/case-studies/icar_stan.html for a different area. However, I am getting these two errors even when I am using adapt_delta=0.97, max_treedepth=15, warmup=10000, iterations=15000:
1: There were 5000 divergent transitions after warmup. Seehttp://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmupto find out why this is a problem and how to eliminate them.
2. Examine the pairs() plot to diagnose sampling problems
I tried using non centered parameters for ‘betas’, however, it throws me 7 common errors of stan.
I am using 10 covariates and an offset. It would be great if anyone can help me understand the reason for this error and ways to fix it.
Here is a stan program for the BYM2 model retrieved from github repository of Mitzi Morris:
data {
int<lower=0> N;
int<lower=0> N_edges;
int<lower=1, upper=N> node1[N_edges]; // node1[i] adjacent to node2[i]
int<lower=1, upper=N> node2[N_edges]; // and node1[i] < node2[i]
int<lower=0> y[N]; // count outcomes
vector<lower=0>[N] E; // exposure
int<lower=1> K; // num covariates
matrix[N, K] x; // design matrix
real<lower=0> scaling_factor; // scales the variance of the spatial effects
}
transformed data {
vector[N] log_E = log(E);
}
parameters {
real beta0; // intercept
vector[K] betas; // covariates
real<lower=0> sigma; // overall standard deviation
real<lower=0, upper=1> rho; // proportion unstructured vs. spatially structured variance
vector[N] theta; // heterogeneous effects
vector[N] phi; // spatial effects
}
transformed parameters {
vector[N] convolved_re;
// variance of each component should be approximately equal to 1
convolved_re = sqrt(1 - rho) * theta + sqrt(rho / scaling_factor) * phi;
}
model {
y ~ poisson_log(log_E + beta0 + x * betas + convolved_re * sigma); // co-variates
// This is the prior for phi! (up to proportionality)
target += -0.5 * dot_self(phi[node1] - phi[node2]);
beta0 ~ normal(0.0, 0.05);
betas ~ normal(0.0, 0.05);
theta ~ normal(0.0, 1.0);
sigma ~ normal(0, 1.0);
rho ~ beta(0.5, 0.5);
// soft sum-to-zero constraint on phi)
sum(phi) ~ normal(0, 0.001 * N); // equivalent to mean(phi) ~ normal(0,0.001)
}
generated quantities {
real logit_rho = log(rho / (1.0 - rho));
vector[N] eta = log_E + beta0 + x * betas + convolved_re * sigma; // co-variates
vector[N] mu = exp(eta);
}
Here is the same model with non centered parameterization:
data {
int<lower=0> N;
int<lower=0> N_edges;
int<lower=1, upper=N> node1[N_edges]; // node1[i] adjacent to node2[i]
int<lower=1, upper=N> node2[N_edges]; // and node1[i] < node2[i]
int<lower=0> y[N]; // count outcomes
vector<lower=0>[N] E; // exposure
int<lower=1> K; // num covariates
matrix[N, K] x; // design matrix
real<lower=0> scaling_factor; // scales the variance of the spatial effects
}
transformed data {
vector[N] log_E = log(E);
}
parameters {
real beta0; // intercept
vector[K] betas_raw; // covariates
real<lower=0> sigma; // overall standard deviation
real<lower=0, upper=1> rho; // proportion unstructured vs. spatially structured variance
vector[N] theta; // heterogeneous effects
vector[N] phi; // spatial effects
vector[K] gamma;
vector<lower=0>[K] tau;
}
transformed parameters {
vector[N] convolved_re;
vector[K] betas;
// variance of each component should be approximately equal to 1
convolved_re = sqrt(1 - rho) * theta + sqrt(rho / scaling_factor) * phi;
betas= gamma + tau .*betas_raw;
}
model {
y ~ poisson_log(log_E + beta0 + x * betas + convolved_re * sigma); // co-variates
// This is the prior for phi! (up to proportionality)
target += -0.5 * dot_self(phi[node1] - phi[node2]);
gamma ~ normal(0,0.05);
tau ~ cauchy(0,2.5);
beta0 ~ normal(0.0, 0.05);
betas_raw ~ normal(0.0, 0.05);
theta ~ normal(0.0, 1.0);
sigma ~ normal(0, 1.0);
rho ~ beta(0.5, 0.5);
// soft sum-to-zero constraint on phi)
sum(phi) ~ normal(0, 0.001 * N); // equivalent to mean(phi) ~ normal(0,0.001)
}
generated quantities {
real logit_rho = log(rho / (1.0 - rho));
vector[N] eta = log_E + beta0 + x * betas + convolved_re * sigma; // co-variates
vector[N] mu = exp(eta);
}