如何在 R 中为包含预测变量的 GLM 生成 AIC 表

Question

我正在使用以下代码运行一系列 glms，在 R Markdown 中运行。我需要 AIC 表包含变量名称（而不仅仅是模型名称）。

# Specify the names of predictor variables
predictor_vars <- c("communitiesinpa", "formalprotection", "attitudewildlife", "competitionlivestock", "fenced", "parkresourcescat", "routinepatrols", "mitigationlivestockbarriers", "mitigationcommunityengagement")

# Create all possible combinations of predictor variables
all_combinations <- unlist(lapply(1:length(predictor_vars), function(x) combn(predictor_vars, x, simplify = FALSE)), recursive = FALSE)

# Create an empty list to store the GLM models
models_list <- list()

# Fit GLMs for all combinations and store them in the models_list
for (i in seq_along(all_combinations)) {
  formula_string <- paste("totalthreatscore ~", paste(all_combinations[[i]], collapse = "+"))
  formula_object <- as.formula(formula_string)
  model <- glm(formula_object, data = glmdata, family = Gamma)
  models_list[[i]] <- model
}

# Print the summary of each model
for (i in seq_along(models_list)) {
  cat("Model with predictors:", paste(all_combinations[[i]], collapse = " + "), "\n")
  print(summary(models_list[[i]]))
  cat("\n")
}

AIC_table <- aictab(models_list)

# Display the AIC table
knitr::kable(AIC_table)

aictab

提供型号（即 1 到 >500），但不包括变量名称。任何人都可以协助修改代码以包含这些吗？

Answer 1

您可以将变量作为

attr

ibute 添加到每个模型。在

mtcars

上进行了演示，代码略有简化。

predictor_vars <- c("cyl", "disp", "hp", "drat", "wt", "qsec", "vs", "am", 
                    "gear", "carb")
outcome_var <- "mpg"

models_list <- lapply(seq_along(predictor_vars), \(m) combn(predictor_vars, m, \(x) {
  do.call('glm', list(fo=reformulate(x, outcome_var), data=quote(mtcars), 
                      family=quote(Gamma))) |> `attr<-`('p.vars', x)
}, simplify=FALSE)) |> unlist(recursive=FALSE)

建议使用 bbmle::AICctab。

AIC_table <- bbmle::AICctab(models_list)
AIC_table <- cbind(as.data.frame(AIC_table),
                   p.vars=sapply(models_list, attr, 'p.vars') |> 
                     sapply(toString))

# Display the AIC table
knitr::kable(AIC_table)

给予

knitr::kable(tail(AIC_table))
# |         |    dAICc| df|p.vars |
# |:--------|--------:|--:|:------|
# |model29  | 0.000000|  4|cyl    |
# |model126 | 1.021748|  5|disp   |
# |model41  | 1.245210|  4|hp     |
# |model129 | 1.340789|  5|drat   |
# |model333 | 1.582143|  6|wt     |
# |model65  | 1.995839|  5|qsec   |

knitr::kable(tail(AIC_table))
# |        |    dAICc| df|p.vars                                            |
# |:-------|--------:|--:|:-------------------------------------------------|
# |model8  | 50.75066|  3|cyl, disp, hp, drat, qsec, vs, am, gear, carb     |
# |model10 | 51.27182|  3|cyl, disp, hp, wt, qsec, vs, am, gear, carb       |
# |model53 | 53.37300|  4|cyl, disp, drat, wt, qsec, vs, am, gear, carb     |
# |model49 | 53.72970|  4|cyl, hp, drat, wt, qsec, vs, am, gear, carb       |
# |model9  | 56.79743|  3|disp, hp, drat, wt, qsec, vs, am, gear, carb      |
# |model6  | 58.92517|  3|cyl, disp, hp, drat, wt, qsec, vs, am, gear, carb |