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

bagging AIC values from mixed effects linear regression models in R

我已经通过替换数据集对数据集进行了 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")

}

在不知道你的 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 来测试它。