Multiple Regression with Indicator and Interaction Variables

In this project we practice building multiple regression models that include indicator and interaction variables. The project has two parts. In Part 1 we work through a complete example together using the Palmer Penguins dataset. In Part 2 you carry out your own analysis from start to finish.

Part 1: Worked Example — Gentoo vs. Non-Gentoo Penguins

In an earlier script we noticed that a scatterplot of flipper length vs. bill depth showed two distinct “clumps” of points. We identified species — and specifically Gentoo penguins — as the likely cause. Here we build a multiple regression model that accounts for this.

1A: The Motivating Scatterplot

Let’s re-examine the scatterplot that showed the two clumps, coloring by species so we can see what is going on.

ggplot(PData, aes(x = flipper_length_mm, y = bill_depth_mm)) +
  geom_point()
## Warning: Removed 2 rows containing missing values or values outside the scale range
## (`geom_point()`).

ggplot(PData, aes(x = flipper_length_mm, y = bill_depth_mm, color = species)) +
  geom_point()
## Warning: Removed 2 rows containing missing values or values outside the scale range
## (`geom_point()`).

The Gentoo penguins are clearly behaving differently from the Adelie and Chinstrap penguins! The two groups have different slopes AND different intercepts. This tells us we need both an indicator variable and an interaction variable in our model.

1B: Creating the Indicator and Interaction Variables

We use mutate() to add two new columns to the data. The indicator variable (which we name GentooIND) equals 1 for Gentoo penguins and 0 for all others. The interaction variable is the product of the indicator and the predictor (in this case: GentooIND times flipper length).

PenguinDataGentoo <- mutate(PData, GentooIND = ifelse(species == "Gentoo", 1, 0))
PenguinDataGentoo <- mutate(PenguinDataGentoo, GentooINT = GentooIND * flipper_length_mm)

1C: Fitting the Multiple Regression Model

We predict bill depth from flipper length, the Gentoo indicator, and the interaction variable.

mlm.BDepth.FLength.Ind.Int <- lm(PenguinDataGentoo$bill_depth_mm ~
                                   PenguinDataGentoo$flipper_length_mm +
                                   PenguinDataGentoo$GentooIND +
                                   PenguinDataGentoo$GentooINT)
mlm.BDepth.FLength.Ind.Int
## 
## Call:
## lm(formula = PenguinDataGentoo$bill_depth_mm ~ PenguinDataGentoo$flipper_length_mm + 
##     PenguinDataGentoo$GentooIND + PenguinDataGentoo$GentooINT)
## 
## Coefficients:
##                         (Intercept)  PenguinDataGentoo$flipper_length_mm  
##                             6.59425                              0.06140  
##         PenguinDataGentoo$GentooIND          PenguinDataGentoo$GentooINT  
##                           -14.83111                              0.04551

Our coefficients are:

  • Intercept = 6.39
  • flipper_length_mm = 0.06
  • GentooIND = -14.50
  • GentooINT = 0.044

So our full equation is: Bill Depth = 6.39 + 0.06(Flipper Length) − 14.50(Indicator) + 0.044(Interaction)

1E: Recovering the Two Group Equations

For non-Gentoo penguins (Indicator = 0):

Bill Depth = 6.39 + 0.06(Flipper Length) − 14.50(0) + 0.044(0 × Flipper Length)

Bill Depth = 6.39 + 0.06(Flipper Length)

For Gentoo penguins (Indicator = 1):

Bill Depth = 6.39 + 0.06(Flipper Length) − 14.50(1) + 0.044(1 × Flipper Length)

Bill Depth = (6.39 − 14.50) + (0.06 + 0.044)(Flipper Length)

Bill Depth = −8.11 + 0.104(Flipper Length)

1E: Plotting the Two Lines of Best Fit

ggplot(PData, aes(x = flipper_length_mm, y = bill_depth_mm, color = species)) +
  geom_point() +
  geom_abline(intercept = 6.39,  slope = 0.06,  color = "red") +
  geom_abline(intercept = -8.11, slope = 0.104, color = "blue")
## Warning: Removed 2 rows containing missing values or values outside the scale range
## (`geom_point()`).

The red line describes non-Gentoo penguins and the blue line describes Gentoo penguins. The two lines fit their respective clusters well, and they are clearly not parallel — confirming that the interaction term was necessary. —

Part 2: Your Turn

Choose one of the following two options. For whichever option you choose, carry out the full analysis described in the sections below.

Option 1 (Penguins — species indicator for Adelie): Use the full PData dataset. Create an indicator variable for Adelie penguins (AdelieIND = 1 if Adelie, 0 otherwise). Predict bill_length_mm from bill_depth_mm, your indicator, and the interaction variable AdelieINT = AdelieIND * bill_depth_mm.

