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, 1-way ANOVA, 2+ way ANOVA, ordinary linear regression, and multiple regression.
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).
During the 2+ way ANOVA (multiple regression) section 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, the same data used in the regression and ANCOVA chapters) to look at these topics.
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
## -21.0129 -2.2501 0.1381 2.3274 19.1092
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 78.7438 6.5133 12.090 < 2e-16 ***
## HandLength 4.2840 0.3427 12.500 < 2e-16 ***
## SexMale 6.8720 0.9121 7.534 4.09e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.93 on 152 degrees of freedom
## Multiple R-squared: 0.8272, Adjusted R-squared: 0.8249
## F-statistic: 363.8 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.8272) and Adjusted R-squared (0.8249).
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 10358.5 10358.5 670.793 < 2.2e-16 ***
## Sex 1 876.6 876.6 56.768 4.087e-12 ***
## Residuals 152 2347.2 15.4
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Here, the column labelled Sum Sq is telling us what the variance EXPLAINED by each of the terms is. The final entry, for Residuals, is the amount not explained (hence “residual”). Therefore we can calculate R-squared from this as (1.03585^{4} + 876.6) / (1.03585^{4} + 876.6 + 2347.2) = 0.8271869. 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 will always increase 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 predictors in the model and the sample size11. 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)))
mod2 <- lm(Height ~ HandLength + Sex + randomVariable, data = morph)
summary(mod2)##
## Call:
## lm(formula = Height ~ HandLength + Sex + randomVariable, data = morph)
##
## Residuals:
## Min 1Q Median 3Q Max
## -21.0371 -2.2358 0.0284 2.2554 19.0136
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 79.1019 6.5570 12.064 < 2e-16 ***
## HandLength 4.2660 0.3449 12.369 < 2e-16 ***
## SexMale 6.8854 0.9144 7.530 4.29e-12 ***
## randomVariable 0.1840 0.3185 0.578 0.564
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.938 on 151 degrees of freedom
## Multiple R-squared: 0.8276, Adjusted R-squared: 0.8241
## F-statistic: 241.6 on 3 and 151 DF, p-value: < 2.2e-16We can ask for the more precise multiple R-squared values like this:
## [1] 0.8271866## [1] 0.8275676If you subtract one from the other, you will see that the Multiple R-squared has “improved” by 0.000381 — 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.0008 (from 0.8249 to 0.8241). This, hopefully, is enough evidence to make you 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 developed in 1970 by a Japanese statistician called Hirotugu Akaike. 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.
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
## mod2 5 868.6055
## mod1 4 869.0918
## mod3 3 916.2803
## mod6 4 918.1688
## mod4 3 976.6764
## mod5 2 1137.2010In the AIC results table, the models are now ordered from best (lowest AIC) to worst (highest AIC).
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 10358.5 10358.5 670.793 < 2.2e-16 ***
## Sex 1 876.6 876.6 56.768 4.087e-12 ***
## Residuals 152 2347.2 15.4
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1So here we already know that the model explains 82.72% 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 1.03585{4}/(1.03585{4} + 876.6) = 0.9219767.
Similarly, the proportion of variance explained by Sex is 876.6/(1.03585^{4} + 876.6) = 0.0780233.
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 (men tend to be both taller and larger-handed) — the share attributed to each term depends on the order in which the terms enter the model. Fitting Height ~ Sex + HandLength instead would give a somewhat different split. Variance partitioning is therefore best treated as a useful summary rather than a unique, definitive attribution.12
In more complicated multiple regression one could use this approach to further group variables into types so one could e.g., one could lump together different types of explanatory variables. 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.
20.4 Conclusion
In conclusion, you now have some tools to understand your models in more detail. R-squared gives a handy summary to tell you how much variation is explained - a high R2 value indicates a good model. It can be used to compare models that use different data. AIC is another measure of model “quality” but can only compare models that use the same data set. 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).
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 predictors in the model.↩︎
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 way of sharing the overlapping variance fairly between correlated predictors.↩︎