3 Level Multilevel Model Specification

Modeling
1

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

2

I think the problem might be that
you have parameter vectors of length N

real a[N]; // page level random effects
vector[N] alpha; // page level random effects

but you’re only initializing the first S or P values.

1 Like
3

Nice catch! Thank you for your response.

I made what I think is the correct update to the Stan model:

real a[P]; // page level random effects
vector[P] alpha; // page level random effects

which are assigned in these components:

// page level
for (i in 1:P){
       alpha[i] <- a[page[i]];
}
a ~ normal(a_mu, sigma_a);

If I’m understanding the syntax correctly, alpha and a are each P parameters and the loop should estimate all P of them. Let me know if that’s wrong.

However, my model is still failing to initialize. Curiously, the error message now indicates that “Random variable[1] is nan, but must not be nan!”, whereas in the last message, it was Random variable[5] - so it’s possible that fixed one issue.

Is there a way to map the Random variable[index] to the associated parameter in the model? That might help tip me off to where the issue is arising…

4

now the problem is that vector alpha is size P but in the generated quantities block, you’re indexing it as if it were size N. similar problem in the model block - the sampling statement is vectorized but the vectors are of different sizes. alas, this can’t be caught by the compiler - it’s a run-time error.
I think the error message is a PyStan thing - agreed, it’s difficult to interpret.

1 Like
5

Thanks, Mitzi!

Your comments have tipped me off to a couple other issues with the indexing in my model. I’m going to go back and rework a couple things.

Thanks for your help!