Is R output reliable (specially IRT package ltm) [duplicate]

This may seem like a foolish question, but it is critical for me to have some expert advice on this.

I am doing some research on students' learning levels in India. I am using R software (ltm package) to calculate student scores (ability estimates as per IRT model).

Features of the data:

No of students: 15,000 Test papers: 3 different forms No of questions in one form: 50

My concern is whether R outputs are reliable enough to publish the results in research community or the results from R software can be questioned, specifically with respect to ltm package. Also if R can handle such large data comfortably without compromising the quality.

I would be very grateful if you could share your inputs and any documentation/research paper/reference you may have.

• R is widely accepted in the statistical/scientific community and many papers used it for the statistical analyses. I can't judge if the ltm package is reliable in all cases, but there has been a peer-reviewed paper about it. Just make sure to cite the R version as well as the package (type citation(package="ltm") to see how). The data size in R is limited by your PC's RAM and Rs numerical accuracy is as good as the one of other packages (as far as I know). – COOLSerdash Jul 10 '13 at 11:52
• Reasons to doubt a package could include failure to produce correct answers that are independently predictable. Failure to match results from quite different software also implies that at least one program is wrong. The point is best addressed through specific enquiries, not by highly general questions. – Nick Cox Jul 10 '13 at 12:05
• "Is R output reliable" is far too broad a question; by including R packages across many repositories and sources, we're looking at tens of thousands of packages and hundreds of thousands of functions either in them or standing on their own. The core R functions are usually highly reliable, the heavily used packages on CRAN generally very solid. The less used packages, and moreso, less used parts of less used packages may still be excellent but you can't assume it. A bit of a skeptical eye and some willingness to engage in learning how to do reasonableness checks is important. – Glen_b Apr 27 '14 at 3:18

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.

R is widely used in scientific circles for published papers. R stores your data in RAM, so either it will be able to process your data set or it won't -- depending on whether the data and processing fit in memory -- there is no degraded mode where you get results but they are less accurate. (Technically, there are packages that let you work with larger data sets than fit in memory, but it's not trivial to use them.)

There are almost always choices of packages that do similar tasks, so if you are very concerned about ltm, you can also look into other packages that do IRT. A quick search on my machine brings up packages MCMCpack, psych, and KernSmoothIRT, and there are probably others if you look on CRAN. Analyze your data with two packages to make sure that the answers are in reasonable agreement.

The beauty of R is that it's free, so you can try it and see. And the packages are free, so you can try more than one. If your data set is too large, or if you're not satisfied with the results, you've only lost a bit of time.

If you, or anyone else, has a question about the results of an analysis in R you/they can always look at the source code to see exactly what computations are being made. With any proprietary software you have to take their word that it is doing the correct things.

> library(fortunes)
> fortune(102)

Mingzhai Sun: When you use it [R], since it is written by so many authors, how
do you know that the results are trustable?
Bill Venables: The R engine [...] is pretty well uniformly excellent code but
you have to take my word for that. Actually, you don't. The whole engine is
open source so, if you wish, you can check every line of it. If people were out
to push dodgy software, this is not the way they'd go about it.
-- Mingzhai Sun and Bill Venables
R-help (January 2004)


I might add that in recent IRT work I am undertaking, when all specifications are made the same, I have found IRT results (2PLM) from the ltm package to map on to the same analyses conducted in the Mplus statistical package extremely well, including model solution indices as well as parameter estimates.