Notes on DataCamp’s course: Multiple and Logistic Regression

Adding a categorical variable to a regression

I am going to run through the example that I gave above about the model of house prices against square footage and number of bathrooms. I’ve used the test.csv from Kaggle here.

Number of bathrooms is a categorical variable which means that the addition of the variable will change the intercept of the model, leading geometrically to parallel lines in two dimensions.

Here is a model of home price as a function of square footage alone. In this model we are not allowing price to differ based on the number of bathrooms in the house.

# load in data
housing_data <- read.csv("test.csv")

# Combine GrLivArea and TotalBsmtSF to get total square footage, add FullBath
#   with BsmtFullBath to ge bathroom_count, and remove houses with 0 or 6
#   bathrooms due to lack of data
housing_data <- housing_data %>% 
  mutate(square_footage = GrLivArea + TotalBsmtSF,
         bathroom_count = as.factor(FullBath + BsmtFullBath)) %>% 
  filter(bathroom_count %in% c(1, 2, 3, 4, 5))

single_var_mod <- lm(SalePrice ~ square_footage, data=housing_data)

ggplot(housing_data, aes(square_footage, SalePrice)) +
  geom_point(alpha=.5) +
  geom_smooth(method = "lm", se = FALSE)

and here is the model that adds the categorical variable bathroom_count to the model. In this case that means adding another two lines to the graphs where each line has the same slope, but differing intercepts.

two_var_mod <- lm(SalePrice ~ square_footage + bathroom_count, data=housing_data)

# augment the model so that I can get .fitted values for the plot
aug_two_var_mod <- augment(two_var_mod)

ggplot(aug_two_var_mod, aes(square_footage, SalePrice, color=bathroom_count)) +
  geom_point(alpha=.5) +
  geom_line(aes(y=.fitted)) +
  facet_wrap(~bathroom_count)

When multiplying explanatory variables together I get what is called an interaction term. The case of linear regression with one numeric variable and one categorical variable leads geometrically to multiple linear regression lines with differing slopes. This means that square_footage is somehow modulated by bathroom_count, this is logical since as the bathroom_count increases we would expect square_footage to increase as well.

interaction_mod <- lm(
  SalePrice ~ square_footage + 
    bathroom_count + 
    square_footage * bathroom_count, 
    data=housing_data
  )

aug_interaction_mod <- augment(interaction_mod)

ggplot(aug_interaction_mod, 
       aes(square_footage, SalePrice, color=bathroom_count)) +
  geom_point(alpha=.5) +
  geom_line(aes(y=.fitted)) +
  facet_wrap(~bathroom_count)

Which model is better?

As one measure of model fit, I can compare the adjusted R² between the models to see which explains more of the total variation in the model. I use the adjusted R² because it penalizes the use of additional varibles as opposed to R² which always increases with the addition of more variables.

model_comparison <- data.frame(
  model =
  c(
    "single variable model",
    "two variable model",
    "model with interaction term"
  ),
  adj_r_squared =
  c(
    glance(single_var_mod)$adj.r.squared, 
    glance(two_var_mod)$adj.r.squared, 
    glance(interaction_mod)$adj.r.squared
  )
)

model_comparison %>% 
  arrange(desc(adj_r_squared))

##                         model adj_r_squared
## 1 model with interaction term     0.6923468
## 2          two variable model     0.6487565
## 3       single variable model     0.6065174

The highest adjusted R² is given by the model with the interaction term followed by the two variable model, indicating that allowing both slopes and intercepts to differ leads to the better model.

Multiple numeric explanatory variables

I can add as many variables to linear models as we want, but the ability to visualize the model breaks down after 3 variables. I will extend my two variable model without interaction to include the square footage of the garage to compare the adjusted R² to those computed previously.

aug_three_var_mod <- lm(
  SalePrice ~ square_footage + GarageArea + bathroom_count, 
  data=housing_data
)

glance(aug_three_var_mod)$adj.r.squared

## [1] 0.6851754

The adjusted R² falls at 0.6851754, which is between the two variable model and the model with the interaction term that I put together.