Error evaluating the log probability at the initial value

Modeling
1

Hi all,

I’m receiving the following error in my code for a simple hierarchical model.

Chain 1: Rejecting initial value:
Chain 1: Error evaluating the log probability at the initial value.
Chain 1: Exception: normal_lpdf: Location parameter is nan, but must be finite! (in ‘model1b61b3dabc9_d0ce4f9540376bb49b8354006ddbdb68’ at line 28)

I did look at the topics and googled the message, but didn’t see a response that helped with this problem. Here is the code.

modelString = "
data {
int<lower=1> n; // number of students
int<lower=1> G; // number of schools
int schid[n]; // school indices
vector[n] ASRREA01; // reading outcome variable
vector[n] ASBG04;
vector[n] ACBG03A;

}

parameters {
vector[G] gamma00;
vector[G] gamma01;
vector[G] alpha;
vector[G] beta1;
vector[n] mu_beta1;
vector[G] mu_alpha;
vector<lower=0>[n] sigma_read;
vector<lower=0>[n] sigma_alpha;
vector<lower=0>[G] sigma_beta1;
}

model {
vector[n] mu;

for (i in 1:n) {
ASRREA01[i] ~ normal(mu[i], sigma_read);
mu[i] = alpha[schid[i]] + beta1[schid[i]]*ASBG04[i];
}

for (g in 1: G) {
alpha[g] ~ normal(mu_alpha, sigma_alpha);
// mu_alpha[g] = gamma00 + gamma01*ACBG03A[g];
beta1[g] ~ normal(mu_beta1, sigma_beta1);

}

// Priors
gamma00 ~ normal(300,100);
gamma01 ~ normal(0, 5);
mu_beta1 ~ normal(0, 5);
mu_alpha ~ normal(100,10);
sigma_read ~ cauchy(1,5);
sigma_alpha ~ cauchy(1,5);
sigma_beta1 ~ cauchy(1,5);

}

"

Notice that when take out the comment in the line that reads
mu_alpha[g] = gamma00 + gamma01*ACBG03A[g];

I get the following error message

error in ‘model1b61ba56473_d6caa3c42c5f7321f5f780652f54c304’ at line 34, column 5

