Hello!

I’m a fairly new Stan user and finding it to be an incredible tool with incredible documentation.

I’m doing some website analysis and am trying to build a 3 level random intercept mixed effects model where the first level is the data level, the second level is the web page level and the third level is the website section level. Whenever I try to run my model I get a “RuntimeError: Initialization failed.” message from PyStan. When I look at the output from the command line I get:

“Exception: normal_lpdf: Random variable[5] is nan, but must not be nan! (in ‘unkown file name’ at line 46)

Initialization between (-2, 2) failed after 100 attempts.”

As a Stan newbie, I’m wondering if my model is somehow misspecified or if I should pursue a different line of inquiry. I’ve tried using a non-centered parameterization and alternative initializations, but ran into the same issue. I should note that I ran a standard multilevel model at the page level on the same dataset without any issues.

Here is the Stan model I built:

```
data {
int<lower=0> N; // number of datapoints
int<lower=0> K; // number of predictors
int<lower=0> S; // number of sections
int<lower=0> P; // number of pages
matrix[N, K] x; // covariates
int y[N]; // successes
int trials[N]; // trials
int<lower=1> page[N]; //page indexes
int<lower=1> section[N]; //section indexes
}
parameters {
vector[K] beta; // fixed effects slopes
real<lower=0> sigma_a; // page level variance
real<lower=0> sigma_sect; // section level variance
real sect_mu; // section level mean
}
transformed parameters {
real sect[S]; // section level random effects
real a_mu[S]; // section level random effects
real a[P]; // page level random effects
vector[P] alpha; // page level random effects
// section level
for (i in 1:S){
a_mu[i] <- sect[section[i]];
}
// page level
for (i in 1:P){
alpha[i] <- a[page[i]];
}
}
model {
// section level hyperpriors
sect_mu ~ normal(0, 10);
sigma_sect ~ cauchy(0, 5);
// section level priors
sect ~ normal(sect_mu, sigma_sect);
// page level priors
sigma_a ~ cauchy(0, 2);
a ~ normal(a_mu, sigma_a);
// fixed effects slope priors
beta ~ normal(0, 1);
// likelihood
y ~ binomial(trials, inv_logit(alpha + x * beta));
}
generated quantities {
vector[N] y_tilde;
for (n in 1:N)
y_tilde[n] = binomial_rng(trials[n], inv_logit(alpha[n] + x[n] * beta));
}
```

Any help or guidance would be greatly appreciated!

Thanks,

Eddie