Hello fellow STANers,
We’ve been working on a hierarchical model (menstrual cycle lengths womanwise and the woman’s characteristics populationwise) for more than six months. We built the model incrementally and tested a lot of different datasets with different characteristics (to learn what the model could fit best and how). This summer, we were happy with our model’s convergence: the Rhat&ESS looked good and resulting interpretable parameters made sense.
Recently, I reran RStan only to find the same MC did not converge anymore. Since then, I’ve gone back to playing with the model, implementing it in cmdStanR, changing the sampler parameters, etc., without any success. It fails to converge at all (the parameters have Rhat between 1.2 and 4) and seems to be hitting the maximum tree depth in most iterations (75%).

Has this “unexplainable” switch of convergence ever been seen before ? The dataset is exactly the same, the code versioning (git) shows no differences in the scripts and there have been no automatic or manual system updates. . The seed hasn’t been changed either.

The two R wrappers, cmdStanR and Rstan give slightly different results (on smaller, easier, converging datasets with which I’m testing my implementation). I imagine this is to be expected?

Do you have any idea of what I could further try? As I’m trying to replicate an existing model (described in a published article), I cannot reparametrize it further unless there is no other option at all.
Thanks so much for any help, all input welcome!
data {
int < lower = 1 > N; // Number of cycles
int < lower = 1 > W; // Number of women
int < lower = 1 > iwoman[N]; // Table: women to whom the cycle belongs
vector[N] lnC; // Cycle lengths
}
parameters {
real < lower = 0 > half_nu; // Degrees of freedom
vector < lower = 0 >[N] phi; // Precision
real mub; // Population mean of women means
real < lower = 0 > psi; // Population variance of women means
real < lower = 0 > a; // Population parameters for women variances
real < lower = 0 > b;
vector[W] mu; // Women means
vector < lower = 0 >[W] sigmasq; // Women variances
}
transformed parameters{
real invpsi;
invpsi = 1/psi;
}
model {
half_nu ~ gamma(1.0e3, 1.0e3);
phi ~ gamma(half_nu, half_nu);
mub ~ normal(0,sqrt(1000));
psi ~ inv_gamma(1.0e3, 1.0e3);
a ~ gamma(1.0e3, 1.0e3);
b ~ gamma(1.0e3, 1.0e3);
mu ~ normal(mub, sqrt(psi));
sigmasq ~ gamma(a, b);
lnC ~ normal( mu[iwoman], sqrt(sigmasq[iwoman] ./ phi) );
}
fit_CL < CL_fullhierarchy$sample(data = CL_train_data,
iter_warmup = 3000, iter_sampling = 3000, thin = 10,
chains = 4, threads_per_chain = 2, seed = 803214053,
adapt_delta = 0.99, step_size = 0.01,
refresh = 300, output_dir = paste0(diagnosticspath, "pi_cycle1/"))