library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.2.1     ✔ readr     2.2.0
## ✔ forcats   1.0.1     ✔ stringr   1.6.0
## ✔ ggplot2   4.0.3     ✔ tibble    3.3.1
## ✔ lubridate 1.9.5     ✔ tidyr     1.3.2
## ✔ purrr     1.2.2     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(pROC)
## Type 'citation("pROC")' for a citation.
## 
## Attaching package: 'pROC'
## 
## The following objects are masked from 'package:stats':
## 
##     cov, smooth, var
epa_data <- read_csv("/Users/danielmedlin/Downloads/epa2012.csv")
## Rows: 1129 Columns: 28
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr  (19): mfr_name, division, carline, mfr_code, transmission_speed, guzzle...
## dbl   (8): model_yr, model_type_index, engine_displacement, no_cylinders, ci...
## date  (1): release_date
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

###Introduction

#Research Question How well does engine displacement, number of cylinders, and combined MPG predict whether a vehicle is classififed as a gas guzzler?

This assignemnt uses the epa2012 dataset from OpenIntro.org and contains information from the EPA abotu vehicles from the 2012 model year. The dataset has 1129 observations and 28 variables, with each observation representing a single vehicle. For this assignement, the main variables being used are the guzzler, engine_displacement, no_cylinders, and comb_mpg. The guzzler variable tells us whehter a vehicle is a guzzler or not using (Y) and (N). The engine_displacement variable describes the engine’s size. The no_cylinders variable tells us the number of cylinders, and comb_mpg tells us the combined gas mileage of each car.

I chose this topic because I was interested in how three other variables of a car determined whether the car was efficient with fuel. I will use logistical regression to examine how the three predictors influcen the probability that a vehicle will be classified as a gas guzzler.

Source: https://www.openintro.org/data/index.php?data=epa2012

###Data Analysis

epa_clean <- epa_data |> select(guzzler, engine_displacement, no_cylinders, comb_mpg)
epa_clean <- epa_clean |> drop_na() |> mutate(guzzler = as.factor(guzzler))
epa_clean |> group_by(guzzler) |> summarise(
  count = n(),
  avg_engine = mean(engine_displacement),
  avg_cylinders = mean(no_cylinders),
  avg_mpg = mean(comb_mpg))
## # A tibble: 2 × 5
##   guzzler count avg_engine avg_cylinders avg_mpg
##   <fct>   <int>      <dbl>         <dbl>   <dbl>
## 1 N        1038       3.32          5.63    22.1
## 2 Y          84       5.52          9.74    14.9

Before I did the logistical regression, I prepared epa2012 by selecting the important variables: guzzler, engine_displacement, no_cylinders, and comb_mpg. I removed observations containing missing values in these variables so th regression model would use only complete data. I also converter the guzzler variable to a factor because it will be the categorical outcome. I also examined the cleaned set and summarized the predictor variables.

###Statistical Analysis

I decided to use logistical regression for this assignment because the guzzler variable is a categorical variable, considering vehicles as either a gas guzzler or not. The variables engine_displacement, no_cylinders, and comb_mpg will be used as predictor variables to find what causes a vehicle year 2012 model vehicle to be considered a gas guzzler.

#Model

guzzler_model <- glm(
  guzzler ~ engine_displacement + no_cylinders + comb_mpg, 
  data = epa_clean,
  family = binomial
)
summary(guzzler_model)
## 
## Call:
## glm(formula = guzzler ~ engine_displacement + no_cylinders + 
##     comb_mpg, family = binomial, data = epa_clean)
## 
## Coefficients:
##                      Estimate Std. Error z value Pr(>|z|)    
## (Intercept)         -14.07973    3.67781  -3.828 0.000129 ***
## engine_displacement  -0.45774    0.26971  -1.697 0.089669 .  
## no_cylinders          2.13496    0.39067   5.465 4.63e-08 ***
## comb_mpg             -0.14241    0.08184  -1.740 0.081841 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 597.01  on 1121  degrees of freedom
## Residual deviance: 247.64  on 1118  degrees of freedom
## AIC: 255.64
## 
## Number of Fisher Scoring iterations: 9
exp(coef(guzzler_model))
##         (Intercept) engine_displacement        no_cylinders            comb_mpg 
##        7.678069e-07        6.327128e-01        8.456696e+00        8.672684e-01

The coefficients show that engine_displacement and comb_mpg are not statistically significant since their p-values are above 0.05, while no_cylinders is statistically significant since its less than 0.05. Calculating the log-odds, the no_clyinders variable shows that when holding for the other variables, each additional cylinder added increases log-odds by 2.135 for the vehicle to be considered a gas-guzzler.

#Diagnostics

epa_clean$prob <- predict(guzzler_model,
                          newdata = epa_clean,
                          type = "response")
epa_clean$predicted <- ifelse(epa_clean$prob > 0.5, "Y", "N")
conf_matrix <- table(
  Predicted = epa_clean$predicted,
  Actual = epa_clean$guzzler
)
conf_matrix
##          Actual
## Predicted    N    Y
##         N 1036   44
##         Y    2   40
TN <- conf_matrix["N", "N"]
FN <- conf_matrix["N", "Y"]
FP <- conf_matrix["Y", "N"]
TP <- conf_matrix["Y", "Y"]

accuracy <- (TP + TN)/ sum(conf_matrix)
sensitivity <- TP / (TP + FN)
specificity <- TN / (TN + FP)

accuracy
## [1] 0.9590018
sensitivity
## [1] 0.4761905
specificity
## [1] 0.9980732

The model shows an accuracy of 95.9%, a sensitivity of 47.6%, and a specificity of 99.8%. So while the model is highly accurate, the model can only identify around half of the gas guzzlers correctly. The model is almost perfect at finding cars that are not gas guzzlers since the specificity is 99.8%. So the model’s accuracy is overstated since there is an imbalace between cars that aren’t gas guzzlers to cars that are, and that the model is almost perfect at finding non-gas guzzlers, which dominate the dataset.

roc_curve <- roc(epa_clean$guzzler,
                 epa_clean$prob)
## Setting levels: control = N, case = Y
## Setting direction: controls < cases
plot(roc_curve,
     main ="ROC Curve")

auc(roc_curve)
## Area under the curve: 0.9542

The model has an AUC of 0.9543, which indicates the model is excellent at distinguishing between gas guzzler and non-gas-guzzlers. While the model has low sensitivity stated before, the AUC suggests the model is still strong overall at distinguishing.

###Conclusion

The logistic regression analysis examined whether engine displacement, number of cylinders, and combined MPG could predict whether a vehicle would be classified as a gas guzzler for 2012 models. The number of cylinders showed to be statistically significant since its p-value was less than 0.05. Each additional cylinder in a car increased the odds by 8.46 for a car to be classified a gas guzzler while holding for engine displacement and MPG. The other variables, engine displacement and combined MPG, were not statistically significant. The model had an overall accuracy of 95.9% and an AUC of 0.9542, meaning the model was strong overall at distinguishing gas guzzlers and non-gas guzzlers. However, the model’s sensitivity of 47.6% meant it was not good at finding gas-guzzlers. For future research, I would look to examine other characteristics like transmission type or vehicle class to see if other variables account for classifying gas-guzzlers.

###References OpenIntro. (n.d.). Vehicle info from the EPA for 2012 [epa2012 dataset]. OpenIntro.org. https://www.openintro.org/data/index.php?data=epa2012