2.1) Features description:

The dataset basically contains several characteristics of cars, an insurance risk rating and normalized losses in use as compared to other cars.

Below we cond find more details about the different features:

1. normalized.losses: numerical - It is the relative average loss payment per insured vehicle year, this variables was normalized
2. make: ordinal with 22 levels - Car brand - “volkswagen”, “volvo”, “nissan”, “porsche”, “honda”, “subaru”, “mazda”, “jaguar”, “dodge”, “mercury”, “toyota”, “chevrolet”, “mercedes-benz”, “peugot”, “mitsubishi”, “plymouth”, “bmw”, “saab”, “isuzu”, “renault”, “alfa-romero” and “audi”
3. fuel.type: ordinal with 2 levels - Type of fuel - “diesel” and “gas”
4. aspiration: ordinal with 2 levels - Type of aspiration, it refers to breathin - “turbo” and “std”
5. num.of.doors: ordinal with 2 levels - Number of door - “four” and “two”
6. body.style: ordinal with 5 levels - Body style of the car (shape) - “hatchback” “sedan” “wagon” “hardtop” and “convertible”
7. drive.wheels: ordinal with 3 levels - The kind of drive wheel - “fwd” “rwd” “4wd”
8. engine.location: ordinal with 2 levels - The location of the engine - “front” “rear”
9. wheel.base: numerical - It is the distance between the centers of the front and rear wheels
10. length: numerical - The length of the car
11. width: numerical - The width of the car
12. height: numerical - The height of the car
13. curb.weight: numerical - The total mass of a vehicle with standard equipment and all necessary operating consumables
14. engine.type: ordinal with 7 levels - The kind of engine - “ohc” “ohcv” “l” “rotor” “ohcf” “dohc” “dohcv”
15. num.of.cylinders: ordinal with 7 levels - The number of cylinders - “four” “eight” “six” “three” “five” “twelve” “two”
16. engine.size: numerical - The size of the engine
17. fuel.system: ordinal with 9 levels - The kind of fuel system - “2bbl” “mfi” “mpfi” “1bbl” “idi” “spdi” “4bbl” “spfi”
18. bore: numerical - The diameter of each cylinder
19. stroke: numerical - The distance travelled by the piston in each cycle
20. compression.ratio: numerical - Measure based on the relative volumes of the combustion chamber and the cylinder
21. horsepower: numerical - The horse power of the engine
22. peak.rpm: numerical - Power band of the engine
23. city.mpg: numerical - It is the related distance traveled by a vehicle and the amount of fuel consumed in city
24. highway.mpg: numerical - It is the related distance traveled by a vehicle and the amount of fuel consumed in highway
25. price: numerical - Price of the car
26. symboling: ordinal with 7 levels - Indicates how safe is the car - -3: Highly secure car, -2: Moderate safe car, -1: Low safe car, 0: Neutral, not safe and not risky, +1 Low risky car, +2 Moderate risky car and +3 Highly risky car

For the purpose of this research, first of all, the target variable symboling will be transformed to three levels:

• Secure, corresponding with symboling levels -3, -2 and -1
• Neutral, neither safe nor risky car, corresponding with symboling levels 0
• Risky, corresponding with symboling levels +3, -+2 and +1

After will be transformed to factor.

cars$symboling <-
  plyr::revalue(as.character(cars$symboling),
                c("-3" = "secure",
                  "-2" = "secure",
                  "-1" = "secure",
                  "0" = "neutral",
                  "1" = "risky",
                  "2" = "risky",
                  "3" = "risky")) %>%
  as.factor()

Below we cand find a barplot of the target variable symboling:

As we can see in the graph, the levels are imbalanced. There are much more cars with symboling level risky than secure, hence we should take further this into consideration, maybe applying some type of resampling.

2.2) Numeric variables:

There are 15 numeric variables.

We can find below a density plot of these variables:

Below we can find the main statistical moments for the numeric variables:

2.3) Categorical variables:

There are 11 categorical variables, one of them is the target variable.

We can find below a barplot of this variables:

Fist of all, we will create crossvalidation tables with the taget variable symboling, and the rest of the ordinal variables. We will be able to see in percentage how the different variables are distributed, along the target variable:

