---
title: "Categorical predictors (Notes)"
subtitle: "STAT 155"
author: "Your Name"
format:
  html:
    toc: true
    toc-depth: 2
    embed-resources: true
---


```{r setup}
#| include: false
knitr::opts_chunk$set(
  collapse = TRUE, 
  warning = FALSE,
  message = FALSE,
  error = TRUE,
  fig.height = 2.75, 
  fig.width = 4.25,
  fig.env = 'figure',
  fig.pos = 'h',
  fig.align = 'center')
```



\


::: {.callout-note title = "Learning goals"}

- Write a model formula for a simple linear regression model with a **categorical predictor** using **indicator variables**.
- Interpret the model coefficients in a simple linear regression model with a categorical predictor.

:::


\




::: {.callout-note title="Additional resources"}

**Required video**

[Simple linear regression: categorical predictor](https://youtu.be/QPpp6Kpk2To) ([slides](https://drive.google.com/file/d/1fbSJW4DhVWxR8rykowPUoAZ6WfzuZ857/view?usp=sharing))


**Optional**

- Watch: [R Code for Categorical Predictors](https://voicethread.com/share/14980390/)
- Read: Section 3.9 in the [STAT 155 Notes](https://mac-stat.github.io/Stat155Notes/) only up through section 3.9.1.

:::




\
\
\
\



# Warm-up {-}

Today we'll explore penguins!

```{r}
# Load packages & data
# The penguins data is stored in R
library(tidyverse)
data(penguins)
penguins <- penguins %>% 
  filter(!is.na(sex), !is.na(species))
```

You can pull up the codebook by typing the following in your **CONSOLE**: ?penguins

```{r}
# The cases / units of observation are penguins
# We have data on 333 penguins, 8 variables on each
head(penguins)
dim(penguins)
```

We'll focus on 3 variables:

- `flipper_len` = the length of the penguin's flippers (arms) in mm
- `species` = Adelie, Chinstrap, or Gentoo
- `sex` = female, male

\
\


**EXAMPLE 1: Get to know the data**

```{r}
# Create a dataset called "chinstraps" that only includes Chinstrap penguins
chinstraps <- penguins %>% 
  ___(species == "Chinstrap")
head(chinstraps)
```

```{r}
# Among Chinstraps, what's the average flipper length (mm)?
chinstraps %>% 
  ___(mean(flipper_len))
```

```{r}
# Among Chinstraps, how many female and male penguins are there?
# (We haven't used this code in a while!)
chinstraps %>% 
  ___(sex)
```


\
\
\
\


**EXAMPLE 2: Flipper length vs sex (bad viz)**

Among *Chinstraps*, let's explore how flipper length might be explained by sex.
Use the results below to explain why a scatterplot is not an effective visualization and a line is not an effective summary of the relationship in this example (or for exploring a quantitative Y variable vs categorical X predictor in general).

```{r}
chinstraps %>% 
  ggplot(aes(y = flipper_len, x = sex)) + 
  geom_point() + 
  geom_smooth(method = "lm", se = FALSE)
```


\
\
\
\


**EXAMPLE 3: Flipper length vs sex (good viz)**

Side-by-side *density plots* provide a better picture of the relationship between flipper length and sex.

```{r}
# Check out a density plot of `flipper_len` alone
chinstraps %>% 
  ggplot(aes(x = flipper_len)) + 
  geom_density()
```

```{r}
# Modify that code to separate the density plot by sex,
# using different line colors to distinguish females and males:
chinstraps %>% 
  ggplot(aes(x = flipper_len)) + 
  geom_density()
```

```{r}
# Modify the code again: use different *fills* to distinguish females and males
chinstraps %>% 
  ggplot(aes(x = flipper_len)) + 
  geom_density()
```

```{r}
# Modify this code to create side-by-side *boxplots* of flipper_len by sex
chinstraps %>% 
  ggplot(aes(y = flipper_len)) + 
  geom_boxplot()
```

```{r}
# Modify this code to create separate *histograms*:
chinstraps %>% 
  ggplot(aes(x = flipper_len)) + 
  geom_histogram(color = "white") + 
  facet_wrap(~ ___)
```


**FOLLOW-UPS**

- Describe the relationship of `flipper_len` with `sex`.
- What do you think was the least effective plot?


\
\
\
\


::: {.callout-note title = "Next up: modeling!"}

We can use linear regression to model the relationship of a quantitative response variable Y with a categorical predictor variable X!
This is a relief because predictors are often categorical.
Examples:

- Model time to campus (Y) by mode of transportation (X) -- walk, bike, bus, etc.
- Model film revenue (Y) by genre (X) -- rom com, sci-fi, etc.
- Model wages (Y) by job industry (X) -- agriculture, construction, etc.
- Model some health outcome (Y) by treatment -- A or B.

Though the linear regression model formula is *not* represented by a line, it still has a familiar form:

$E[Y | X]$ = $\beta_0$ + $\beta_1$ (something related to $X$)

:::





\
\
\
\





# Exercises {-}



## Exercise 1: Summaries by category {-}

In exploring the relationship of `flipper_len` with `sex`, we've explored some visual summaries.
Let's follow these up with some numerical summaries.
Recall that among all Chinstraps, males and females *combined*, the average flipper length was 195.8mm:

```{r}
chinstraps %>% 
  summarize(mean(flipper_len))
```

We can calculate separate averages for female and male penguins by adding the `group_by()` function:

```{r}
# PAUSE: Note where the group_by() function is, and what information it uses!
chinstraps %>% 
  group_by(sex) %>% 
  summarize(mean(flipper_len))
```

### Part a {-}

Examine the group means, store them in your brain for later, and make sure that you can match these numbers up with what we observed in the plots.

### Part b {-}

In the chunk below, calculate the *difference* between the mean flipper length for male and female penguins (mean for males - mean for females).
Check out the actual number and keep this in your brain!

```{r}

```

### Part c {-}

In general, the `group_by()` is a unique function in that we'd almost never use it on its own.
All it does is change the background structure of the dataset:

```{r}
# Check out the head of the data set
chinstraps %>% 
  head()
```

```{r}
# Check out the head of the data AFTER using group_by
# How did this change?
chinstraps %>% 
  group_by(sex) %>% 
  head()
```

So we *almost always* follow up `group_by()` with `summarize()`:

```{r}
# Example: calculate the minimum flipper length by sex
chinstraps %>% 
  group_by(sex) %>% 
  summarize(min(flipper_len))
```




\
\
\
\


## Exercise 2: modeling `flipper_len` by `sex` {-}

Now that we've visually and numerically explored the relationship of `flipper_len` and `sex` for Chinstrap penguins, let's build a simple linear regression model of this relationship:

```{r}
# Construct the model
flipper_model_1 <- lm(flipper_len ~ sex, data = chinstraps)

# Summarize the model
summary(flipper_model_1)
```

Fill in the resulting model formula:

E[`flipper_len` | `sex`] = ___ + ___ `sexmale`


IMPORTANT: Discuss the following before moving on.

- Where have you encountered the 2 model coefficients before?!

- Notice that `sex` shows up as `sexmale` in the model! This is a 0/1 **indicator variable** of whether a penguin is *male*:    
    - `sexmale` = 1 if the penguin is male
    - `sexmale` = 0 otherwise (i.e the penguin is female)

- In contrast, we do not observe `sexfemale` in the model output! This is because `female` is what we call the **reference level** or **reference category** of the `sex` variable. It was chosen to be the reference level simply because, *alphabetically*, it's the first `sex` category. You’ll observe below that female penguins are, indeed, still in the model. You’ll also observe why the term "reference level" makes sense!



\
\
\
\


## Exercise 3: Making sense of the model {-}

Recall our model formula:

E[`flipper_len` | `sex`] = 191.735 + 8.176 `sexmale`

### Part a {-}

Use the above model formula to calculate the expected flipper length for a *female* penguin. HINT: What is `sexmale` for female penguins, 0 or 1?

```{r}
191.735 + 8.176 * ___
```

### Part b {-}

Similarly, calculate the expected flipper length for a *male* penguin.

```{r}
191.735 + 8.176 * ___
```

### Part c {-}

Re-examine these 2 calculations. Where have you come across these numbers before?!





\
\
\
\


## Exercise 4: Interpreting coefficients {-}

Recall that our model formula is *not a formula for a line*:

E[`flipper_len` | `sex`] = 191.735 + 8.176 `sexmale`

Thus we can't interpret the coefficients as we have before when using a quantitative predictor.
Taking this into account and reflecting upon your calculations above...   

### Part a {-}

Interpret the intercept coefficient (`191.735`) in the penguin context.
*Make sure to use non-causal language, include units, and talk about averages rather than individual cases.*


### Part b {-}

Interpret the `sexmale` coefficient (`8.176`) in context.
Hint: How did you use this value in the prediction calculations above?





\
\
\
\


## Exercise 5: The least squares criterion {-}

The coefficient estimates in this model were selected in the same way as when our predictor is quantitative: by minimizing the sum of the squared residuals.

E[`flipper_len` | `sex`] = 191.735 + 8.176 `sexmale`

Use this model to calculate the residual for the first penguin in the dataset:

```{r}
# Observed data on the first penguin
chinstraps %>% 
  select(flipper_len, sex) %>% 
  head(1)
```



\
\
\
\


## Exercise 6: Categorical predictors with >2 categories {-}

Let's zoom back out to all `penguins`, not just Chinstraps, and explore how `flipper_len` can be modeled by `species`.

```{r}
# Count up the number of penguins of each species
# (How many species are there?)
penguins %>% 
  
```

```{r}
# Construct a side-by-side density plot of the relationship between flipper_len and species

```

```{r}
# Calculate the average flipper_len for each species
# Do this without using any filter()!!!

```


**FOLLOW-UP:**

- Describe what you've learned so far about the relationship of `flipper_len` with `species`. Be sure to comment on all 3 species!
- What is the *reference level* of the `species` variable?


\
\
\
\


## Exercise 7: Build & interpret the model {-}

Let's now model this relationship:

```{r}
# Construct the model
flipper_model_2 <- lm(flipper_len ~ species, data = penguins)

# Summarize the model
summary(flipper_model_2)
```

### Part a {-}

Complete the model formula below.
THINK: Why don't Adelie penguins show up here? 

E[`flipper_len` | `species`] = ___ + ___ speciesChinstrap + ___ speciesGentoo

The `species` related predictors are, again, 0/1 indicator variables: 

- `speciesChinstrap` = 1 for Chinstrap penguins and 0 otherwise (for Adelie or Gentoo)
- `speciesGentoo` = 1 for Gentoo penguins and 0 otherwise (for Adelie or Chinstrap)


### Part b {-}

Use the model to make some predictions.
Use your model formula above, not the `predict()` function.

```{r}
# Predict flipper_len for an Adelie penguin


# Predict flipper_len for a Chinstrap penguin


# Predict flipper_len for a Gentoo penguin


```


### Part c {-}

Interpret all 3 model coefficients: 190.1027, 5.7208, 27.1326.
*Make sure to frame this is context, use non-causal language, include units, and talk about averages rather than individual cases.*


### Part d {-}

How *strong* is this model?
Since `species` is a *categorical* predictor, we *cannot* use correlation to answer this question -- correlation can only be used to measure the strength and direction of a linear relationship between 2 *quantitative* variables.
BUT we *can* use R-squared.

1. Report and interpret the R-squared value for `flipper_model_2`.


2. Based on this and the fact that `flipper_model_1` had an R-squared of 0.33, comment on which is the stronger predictor of `flipper_len`: `species` or `sex`



\
\
\
\


## Exercise 8: Really understand how the coefficients work {-}

Next, consider the model of `flipper_len` by `island`.
The average `flipper_len` of penguins on each `island` is calculated below:

```{r}
penguins %>% 
  group_by(island) %>% 
  summarize(mean(flipper_len))
```

### Part a {-}

Using just the above averages and your understanding of categorical predictors, fill in the model coefficients and indicator variables below.
Do NOT use `lm()` yet!!!

E[`flipper_len` | `island`] = ___ +/- ___ island??? +/- ___ island???


### Part b {-}

Check your work to Part a.

```{r}
flipper_model_3 <- lm(flipper_len ~ island, penguins)
summary(flipper_model_3)
```



\
\
\
\




## Exercise 9: Diamond price vs color {-}

Load data on `diamonds` that's stored in the `tidyverse` package:

```{r}
data("diamonds")

# Wrangle the data (don't worry about the code)
diamonds <- diamonds %>% 
    mutate(
        cut = factor(cut, ordered = FALSE),
        color = factor(color, ordered = FALSE),
        clarity = factor(clarity, ordered = FALSE)
    )

head(diamonds)
```

We'll explore the relationship between 2 variables:

- `price` of a diamond in USD
- `color` of a diamond, categorized as D (best) through J (worst)

Practice the below steps in order to explore the relationship of diamond `price` with `color`.

- Before exploring any data, write down what you expect a plot of `price` vs `color` to look like. Then construct and interpret an appropriate plot.

- Compute the average price for each color. (Remembering the color categories, why is this surprising? What do you think is going on here?)

- Fit an appropriate linear model with `lm()` and display a model summary.

- Write out the model formula from the above summary.

- Which color is the reference level? How can you tell from the model summary?

- Interpret the intercept and two other coefficients from the model in terms of the data context.




\
\
\
\



# Reflection {-}

Through the exercises above, you learned how to build and interpret models that incorporate a categorical predictor variable.
For the benefit of your future self, summarize how one can interpret the coefficients for a categorical predictor.







## Done!

- Finalize your notes: (1) Render your notes to an HTML file; (2) Inspect this HTML in your Viewer -- check that your work translated correctly; and (3) Outside RStudio, navigate to your 'Activities' subfolder within your 'STAT155' folder and locate the HTML file. Open and inspect this in a web browser.
- Clean up your RStudio session: End the rendering process by clicking the 'Stop' button in the 'Background Jobs' tab in the lower left pane.
- Check the wrap-up announcements and solutions in the online course manual, at the bottom of the corresponding chapter.
- Work on homework!


