graded irt model

[Teddy Groves]

Hi Stan users,

I've been trying to fit a graded item response model and have some questions. 

I started with the code that David Furr posted here and made the following changes. First, each of my items has the same number of possible grades, so I put the difficulty parameters in a matrix rather than a vector. Second, I haven't tried to fit any discrimination parameters. Finally, I've attempted to add a subject-level predictor.

My data doesn't quite fit the standard irt paradigm as the 'items' represent types of graded actions rather than exam questions. As a result, many items are attempted multiple times by the same subject. I think this situation still fits within the standard IRT framework but I am no expert.

The STAN code I've been using is as follows:

graded_irt_code_no_discrimination = """
data {
    int<lower=1> I;                // number of types of actions
    int<lower=1> J;                // number of subjects
    int<lower=1> N;                // number of actions
    int<lower=1> G;                // number of grades
    int<lower=1,upper=I> ii[N];    // type of action n
    int<lower=1,upper=J> jj[N];    // subject for action n
    vector[J] x;                   // subject-level predictor
    int<lower=1,upper=G+1> y[N];   // grade for action n, i.e. 1 (max loss), ..., G+1(loss 0)
}
parameters {
    vector[I] kappa_one;                // First difficulties
    matrix<lower=0>[I, G-1] delta;      // Difficulty step-changes
    real theta[J];                      // Ability
    real beta_x;                        // coefficient for x
    real<lower=0> sigma_x;
    }
transformed parameters{
    ordered[G] kappa[I];    // difficulties
    for (i in 1:I){
        kappa[i, 1] <- kappa_one[i];
        for (g in 2:G){
            kappa[i, g] <- kappa_one[i] + cumulative_sum(delta[i])[g-1];
        }
    }
}
model {
    theta ~ normal(x * beta_x, sigma_x);
    for (i in 1:I){
        for (g in 1:G-1){
            delta[i, g] ~ normal(0, 5);
        }
    }
    kappa_one ~ normal(0, 5);
    for (n in 1:N) {
        y[n] ~ ordered_logistic(theta[jj[n]], kappa[ii[n]]);
    }
}
"""

Unfortunately I can't share the data I used as it's commercially sensitive. There are about 1600 subjects, 11 grades, 40 items and 18176 responses.

The model runs as written above, though it seems quite slow to converge. Here are the highlights of the fit summary after 300 iterations in five chains, which took my computer about 1700 seconds:

Inference for Stan model: anon_model_e394c55777111f927c5e15a6c0afa0a2.
5 chains, each with iter=300; warmup=150; thin=1; 
post-warmup draws per chain=150, total post-warmup draws=750.

                mean se_mean     sd   2.5%    25%    50%    75%  97.5%  n_eff   Rhat
kappa_one[0]  -13.86    0.24   2.75 -19.49 -15.78 -13.81 -11.84  -8.84    127   1.01
kappa_one[1]  -10.79    0.19   2.45 -16.31 -12.39 -10.49  -8.98  -6.73    163    1.0
...
delta[0,0]      1.39    0.08   1.29   0.04   0.39   1.05   1.96   4.66    251    1.0
delta[1,0]      1.63    0.09   1.43   0.06   0.58   1.29   2.25   5.74    251    1.0
...
theta[0]        0.48    0.02   0.26  -0.03   0.31   0.48   0.64   1.05    251    1.0
theta[1]        0.26    0.02    0.3  -0.33   0.05   0.25   0.45   0.85    251    1.0
...
beta_x          0.37    0.01   0.04    0.3   0.35   0.37    0.4   0.43      9   1.17
sigma_x         0.32  6.3e-3   0.03   0.27    0.3   0.32   0.34   0.37     17   1.32
...
kappa[20,8]    -0.25  7.3e-3   0.07  -0.39   -0.3  -0.25   -0.2  -0.12     95   1.04
kappa[21,8]     0.92  7.9e-3   0.13   0.67   0.84   0.92    1.0   1.19    251   1.02
...
lp__          -2.7e4   27.09 108.34 -2.7e4 -2.7e4 -2.7e4 -2.7e4 -2.7e4     16   1.36

I also got quite a few of the following warning message in the first 30 iterations, and a couple more later:

Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:

validate transformed params: kappa[k0__] is not a valid ordered vector. The element at 2 is -7.65643, but should be greater than the previous element, -7.65643

I'd really like to know:

• Is there a problem with fitting this kind of model to data where the subjects attempt the same item more than once?
• Why is the warning triggered even though the delta parameter is constrained to be greater than or equal to 0?
  
• Is the speed of the model inevitable given the amount of data or is there a way to make it faster? 
• Is there is a neater or more efficient way to obtain the kappa vectors in the transformed parameters block?
  
