Executive Summary

The data was extracted from the 1974 Motor Trend US magazine, and comprises fuel consumption and 10 aspects of automobile design and performance for 32 automobiles (1973–74 models). The goal here is to answer the question; Is an automatic or manual transmission better for MPG? The other goal is to quantify the MPG difference between automatic and manual transmissions

Very quickly it became apparent that manual transmission on average get better MPG. This was shown easily with a priliminary t test that gave a predicted true difference of 7.244. Through modeling can we get a more accurate number?

The best model I came up with is {r}lm(mpg ~ am + hp, mtcars). This model provides that the true difference in means 5.277. This takes into account variance inflation, variable correlation, analysis of variance, and leverage points. I walk through my process below.

Data Manipulation

data(mtcars)

library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(ggplot2)
library(plyr)
## -------------------------------------------------------------------------
## You have loaded plyr after dplyr - this is likely to cause problems.
## If you need functions from both plyr and dplyr, please load plyr first, then dplyr:
## library(plyr); library(dplyr)
## -------------------------------------------------------------------------
## 
## Attaching package: 'plyr'
## The following objects are masked from 'package:dplyr':
## 
##     arrange, count, desc, failwith, id, mutate, rename, summarise,
##     summarize
library(car)

str(mtcars)
## 'data.frame':    32 obs. of  11 variables:
##  $ mpg : num  21 21 22.8 21.4 18.7 18.1 14.3 24.4 22.8 19.2 ...
##  $ cyl : num  6 6 4 6 8 6 8 4 4 6 ...
##  $ disp: num  160 160 108 258 360 ...
##  $ hp  : num  110 110 93 110 175 105 245 62 95 123 ...
##  $ drat: num  3.9 3.9 3.85 3.08 3.15 2.76 3.21 3.69 3.92 3.92 ...
##  $ wt  : num  2.62 2.88 2.32 3.21 3.44 ...
##  $ qsec: num  16.5 17 18.6 19.4 17 ...
##  $ vs  : num  0 0 1 1 0 1 0 1 1 1 ...
##  $ am  : num  1 1 1 0 0 0 0 0 0 0 ...
##  $ gear: num  4 4 4 3 3 3 3 4 4 4 ...
##  $ carb: num  4 4 1 1 2 1 4 2 2 4 ...
summary(mtcars)
##       mpg             cyl             disp             hp       
##  Min.   :10.40   Min.   :4.000   Min.   : 71.1   Min.   : 52.0  
##  1st Qu.:15.43   1st Qu.:4.000   1st Qu.:120.8   1st Qu.: 96.5  
##  Median :19.20   Median :6.000   Median :196.3   Median :123.0  
##  Mean   :20.09   Mean   :6.188   Mean   :230.7   Mean   :146.7  
##  3rd Qu.:22.80   3rd Qu.:8.000   3rd Qu.:326.0   3rd Qu.:180.0  
##  Max.   :33.90   Max.   :8.000   Max.   :472.0   Max.   :335.0  
##       drat             wt             qsec             vs        
##  Min.   :2.760   Min.   :1.513   Min.   :14.50   Min.   :0.0000  
##  1st Qu.:3.080   1st Qu.:2.581   1st Qu.:16.89   1st Qu.:0.0000  
##  Median :3.695   Median :3.325   Median :17.71   Median :0.0000  
##  Mean   :3.597   Mean   :3.217   Mean   :17.85   Mean   :0.4375  
##  3rd Qu.:3.920   3rd Qu.:3.610   3rd Qu.:18.90   3rd Qu.:1.0000  
##  Max.   :4.930   Max.   :5.424   Max.   :22.90   Max.   :1.0000  
##        am              gear            carb      
##  Min.   :0.0000   Min.   :3.000   Min.   :1.000  
##  1st Qu.:0.0000   1st Qu.:3.000   1st Qu.:2.000  
##  Median :0.0000   Median :4.000   Median :2.000  
##  Mean   :0.4062   Mean   :3.688   Mean   :2.812  
##  3rd Qu.:1.0000   3rd Qu.:4.000   3rd Qu.:4.000  
##  Max.   :1.0000   Max.   :5.000   Max.   :8.000
plot(mtcars)

