Hi everyone!

I have been working with @lizzieinvancouver to understand why simple intercept-only models with two grouping factors would result in poorly mixing chains and low neffective values.

This issue came to our attention when working on a more complex model, but in running a similar model on a much smaller dataset we found the same thing.

The data we used for this example was of scores from the 1932 Olympic Games in Lake Placid. OlympicGames_1932.csv (1.8 KB)

The stan code for our intercept-only model is the following:

```
data {
int<lower=0> N;//number of observation
int<lower=0> n_pair;//number of pairs
int<lower=0> n_judge;//number of judges
int<lower=1,upper = n_pair> pair[N];//vector of pairs
int<lower=1,upper = n_judge> judge[N];//vector of judges
vector[N] y;
}
parameters {
vector[n_pair] gamma;
vector[n_judge] delta;
real<lower=0> mu;
real<lower=0> sigma_gamma;
real<lower=0> sigma_delta;
real<lower=0> sigma_y;
}
transformed parameters {
vector[N] y_hat;
for (i in 1:N)
y_hat[i] = mu + gamma[pair[i]] + delta[judge[i]];
}
model {
mu ~ normal(1,10);
sigma_gamma ~ normal(1,10);
sigma_delta ~ normal(1,10);
sigma_y ~ normal(1,1);
gamma ~ normal(0, sigma_gamma);
delta ~ normal(0, sigma_delta);
//likelihood
y ~ normal(y_hat, sigma_y);
}
generated quantities {
real y_rep[N];
for(n in 1:N)
y_rep[n] = normal_rng(y_hat[n],sigma_y);
}
```

The R code to run the model is:

```
prog <- subset(dat, criterion == "Program")
prog_data <- list(y = prog$score,
N = dim(prog)[1],
pair = prog$pair,
judge = prog$judge,
n_pair = length(unique(prog$pair)),
n_judge = length(unique(prog$judge)))
mdl<- stan('stan/OlympicGames_1932mdl.stan',
data = prog_data)
```

Even with this simple example, we observe poorly miking chains and neffective values as low as 465:

Any thoughts on what processes are underlying this result and how they might be addressed would be much appreciated.

Thanks!