# Load needed libraries
library(tidyverse)
-library(lmerTest)# Note that these lines only need to be run once per computer
## So you can skip this step if you've installed these before
install.packages("tidyverse")
-install.packages("lmerTest")We’ll go ahead and load some of these libraries as well to be able to better demonstrate these concepts.
# Load needed libraries
library(tidyverse)
-library(lmerTest)We might also fit other candidate models for pairs of X, W, and Z but for the sake of simplicity in this hypothetical we’ll skip those. Note that for this method to be appropriate you need to fit the same type of model in all cases!
-Once we’ve fit all of our models and assigned them to objects, we can use the AIC function included in base R to compare the AIC score of each model. “AIC” stands for Akaike (AH-kuh-ee-kay) Information Criterion and is one of several related information criteria for summarizing a model’s explanatory power. The lowest AIC best explains the data and models with more parameters are penalized to make it mathematically possible for a model with fewer explanatory variables to still do a better job capturing the variation in the data. Technically any difference in AIC indicates model improvement but many scientists use a rule of thumb of a difference of 2. So, if two models have AIC scores that differ by less than 2, you can safely say that they have comparable explanatory power. That is definitely a semi-arbitrary threshold but so is the 0.05 threshold for p-value “significance”.
Once we’ve fit all of our models and assigned them to objects, we can use the AIC function included in base R to compare the AIC score of each model. “AIC” stands for Akaike (AH-kuh-ee-kay) Information Criterion and is one of several related information criteria for summarizing a model’s explanatory power. Models with more parameters are penalized to make it mathematically possible for a model with fewer explanatory variables to still do a better job capturing the variation in the data.
The model with the lowest AIC best explains the data. Technically any difference in AIC indicates model improvement but many scientists use a rule of thumb of a difference of 2. So, if two models have AIC scores that differ by less than 2, you can safely say that they have comparable explanatory power. That is definitely a semi-arbitrary threshold but so is the 0.05 threshold for p-value “significance”.
Let’s check out an example using AIC to compare the strengths of several models. Rather than using simulated data–as we did earlier in the mixed-effect model section–we’ll use some real penguin data included in the palmerpenguins package.
This dataset includes annual data on three penguin species spread across several islands. The sex of the penguins was also recorded in addition to the length of their flippers, body mass, and bill length and depth.
+For the purposes of this example, our research question is as follows: what factors best explain penguin body mass?
+# Load the penguins data from the `palmerpenguins` package
+data(penguins)
+
+# Make a version where no NAs are allowed
+peng_complete <- penguins[complete.cases(penguins), ]
+
+# Check the structure of it
+dplyr::glimpse(peng_complete)NA values in any column. It is better to identify and handle NAs more carefully but for this context we just want to have the same number of observations in each model
+Rows: 333
+Columns: 8
+$ species <fct> Adelie, Adelie, Adelie, Adelie, Adelie, Adelie, Adel…
+$ island <fct> Torgersen, Torgersen, Torgersen, Torgersen, Torgerse…
+$ bill_length_mm <dbl> 39.1, 39.5, 40.3, 36.7, 39.3, 38.9, 39.2, 41.1, 38.6…
+$ bill_depth_mm <dbl> 18.7, 17.4, 18.0, 19.3, 20.6, 17.8, 19.6, 17.6, 21.2…
+$ flipper_length_mm <int> 181, 186, 195, 193, 190, 181, 195, 182, 191, 198, 18…
+$ body_mass_g <int> 3750, 3800, 3250, 3450, 3650, 3625, 4675, 3200, 3800…
+$ sex <fct> male, female, female, female, male, female, male, fe…
+$ year <int> 2007, 2007, 2007, 2007, 2007, 2007, 2007, 2007, 2007…With our data in hand and research question in mind, we can fit several candidate models that our scientific intuition and the published literature support as probable then compare them with AIC.
+# Species and sex
+mod_spp <- lm(body_mass_g ~ species + sex, data = peng_complete)
+
+# Island alone
+mod_isl <- lm(body_mass_g ~ island, data = peng_complete)
+
+# Combination of species and island
+mod_eco <- lm(body_mass_g ~ island + species + sex, data = peng_complete)
+
+# Body characteristics alone
+mod_phys <- lm(body_mass_g ~ flipper_length_mm + bill_length_mm + bill_depth_mm,
+ data = peng_complete)
+
+# Global model
+mod_sink <- lm(body_mass_g ~ island + species + sex +
+ flipper_length_mm + bill_length_mm + bill_depth_mm,
+ data = peng_complete)Once we’ve fit all of these models, we can use the AIC function from base R (technically from the stats package included in base R).
# Compare models
+AIC(mod_spp, mod_isl, mod_eco, mod_phys, mod_sink) %>%
+ dplyr::arrange(AIC)AIC function doesn’t sort by AIC score automatically so we’re using the arrange function to make it easier for us to rank models by their AIC scores
+ df AIC
+mod_sink 10 4727.242
+mod_spp 5 4785.594
+mod_eco 7 4789.480
+mod_phys 5 4929.554
+mod_isl 4 5244.224Interestingly, it looks like the best model (i.e., the one that explains most of the data) is the global model that included most of the available variables. As stated earlier, it is not always the case that the model with the most parameters has the lowest AIC so we can be confident this is a “real” result. The difference between that one and the next (incidentally the model where only species and sex are included as explanatory variables) is much larger than 2 so we can be confident that the global model is much better than the next best.
+With this result your interpretation would be that penguin body mass is better explained by a combination of species, sex, physical characteristics of the individual penguin, and the penguin’s home island than it is by any of the other candidate models. In a publication you’d likely want to report this entire AIC table (either parenthetically or in a table) so that reviewers could evaluate your logic.