mtcars$cyl <- factor(mtcars$cyl)
mtcars$vs <- factor(mtcars$vs)
mtcars$gear <- factor(mtcars$gear)
mtcars$carb <- factor(mtcars$carb)
mtcars$am <- factor(mtcars$am)
mtcars$CarModel <- rownames(mtcars)

Initial Exploration

Here I do the initial t test and boxplot. There apprears to be a significant difference however other confounding variables are not taken into account.

t.test(mpg ~ am, mtcars)
## 
##  Welch Two Sample t-test
## 
## data:  mpg by am
## t = -3.7671, df = 18.332, p-value = 0.001374
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -11.280194  -3.209684
## sample estimates:
## mean in group 0 mean in group 1 
##        17.14737        24.39231
ggplot(mtcars, aes(x = am, y = mpg, colour = factor(am),label = CarModel))+
        geom_boxplot() +
        geom_point() +
        geom_text(data  = subset(mtcars, mpg < 12 | mpg > 32), vjust = "inward", nudge_x = 0.01, hjust ="inward")

Further Exploration

The goal here is to look for the most confounding variables so that they can be accunted for in the model. Inital observation was very interesting. Many variables seemed to impact mpg directly. The two variables that seem the most influential are hp and wt. These will be my focus as I continue into the modeling process.

notes on specific variables

  • Cyl: seems to have effect but is a factor variable and the smaller less cyl cars tend to be standard. Definate contender for t=further investigation. Also higly related to hp.( use either or)
  • V shape: engine seem to have a siginificant effect as well. Another contender for investigation
  • qsec: is a very interesting figure. From it there appears to be very little relation with am and qsec. (will not be considered in further investigation)
  • drat: very variable relationship, not a good regressor
  • gear: seems to be relationship but small amount of data for each factor, going to cause high variable inflation

Best regressors: Cyl, V, hp, and wt

ggplot(mtcars, aes(x = hp, y = mpg, colour = factor(am), label = CarModel)) + 
        geom_point() + 
        geom_text(data  = subset(mtcars, hp > 225 | mpg > 25), vjust = "inward", nudge_x = 0.01, hjust ="inward")

ggplot(mtcars, aes(x = wt, y = mpg, colour = factor(am), label = CarModel)) +
        geom_point() + 
        geom_text(data  = subset(mtcars, wt > 4 | mpg > 25), vjust = "inward", nudge_x = 0.01, hjust ="inward")

ggplot(mtcars, aes(x = cyl, y = mpg, colour = factor(am), label = CarModel)) +
        geom_point() +
        geom_text(data  = subset(mtcars, wt > 4 | mpg > 25), vjust = "inward", nudge_x = 0.01, hjust ="inward")

ggplot(mtcars, aes(x = vs, y = mpg, colour = factor(am), label = CarModel)) + geom_point() +
        geom_text(data  = subset(mtcars, wt > 4 | mpg > 25), vjust = "inward", nudge_x = 0.01, hjust ="inward")

ggplot(mtcars, aes(x = qsec, y = mpg, colour = factor(am))) + geom_point() + geom_smooth(method = lm)

ggplot(mtcars, aes(x = drat, y = mpg, colour = factor(am), label = CarModel)) + geom_point() + geom_text(data  = subset(mtcars, wt > 4 | mpg > 25), vjust = "inward", nudge_x = 0.01, hjust ="inward")

ggplot(mtcars, aes(x = gear, y = mpg, colour = factor(am))) + geom_point()

Constucing Multivariable Linear Model