32:   for (g in 1: G) {
33:      alpha[g] ~ normal(mu_alpha, sigma_alpha);
34:      mu_alpha[g] = gamma00 + gamma01*ACBG03A[g];
        ^
35:      beta1[g] ~ normal(mu_beta1, sigma_beta1);

PARSER EXPECTED: “}”
Error in stanc(file = file, model_code = model_code, model_name = model_name, :
failed to parse Stan model ‘d6caa3c42c5f7321f5f780652f54c304’ due to the above error.

Thanks in advance,

David

2

Hi David,

The error here is that you’re defining a log-probability conditional on mu_alpha before you define what mu_alpha is. Because you haven’t assigned a value to mu_alpha yet, it is automatically initialised to nan, which causes your errors. You’ll likely have similar issues with mu.

To fix this you just need to swap the order of the declarations:

for (g in 1: G) {
mu_alpha[g] = gamma00 + gamma01*ACBG03A[g];
alpha[g] ~ normal(mu_alpha, sigma_alpha);
beta1[g] ~ normal(mu_beta1, sigma_beta1);

}

Another thing to be aware of is that you’ve both defined a prior for mu_alpha and are assigning values to it. While this is completely legal (see this blog post for more information) it’s important to check that the priors are consistent with the range of values expected from the transformation.

Something else that will cause some issues for you is that you’re mixing different variable sizes in a way that I don’t think you want. For example, the construction of mu_alpha:

mu_alpha[g] = gamma00 + gamma01*ACBG03A[g];

mu_alpha[g] (the g-th element) is a real, but the construction: gamma00 + gamma01*ACBG03A[g] is a vector (i.e., vector + vector*real). You’ll run into a similar issue with some of your priors:

alpha[g] ~ normal(mu_alpha, sigma_alpha);

Where alpha[g] is a real, mu_alpha is a vector of size G, and sigma_alpha is a vector of size n.

Let me know if you have any questions about sizes and indexing and I’ll be happy to help.

Cheers,
Andrew

3

Hi Andrew,

Thanks very much for your assistance. I changed everything to real but I am still getting an error, which I believe is due to my lack of understanding of the order of things. The code I am working with is based on rjags code and there clearly are some differences. Any advice would help here.

Thanks

David

SYNTAX ERROR, MESSAGE(S) FROM PARSER:
Cannot assign to variable outside of declaration block; left-hand-side variable origin=parameter
Illegal statement beginning with non-void expression parsed as
mu_alpha[g]
Not a legal assignment, sampling, or function statement. Note that

  • Assignment statements only allow variables (with optional indexes) on the left;
  • Sampling statements allow arbitrary value-denoting expressions on the left.
  • Functions used as statements must be declared to have void returns

error in ‘model2702fed3ddc_e943ae16db983eb9b320834178fc1a82’ at line 33, column 5

31:      
32:   for (g in 1: G) {
33:      mu_alpha[g] = gamma00 + gamma01*ACBG03A[g];
        ^
34:      alpha[g] ~ normal(mu_alpha, sigma_alpha);

PARSER EXPECTED: “}”
Error in stanc(file = file, model_code = model_code, model_name = model_name, :
failed to parse Stan model ‘e943ae16db983eb9b320834178fc1a82’ due to the above error.

—code—
modelString = "
data {
int<lower=1> n; // number of students
int<lower=1> G; // number of schools
int schid[n]; // school indices
vector[n] ASRREA01; // reading outcome variable
vector[n] ASBG04;
vector[n] ACBG03A;

}

parameters {
real gamma00[G];
real gamma01[G];
real alpha[G];
real beta1[G];
real mu_beta1[G];
real mu_alpha[G];
real<lower=0> sigma_read[n];
real<lower=0> sigma_alpha[G];
real<lower=0>sigma_beta1[G];
}

model {
real mu[n];

for (i in 1:n) {
ASRREA01[i] ~ normal(mu[i], sigma_read);
mu[i] = alpha[schid[i]] + beta1[schid[i]]*ASBG04[i];
}

for (g in 1: G) {
mu_alpha[g] = gamma00 + gamma01*ACBG03A[g];
alpha[g] ~ normal(mu_alpha, sigma_alpha);
beta1[g] ~ normal(mu_beta1, sigma_beta1);

}

// Priors
gamma00 ~ normal(300,100);
gamma01 ~ normal(0, 5);
mu_beta1 ~ normal(0, 5);
mu_alpha ~ normal(100,10);
sigma_read ~ cauchy(1,5);
sigma_alpha ~ cauchy(1,5);
sigma_beta1 ~ cauchy(1,5);

}

"

4

Ah, I should have clarified that sorry. When you want to assign transformations of parameters to a new variable (as well as give those transformations priors), you need to use the transformed parameters block, like so:

transformed parameters {
  vector[G] mu_alpha;

  for(g in 1:G) {
    mu_alpha[g] = gamma00 + gamma01*ACBG03A[g];
  }
}

Also, while you’ve changed the types of some the variables to real, you still have issues with mismatched sizes. This is because real gamma00[G] is now an array of size G, rather than a vector. For more background on the typing in Stan, see this section of the manual.

Would you be able to describe your model, or give the notation for which parameters are at which level? Then I can give you some more concrete direction on the right way to express it in Stan

5

So, this is very basic multilevel model with predictors at both the within school and between school level. The model looks like this

Level 1 (students within schools)
ASRRREA01_(ig) = alpha_g + beta1_g * ASBG04_(ig) + r_(ig): the variance of r_(ig) is sigma_read

Level 2 (between schools)
alpha_(g) = gamm00 + gamm01*ACBG03A + u_alpha_g: The variance of u__alpha(g) is
sigma_alpha

beta1_(g) = mu_beta1 + u_beta1_g: The variance of u_beta1(g) is sigma_beta1

I know this can be set up in a hierarchical Bayesian framework, of course, and I did this successfully in jags, but now just learning the syntax of Stan.

Let me know if you need any more information.

Thanks,

David

6

Hi David,

If I’m reading that correctly, your model is:

ASRRREA01_{ig} = \alpha_g + beta_{1g} * ASBG04_{ig} + r_{ig}
r_{ig} \sim N(0,\sigma_{read})
\alpha_g = \gamma_{00} + \gamma_{01} * ACBG03A_g + u_{\alpha g}
u_{\alpha g} \sim N(0,\sigma_{\alpha})
\beta_g = \mu_{beta} + u_{\beta g}
u_{\beta g} \sim N(0,\sigma_{\beta})

Before I run through the Stan syntax, I also need to point out that you’ll need to change how ACBG03A is being passed. At the moment, it’s a vector of size n:

vector[n] ACBG03A;

Which I’m assuming contains the Level-2 covariate repeated for each individual in the same school. This instead needs to be passed as a vector of size G:

vector[G] ACBG03A;

So that you’re only passing the covariate value once per school.

With writing Stan models, it’s generally best to write “top down”, so you specify the higher levels before the lower levels (which is the opposite for JAGS/BUGS, if I recall correctly). This helps you avoid the issues of using variables before defining them that you ran into earlier.

As for the model, an initial (JAGS-esque) syntax using loops would like so:

data {

  int<lower=1> n; // number of students
  int<lower=1> G; // number of schools
  int schid[n]; // school indices
  vector[n] ASRREA01; // reading outcome variable
  vector[n] ASBG04;
  vector[G] ACBG03A;

}

parameters {
  vector[G] alpha;
  vector[G] beta1;
  real gamma00;
  real gamma01;
  real mu_beta1;
  real<lower=0> sigma_alpha;
  real<lower=0> sigma_beta1;
  real<lower=0> sigma_read;
}

model{

  gamma00 ~ normal(300,100);
  gamma01 ~ normal(0, 5);
  mu_beta1 ~ normal(0, 5);
  sigma_read ~ cauchy(1,5);
  sigma_alpha ~ cauchy(1,5);
  sigma_beta1 ~ cauchy(1,5);
  
  for(g in 1:G) {
    alpha[g] ~ normal(gamma00 + gamma01 * ACBG03A[g], sigma_alpha);
    beta1[g] ~ normal(mu_beta1, sigma_beta1);
  }

  for(i in 1:n) {
    ASRREA01[i] ~ normal(alpha[schid[i]] + beta1[schid[i]] * ASBG04[i], sigma_read);
  }
}

Note that only the alpha and beta1 parameters are given sizes, the other parameters are just declared as real (with no size), since they should only be a single parameter. However, this model’s reliance on looping is inefficient. One of Stan’s strengths is that most functions are vectorised, so we can mostly do away with loops. This allows the c++ on the back-end to take advantage of more efficient matrix arithmetic or to avoid repeating calculations.

The vectorised version of your model would look like so:

data {

  int<lower=1> n; // number of students
  int<lower=1> G; // number of schools
  int schid[n]; // school indices
  vector[n] ASRREA01; // reading outcome variable
  vector[n] ASBG04;
  vector[G] ACBG03A;
}

parameters {
  real gamma00;
  real gamma01;
  real mu_beta1;
  real<lower=0> sigma_alpha;
  real<lower=0> sigma_beta1;
  real<lower=0> sigma_read;
  vector[G] alpha;
  vector[G] beta1;
}

model{

  gamma00 ~ normal(300,100);
  gamma01 ~ normal(0, 5);
  mu_beta1 ~ normal(0, 5);
  sigma_read ~ cauchy(1,5);
  sigma_alpha ~ cauchy(1,5);
  sigma_beta1 ~ cauchy(1,5);
  
  alpha ~ normal(gamma00 + gamma01 * ACBG03A, sigma_alpha);
  beta1 ~ normal(mu_beta1, sigma_beta1);

  ASRREA01 ~ normal(alpha[schid] + beta1[schid] .* ASBG04, sigma_read);
}

The only changes to the construction of alpha and beta1 are the removal of the indexing and enclosing loop declaration. The declaration for ASRREA01 has two changes to note here. First, rather than iterating through the schid indexes, you can simply pass the whole array of indexes (i.e., alpha[schid]) and the result will be a vector of length n with the appropriate alpha for each individual. Next is the use of the elementwise multiplication operator: .* so that each respective element of the two vectors beta1[schid] and ASBG04 are multiplied together

7

Something to also be aware of with hierarchical models, is that you can often run into many divergences or long run times. This is because the hierarchical nature of the parameters can induce some tricky posterior geometry that the sampler will struggle with. For more information on this and how to resolve it using the non-centered parameterisation, see this section of the manual

8

Hi Andrew, sorry I didn’t get back to you. Yes, I just tried the syntax you provided in JAGS-esque form and got the following error (not the same as the last one). Thanks also for the pointers on syntax. David

Error in new_CppObject_xp(fields$.module, fields$.pointer, …) :
Exception: mismatch in dimension declared and found in context; processing stage=data initialization; variable name=ACBG03A; position=0; dims declared=(72); dims found=(896) (in ‘model5fa5a112f48_ea985ca51bf54f9d294f35274780b4bb’ at line 9)

failed to create the sampler; sampling not done

Here is the code, which I copied in directly from what you sent me.

data {

int<lower=1> n; // number of students
int<lower=1> G; // number of schools
int schid[n]; // school indices
vector[n] ASRREA01; // reading outcome variable
vector[n] ASBG04;
vector[G] ACBG03A;

}

parameters {
vector[G] alpha;
vector[G] beta1;
real gamma00;
real gamma01;
real mu_beta1;
real<lower=0> sigma_alpha;
real<lower=0> sigma_beta1;
real<lower=0> sigma_read;
}

model{

gamma00 ~ normal(300,100);
gamma01 ~ normal(0, 5);
mu_beta1 ~ normal(0, 5);
sigma_read ~ cauchy(1,5);
sigma_alpha ~ cauchy(1,5);
sigma_beta1 ~ cauchy(1,5);

for(g in 1:G) {
alpha[g] ~ normal(gamma00 + gamma01 * ACBG03A[g], sigma_alpha);
beta1[g] ~ normal(mu_beta1, sigma_beta1);
}

for(i in 1:n) {
ASRREA01[i] ~ normal(alpha[schid[i]] + beta1[schid[i]] * ASBG04[i], sigma_read);
}
}

9

That error is related to this section here:

What I was saying there was that you need to change the data that is being passed to Stan to match the way that it is being used in the model. Because ACBG03A is a school-level covariate, there should only be one value per school. In other words, you should only pass G values of ACBG03A to Stan.

The error:

variable name=ACBG03A; position=0; dims declared=(72); dims found=(896)

Indicates that the Stan syntax declared that ACBG03A should be length G, but what was passed was a vector of length n