How to generate an AIC table for GLMs that includes the names of predictor variables

177 views Asked by At

I am running a series of glms using the below code, run in R Markdown. I need the AIC table to include the variable names (not just the name of the model).

# 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 provides the model numbers (i.e., 1 to >500), but does not include the variable names. Can anyone assist with modifying the code to include these?

1

There are 1 answers

5
jay.sf On

You could add the variables as attribute to each model. Demonstrated on mtcars with slightly simplified code.

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)

(This answer elaborates on the use of do.call here.)

The use of bbmle::AICctab was recommended1. We can produce character strings out of the p.vars using toString() and cbind.

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

knitr::kable(AIC_table)
# |         |    logLik|     AICc|  dLogLik|    dAICc| df|    weight|p.vars |
# |:--------|---------:|--------:|--------:|--------:|--:|---------:|:------|
# |model29  | -67.10019| 143.6819| 30.77475| 0.000000|  4| 0.0701785|cyl    |
# |model126 | -66.19796| 144.7036| 31.67698| 1.021748|  5| 0.0421050|disp   |
# |model41  | -67.72280| 144.9271| 30.15215| 1.245210|  4| 0.0376539|hp     |
# |model129 | -66.35748| 145.0227| 31.51746| 1.340789|  5| 0.0358968|drat   |
# |model333 | -64.95201| 145.2640| 32.92294| 1.582143|  6| 0.0318160|wt     |
# |model65  | -66.68501| 145.6777| 31.18994| 1.995839|  5| 0.0258710|qsec   |
# ...
# |model8  | -93.78769| 194.4325| 4.087253| 50.75066|  3|      0|cyl, disp, hp, drat, qsec, vs, am, gear, carb     |
# |model10 | -94.04827| 194.9537| 3.826675| 51.27182|  3|      0|cyl, disp, hp, wt, qsec, vs, am, gear, carb       |
# |model53 | -93.78669| 197.0549| 4.088254| 53.37300|  4|      0|cyl, disp, drat, wt, qsec, vs, am, gear, carb     |
# |model49 | -93.96504| 197.4116| 3.909904| 53.72970|  4|      0|cyl, hp, drat, wt, qsec, vs, am, gear, carb       |
# |model9  | -96.81108| 200.4793| 1.063868| 56.79743|  3|      0|disp, hp, drat, wt, qsec, vs, am, gear, carb      |
# |model6  | -97.87495| 202.6070| 0.000000| 58.92517|  3|      0|cyl, disp, hp, drat, wt, qsec, vs, am, gear, carb |