1. For body.style and drive.wheels:

• Most total common cars:
  • For body.style - sedan
  • For drive.wheels - fwd
• Most secure and risky cars by features:
  • For body.style - convertible is the most secure car
  • For body.style - hatchback is the most risky car
  • For drive.wheels - 4wd is the most secure car
  • For drive.wheels - fwd is the most risky car

2. For fuel.type, aspiration, num.of.doors and engine.location:

• Most total common cars:
  • For fuel.type - gas
  • For aspiration - std
  • For num.of.doors - four
  • For engine.location - front
• Most secure and risky cars by features:
  • For fuel.type - gas is more likely to be risky or secure
  • For fuel.type - diesel is the most neutral car
  • For aspiration - std is the most secure car
  • For aspiration - turbo is the most risky car
  • For num.of.doors - four is the most secure car
  • For num.of.doors - two is the most risky car
  • For engine.location - front is more likely to be risky or secure
  • For engine.location - rear is the most neutral

3. For engine.type:

• Most total common cars:
  • For engine.type - ohc
• Most secure and risky cars by features:
  • For engine.type - ohxf is the most secure car
  • For engine.type - dohc is the most risky car

4. For num.of.cylinders:

• Most total common cars:
  • For num.of.cylinders - four
• Most secure and risky cars by features:
  • For num.of.cylinders - six is the most secure car
  • For num.of.cylinders - three is the most risky car

5. For fuel.system:

• Most total common cars:
  • For fuel.system - mpfi
• Most secure and risky cars by features:
  • For fuel.system - mpfi is the most secure car
  • For fuel.system - idi is the most risky car

6. For make:

• Most total common cars:
  • For make - toyota
• Most secure and risky cars by features:
  • For make - volvo is the most secure car
  • For make - saab is the most risky car

3.1) Var. transformation:

any(sapply(cars, is.character)) %>% 
  knitr::knit_print()

cars_preProces <-
  preProcess(cars, method = c("range"))
cars <- predict(cars_preProces,
               cars)

3.2) Missing values:

any(is.na(cars)) %>% 
  knitr::knit_print()

[1] FALSE

There is not any missing value, so we can continue further.

3.3) Unique variables:

Finally, we should check the numeric variables, in order to be sure tha they are unique.

To consider one variable as unique, a threshold of 500k was selected, hence half of the possible maximum number of unique values, 0.5 * 1.000.000 = 500.000

In order to be able to check this, the function find_if_unique_length was created:

find_if_unique_length <- function(x) {
  cars_numeric_vars_index <- c()
  for (i in cars_numeric_vars) {
    cars_numeric_vars_index <-
      c(cars_numeric_vars_index, grep(i, colnames(x)))
  }
  cars_numeric_vars_unique <- c()
  for (i in cars_numeric_vars_index) {
    a <- length(unique(x[, i]))
    if (a < (0.5) * dim(cars)[1]) {
      cars_numeric_vars_unique <-
        c(cars_numeric_vars_unique,
          paste(colnames(cars[i]), "NOT unique"))
    }
    else {
      cars_numeric_vars_unique <-
        c(cars_numeric_vars_unique, paste(colnames(cars[i]), "UNIQUE"))
    }
  }
  return(cars_numeric_vars_unique)
}

Now we will pass the function to our dataset:

data.frame(Results= find_if_unique_length(cars)) %>% 
 kable(align = "l",digits = 2) 

As per the selected threshold of 500k, we can conclude that all the variables are unique, so we are ready to divide our data in two parts.

4.1) For models that can manage ordinal variables:

We are ready to divide the data in two samples, train and test. The train sample is used to train the model, and the test sample is used to make the prediction and verify the performance of the model.

The data will be divided in 70% training, and 30% test, with the help of the function createDataPartition of the caret package:

set.seed(16)
cars_which_train <-
  createDataPartition(cars$symboling,
                      p = 0.7, 
                      list = FALSE)

cars_train <- cars[cars_which_train,]

cars_test <- cars[-cars_which_train,]

4.2) For the modesl sensible to data distribution:

