Course 7 project
Feng Li
Executive Summary
This project assumes that you work for Motor Trend, a magazine about the automobile industry. Looking at a data set (mtcars from datasets library) of a collection of cars, they are interested in exploring the relationship between a set of variables and miles per gallon (MPG) (outcome). They are particularly interested in the following two questions:
- “Is an automatic or manual transmission better for MPG”
- “Quantify the MPG difference between automatic and manual transmissions”
The purpose of this project report is to address above two questions use the knowledge from course 7 regression models.
Load data and have a quick look
First load the data.
library(datasets)
data(mtcars)Question 1: “Is an automatic or manual transmission better for MPG”
Recall exploratory analysis approaches, we’ll do a data summary and two plots - boxplot and scatter plot to miles per gallon(mpg) versus transmission type(am).
library(ggplot2)
table(mtcars$am)##
## 0 1
## 19 13summary(mtcars$mpg[mtcars$am ==0])## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 10.40 14.95 17.30 17.15 19.20 24.40summary(mtcars$mpg[mtcars$am ==1])## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 15.00 21.00 22.80 24.39 30.40 33.90mpg.a <- mtcars$mpg[mtcars$am == 0]
mpg.m <- mtcars$mpg[mtcars$am == 1]
mtcars$am <- factor(mtcars$am,labels=c("Automatic","Manual"))
g <- ggplot(mtcars, aes(x=am, y=mpg, fill=am))
g1 <- g + geom_boxplot() +
geom_hline(aes(yintercept=mean(mpg.a)), linetype="dashed") +
geom_hline(aes(yintercept=mean(mpg.m)), linetype="dashed") +
ylab("Miles Per Gallon") +
xlab("Transmission Type")
g1
From above analysis, we know in this dataset, we have 19 automatic transmission cars and 13 manual transmission cars. the mpg mean value of automatic transmisison cars is 17.15mile/gallon which is less than manual transmission cars 24.39mile/gallon.
So to the first question, overall automatic transmission is better for MPG. From boxplot, we can clearly see this comparison.
Question 2: “Quantify the MPG difference between automatic and manual transmissions”
1. Statistical inference with a t test
Do a t test to mpg for automatic cars and mpg for manual cars. From this, we’ll know if there is a significant change to mpg when transmission changes from automatic to manual. So the H0 is mu =0 which means no change. We’ll check if we can reject.
t.test(mpg.a, mpg.m, paired=F, var.equal = F)##
## Welch Two Sample t-test
##
## data: mpg.a and mpg.m
## 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 of x mean of y
## 17.14737 24.39231Ok, so from p-value (0.001374) is far less than 0.05, we can tell that we are able to reject null hypothesis which means mpg has significant difference between automatic and manual transmissions.
2. Regression models analysis.
2.1 First, single variable model. We’ll apply linear models to mpg as outcome and am as predictor like following
fit.one <- lm(mpg ~ am, data=mtcars)
summary(fit.one)$coefficients## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 17.147368 1.124603 15.247492 1.133983e-15
## amManual 7.244939 1.764422 4.106127 2.850207e-04So we only see amManual varialbel listed in coefficient list. This is because am is categorical variable. In this case, R automatically takes one level of am as reference line or base line for coefficient analysis. Here amAutomatic is the baseline. the coefficient and pvalue of amManual is a comparison to amAutomatic.
This is exactly we want to answer the second question about the mpg difference between amAutomatic and amManual.
From output, we can tell following information:
1. The intercept 17.147 is the mean of of mpg from automatic transmission. This is reference line.
2. The coefficient of amManual is in a comparison to amAutomatic of the means from the two group. The 7.245 is the change of the mpg mean between amAutomatic and amManual (amManual - amAutomatic). If you want a mean of amManual, it’ll be 17.147 + 7.245 = 24.392.
3. p-value of amManual 0.000285 is to test whether or not amManual is different from amAutomatic. We can tell a significant difference as it’s much less than 0.05.
As we know from model selection, if we omitted the variables that should have been included in the model, out model ends up bias. So we’ll try to add more variables to the single variable model. But including unnecessary variables will cause variance inflation. We’ll need to find a best way through statistical analysis.
2.2 Second, consider more variales for better models.
As we know, a third variable can distort, or confound the relationship amonge two existing ones. What we want is to add those who would contribute instead those are unnecessary or distracted. This is a dynamic process and would cost lots of effort to do residual plot and diagnostics analysis to choose the best models.
In this project report, we’ll do analysis to two models as following.
- First, we’ll use nested likelihood ratio test to do multiple models with additional predictors. We try to find a best model . There are multiple ways to do this. The following is one of my tests:
fit.two <- lm(mpg ~ am + hp, data=mtcars)
fit.three <- lm(mpg ~ am + hp + wt, data=mtcars)
fit.four <- lm(mpg ~ am + hp + wt + qsec, data=mtcars)
fit.five <- lm(mpg ~ am + hp + wt + qsec + cyl, data=mtcars)
anova(fit.one, fit.two, fit.three, fit.four, fit.five)## Analysis of Variance Table
##
## Model 1: mpg ~ am
## Model 2: mpg ~ am + hp
## Model 3: mpg ~ am + hp + wt
## Model 4: mpg ~ am + hp + wt + qsec
## Model 5: mpg ~ am + hp + wt + qsec + cyl
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 30 720.90
## 2 29 245.44 1 475.46 77.3500 2.862e-09 ***
## 3 28 180.29 1 65.15 10.5987 0.003138 **
## 4 27 160.07 1 20.22 3.2903 0.081250 .
## 5 26 159.82 1 0.25 0.0405 0.842062
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1fit.chosen <- fit.fourWe can see from p-value that fit.two makes significant difference from fit.one and fit.three makes significant difference from fit.two. So we are getting better addming variables. fit.four has 0.081250 p-value compared to fit.three, it’s larger than 0.05 but really close. so we try one more fit.five. But we can see fit.five is not improving much anymore.
So we tend to use fit.four as our best model.
Now let’s do residuals plot and dfbetas and hatvalues checking to model fit.four.
summary(fit.chosen)##
## Call:
## lm(formula = mpg ~ am + hp + wt + qsec, data = mtcars)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.4975 -1.5902 -0.1122 1.1795 4.5404
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 17.44019 9.31887 1.871 0.07215 .
## amManual 2.92550 1.39715 2.094 0.04579 *
## hp -0.01765 0.01415 -1.247 0.22309
## wt -3.23810 0.88990 -3.639 0.00114 **
## qsec 0.81060 0.43887 1.847 0.07573 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.435 on 27 degrees of freedom
## Multiple R-squared: 0.8579, Adjusted R-squared: 0.8368
## F-statistic: 40.74 on 4 and 27 DF, p-value: 4.589e-11par(mfrow=c(2,2))
plot(fit.chosen)
par(mfrow=c(1,1))
# check out liers with high residuals
sort(round(dfbetas(fit.chosen)[,2], 3))## Volvo 142E Datsun 710 Pontiac Firebird
## -0.428 -0.304 -0.223
## Mazda RX4 Wag Hornet Sportabout Lotus Europa
## -0.199 -0.176 -0.147
## Mazda RX4 Merc 240D Cadillac Fleetwood
## -0.141 -0.118 -0.118
## Ford Pantera L Ferrari Dino Fiat X1-9
## -0.100 -0.057 -0.051
## Merc 450SL Merc 230 Merc 450SE
## -0.049 -0.047 -0.037
## Hornet 4 Drive Lincoln Continental Merc 280
## -0.034 -0.018 -0.016
## Camaro Z28 Porsche 914-2 Duster 360
## -0.014 0.010 0.013
## Honda Civic Valiant Maserati Bora
## 0.013 0.051 0.053
## Merc 450SLC Merc 280C Dodge Challenger
## 0.061 0.088 0.164
## AMC Javelin Toyota Corona Toyota Corolla
## 0.246 0.345 0.362
## Fiat 128 Chrysler Imperial
## 0.476 0.549# check leverage points
sort(round(hatvalues(fit.chosen),3))## Merc 450SL Merc 450SLC Merc 280C
## 0.060 0.061 0.063
## Merc 450SE Pontiac Firebird Merc 280
## 0.065 0.075 0.077
## Hornet 4 Drive Hornet Sportabout AMC Javelin
## 0.078 0.093 0.093
## Datsun 710 Valiant Fiat X1-9
## 0.096 0.098 0.106
## Ferrari Dino Volvo 142E Dodge Challenger
## 0.117 0.125 0.126
## Fiat 128 Honda Civic Merc 240D
## 0.128 0.128 0.130
## Porsche 914-2 Mazda RX4 Wag Mazda RX4
## 0.152 0.153 0.157
## Toyota Corolla Camaro Z28 Lotus Europa
## 0.160 0.162 0.164
## Duster 360 Toyota Corona Ford Pantera L
## 0.194 0.194 0.224
## Chrysler Imperial Cadillac Fleetwood Lincoln Continental
## 0.232 0.237 0.274
## Merc 230 Maserati Bora
## 0.460 0.520The following are the understanding to above plot and brief dignostic analysis:
1. From residuals-fitted value plot, we don’t significant pattern with only a bit curve. Which is good. 2. From QQ plot, it looks a normal distribution. It’s good too. 3. From scale-residual, all the points kind of randomly scattered. a bit linear pattern but not too bad. 4. From residual-leverage plot, we don’t see influential points which means not many outliers or leverage points. There maybe some points in upper right corner and lower right coner but they are far away from cook’s distance which means they are not influential to our model.
dfbetas and hatvalues outputs are also consistent to above plot4. From the outputs, we don’t see outliers with high residuals and we do see points that are suspicious to be leverage points.
- Last, we’ll do a model using all the variables and have a comparison to our fit.chosen.
fit.all <- lm(mpg ~ ., data=mtcars)
anova(fit.chosen, fit.all)## Analysis of Variance Table
##
## Model 1: mpg ~ am + hp + wt + qsec
## Model 2: mpg ~ cyl + disp + hp + drat + wt + qsec + vs + am + gear + carb
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 27 160.07
## 2 21 147.49 6 12.572 0.2983 0.9308So we see that our chosen model is not significant different from fit.call with all the variables. That means that we did not omitt necessary variables causing bias. At the same time, we did not included unnecessary variable causing model variance inflation.
In one word, we may find a “useful”" model so far. (As Brain said, all models are wrong, some models are useful…)