In the previous tutorial, we practised fitting a linear model with a single predictor to see if we can say anything about the length of a penguin’s flipper if we know something about its body mass. This time, we will build more sophisticated models with more than one predictor. To do that, we’ll make use of the understanding of multiple-predictor linear models that you gained in lecture 9.

Warm-up quiz

Before we start, let’s warm up by answering a few quiz questions!

As you know, the general equation for the linear model with n predictors is

\[y_i = b_0 + b_1\times x_{1_i} + b_2\times x_{2_i} + b_3\times x_{3_i} + \dots + b_{n-1}\times x_{n-1_i} + b_n\times x_{n_i} + e_i\]

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 1punk!

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 2punk!

In this tutorial, just like in the previous one, we’ll continue using the penguins data set from the palmerpenguins package. You will also need the packages tidyverse, broom, and parameters (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(parameters)

Data

Task 3punk!

Just like last time, 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 4punk!

Remind yourself of what’s in the data 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                                 

Single-predictor model

Last time we looked at how we can use information about a penguin’s body mass to make a more accurate prediction about the length of its flippers than simply guessing the mean flipper length for any given specimen.

Task 5

Let’s take a look at this model again.

Task 5.1punk!

Fit the model to the data using the lm(), saving it’s output as m_mass

m_mass <- lm(flipper_length_mm ~ body_mass_g, peng_data)

# alternatively
m_mass <- peng_data %>%
  lm(formula = flipper_length_mm ~ body_mass_g, data = .)

Task 5.2punk!

Use broom::tidy() to get a neat tibble of model summary statistics for this model. Save it as m_mass_summary

m_mass_summary <- m_mass %>%
  broom::tidy()

Task 5.3Prog-rocK

Inspect the results and answer the following questions.

Task 5.4punk!

Use broom::glance() to get the model fit statistics for m_mass and answer the quiz question.

m_mass %>% broom::glance()

# A tibble: 1 × 12
  r.squared adj.r.squared sigma statistic   p.value    df logLik   AIC
      <dbl>         <dbl> <dbl>     <dbl>     <dbl> <dbl>  <dbl> <dbl>
1     0.759         0.758  6.91     1071. 4.37e-107     1 -1146. 2297.
# … with 4 more variables: BIC <dbl>, deviance <dbl>,
#   df.residual <int>, nobs <int>

m_mass %>%
  broom::glance() %>%
  dplyr::mutate(rsq_perc = r.squared * 100) %>%
  dplyr::pull(rsq_perc) %>%
  round(1)

[1] 75.9

Now that you’ve refreshed your memory of the single-predictor model we built last time, let’s extend it by adding another predictor.

Multiple-predictor model

Imagine we’re also interested in whether or not bill length is related to flipper length. We can extend out m_mass model to explore this research question.

Before we do that, let’s look at what a simple model of flipper length as predicted by bill length tells us.

Task 6Prog-rocK

Fit this model, saving it as m_bill, inspect the results and answer the following questions.

m_bill <- lm(flipper_length_mm ~ bill_length_mm, peng_data)
m_bill_smry <- m_bill %>% broom::tidy()

m_bill_smry

# A tibble: 2 × 5
  term           estimate std.error statistic  p.value
  <chr>             <dbl>     <dbl>     <dbl>    <dbl>
1 (Intercept)      127.       4.67       27.2 3.65e-87
2 bill_length_mm     1.69     0.105      16.0 1.74e-43

m_bill %>%
  broom::glance() %>%
  dplyr::mutate(rsq_perc = r.squared * 100) %>%
  dplyr::pull(rsq_perc) %>%
  round(2)

[1] 43.06

OK, so we know that bill length is a significant predictor or flipper length. However, it could well be the case, that both of these variables are just reflections of body mass and, if we take body mass into account, the relationship between bill and flipper length disappears. In other words, it’s possible that larger penguins have both longer bills and longer flippers but once you only consider birds of a certain body mass, the ones with longer bills do not have longer flippers. Putting both of the predictors in the model allows us to see whether or not this is actually the case.

Let’s do it then!

Task 7

Let’s fit and interpret a two-predictor model of flipper length one step at a time.

Task 7.1punk!

First of all, let’s fit the model using lm() and save it as m_mass_bill.

m_mass_bill <- lm(flipper_length_mm ~ body_mass_g + bill_length_mm, peng_data)

Task 7.2Prog-rocK

Get the summary of the model and have a look at it. Are the model coefficients the same as the ones from m_mass and m_bill? If not, how are they different?

m_mass_bill_smry <- m_mass_bill %>% broom::tidy()
m_mass_bill_smry

# A tibble: 3 × 5
  term           estimate std.error statistic   p.value
  <chr>             <dbl>     <dbl>     <dbl>     <dbl>
1 (Intercept)    122.      2.86         42.7  1.70e-138
2 body_mass_g      0.0131  0.000545     23.9  7.56e- 75
3 bill_length_mm   0.549   0.0801        6.86 3.31e- 11

Task 8jazz...

Given the coefficients of our model – 121.956, 0.013, and 0.549 – give the predicted value of flipper length for the following values of body mass and bill length. (2 decimal places)

A good thing to notice is the relationships between the answers to each of the questions and the \(b\) coefficients.

The values of body mass in questions 14 and 15 are the same and the values of bill length differ by 15 mm. That means that the difference between the results is going to be 15 × 0.549 = 8.24. That’s going to be the same as the difference between answers to questions 16 and 17 because the difference in the respective values of bill length in these two questions are also 15 mm, while the values of body mass are, again, the same.

Questions 14 and 16 on the one hand and 15 and 17 on the other have the same values of bill length but, in both cases, they differ in body mass by 400 g. That means that the differences between both pairs of results is going to be 400 × 0.013 = 5.2.

Finally, the differences between questions 14 and 17 and 15 and 16, respectively, is going to be the sum of the previous two differences: 400 × 0.013 + 15 × 0.549 = 13.44.

m_mass_bill_smry

# A tibble: 3 × 5
  term           estimate std.error statistic   p.value
  <chr>             <dbl>     <dbl>     <dbl>     <dbl>
1 (Intercept)    122.      2.86         42.7  1.70e-138
2 body_mass_g      0.0131  0.000545     23.9  7.56e- 75
3 bill_length_mm   0.549   0.0801        6.86 3.31e- 11

m_mass_bill_smry %>%
  dplyr::mutate(x = c(1, 4200, 45), # step 1
                bx = estimate * x) %>% # step 2
  dplyr::pull(bx) %>% # step 3
  sum() %>% # step 4
  round(2) # step 5

[1] 201.48

Task 9punk!

Look at the model fit statistic, R² in particular.

m_mass_bill %>% broom::glance()

# A tibble: 1 × 12
  r.squared adj.r.squared sigma statistic   p.value    df logLik   AIC
      <dbl>         <dbl> <dbl>     <dbl>     <dbl> <dbl>  <dbl> <dbl>
1     0.788         0.787  6.49      631. 4.88e-115     2 -1123. 2255.
# … with 4 more variables: BIC <dbl>, deviance <dbl>,
#   df.residual <int>, nobs <int>

Transforming variables in the model

When we were reporting our body mass model in the last tutorial, we changed the scale of the predictor from grammes to kilogrammes. This is by no means a problematic thing to do. After all, given the mean and SD of the body_mass_g variable (4201.75 and 801.95, respectively), it is not really reasonable to expect every 1 g increase in body mass to have a substantial effect on flipper length and so working with body mass in kg is fine.

Task 10Prog-rocK

Rescale the body mass variable, refit the model and inspect the results to see how the coefficients change.

# first, create the re-scaled variable
peng_data <- peng_data %>%
  dplyr::mutate(mass_kg = body_mass_g / 1000)
# re-fti the model
m_massKg_bill <- lm(flipper_length_mm ~ mass_kg + bill_length_mm, peng_data)
# see resutls
m_massKg_bill_smry <- m_massKg_bill %>% broom::tidy()
m_massKg_bill_smry

# A tibble: 3 × 5
  term           estimate std.error statistic   p.value
  <chr>             <dbl>     <dbl>     <dbl>     <dbl>
1 (Intercept)     122.       2.86       42.7  1.70e-138
2 mass_kg          13.1      0.545      23.9  7.56e- 75
3 bill_length_mm    0.549    0.0801      6.86 3.31e- 11

You should be able to see that both the intercept and the slope for bill length stayed the same as in the original model but the slope for mass (kg) scaled up by 1,000, compared to the one in the original model. Other than that, it’s the same number.

So as you can see, scaling the predictors scales the coefficients accordingly. If we wanted, we could scale flipper or bill length to some other units (such as centimetres) but, in this case, millimetres work well.

What doesn’t work as well, however, is the intercept. Not that there’s something wrong with it but, in the current model, it’s not very meaningful. After all, the flipper length of a hypothetical massless bill-less penguin is not that interesting. If may be more interesting to know what the expected flipper length is for a penguin of average size and average bill length.

Task 11Prog-rocK

Re-fit the model with the predictors centred around their respective means.

peng_data <- peng_data %>%
  mutate(
    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)
# see resutls
m_massKg_bill_cntr_smry <- m_massKg_bill_cntr %>% broom::tidy()
m_massKg_bill_cntr_smry

# A tibble: 3 × 5
  term             estimate std.error statistic  p.value
  <chr>               <dbl>     <dbl>     <dbl>    <dbl>
1 (Intercept)       201.       0.351     573.   0       
2 mass_kg_cntr       13.1      0.545      23.9  7.56e-75
3 bill_length_cntr    0.549    0.0801      6.86 3.31e-11

Standardised B

The centred model we just made is pretty informative. However, it’s not immediately obvious which of our two predictors has a stronger effect on the outcome. That’s because the two predictors are measured on incomparable scales (kg vs mm). This is where standardised coefficients (Bs) come in handy: because they are both on the scale of standard deviations, they allow us to compare effects of very different variables.

Task 12

Use parameters::model_parameters() with the standardize = "refit" argument to get B coefficients for the predictors in our model.

m_massKg_bill_cntr %>% parameters::model_parameters(standardize = "refit")

Parameter        | Coefficient |   SE |        95% CI |    t(339) |      p
--------------------------------------------------------------------------
(Intercept)      |   -8.96e-16 | 0.02 | [-0.05, 0.05] | -3.59e-14 | > .999
mass kg cntr     |        0.74 | 0.03 | [ 0.68, 0.81] |     23.94 | < .001
bill length cntr |        0.21 | 0.03 | [ 0.15, 0.27] |      6.86 | < .001

Reporting

Task 13jazz...

Report the findings of the model. At the very least, you should report R² and the individual model coefficients along with their test of significance.

That’s all for today. Good job!