Stan throws error in population-averaged longitudinal GR IRT model.

[Shamil Sadigov]

Hi All,

 

I got the graded response IRT model to run (with a lot of help from Bob - thanks a lot) for only one time point. Now I am modelling a longitudinal GR IRT model using a population-averaged approach with a covariance structure on the ability parameters. I have tried the AR(1), compound symmetry and even independence covariances in Stan, but I get errors.

I fixed the mean of the ability parameter at visit 1 to 0, and it's variance to 1 for identifiability, and the discrimination parameter is fixed to be positive.

 

There are 7 variables measured at 4 time points. Each variable is on the 1-7 scale.

 

I get the following error:

 

Error : Error during user-specified initialization:
Error in function stan::prob::normal_log(N4stan5agrad3varE): Random variable is 1.#QNAN:0, but must not be nan!
error occurred during calling the sampler; sampling not done

 

Here is the model with the CS covariance structure.

 

 
 data {
  real m_mu_theta;               // prior for mean of theta
  real var_mu_theta;             // prior for variance of mean of theta
  int<lower=0> T;                // number of categorical time points
  int<lower=0> N;                // number of subjects
  int<lower=0> J;                // number of items
  int<lower=1, upper=7> K[J];    // max level for each item
  int Y[N, J, T];                // data, 3D array of integers
}

parameters {
    ordered[6] kappa[J];                    // intercept : array of ordered vectors
    vector[T] theta[N];                     // ability : theta[N,T] - array of vectors
    real<lower=0.0, upper=10.0> alpha[J];   // discrimination
    real<lower=0.0, upper=1.0> rho;         // Correlation parameter

}

model{
   real sigsq_theta;
   vector[T] mu_theta;
   matrix[T, T] var_theta;
   matrix[T, T] L_theta;
   real p[N, J, T, 7];
   real Q[N, J, T, 6];
 
// prior for mu_theta
   mu_theta[1] <- 0.0;  // fix to 0 for identifiability
   for (t in 2:T){
      mu_theta[t] ~ normal(m_mu_theta, var_mu_theta);
   }

// Compound symmetry structure for the covariance of thetas
   sigsq_theta <- 1.0;  // fix to 1 for identifiability
   var_theta[1,1] <- sigsq_theta;
   for (t in 2:T) {
      var_theta[t,t] <- sigsq_theta;
      for (j in 1:(t-1)) {
        var_theta[t,j] <- sigsq_theta*rho;
        var_theta[j,t] <- var_theta[t,j];
      }
   }
   L_theta <- cholesky_decompose(var_theta);
   theta ~ multi_normal_cholesky(mu_theta, L_theta);

  //  alpha ~ uniform(0.0, 10.0);
   for (j in 1:J) {
      for (k in 1:(K[j] - 1)) {
         kappa[j,k] ~ normal(0.0, 5.0);
      }
   }
  for (t in 1:T){
    for (i in 1:N) {
      for (j in 1:J) {
        for (k in 1:(K[j] - 1)) {
          Q[i, j, t, k] <- inv_logit(kappa[j,k] - alpha[j]*theta[i,t]);
          }
          p[i, j, t, 1] <-  Q[i, j, t, 1];
        for (k in 2:(K[j] - 1)){
          p[i, j, t, k] <- Q[i, j, t, k] - Q[i, j, t, k-1];
        }
        p[i, j, t, K[j]] <- 1 - Q[i, j, t, K[j] - 1];

      // incement log probability directly because Y[i, j]
      // has categorical distribution with varying dimension
        increment_log_prob(log(p[i, j, t, Y[i, j, t]]));
      }
    }
  }
}

[Daniel Lee]

Hi Shamil,

There are two issues with your use of mu_theta (local variable in the model block). I stopped there. Hopefully this will point you in the right direction.

1. Reading the error message.

