从 R 中的混合效应线性回归模型中装袋 AIC 值

Question

我已经用替换数据集重新采样了 1000 次，现在想要将三个模型拟合到这 1000 个数据集中的每一个数据集中，并收集它们的 AIC 分数。此过程的最终目标是获得所有模型中每个模型的平均 AIC 分数以及它们的 95% 置信区间。下面的代码有问题，我不知道我在哪里犯了错误。发生的情况是最终矩阵仅包含前几次迭代的 AIC 得分向量（即不是全部 1000）。我在每次迭代中初始化主矩阵或向量的方式是否有错误？或者我的行添加程序有缺陷？或者，如果代码是正确的，它是否与输入该代码的数据集有关？如果是后者，那么当代码读取这些数据集并跳过它们时，为什么我没有收到错误消息？我已经为此苦苦挣扎了好几天，并且非常困惑，因此将不胜感激。

require(lme4)
require(lmerTest)

# initializing an empty matrix for storing each vector of AIC scores from each iteration
# the matrix has width 3 because three models are fitted at each iteration
AIC.scores = data.frame(matrix(, nrow = 0, ncol = 3))

#fit regression models to each of 1000 datasets
for(iter in 1:1000){

  #retrieving the data set, named accordingly, for the current iteration
  data = read.csv(paste("data_set_", iter,".csv", sep=""), header=TRUE)

  #initializing vector of AICs from models in current iteration
  AIC.score = vector(mode="numeric", length=3)

  mod1 = lmer(RT.log ~ crit.var1.log.std +
                       (1|Subject) +
                       (1|Item),
                          data = data,
                          REML=FALSE)

  AIC.score[1] = summary(mod1)$AIC[1]

  mod2 = lmer(RT.log ~ crit.var2.log.std +
                       (1|Subject) +
                       (1|Item),
                          data = data,
                          REML=FALSE)

  AIC.score[2] = summary(mod2)$AIC[1]

  mod3 = lmer(RT.log ~ crit.var3.log.std +
                       (1|Subject) +
                       (1|Item),
                          data = data,
                          REML=FALSE)

  AIC.score[3] = summary(mod3)$AIC[1]

  #adding vector of AICs scores from current iteration to main matrix
  AIC.scores = rbind(AIC.scores, t(AIC.score))

  cat("bagging iteration", iter, "completed!\n")

 }

#renaming column names in AIC score matrix
colnames(AIC.scores) = c("model1", "model2", "model3")

# function for calculating mean AIC and 95% C.I.s for each model across all  iterations
norm.interval = function(data, z=1.96) {
  mean = mean(data)
  variance = var(data)
  sd = sqrt(variance/length(data))
  c(mean, mean - z * sd, mean + z * sd)
}

for (i in 1:3) {

  cat("The mean, lCI, uCI for model", i, "are:", norm.interval(AIC.scores[,i]), "\n")

}

Answer 1

在不知道你的 lmer 模型是什么或你的数据是什么的情况下在黑暗中拍摄。

将所有 data.frame 作为单个列表读取：

all.data <- lapply(paste0("data_set_", 1:1000, ".csv"), read.csv, header=TRUE)

根据上面的数据列表，以文本形式生成 lmer 公式：

all.form <- lapply(paste0("all.data[[", 1:1000, "]]"), function(x) list(
  mod1 = paste0("lmer(RT.log ~ crit.var1.log.std +  (1|Subject) + (1|Item), REML=FALSE, data =", x, ")"),
  mod2 = paste0("lmer(RT.log ~ crit.var2.log.std +  (1|Subject) + (1|Item), REML=FALSE, data =", x, ")"),
  mod3 = paste0("lmer(RT.log ~ crit.var3.log.std +  (1|Subject) + (1|Item), REML=FALSE, data =", x, ")")
  )) 

执行 lmer 公式并将模型写到一个列表中：

all.lmer.mod <- lapply(all.form. function(x) lapply(x, function(y) eval(parse(text=y))))

提取lmer模型的所有AIC值：

all.AIC <- lapply(all.lmer.mod, function(x) lapply(x, AIC))

请注意，如果您有大型 data.frame 并且您有许多主题和项目级别，则该过程将需要很长时间。通过将顶部的1:1000 更改为首先说1:2 来测试它。