In tutorial 10, we extended our model of penguin flipper length by including bill length as an additional predictor and learnt how to centre, scale, and standardise variables in the linear model to change the interpretation of the individual \(b\) coefficients. The last tutorial of this module will guide you through model comparison using the F-test with R’s anova() function. Before you start this tutorial, make sure you are comfortable with the concepts discussed in lecture 10.

Warm-up quiz

As usual, here’s a little quiz to get you in the mood. Consider this linear model when answering the questions below.

Call:
lm(formula = live ~ love + laugh, data = yolo)

Residuals:
    Min      1Q  Median      3Q     Max 
-49.382 -12.559  -0.789  12.820  50.012 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   48.689      1.578  30.850  < 2e-16 ***
love          10.202      1.212   8.419 1.76e-15 ***
laugh         -3.306      2.373  -1.393    0.165    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 19.14 on 291 degrees of freedom
Multiple R-squared:  0.201, Adjusted R-squared:  0.1955 
F-statistic:  36.6 on 2 and 291 DF,  p-value: 6.609e-15

Setting Up

As usual, all you need is this tutorial and RStudio. Remember that you can easily switch between windows with the Alt + ↹ Tab (Windows) and ⌘ Command + ⇥ Tab (Mac OS) shortcuts.

Task 1

Create a new week_10 project and open it in RStudio. Then, create two new folders in the new week_10 folder: r_docs and data. Finally, open a new R Markdown file and save it in the r_docs folder. Get into the habit of creating new code chunk in the file for each of the tasks below.

Task 2

In this tutorial, just like in the previous two, we’ll continue using the penguins data set from the palmerpenguins package. You will also need the packages tidyverse, broom, and QuantPsyc (install it if you don’t have it yet; watch out for spelling and upper/lower case letters!). Put the code that loads all of these in the setup chunk of your Rmd file and run it.

library(palmerpenguins)
library(tidyverse)
library(broom)
library(QuantPsyc)

Data

Task 3

Just like the last couple of times, save the penguins data set in your environment as peng_data to avoid having to type out palmerpenguins::penguins every time you want to use it.

peng_data <- palmerpenguins::penguins

Task 4

It’s always good to refresh your memory about the dataset by looking at a rough summary.

summary(peng_data)

      species          island    bill_length_mm  bill_depth_mm  
 Adelie   :152   Biscoe   :168   Min.   :32.10   Min.   :13.10  
 Chinstrap: 68   Dream    :124   1st Qu.:39.23   1st Qu.:15.60  
 Gentoo   :124   Torgersen: 52   Median :44.45   Median :17.30  
                                 Mean   :43.92   Mean   :17.15  
                                 3rd Qu.:48.50   3rd Qu.:18.70  
                                 Max.   :59.60   Max.   :21.50  
                                 NA's   :2       NA's   :2      
 flipper_length_mm  body_mass_g       sex           year     
 Min.   :172.0     Min.   :2700   female:165   Min.   :2007  
 1st Qu.:190.0     1st Qu.:3550   male  :168   1st Qu.:2007  
 Median :197.0     Median :4050   NA's  : 11   Median :2008  
 Mean   :200.9     Mean   :4202                Mean   :2008  
 3rd Qu.:213.0     3rd Qu.:4750                3rd Qu.:2009  
 Max.   :231.0     Max.   :6300                Max.   :2009  
 NA's   :2         NA's   :2                                 

Adding a categorical predictor

Let’s start by rebuilding one of the models we fitted last time.

Task 5

Build a linear model predicting flipper length by body mass in kilogrammes and bill length in millimetres, with both predictors mean-centred.

peng_data <- peng_data %>%
  mutate(
    mass_kg = body_mass_g / 1000,
    mass_kg_cntr = mass_kg - mean(mass_kg, na.rm = T), # using mean
    bill_length_cntr = scale(bill_length_mm, scale = FALSE)) # using scale

m_massKg_bill_cntr <- lm(flipper_length_mm ~ mass_kg_cntr + bill_length_cntr, peng_data)

Task 6

Create a new model m_mass_bill_sex that adds the sex variable in addition to the two predictors in the previous model.

