wampus.data <- read.csv("data/wampus.csv")7 Model selection
7.1 Model selection
Ready your best catwalk, it is time for model selection!
A typical goal of model fitting, no matter the type of GLM used, is to determine which combination of explanatory variables provides the best fit to the data. We might also be interested in figuring out in a set of ecological hypotheses which ones are the most relevant to explain our data.
Model selection is an essential tool to answer those questions, and to ensure that we are ending up with a model that leads to a good trade-off between the number of parameters used and the fit.
A full review of model selection is beyond the scope of this introduction to GLMs in R; there are literal books written solely on this topic. (And frankly, this can be a messy and controversial subject.) For now, I will just introduce you to some methods to start playing around with it in R.
We have already learned about the likelihood ratio test, permitting us to compare two models, as long as they are nested one within the other (which is already a pretty big limitation). But, it also does not give us a direct idea of how valuable adding/removing variables is (i.e., is it really useful to add a variable in the model considering the fact that it makes the model that much more complex?).
7.1.1 Using AIC to compare models
One model selection approach is to use a metric to compare models, based on fit and model complexity. One such metric is the Akaike’s Information Criterion (AIC). It provides a score for each model based on how well it fits the data (the higher the likelihood, the lower the AIC, the better), but penalizes that score by the number of parameters used to achieve that fit (more parameters = higher AIC = bad). Everything else constant, we will favor the model with the lowest AIC.
Tip
Models within 2 AIC points of each other are usually considered comparable and equally useful!
Let’s take a look on how to conduct model selection in practice by analyzing sightings of the six-legged wampus cat in the Appalachians.
The first step is to fit the models we want to compare. We can then extract the AIC values using the AIC() on the model object. We start by running a Poisson regression on Wampus cat counts in function of blood alcohol concentration of the observer:
wampus.mod.1=glm(Y.wampus~BloodAlcoholConcentration,
family="poisson",
data=wampus.data)
AIC(wampus.mod.1)[1] 110.3794Without anything else, this value is not really informative. So, let’s look at how counts relate to other variables (luminosity at time of observation, vegetation coverage, and temperature). We can run all the models, and feed all of them to our AIC() function which is going to extract all of them at once.
wampus.mod.2=glm(Y.wampus~Luminosity,
family="poisson",
data=wampus.data)
wampus.mod.3=glm(Y.wampus~VegCover,
family="poisson",
data=wampus.data)
wampus.mod.4=glm(Y.wampus~Temperature,
family="poisson",
data=wampus.data)
wampus.mod.null=glm(Y.wampus~1,
family="poisson",
data=wampus.data)
Important
In order to be able to compare models using AIC, all models need to be fitted to the SAME dataset!
AIC(wampus.mod.1,
wampus.mod.2,
wampus.mod.3,
wampus.mod.4,
wampus.mod.null) df AIC
wampus.mod.1 2 110.3794
wampus.mod.2 2 158.5802
wampus.mod.3 2 179.2260
wampus.mod.4 2 139.9127
wampus.mod.null 1 178.3202Before we go further, I should note that AIC can overfit data by selecting too many parameters, especially when sample size is small. Because of this, the AICc (AIC corrected for small samples size) was developed. The AICc applies a stronger penalty to the number of parameters for low sample sizes. As sample sizes increase, AICc converges to the ‘regular’ AIC. Because of this, it is best practice to use AICc by default. There are several packages to do that in R, but for now, we will use the MuMIn package and its AICc() function.
library("MuMIn")
AICc.tab = MuMIn::AICc(wampus.mod.1,
wampus.mod.2,
wampus.mod.3,
wampus.mod.4,
wampus.mod.null)
AICc.tab[order(AICc.tab$AICc),] # To make things easier to read, we can reorder the table by ascending AICc values df AICc
wampus.mod.1 2 111.0852
wampus.mod.4 2 140.6186
wampus.mod.2 2 159.2861
wampus.mod.null 1 178.5424
wampus.mod.3 2 179.9318Out of our 5 tested models, it appears that the first one is the best (the one relating Wampus cat sightings to the blood alcohol concentration of the observer, go figure…).
TipUpdating models
You can use the update function to update and rerun a model with an updated formula.
7.1.2 Stepwise selection
“Congratulations, you’re still in the running for becoming America’s Next Top Model”
— Tyra Banks
The approach above is great if we have a small number of pre-determined models (preferably based on clear ecological hypotheses) to compare. However, when we want have several independent variables we want to consider, we can build up our model incrementally. Two common options are to start with a full model containing all the variables, and eliminate them sequentially one by one until AIC/AICc is not improved anymore (backward stepwise model selection). Or we can start from the null model, and add variables one by one, until AIC/AICc can’t be improved further (forward stepwise model selection). A combination of both methods can also be used. Base R has the “step()” function, but I would recommend instead using the stepAIC() function in the MASS packages: it is better and more versatile.
Let’s start with our full model:
wampus.mod.full=glm(Y.wampus ~ BloodAlcoholConcentration + Luminosity + VegCover + Temperature,
family="poisson",
data=wampus.data)
summary(wampus.mod.full)
Call:
glm(formula = Y.wampus ~ BloodAlcoholConcentration + Luminosity +
VegCover + Temperature, family = "poisson", data = wampus.data)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.68046 0.79103 0.860 0.389664
BloodAlcoholConcentration 0.83696 0.18092 4.626 3.72e-06 ***
Luminosity -0.45405 0.11904 -3.814 0.000136 ***
VegCover 0.03066 0.14990 0.205 0.837919
Temperature 0.01220 0.01415 0.863 0.388381
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for poisson family taken to be 1)
Null deviance: 126.707 on 19 degrees of freedom
Residual deviance: 41.527 on 15 degrees of freedom
AIC: 101.14
Number of Fisher Scoring iterations: 5And run it through the stepAIC() function…
library(MASS)
step.mod = stepAIC(wampus.mod.full)Start: AIC=101.14
Y.wampus ~ BloodAlcoholConcentration + Luminosity + VegCover +
Temperature
Df Deviance AIC
- VegCover 1 41.569 99.182
- Temperature 1 42.281 99.895
<none> 41.527 101.140
- Luminosity 1 56.601 114.214
- BloodAlcoholConcentration 1 67.765 125.378
Step: AIC=99.18
Y.wampus ~ BloodAlcoholConcentration + Luminosity + Temperature
Df Deviance AIC
- Temperature 1 42.300 97.914
<none> 41.569 99.182
- Luminosity 1 56.726 112.339
- BloodAlcoholConcentration 1 68.789 124.403
Step: AIC=97.91
Y.wampus ~ BloodAlcoholConcentration + Luminosity
Df Deviance AIC
<none> 42.300 97.914
- Luminosity 1 56.766 110.379
- BloodAlcoholConcentration 1 104.967 158.580The output shows us the path taken to get to the selected model (this can be hidden by setting the argument trace to 0, if you’re more of a “destination over journey” kind of person), and a quick summary of the final model. The function also returns the final model in an object that can be accessed and handled as shown in previous chapters.
summary(step.mod)
Call:
glm(formula = Y.wampus ~ BloodAlcoholConcentration + Luminosity,
family = "poisson", data = wampus.data)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 1.3530 0.1306 10.356 < 2e-16 ***
BloodAlcoholConcentration 0.7010 0.0851 8.237 < 2e-16 ***
Luminosity -0.4382 0.1167 -3.754 0.000174 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for poisson family taken to be 1)
Null deviance: 126.71 on 19 degrees of freedom
Residual deviance: 42.30 on 17 degrees of freedom
AIC: 97.914
Number of Fisher Scoring iterations: 5This is an overall great and quick approach that can be really useful if there are a lot of variables involved. However, because the selection process is done one step at a time, the path taken can lead us to local optima that won’t allow for the selection of the best model. For that, we would have to test all possible combinations, which can quickly become very resource consuming. Every time we add one possible variable, we multiply the total number of possible model by 2! For ‘n’ variables, we have 2n possible models!
7.1.3 Automated model selection using dredge()
If we want to test every possible combinations of parameters, we can rely on the MuMin package and its dredge() function. This function takes a full model and generates a model selection table of all combinations of variables in the full model.
Important
The dredge() function requires that NAs in datasets trigger an error and stop the computations. This can be done globally, but I would recommend instead doing it at the level of the initial full model by specifying the following argument in the model call na.action = "na.fail".
wampus.mod.full=glm(Y.wampus ~ BloodAlcoholConcentration + Luminosity + VegCover + Temperature,
family="poisson",
data=wampus.data,
na.action = "na.fail")
model_selection_table <- dredge(wampus.mod.full)
model_selection_tableGlobal model call: glm(formula = Y.wampus ~ BloodAlcoholConcentration + Luminosity +
VegCover + Temperature, family = "poisson", data = wampus.data,
na.action = "na.fail")
---
Model selection table
(Intrc) BldAC Lmnst Tmprt VgCvr df logLik AICc delta weight
4 1.3530 0.7010 -0.4382 3 -45.957 99.4 0.00 0.642
8 0.7029 0.8297 -0.4550 0.011970 4 -45.591 101.8 2.44 0.190
12 1.3460 0.7043 -0.4369 0.02079 4 -45.947 102.6 3.15 0.133
16 0.6805 0.8370 -0.4541 0.012200 0.03066 5 -45.570 105.4 6.01 0.032
2 1.2190 0.7446 2 -53.190 111.1 11.67 0.002
10 1.2040 0.7510 0.05071 3 -53.128 113.8 14.34 0.000
6 1.0820 0.7701 0.002527 3 -53.170 113.8 14.43 0.000
14 1.0620 0.7776 0.002622 0.05116 4 -53.107 116.9 17.47 0.000
7 4.1160 -0.4421 -0.050430 3 -59.201 125.9 26.49 0.000
15 4.0730 -0.4669 -0.048720 -0.14490 4 -58.689 128.0 28.63 0.000
5 4.2330 -0.055520 2 -67.956 140.6 41.21 0.000
13 4.2230 -0.055110 -0.04091 3 -67.911 143.3 43.91 0.000
11 1.7350 -0.4822 -0.24040 3 -75.485 158.5 59.06 0.000
3 1.6700 -0.4362 2 -77.290 159.3 59.87 0.000
1 1.5370 1 -88.160 178.5 79.13 0.000
9 1.5650 -0.12390 2 -87.613 179.9 80.52 0.000
Models ranked by AICc(x) The function returns a full model selection table, already sorted by AICc for our convenience. Each row corresponds to one model. It also shows for each variable the corresponding parameter estimate, and comes with other useful information such as AICc weights. AICc weights are normalized values derived from the AICc scores that represent the relative probability of a model being the best among the set of candidates. AICc weights sum to 1 over all models in the candidate set. They give a direct and interpretable measure of support for each model. Higher weight indicates a better model.
We can easily extract information from the top models (i.e., the ones with a AICc within 2 AICc of the best model):
# View the top models
subset(model_selection_table, delta < 2)Global model call: glm(formula = Y.wampus ~ BloodAlcoholConcentration + Luminosity +
VegCover + Temperature, family = "poisson", data = wampus.data,
na.action = "na.fail")
---
Model selection table
(Intrc) BldAC Lmnst df logLik AICc delta weight
4 1.353 0.701 -0.4382 3 -45.957 99.4 0 1
Models ranked by AICc(x) We can also directly extract the top model using get.models() for further exploration and visualization!
# Extract the best model
best_model <- get.models(model_selection_table, 1)[[1]]
summary(best_model)
Call:
glm(formula = Y.wampus ~ BloodAlcoholConcentration + Luminosity +
1, family = "poisson", data = wampus.data, na.action = "na.fail")
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 1.3530 0.1306 10.356 < 2e-16 ***
BloodAlcoholConcentration 0.7010 0.0851 8.237 < 2e-16 ***
Luminosity -0.4382 0.1167 -3.754 0.000174 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for poisson family taken to be 1)
Null deviance: 126.71 on 19 degrees of freedom
Residual deviance: 42.30 on 17 degrees of freedom
AIC: 97.914
Number of Fisher Scoring iterations: 5
Tip
stepAIC() vs dredge()
The stepwise regression approach using stepAIC() is fast, making it suitable for large datasets with many variables. But, it can get stuck in local optima and ignore the “best” model if it is not along the immediate path.
Testing all possible subsets with dredge() through an exhaustive search on the other hand can be slow and computationally intensive for models with many parameters. It will identify the global optimum through all comparable models, but is more suited for small to medium-sized datasets, and where finding the “true” best model is critical.
7.2 Model validation
“All models are wrong, but some are useful.”
— George Box
CautionA note on model selection
“Best model” does not necessarily mean “good model”!
“Best model” does not necessarily mean “good model”… Sad face. I hate having to say this to you after all we’ve been through. But, it is true. And this would have been too easy anyway, right?
Because we are using an arbitrary set of models (hopefully based on ecological hypotheses), we cannot guarantee that whatever is coming out of it is actually meaningful without at least running some checks.
This is where model validation comes in: we always need to test the selected model(s) against independent data to confirm its/their predictive power. There are several model validation techniques available, such as:
Holdout Validation (Train-Test Split)
K-fold cross validation
Leave-One-Out Cross-Validation (LOOCV)
Bootstrapping
We are going to just look at one of them today, the K-fold cross validation. The concept behind it is to split the original dataset in ‘K’ subsets, train the model on all subsets but for one (our test data), and see how well they perform at predicting what is seen in the test data (which can be assessed using a variety of metrics).
The easiest option to do this in R is probably with the caret package, and its train() function. We first need to set up our training control parameters, and then run our selected model through the training function:
library(caret)
# Define the training control with 10-fold cross-validation
train_control <- trainControl(method = "cv", # 'cv' stands for cross-validation
number = 5) # K is typically set at somewhere between 5 and 10
# Train the model using the train function, applying the defined cross-validation
training.mod <- train(Y.wampus~BloodAlcoholConcentration+Luminosity ,
data = wampus.data,
method = "glm",
family="poisson",
trControl = train_control)
# Print the results, including average performance metrics like RMSE and R-squared
print(training.mod)Generalized Linear Model
20 samples
2 predictor
No pre-processing
Resampling: Cross-Validated (5 fold)
Summary of sample sizes: 16, 15, 17, 16, 16
Resampling results:
RMSE Rsquared MAE
3.076128 0.696879 2.079233Ok, so let’s decode the output, shall we?
First, we are reminded of the kind of model we fitted, the sample size, and the number of predictors. So far, nothing we haven’t seen before.
Second, we get some processing information:
We did not pre-process the data (e.g., no scaling, no centering, no filtering, etc.). You can check the
preProcOptionsargument of thetrainControlfunction, and the associated help file? caret::preProcessfor more information.We are reminded of the validation resampling method used (here a 5-fold cross-validation).
We are kindly informed of the sample sizes of each ‘fold’.
Finally, and this is the part we’ve all been waiting for, the resampling results. train provides us with the following metrics: RMSE, R-squared, and MAE. These metrics give us an idea of how well the model performed on previously unseen data (i.e., the training data). How about we take a closer look at each of them?
RMSE: This is the Root Mean Squared Error. This measure combines bias (difference between expected value of the estimator and the parameter) and precision in one overall measure of how “close” the predicted values are to the observed data on average. The lower the RMSE, the more closely a model can predict the actual observations.
R-squared: A familiar sight! This is a measure of the correlation between the model predictions made and the data. The higher the R2, the better the model is at predicting observations.
MAE: My favorite. The Mean Absolute Error. (First, it’s a pretty cool name, you gotta admit!) It measures the average absolute difference between the predictions and the observations. Both RMSE and MAE are expressed in the same unit as the response variable, but MAE is less sensitive to outliers, and frankly more directly interpretable. It simply tells you how far, on average, the predicted values are from the true values. Unsurprisingly, the lower the MAE, the more accurately a model can predict the actual observations.
With that, we can determine if we are content with how well the selected model performs with predictions. If not, time to get back to the drawing board. But in this case, I am pretty happy with the results!
The caret package can also handle several other validation methods. You can select them with the method argument in trainControl (e.g., “boot”, “cv”, “repeatedcv”, “LOOCV”, …). Feel free to explore!
7.3 Words of wisdom: of the importance of putting hypotheses first.
The tools described above for model selection and validation are powerful, but are only that: tools. They need to be used mindfully. As scientists, our job is to put our hypotheses first, lest we want to get lost in the treacherous waters of technical stats. Statistics are here to help use make sense of the data, but this is only possible if we put the data first. We need to think carefully of the processes that led to these data, and how the statistical approach we use to analyze them is going to get us to our answer.
While model dredging can seem tempting and can be useful, it should be used cautiously to avoid selecting models with no ecological meaning.
Interpretability vs. Fit: we have to balance the needs for flexible and well-fitting models against how interpretable they are. What is the point of a good model if it actually does not shed light on the biological processes we are attempting to understand?
Put your hypotheses first!