Error evaluating the log probability at the initial value

Modeling
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