ML person parameters with fscores seem to be slightly off in a GRM

[philipp...@googlemail.com]

Dear Philip,

I have implemented a kind of NR person parameter estimation scheme and to validate it, I compared it to mirt:fscores. The functions agree for 1 dim 2PL, 1 dim GRM, but unfortunately, there is a slight discrepancy for a 2 dim GRM, and I seem to be able to obtain a slightly better likelihood (2nd decimal place). Below is a code example (excuse the messy code for the log-likelhood).

From what I read in the source code of fscores.internal, it could just be that nlm converges too early, but this is just a wild guess.

Kind wishes,

  Philipp

# log lik of a GRM in a slightly different parametrization as in mirt (the thresholds have a negative sign) 

log_lik <- function(loadings, # D columns of loadings

                    thresh, # matrix of thresholds of same number of rows as loadings. 

                    x, # response pattern, coded from 0, ..., K-1.

                    theta){

  I <- nrow(loadings)

  stopifnot(I == nrow(thresh))

  D <- ncol(loadings)

  thresh <- cbind(rep(-Inf,I), thresh, rep(Inf,I))

  idx <- x +1 # easier for indexing

  delta <- matrix(thresh[cbind(rep(1:I,2),c(idx,idx+1))],ncol = 2)

  

  nu <- matrix(loadings %*% theta, nrow=I,ncol =2) - delta

  

  u <- 1/(1+exp(-nu))

  

  p <-  u[,1]-u[,2]

  

  sum(log(p))

}

# try two dim model (is this too complex a model for 2 dimensons?)

library(mirt)

pmod2 <- mirt(Science, 2)

loadings <- coef(pmod2,simplify = TRUE)$items[,1:2]

thresh <- -coef(pmod2,simplify = TRUE)$items[,3:5] # sign change needed, since mirt has intercept parametrization

X <- as.matrix(Science) - 1 # code from zero to K-1

scores <- fscores(pmod2,"ML",rotate = "none")

log_lik(loadings, thresh, X[2,], scores[2,])

log_lik(loadings, thresh, X[2,], c(0.03681344, 0.31167518)) # slightly better (values obtained with other method)

# optim's value is identical apart from the 8th decimal place

optim(c(0.03681344, 0.31167518),function(theta)log_lik(loadings,thresh,X[2,],theta), lower = c(0,0), upper = c(0.5,0.5), method = "L-BFGS-B")

[Phil Chalmers]

For this example I would agree the nlm() just didn't converge well for the ML estimate, largely because the likelihood surface is so flat near the ML (which of course happens because the model is complex and the likelihood function so non-informative; hence why the SEs are so large). Take a look (I added this to your previous code):

Th <- expand.grid(seq(0,.5, length.out=200), seq(-.7,.7, length.out=200))
LL <- apply(Th, 1, function(theta)log_lik(loadings,thresh,X[2,],theta))
rgl::plot3d(Th[,1], Th[,2], LL)
rgl::spheres3d(scores[2,1], scores[2,2], log_lik(loadings,thresh,X[2,],scores[2,]), col='red', radius = .01)
rgl::spheres3d(0.03681344,0.31167518, log_lik(loadings,thresh,X[2,],c(0.03681344, 0.31167518)), col='blue', radius = .01)

When I checked the convergence info nlm did indeed return a suboptimal code (code = 2 rather than 1), so the optimizer identified a problem. I think in future releases fscores() should print warning messages for problematic estimates automatically so that users are aware of potential convergence issues, and return the convergence information as an attribute to be further inspected. Thanks for bringing this to my attention.

Phil

--
You received this message because you are subscribed to the Google Groups "mirt-package" group.
To unsubscribe from this group and stop receiving emails from it, send an email to mirt-package...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/mirt-package/48b64d29-6158-4614-a293-30ac70168a0dn%40googlegroups.com.