(FYI, NaN is covered it's in the manual in Appendix B.1., pg 372.)

This error:

Error : Error during user-specified initialization:
Error in function stan::prob::normal_log(N4stan5agrad3varE): Random variable is 1.#QNAN:0, but must not be nan!
error occurred during calling the sampler; sampling not done

comes from this line:

  mu_theta[t] ~ normal(m_mu_theta, var_mu_theta);

The error message is a little hard to parse, but hopefully this will make sense. A sampling statement like "y ~ normal(mu, sigma)" is actually called like "increment_log_prob(normal_log(y, mu, sigma))".

We've called the first argument, the left side of the the sampling statement, the "Random variable." The error message is saying that the left hand side, "mu_theta[t]" in your code, is "1.#QNAN:0".

"1.#QNAN" is how Windows prints NaN. Under Mac or Linux, this is typically "nan", but isn't guaranteed to be anything reasonable.

":0" that's at the end is how Stan prints automatic differentiation variables. I'll skip the details, but that's where it's coming from.

So, if you actually want to use the sampling statement over mu_theta, you have to set all the values of mu_theta. You've only set the first value to 0.0. The rest are NaN. Stan is an imperative language and won't automatically do anything for you between the time of declaring mu_theta and using mu_theta.

2. What should you be doing?

I think what you're trying to do is use mu_theta as a parameter. If this is the case, you should move it to the parameters block.

I think you might have moved this to the model block due to the non-identifiability and trying to set the first value as 0. There are multiple ways to do this. One way is:

within the parameters block, defining a T-1 vector of parameters and leaving mu_vector in the model block. So it'd look like:

parameters {

   vector[T-1] mu_theta_raw;

}

model {

   vector[T] mu_theta;
...

   mu_theta[1] <- 0.0;

   for (t in 1:(T-1))

      mu_theta[t+1] <- mu_theta_raw[t];
...
}

(leave the rest the same)

By the time the model block gets executed, the parameters have already been set to values. Transformed parameters block has been executed.

Hopefully that helps.

Daniel

--
You received this message because you are subscribed to the Google Groups "Stan users mailing list" group.
To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[Shamil Sadigov]

Hi Daniel,

 

Thank you for the solution. The trick works -  the model is now sampling. I have a question about the "mu_theta_raw" parameter - what is its implied pdf according to the code?

 

I have also been getting the following error:

 

 

Warning (non-fatal): Left-hand side of sampling statement (~) contains a non-linear transform of a parameter or local variable. 
You must call increment_log_prob() with the log absolute determinant of the Jacobian of the transform. 
Sampling Statement left-hand-side expression: 
get_base1(mu_theta,t,"mu_theta",1) ~ normal_log(...) 
 

But I am not sure what non-linear transformation I am actually doing. Inverse logit perhaps?

 

-Shamil.

[Bob Carpenter]

Which version of Stan are you using? I thought we'd gotten
rid of that erroneous error message recently. If it's just indexing,
you don't need a Jacobian. It's just a projection.

- Bob

[Shamil Sadigov]

Hi Bob,

 

I am using rstan_2.4.0.

R version 3.1.0 (2014-04-10)
Platform: x86_64-w64-mingw32/x64 (64-bit)

G__~1.EXE (GCC) 4.6.3 20111208 (prerelease)

 

Shamil

[Bob Carpenter]

Ah, I see what's going on. You always get the warning
with a local variable because there's no way to infer all that
can happen to it before it gets to sampling statement.

In your case

parameters {
> vector[T-1] mu_theta_raw;
> }
> model {
> vector[T] mu_theta;
> ...
>
> mu_theta[1] <- 0.0;
> for (t in 1:(T-1))
> mu_theta[t+1] <- mu_theta_raw[t];
> ...
> }

Do you really want to assign mu_theta[1] a prior? It's set to 0.
You should be able to just vectorize and replace

> > mu_theta[t] ~ normal(m_mu_theta, var_mu_theta);

with

mu_theta_raw ~ normal(m_mu_theta, var_mu_theta);

You can still assign to mu_theta. The parameters value are set before
entering the model block, so the above will be well defined.

- Bob

[Shamil Sadigov]

Hi Bob,

 

Would it be better to declare mu_theta in the transformed parameters {} block?

 

I didn't want to assign mu_theta[1] a prior, since it is fixed to 0.

 

I updated the code and ran it -  the mixing of the chains was really poor, and the effective sample size was =2 for most parameters and the R_hats were quite high. The same model using only one time point looked quite good in terms of mixing.

 

Here is the new code:

 

 data {
  real m_mu_theta;               // prior for mean of theta
  real var_mu_theta;             // prior for variance of mean of theta
  int<lower=0> T;                // number of categorical time points
  int<lower=0> N;                // number of subjects
  int<lower=0> J;                // number of items
  int<lower=1, upper=7> K[J];    // max level for each item
  int Y[N, J, T];                // data, 3D array of integers
}

parameters {
    ordered[6] kappa[J];                    // intercept : array of ordered vectors
    vector[T] theta[N];                     // ability : theta[N,T] - array of vectors
    real<lower=0.0, upper=10.0> alpha[J];   // discrimination
    real<lower=0.0, upper=1.0> rho;         // Correlation parameter

    vector[T-1] mu_theta_raw;
}

model{
   real sigsq_theta;
   vector[T] mu_theta;
   matrix[T, T] var_theta;
   matrix[T, T] L_theta;
   real p[N, J, T, 7];
   real Q[N, J, T, 6];
 
// prior for mu_theta
   mu_theta[1] <- 0.0;  // fix to 0 for identifiability

   for (t in 1:(T-1))
      mu_theta[t+1] <- mu_theta_raw[t];

   mu_theta_raw ~ normal(m_mu_theta, var_mu_theta);
  

 

 

Is there any evident mistake in the code or are the data and priors/initial values the cause of poor mixing?

 

Shamil.

[Bob Carpenter]

On Aug 20, 2014, at 6:01 AM, Shamil Sadigov <sha...@gmail.com> wrote:

> Hi Bob,
>

> Would it be better to declare mu_theta in the transformed parameters {} block?
>
> I didn't want to assign mu_theta[1] a prior, since it is fixed to 0.

The only difference between declaring in transformed parameters and
the model is whether you want it saved to the output each iteration.

If mu_theta[1] is fixed to 0 and the prior is constant, then having
it there won't do anything (it's just a constant term that drops out).

Looking at the rest of the code I'd try to figure out if you can
use our ordered-logit rather than building something like it yourself.
We have the gradients built in.

I'm not sure about all the indexing --- always worth checking

You can also vectorize the sampling statements. For example,

> for (j in 1:J) {
> for (k in 1:(K[j] - 1)) {
> kappa[j,k] ~ normal(0.0, 5.0);
> }
> }

should be:

for (j in 1:J)

kappa[j] ~ normal(0,5);

You don't need the .0 --- integers are automatically promoted
to real values in real argument positions.

I assume that the size of kappa[j] is K[j] - 1. If not,
you have an unconstrained variable with an improper prior,
and that can be a problem if the posterior's improper.

For more speedup, you can try the non-centered parameterization
for mu_theta (what we used to call the "Matt trick") --- see the manual chapter
on reparameterization and changes of variables.

The reparameterization will affect the convergence behavior, and is
particularly useful if the data sizes are small. The vectorization
just makes iterations go faster.

- Bob

[Shamil Sadigov]

I have been wondering how to specify kappa[j] to be of varying length K[j]-1 in an array of size J (syntax-wise).

In the likelihood statments I  have the index for kappas go until K[j]-1, but in the parameter declaration, I am not sure how to define an array of length J made of vectors of (varying) length K[j] - 1. Right now the code declares all kappa[j]s to be of size max(K[j] - 1), which is =6. Some priors for kappa then don't have the corresponding likelihood, I can see how this might cause an improper posterior.

 

At different time points the max realized response is different, although the possible max is 7 for all items X1-X7.

 

                                                                                                                                   
    T1 T2 T3 T4 T5 T6 T7                                                                                                                                  
X1  6  6  6  5  5  6  6                                                                                                                                  
X2  6  6  6  6  6  6  6                                                                                                                                  
X3  7  6  6  6  7  6  6                                                                                                                                  
X4  7  6  6  6  6  6  6                                                                                                                                  
X5  6  6  6  7  6  6  6                                                                                                                                  
X6  7  7  6  7  6  6  6                                                                                                                                  
X7  6  6  6  6  6  6  6                                                                                                                                  

 

- I don't want to model a separate intercept kappa_{ijt} for each time point - only for each level of each variable kappa_{ij}. Would that still be possible with varying max level?

 

-Shamil

[Shamil Sadigov]

Update on chain mixing:

 

With more informative priors the mixing and n_eff look quite good now. R_hats are ~=1.

 

  alpha ~ normal(1.0, 3.0);                   
   for (j in 1:J) {                   
         kappa[j] ~ normal(0.0, 2.5);                   
      }            

[Bob Carpenter]

I'm guessing it's mixing better because the prior will ensure
the posterior's proper.

To answer the earlier questions about how to code ragged arrays more
efficiently in Stan (that is, without dummy parameters filling
in the gaps), it's a difficult problem.

You need a vector of the size of the sum of the lengths.

And you can give it one big prior (which is more efficient, to boot),
or you can assign a prior that varies by j in the following loop.

Then you need to keep track of where you're at in the array, starting
from zero, and make a local variable of the size K[j] - 1 in the loop
out of the next variables in the array, then update the array position.
There are vector slicing operations that can make this code a bit more
compact, but it's not going to be pretty.

- Bob