Option 2 (Penguins — species indicator for Chinstrap): Use the full PData dataset. Create an indicator variable for Chinstrap penguins (ChinstrapIND = 1 if Chinstrap, 0 otherwise). Predict bill_length_mm from body_mass_g, your indicator, and the interaction variable ChinstrapINT = ChinstrapIND * body_mass_g.

State which option you chose:

Option 1. Option 1 (Penguins — species indicator for Adelie): Use the full PData dataset. Create an indicator variable for Adelie penguins (AdelieIND = 1 if Adelie, 0 otherwise). Predict bill_length_mm from bill_depth_mm, your indicator, and the interaction variable AdelieINT = AdelieIND * bill_depth_mm. —

2A: The Motivating Scatterplot

Make a scatterplot of your two quantitative variables (x = predictor, y = response), coloring the points by your grouping variable. Do the two groups appear to have different slopes, different intercepts, or both?

ggplot(PData, aes(x =  bill_depth_mm, y =  bill_length_mm)) +
  geom_point()
## Warning: Removed 2 rows containing missing values or values outside the scale range
## (`geom_point()`).

ggplot(PData, aes(x = bill_depth_mm, y = bill_length_mm, color = species)) +
  geom_point()
## Warning: Removed 2 rows containing missing values or values outside the scale range
## (`geom_point()`).

Describe what you see: The Adelie penguins seem to have different data than the Gentoo and Chinstrap data sets. The longer the bill depth the shorter the bill length is. Gentoo have a shorter bill depth but a larger bill length and Chinstrap have a bill depth and bill length thats between the Gentoo and Adelie data sets.*

2B: Creating the Indicator and Interaction Variables

Option 1 (Penguins — species indicator for Adelie): Use the full PData dataset. Create an indicator variable for Adelie penguins (AdelieIND = 1 if Adelie, 0 otherwise). Predict bill_length_mm from bill_depth_mm, your indicator, and the interaction variable AdelieINT = AdelieIND * bill_depth_mm. —

Use mutate() and ifelse() to create your indicator variable and your interaction variable.

PenguinDataAdelie <- mutate(PData, AdelieIND = ifelse(species == "Adelie", 1, 0))
PenguinDataAdelie <- mutate(PenguinDataAdelie, AdelieINT = AdelieIND * bill_depth_mm)

2C: Fitting the Multiple Regression Model

Fit your model using lm(). Report the four coefficients (intercept, predictor slope, indicator, interaction).

mlm.BLength.BDepth.Ind.Int <- lm(PenguinDataAdelie$bill_length_mm ~
                                   PenguinDataAdelie$bill_depth_mm+
                                   PenguinDataAdelie$AdelieIND +
                                   PenguinDataAdelie$AdelieINT)
mlm.BLength.BDepth.Ind.Int
## 
## Call:
## lm(formula = PenguinDataAdelie$bill_length_mm ~ PenguinDataAdelie$bill_depth_mm + 
##     PenguinDataAdelie$AdelieIND + PenguinDataAdelie$AdelieINT)
## 
## Coefficients:
##                     (Intercept)  PenguinDataAdelie$bill_depth_mm  
##                        34.42413                          0.83634  
##     PenguinDataAdelie$AdelieIND      PenguinDataAdelie$AdelieINT  
##                       -11.35605                          0.02069

Write out the full regression equation:

y-hat = .83634(predictor) + -11.35605(indicator) + .02069 (interaction) + 34.42413

2D: Recovering the Two Group Equations

Use the four coefficients to work out the prediction equation for each group, the same way we did in Section 1E. Show your algebra.

Group 1 (reference category, indicator = 0):

y-hat = .83634(PenguinDataAdelie\(bill_depth_mm) + -11.35605(0) + .02069 (0*PenguinDataAdelie\)bill_depth_mm) + 34.42413

y-hat= .83634(PenguinDataAdelie$bill_depth_mm) + 34.42413

Group 2 (indicator = 1):

y-hat = .83634(PenguinDataAdelie\(bill_depth_mm) + -11.35605(1) + .02069 (1*PenguinDataAdelie\)bill_depth_mm) + 34.42413

y-hat = .85703(PenguinDataAdelie$bill_depth_mm) + 23.06808

2E: Plotting the Two Lines of Best Fit

Use geom_abline() to draw both lines on top of your color-coded scatterplot. Use the group equations from 2D for the slope and intercept of each line.

ggplot(PData, aes(x = bill_depth_mm, y = bill_length_mm, color = species)) +
  geom_point() +
  geom_abline(intercept = 34.42413,  slope = .83634,  color = "hotpink") +
  geom_abline(intercept = 23.06808, slope = .85703, color = "purple")
## Warning: Removed 2 rows containing missing values or values outside the scale range
## (`geom_point()`).

Do the lines fit their respective groups well? Are they parallel or do they cross? Does this match what you expected from your scatterplot in 2A?

Yes! The lines fit their respective groups well, they seem to cutt across the data sets accordingly at a mid point. The lines of best fit are parallel and match the scatterplot in 2A!