I was able to fit the Longitudinal IRT model in Winbugs for an ordinal response by extending the BUGS code I took from the paper by Curtis in JSS http://www.jstatsoft.org/v36/c01/paper/

However, I am having difficulty introducing a group effect into the code. The data are 7-category ordinal responses coded 1-7, and measured at multiple time points for each subject. I would like to compare 3 groups of subjects on their latent variable (ability parameter). Here is the BUGs code (I am using R2WinBUGS) without the group effect and it runs as expected:

lgrmar1 <- function(){
        #  Graded response longitudinal model
        #  Unconstrained - discrimination parameters alpha freely estimated    
for (t in 1:T){
    for (i in 1:n){
        for (j in 1:p){
            Y[i, j, t] ~ dcat(prob[i, j, t, 1:K[j]])

            for (k in 1:(K[j] - 1)) {
            logit(P[i, j, t, k]) <- kappa[j, k] - alpha[j]*theta[i, t]

            P[i, j, t, K[j]] <- 1.0

        for (j in 1:p){
            prob[i, j, t, 1] <- P[i, j, t, 1]
            for (k in 2:K[j]) {
                prob[i, j, t, k] <- P[i, j, t, k] - P[i, j, t, k - 1]

for (i in 1:n){
    theta[i, 1:T] ~ dmnorm(mu.theta[], Pr.theta[,])

# Prior for mu.theta
mu.theta[1] <- 0.0
for (t in 2:T){
    mu.theta[t] ~ dnorm(m.mu.theta, pr.mu.theta)

pr.mu.theta <- pow(s.mu.theta, -2)

# AR(1) structure for Sigma.theta
sigsq.theta <- 1.0
Sigma.theta[1, 1] <- sigsq.theta
for (t in 2:T){
    Sigma.theta[t, t] <- sigsq.theta
    for (j in 1:(t - 1)){
        Sigma.theta[t, j] <- sigsq.theta*pow(rho, t - j)
        Sigma.theta[j, t] <- Sigma.theta[t, j]

Pr.theta[1:T, 1:T] <- inverse(Sigma.theta[,])
rho ~ dunif(-1.0, 1.0)

# Priors on item parameters
for (j in 1:p){
    alpha[j] ~ dnorm(m.alpha, pr.alpha) %_% I(0, )
pr.alpha <- pow(s.alpha, -2)

# thresholds need to be ranked
for (j in 1:p){
    for (k in 1:K[j] - 1){
        kappa.star [j, k] ~ dnorm(m.kappa, pr.kappa)
        kappa[j, k] <- ranked(kappa.star[j, 1:(K[j] - 1)], k)
pr.kappa <- pow(s.kappa, -2)       

Constants are: 
T <- dim(Y)[3] 
p <- dim(Y)[2] 
sample size 
n <- dim(Y)[1]  
K - max level for ordinal variable for each item 
K <- apply(apply(Y, c(2,3), max), 1, max)

Prior parameters: 
m.alpha <- 1.0 
s.alpha <- 2.5 
m.kappa <- 0.0
s.kappa <- 2.5 
m.mu.theta <- 0 
s.mu.theta <- 1

My question is, how can I code a group effect on the latent variable scale? Should I model the group effect directly in the likelihood statement (say, logit(P[i, j, t, k]) <- kappa[j, k] - alpha[j]*theta[i, t] + coef1*grp1 + coef2*grp2, there are 3 groups), or should the group effect enter though the latent variable theta? For the latter case, what would be the modification to the code?

  • $\begingroup$ I came back to answer my own question. The answer: both ways can be used to model the group effect. Modeling theta allows for an easier interpretation. We just have to fix the distribution of one of the group to, say N(0,1) via a prior for identifiability. For the categorical time models it is more straightforward to model the Theta, and for the random coefficients continuous time model it looks like modeling the latent response is the way to go, since the variance of Theta can be controlled more easily for identifiability. $\endgroup$ – user151310 Nov 19 '15 at 15:51

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.