Load Libraries

library(caTools)
library(xgboost)
library(caret)
## Loading required package: ggplot2
## Loading required package: lattice
library(Metrics)
## 
## Attaching package: 'Metrics'
## The following objects are masked from 'package:caret':
## 
##     precision, recall
library(readr)
library(fastDummies)
## Thank you for using fastDummies!
## To acknowledge our work, please cite the package:
## Kaplan, J. & Schlegel, B. (2023). fastDummies: Fast Creation of Dummy (Binary) Columns and Rows from Categorical Variables. Version 1.7.1. URL: https://github.com/jacobkap/fastDummies, https://jacobkap.github.io/fastDummies/.
library(stringr)
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following object is masked from 'package:xgboost':
## 
##     slice
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(xgboost) # for xgboost
library(tidyverse) 
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ forcats   1.0.0     ✔ tibble    3.2.1
## ✔ lubridate 1.9.3     ✔ tidyr     1.3.0
## ✔ purrr     1.0.2
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ✖ purrr::lift()   masks caret::lift()
## ✖ dplyr::slice()  masks xgboost::slice()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(tidyr)
library(lightgbm)
## 
## Attaching package: 'lightgbm'
## 
## The following object is masked from 'package:dplyr':
## 
##     slice
## 
## The following object is masked from 'package:xgboost':
## 
##     slice
library(data.table)
## 
## Attaching package: 'data.table'
## 
## The following objects are masked from 'package:lubridate':
## 
##     hour, isoweek, mday, minute, month, quarter, second, wday, week,
##     yday, year
## 
## The following object is masked from 'package:purrr':
## 
##     transpose
## 
## The following objects are masked from 'package:dplyr':
## 
##     between, first, last
library(ggplot2)
library(Matrix)
## 
## Attaching package: 'Matrix'
## 
## The following objects are masked from 'package:tidyr':
## 
##     expand, pack, unpack
library(tm) # for text mining
## Loading required package: NLP
## 
## Attaching package: 'NLP'
## 
## The following object is masked from 'package:ggplot2':
## 
##     annotate
library(SnowballC) # for text stemming
library(mice)
## 
## Attaching package: 'mice'
## 
## The following object is masked from 'package:stats':
## 
##     filter
## 
## The following objects are masked from 'package:base':
## 
##     cbind, rbind
library(rpart)
library(rpart.plot)

Data Collection

train = read.csv("/Users/Joyce630/Downloads/usedcars2023 2/analysisData.csv")
train2 = read.csv("/Users/Joyce630/Downloads/usedcars2023 2/scoringData.csv")
train <-bind_rows(train, train2)
#Chech missing value, including NULL, "", empty value
apply(train, MARGIN = 2, FUN = function(x) sum(is.na(x)))
##                       id                make_name               model_name 
##                        0                        0                        0 
##                trim_name                body_type fuel_tank_volume_gallons 
##                        0                        0                        4 
##                fuel_type     highway_fuel_economy        city_fuel_economy 
##                        0                     5398                     5398 
##                    power                   torque             transmission 
##                        0                        0                        0 
##     transmission_display             wheel_system     wheel_system_display 
##                        0                        0                        0 
##         wheelbase_inches      back_legroom_inches     front_legroom_inches 
##                        0                        0                        0 
##            length_inches             width_inches            height_inches 
##                        0                        0                        0 
##              engine_type      engine_displacement               horsepower 
##                        0                      389                      389 
##             daysonmarket              description           exterior_color 
##                        0                        0                        3 
##           interior_color            major_options          maximum_seating 
##                        0                        0                        0 
##                     year                    fleet            frame_damaged 
##                        0                        0                        0 
##         franchise_dealer           franchise_make            has_accidents 
##                        0                        0                        0 
##                    isCab                   is_cpo                   is_new 
##                        0                        0                        0 
##              listed_date            listing_color                  mileage 
##                        0                        0                     2245 
##              owner_count                  salvage            seller_rating 
##                    24741                        0                      677 
##                    price 
##                    10000
#Check ""  total
apply(train,
      MARGIN = 2, 
      FUN = function(x) sum(x == ""))
