如何使用R语言mirt包固定部分项目的参数？

需求是仅固定原有项目的参数，新增项目的参数允许自由估计。用原有数据拟合模型一切正常，但加入新列后运行报错，提示itemtype specification is not the correct length，同时出现两条警告信息。

复现代码 #

# 模拟数据
set.seed(1)
dat <- matrix(sample(c(0,1), 100, TRUE), ncol = 5)
colnames(dat) <- paste0('item', 1:5)

# 拟合模型并提取项目参数
library(mirt)
mod <- mirt(dat, model = 1, itemtype = '2PL')
param <- mod2values(mod)
param$est <- FALSE

# 原有数据拟合正常
mod2 <- mirt(dat, model = 1, itemtype = '2PL', pars = param)

# 新增列后报错
new.data <- cbind(dat, item6 = sample(c(0,1), 20, TRUE))
mod3 <- mirt(new.data, model = 1, itemtype = '2PL', pars = param)

# 错误信息：itemtype specification is not the correct length
# 警告信息：
# 1: EM cycles terminated after 500 iterations. 
# 2: In itemtype %in% c("2PLNRM", "3PLNRM", "3PLuNRM", "4PLNRM") & K <  :
#   longer object length is not a multiple of shorter object length

问题根源 #

报错的核心是：你手里的param只包含原来5个项目的参数，但new.data新增了第6个项目，此时参数表的长度和数据列数不匹配；另外，虽然itemtype = '2PL'是单个值时mirt通常会自动循环匹配，但搭配pars参数时需要严格对应每个项目的配置，否则就会触发长度不匹配的错误。

解决方法 #

要实现“固定旧项目参数，新项目自由估计”，得手动扩展参数表，把新项目的初始参数加进去，同时设置好哪些参数要固定、哪些要估计：

方法1：通过临时模型获取新项目参数 #

先给单独的新项目拟合一个临时模型，提取它的初始参数，再和旧参数合并：

# 给新项目单独拟合临时模型，获取初始参数
temp_mod <- mirt(new.data[, "item6", drop=FALSE], model = 1, itemtype = '2PL')
temp_param <- mod2values(temp_mod)
temp_param$est <- TRUE  # 新项目参数允许自由估计

# 合并旧参数和新项目参数
new_param <- rbind(param, temp_param)

# 拟合新模型，itemtype指定为和项目数一致的向量（6个项目就6个'2PL'）
mod3 <- mirt(new.data, model = 1, itemtype = rep('2PL', 6), pars = new_param)

# 查看拟合结果
coef(mod3, simplify = TRUE)

方法2：手动构造新项目参数 #

如果不想拟合临时模型，直接手动造新项目的参数行就行（2PL模型每个项目有区分度a1和难度d两个参数）：

# 手动构造新项目的参数行，初始值设区分度1、难度0
new_item_param <- data.frame(
  item = rep("item6", 2),
  name = c("a1", "d"),
  group = "all",
  value = c(1, 0),
  est = TRUE,
  stringsAsFactors = FALSE
)

# 合并参数表
new_param <- rbind(param, new_item_param)

# 拟合模型
mod3 <- mirt(new.data, model = 1, itemtype = rep('2PL', 6), pars = new_param)

关键注意点 #

• 使用pars参数时，参数表的总行数必须和所有项目的总参数数匹配（5个旧项目×2 + 1个新项目×2 = 12行）。
• est列是开关：FALSE表示固定该参数，TRUE表示允许模型估计该参数。
• itemtype最好指定为和项目数一致的向量，避免自动匹配时出问题，所有项目都是2PL的话，用rep('2PL', ncol(new.data))就行。

内容的提问来源于stack exchange，提问作者Julien