Setup

library(tidyverse)

## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.1 ──

## ✓ ggplot2 3.3.5     ✓ purrr   0.3.4
## ✓ tibble  3.1.5     ✓ dplyr   1.0.7
## ✓ tidyr   1.1.4     ✓ stringr 1.4.0
## ✓ readr   2.0.2     ✓ forcats 0.5.1

## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()

Look

Use glimpse() to look at mpg.

Solution

glimpse(mpg)

## Rows: 234
## Columns: 11
## $ manufacturer <chr> "audi", "audi", "audi", "audi", "audi", "audi", "audi", "…
## $ model        <chr> "a4", "a4", "a4", "a4", "a4", "a4", "a4", "a4 quattro", "…
## $ displ        <dbl> 1.8, 1.8, 2.0, 2.0, 2.8, 2.8, 3.1, 1.8, 1.8, 2.0, 2.0, 2.…
## $ year         <int> 1999, 1999, 2008, 2008, 1999, 1999, 2008, 1999, 1999, 200…
## $ cyl          <int> 4, 4, 4, 4, 6, 6, 6, 4, 4, 4, 4, 6, 6, 6, 6, 6, 6, 8, 8, …
## $ trans        <chr> "auto(l5)", "manual(m5)", "manual(m6)", "auto(av)", "auto…
## $ drv          <chr> "f", "f", "f", "f", "f", "f", "f", "4", "4", "4", "4", "4…
## $ cty          <int> 18, 21, 20, 21, 16, 18, 18, 18, 16, 20, 19, 15, 17, 17, 1…
## $ hwy          <int> 29, 29, 31, 30, 26, 26, 27, 26, 25, 28, 27, 25, 25, 25, 2…
## $ fl           <chr> "p", "p", "p", "p", "p", "p", "p", "p", "p", "p", "p", "p…
## $ class        <chr> "compact", "compact", "compact", "compact", "compact", "c…

Tables

There are several categorical variables which might be useful in predicting highway gas mileage. Get tables of class, trans, fl and year.

Answer

table(mpg$class)

## 
##    2seater    compact    midsize    minivan     pickup subcompact        suv 
##          5         47         41         11         33         35         62

table(mpg$trans)

## 
##   auto(av)   auto(l3)   auto(l4)   auto(l5)   auto(l6)   auto(s4)   auto(s5) 
##          5          2         83         39          6          3          3 
##   auto(s6) manual(m5) manual(m6) 
##         16         58         19

table(mpg$fl)

## 
##   c   d   e   p   r 
##   1   5   8  52 168

table(mpg$year)

## 
## 1999 2008 
##  117  117

Objective

Build a model to predict hwy (highway mpg)

Our first model was an equation of the form

\[hwy = intercept + slope * displ\] ## Run Again

To refresh our minds, run the basic model. Then put the residuals from the model into the mpg dataframe and do a scatterplot of residuals against displ.

Solution

model1 = lm(hwy~displ,data = mpg)
mpg$residuals = model1$residuals
mpg %>% 
  ggplot(aes(x = displ,y=residuals)) +
  geom_point()

Other Possibilities

Non-Linearity

The graphical analysis of residuals shows that a quadratic model might be more appropriate than a linear model. Create the variable displ2, which is the square of displ. Add this variable to the dataframe. Then run a second model with the quadratic term. Compare the two on the basis of the standard error of the residuals (the RMSE).

Solution

mpg = mpg %>% 
  mutate(displ2 = displ^2)


model2 = lm(hwy~displ + displ2,data = mpg)

summary(model1)

## 
## Call:
## lm(formula = hwy ~ displ, data = mpg)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -7.1039 -2.1646 -0.2242  2.0589 15.0105 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  35.6977     0.7204   49.55   <2e-16 ***
## displ        -3.5306     0.1945  -18.15   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.836 on 232 degrees of freedom
## Multiple R-squared:  0.5868, Adjusted R-squared:  0.585 
## F-statistic: 329.5 on 1 and 232 DF,  p-value: < 2.2e-16

summary(model2)

## 
## Call:
## lm(formula = hwy ~ displ + displ2, data = mpg)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -6.6258 -2.1700 -0.7099  2.1768 13.1449 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  49.2450     1.8576  26.510  < 2e-16 ***
## displ       -11.7602     1.0729 -10.961  < 2e-16 ***
## displ2        1.0954     0.1409   7.773 2.51e-13 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.423 on 231 degrees of freedom
## Multiple R-squared:  0.6725, Adjusted R-squared:  0.6696 
## F-statistic: 237.1 on 2 and 231 DF,  p-value: < 2.2e-16

Adding the quadratic term reduces the RMSE from about 3.8 to about 3.4. This is good, but what does this model imply?

Predictions

Create predicted values of hwy for values of displ from .5 liters to 10 liters. Create predictions from these values using model2 and graph them.

Solution

displ = seq(from = .5, to = 10, by = .5)
displ2 = displ^2
new = data.frame(displ,displ2)
preds = predict(model2,new)
results = cbind(new,preds)
results %>% 
  ggplot(aes(x = displ,y = preds)) +
  geom_point()

Hmmmm? A 10 liter engine would get 40 mpg??

We got a better fit to the data we had, but with new data outside the range we had, the model fails. We should look at other ways to improve our model.

Categorical Variables

What are some other features that might influence hwy?

What about vehicle class, drive type, or fuel type?

Use geom_jitter to see how class effects hwy.

Solution

mpg %>% 
  ggplot(aes(x = class, y = hwy)) + geom_jitter()

What do you see?

What about fuel type? Repeat the exercise with fl.

Answer

mpg %>% 
  ggplot(aes(x = fl, y = hwy)) + geom_jitter()

What do you see?

Transmission

Repeat the exercise for trans.

Solution

mpg %>% 
  ggplot(aes(x = trans, y = hwy)) + geom_jitter()

Year

What about year?

Solution

mpg %>% 
  ggplot(aes(x = year, y = hwy)) + geom_jitter()

What do you see?

Drivetrain

What about drv?

Solution

mpg %>% 
  ggplot(aes(x = drv, y = hwy)) + geom_jitter()

The three values are significantly different.

Categorical Variables and Residuals

How could we see the relationship between a categorical variable and the residuals?

Use color in a scatterplot of residuals against displ in model1.

Do this for drv.

Solution

mpg %>% 
  ggplot(aes(x = displ,y=residuals, color = drv)) +
  geom_point()

Class

Do the same for class.

Solution

mpg %>% 
  ggplot(aes(x = displ,y=residuals, color = class)) +
  geom_point()