library(readr)
## Warning: package 'readr' was built under R version 3.2.5
C3.Cars <- read_csv("P:/SUNY-Brockport/SUNY Brockport_1/teaching/Staitstical Methods/Spring 2017/R_Projects_in_class_spring17/Ch3/C3 Cars.csv")
## Parsed with column specification:
## cols(
##   Mileage = col_integer(),
##   Price = col_double(),
##   Make = col_character(),
##   Model = col_character(),
##   Trim = col_character(),
##   Type = col_character(),
##   Cyl = col_integer(),
##   Liter = col_double(),
##   Doors = col_integer(),
##   Cruise = col_integer(),
##   Sound = col_integer(),
##   Leather = col_integer()
## )
View(C3.Cars)

require(mosaic)
## Loading required package: mosaic
## Warning: package 'mosaic' was built under R version 3.2.5
## Loading required package: dplyr
## Warning: package 'dplyr' was built under R version 3.2.5
## 
## 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
## Loading required package: lattice
## Loading required package: ggplot2
## Warning: package 'ggplot2' was built under R version 3.2.5
## Loading required package: mosaicData
## Warning: package 'mosaicData' was built under R version 3.2.5
## Loading required package: Matrix
## 
## The 'mosaic' package masks several functions from core packages in order to add additional features.  
## The original behavior of these functions should not be affected by this.
## 
## Attaching package: 'mosaic'
## The following object is masked from 'package:Matrix':
## 
##     mean
## The following objects are masked from 'package:dplyr':
## 
##     count, do, tally
## The following objects are masked from 'package:stats':
## 
##     binom.test, cor, cov, D, fivenum, IQR, median, prop.test,
##     quantile, sd, t.test, var
## The following objects are masked from 'package:base':
## 
##     max, mean, min, prod, range, sample, sum
require(leaps)
## Loading required package: leaps
## Warning: package 'leaps' was built under R version 3.2.5
require(graphics)

##Packages lattice and dplyr are automatically added with mosaic
attach(C3.Cars)

#Create a scatterplot between Mileage and Price:
xyplot( Price ~ Mileage, data=C3.Cars, main="")

#fitting LM line or any arbitrary line using abline, 

#abline(h=4000)
plot(Mileage,Price)
abline(lm(Mileage~Price),col="red")
#fitting a lowess curve
lines(lowess(C3.Cars),col="green")
## Warning in data.matrix(x): NAs introduced by coercion
## Warning in data.matrix(x): NAs introduced by coercion

## Warning in data.matrix(x): NAs introduced by coercion

## Warning in data.matrix(x): NAs introduced by coercion

#The scatterplot( ) function in the car package 
#offers many enhanced features, including fit lines, marginal box plots,
#conditioning on a factor, and interactive point identification. 

# Enhanced Scatterplot of MPG vs. price
# by Number of Car Cylinders
library(car)
## Warning: package 'car' was built under R version 3.2.5
## 
## Attaching package: 'car'
## The following objects are masked from 'package:mosaic':
## 
##     deltaMethod, logit
## The following object is masked from 'package:dplyr':
## 
##     recode
scatterplot(Mileage~Price |Cyl , data=C3.Cars,
            xlab="Mileage", ylab="Price",
            main="Scatter Plot",
            labels=row.names(C3.Cars)) 

# for more on scatter plots refer to the below url
# http://www.statmethods.net/graphs/scatterplot.html