Note that construcing and choosing a linear model was not a linear process I went back and forth looking at the impact of changes in the uncertainty in the model. Ultimatly the goal was to reduce uncertainty in this section.

Observing Uncertainty

Going into this process I knw that hp and wt where suspect. This was confirmed when fitting different model variable and checking for significance. Many were signifiacnt but the most was hp. Then adding another variable wt was the next most significant. This matches my intial investigations!

Two other variables that were investigated were vs and cyl. Both were significant but less so then hp or wt. Cyl remained signifant and could be included in the model but under further investigation causes high inflation.

fit <- lm(mpg ~ am, mtcars)
fit2 <- lm(mpg ~ am + hp, mtcars)
fit3 <- lm(mpg ~ am + wt + hp, mtcars)
fit4 <- lm(mpg ~ am + wt + hp + cyl, mtcars)
fit5 <- lm(mpg ~ am + wt + cyl + hp + vs, mtcars)

anova(fit, fit2, fit3, fit4, fit5)
## Analysis of Variance Table
## 
## Model 1: mpg ~ am
## Model 2: mpg ~ am + hp
## Model 3: mpg ~ am + wt + hp
## Model 4: mpg ~ am + wt + hp + cyl
## Model 5: mpg ~ am + wt + cyl + hp + vs
##   Res.Df    RSS Df Sum of Sq       F    Pr(>F)    
## 1     30 720.90                                   
## 2     29 245.44  1    475.46 82.7287 2.096e-09 ***
## 3     28 180.29  1     65.15 11.3357  0.002462 ** 
## 4     26 151.03  2     29.27  2.5461  0.098534 .  
## 5     25 143.68  1      7.35  1.2782  0.268968    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual Plotting and Diagnostics

This is was done simutaneosly with constructing the linear model. Anova testing narrowed the model choice down to two models, lm(mpg ~ am + hp, mtcars) and lm(mpg ~ am + wt + hp, mtcars). Now residual plotting. Ignore the warning message it is refering to factors not being sensible for vifs but the variables we are observing are not factors ( i.e hp and wt).

Through further investigation, the additional wt factor has a big impact on the variable inflation. This is why I choose lm(mpg ~ am + wt + hp, mtcars) as the best model. This also makes sense. Hp back in the 70’s was very colsely related to the size of the engine which contributes most to the weight of the car. That is why it is not sensible to inclue both wt and hp in the model because that would be observing the same varialbe twice. Therefore hp can be a good predictor of car weight.

vif(fit2)
##       am       hp 
## 1.062867 1.062867
vif(fit3)
##       am       wt       hp 
## 2.271082 3.774838 2.088124

Looking at Leverage Points and Residual Plotting

Observing leverage points throughout the analysis is interesing. In the exploration I tried to keep track of the outliers by labeling them in the plots. These where then observed again through residual ploting and beta testing. The Maserati Bora is shown to have lots of leverage and is a possible removal. Depening on the study the Bora is a completly different class of car than the others.

betas <- dfbetas(fit3)
betas <- as.data.frame(betas)
betas$CarModel <- rownames(betas)

plot(fit2)

plot(fit3)

am1.betas <- arrange(betas, desc(abs(am1)))[1:5,c(2,5)]
wt.betas <- arrange(betas, desc(abs(wt)))[1:5,c(3,5)]
hp.betas <- arrange(betas, desc(abs(hp)))[1:5,c(4,5)]

hp.betas
##           hp          CarModel
## 1  0.5809822     Maserati Bora
## 2 -0.4329640          Fiat 128
## 3 -0.3350091 Chrysler Imperial
## 4 -0.2400377    Toyota Corolla
## 5  0.2364049      Lotus Europa
wt.betas
##           wt          CarModel
## 1  0.9726881 Chrysler Imperial
## 2 -0.3897100      Lotus Europa
## 3 -0.2552249     Mazda RX4 Wag
## 4  0.2481860     Toyota Corona
## 5  0.2100033          Fiat 128