There are some models sensible to data distribution based on euclidean distance.

In our dataset there are many factors and this can be time consuming for some models, hence all these factors will be trasnformed to numeric:

cars2 <- cars[-26]
indx <- sapply(cars2, is.factor)
cars2[indx] <- lapply(cars2[indx], function(x) as.numeric(as.character(x)))
cars2 <- cbind(cars2, cars[26])
cars_preProces1 <-
  preProcess(cars2, method = c("range"))
cars2 <- predict(cars_preProces1,
               cars2)

set.seed(16)
cars_which_train1 <-
  createDataPartition(cars2$symboling,
                      p = 0.7, 
                      list = FALSE)

cars_train1 <- cars2[cars_which_train1,]

cars_test1 <- cars2[-cars_which_train1,]

5.1) Correlation between features:

We will check the correlation between all the numerical variables, so we will use the list cars_numeric_vars.

Firstly, the correlation will be checked by the graph corplot:

It seems that the variables are not highly correlated, there are not values close to dark blue or dark red.

To be sure, we will check the maximum and minimum correlations, but in the maximum the value one should be removed, because each variable is evaluated against itself:

The maximum correlation was 0.17139, and the minimum -0.16221, hence there is not any variable to be omitted.

5.2) Relationship with the target variable:

Now we will check if the characteristic variables have relationship with the target variable symboling, so we will use the list cars_mult_bin_vars.

In order to be able to check this, we will use ANOVA, with the created function result_aov_pvalue:

result_aov_pvalue <- function(data, var) {
  result <- c()
  for (i in var) {
    if (summary(aov(data[, 10] ~ data[, i]))[[1]][["Pr(>F)"]][1] < 0.05) {
      result <-
        c(result,
          paste("Reject H0 -", i, "impact in symboling"))
    }
    else {
      result <-
        c(result,
          paste("NO reject H0 -", i, "has not impact in symboling"))
    }
  }
  return(result)
}

We will apply the function to our data:

The null hypothesis means that the variable has not impact in the taget variable (symboling).

For all the cases we reject this null hypothesis. Hence, considering 5% as the level of significance, we can conclude that all the characteristic variables have impact in the target variable.

5.3) Variables with near zero variance:

The variables with zero or near zero variance can have negative impact on the final result of the applied algorithm, for this reason is important to be checked.

The function nearZeroVar from the caret package will be used:

cars_nzv_stats <- nearZeroVar(cars_train,
                              saveMetrics = TRUE)
saveRDS(cars_nzv_stats, "cars_nzv_stats.rds")

cars_nzv_stats <- readRDS("cars_nzv_stats.rds")
cars_nzv_stats_res <- cars_nzv_stats %>%
  rownames_to_column("variable") %>%
  arrange(-zeroVar, -nzv, -freqRatio)
cars_nzv_stats_res[c(1, 4:5)] %>% 
kable(align = "l",digits = 2)

cars_nzv_unsel <- c(cars_nzv_stats_res[1][cars_nzv_stats_res[5] == TRUE])
sort(cars_nzv_unsel) %>% 
  knitr::knit_print()

character(0)

As we can see there is not any variable with TRUE for nzv and zeroVar, so the model does not suggest to omit any feature.

It can be concluded that no variables should be omitted considering near zero variance method.

5.4) Linear regressions in dataset:

We will check if in the dataset carsthere is linear regression with the function findLinearCombos:

cars_linearCombos <- findLinearCombos(cars_train[, cars_numeric_vars])  %>% 
  knitr::knit_print()

$linearCombos list()

$remove NULL

The function returns NULL, so there is not linear regression in our data.

5.5) Rank features - Learning Vector Quantization:

Learning Vector Quantization will be applied in order to rank the features by importance in relationship with symboling.

This algorithm is being applied with down-sampling and without cross-validation, and it can give an idea of which feauture could be omitted.

set.seed(16)
ctrl_cvnone <- trainControl(method = "none",
                            sampling = "down")
rank_features <-
  train(
    symboling ~ .,
    data = cars_train,
    method = "lvq",
    trControl = ctrl_cvnone
  )
saveRDS(rank_features, "rank_features.rds")

