R语言循环模拟报错及结果聚合问题求助

核心问题点 #

1. 聚合逻辑位置错误：原代码将均值计算放在内层j循环中，每次模拟都会重复计算，且由于xcoef/zcoef是逐步填充的，中间计算的均值并不准确，最终仅保留最后一次循环的错误结果。
2. 结果存储维度不匹配：初始空数据框results直接通过results$xcoef[i]赋值会引发报错，因为数据框没有足够的行数支撑索引赋值。
3. 缺失一类错误率（t1er）计算：原代码未利用int_pvalue存储交互项的p值，也未实现错误拒绝交互项原假设的比例统计。

修正后的代码 #

set.seed(2024)
# 设定条件
b1 <- c(-0.3, 0.3)
b2 <- c(-0.3, 0.3)
conditions <- expand.grid(b1 = b1, b2 = b2)

# 其他参数
sample_size <- 200
nsim <- 10  # 模拟次数
alpha <- 0.05  # 显著性水平，用于计算一类错误率

# 预初始化结果数据框，行数匹配条件数量
results <- data.frame(
  condition_id = 1:nrow(conditions),
  b1_true = conditions$b1,
  b2_true = conditions$b2,
  xcoef_mean = numeric(nrow(conditions)),
  zcoef_mean = numeric(nrow(conditions)),
  t1er = numeric(nrow(conditions))
)

# 嵌套循环执行模拟
for (i in 1:nrow(conditions)){
  # 为当前条件初始化模拟结果存储向量
  xcoef <- numeric(nsim)
  zcoef <- numeric(nsim)
  int_pvalue <- numeric(nsim)
  
  for(j in 1:nsim){
    # 生成真实数据
    sigma <- matrix(c(1,0,0,1), nrow=2)
    XZ <- as.data.frame(MASS::mvrnorm(n=sample_size, mu=c(0,0), Sigma=sigma))
    colnames(XZ) <- c("x", "z")
    XZ$err <- rnorm(n=sample_size, mean=0, sd=1)
    
    # 真实模型无交互项（系数为0），用于计算一类错误率
    XZ$True_y <- conditions$b1[i] * XZ$x + conditions$b2[i] * XZ$z + XZ$err
    
    # 拟合模型并提取结果
    lmodel <- summary(lm(True_y ~ x * z, data=XZ))
    xcoef[j] <- lmodel$coefficients["x", "Estimate"]
    zcoef[j] <- lmodel$coefficients["z", "Estimate"]
    int_pvalue[j] <- lmodel$coefficients["x:z", "Pr(>|t|)"]
  }
  
  # 内层循环结束后，聚合当前条件的所有模拟结果
  results$xcoef_mean[i] <- mean(xcoef)
  results$zcoef_mean[i] <- mean(zcoef)
  # 一类错误率：交互项p值小于alpha的比例（错误拒绝原假设的次数占比）
  results$t1er[i] <- mean(int_pvalue < alpha)
}

# 查看最终结果
print(results)

关键修改说明 #

1. 调整聚合逻辑位置：将均值计算和一类错误率统计移到内层j循环之外，确保所有模拟完成后再进行跨模拟的结果聚合，得到准确的统计量。
2. 预初始化结果数据框：根据条件数量预先设定results的行数，避免维度不匹配报错，同时添加真实系数列方便对比模拟估计值与真实值。
3. 完善一类错误率计算：提取交互项的p值，计算其小于显著性水平alpha的比例——由于真实模型中交互项不存在，拒绝原假设即为错误，该比例就是一类错误率。
4. 简化系数提取语法：使用["x", "Estimate"]的索引方式，比双括号更直观且易维护。

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