R语言for循环回退旧值：模拟数据与Omega系数计算异常

问题概述 #

使用simstudy生成多层有序分类数据，通过psych::omega()（调用lavaan）计算单级Omega系数，采用三重循环遍历组内变异值、项目间相关系数rho及多轮迭代时，出现以下问题：

• 测试3-10次迭代时，频繁触发lavaan模型不收敛错误
• 循环会回退至之前的参数组合后提前终止，例如var=1.5 + rho=0.1时，7次成功迭代后，剩余3次自动复用var=1.0 + rho=0.6的参数，随后整个循环停止
• 迭代次数设为目标的1000次时，输出结果全为NA
• 已尝试更新R/RStudio及相关包、添加/移除tryCatch、调整迭代次数，但未解决核心问题

代码核心问题分析 #

1. 结果数据框赋值逻辑错误：
   原代码中results_df[i,] <- (c(i, v, r, SL_Omega_values[1,1]))每次迭代都覆盖数据框的第i行，而非追加新行。当进入不同var/rho组合的循环时，会覆盖之前组的同一行数据，造成"循环回退参数"的假象。
2. 数据类型不统一：
   SL_Omega_values初始为字符型（"Convergence failed"），成功计算时为数值型，直接赋值会导致数据框列类型混乱，最终出现大量NA。
3. 有序分类变量处理不当：
   psych::omega()中设置poly=FALSE，但模拟生成的是有序分类项目，应使用polychoric相关而非Pearson相关，这会增加模型拟合难度与收敛失败概率。
4. 文件写入时机错误：
   每次迭代都写入文件，但因数据框未正确追加，文件内容会重复覆盖，导致结果混乱。

修复后的代码 #

library(simstudy)
library(dirmult)
library(tidyverse)
library(psych)
library(GPArotation)
library(multilevel)
library(lavaan)

# 初始化结果数据框，统一列类型
results_df <- tibble(
  iteration = integer(),
  within_cluster_variance = numeric(),
  inter_item_rho = numeric(),
  SL_Omega_raw = numeric()  # 用NA表示收敛失败，统一数值型
)

set.seed(123)

var_value <- c(0.0, 0.5, 1.0, 1.5, 2.0, 2.5)
rho_value <- c(0.1, 0.2, 0.3, 0.4, 0.5, 0.6)

for(v in var_value){
  for(r in rho_value){
    # 为当前var-rho组合创建临时结果存储
    temp_results <- tibble()
    for(i in 1:5){
      # ---------------------- 生成数据 ----------------------
      # 班级水平数据
      class.level <- defData(varname = "class_zscore", dist = "normal", 
                             formula = 0, variance = 1, id = "Class_ID")
      class.level <- defData(class.level, varname = "Student_Count", dist = "noZeroPoisson", formula = 20)
      class.data <- genData(100, class.level)
      
      # 学生水平数据
      gen.student <- defDataAdd(varname = "student_EXTERNALIZE_score", dist = "normal", 
                                formula = "class_zscore", variance = v)
      dtClass <- genCluster(class.data, "Class_ID", numIndsVar = "Student_Count", level1ID = "Student_ID")
      dtClass <- addColumns(gen.student, dtClass)
      
      # 生成有序分类项目
      baseprobs <- matrix(c(
        0.973, 0.018, .006, 0.003,
        0.829, 0.095, .050, 0.026,
        0.765, 0.115, .069, 0.051,
        0.882, 0.068, .032, 0.018,
        0.717, 0.106, .081, 0.096,
        0.880, 0.062, .038, 0.020,
        0.905, 0.045, .034, 0.016),
        nrow = 7, byrow = TRUE)
      
      student.items <- genOrdCat(dtClass , adjVar = "student_EXTERNALIZE_score", 
                                 baseprobs, prefix = "Item", 
                                 asFactor=FALSE, idname = "Student_ID",
                                 corstr = "cs", rho = r)
      
      ITEMS <- student.items[,6:12]
      
      # ---------------------- 计算Omega系数 ----------------------
      sl_omega <- NA
      tryCatch({
        # 针对有序分类变量，设置poly=TRUE，使用polychoric相关
        sl_omega_fit <- psych::omega(ITEMS, nfactors = 1, poly = TRUE, plot = FALSE, lavaan = TRUE,
                                     control = list(iter.max = 10000))  # 增加lavaan迭代次数
        sl_omega <- sl_omega_fit$omega.tot
      }, error = function(e) {
        # 收敛失败时保留NA，同时输出错误信息便于排查
        message(sprintf("Iteration %d (var=%.1f, rho=%.1f): %s", i, v, r, e$message))
      })
      
      # 追加当前迭代结果到临时数据框
      temp_results <- bind_rows(temp_results, tibble(
        iteration = i,
        within_cluster_variance = v,
        inter_item_rho = r,
        SL_Omega_raw = sl_omega
      ))
    }
    # 将当前var-rho组合的结果追加到总数据框
    results_df <- bind_rows(results_df, temp_results)
    # 写入当前var-rho组合的结果文件
    write_csv(results_df %>% filter(within_cluster_variance == v, inter_item_rho == r),
              paste0("results_df_", v, "_", r, ".csv"))
  }
}
# 写入所有结果
write_csv(results_df, "results_df_all.csv")

收敛问题额外优化建议 #

• 调整lavaan拟合参数：通过control参数传递更多选项，例如更换优化器optimizer = "nlminb"，或增加迭代次数iter.max = 20000
• 优化项目分布：部分项目的类别概率极端（如第一个项目97.3%集中在第一类），可能导致协方差矩阵奇异，可适当调整baseprobs，或在模拟后检查协方差矩阵是否可逆
• 增加样本量：当前总样本量约2000，当rho较小或组内变异较大时，样本量可能不足以稳定估计，可尝试增加班级数量（如genData(200, class.level)）或班级人数（如formula = 30）
• 考虑多层模型：原数据为多层结构，单级Omega忽略聚类效应可能导致模型拟合偏差，可尝试使用multilevel包的多层可靠性分析函数

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