#subsetting data to look further into it
car.subset <- C3.Cars[C3.Cars$Price> 52000,]
head(car.subset)
## # A tibble: 6 × 12
##   Mileage    Price     Make  Model            Trim        Type   Cyl Liter
##     <int>    <dbl>    <chr>  <chr>           <chr>       <chr> <int> <dbl>
## 1     583 70755.47 Cadillac XLR-V8 Hardtop Conv 2D Convertible     8   4.6
## 2    6420 68566.19 Cadillac XLR-V8 Hardtop Conv 2D Convertible     8   4.6
## 3    7892 69133.73 Cadillac XLR-V8 Hardtop Conv 2D Convertible     8   4.6
## 4   12021 66374.31 Cadillac XLR-V8 Hardtop Conv 2D Convertible     8   4.6
## 5   15600 65281.48 Cadillac XLR-V8 Hardtop Conv 2D Convertible     8   4.6
## 6   18200 63913.12 Cadillac XLR-V8 Hardtop Conv 2D Convertible     8   4.6
## # ... with 4 more variables: Doors <int>, Cruise <int>, Sound <int>,
## #   Leather <int>
#Create a "regression" object, named cars.lm, that will contain 
#summary information about the fitted
#regression model of Price on Mileage:
#subset model
C3.Cars.subset <- lm(Mileage~Price,data=car.subset)
summary(C3.Cars.subset)
## 
## Call:
## lm(formula = Mileage ~ Price, data = car.subset)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2216.47 -1204.31   -72.41  1544.19  1779.85 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.473e+05  5.448e+03   27.04 3.77e-09 ***
## Price       -2.043e+00  8.617e-02  -23.71 1.07e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1638 on 8 degrees of freedom
## Multiple R-squared:  0.986,  Adjusted R-squared:  0.9842 
## F-statistic:   562 on 1 and 8 DF,  p-value: 1.067e-08
#full data model

cars.lm = lm(C3.Cars$Price~C3.Cars$Mileage,data=C3.Cars)

#The estimated coefficients, R^2, t-statistics, and other results 
#about the model can be extracted from cars.lm
#using the summary()command:
summary(cars.lm)
## 
## Call:
## lm(formula = C3.Cars$Price ~ C3.Cars$Mileage, data = C3.Cars)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -13905  -7254  -3520   5188  46091 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      2.476e+04  9.044e+02  27.383  < 2e-16 ***
## C3.Cars$Mileage -1.725e-01  4.215e-02  -4.093 4.68e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 9789 on 802 degrees of freedom
## Multiple R-squared:  0.02046,    Adjusted R-squared:  0.01924 
## F-statistic: 16.75 on 1 and 802 DF,  p-value: 4.685e-05
#ANOVA 
anova(cars.lm)  # anova table also is useful
## Analysis of Variance Table
## 
## Response: C3.Cars$Price
##                  Df     Sum Sq    Mean Sq F value    Pr(>F)    
## C3.Cars$Mileage   1 1.6056e+09 1605590375  16.755 4.685e-05 ***
## Residuals       802 7.6856e+10   95830165                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#Create a vector of residuals from the cars.lm model:
resid1 <- cars.lm$res

#residuals can be found using cars.lm$res or residuals(cars.lm)
#fitted/predicted values can be found using fitted(cars.lm)
fit<-fitted(cars.lm)
#Residual for the first car (Buick Century):
resid1[1]  # to show the first element of the vector of residuals
##         1 
## -6032.165
#It is often helpful to attach the residuals and fitted values to data frame.
Cars<-mutate(C3.Cars, resid1,fit)
Cars[1,c(1:3,13:14)] #to show the first row and selected columns of the Cars data frame
## # A tibble: 1 × 5
##   Mileage   Price  Make    resid1      fit
##     <int>   <dbl> <chr>     <dbl>    <dbl>
## 1    8221 17314.1 Buick -6032.165 23346.27
#############Stepwise regression########
#Note that the stepwise regression method implemented in R does not use 
#the improvement in R^2 as the basis for adding a new variable to 
#the model. 
#It uses a criterion called AIC. 
#Models with smaller AIC values are preferred. 
#To perform stepwise regression for the cars data in R, 
#first create a regression object that contains all relevant explanatory variables
################################################################
full.model<-lm(Price~Cyl+Liter+Doors+Cruise+Sound+Leather+Mileage, data = Cars)