• Have I coded the subject-level predictor correctly? (I'm worried that I might have inadvertently caused scale or identifiability issues.)
  
• What would be the best way to add item-level predictors?
  

Any help would be greatly appreciated!

[Bob Carpenter]

See inline.

This is going to be super slow as you're continually recomputing
terms in those cumulative sums.

You can just declare kappa as a parameter directly (as
you've declared it) and put a prior on the differences.

for (i in 1:I)

kappa[i, 2:G] - kappa[i, 1:(G - 1)] ~ normal(0, 5);

The Jacobian has a constant diagonal, so you don't need
to adjust for it.

It should also be a lot faster.

> model {
> theta ~ normal(x * beta_x, sigma_x);
> for (i in 1:I){
> for (g in 1:G-1){
> delta[i, g] ~ normal(0, 5);
> }
> }
> kappa_one ~ normal(0, 5);
> for (n in 1:N) {
> y[n] ~ ordered_logistic(theta[jj[n]], kappa[ii[n]]);
> }
> }
> """

The ordered logistic is slow and hasn't been vectorized yet.

The bigger problem for statistical mixing may be that the
ordered logistic isn't identified if you have an intercept,
because it allows all the cutpoints and the intercept to
move together.

>
>
> Unfortunately I can't share the data I used as it's commercially sensitive. There are about 1600 subjects, 11 grades, 40 items and 18176 responses.

It's always good to start with simulated data to make sure
your model's working as intended.

> The model runs as written above, though it seems quite slow to converge. Here are the highlights of the fit summary after 300 iterations in five chains, which took my computer about 1700 seconds:

Was that running in parallel? How many cores do you have?
We all have four-core CPUs and run four chains in parallel.
Doesn't save compute cycles, but it can be a huge win in
wall time.

> Inference for Stan model: anon_model_e394c55777111f927c5e15a6c0afa0a2.
> 5 chains, each with iter=300; warmup=150; thin=1;
> post-warmup draws per chain=150, total post-warmup draws=750.
>
> mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
> kappa_one[0] -13.86 0.24 2.75 -19.49 -15.78 -13.81 -11.84 -8.84 127 1.01
> kappa_one[1] -10.79 0.19 2.45 -16.31 -12.39 -10.49 -8.98 -6.73 163 1.0
> ...
> delta[0,0] 1.39 0.08 1.29 0.04 0.39 1.05 1.96 4.66 251 1.0
> delta[1,0] 1.63 0.09 1.43 0.06 0.58 1.29 2.25 5.74 251 1.0
> ...
> theta[0] 0.48 0.02 0.26 -0.03 0.31 0.48 0.64 1.05 251 1.0
> theta[1] 0.26 0.02 0.3 -0.33 0.05 0.25 0.45 0.85 251 1.0
> ...
> beta_x 0.37 0.01 0.04 0.3 0.35 0.37 0.4 0.43 9 1.17
> sigma_x 0.32 6.3e-3 0.03 0.27 0.3 0.32 0.34 0.37 17 1.32
> ...
> kappa[20,8] -0.25 7.3e-3 0.07 -0.39 -0.3 -0.25 -0.2 -0.12 95 1.04
> kappa[21,8] 0.92 7.9e-3 0.13 0.67 0.84 0.92 1.0 1.19 251 1.02
> ...
> lp__ -2.7e4 27.09 108.34 -2.7e4 -2.7e4 -2.7e4 -2.7e4 -2.7e4 16 1.36

The problem here is thta the model hasn't converged in the beta and sigma parameters.

> I also got quite a few of the following warning message in the first 30 iterations, and a couple more later:
>
> Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
>
> validate transformed params: kappa[k0__] is not a valid ordered vector. The element at 2 is -7.65643, but should be greater than the previous element, -7.65643

normal(0, 5) with truncation below at zero isn't a zero-avoiding
prior---there's a non-zero density at zero. If you really want them
to be non-zero, you should use a lognormal or a gamma of the right shape.

> I'd really like to know:
> • Is there a problem with fitting this kind of model to data where the subjects attempt the same item more than once?

Not in principle, but it makes a probably unwarranted independence
assumption across the multiple tries (maybe the later tries will
be more likely to succeed, or maybe the student has the same problem
all the time).

> • Why is the warning triggered even though the delta parameter is constrained to be greater than or equal to 0?

See above. It's the equal that's the problem.

> • Is the speed of the model inevitable given the amount of data or is there a way to make it faster?

See above.

> • Is there is a neater or more efficient way to obtain the kappa vectors in the transformed parameters block?

See above---just don't do it, or combine terms in
the cumulative sum claculations.

> • Have I coded the subject-level predictor correctly? (I'm worried that I might have inadvertently caused scale or identifiability issues.)

I'm not sure why you're sampling here. I'm not an expert
on these models, though.

> • What would be the best way to add item-level predictors?

Just add them like in any other GLM.

> Any help would be greatly appreciated!

Hope that helps.

- Bob