Graded Response IRT model initialization failed error

[Shamil Sadigov]

Dear Stan Developers,

 

I am trying to fit a 2PL longitudinal graded reponse IRT model in Stan. I started off with just the GRM IRT part, without the longitudinal component. I am getting an error related to the initial value generation. I tightened the bounds and variances on the priors, but still get the same error (Initialization failed after 100 attempts....).

 

The model has 7 ordinal variables, with a variable maxium level for each dimension.

 

data {
  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];                                  // data
}

parameters {
    real kappa[J, 6];                                       // intercept
    real theta[N];                                             // ability
    real<lower=0.0, upper=2.0> alpha[J];    // discrimination

}

model{
   real p[N, J, 7];
   real Q[N, J, 6];
   theta ~ normal(0.0, 1.0);
  //  alpha ~ uniform(0.0, 2.0);
   for (j in 1:J) {
      for (k in 1:(K[j] - 1)) {
         kappa[j,k] ~ normal(0.0, 2.0);
      }
   }
   
  for (i in 1:N) {
    for (j in 1:J) {
      for (k in 1:(K[j] - 1)) {
        Q[i, j, k] <- inv_logit(kappa[j,k] - alpha[j]*theta[i]);
        }
        p[i, j, 1] <-  Q[i, j, 1];
      for (k in 2:(K[j] - 1)){
        p[i, j, k] <- Q[i, j, k] - Q[i, j, k-1];
      }
      p[i, j, K[j]] <- 1 - Q[i, j, K[j] - 1];

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

 

I have attached the simulated data and the error message.

 

Data set name is: Y

Input parameters are:

# n.variables   
J <- dim(Y)[2]
# sample size
N <- dim(Y)[1]
# K - max level for ordinal variable for each item
K <- apply(Y, 2, max)

 

Data are in csv format.

 

 

 

[Marco Inacio]

You are using random  init = "random" (the default) and Stan is initilizing your parameter alpha (which has constraints) with invalid values. Try manually setting the initial values of your parameters using the init argument (run ?rstan::stan to see an example).

--
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.

[Bob Carpenter]

You didn't give us a complete program for simulation. The
data is a CSV file with headers "V1","V2","V3","V4","V5","V6","V7".

Taking a quick look at your program, though, I don't see anything
that guarantees p[i,j,k] is positive. You have:

> for (k in 1:(K[j] - 1)) {
> Q[i, j, k] <- inv_logit(kappa[j,k] - alpha[j]*theta[i]);
> }
> p[i, j, 1] <- Q[i, j, 1];
> for (k in 2:(K[j] - 1)){
> p[i, j, k] <- Q[i, j, k] - Q[i, j, k-1];
> }

Can't those differences result in negative p? There are no
constraints on kappa, alpha, and theta.

- Bob

> --
> 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.

> <simdata.txt><stan_error1.txt>

[Shamil Sadigov]

Hi All,

 

Thank you for looking into this. The threshold parameter kappa does indeed need a constrain - kappas need to be ordered to ensure the probabilities are positive. I didn't know how to code it in Stan at the time. I managed to code the ordering of the thresholds and now the program works. .

Alpha has a constrain to be positive. I think putting a lower bound of 0 in the parameter declaration takes care of it.

I posted the code below.

 

 
 data {
  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];                   // data
}

parameters {
    matrix[J, 6] kappastar;                 // intercept

    real theta[N];                          // ability

    real<lower=0.0, upper=10.0> alpha[J];   // discrimination

}

model{
   matrix[J, 6] kappa;

   real p[N, J, 7];
   real Q[N, J, 6];
   theta ~ normal(0.0, 1.0);

  //  alpha ~ uniform(0.0, 10.0);

   for (j in 1:J) {
      for (k in 1:(K[j] - 1)) {

         kappastar[j,k] ~ normal(0.0, 50.0);
      }
 kappa[j] <- sort_asc(kappastar[j]);     

[Bob Carpenter]

A better solution might be:

ordered[6] kappastar[J]

which makes J an array of ordered vectors. That is,
kappastar[j] will be an ordered vector of size 6.

You can use the same priors, but normal(0,50) won't do
much beyond just using an improper uniform prior.

You should be careful about identifiability with this model.

We have an ordinal logistic density built in, which might save
you all the index fiddling.

- Bob