Question on Modeling and Interpreting a Polytomous Rasch Testlet

[Joseph Miller]

Hi Phil and Everyone, 

I am a bit new to the mirt package and simulation, and I am trying to figure out how to obtain the correct overall estimates for a polytomous Rasch testlet model in mirt-- specifically, the correct overall difficulty estimate and the correct covariances from the model. After doing some searching and correctly simulating and estimating a dichotomous Rasch testlet, I believe I found the correct Polytomous Rasch Testlet model to be: 

I then simulated the data such that: 

set.seed(1234)

N <- 10000 # number of persons
I <- 12 # number of items
K <- 4 # number of response categories (0,1,2,3)
TT <- 3 # number of testlets
items_per_testlet <- I / TT

# Assign items to testlets
testlet_id <- rep(1:TT, each = items_per_testlet)

# Simulate person abilities
theta <- rnorm(N, mean = 0, sd = 1)

# Simulate testlet effects (gamma)
sdt <- seq(0.5, 1.5, length.out = TT)
gamma <- matrix(0, nrow = N, ncol = TT)
for (tt in 1:TT) {
  gamma[, tt] <- rnorm(N, mean = 0, sd = sdt[tt])
}

# Simulate item difficulties
delta <- rnorm(I, mean = 0, sd = 1)

# Simulate thresholds (step parameters)
tau <- matrix(NA, nrow = I, ncol = K-1)
for (i in 1:I) {
  tau[i, ] <- sort(rnorm(K-1, mean = 0, sd = 1))
}

# Set discrimination (a_i = 1 for PCM/Rasch)
a_i <- 1

resp <- matrix(NA, nrow = N, ncol = I)
for (n in 1:N) {
  for (i in 1:I) {
    d <- testlet_id[i]
    logits <- a_i * (theta[n] + gamma[n, d] - delta[i] - c(0, tau[i, ]))
    exp_logits <- exp(logits)
    cum_probs <- exp_logits / sum(exp_logits)
    cum_probs <- cumsum(cum_probs)
    u <- runif(1)
    resp[n, i] <- findInterval(u, cum_probs)
  }
}
colnames(resp) <- paste0("X", 1:I)

I then estimate the model using the following equation: 

mod_z <- mirt::bfactor(resp, testlet_id, itemtype = "Rasch", verbose= F)

coef(mod_z, simplify = T)

Overall, I can obtain the correct tau / step parameters, but I can't seem to get an overall difficulty parameter (unable to find it or calculate it) nor can I get the correct covariances of each testlet (all zeros in coef). 

1. Is there a recommended approach in mirt to obtain an overall test difficulty estimate in this context? or is this being simulated correctly?

2. How can I extract the estimated covariances between the testlet factors from the fitted model? or is this being simulated correctly?

Thanks in advance to any and all help!

Joe