Lecture 26 Comparison of general linear models

In class version

26.1 Learning objectives

By the end of this lecture you should be able to:

  1. Identify when one linear model is nested within another.
  2. State the hypothesis tested by comparing two nested models.
  3. Use a partial F test to decide whether the additional terms in a larger model improve the fit.
  4. Interpret the change in RSS, degrees of freedom, F statistic and P-value in a nested-model comparison.
  5. Explain the difference between anova(model) and anova(reduced_model, full_model).
  6. Apply nested-model comparisons to the Samara and body-fat examples introduced in the course.

26.2 Quick recall of Lecture 25

Lecture 25 gave us three increasingly flexible ways to combine a factor A with a numerical covariate x. The factor can do nothing, shift the intercept, or also change the slope.

The progression was:

Y ~ x
Y ~ A + x
Y ~ A * x

Now we ask the next question: is that extra flexibility actually supported by the data?

26.3 Nested linear models

A model is nested within another model if it can be obtained by placing restrictions on the parameters of the larger model.

For one factor A and one numerical covariate x, the nesting relationships are:

Nesting relationships for one factor and one numerical covariate.
Nesting relationships for one factor and one numerical covariate.

Both Y ~ A and Y ~ x are nested within Y ~ A + x. However, Y ~ A and Y ~ x are not generally nested within each other.

The separate-lines model Y ~ A * x contains the parallel-lines model Y ~ A + x, because setting all A:x interaction coefficients to zero produces the parallel-lines model.

26.4 Recall from Lecture 14: the partial F test

We do not need a new test here. Lecture 14’s notes already gave us nested-model comparisons and the partial F test. We now use the same machinery for the factor-and-covariate models from Lecture 25, where one scientific effect can correspond to several coefficients.

Adding terms to a linear model cannot increase its residual sum of squares. A larger model has more flexibility, so its RSS will be the same or smaller. The question is whether the decrease in RSS is large enough to justify the additional parameters, relative to the residual variation that remains in the larger model.

Let \(M_R\) be a reduced model and \(M_F\) be a full model, with \(M_R\) nested within \(M_F\). The hypotheses are

\[ H_0:\text{the additional coefficients in }M_F\text{ are all zero} \]

versus

\[ H_1:\text{at least one of the additional coefficients is non-zero}. \]

The partial F statistic is

\[ F = \frac{ (RSS_R-RSS_F)/(df_R-df_F) }{ RSS_F/df_F }. \]

The numerator is the reduction in RSS per additional parameter.

The denominator is the residual mean square from the full model, which estimates the remaining error variance.

A large F statistic means that the additional terms have removed a large amount of residual variation relative to the background residual variation.

Under \(H_0\), the statistic follows an F distribution with \(df_R-df_F\) and \(df_F\) degrees of freedom.

The models must be nested and fitted to the same observations. The usual linear-model assumptions also apply.

26.5 The Samara data revisited

Recall that we model the mean speed of fall of samara as a function of disk loading and the tree from which each fruit fell.

Grab the data here: samara.csv

We use the same three models as Lecture 25:

m1 <- lm(Velocity ~ Load, data = Samara)
m2 <- lm(Velocity ~ Load + TreeF, data = Samara)
m3 <- lm(Velocity ~ Load * TreeF, data = Samara)

Here m1 is one common regression line, m2 gives parallel lines with different intercepts but a common slope, and m3 gives separate lines with different intercepts and potentially different slopes.

26.6 Comparing the Samara models

First look at what happens as we make the model more flexible: residual degrees of freedom decrease and RSS gets smaller.

Table 26.1: Table 26.2: Residual information for the three Samara models
Model Formula Residual df RSS
m1: one line Velocity ~ Load 33 0.215
m2: parallel lines Velocity ~ Load + TreeF 31 0.203
m3: separate lines Velocity ~ Load * TreeF 29 0.165

Before looking at the tests, compare the fitted lines. The observations, colours, axes, and scales are the same in all three panels.

unlabelled

Figure 26.1: The same Samara observations with fitted lines from the three candidate models.

Moving from m1 to m2 allows the trees to move vertically apart. Moving from m2 to m3 also allows their slopes to differ.

Model Formula What changed?
m1 Velocity ~ Load One common line
m2 Velocity ~ Load + TreeF Adds TreeF: 2 additional df, allowing vertical shifts
m3 Velocity ~ Load * TreeF Adds TreeF:Load: 2 additional df, allowing slope differences

Now we can ask the two questions we actually care about:

  • Does TreeF shift the lines vertically? Compare m1 with m2.
  • Do the slopes differ by Tree? Compare m2 with m3.

26.7 Comparison 1: Do the slopes need to differ?

m3 adds the TreeF:Load interaction to m2. The null hypothesis is therefore that the interaction coefficients are all zero.

So, once we allow the trees to have different intercepts, is there evidence that the relationship between Load and Velocity also differs among trees?