##                       id                make_name               model_name 
##                        0                        0                        0 
##                trim_name                body_type fuel_tank_volume_gallons 
##                        9                        0                       NA 
##                fuel_type     highway_fuel_economy        city_fuel_economy 
##                      669                       NA                       NA 
##                    power                   torque             transmission 
##                     5538                     6249                     1526 
##     transmission_display             wheel_system     wheel_system_display 
##                     1526                      350                      350 
##         wheelbase_inches      back_legroom_inches     front_legroom_inches 
##                        0                        0                        0 
##            length_inches             width_inches            height_inches 
##                        0                        0                        0 
##              engine_type      engine_displacement               horsepower 
##                      672                       NA                       NA 
##             daysonmarket              description           exterior_color 
##                        0                      983                       NA 
##           interior_color            major_options          maximum_seating 
##                        0                     2655                        0 
##                     year                    fleet            frame_damaged 
##                        0                    23517                    23517 
##         franchise_dealer           franchise_make            has_accidents 
##                        0                    12104                    23517 
##                    isCab                   is_cpo                   is_new 
##                    23517                    46983                        0 
##              listed_date            listing_color                  mileage 
##                        0                        0                       NA 
##              owner_count                  salvage            seller_rating 
##                       NA                    23517                       NA 
##                    price 
##                       NA
#To facilitate filling missing value, convert to all "" to NA
#Fill all empty strings with NA
train <- lapply(train, function(x) ifelse(x == "", NA, x))
train <- as.data.frame(train)

#Now check again
apply(train,
      MARGIN = 2, 
      FUN = function(x) sum(is.na(x)))
##                       id                make_name               model_name 
##                        0                        0                        0 
##                trim_name                body_type fuel_tank_volume_gallons 
##                        9                        0                        4 
##                fuel_type     highway_fuel_economy        city_fuel_economy 
##                      669                     5398                     5398 
##                    power                   torque             transmission 
##                     5538                     6249                     1526 
##     transmission_display             wheel_system     wheel_system_display 
##                     1526                      350                      350 
##         wheelbase_inches      back_legroom_inches     front_legroom_inches 
##                        0                        0                        0 
##            length_inches             width_inches            height_inches 
##                        0                        0                        0 
##              engine_type      engine_displacement               horsepower 
##                      672                      389                      389 
##             daysonmarket              description           exterior_color 
##                        0                      983                        3 
##           interior_color            major_options          maximum_seating 
##                        0                     2655                        0 
##                     year                    fleet            frame_damaged 
##                        0                    23517                    23517 
##         franchise_dealer           franchise_make            has_accidents 
##                        0                    12104                    23517 
##                    isCab                   is_cpo                   is_new 
##                    23517                    46983                        0 
##              listed_date            listing_color                  mileage 
##                        0                        0                     2245 
##              owner_count                  salvage            seller_rating 
##                    24741                    23517                      677 
##                    price 
##                    10000

Data Cleaning

#Fill missing values with NA and count NA within each column
train <- lapply(train, function(x) ifelse(x == "", NA, x))
train <- as.data.frame(train)

apply(train,
      MARGIN = 2, 
      FUN = function(x) sum(is.na(x)))
##                       id                make_name               model_name 
##                        0                        0                        0 
##                trim_name                body_type fuel_tank_volume_gallons 
##                        9                        0                        4 
##                fuel_type     highway_fuel_economy        city_fuel_economy 
##                      669                     5398                     5398 
##                    power                   torque             transmission 
##                     5538                     6249                     1526 
##     transmission_display             wheel_system     wheel_system_display 
##                     1526                      350                      350 
##         wheelbase_inches      back_legroom_inches     front_legroom_inches 
##                        0                        0                        0 
##            length_inches             width_inches            height_inches 
##                        0                        0                        0 
##              engine_type      engine_displacement               horsepower 
##                      672                      389                      389 
##             daysonmarket              description           exterior_color 
##                        0                      983                        3 
##           interior_color            major_options          maximum_seating 
##                        0                     2655                        0 
##                     year                    fleet            frame_damaged 
##                        0                    23517                    23517 
##         franchise_dealer           franchise_make            has_accidents 
##                        0                    12104                    23517 
##                    isCab                   is_cpo                   is_new 
##                    23517                    46983                        0 
##              listed_date            listing_color                  mileage 
##                        0                        0                     2245 
##              owner_count                  salvage            seller_rating 
##                    24741                    23517                      677 
##                    price 
##                    10000
#Assign following new car NA to 0
train$owner_count[is.na(train$owner_count)] <- 0
train$mileage[is.na(train$mileage)] <- 0
train$price[is.na(train$price)] <- 0

#Extract torque and power from its units
extract_values <- function(str) {
  parts <- strsplit(gsub(",", "", str), " ")[[1]]
  return(as.numeric(parts[c(1, 4)]))
}

#Apply the function to the 'power' and 'torque' columns
train <- train %>%
  rowwise() %>%
  mutate(Power_HP = extract_values(power)[1],
         Power_RPM = extract_values(power)[2],
         Torque_lbft = extract_values(torque)[1],
         Torque_RPM = extract_values(torque)[2]) %>%
  ungroup() %>%
  select(-power, -torque)