Learning Vector Quantization

700001 samples 25 predictor 3 classes: ‘neutral’, ‘risky’, ‘secure’

No pre-processing Resampling: None Addtional sampling using down-sampling

• Most important feautures by level:
  • neutral: engine.type
  • risky: num.of.doors
  • secure: engine.type
• Less important feautures by level:
  • neutral: compression.ratio
  • risky: compression.ratio
  • secure: compression.ratio

As we can see, the feature engine.type is one of the most important for all the levels.

On the other hand, in case we need to ommit some feautures, probably we will choose compression.ratio, as it is the least important for all the levels.

5.6) Conclusions:

The previous analysis did not suggest to ommit any feauture, but the variable make has too many levels. Cosequently, it can be a problem in terms of computational time when training the different models.

For this reason, in the less demanding models will be tested whether omitting the feature make gives much worst results.

6.1) Functions:

Two functions were created in order to evaluate the results of the models:

accuracy_multinom <- function(predicted, real) {
  ctable_m <- table(predicted,
                    real)
  accuracy <- (100 * sum(diag(ctable_m)) / sum(ctable_m))
  base_ <- diag(ctable_m) / colSums(ctable_m)
  balanced_accuracy <- mean(100 * ifelse(is.na(base_), 0, base_))
  base_2 <- diag(ctable_m) / rowSums(ctable_m)
  correctly_predicted <-
    mean(100 * ifelse(is.na(base_2), 0, base_2))
  return(
    data.frame(
      accuracy = accuracy,
      balanced_accuracy = balanced_accuracy,
      balanced_correctly_predicted = correctly_predicted
    )
  )
}

plot_model_fitted <- function(fitt) {
  require(scales)
  a <- data.frame(Symboling = fitt)
  ggplot(a, aes(Symboling)) +
    geom_bar(fill = "#9F1E42") +
    theme_minimal() +
    scale_y_continuous(labels = comma) 
    
}

6.2) MLR - Multinomial Logistic Regression:

The first model to be applied will be multinomial logistic regression. It is a classification method that generalize logistic regression to multiclass problems, in this case three classes or levels.

set.seed(16)
mlr_multinomial <- multinom(symboling ~ .,
                    data = cars_train,
                      maxit = 1000)
saveRDS(mlr_multinomial, "mlr_multinomial.rds")

set.seed(16)
mlr_multinomial_fitted <- predict(mlr_multinomial,
                              cars_test)
saveRDS(mlr_multinomial_fitted, "mlr_multinomial_fitted.rds")

set.seed(16)
mlr_multinomial_fitted_prob <- predict(mlr_multinomial,
                              cars_test,
                              type= "prob")
saveRDS(mlr_multinomial_fitted_prob, "mlr_multinomial_fitted_prob.rds")

set.seed(16)
mlr_multinomial_sel <- multinom(symboling ~ .,
                    data = cars_train %>% 
                      dplyr::select(-make),
                      maxit = 1000)
saveRDS(mlr_multinomial_sel, "mlr_multinomial_sel.rds")

mlr_multinomial_sel <- readRDS("mlr_multinomial_sel.rds")
data.frame(Residual.Deviance = round(mlr_multinomial_sel[["deviance"]], 2), AIC = 
             round(mlr_multinomial_sel[["AIC"]], 2)) %>%
  kable(align = "l",digits = 2) 

set.seed(16)
mlr_multinomial_sel_fitted <- predict(mlr_multinomial_sel,
                              cars_test)
saveRDS(mlr_multinomial_sel_fitted, "mlr_multinomial_sel_fitted.rds")

set.seed(16)
mlr_multinomial_sel_fitted_prob <- predict(mlr_multinomial_sel,
                              cars_test,
                              type= "prob")
saveRDS(mlr_multinomial_sel_fitted_prob, "mlr_multinomial_sel_fitted_prob.rds")

6.3) PMR - Penalized Multinomial Regression:

Caret package brings several models, one of them is Penalized Multinomial Regression. It is the closest one to the normal multinomial model.

The benefits of applying this algorithm with caret package are that we will be able to apply resampling and cross-validation. Adittionally, as the previous normal multinomial model, we will apply 1.000 as the maximum number of iterations.

