如何在R中对PoLCA潜类别模型执行bootstrap似然比检验

报错原因 #

boot.LR()是mirt包专属的似然比检验工具，仅支持mirt包输出的SingleGroupClass类模型对象。你传入的是poLCA包生成的poLCA类模型，不符合mirt的输入格式要求，因此触发了校验报错。

解决方法 #

你可以选择两种方案解决：

方案1：手动编写适配poLCA的bootstrap似然比检验代码 #

直接自己实现检验逻辑，不需要跨包调用，示例代码如下：

# 定义函数：比较k类和k+1类poLCA模型的bootstrap LR检验
boot_lr_poLCA <- function(mod_k, mod_k1, nboot = 200) {
  # 计算观测值的LR统计量
  obs_lr <- -2 * (mod_k$llik - mod_k1$llik)
  # 自由度：两个模型的参数量之差
  df <- mod_k1$npar - mod_k$npar
  boot_lr <- numeric(nboot)
  
  for (i in 1:nboot) {
    # 从简约模型（k类）生成bootstrap模拟数据
    sim_dat <- poLCA.simdata(N = mod_k$Nobs, probs = mod_k$probs, 
                            P = mod_k$P, x = mod_k$x)
    # 用k类和k+1类分别拟合模拟数据
    fit_k <- poLCA(mod_k$formula, data = sim_dat$dat, nclass = mod_k$nclass,
                  maxiter = 10000, nrep = 50, verbose = FALSE)
    fit_k1 <- poLCA(mod_k1$formula, data = sim_dat$dat, nclass = mod_k1$nclass,
                   maxiter = 10000, nrep = 50, verbose = FALSE)
    # 存储模拟LR值
    boot_lr[i] <- -2 * (fit_k$llik - fit_k1$llik)
    cat("完成第", i, "次bootstrap\n")
  }
  # 计算p值：观测LR大于模拟LR的比例
  p_val <- mean(boot_lr >= obs_lr, na.rm = TRUE)
  return(list(obs_lr = obs_lr, df = df, boot_lr = boot_lr, p_value = p_val))
}

# 调用示例，比较1类和2类模型
set.seed(421)
lr_test_result <- boot_lr_poLCA(LCA_1, LCA_2, nboot = 200)
print(lr_test_result$p_value)

方案2：改用mirt包做潜类别分析，直接调用boot.LR() #

如果你不想自己写bootstrap逻辑，可以直接用mirt包的潜类别建模功能，生成的模型可以直接适配boot.LR()，示例代码如下：

library(mirt)
# 拟合1类潜类别模型
mirt_LCA1 <- mirt(data = df, model = 1, itemtype = "graded", 
                 mixture = 1, verbose = FALSE)
# 拟合2类潜类别模型
mirt_LCA2 <- mirt(data = df, model = 1, itemtype = "graded", 
                 mixture = 2, verbose = FALSE)
# 直接调用bootstrap似然比检验
set.seed(421)
boot_lr_result <- boot.LR(mirt_LCA1, mirt_LCA2, R = 200)
print(boot_lr_result)

注意事项 #

• bootstrap的迭代次数nboot/R建议至少设置为200，结果才可靠；如果数据量大会耗时较长，可以开启并行加速。
• 拟合潜类别模型时nrep参数要设置足够大，避免陷入局部最优解影响检验结果。

内容的提问来源于stack exchange，提问作者E. Williamson