# using update()
m_mass_bill_sex <- update(m_massKg_bill_cntr, ~ . + sex)
# or build it from scratch
m_mass_bill_sex <- lm(flipper_length_mm ~ mass_kg_cntr + bill_length_cntr + sex, peng_data)

Evaluating results

Task 7

Look at the summary of this model and answer the following questions.

summary(m_mass_bill_sex)

Call:
lm(formula = flipper_length_mm ~ mass_kg_cntr + bill_length_cntr + 
    sex, data = peng_data)

Residuals:
    Min      1Q  Median      3Q     Max 
-23.725  -3.917   0.330   4.199  14.205 

Coefficients:
                  Estimate Std. Error t value Pr(>|t|)    
(Intercept)      203.22727    0.49903 407.248  < 2e-16 ***
mass_kg_cntr      14.01511    0.53535  26.179  < 2e-16 ***
bill_length_cntr   0.60591    0.07599   7.973 2.57e-14 ***
sexmale           -4.71265    0.74065  -6.363 6.64e-10 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 6.067 on 329 degrees of freedom
  (11 observations deleted due to missingness)
Multiple R-squared:  0.8143,    Adjusted R-squared:  0.8126 
F-statistic:   481 on 3 and 329 DF,  p-value: < 2.2e-16

Task 8

What is the unique proportion of variance in the flipper length variable that’s attributable to sex? In other words, how much variance in the outcome is accounted for by sex?

rsq_m1 <- m_massKg_bill_cntr %>% broom::glance() %>% dplyr::pull(r.squared)
rsq_m2 <- m_mass_bill_sex %>% broom::glance() %>% dplyr::pull(r.squared)

# difference in % rounded to 1 dp
round((rsq_m2 - rsq_m1) * 100, 1)

[1] 2.6

Model comparison

Now that we know all about the effect of sex on flipper length in penguins (this is going to win you a pub quiz one day!), let’s look at a formal comparison between a model that includes the variable and one that does not.

Task 9

Let’s formally compare m_massKg_bill_cntr and m_mass_bill_sex.

Task 9.1

Enter the objects containing the two linear models as arguments to the anova() function.

anova(m_massKg_bill_cntr, m_mass_bill_sex)

Well, that didn’t work, did it?

The error message we get from R for once actually tells us what the problem is in plain English:

Error in anova.lmlist(object, ...) : models were not all fitted to the same size of dataset

So the problem is that the two models we are trying to compare are based on datasets of different sizes.

Task 9.2

Pause to think why this might be the case…

So now, to compare our two models, we need to refit m_massKg_bill_cntr using the same subset of peng_data that was used to fit m_mass_bill_sex.

Task 9.3

Create an object in the environment called peng_subset that will allow you to refit m_massKg_bill_cntr such that it is comparable with the other model.

peng_subset <- peng_data %>%
  filter(!is.na(sex))

Task 9.4

Create a model called m_massKg_bill_cntr_compare that is the same as m_massKg_bill_cntr but, this time round, fitted to peng_subset.

m_massKg_bill_cntr_compare <- lm(flipper_length_mm ~ mass_kg_cntr + bill_length_cntr, peng_subset)

Task 9.5

Have a look at coefficients of this new model and compare them to the ones from the original model

coef(m_massKg_bill_cntr_compare)

     (Intercept)     mass_kg_cntr bill_length_cntr 
     200.8593884       13.0173337        0.5440353 

coef(m_massKg_bill_cntr)

     (Intercept)     mass_kg_cntr bill_length_cntr 
     200.9152047       13.0507787        0.5492244 

As you can see, the coefficients are ever so slightly different. This is not that surprising because the data these two versions of the model use are themselves different. What is reassuring is that the discrepancies are very small.

Now that the first model is fitted to the same data as the second one, we can finally compare them.

Task 9.6

Compare the m_massKg_bill_cntr_compare and m_mass_bill_sex models using the anova() function.

anova(m_massKg_bill_cntr_compare, m_mass_bill_sex)

Analysis of Variance Table

Model 1: flipper_length_mm ~ mass_kg_cntr + bill_length_cntr
Model 2: flipper_length_mm ~ mass_kg_cntr + bill_length_cntr + sex
  Res.Df   RSS Df Sum of Sq      F    Pr(>F)    
1    330 13598                                  
2    329 12108  1      1490 40.486 6.635e-10 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

