A good way to evaluate the quality of software is to perform simulations with known population parameters and observe how well these values can be recovered. Better yet, comparing parameter recovery to other known software is also a good idea since then you will have a general idea of what's happening if there are peculiarities in the results.
In my experience, ltm seems to work pretty well with dichotomous datasets but doesn't always behave nicely with polytomous models. For instance, consider these two randomly simulated datasets which are estimated with both ltm (version 1.0-0) and mirt (version 1.17.1):
library(ltm)
library(mirt)
#this seed agrees just fine...
set.seed(1)
a <- matrix(rlnorm(20, .2, .3))
d <- matrix(rnorm(20*4, 0, 2), 20)
d <- t(apply(d, 1, sort, decreasing=TRUE))
dat <- simdata(a, d, 500, itemtype = 'graded')
ltmmod <- grm(dat, IRT.param=FALSE)
logLik(ltmmod)
'log Lik.' -11827.25 (df=100)
mirtmod <- mirt(dat, 1)
Iteration: 17, Log-Lik: -11826.036, Max-Change: 0.00006
mirtmod@logLik
[1] -11826.04
The first seed seems to be estimated well enough by both packages, resting on similar maximum likelihood locations. However, for other random seeds:
# this one is way different
set.seed(1234)
a <- matrix(rlnorm(20, .2, .3))
d <- matrix(rnorm(20*4, 0, 2), 20)
d <- t(apply(d, 1, sort, decreasing=TRUE))
dat <- simdata(a, d, 500, itemtype = 'graded')
ltmmod <- grm(dat, IRT.param=FALSE)
logLik(ltmmod)
'log Lik.' -13462.7 (df=100)
mirtmod <- mirt(dat, 1)
Iteration: 18, Log-Lik: -11913.724, Max-Change: 0.00002
mirtmod@logLik
[1] -11913.72
Clearly, ltm is getting hung up somewhere when estimating these models, perhaps resting on local minimums. How often this occurs, and what it means for inferring population parameters, can be checked out by running a quick simulation:
library(SimDesign)
# SimFunctions(summarise=FALSE)
mbias <- function(sample, pop) mean(sample - pop)
mRMSE <- function(sample, pop) sqrt(mean((sample - pop)^2))
Design <- data.frame(N=500)
#-------------------------------------------------------------------
Generate <- function(condition, fixed_objects = NULL) {
while(TRUE){
a <- matrix(rlnorm(20, .2, .3))
d <- matrix(rnorm(20*4, 0, 2), 20)
d <- t(apply(d, 1, sort, decreasing=TRUE))
dat <- simdata(a, d, condition$N, itemtype = 'graded')
ncats <- apply(dat, 2, function(x) length(unique(x)))
if(all(ncats == 5)) break
}
ret <- list(data=dat, a=a, d=d)
ret
}
Analyse <- function(condition, dat, fixed_objects = NULL) {
Attach(dat)
ltmmod <- grm(data, IRT.param=FALSE)
mirtmod <- mirt(data, 1, verbose = FALSE)
if(ltmmod$convergence != 0) stop('ltm failed to converge')
if(!extract.mirt(mirtmod, 'converged')) stop('mirt failed to converge')
cfs <- coef(ltmmod)
ltm_as <- cfs[,'beta', drop=FALSE]
ltm_ds <- -cfs[,1:4]
mv <- mod2values(mirtmod)
mirt_as <- mv$value[mv$name == 'a1']
mirt_ds <- matrix(mv$value[mv$name %in% c('d1', 'd2', 'd3', 'd4')], 20,
byrow=TRUE)
ret <- data.frame(ltm_logLik = logLik(ltmmod),
mirt_logLik = logLik(mirtmod),
ltm_bias_a = mbias(ltm_as, a),
mirt_bias_a = mbias(mirt_as, a),
ltm_bias_d = mbias(ltm_ds, d),
mirt_bias_d = mbias(mirt_ds, d),
ltm_RMSE_a = mRMSE(ltm_as, a),
mirt_RMSE_a = mRMSE(mirt_as, a),
ltm_RMSE_d = mRMSE(ltm_ds, d),
mirt_RMSE_d = mRMSE(mirt_ds, d))
ret
}
#-------------------------------------------------------------------
results <- runSimulation(design=Design, replications=100,
generate=Generate, analyse=Analyse,
packages=c('mirt', 'ltm'), parallel=TRUE)
The results object is a data.frame here containing all the replications. Summarising the results gives
# post analysis
round(apply(results[,3:10], 2, function(x)
c(mean=mean(x), sd=sd(x), max=max(x))), 3)
ltm_bias_a mirt_bias_a ltm_bias_d mirt_bias_d
mean 0.051 0.010 -0.018 0.000
sd 0.182 0.044 0.168 0.063
max 0.923 0.138 0.442 0.185
ltm_RMSE_a mirt_RMSE_a ltm_RMSE_d mirt_RMSE_d
mean 0.217 0.128 0.262 0.173
sd 0.228 0.027 0.234 0.033
max 1.118 0.202 1.261 0.286
From this very small simulation we can see that ltm potentially can be very off from the population parameters (it's especially prevalent by looking at the max statistic above). The same is seen in the log-likelihoods comparing the models, which should in theory be very close.
plot(results$mirt_logLik, results$ltm_logLik,
ylab = 'ltm log-lik', xlab = 'mirt log-lik')
For the most part the log-likelihoods are comparable, though ltm occasionally had much higher minimums than mirt. This problem has been discussed elsewhere as well (https://groups.google.com/forum/#!topic/mirt-package/uK3W4XAMQ9Q). Because this appears to happen so often whenever running the ltm::grm() function I would highly recommend re-estimating with different starting values for every model you fit to see if the ML estimates are consistent. Or, choose alternative IRT software which do not have this property.
Ris widely accepted in the statistical/scientific community and many papers used it for the statistical analyses. I can't judge if theltmpackage is reliable in all cases, but there has been a peer-reviewed paper about it. Just make sure to cite theRversion as well as the package (typecitation(package="ltm")to see how). The data size inRis limited by your PC's RAM andRs numerical accuracy is as good as the one of other packages (as far as I know). $\endgroup$