anova(m2, m3)
Analysis of Variance Table

Model 1: Velocity ~ Load + TreeF
Model 2: Velocity ~ Load * TreeF
  Res.Df     RSS Df Sum of Sq     F  Pr(>F)  
1     31 0.20344                             
2     29 0.16549  2  0.037949 3.325 0.05011 .
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

The evidence for different slopes is borderline. At a strict 5% significance level we would not reject the parallel-slopes model, but a P-value this close to 0.05 should not be interpreted as evidence that the slopes are identical.

26.8 Comparison 2: If slopes are parallel, do the trees need different intercepts?

m2 adds TreeF to the model that already contains Load. This tests whether Tree provides additional information about Velocity after accounting for Load.

anova(m1, m2)
Analysis of Variance Table

Model 1: Velocity ~ Load
Model 2: Velocity ~ Load + TreeF
  Res.Df     RSS Df Sum of Sq      F Pr(>F)
1     33 0.21476                           
2     31 0.20344  2  0.011322 0.8626 0.4319

There is little evidence that Tree requires an additional intercept shift after accounting for Load.

26.9 Comparison 3: Joint comparison

This comparison adds both the Tree main effect and the Tree-by-Load interaction at once. It tests whether allowing Tree to affect either the intercept or the slope improves the model compared with one common regression line.

anova(m1, m3)
Analysis of Variance Table

Model 1: Velocity ~ Load
Model 2: Velocity ~ Load * TreeF
  Res.Df     RSS Df Sum of Sq      F  Pr(>F)  
1     33 0.21476                              
2     29 0.16549  4  0.049272 2.1585 0.09885 .
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

This is a different question from carrying out the previous two tests separately. The joint comparison has a P-value of approximately 0.099.

26.10 Why explicit comparisons matter

This is where the orthogonality point from Lecture 23 matters. In an orthogonal design, changing the order of terms does not change the sums of squares attributed to them. The Samara data are not orthogonal: the trees were sampled over different Load distributions, so Tree and Load overlap in the information they carry.

unlabelled

Figure 26.2: The Samara data have different Load distributions across trees.

That means the variation attributed to Tree depends on whether Load has already been accounted for. The two uses of anova() below therefore answer different questions.

samara_tree_first <- lm(Velocity ~ TreeF + Load, data = Samara)
samara_load_first <- lm(Velocity ~ Load + TreeF, data = Samara)

anova(samara_tree_first)
Analysis of Variance Table

Response: Velocity
          Df  Sum Sq  Mean Sq F value    Pr(>F)    
TreeF      2 0.53942 0.269708  41.098 1.913e-09 ***
Load       1 0.31554 0.315542  48.082 8.884e-08 ***
Residuals 31 0.20344 0.006563                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(samara_load_first)
Analysis of Variance Table

Response: Velocity
          Df  Sum Sq Mean Sq  F value    Pr(>F)    
Load       1 0.84364 0.84364 128.5517 1.471e-12 ***
TreeF      2 0.01132 0.00566   0.8626    0.4319    
Residuals 31 0.20344 0.00656                       
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

anova(model) produces a sequential, or Type I, ANOVA table. Terms are added to the model in sequence, so in a non-orthogonal design the sums of squares can depend on the order of the terms.

In contrast, anova(reduced_model, full_model) directly compares two fitted nested models. The hypothesis is determined by the terms that are present in the full model but absent from the reduced model.

If our scientific question is specifically whether Tree adds information after accounting for Load, the comparison should be written explicitly:

anova(m1, m2)
Analysis of Variance Table

Model 1: Velocity ~ Load
Model 2: Velocity ~ Load + TreeF
  Res.Df     RSS Df Sum of Sq      F Pr(>F)
1     33 0.21476                           
2     31 0.20344  2  0.011322 0.8626 0.4319

Here m1 already contains Load, and m2 adds TreeF. The test therefore asks whether Tree improves the model after Load has already been accounted for.

Type II and Type III sums of squares provide other conventions for adjusted tests in non-orthogonal designs. We do not need them here because we can state the scientific comparison directly using reduced and full models.

26.11 A small R detail

When several nested models are supplied in one anova() call, R uses the largest model’s residual mean square as the common denominator. For a specific scientific hypothesis, use anova(reduced, full) so the comparison is explicit. This is why the m1 versus m2 result can differ when m3 is included in the same call.

anova(m1, m2, m3)
Analysis of Variance Table

Model 1: Velocity ~ Load
Model 2: Velocity ~ Load + TreeF
Model 3: Velocity ~ Load * TreeF
  Res.Df     RSS Df Sum of Sq     F  Pr(>F)  
1     33 0.21476                             
2     31 0.20344  2  0.011322 0.992 0.38306  
3     29 0.16549  2  0.037949 3.325 0.05011 .
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Model comparison does not replace model checking. Once we have chosen the mean structure, the usual LINE diagnostics still need to be checked.

26.12 Body Fat case study