For this model we will apply cross-validation with 5 fold:

set.seed(16)

ctrl_cv <- trainControl(method = "cv",
                          number = 5)
options(scipen=999)
pmr_multinomial <- train(
  symboling ~ .,
  data = cars_train,
  method = "multinom",
  trControl = ctrl_cv,
  maxit = 1000
)
saveRDS(pmr_multinomial, "pmr_multinomial.rds")

set.seed(16)
pmr_multinomial_fitted <- predict(pmr_multinomial,
                              cars_test)
saveRDS(pmr_multinomial_fitted, "pmr_multinomial_fitted.rds")

set.seed(16)
pmr_multinomial_fitted_prob <- predict(pmr_multinomial,
                              cars_test,
                              type= "prob")
saveRDS(pmr_multinomial_fitted_prob, "pmr_multinomial_fitted_prob.rds")

6.4) Test - Choosing the kind of sampling:

As we saw the multinom model of caret package is not really demanding in terms of computational time, so we will use this model to test which sampling method should be applied in KNN and SVM.

We will see if the accuracy changes if we apply no-sampling and down-sampling.

set.seed(16)

ctrl_cv <- trainControl(method = "cv",
                        number = 2)
options(scipen=999)
pmr_multinomial_nosam <- train(
  symboling ~ .,
  data = cars_train,
  method = "multinom",
  trControl = ctrl_cv,
  maxit = 1000
)
saveRDS(pmr_multinomial_nosam, "pmr_multinomial_nosam.rds")

set.seed(16)
pmr_multinomial_nosam_fitted <- predict(pmr_multinomial_nosam,
                              cars_test)
saveRDS(pmr_multinomial_nosam_fitted, "pmr_multinomial_nosam_fitted.rds")

set.seed(16)
pmr_multinomial_nosam_fitted_prob <- predict(pmr_multinomial_nosam,
                              cars_test,
                              type= "prob")
saveRDS(pmr_multinomial_nosam_fitted_prob, "pmr_multinomial_nosam_fitted_prob.rds")

set.seed(16)
ctrl_cv <- trainControl(method = "cv",
                          number = 2,
                          sampling = "down",
                        classProbs = TRUE,
                        summaryFunction = fiveStats)
options(scipen=999)
pmr_multinomial_down <- train(
  symboling ~ .,
  data = cars_train,
  method = "multinom",
  trControl = ctrl_cv,
  maxit = 1000
)
saveRDS(pmr_multinomial_down, "pmr_multinomial_down.rds")

set.seed(16)
pmr_multinomial_down_fitted <- predict(pmr_multinomial_down,
                              cars_test)
saveRDS(pmr_multinomial_down_fitted, "pmr_multinomial_down_fitted.rds")

set.seed(16)
pmr_multinomial_down_fitted_prob <- predict(pmr_multinomial_down,
                              cars_test,
                              type= "prob")
saveRDS(pmr_multinomial_down_fitted_prob, "pmr_multinomial_down_fitted_prob.rds")

6.5) KNN - K-Nearest Neighbors:

k-nearest neighbors algorithm is one of the most famous classification methods. In this case the input consists of the k closest training samples in the feature space.

It uses euclidean distance in order to measure the distance between neighbors, for this reason, for this algorithm is important to have scaled data, in our case the data was scaled before, with the range method, in the range [0, 1].

The algorithm will be applied with cross-validation, 3 folds, and down-sampling.

Adittionally, we will apply tuneGrid, in order to personalize more parameters, in this case several values of k are going to be tested, corresponding to the sequence seq(5, 89, 14))). Hence, 18 values will be tested 5, 19.. till 80. The algorithm will return the finalmodel with the greatest accuracy.

 set.seed(16) 
 k_value <- data.frame(k = c(seq(5, 145, 28)))
 ctrl_cv3 <- trainControl(method= "cv",
                             number= 3,
                          sampling= "down") 
 knn_model <- train(symboling ~ ., 
                      cars_train1 %>% 
                      dplyr::select(-make), 
                      method = "knn", 
                      trControl =  ctrl_cv3,
                      tuneGrid = k_value) 
 saveRDS(knn_model, "knn_model.rds") 

 set.seed(16) 
 knn_model_fitted <- predict(knn_model, 
                               cars_test1) 
 saveRDS(knn_model_fitted, "knn_model_fitted.rds") 

 set.seed(16) 
 knn_model_fitted_prob <- predict(knn_model, 
                               cars_test1, 
                               type= "prob") s
 saveRDS(knn_model_fitted_prob, "knn_model_fitted_prob.rds") 

 knn_model$finalModel$k %>% 
  knitr::knit_print()