#Now use the step function, and store the results in the object step.model:
step.model <- step(full.model)
## Start:  AIC=14331.39
## Price ~ Cyl + Liter + Doors + Cruise + Sound + Leather + Mileage
## 
##           Df  Sum of Sq        RSS   AIC
## - Liter    1   44992316 4.3492e+10 14330
## <none>                  4.3447e+10 14331
## - Sound    1  663674405 4.4111e+10 14342
## - Doors    1 1265022645 4.4712e+10 14352
## - Mileage  1 1548003060 4.4995e+10 14358
## - Cyl      1 1681894212 4.5129e+10 14360
## - Leather  1 1714076275 4.5161e+10 14360
## - Cruise   1 4986166639 4.8433e+10 14417
## 
## Step:  AIC=14330.22
## Price ~ Cyl + Doors + Cruise + Sound + Leather + Mileage
## 
##           Df  Sum of Sq        RSS   AIC
## <none>                  4.3492e+10 14330
## - Sound    1 6.8659e+08 4.4178e+10 14341
## - Doors    1 1.2297e+09 4.4722e+10 14351
## - Mileage  1 1.5632e+09 4.5055e+10 14357
## - Leather  1 1.6943e+09 4.5186e+10 14359
## - Cruise   1 4.9514e+09 4.8443e+10 14415
## - Cyl      1 1.3563e+10 5.7055e+10 14546
#To see the suggested final model, apply the summary command to step.model:
summary(step.model)
## 
## Call:
## lm(formula = Price ~ Cyl + Doors + Cruise + Sound + Leather + 
##     Mileage, data = Cars)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -13104  -5566  -1544   3877  33349 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  7.323e+03  1.771e+03   4.135 3.92e-05 ***
## Cyl          3.200e+03  2.030e+02  15.765  < 2e-16 ***
## Doors       -1.463e+03  3.083e+02  -4.747 2.45e-06 ***
## Cruise       6.206e+03  6.515e+02   9.525  < 2e-16 ***
## Sound       -2.024e+03  5.707e+02  -3.547 0.000412 ***
## Leather      3.327e+03  5.971e+02   5.572 3.45e-08 ***
## Mileage     -1.705e-01  3.186e-02  -5.352 1.14e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 7387 on 797 degrees of freedom
## Multiple R-squared:  0.4457, Adjusted R-squared:  0.4415 
## F-statistic: 106.8 on 6 and 797 DF,  p-value: < 2.2e-16
#Note that the final model left out one of the original seven variables.

######Best Subset Regression###################
#Requires loading the R package leaps. 
#The following code can be used to select the "best" models based on the lowest value of Cp. 
#You can find additional information about various options for the leaps function by using the help function
##################################################
require(leaps)
help(leaps)
## starting httpd help server ...
##  done
#First, create a "design matrix" whose columns contain the values of the explanatory variables of interest

full.data=cbind(Cars$Mileage,Cars$Cyl,Cars$Liter,Cars$Doors,Cars$Cruise,Cars$Sound,Cars$Leather)
head(full.data)
##       [,1] [,2] [,3] [,4] [,5] [,6] [,7]
## [1,]  8221    6  3.1    4    1    1    1
## [2,]  9135    6  3.1    4    1    1    0
## [3,] 13196    6  3.1    4    1    1    0
## [4,] 16342    6  3.1    4    1    0    0
## [5,] 19832    6  3.1    4    1    0    1
## [6,] 22236    6  3.1    4    1    1    0
#Now use the leaps function, and store the results in the object best.lm:

best.lm <- leaps(full.data,Cars$Price, method="Cp", names = c('Mileage','Cyl','Liter','Doors','Cruise','Sound','Leather'), nbest=3)
#method=c("Cp", "adjr2", "r2") # u can pick any of them
names(best.lm)  
## [1] "which" "label" "size"  "Cp"
data1<-cbind(best.lm$which,best.lm$Cp)
head(data1)
##   Mileage Cyl Liter Doors Cruise Sound Leather          
## 1       0   1     0     0      0     0       0 171.95875
## 1       0   0     1     0      0     0       0 189.68653
## 1       0   0     0     0      1     0       0 370.65991
## 2       0   1     0     0      1     0       0  87.57869
## 2       0   0     1     0      1     0       0 110.43981
## 2       0   1     0     1      0     0       0 145.78142
#####Cp statistic########## 
# is a measure of error in the best subset model, relative to incorporating all variables
#Adequate models have Cp roughly = parameters i nthe model( including the constant/intercept)
#and/or Cp is at its minimum
#####################################
#To identify the best models, examine the output of the output matrix we have called data1.
#Each row of the matrix provides a value of 1 if the explanatory variable 
#was included, and a value of 0 if it was not. For example, 
#the best model with one explanatory variable contains the predictor Cyl, and has an adjusted R-squared value of 0.323.


#We will use the regression model that contains all explanatory variables except Liter:
cars.lm2 <- lm(Price~Cyl+Doors+Cruise+Sound+Leather+Mileage)

#Create a vector of predicted (fitted) values of Price:
fit2 <- cars.lm2$fitted

#Create a vector of residuals:
resid2 <- cars.lm2$res

Cars<-mutate(Cars, resid2, fit2)

