使用foreach和doParallel并行化模拟时遇下标越界错误求助

错误根源分析 #

你遇到的subscript out of bounds错误来自于lmer_function中对t()函数的调用。代码里大量使用t(coefficients(summary(model))["变量名", c("Estimate", "Std. Error", "Pr(>|t|)")])提取模型系数，但当某个变量因完全共线性、方差为0等原因被lmer自动从模型中剔除时，coefficients(summary(model))["变量名", ...]会返回NULL，此时调用t()就会触发下标越界错误。

你仅对z1j、z2j做了存在性检查，但x1ij、x2ij未做检查，一旦这些变量在某次模拟中被剔除，就会触发错误。另外并行环境下的随机数生成可能导致某些模拟场景出现串行时未碰到的共线性问题，因此你无法复现错误。

具体修复步骤 #

1. 统一变量存在性检查逻辑 #

修改lmer_function，对所有要提取的变量（x1ij、x2ij、z1j、z2j）都做存在性判断，避免直接提取不存在的行：

lmer_function <- function(dgp_grid, model){
  # 定义通用系数提取函数
  extract_coef <- function(cc, var_name) {
    if (!var_name %in% rownames(cc)) {
      return(t(data.frame(Estimate = NA, `Std. Error` = NA, `Pr(>|t|)` = NA)))
    } else {
      return(t(cc[var_name, c("Estimate", "Std. Error", "Pr(>|t|)")]))
    }
  }
  
  if (model == 1) { # lmer模型（两个预测变量）
    model_fit <- lmer(yij_2_p ~ x1ij + z1j + (1|nj), data = dgp_grid)
    cc <- coefficients(summary(model_fit))
    
    coefs_x1ij <- extract_coef(cc, "x1ij")
    coefs_z1j <- extract_coef(cc, "z1j")
    
    res_lmer <- cbind(data.frame(coefs_x1ij), data.frame(coefs_z1j))
    colnames(res_lmer) <- c(
      "estimate_x1ij_lmer_model",
      "std.err_x1ij_lmer_model",
      "p_x1ij_lmer_model",
      "estimate_z1j_lmer_model",
      "std.err_z1j_lmer_model",
      "p_z1j_lmer_model"
    )
    return(res_lmer)
    
  } else if (model == 2) { # lmer模型（四个预测变量）
    model_fit <- lmer(yij_4_p ~ x1ij + z1j + x2ij + z2j + (1|nj), data = dgp_grid)  
    cc <- coefficients(summary(model_fit))
    
    coefs_x1ij <- extract_coef(cc, "x1ij")
    coefs_z1j <- extract_coef(cc, "z1j")
    coefs_x2ij <- extract_coef(cc, "x2ij")
    coefs_z2j <- extract_coef(cc, "z2j")
    
    res_lmer <- cbind(
      data.frame(coefs_x1ij),
      data.frame(coefs_z1j),
      data.frame(coefs_x2ij),
      data.frame(coefs_z2j)
    )
    
    colnames(res_lmer) <- c(
      "estimate_x1ij_lmer_model_4p",
      "std.err_x1ij_lmer_model_4p",
      "p_x1ij_lmer_model_4p",       
      "estimate_z1j_lmer_model_4p",
      "std.err_z1j_lmer_model_4p",
      "p_z1j_lmer_model_4p",        
      "estimate_x2ij_lmer_model_4p",
      "std.err_x2ij_lmer_model_4p",
      "p_x2ij_lmer_model_4p",       
      "estimate_z2j_lmer_model_4p",
      "std.err_z2j_lmer_model_4p",
      "p_z2j_lmer_model_4p"
    )
    return(res_lmer) 
    
  } else {
    stop("model not specified")
  }
}

2. 优化并行任务拆分（避免重复计算） #

当前并行代码让每个核跑完整的manyimulations_gamma，相当于重复计算了4次相同的模拟，属于算力浪费。正确做法是把设计网格拆分成多个子集，每个核跑一部分任务：

# 先创建完整设计网格
design_grid <- expand.grid(ni = ni,
                           nj = nj, 
                           RI_sd= RI_sd,
                           gamma02 = gamma02,
                           gamma20 = gamma20,
                           replication = 1:niter
)
design_grid$ID <- seq.int(nrow(design_grid))

# 拆分设计网格为number_cores份
split_design <- split(design_grid, seq_len(number_cores))

# 并行运行
doParallel::registerDoParallel(number_cores)
start <- Sys.time()
temp_result <- foreach (subset = split_design, .packages=c("lme4")) %dopar% {
  manyimulations_gamma(subset)
}
Sys.time()-start
doParallel::stopImplicitCluster()

# 合并结果
res <- do.call(rbind, temp_result)
save(res, file = "sim_res_gamma_test.rda")  

3. 确保并行环境的随机数独立性 #

并行模拟时需设置独立的随机种子，避免不同核的随机数重复：

set.seed(123)
cl <- makeCluster(number_cores, type = "PSOCK")
registerDoParallel(cl)

temp_result <- foreach (subset = split_design, .packages=c("lme4"), .seed=123) %dopar% {
  manyimulations_gamma(subset)
}
stopCluster(cl)

额外注意事项 #

• 检查数据生成逻辑：比如z1j是rep(rbinom(nj,1,0.5), each=ni)，若nj过小可能出现所有z1j为0或1的情况，导致变量方差为0被模型剔除，可调整dgp逻辑或添加判断。
• 添加错误捕获：在one_simulation_gamma中用tryCatch包裹，避免单个模拟失败导致整个并行任务终止：

one_simulation_gamma <- function(ni,nj,gamma02,gamma20,RI_sd) {
  tryCatch({
    dgp_grid <-dgp_pilot_study_gamma(ni,nj,RI_sd,gamma02,gamma20)
    res_lmer_function_model_2 <- lmer_function(dgp_grid, model = 2)
    return(res_lmer_function_model_2)
  }, error = function(e) {
    return(as.data.frame(matrix(NA, nrow=1, ncol=12, 
                                dimnames=list(NULL, c(
                                  "estimate_x1ij_lmer_model_4p",
                                  "std.err_x1ij_lmer_model_4p",
                                  "p_x1ij_lmer_model_4p",       
                                  "estimate_z1j_lmer_model_4p",
                                  "std.err_z1j_lmer_model_4p",
                                  "p_z1j_lmer_model_4p",        
                                  "estimate_x2ij_lmer_model_4p",
                                  "std.err_x2ij_lmer_model_4p",
                                  "p_x2ij_lmer_model_4p",       
                                  "estimate_z2j_lmer_model_4p",
                                  "std.err_z2j_lmer_model_4p",
                                  "p_z2j_lmer_model_4p"
                                )))))
  })
}

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