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.


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,

#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[1] <- 0
  a1[1] ~ dnorm(0,2)
  d[1] ~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) 

1 Answer 1


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[1] <- 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.

  • $\begingroup$ Thanks, Phil--your package notes have already been a great source for my interest in MIRT. I went ahead and tried this, only now it seems like none of the a/theta parameters have converged on either dimension (as evidenced by all Rhats and traceplots). In general, does rotation of these parameters have an impact on these convergence statistics? Parameter recovery is also off for the a parameters, though the d's look pretty good at a glance. Perhaps I need to look into the SEM/CFA MCMC literature. $\endgroup$
    – metrikrn
    Jul 9, 2015 at 17:52
  • $\begingroup$ 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. $\endgroup$ Jul 9, 2015 at 18:16

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.