Problems with non-centered variance parameters (Matt tick)

Modeling
1

Hi folks,

I know that this is a much-discussed issue and I have been reading through a large number of posts, but after trying to reparameterize and test my model for over a week I would appreciate some help to see where I am going wrong.

In brief: I am looking to estimate a hierarchical model with four levels–repeated observations within individuals in a given year (on a unit scale), within individuals over the years, within geographical districts (called woredas). The n_eff on the variance parameters is low when I write the model. So far so familiar. However, when I try to reparameterize this, the problem actually becomes worse, with n_eff on the critical parameters cut roughly in half. Here are my questions:

  1. may the issue stem from the fact that I currently do not include the means in the added declarations in the transformed parameter block? Or is there anything else wrong here?

  2. If the means need to be added, how can I do that? I have tried something along the lines of
    rt_w = a + sigma_wrt_w_tilde
    rt_s = rt_s + sigma_w    rt_t_tilde
    rt_t = rt_s + sigma_t*rt_t_tilde
    but then I am left with rt ~ normal(rt_t[id_year],sigma) in the model block, which does not match in dimensions. It seems like I am missing something obvious here, but I just cannot figure it out and it is driving me insane.

  3. if and when I can get this to work, might it be possible/helpful to use a de-centered parametrization also for the residual variance?

Any insights much appreciated!

Here is the original, centered model:

data{
    int<lower=1> N;
    int<lower=1> N_id_year;
    int<lower=1> N_id;
    int<lower=1> N_woreda;
    vector[N] rt;
    int woreda[N];
    int id[N];
    int id_year[N];
}
parameters{
    real<lower=0> tau;
    vector[N_woreda] rt_w;
    vector[N_id] rt_s;
    vector[N_id_year] rt_t;
    real<lower=0,upper=1> a;
    vector[N_id] tau_s;
    vector[N_woreda] tau_w;
    real<lower=0> sigma_w;
    real<lower=0> sigma_s;
    real<lower=0> sigma_t;
    real<lower=0> sigma_tau_w;
    real<lower=0> sigma_tau_s;
}
model{
    vector[N] sigma;
    vector[N] mu;
    tau ~ cauchy( 0.25 , 0.25 );
    sigma_tau_s ~ cauchy( 0.25 , 0.25 );
    sigma_tau_w ~ cauchy( 0.25 , 0.25 );
    sigma_w ~ cauchy( 0.25 , 0.25 );
    sigma_s ~ cauchy( 0.25 , 0.25 );
    sigma_t ~ cauchy( 0.25 , 0.25 );
    tau_w ~ normal( 0 , sigma_tau_w );
    tau_s ~ normal( 0 , sigma_tau_s );
    rt_w ~ normal( 0 , sigma_w );
    rt_s ~ normal( 0 , sigma_s );
    rt_t ~ normal( 0 , sigma_t );
    a ~ normal( 0.5 , 0.25 );
    sigma = tau + tau_s[id] + tau_w[woreda];
    mu = a + rt_w[woreda] + rt_s[id] + rt_t[id_year];
    rt ~ normal( mu , sigma );
}
generated quantities{
    vector[N] sigma;
    vector[N] mu;
    real dev;
    dev = 0;
    sigma = tau + tau_w[woreda] + tau_s[id];
    mu = a + rt_w[woreda] + rt_s[id] + rt_t[id_year];
    dev = dev + (-2)*normal_lpdf( rt | mu , sigma );
}

And here is the non-centered version I have been trying:

data{
    int<lower=1> N;
    int<lower=1> N_id_year;
    int<lower=1> N_id;
    int<lower=1> N_woreda;
    vector[N] rt;
    int woreda[N];
    int id[N];
    int id_year[N];
}
parameters{
    real<lower=0> tau;
    vector[N_woreda] rt_w_tilde;
    vector[N_id] rt_s_tilde;
    vector[N_id_year] rt_t_tilde;
    real<lower=0,upper=1> a;
    vector[N_id] tau_s;
    vector[N_woreda] tau_w;
    real<lower=0> sigma_w;
    real<lower=0> sigma_s;
    real<lower=0> sigma_t;
    real<lower=0> sigma_tau_w;
    real<lower=0> sigma_tau_s;
}
transformed parameters{
    vector[N_woreda] rt_w;
    vector[N_id] rt_s;
    vector[N_id_year] rt_t;
    rt_w = sigma_w*rt_w_tilde;
    rt_s = sigma_s*rt_s_tilde;
    rt_t = sigma_t*rt_t_tilde;
}
model{
    vector[N] sigma;
    vector[N] mu;
    tau ~ cauchy( 0.25 , 0.25 );
    sigma_tau_s ~ cauchy( 0.25 , 0.25 );
    sigma_tau_w ~ cauchy( 0.25 , 0.25 );
    sigma_w ~ cauchy( 0.25 , 0.25 );
    sigma_s ~ cauchy( 0.25 , 0.25 );
    sigma_t ~ cauchy( 0.25 , 0.25 );
    tau_w ~ normal( 0 , sigma_tau_w );
    tau_s ~ normal( 0 , sigma_tau_s );
    rt_w_tilde ~ normal( 0 , 1 );
    rt_s_tilde ~ normal( 0 , 1 );
    rt_t_tilde ~ normal( 0 , 1 );
    a ~ normal( 0.5 , 0.25 );
    sigma = tau + tau_s[id] + tau_w[woreda];
    mu = a + rt_w[woreda] + rt_s[id] + rt_t[id_year];
    rt ~ normal( mu , sigma );
}
generated quantities{
    vector[N] sigma;
    vector[N] mu;
    real dev;
    dev = 0;
    sigma = tau + tau_w[woreda] + tau_s[id];
    mu = a + rt_w[woreda] + rt_s[id] + rt_t[id_year];
    dev = dev + (-2)*normal_lpdf( rt | mu , sigma );
}