#Create a custom function to calculate the mode
#First, calculate the overall median or mean, fill hp
overall_median <- median(train$horsepower, na.rm = TRUE)
train <- train %>%
  group_by(trim_name) %>%
  mutate(horsepower = ifelse(is.na(horsepower), mean(horsepower, na.rm = TRUE), horsepower)) %>%
  mutate(horsepower = ifelse(is.na(horsepower), overall_median, horsepower)) %>%
  ungroup()

#calculate mean/median by grouping the same trim model 
library(dplyr)
impute_by_group <- function(data, group_var, target_var, use_median=FALSE) {
  if (use_median) {
    overall_value <- median(data[[target_var]], na.rm = TRUE)
  } else {
    overall_value <- mean(data[[target_var]], na.rm = TRUE)
  }
  data <- data %>%
    group_by(!!sym(group_var)) %>%
    mutate(!!sym(target_var) := ifelse(is.na(.data[[target_var]]),
                                       if (use_median) median(.data[[target_var]], na.rm = TRUE) 
                                       else mean(.data[[target_var]], na.rm = TRUE),
                                       .data[[target_var]])) %>%
    mutate(!!sym(target_var) := ifelse(is.na(.data[[target_var]]), overall_value, .data[[target_var]])) %>%
    ungroup()
  
  
  return(data)
}
train <- impute_by_group(train, 'trim_name', 'engine_displacement', use_median=FALSE)
train <- impute_by_group(train, 'trim_name', 'highway_fuel_economy', use_median=FALSE)
train <- impute_by_group(train, 'trim_name', 'city_fuel_economy', use_median=FALSE)
train <- impute_by_group(train, 'trim_name', 'Power_HP', use_median=FALSE)
train <- impute_by_group(train, 'trim_name', 'Power_RPM', use_median=FALSE)
train <- impute_by_group(train, 'trim_name', 'Torque_lbft', use_median=FALSE)
train <- impute_by_group(train, 'trim_name', 'Torque_RPM', use_median=FALSE)
train <- impute_by_group(train, 'trim_name', 'fuel_tank_volume_gallons', use_median=FALSE)


#calculate mode by grouping the same trim model 
calculate_mode <- function(x) {
  unique_x <- unique(x[!is.na(x)])
  tabulate_match <- tabulate(match(x, unique_x))
  unique_x[which.max(tabulate_match)]
}

# Function to impute missing values by group using the mode
impute_by_group_mode <- function(data, group_var, target_var) {
  overall_mode <- calculate_mode(data[[target_var]])
  data <- data %>%
    group_by(!!sym(group_var)) %>%
    mutate(!!sym(target_var) := ifelse(is.na(.data[[target_var]]), 
                                       calculate_mode(.data[[target_var]]), 
                                       .data[[target_var]])) %>%
    mutate(!!sym(target_var) := ifelse(is.na(.data[[target_var]]), overall_mode, .data[[target_var]])) %>%
    ungroup()
  
  return(data)
}

train <- impute_by_group_mode(train, 'trim_name', 'fuel_type')
train <- impute_by_group_mode(train, 'trim_name', 'transmission')
train <- impute_by_group_mode(train, 'trim_name', 'transmission_display')
train <- impute_by_group_mode(train, 'trim_name', 'wheel_system')
train <- impute_by_group_mode(train, 'trim_name', 'engine_type')

EDA

ggplot(train, aes(x = make_name, y = price)) +
  geom_boxplot() +
  theme(axis.text.x = element_text(angle = 90, hjust = 1)) + 
  labs(title = "Price by Car Make", x = "Car Make", y = "Price")

library(dplyr); library(tidyr)
train |> 
  select(-price)|>
  select_if(is.numeric)|>
  pivot_longer(cols = 1:22,names_to = 'numeric_predictor', values_to = 'values'  )|>
  ggplot(aes(x = values))+
  geom_histogram()+
  facet_wrap(numeric_predictor~., scales = 'free')+
  theme_bw()
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## Warning: Removed 677 rows containing non-finite outside the scale range
## (`stat_bin()`).

One hot Encoding for Make name since we see the price distribution of different models vary

# Convert 'make_name' column in the 'train' dataframe to a factor (categorical variable)
train$make_name <- factor(train$make_name)
# Perform one-hot encoding on the 'make_name' column, excluding the intercept term
encoded_data <- model.matrix(~ make_name - 1, data = train)
encoded_df <- as.data.frame(encoded_data)
# Combine the original 'train' dataframe with the new one-hot encoded dataframe
train <- cbind(train, encoded_df)
# Remove the original 'make_name' column 
train$make_name <- NULL

