Chapter 20 Evaluating linear models
We have now covered a range of linear models that are all fitted using the same tool in R (lm): t-tests, one-way ANOVA, linear regression, ANCOVA and n-way ANOVA (the last two being models with more than one explanatory variable).
The models are all fitted in the same way and have the same assumptions. We have already covered the four key diagnostic plots (see the ANOVA and linear regression sections, and the GSWR textbook), how to evaluate the significance of parameters, and the meaning of the coefficients.
There are some additional useful points to consider: proportion of variance explained (R-squared value), proportion of variance explained by different variables in the model, and Akaike’s Information Criterion (and likelihood).
In the ANCOVA and n-way ANOVA chapters you may have realised that there may be multiple ways to fit a model. For example, you may have a choice of parameters to include - should you include them or not? Which ones should you include? Would a log-transformed explanatory variable be better?
We will use the morphometry.csv data (body measurements from a sample of people, the same data used in the regression and ANCOVA chapters) to look at these topics.
In this chapter we use the dplyr package:
morph <- read.csv("CourseData/morphometry.csv") %>%
mutate(Height = Height / 10, HandLength = HandLength / 10) # mm -> cm20.1 R-squared value
Let’s fit a simple model: Height ~ HandLength + Sex
##
## Call:
## lm(formula = Height ~ HandLength + Sex, data = morph)
##
## Residuals:
## Min 1Q Median 3Q Max
## -10.7529 -2.6604 -0.3563 2.9424 9.7437
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 85.6144 7.5536 11.33 < 2e-16 ***
## HandLength 3.9115 0.3967 9.86 < 2e-16 ***
## SexMale 7.4154 1.0159 7.30 1.5e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.064 on 152 degrees of freedom
## Multiple R-squared: 0.8057, Adjusted R-squared: 0.8032
## F-statistic: 315.2 on 2 and 152 DF, p-value: < 2.2e-16The model summary here shows us the R-squared value, which is a measure of the proportion of variation in the response variable that is explained by variation in the explanatory variable(s). In a linear regression, an r-squared value of 1 (100%) would mean that all data points fall on the line. As the r-squared value declines, there exists more noise in the relationship (i.e. the points become more spread out around the line).
There are two types of R-squared value shown in this summary: Multiple R-squared (0.8057) and Adjusted R-squared (0.8032).
We’ll look at multiple R-squared first. This value is calculated as the amount of explained variation divided by the total amount of variation. Take a look at the anova summary table:
## Analysis of Variance Table
##
## Response: Height
## Df Sum Sq Mean Sq F value Pr(>F)
## HandLength 1 9535.3 9535.3 577.198 < 2.2e-16 ***
## Sex 1 880.3 880.3 53.285 1.501e-11 ***
## Residuals 152 2511.0 16.5
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Here, the column labelled Sum Sq tells us how much of the variation each term EXPLAINS, measured as a sum of squares. The final entry, for Residuals, is the amount not explained (hence “residual”). Therefore we can calculate R-squared from this as (9535.3 + 880.3) / (9535.3 + 880.3 + 2511) = 0.8057. You can check that this matches the figure indicated by summary(mod1).
But what is the Adjusted R-squared?
Multiple R-squared is a measure of R-squared value for models that can have multiple predictor variables. It measures the amount of variation in the response variable that can be explained by the predictor variables. When additional terms are added to the model, the multiple R-squared can never decrease (and in practice almost always increases) because terms will always explain some portion of the variance, even if it is very small. This behaviour can be a bit annoying, so adjusted R-squared controls against this increase by adding a penalty based on the number of coefficients in the model and the sample size14. When reporting R-squared values for models with >1 term you should report the adjusted R-squared value.
You can test this by adding terms to the model. Let’s start with something silly - we’ll add a variable that is simply a vector of random numbers to the model. By definition this cannot have any meaningful explanatory power, but what will it do to the multiple R-squared value?
morph <- morph %>%
mutate(randomVariable = rnorm(nrow(morph)))
mod1_random <- lm(Height ~ HandLength + Sex + randomVariable, data = morph)
summary(mod1_random)##
## Call:
## lm(formula = Height ~ HandLength + Sex + randomVariable, data = morph)
##
## Residuals:
## Min 1Q Median 3Q Max
## -10.455 -2.656 -0.355 3.010 9.501
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 86.0710 7.5998 11.325 < 2e-16 ***
## HandLength 3.8864 0.3993 9.733 < 2e-16 ***
## SexMale 7.4893 1.0240 7.314 1.42e-11 ***
## randomVariable -0.2159 0.3295 -0.655 0.513
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.072 on 151 degrees of freedom
## Multiple R-squared: 0.8063, Adjusted R-squared: 0.8024
## F-statistic: 209.5 on 3 and 151 DF, p-value: < 2.2e-16We can ask for the more precise multiple R-squared values like this:
## [1] 0.8057465## [1] 0.8062974If you subtract one from the other, you will see that the Multiple R-squared has “improved” by 0.0005509 : an increase that happens automatically whenever any term is added to a model, even a completely useless one. The Adjusted R-squared, however, is not so easily fooled: because it penalises each extra parameter, the random variable earns essentially no reward. Here the Adjusted R-squared actually decreases slightly, changing by just -0.0007 (from 0.8032 to 0.8024). Be aware, though, that adjusted R-squared is not a noise detector: it rises whenever a new term improves the fit by more than an average noise term would, so with other random numbers (try a different seed) it rises roughly one time in three. Even so, it is a better summary than multiple R-squared, so you should favour reporting adjusted rather than multiple R-squared values.
The adjusted R-squared value can be used as a single-number summary of model explanatory power.
20.2 Akaike Information Criterion (AIC)
The Akaike information criterion (AIC) is an estimator of prediction error in a statistical model introduced in the early 1970s by the Japanese statistician Hirotugu Akaike (presented in 1971, published in 1973 and 1974). In a nutshell, it is a measure of the relative quality of statistical models for a given set of data. This last part is important. AIC is only comparable among statistical models that use the same data (and which have the same response variable). In other words, given a collection of plausible models that use the same data set, AIC estimates the quality of each model relative to the others. If you are interested in the details you can read further on Wikipedia or a more advanced statistics book; otherwise you can simply trust that AIC estimates the relative quality of models, with lower values being better.
The formula is \(AIC = -2 \ln(L) + 2k\), where \(L\) is the likelihood (how probable the observed data are under the fitted model) and \(k\) is the number of estimated parameters. The first part rewards a good fit; the second part penalises complexity, adding 2 for each extra parameter. For lm models, \(k\) includes the residual standard deviation as well as the coefficients, which is why the df column in the table below is one larger than the number of coefficients.
You can get R to tell you the AIC value for a model using the function AIC() e.g. AIC(mod1).
Here’s a simple example of use in practice:
mod1 <- lm(Height ~ HandLength + Sex, data = morph)
mod2 <- lm(Height ~ HandLength * Sex, data = morph)
mod3 <- lm(Height ~ HandLength, data = morph)
mod4 <- lm(Height ~ Sex, data = morph)
mod5 <- lm(Height ~ 1, data = morph)
mod6 <- lm(Height ~ HandLength + randomVariable, data = morph)
(AICtable <- AIC(mod1, mod2, mod3, mod4, mod5, mod6) %>%
arrange(AIC))## df AIC
## mod1 4 879.5499
## mod2 5 881.3712
## mod3 3 924.1299
## mod6 4 926.1124
## mod4 3 954.1882
## mod5 2 1129.5315In the AIC results table, the models are now ordered from best (lowest AIC) to worst (highest AIC).
AIC values only mean something relative to each other, so it is usual to look at the difference between each model and the best one (\(\Delta\)AIC):
## df AIC deltaAIC
## mod1 4 879.5499 0.000000
## mod2 5 881.3712 1.821322
## mod3 3 924.1299 44.580005
## mod6 4 926.1124 46.562507
## mod4 3 954.1882 74.638320
## mod5 2 1129.5315 249.981632Common rules of thumb (Burnham and Anderson 2002) are that models within about 2 AIC units of the best have similar support, models 4 to 7 units worse have considerably less support, and models more than 10 units worse have essentially none. When a more complex model is within 2 units of a simpler one, as the interaction model mod2 is here (1.8 units worse than mod1), the extra term adds essentially nothing and the simpler model is preferred (Arnold 2010). Finally, choose a small set of sensible candidate models before comparing them: fitting every possible combination of variables and picking the lowest AIC will find patterns that are there by chance.
20.3 Variance partitioning
When you have a model with numerous terms (e.g. a multiple regression model, or an 2-way ANOVA for example) it is often useful to express the results in term of variance explained.
We can do this using variance partitioning.
Consider our earlier model lm(Height ~ HandLength + Sex). What proportion of the variance in height is explained by hand length? And what proportion by Sex? (and so on, for more complicated models…)
This is done by examining the anova summary (e.g. anova(mod1)), using the Sum Sq column.
## Analysis of Variance Table
##
## Response: Height
## Df Sum Sq Mean Sq F value Pr(>F)
## HandLength 1 9535.3 9535.3 577.198 < 2.2e-16 ***
## Sex 1 880.3 880.3 53.285 1.501e-11 ***
## Residuals 152 2511.0 16.5
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1So here we already know that the model explains 80.6% of variation in Height. We can partition this among the terms by using the Sums of Squares values. The proportion of variance explained by HandLength is 9535.3/(9535.3 + 880.3) = 0.92.
Similarly, the proportion of variance explained by Sex is 880.3/(9535.3 + 880.3) = 0.08.
A word of caution. These proportions come from the Sum Sq column of anova(), which uses sequential (Type I) sums of squares: each term is credited with the variation it explains after the terms listed before it in the model formula. When the predictors are correlated (as HandLength and Sex are here, because men tend to have larger hands than women), the share attributed to each term depends on the order in which the terms enter the model. Fitting Height ~ Sex + HandLength instead gives a very different split: the share credited to HandLength falls from about 92% to about 15%. This happens because most of the explained variation can be explained by either predictor; that shared part goes to whichever term enters first. Variance partitioning is therefore best treated as a useful summary rather than a unique, definitive attribution.15
In more complicated models, this approach can be taken further by grouping explanatory variables into types. Imagine you had data on human cholesterol level, for example. You might have explanatory variables including various genotypes, morphology (height/weight), various dietary factors and so on. After partitioning variance among the many variables, it could then be useful to group these variables into a smaller number of “types”, such as “genetic”, “morphological” and “diet”. Thus, variance partitioning can help make sense of complex data and can improve how such results are communicated.
The logical process is the same for Generalised Linear Models (GLM), which we will cover soon, except we use an analogous quantity called Deviance rather than Sum of Squares.
Try it: Overfitting and AIC Lab (open the lab)
Fit polynomials of increasing degree to a sample of daily temperatures.
- With 20 sampled days, which polynomial degree predicts the unseen days best? Does AIC pick the same degree?
- Draw a new sample of days, then try 80 days. Does the best degree change? Why do small samples overfit sooner?
20.4 Conclusion
In conclusion, you now have some tools to understand your models in more detail. R-squared gives a handy summary of how much of the variation in the response is explained (report adjusted R-squared for models with more than one term). It is not a measure of effect size, and a low value does not mean a bad model: noisy biological systems often give low R-squared values even when an effect is clear and precisely estimated. Because R-squared depends on the range and noise of each data set, do not use it to compare models fitted to different data. AIC is another measure of model “quality” that can only compare models fitted to the same data set, with the same response variable. Low AIC values are better than high ones. Variance partitioning can be used as a handy way to sum up your model (in addition to significance and coefficient values).
References
The formula is: \(R^2_{adj} = 1 - \left( \frac{(1 - R^2)(n - 1)}{n - k - 1} \right)\), where n is the sample size and k is the number of estimated coefficients excluding the intercept. Note that k counts coefficients, not variables: an interaction adds coefficients, and a categorical variable with L levels adds L - 1.↩︎
If you need an order-independent partition where the contributions still add up to the model’s R-squared, the standard tool is the
lmgmetric in therelaimpopackage (relaimpo::calc.relimp(mod1, type = "lmg")). It averages each predictor’s contribution over every possible order of entry: for two predictors this is simply the average of the two orderings (HandLength + SexandSex + HandLength). This is equivalent to the “Shapley value” from game theory: a principled convention for dividing the overlapping variance between correlated predictors (a convention, because the shared part cannot truly be attributed to either).↩︎