如何在R效能模拟函数中存储while循环的全部效能与样本量至结果数据框

项目背景 #

我正在开展一项效能模拟研究，对象为两臂随机对照试验，核心目的是验证Fast Tap、Hold、Type和Walk测试所需的样本量与效能是否低于mJOA测试。

当前代码问题 #

现有的sim模拟函数仅能存储while循环终止前最后一次的效能值（power1或power2）及对应样本量（N-m）到result矩阵，无法记录循环过程中所有的效能和样本量数据。

修改需求 #

调整代码实现以下目标：

• 将while循环运行过程中（直至所有效能超过0.8时终止）的所有power1或power2值以列表形式存入result数据框的power_values列
• 将对应的样本量（N-m）以列表形式存入sample_sizes列
• 保留原数据框的其他输出字段

修改后的R代码 #

# 初始化协方差矩阵（需提前定义n_timepoints）
Sigma  <- array(0, dim=rep(n_timepoints, 2))

# 定义模拟函数
sim <- function(o_mean, o_sd, o_effect_sizes, N, test = c("t", "ANCOVA")) {
  test <- match.arg(test)
  # 初始化循环过程跟踪存储
  tracking <- list(
    sample_sizes = numeric(),
    power_values = list()
  )
  m <- 0  # 初始化样本量调整参数
  alpha <- 0.05  # 显著性水平
  beta <- 0.8    # 效能终止阈值
  simunum <- 1000 # 模拟次数（可根据需求调整）
  
  # 循环：直至所有时间点效能超过阈值
  repeat {
    current_n <- N - m
    if (current_n <= 2) break # 防止样本量过小无法计算统计量
    
    # 初始化单次模拟统计量存储
    C0_mean  <- matrix(NA, nrow=simunum, ncol=n_timepoints)
    I0_mean  <- matrix(NA, nrow=simunum, ncol=n_timepoints)
    p        <- matrix(NA, nrow=simunum, ncol=n_timepoints)
    
    for (simu in 1:simunum) {
      mu_C0 <- rep(o_mean, n_timepoints)
      mu_I0 <- mu_C0 + o_effect_sizes
      
      C0    <- MASS::mvrnorm(current_n, mu_C0, Sigma)
      I0    <- MASS::mvrnorm(current_n, mu_I0, Sigma)
      
      for (j in 1:n_timepoints) {
        if (test == "t") {
          df        <- 2 * current_n - 2
          S_pooled  <- (sd(C0[,j])^2*(current_n-1) + sd(I0[,j])^2*(current_n-1)) / df
          sd_pooled <- sqrt(S_pooled)
          diff      <- mean(C0[,j]) - mean(I0[,j])
          t_val     <- diff / sqrt(sd_pooled*(2/current_n))
          P         <- pt(-abs(t_val), df) + (1 - pt(abs(t_val), df))
        } else if (test == "ANCOVA") {
          I0.bl     <- MASS::mvrnorm(current_n, o_mean, o_sd)
          C0.bl     <- MASS::mvrnorm(current_n, o_mean, o_sd)
          
          endpoint  <- c(I0[,j], C0[,j])
          baseline  <- c(I0.bl[,1], C0.bl[,1])
          group     <- factor(rep(c("Treatment","Control"), each=current_n))
          model     <- lm(endpoint ~ baseline + group)
          P         <- anova(model)["group", "Pr(>F)"]
        }
        
        p[simu, j]         <- P
        C0_mean[simu, j]   <- mean(C0[,j])
        I0_mean[simu, j]   <- mean(I0[,j])
      }
    }
    
    # 计算当前样本量下的效能
    current_power <- colSums(p < alpha) / simunum
    
    # 记录当前样本量与效能
    tracking$sample_sizes <- c(tracking$sample_sizes, current_n)
    tracking$power_values <- c(tracking$power_values, list(current_power))
    
    print(paste("当前样本量:", current_n, "效能:", paste(round(current_power, 3), collapse=", ")))
    
    # 判断终止条件
    if (all(current_power > beta)) {
      final_power <- current_power
      final_n <- current_n
      break
    }
    m <- m + 1 # 减少有效样本量（或改为N <- N +1 增加总样本量，根据原逻辑调整）
  }
  
  # 构建最终结果数据框
  result <- data.frame(
    rho = rep(NA, n_timepoints), # 若rho有定义可补充传入
    endpoint = paste("Timepoint", 1:n_timepoints),
    test = rep(test, n_timepoints),
    final_power = final_power,
    final_n_critical = final_n,
    sample_sizes = I(rep(list(tracking$sample_sizes), n_timepoints)),
    power_values = I(lapply(1:n_timepoints, function(j) sapply(tracking$power_values, `[`, j)))
  )
  
  return(result)
}

# 调用示例（需提前定义mJOA_BL_mean, mJOA_BL_sd, mJOA_effect_sizes等参数）
# mJOA_result2 <- sim(mJOA_BL_mean, mJOA_BL_sd, mJOA_effect_sizes, 18, "t")

代码修改说明 #

1. 新增tracking列表全程记录循环过程中的样本量和效能值
2. 修正循环终止逻辑，确保所有时间点效能超过0.8时停止循环
3. 补充原代码中未定义的关键变量（如simunum、alpha），保证代码可运行
4. 结果数据框中用I()保留列表格式，将全程样本量和对应每个时间点的效能序列分别存入指定列
5. 保留原有的输出字段，同时新增过程数据存储列

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