#same for the body_type
train$body_type <- factor(train$body_type)
encoded_data <- model.matrix(~ body_type - 1, data = train)
encoded_df <- as.data.frame(encoded_data)
train <- cbind(train, encoded_df)
train$body_type <- NULL

#same for the fuel_type
train$fuel_type <- factor(train$fuel_type)
encoded_data <- model.matrix(~ fuel_type - 1, data = train)
encoded_df <- as.data.frame(encoded_data)
train <- cbind(train, encoded_df)
train$fuel_type <- NULL

#same for the transmission_display
train$transmission_display <- factor(train$transmission_display)
encoded_data <- model.matrix(~ transmission_display - 1, data = train)
encoded_df <- as.data.frame(encoded_data)
train <- cbind(train, encoded_df)
train$transmission_display <- NULL

#same for the wheel_system
train$wheel_system <- factor(train$wheel_system)
encoded_data <- model.matrix(~ wheel_system - 1, data = train)
encoded_df <- as.data.frame(encoded_data)
train <- cbind(train, encoded_df)
train$wheel_system <- NULL

#same for the engine_typ
train$engine_type <- factor(train$engine_type)
encoded_data <- model.matrix(~ engine_type - 1, data = train)
encoded_df <- as.data.frame(encoded_data)
train <- cbind(train, encoded_df)
train$engine_type <- NULL

#convert major options to sparse matrix
split_options <- strsplit(gsub("\\[|\\]|'|\"", "", train$major_options), ",\\s*")
all_options <- unlist(split_options)
options_matrix <- table(rep(1:length(split_options), sapply(split_options, length)), all_options)
options_df <- as.data.frame.matrix(options_matrix)
train <- cbind(train, options_df)
train$major_options <- NULL

True and False value missing

#FIll true= 1, false = -1 and NA = 0
fill_and_replace_values <- function(data, column_name) {
  # Convert to lowercase for uniform comparison
  data[[column_name]] <- tolower(as.character(data[[column_name]]))
  data[[column_name]] <- ifelse(data[[column_name]] %in% c('true', 'yes', '1'), 1,
                                ifelse(data[[column_name]] %in% c('false', 'no', '0'), -1, 
                                       data[[column_name]]))
  data[[column_name]][data[[column_name]] == "" | is.na(data[[column_name]])] <- 0
  data[[column_name]] <- as.numeric(data[[column_name]])
  if (any(is.na(data[[column_name]]), na.rm = TRUE)) {
    warning("Some non-numeric values were coerced to NA as they are not recognized boolean values.")
  }
  
  return(data)
}

# Usage:all true false columns
train <- fill_and_replace_values(train, 'isCab')
train <- fill_and_replace_values(train, 'has_accidents')
train <- fill_and_replace_values(train, 'is_cpo')
train <- fill_and_replace_values(train, 'salvage')
train <- fill_and_replace_values(train, 'frame_damaged')
train <- fill_and_replace_values(train, 'fleet')
train <- fill_and_replace_values(train, 'franchise_dealer')
train <- fill_and_replace_values(train, 'is_new')

Feature Engineering

#convert 2023...etc years to how many years from now
convert_year_to_age <- function(dataframe, year_column_name) {
  # Ensure the column exists in the dataframe
  if (!(year_column_name %in% names(dataframe))) {
    stop("Column not found in the dataframe")
  }
  current_year <- as.numeric(format(Sys.Date(), "%Y"))
  dataframe[[year_column_name]] <- current_year - dataframe[[year_column_name]]
  return(dataframe)
}
train <- convert_year_to_age(train, "year")

Develop Models

char_columns <- sapply(train, is.character)
train <- train[, !char_columns]
trainData <- train[train$price != 0, ]
testData <- train[train$price == 0, ]
target_column <- 'price'
x_train <- trainData %>% select(-all_of(target_column))
y_train <- trainData %>% select(all_of(target_column))
x_test <- testData %>% select(-all_of(target_column))
y_test <- testData %>% select(all_of(target_column))

GBM Model (Omitted parameter tuning due to time constraint)

hyperparam <- expand.grid(
  n.trees = c(3000,3500,4000,4500,5000,5500,6000),
  interaction.depth = 1:6,
  shrinkage = seq(0.01,0.1,0.01),
  n.minobsinnode = c(5,10,15)
)
library(caret)
# Create a train control object for cross-validation
trctrl <- trainControl(
  method = "cv",
  number = 5  # Adjust the number of folds as needed
)