6.6) LDA - Linear Discriminant Analysis:

Now Linear discriminant analysis will be applied. This algorithm is popular when dealing with a multiclass dependent variable, like our variable symboling.

It finds a linear combination of features that characterizes or separates two or more classes of objects.

In this case it will be applied with repeat cross-validation, 10 folds, repeated 3 times and down sampling In this case as it is not based on euclidean distance, it can be applied to train_cars instead of train_cars1. Aditionally, the feature make will not be omitted, because this model can manage very well the ordinal variables with several levels.

 set.seed(16) 
 ctrl_cv10 <- trainControl(method= "repeatedcv", 
                           number = 10, 
                           repeats = 3,
                           sampling= "down") 
lda_model <- train( 
   symboling ~ ., 
   data = cars_train, 
   method = "lda", 
   trControl = ctrl_cv10) 
 saveRDS(lda_model , "lda_model.rds") 

 set.seed(16) 
 lda_model_fitted <- predict(lda_model, 
                               cars_test) 
 saveRDS(lda_model_fitted, "lda_model_fitted.rds") 

 set.seed(16) 
 lda_model_fitted_prob <- predict(lda_model, 
                               cars_test, 
                               type= "prob") 
 saveRDS(lda_model_fitted_prob, "lda_model_fitted_prob.rds") 

6.7) QDA - Quadratic Discriminant Analysis:

Similarly to the case of SVM there are other similar algorithms that apply other kind of combination rather than linear, for example quadratic, this last will be applied in this case.

In this case it will be applied with repeat cross-validation, 10 folds, repeated 3 times and down sampling, in this case as it is not based on euclidean distance, it can be applied to train_cars instead of train_cars1. Aditionally, the feature make will not be omitted, because this model can manage well the ordinal variables with several levels.

 set.seed(16) 
 ctrl_cv10 <- trainControl(method= "repeatedcv", 
                           number = 10, 
                           repeats = 3,
                           sampling= "down") 
qda_model <- train( 
   symboling ~ ., 
   data = cars_train, 
   method = "qda", 
   trControl = ctrl_cv10) 
 saveRDS(qda_model , "qda_model.rds") 

 set.seed(16) 
 qda_model_fitted <- predict(qda_model, 
                               cars_test) 
 saveRDS(qda_model_fitted, "qda_model_fitted.rds") 

 set.seed(16) 
 qda_model_fitted_prob <- predict(qda_model, 
                               cars_test, 
                               type= "prob") 
 saveRDS(qda_model_fitted_prob, "qda_model_fitted_prob.rds") 

7.1) Summary:

7.2) Conclusions:

• For the dataset cars the best classification model is Quadratic Discriminant with repeated cross validation and down-sampling
• In imbalanced data one cannot trust the accuracy to compare models, the model tends to predict the most common class, hence other kind of measures are recomended, in our case ROC was used
• K nearest neighbors and support vector machine are models that do not work very well with multiclass dependent variable, it would be better to use them in case of binomial depdent varaible
• Before applying demanding models, it is convenient to analize ceteris paribus what is expected to happen, in our case it was analized what will happen when down-sampling is applied
• Computational time matters, and it should be taken into consideration when applying machine learning models, hence the data preparation is really important in these cases. In certain instances it is recommended to parallel the process, when the algorithm allows to do it (doParallel and doMC packages can be used)
• In case of having a big data and imbalance clases it is really productive to use down-sampling, in terms of computational time and performance of the model
• When performing a Machine Learning analysis it is important to have enough memory in the system which would allow one to save the results