R语言用for循环替代parSim做多场景仿真仅返回平均结果求助

代码核心问题 #

• foreach循环未提取当前场景i的对应参数，直接使用全局参数向量，计算逻辑错误
• 内层for循环第一次迭代就执行return，仅返回1次重复结果，无法完成全部100次仿真
• 无结果存储容器，每次迭代的结果被覆盖，无法保留全量数据
• family参数错误放置在summary()函数中，应该属于glm()的入参，会导致模型拟合错误

修正后完整代码 #

library(foreach)
library(dplyr)

# 基础参数设置
fail = 0
reps = 100
sampleSize = c(120, 250, 350, 500, 730)
b0 = c(-1.09,-0.3962,-0.51,-0.015,-0.5965)
b1 = c(0, 0, 0, 0, 0)
b2 = c(0.063626, 0.0321, 0.2031, 0.52031, 0.8101)
b3 = c(0.3954, 0.5954, 0.8012, 0.9547, 0.9801)

# 生成所有仿真场景
df.scenarios = expand.grid(sampleSize, b0, b1, b2, b3)
colnames(df.scenarios) = c("sampleSize", "b0", "b1", "b2", "b3")

# 仿真执行，.combine = rbind 自动合并所有场景的结果
all_results <- foreach::foreach (i = 1:nrow(df.scenarios), .combine = rbind)  %do% {
  # 提取当前场景的参数
  cur_n <- df.scenarios$sampleSize[i]
  cur_b0 <- df.scenarios$b0[i]
  cur_b1 <- df.scenarios$b1[i]
  cur_b2 <- df.scenarios$b2[i]
  cur_b3 <- df.scenarios$b3[i]
  
  # 初始化当前场景的结果列表
  scen_results <- vector("list", reps)
  
  for (r in 1:reps)  {
    # 生成模拟数据，使用当前场景的样本量
    sex = rbinom(cur_n, size = 1, p = 0.5)
    treatmentgroup = rbinom(cur_n, size = 1, p = 0.5)
    age_cat  = rbinom(cur_n, size = 1, p = 0.5)
    prob = 1 / (1 + exp(-cur_b0 + cur_b1 * treatmentgroup + cur_b2 * sex + cur_b3 * age_cat))
    mortality = rbinom(cur_n, size = 1, p = prob)
    
    dat = data.frame(mortality, sex, treatmentgroup, age_cat)
    dat <- as.data.frame(lapply(dat, as.numeric))
    dat = dat %>% mutate(id = row_number())
    dat$id <- paste0("S0_", dat$id)
    dat = dat %>%  relocate(id)
    
    # 拟合模型，family参数放在glm中
    or_fit <- glm(mortality ~ treatmentgroup, data = dat, family = binomial(link = "logit"))
    or <- summary(or_fit)
    
    b <- or$coef[2]
    se <- or$coef[2,2]
    t <- or$coef[2,3]/or$coef[2,4]
    p.val <- or$coef[2,4]
    
    OR <- exp(b)
    LCI <- exp(b-1.96*se)
    UCI <- exp(b+1.96*se)
    
    # 存储结果，追加场景和重复标识
    scen_results[[r]] <- data.frame(
      scenario_id = i,
      rep_id = r,
      sampleSize = cur_n,
      b0 = cur_b0,
      b1 = cur_b1,
      b2 = cur_b2,
      b3 = cur_b3,
      estimate = b, 
      std.error = se, 
      z.value = t, 
      p.value = p.val, 
      OR = OR, 
      lower_bound = LCI, 
      upper_bound = UCI, 
      method =  "Unadjusted.logit", 
      fail = fail
    )
  }
  # 合并当前场景的所有重复结果后返回
  do.call(rbind, scen_results)
}

# 查看最终结果维度，应该是 场景数*reps 行
dim(all_results)

如果需要进一步提速，可以把%do%改为%dopar%，提前注册并行后端即可实现多线程运行，仿真效率会大幅提升。

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