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