R输出是否可靠(特别是IRT包ltm)

机器算法验证 r 项目反应理论
2022-03-21 17:36:27

这似乎是一个愚蠢的问题,但对我来说,获得一些专家建议是至关重要的。

我正在研究印度学生的学习水平。我正在使用 R 软件(ltm包)来计算学生分数(根据 IRT 模型的能力估计)。

数据特点:

学生人数:15,000 试卷:3 种不同的形式 一种形式的问题数量:50

我担心的是 R 输出是否足够可靠以在研究社区中发布结果,或者 R 软件的结果是否会受到质疑,特别是在ltm软件包方面。此外,如果 R 可以在不影响质量的情况下轻松处理如此大的数据。

如果您能分享您的意见以及您可能拥有的任何文档/研究论文/参考资料,我将不胜感激。

4个回答

评估软件质量的一个好方法是使用已知的总体参数执行模拟,并观察这些值的恢复情况。更好的是,将参数恢复与其他已知软件进行比较也是一个好主意,因为如果结果有特殊性,您将大致了解正在发生的事情。

以我的经验,ltm似乎可以很好地处理二分数据集,但并不总是能很好地处理多分模型。ltm例如,考虑使用(版本 1.0-0)和mirt(版本 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

两个包似乎都很好地估计了第一个种子,基于相似的最大似然位置。但是,对于其他随机种子:

# 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

显然,ltm在估计这些模型时,可能会停留在局部最小值上。可以通过运行快速模拟来检查这种情况发生的频率以及推断总体参数的含义:

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)

results对象是一个data.frame包含所有复制的此处。总结结果给出

# 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

从这个非常小的模拟中,我们可以看到ltm可能与总体参数相差很大(通过查看max上面的统计数据尤其普遍)。在比较模型的对数似然中也可以看到同样的情况,理论上应该非常接近。

plot(results$mirt_logLik, results$ltm_logLik, 
     ylab = 'ltm log-lik', xlab = 'mirt log-lik')

对数似然

在大多数情况下,对数似然是可比较的,尽管ltm偶尔有比 高得多的最小值mirt这个问题也在其他地方讨论过(https://groups.google.com/forum/#!topic/mirt-package/uK3W4XAMQ9Q)。因为这似乎在运行该ltm::grm()函数时经常发生,所以我强烈建议为您适合的每个模型使用不同的起始值重新估计,以查看 ML 估计是否一致。或者,选择不具有此属性的替代 IRT 软件。

R 在科学界被广泛用于发表论文。R 将您的数据存储在 RAM 中,因此它要么能够处理您的数据集,要么不能 - 取决于数据和处理是否适合内存 - 没有降级模式,您可以获得结果,但它们更少准确的。(从技术上讲,有些软件包可以让您处理比内存更大的数据集,但使用它们并非易事。)

几乎总是可以选择执行类似任务的包,所以如果您非常关心ltm,您还可以查看其他执行 IRT 的包。在我的机器上快速搜索会显示包MCMCpack、、psychKernSmoothIRT,如果您查看CRAN ,可能还有其他包。使用两个包分析您的数据,以确保答案合理一致。

R 的美妙之处在于它是免费的,所以你可以试试看。而且这些包是免费的,所以你可以尝试多个。如果您的数据集太大,或者您对结果不满意,您只会浪费一点时间。

如果您或其他任何人对 R 中的分析结果有疑问,您/他们总是可以查看源代码以准确了解正在进行的计算。对于任何专有软件,您都必须相信他们所做的事情是正确的。

> 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)

我可能会补充一点,在我最近进行的 IRT 工作中,当所有规范都相同时,我发现 ltm 包中的 IRT 结果 (2PLM) 可以很好地映射到 Mplus 统计包中进行的相同分析,包括模型解指数以及参数估计。