# Capture the console output during model training
tuned_model <- train(
  price ~ .,
  data = trainData,
  method = "gbm",
  trControl = trctrl,
  tuneGrid = hyperparam
)

# View the best hyperparameters
print(tuned_model)

# Train the final model with the best hyperparameters
final_gbm_model <- gbm(
  price ~ .,
  data = trainData,
  distribution = "gaussian",
  n.trees = tuned_model$bestTune$n.trees,
  interaction.depth = tuned_model$bestTune$interaction.depth,
  shrinkage = tuned_model$bestTune$shrinkage,
  n.minobsinnode = tuned_model$bestTune$n.minobsinnode
)

After tuning, the best parameters are shown below

library(gbm)
## Loaded gbm 2.1.8.1
set.seed(617)
boost = gbm(price~.,
            data=trainData,
            distribution="gaussian",
            n.trees = 1000,
            interaction.depth = 5,
            shrinkage = 0.08)
pred_train = predict(boost, n.trees=1000)
rmse_train_boost = sqrt(mean((pred_train - trainData$price)^2)); rmse_train_boost
## [1] 3134.479

XGBoost Model (Tune Hyperparameters)

xgb_nround = 5000
xgb_early_stopping_rounds = 3000
xgb_params = list(
  objective = "reg:squarederror", 
  eval_metric = "rmse",
  eta <- c(0.04, 0.06, 0.08, 0.1),
  max_depth <- c(5,6,7,8,9,10),
  subsample <- c(0.5, 0.6, 0.7, 0.8, 0.9, 1),
  colsample_bytree <- 0.5,
  gamma <- c(8, 9, 10),
  alpha <- c(8, 9, 10),
  lambda <- 0)

train_index <- sample.split(trainData$price, SplitRatio = 35000, group = NULL)
train_x <- trainData[train_index, !names(trainData)=='price']
train_y <- trainData[train_index, 'price']
x_test <- testData %>% select(-all_of(target_column))
y_test <- testData %>% select(all_of(target_column))

dtrain <- xgb.DMatrix(data = as.matrix(train_x), label = as.matrix(train_y))
dtest <- xgb.DMatrix(data = as.matrix(test_x), label = as.matrix(test_y))

xgb_model <- xgb.train(
  params = xgb_params,
  data = dtrain,
  nrounds = xgb_nround,
  early_stopping_rounds = xgb_early_stopping_rounds,
  watchlist = list(train = dtrain, test = dtest),
  verbose =2,
  nthread = -1
)
#Initialize the list      
xgb_model <- list()
xgb_rmse <- list()     
cat(
  'best RMSE: ', min(unlist(xgb_rmse)), ' ',
  'eta',xgb_params$eta, ' ',
  'max_depth',xgb_params$max_depth, ' ',
  'subsample',xgb_params$subsample, ' ',
  'colsample_bytree' = xgb_params$colsample_bytree, ' ',
  'gamma',xgb_params$gamma, ' ',
  'alpha',xgb_params$alpha, ' ',
  'lambda',xgb_params$lambda, '\n ',
  '|nround:',unlist(nround_ls), '\n ',
)

XGBoost After tuning

k <- 5
xgb_nround = 2500
xgb_early_stopping_rounds = 10
xgb_params = list(
  objective = "reg:squarederror", 
  eval_metric = "rmse",
  eta = 0.07, # Learning rate
  min_child_weight = 1, # Minimum sum of instance weight needed in a child
  max_depth = 7, # Maximum depth of trees
  subsample = 0.5,  # Subsample ratio of the training instance
  colsample_bytree = 0.5, # Subsample ratio of columns when constructing each tree
  gamma = 8, # Minimum loss reduction required to make a further partition on a leaf node
  alpha = 8  # L1 regularization term on weights
)
x_test <- testData %>% select(-all_of(target_column))
y_test <- testData %>% select(all_of(target_column))
##########Prediction##############
dtrain <- xgb.DMatrix(data = as.matrix(x_train), label = as.matrix (y_train))
dpredict <- xgb.DMatrix(data = as.matrix(x_test))


final_model <- xgb.train(
  params = xgb_params,
  data = dtrain,
  nrounds = xgb_nround,
  early_stopping_rounds = xgb_early_stopping_rounds,
  watchlist = list(train = dtrain),
  verbose = 0,
  nthread = -1
)
pred_train_xgb = predict(final_model, dtrain)
rmse_train_xgboost = sqrt(mean((pred_train_xgb - trainData$price)^2)); rmse_train_xgboost
## [1] 675.2887