OK, that should have worked now. Let’s break down the output of the anova() function.

First of all and least importantly, we have the heading Analysis of Variance Table. This is slightly confusing terminology as the term analysis of variance (or ANOVA) is oftern reserved for a different context. Let’s just agree to call this the F-test and never mind the heading (or the name of the function, for that matter).

Next, we have a re-statement of the models we’re comparing along with labels Model 1 and Model 2. These labels correspond to the rows of the table below:

  Res.Df   RSS Df Sum of Sq      F    Pr(>F)    
1    330 13598                                  
2    329 12108  1      1490 40.486 6.635e-10 ***

The first row of this table gives us information about Model 1 (in our case m_massKg_bill_cntr_compare), while the second row contains info about Model 2 (m_mass_bill_sex), plus the results of the F-test itself. Remember from lecture 10, that the F-test compares the amount of variance in the outcome variable explained by a model to the variance not explained by the model. This table contains the value of the F-statistic (40.49) and its associated p-value (6.64 × 10⁻¹⁰, which is <.001). On top of that, it also contains all the information to actually calculate F and p. Let’s go over this calculation together.

In the lecture, you learnt that the formula for F is

\[F = \frac{\frac{SS_M}{df_M}}{\frac{SS_R}{df_R}}\] where the model sum of squares (\(SS_M\), Sum of Sq in the table above) is just the total SS minus the residual SS (the latter labelled RSS in the table above). You also learnt that the total sum of squares for a model is just the residual SS of the model that we are comparing it to. So, in our case, \(SS_T\) of our model m_mass_bill_sex is the same as \(SS_R\) of m_massKg_bill_cntr_compare. Reading the table, it all checks out:

     RSS Sum of Sq
1  13598
2  12108      1490

Sum of Sq(row 2) = RSS(1) − RSS(2).

The difference between our two models is just one extra coefficient (for sex) and so the model degrees of freedom is 1. The residual degrees of freedom, as mentioned earlier is \(N - p - 1\), where \(p\) is the total number of predictors in the model: 333 − 3 − 1 = 329.

Now that we have all the numbers we need, we can just plug them into the formula:

\[F = \frac{\frac{SS_M}{df_M}}{\frac{SS_R}{df_R}} = \frac{\frac{1,490}{1}}{\frac{12,108}{329}} \approx\frac{1,490}{36.8} \approx40.4\]

Nice!

So what does the F-test tell us? As explained in the lecture, the F-test tests the null hypothesis that the new model is NOT a significant improvement on the old model. If the probability of getting the F-value we got from the test, or an even larger one, if the null is true is small enough – smaller than our chosen \(\alpha\) level – we reject the null hypothesis.

In our case, the p-value is far less than any conventional \(\alpha\) so we can say that:

Task 10

For your last task, perform an F-test comparing m_massKg_bill_cntr_compare to the null model.

m_null <- lm(flipper_length_mm ~ 1, peng_subset)

# You can also compare all 3 models
anova(m_null, m_massKg_bill_cntr_compare)

Analysis of Variance Table

Model 1: flipper_length_mm ~ 1
Model 2: flipper_length_mm ~ mass_kg_cntr + bill_length_cntr
  Res.Df   RSS Df Sum of Sq      F    Pr(>F)    
1    332 65219                                  
2    330 13598  2     51620 626.35 < 2.2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Well done finishing your last Analysing Data tutorial!

summary(m_massKg_bill_cntr_compare)

Call:
lm(formula = flipper_length_mm ~ mass_kg_cntr + bill_length_cntr, 
    data = peng_subset)

Residuals:
     Min       1Q   Median       3Q      Max 
-20.9869  -4.5868   0.4285   4.8784  14.0977 

Coefficients:
                  Estimate Std. Error t value Pr(>|t|)    
(Intercept)      200.85939    0.35180 570.940  < 2e-16 ***
mass_kg_cntr      13.01733    0.54163  24.034  < 2e-16 ***
bill_length_cntr   0.54404    0.07975   6.822 4.31e-11 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 6.419 on 330 degrees of freedom
Multiple R-squared:  0.7915,    Adjusted R-squared:  0.7902 
F-statistic: 626.3 on 2 and 330 DF,  p-value: < 2.2e-16