Skip to main content

Reproducible Research Using R An Open Educational Resource

Christian Martinez

Section 8.9 Adding more variables

Throughout this chapter, we have only looked at 2 variables: price and volume to see if the former can predict the latter. But, what if we have more variables? Can we also utilize them in our regression model?
🔗
Let’s find out!
🔗
Let us create two new variables to add to our data: one where David Portnoy visited the pizza shop and another on how much was spent on advertising.
🔗
# Now, let's add some other variables to our data
set.seed(42)

# 1) Dave Portnoy visit.
pizza_sales <- pizza_sales |>
  mutate(portnoy_visit = rbinom(n(), size = 1, prob = 0.06))  # about 3-6% of weeks

# 2) Ad spend
pizza_sales <- pizza_sales |>
  mutate(ad_spend = round(runif(n(), 800, 6000), 0))
Now, since we have multiple variables, we can run what’s called multivariate regression model. To add more predictors, we simply use the + sign in our lm command after our first predictor.
🔗
m1 <- lm(volume ~ price, data = pizza_sales)                      # baseline

m2 <- lm(volume ~ price + portnoy_visit, data = pizza_sales)      # add portnoy visit

m3 <- lm(volume ~ price + ad_spend, data = pizza_sales)           # add ad spending

m4 <- lm(volume ~ price + ad_spend + portnoy_visit, data = pizza_sales)  # both
Boom! In the code above, we have now run 4 regression models, with every combination included.
🔗
Now, we can run the summary or glance command on each model individually, compare the models manually, but that would require a lot of working memory power on our end. Or we could utilize some commands in R to make our lives much easier.
🔗
# How do we know which model is the best?

AIC(m1, m2, m3, m4) |> arrange(AIC)
# This code provides us with an AIC (Akaike Information Criterion) value.
# The smaller the AIC, the better the model's trade-off between complexity and fit.

df      AIC
m2  4 3486.244
m4  5 3488.239
m1  3 3491.395
m3  4 3493.327

# Let's say we want to use the broom function again to do a model comparison
models <- list(m1, m2, m3, m4)
# In the code above, we are taking all of our models and putting them into a "list" of models

names(models) <- c("m1","m2","m3","m4")
# In the code above, we are just giving each item in the list a name

# In the code below, we are saying "Hey, run glance on every single item on the models list and create a data frame

map_df(models, ~glance(.x), .id = "model") |> select(model, r.squared,adj.r.squared,p.value,AIC)

# A tibble: 4 × 5
  model r.squared adj.r.squared  p.value   AIC
  <chr>     <dbl>         <dbl>    <dbl> <dbl>
1 m1        0.431         0.427 1.36e-20 3491.
2 m2        0.456         0.449 5.57e-21 3486.
3 m3        0.431         0.424 1.79e-19 3493.
4 m4        0.456         0.446 5.07e-20 3488.
With both of these, we now have:
🔗
  1. R²: A higher R squared means the model is a better predictor. Remember, R^2 always goes up with each added predictor, no matter how salient they actually are.
    🔗
  2. Adjusted R²: A higher adjusted R squared means the model is a better predictor. Adjusted R² always takes into account the number of predictors and does not inherently increase with more predictors.
    🔗
  3. p-value: If these models are statistically significant.
    🔗
  4. AIC: tells us how good our models are. A lower number means a better model.
    🔗
Overall, the goal is to create the most parsimonious regression model as possible. That is, the model with the fewest necessary predictors. As such, the model that is most parsimonious is m1, the model with just price. Yes, other models have higher adjusted R^2 values and lower AIC numbers, but only slightly. In this case, m1 is the model to go with.
🔗

Subsection 8.9.1 Bonus code

If you are looking to compare multiple regression models in a faster way, you could also utilize the code below.
🔗
full_model <- lm(volume ~ price + ad_spend + portnoy_visit, data = pizza_sales)

step_model <- step(full_model, direction = "both")

Start:  AIC=3043.53
volume ~ price + ad_spend + portnoy_visit

                Df  Sum of Sq        RSS    AIC
- ad_spend       1 1.4575e+06 4.4042e+10 3041.5
<none>                        4.4041e+10 3043.5
- portnoy_visit  1 2.0473e+09 4.6088e+10 3048.6
- price          1 3.1199e+10 7.5240e+10 3125.1

Step:  AIC=3041.54
volume ~ price + portnoy_visit

                Df  Sum of Sq        RSS    AIC
<none>                        4.4042e+10 3041.5
+ ad_spend       1 1.4575e+06 4.4041e+10 3043.5
- portnoy_visit  1 2.0659e+09 4.6108e+10 3046.7
- price          1 3.1252e+10 7.5294e+10 3123.2
🔗
🔗
PrevTopNext