tail(Cars,10)  # shows the 10 rows with the largest residuals. 
## # A tibble: 10 × 16
##    Mileage    Price   Make    Model          Trim  Type   Cyl Liter Doors
##      <int>    <dbl>  <chr>    <chr>         <chr> <chr> <int> <dbl> <int>
## 1     5826 18173.98 Saturn L Series L300 Sedan 4D Sedan     6     3     4
## 2     7755 18490.98 Saturn L Series L300 Sedan 4D Sedan     6     3     4
## 3    10102 17322.08 Saturn L Series L300 Sedan 4D Sedan     6     3     4
## 4    10986 17978.36 Saturn L Series L300 Sedan 4D Sedan     6     3     4
## 5    14242 16425.17 Saturn L Series L300 Sedan 4D Sedan     6     3     4
## 6    16229 16507.07 Saturn L Series L300 Sedan 4D Sedan     6     3     4
## 7    19095 16175.96 Saturn L Series L300 Sedan 4D Sedan     6     3     4
## 8    20484 15731.13 Saturn L Series L300 Sedan 4D Sedan     6     3     4
## 9    25979 15118.89 Saturn L Series L300 Sedan 4D Sedan     6     3     4
## 10   35662 13585.64 Saturn L Series L300 Sedan 4D Sedan     6     3     4
## # ... with 7 more variables: Cruise <int>, Sound <int>, Leather <int>,
## #   resid1 <dbl>, fit <dbl>, resid2 <dbl>, fit2 <dbl>
# which Make and Type correspond to these values?

####################################
# model diagnostics using library car
#####################################

library(car)

# Assessing Outliers
outlierTest(cars.lm2) # Bonferonni p-value for most extreme obs
##     rstudent unadjusted p-value Bonferonni p
## 151 4.612092         4.6435e-06    0.0037334
## 153 4.549033         6.2291e-06    0.0050082
## 152 4.434826         1.0510e-05    0.0084498
## 154 4.255602         2.3331e-05    0.0187580
## 155 4.186102         3.1543e-05    0.0253610
## 156 4.055701         5.4902e-05    0.0441410
# Normality of Residuals
# qq plot for studentized resid
qqPlot(cars.lm2, main="QQ Plot")

# distribution of studentized residuals

library(MASS)
## 
## Attaching package: 'MASS'
## The following object is masked from 'package:dplyr':
## 
##     select
sresid <- studres(cars.lm2)
hist(sresid, freq=FALSE,
     main="Distribution of Studentized Residuals")
xfit<-seq(min(sresid),max(sresid),length=40)
yfit<-dnorm(xfit)
lines(xfit, yfit) 

# Evaluate homoscedasticity
# non-constant error variance test
ncvTest(cars.lm2)
## Non-constant Variance Score Test 
## Variance formula: ~ fitted.values 
## Chisquare = 139.5819    Df = 1     p = 3.285834e-32
# plot studentized residuals vs. fitted values
spreadLevelPlot(cars.lm2)

## 
## Suggested power transformation:  -0.1152622
#pearsons residuals vs fitted value plus curvature test 
residualPlot(cars.lm2)

##############################
# Evaluate MultiCollinearity
###############################
#initial diagnostics 
#Plot Cyl versus Liter:
par(mfrow=c(1,1))
plot(C3.Cars$Liter,C3.Cars$Cyl)

#Correlation between Cyl and Liter:
cor(C3.Cars$Liter,C3.Cars$Cyl)
## [1] 0.9578966
vif(cars.lm2) # variance inflation factors
##      Cyl    Doors   Cruise    Sound  Leather  Mileage 
## 1.167272 1.010767 1.164616 1.045804 1.049990 1.003421
sqrt(vif(cars.lm2)) > 2 # problem?
##     Cyl   Doors  Cruise   Sound Leather Mileage 
##   FALSE   FALSE   FALSE   FALSE   FALSE   FALSE
#############################
# Evaluate Nonlinearity
# component + residual plot
##############################
crPlots(cars.lm2)

###############################
# Test for Autocorrelated Errors
#################################
durbinWatsonTest(cars.lm2)
##  lag Autocorrelation D-W Statistic p-value
##    1       0.8462678     0.3042106       0
##  Alternative hypothesis: rho != 0
#another way of getting diagnostic curves
par(mfrow=c(2,2))
plot(cars.lm2)