# MCMC Estimation of Multidimensional 2PL IRT Model Using JAGS

I'm trying to prepare for some more advanced work involving MIRT models I'll be doing later this year by fitting a very simple multidimensional 2PL model to some simulated data using MCMC methods in JAGS. I'm essentially replicating Bolt & Lall (2003) using the WinBUGS code they provide as a foundation (I've included my code below for those interested). Note: Unlike B&L, I used the mirt package to simulate 2D dichotomous data.

I'm fairly new to this area, and certainly new to estimating multidimensional models using MCMC. The code below runs, but the estimates related to the second dimension (discrimination/ability) are awful (several Rhats for the a's on the second dimension are upwards of 12.00 and pitiful mixing in the traceplots). This is the case even after running 100K+ iterations and playing around with different simulating parameters. The parameters for the first dimension and d parameters all seem to have converged nicely.

I've spent a lot of time trying to read up on this, but can't find good references. My questions are whether this has something to do with the data-generating process (perhaps the parameters I'm using), my constraints (it seems as though I'm following the minimum requirement for identifiability and again, mine replicate those of B&L), rotation (though, this doesn't seem to have anything to do with orphaned chains/non-convergence (?)). Let me know if there's something I'm being ignorant about--I feel like I've hit a wall on such an easy model! Also do let me know if there's anything else I can provide. Thanks.

library(R2jags)
library(mirt)

a <- matrix(c(
0.79, 1.36,
0.93, 1.38,
0.58, 0.38,
0.87, 0.87,
0.83, 0.79,
0.31, 0.99,
0.60, 0.48,
0.60, 0.87,
1.64, 0.15,
1.11, 1.30,
0.53, 0.97,
1.26, 0.39,
2.37, 0.00,
1.17, 1.76,
0.96, 1.26,
0.56, 0.46,
1.17, 0.20,
0.63, 0.26,
1.01, 0.47,
0.81, 0.77),20,2,byrow=T)
d <- matrix(c(-0.90, -1.20, 1.00, -0.97, -1.08,
-1.53, -0.61, -0.60, 1.24, -0.69, -1.31, 0.92,
2.49, -0.06, -0.48, -0.82, 1.11, 0.66, -0.15,
-1.08),ncol=1)

#simulating data
y <- simdata(a,d,1000,"dich")

m2pl.model <- function() {
for(i in 1:N) {
for(j in 1:J) {
y[i,j] ~ dbern(p[i,j])
logit(p[i,j]) <- a1[j] * theta[i,1] + a2[j] * theta[i,2] + d[j]
}
theta[i,1:2] ~ dmnorm(mu[1:2], tau[1:2,1:2])
}

a2 <- 0
a1 ~ dnorm(0,2)
d ~dnorm(0,2)

for(j in 2:J) {
a1[j] ~ dnorm(0,2)
a2[j] ~ dnorm(0,2)
d[j] ~ dnorm(0,2)
}}

N <- nrow(y)
J <- ncol(y)
mu <- c(0,0)
tau <- matrix(c(1,0,0,1),2,2,byrow=T)
m2pl.data <- list("y","N","J","mu","tau")

m2pl.params <- c("a1","a2","d","theta")

m2pl.inits <- function(){list(
a1 = rnorm(J,1,.5),
a2 = c(NA,rnorm(J-1,1,.5)),
d  = rnorm(J,0,.5))}

m2pl.model.par <- jags.parallel(data=m2pl.data, inits = m2pl.inits,
m2pl.params, n.iter = 20000,model.file = m2pl.model, n.burnin = 10000,
n.thin = 15)


If you want the orientation to be the same as the simulated data, be sure to set a 0 in your a object which matches the fixed value during estimation. Otherwise, the raw coefficients will be in a different orientation and therefore not easily comparable upon convergence without some rotation.
So in your case, it looks like you set a2 <- 0, so do the same to first row, second column in a when generating your data. Note that upon convergence you may have to manually multiply columns by -1 to flip the signs if the majority of coefficients are negative (easy to automate), because even the signs of the loadings are not invariant in exploratory IRT models. Hope that helps.
• That's very strange, because setting up what I just described works fine in mirt. It's possible that you need to force one or more parameters to be positive throughout the draws, because they could be jumping back and forth between the positive and negative counterparts which will fit equally well. – philchalmers Jul 9 '15 at 18:16