Extracting and plotting AIC weights

AIC weights

# Nicholas Judd | Donders Institute
# 26-12-2022
# njudd.com

# replicating https://twitter.com/rogierK/status/1004443199031607296/photo/2
# Wagenmakers, E.-J., & Farrell, S. (2004). 
# AIC model selection using Akaike weights. 
# Psychonomic Bulletin & Review, 11(1), 192–196.

if (!require(pacman)){
  install.packages('pacman')
}
set.seed(42) # a seed where the noise is uninformative; also a magical number
pacman::p_load(qpcR, AICcmodavg)

df <- iris # using the iris dataset
df$noise <- rnorm(dim(iris)[1])

m1 <- lm(Sepal.Width ~ Petal.Length, data = df)
m2 <- lm(Sepal.Width ~ Petal.Length + Petal.Width, data = df)
m3 <- lm(Sepal.Width ~ Petal.Length + Petal.Width + Species, data = df)
m4 <- lm(Sepal.Width ~ Petal.Length + Petal.Width + Species + noise, data = df)

mod_list <- list(m1, m2, m3, m4)
mod_list.names <- c('m1', 'm2', 'm3', 'm4')
# manually with qpcR::akaike.weights()

AICvec <- map_dbl(mod_list, AIC) # AICvec <- c(AIC(m1), AIC(m2), AIC(m3), AIC(m4))

AICvec_weights <- akaike.weights(AICvec)$weights

# replicating Rogier's graph
plot_df <- data.frame(model = mod_list.names,
                      variable = rep("AIC weight", length(mod_list.names)),
                      values_qpcR = AICvec_weights)
ggplot(plot_df, aes(model, values_qpcR, fill = model)) +
  geom_bar(stat = "identity", position = "dodge") +
  ylim(0, 1) +
  xlab("Akaike weights") +
  theme_classic(base_size = 15)

# now with aictable
aictab(mod_list, modnames = mod_list.names)

# adding the aic weights from aictab to the df
plot_df$values_AICcmodavg <- aictab(mod_list, sort =F)$AICcWt

ggplot(plot_df, aes(model, values_AICcmodavg, fill = model)) +
  geom_bar(stat = "identity", position = "dodge") +
  ylim(0, 1) +
  xlab("Akaike weights") +
  theme_classic(base_size = 15)


# the values differ slightly... 
# because aictab uses the second-order Akaike information criterion
# which is good when sample size is small
# (Burnham & Anderson 2002 recommend its use when n / K < 40)
aictab(mod_list, modnames = mod_list.names, second.ord = FALSE)
# now the AIC values and weights are identical to the akaike.weights function