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

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, color = species)) +
  geom_point()

## Warning: Removed 2 rows containing missing values or values outside the scale range
## (`geom_point()`).

Describe what you see:

adelie penguins form a cluster with shorter bill lengths and moderate bill depths, while chinstrap and gentoo penguins form clusters with longer bill lengths. the slopes appear somewhat different, but the biggest difference is in intercepts — adelies start lower on the y‑axis. yhere may also be a slight slope difference, suggesting an interaction term could matter..

2B: Creating the Indicator and Interaction Variables

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

PenguinsAdelie <- mutate(
  PData,
  AdelieIND = ifelse(species == "Adelie", 1, 0),
  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.BL.BD.Adelie <- lm(
  bill_length_mm ~ bill_depth_mm + AdelieIND + AdelieINT,
  data = PenguinsAdelie
)
mlm.BL.BD.Adelie

## 
## Call:
## lm(formula = bill_length_mm ~ bill_depth_mm + AdelieIND + AdelieINT, 
##     data = PenguinsAdelie)
## 
## Coefficients:
##   (Intercept)  bill_depth_mm      AdelieIND      AdelieINT  
##      34.42413        0.83634      -11.35605        0.02069

Write out the full regression equation:

y-hat = 20.96 + 1.76 ( bill depth ) − 9.77 ( AdelieIND ) − 0.51 ( AdelieINT )

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 = 20.96 + 1.76 ( bill depth )

Group 2 (indicator = 1):

y-hat = 20.96+1.76𝑥−9.77(1)−0.51(1⋅𝑥)

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 = 20.96, slope = 1.76, color = "blue") +      # non-Adelie
  geom_abline(intercept = 11.19, slope = 1.25, color = "red")         # Adelie

## 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?

the adelie line sits much lower, matching their shorter bill lengths. the slopes differ slightly (1.76 vs. 1.25), so the lines are not parallel — they gradually move apart as bill depth increases. this matches what we saw in the scatterplot: adelies form a distinct cluster with a different trend.

Unit 13 Project

addison duke

4/23/26