Now do the same thing with the Body Fat data. The response is percentage body fat, Age is the numerical covariate, and Gender is the factor.

Grab the data here: fat.csv

Use the same nested-model logic as for the Samara data. Does Gender need to change the Age slope, the intercept, or neither?

F0 <- lm(Percent.Fat ~ Age, data = Fat)
F1 <- lm(Percent.Fat ~ Gender + Age, data = Fat)
F2 <- lm(Percent.Fat ~ Gender * Age, data = Fat)

Before looking at the tests, predict which comparison tests the slope difference and which tests the intercept difference.

unlabelled

Figure 26.3: The same Body Fat observations with fitted lines from the three candidate models.

26.13 Does Gender need to change the slope?

F2 adds one Gender-by-Age interaction coefficient. This tests whether the Age slope differs between females and males.

anova(F1, F2)
Analysis of Variance Table

Model 1: Percent.Fat ~ Gender + Age
Model 2: Percent.Fat ~ Gender * Age
  Res.Df    RSS Df Sum of Sq      F Pr(>F)  
1     15 360.88                             
2     14 282.02  1    78.853 3.9144 0.0679 .
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

The P-value is approximately 0.068, so there is not strong evidence at the 5% level that separate slopes are required. We therefore continue with the parallel-lines model F1.

26.14 Does Gender need to change the intercept?

F1 adds the Gender intercept shift to a model that already contains Age.

anova(F0, F1)
Analysis of Variance Table

Model 1: Percent.Fat ~ Age
Model 2: Percent.Fat ~ Gender + Age
  Res.Df    RSS Df Sum of Sq      F  Pr(>F)  
1     16 529.66                              
2     15 360.88  1    168.79 7.0157 0.01824 *
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

The P-value is approximately 0.018, providing evidence that Gender contributes to expected percentage body fat after accounting for Age.

Of these three prespecified models, F1 is the simplest model supported by these comparisons: the fitted lines have different intercepts but a common Age slope.

When the larger model adds exactly one coefficient, the partial F test and the corresponding two-sided t test are equivalent:

\[F=t^2.\]

This is why the P-value from anova(F1, F2) matches the P-value for the Gender:Age coefficient in summary(F2), and why anova(F0, F1) matches the Gender coefficient test in summary(F1).

26.15 Fitted equations for the retained Body Fat model

summary(F1)

Call:
lm(formula = Percent.Fat ~ Gender + Age, data = Fat)

Residuals:
   Min     1Q Median     3Q    Max 
-6.638 -3.455 -1.103  3.297  8.952 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)  
(Intercept)  15.0708     6.2243   2.421   0.0286 *
GenderM      -9.7914     3.6966  -2.649   0.0182 *
Age           0.3392     0.1196   2.835   0.0125 *
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 4.905 on 15 degrees of freedom
Multiple R-squared:  0.7461,    Adjusted R-squared:  0.7123 
F-statistic: 22.04 on 2 and 15 DF,  p-value: 3.424e-05

With females as the reference level:

\[ E(\text{Fat}\mid\text{Female}) = \beta_0+\beta_{\text{Age}}\text{Age} \]

and

\[ E(\text{Fat}\mid\text{Male}) = (\beta_0+\beta_{\text{GenderM}}) + \beta_{\text{Age}}\text{Age}. \]

26.16 Body Fat fitted equations

The common Age coefficient gives the slope of both lines. GenderM gives the vertical difference between the male and female fitted lines.

co <- coef(F1)
beta0 <- unname(co["(Intercept)"])
beta_age <- unname(co["Age"])
beta_male <- unname(co["GenderM"])
tibble(model = c("Females", "Males"), equation = c(sprintf("Fat = %.2f + %.2f * Age",
    beta0, beta_age), sprintf("Fat = %.2f + %.2f * Age", beta0 + beta_male, beta_age))) |>
    kable(caption = "Fitted equations by Gender (parallel-slopes model)") |>
    kable_styling(full_width = FALSE)
Table 26.3: Table 26.4: Fitted equations by Gender (parallel-slopes model)
model equation
Females Fat = 15.07 + 0.34 * Age
Males Fat = 5.28 + 0.34 * Age

26.17 Reporting a nested-model comparison

A good write-up of a nested-model comparison needs five things:

  1. State the reduced and full models.
  2. State what additional terms are being tested and translate that into the scientific question.
  3. Report the change in RSS and degrees of freedom, the F statistic and the P-value.
  4. State whether the data provide evidence that the additional model structure is needed.
  5. Interpret the retained model in terms of the response and explanatory variables.

A non-significant comparison does not prove that the smaller model is true. It means that the data do not provide sufficient evidence that the additional structure improves the model.

Everything in this lecture was manageable because the candidate models were specified in advance. The scientific question told us which reduced model to compare with which full model.

The next problem is harder. What if there are many plausible predictors and no single obvious reduced-versus-full comparison? That takes us from model comparison to model selection, which is where we go next.