Recall the general linear model, y = \beta_0 + \beta_1 x_1 + ... + \beta_k x_k + \varepsilon
In the last lecture, we introduced the concept of categorical predictors.
We also discussed how to deal with categorical variables either as a factor variable or as a set of indicator variables.
duck_incidents <- read_csv("https://raw.githubusercontent.com/samanthaseals/SDSII/refs/heads/main/files/data/lectures/W2_duck_incidents.csv") %>%
mutate(loc_yard = if_else(location == "Backyard", 1, 0),
loc_garage = if_else(location == "Garage", 1, 0),
loc_kitchen = if_else(location == "Kitchen", 1, 0),
loc_living = if_else(location == "Living Room", 1, 0))
duck_incidents %>% head()Recall that for a categorical variable with c categories, we include c-1 indicator variables in the model.
The category that is not included in the model is called the reference group.
How does R handle reference groups when modeling?
If we include a factor variable in the model, R automatically creates the indicator variables and selects the reference group (default is the “first” level).
If we include indicator variables in the model, we must ensure that we include only c-1 of them.
The cth one is the reference group.
If we include all c, R will not return an estimate for one of the categorical terms.
| Nephew | x_{\text{H}} | x_{\text{D}} | x_{\text{L}} |
|---|---|---|---|
| Huey | 1 | 0 | 0 |
| Dewey | 0 | 1 | 0 |
| Louie | 0 | 0 | 1 |
We only need two indicators to represent the three categories.
Suppose Louie is the reference group. That means we include x_{\text{H}} and x_{\text{D}} in the model.
When x_{\text{H}} = 1 and x_{\text{D}} = 0, we know nephew = \text{Huey}.
When x_{\text{H}} = 0 and x_{\text{D}} = 1, we know nephew = \text{Dewey}.
When x_{\text{H}} = 0 and x_{\text{D}} = 0, we know nephew = \text{Louie}.
Now, let’s consider this from the mathematical standpoint.
If we included all three indicators, we would have the following in our model:
y = \beta_0 + \beta_1 x_{\text{H}} + \beta_2 x_{\text{D}} + \beta_3 x_{\text{L}} + \varepsilon
y = \beta_0 + \beta_1 x_{\text{H}} + \beta_2 x_{\text{D}} + \beta_3 (1 - x_{\text{H}} - x_{\text{D}}) + \varepsilon
y = (\beta_0 + \beta_3) + (\beta_1 - \beta_3) x_{\text{H}} + (\beta_2 - \beta_3) x_{\text{D}} + \varepsilon
How do we change reference groups for modeling/testing purposes?
If using a factor variable, we can relevel the factor using the factor() function with the levels argument.
If using indicator variables, we simply update the set of c-1 indicators in the model.
Prepare yourself to see a lot of “if you are using a factor variable… but if you are using indicator variables…” throughout the course.
What is my approach (in practice) to choosing a reference group?
If there is some natural order, I select the lowest as the reference.
If there is not an order, I will select one that “makes sense” as the comparison group when thinking of the context surrounding the question being asked of the data.
When I am communicating results to collaborators, I will explain my choice and then ask if they would like a different reference group.
To actually model with categorical predictors in R, we use the appropriate modeling function as usual.
If we are using stored-as-factor variables, we simply include them,
as.factor(),m1 <- glm(damage_cost ~ location,
family = "gaussian",
data = duck_incidents)
tidy(m1, conf.int = TRUE)m2 <- glm(damage_cost ~ loc_garage + loc_kitchen + loc_living,
family = "gaussian",
data = duck_incidents)
tidy(m2, conf.int = TRUE) (Intercept) locationGarage locationKitchen locationLiving Room
202.77983 64.34182 131.84687 58.69676 \hat{\text{cost}} = 202.78 + 64.34 \text{ garage} + 131.85 \text{ kitchen} + 58.70 \text{ living}
From here, we can easily create the estimated costs for each location:
We know that for continuous predictors, the interpretation is “for a one unit increase in x_i, we expect the outcome to [increase or decrease] by \beta_i units, holding all other predictors constant.”
For categorical predictors, what is a one unit increase in x_j?
For an indicator variable, a one unit increase means going from not being in that category to being in that category.
For a factor variable, a one unit increase means going from the reference group to the category represented by that term.
Thus, the interpretation for categorical predictors becomes “for those in category j, as compared to the reference group, we expect the outcome to [increase or decrease] by \beta_j units, holding all other predictors constant.”
(Intercept) locationGarage locationKitchen locationLiving Room
202.77983 64.34182 131.84687 58.69676 We can interpret the model estimates as follows:
The average cost of damage for incidents that occurred in the Backyard is $202.78.
For incidents that occurred in the Garage, we expect the damage cost to increase by $64.34.
For incidents that occurred in the Kitchen, we expect the damage cost to increase by $131.85.
For incidents that occurred in the Living Room, we expect the damage cost to increase $58.70.
We have said that the choice of reference group is arbitrary mathematically. Let’s prove this to ourselves.
Let’s now model damage cost using location but with “Kitchen” as the reference group.
m3 <- glm(damage_cost ~ loc_yard + loc_garage + loc_living,
family = "gaussian",
data = duck_incidents)
tidy(m3, conf.int = TRUE)We see that the estimates themselves have changed.
(Intercept) loc_garage loc_kitchen loc_living
202.77983 64.34182 131.84687 58.69676 (Intercept) loc_yard loc_garage loc_living
334.62670 -131.84687 -67.50505 -73.15011 loc_kitchen was $131.85 – this is the difference between the backyard and the kitchen averages.
(Intercept) locationGarage locationKitchen locationLiving Room
139.288148 51.491118 113.783352 66.873723
nephewHuey nephewLouie
-7.974478 206.606636 \begin{align*} \hat{\text{cost}} = 139.29 &+ 51.49 \text{ garage} + 113.78 \text{ kitchen} + 66.87 \text{ living} \\ & - 7.97 \text{ Huey} + 206.61 \text{ Louie} \end{align*}
We now have two reference groups!
\begin{align*} \hat{\text{cost}} = 139.29 &+ 51.49 \text{ garage} + 113.78 \text{ kitchen} + 66.87 \text{ living} \\ & - 7.97 \text{ Huey} + 206.61 \text{ Louie} \end{align*}
The intercept represents the expected damage cost for incidents that occurred in the backyard with Dewey present.
Locations:
Nephews:
(Intercept) locationGarage locationKitchen locationLiving Room
62.162513 54.773465 92.647154 58.373882
nephewHuey nephewLouie sugar_grams
-0.899301 198.792202 1.835897 \begin{align*} \hat{\text{cost}} = 62.16 &+ 54.77 \text{ garage} + 92.65 \text{ kitchen} + 58.37 \text{ living} \\ & - 0.90 \text{ Huey} + 198.79 \text{ Louie} \\ & + 1.84 \text{ sugar} \end{align*}
\begin{align*} \hat{\text{cost}} = 62.16 &+ 54.77 \text{ garage} + 92.65 \text{ kitchen} + 58.37 \text{ living} \\ & - 0.90 \text{ Huey} + 198.79 \text{ Louie} \\ & + 1.84 \text{ sugar} \end{align*}
The interpretation of predictors remains the same as we have previously learned.
punished binary indicator we created that identifies if Donald punishes the nephews (either “Assigns Chores” or “Grounds”) or not (either “Laughs” or “Yells”).donald_action factor variable instead,\hat{\text{cost}} = 86.61 + 3.80 \text{ sugar} + 8.28 \text{ punished}
Interpreting our model coefficients,
The average cost of damage when Donald does not punish the nephews and no sugar is involved is $86.61.
For each additional gram of sugar intake, we expect the damage cost to increase by $3.80.
When Donald punishes the nephews, we expect the damage cost to increase by $8.28.
In this lecture, we have introduced modeling with categorical predictors.
Always remember that if there are c categories, there will be c-1 predictor terms in the model.
Whether we use a true factor variable or we use a set of indicator variables, estimation, model fit, and omnibus test results will remain the same.
The choice of reference group is arbitrary mathematically, but important for interpretation.
In the next lecture, we will discuss statistical inference with categorical predictors.