e-commerce-shipping-delay-prediction.Rmd

---
title: "Data_Mining_project"
author: "Pratiksha Thakkar"
date: "`r Sys.Date()`"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

```{r, include=FALSE}
library(viridis)
library(dplyr)
library(knitr)
library(kableExtra)

url<-"/Users/newuser/Downloads/E Commerce Shipping Dataset.csv"

library(readr)

dmdata<-read_delim(file=url,delim=",")

```

### Analyzing the historical data of E-Commerce Shipping to identify patterns and factors contributing to late deliveries 


```{r}
names(dmdata)
```




```{r}
dmdatadf<- as.data.frame(dmdata)
head(dmdatadf)
```
```{r}

total_observations <- nrow(dmdatadf)
total_observations

```

```{r}
str(dmdatadf)
```

```{r}
info_data <- data.frame(missing = colSums(is.na(dmdatadf)))
info_data$percent_missing <- (info_data$missing/nrow(dmdatadf)) * 100
info_data$unique <- lapply(dmdatadf, function(column) length(unique(column)))
info_data$dtype <- sapply(dmdatadf, class)
info_data
```

As we can see from the table above, there are no null values in the entire dataset


```{r}
warehouse_summary <- dmdatadf %>%
  group_by(Warehouse_block) %>%
  summarize(count = n(), percentage = round(n() / nrow(dmdatadf), 3))

# Display summary using kable
warehouse_summary %>%
  kable(caption = "Warehouse Block Summary") %>%
  kable_styling(full_width = FALSE) %>%
  column_spec(1:3, bold = TRUE, color = viridis(7), background = "white")
```




```{r}
reached_on_time_summary <- dmdatadf %>%
  group_by(Reached.on.Time_Y.N) %>%
  summarize(count = n(), percentage = round(n() / nrow(dmdatadf), 3))

# Display summary using kable
reached_on_time_summary %>%
  kable(caption = "Reached on Time Summary") %>%
  kable_styling(full_width = FALSE) %>%
  column_spec(1:3, bold = TRUE, color = viridis(7), background = "white")

```



```{r}

shipment_summary <- dmdatadf %>%
  group_by(Mode_of_Shipment) %>%
  summarize(count = n(), percentage = round(n() / nrow(dmdatadf), 3))

# Display summary using kable
shipment_summary %>%
  kable(caption = "Mode of Shipment Summary") %>%
  kable_styling(full_width = FALSE) %>%
  column_spec(1:3, bold = TRUE, color = viridis(7), background = "white")
```







```{r}

library(kableExtra)

# Assuming data_dict is your data frame
dmdatadf %>% select(c(Customer_rating,Cost_of_the_Product,Discount_offered,Weight_in_gms,Discount_offered)) %>% summary() %>%
  kable(caption = "data_summary") %>%
  kable_styling(full_width = FALSE) %>%
   column_spec(1:5, bold = TRUE, color = viridis(1), background = "white")
```




```{r}
data_dict <- data.frame(
  Variable = names(dmdata),
  Description = c(
    "customer ID", 
    "Warehouse_block", 
    "Mode of Shipment", 
    "No. of calls made to customer care", 
    "Customer ratings out of 5", 
    "Cost of the product", 
    "Number of Prior purchases ",
    "Product_importance",
    "Gender ",
    "Discount offered ",
    "Weight of the product in grams",
    "If the product reached on time or not"),
   Example = c("1,2,3 ...", "A,B,C,","Ship,Flight,road" ," <2 low , 3-4 medium, >5 high", "0, 1, 2 ...", "1 to 5", "96 -310 USD  ..", "1,2,3...", "Male/Female", "1-65 USD", "1001 gms to 7846 gms", "Yes/No")
)



library(kableExtra)

# Assuming data_dict is your data frame
data_dict %>%
  kable(caption = "Data Dictionary") %>%
  kable_styling(full_width = FALSE) %>%
   column_spec(1:3, bold = TRUE, color = viridis(1), background = "white")
```


#### data visualization


```{r}
library(ggplot2)
library(viridis)

ggplot(dmdatadf, aes(x = Warehouse_block)) +
  geom_bar(fill =magma(5)) +
  labs(title = "Distribution of Warehouse Blocks")

```


```{r}
ggplot(dmdatadf, aes(x = Mode_of_Shipment)) +
  geom_bar(fill = magma(3)) +
  labs(title = "Distribution of Shipment Modes")
```


```{r}
ggplot(dmdatadf, aes(x =Customer_care_calls )) +
  geom_histogram(fill = "salmon", bins = 5) +
  labs(title = "Distribution of Customer Care Calls")
```



```{r}
ggplot(dmdatadf, aes(x = Customer_care_calls)) +
  geom_bar(fill = "salmon") +
  labs(title = "Distribution of Customer Care Calls")
```


```{r}


```
```{r}
ggplot(dmdatadf, aes(x = Customer_rating)) +
  geom_bar(fill = "salmon") +
  labs(title = "Distribution of Customer Ratings")
```


```{r}
ggplot(dmdatadf, aes(x = Cost_of_the_Product)) +
  geom_histogram(fill = "salmon", bins = 15) +
  labs(title = "Distribution of Product Costs")

```

```{r}
ggplot(dmdatadf, aes(x = Prior_purchases)) +
  geom_histogram(fill = "salmon", bins = 10) +
  labs(title = "Distribution of Prior Purchases")

```

```{r}
ggplot(dmdatadf, aes(x = Prior_purchases)) +
  geom_bar(fill = "salmon", bins = 10) +
  labs(title = "Distribution of Prior Purchases")
```


```{r}
ggplot(dmdatadf, aes(x = Product_importance)) +
  geom_bar(fill = "salmon") +
  labs(title = "Distribution of Product Importance")

```
```{r}
ggplot(dmdatadf, aes(x = Gender)) +
  geom_bar(fill = "salmon") +
  labs(title = "Distribution of Gender")

```


```{r}
ggplot(dmdatadf, aes(x = Discount_offered)) +
  geom_histogram(fill = "salmon", bins = 10) +
  labs(title = "Distribution of Discounts Offered")

```

```{r}
ggplot(dmdatadf, aes(x = Weight_in_gms)) +
  geom_histogram(fill = "salmon", bins = 40) +
  labs(title = "Distribution of Product Weights")

```

```{r}
library(ggplot2)


ggplot(dmdatadf, aes(x = factor(Reached.on.Time_Y.N), fill = factor(Reached.on.Time_Y.N))) +
  geom_bar() +
  scale_x_discrete(labels = c("0" = "On Time", "1" = "Not On Time")) +
  labs(title = "Distribution of On-Time Delivery",
       x = "Reached on Time", y = "Frequency")


```




```{r}
library(ggplot2)

total_reviews <- aggregate(Customer_rating ~ Reached.on.Time_Y.N, data = dmdatadf, FUN = length)

ggplot(dmdatadf, aes(x = Reached.on.Time_Y.N, fill = factor(Customer_rating))) +
  geom_bar(position = "stack", alpha = 0.7) +
  geom_text(data = total_reviews, aes(label = Customer_rating, y = Customer_rating, group = Reached.on.Time_Y.N),
            position = position_stack(vjust = 0.5), color = "black", size = 3) +
  labs(title = "Distribution of Customer Ratings by On-Time Delivery",
       x = "Reached on Time", y = "Frequency", fill = "Customer Rating")
```

which mode of shipment is used as per the weight of product and does it affect the time it reaches.
```{r}
ggplot(dmdatadf, aes(x = factor(Mode_of_Shipment), fill = factor(Reached.on.Time_Y.N))) +
  geom_bar() +
  scale_x_discrete(labels = c("0" = "Flight", "1" = "Road" ,"2" = "Ship")) +
  labs(title = "Distribution of On-Time Delivery",
       x = "Mode of shipment", y = "Weight of product in gms")
```


```{r}
min_rating <- min(dmdatadf$Customer_rating)
max_rating <- max(dmdatadf$Customer_rating)

# Create a histogram with dynamically calculated breaks
hist(dmdatadf$Customer_rating,
     main = "Histogram of Customer Ratings",
     xlab = "Rating",
     ylab = "Frequency",
     col = "skyblue",
     border = "black",
     breaks = seq(min_rating - 0.5, max_rating + 0.5, by = 1),  # Adjust breaks to cover the range of values
     ylim = c(0, max(table(dmdatadf$Customer_rating)) + 1)  # Set y-axis limits
)
str(dmdatadf)
```

##### labelling.  
```{r}
library(caret)
library(fastDummies)
categorical_cols <- c("Warehouse_block", "Mode_of_Shipment")


dmdatadf <- dummy_cols(dmdatadf,select_columns = categorical_cols)
#dmdatadf = data.frame(predict(dmy, newdata=dmdatadf))
#data_combined <- cbind(dmdatadf[, setdiff(names(dmdatadf), categorical_cols), drop = TRUE])
dmdatadf <- dmdatadf[, !names(dmdatadf) %in% categorical_cols]
data_combined <- subset(dmdatadf, select=-c(ID))

data_combined[, grep(paste(categorical_cols, collapse="|"), names(data_combined))] <- lapply(data_combined[, grep(paste(categorical_cols, collapse="|"), names(data_combined))], factor)




data_combined$Product_importance = factor(ifelse(data_combined$Product_importance == "low", 0,ifelse (data_combined$Product_importance == "medium", 1, ifelse (data_combined$Product_importance == "high",2,3))))

data_combined$Prior_purchases = factor(ifelse(data_combined$Prior_purchases == 2, 0,ifelse (data_combined$Prior_purchases == 3, 1, 2 )))


data_combined$Weight_in_gms <- cut(data_combined$Weight_in_gms, 
                                    breaks = c(-Inf, 2000, 4000, Inf), 
                                    labels = c("0", "1", "2"), 
                                   right = FALSE)
str(data_combined)

```


```{r}
library(ggplot2)

# Apply logarithmic transformation to Discount_offered variable
data_combined$log_discount <- log(data_combined$Discount_offered + 1)  # Adding 1 to handle zero values if present

# Ploing histogram of the transformed variable
ggplot(data_combined, aes(x = log_discount)) +
  geom_histogram(fill = "salmon", bins = 10) +
  labs(title = "Distribution of Log-transformed Discounts Offered")
```


```{r}
# Load necessary libraries
library(caret)
library(lattice)
# Set seed for reproducibility
set.seed(123)

# Split the data into partition (50%), validation (30%), and holdout (20%)
partitionIndex <- createDataPartition(data_combined$Reached.on.Time_Y.N, p = 0.5, list = FALSE)
partition_data <- data_combined[partitionIndex, ]
remaining_data <- data_combined[-partitionIndex, ]

validationIndex <- createDataPartition(data_combined$Reached.on.Time_Y.N, p = 0.3, list = FALSE)
validation_data <- remaining_data[validationIndex, ]
holdout_data <- remaining_data[-validationIndex, ]

# Print the dimensions of each subset
cat("Partition data:", dim(partition_data), "\n")
cat("Validation data:", dim(validation_data), "\n")
cat("Holdout data:", dim(holdout_data), "\n")

```


```{r}
# Load the rpart package
library(rpart)
library(partykit)
library(rpart.plot)
partition_data$Reached.on.Time_Y.N <- factor(partition_data$Reached.on.Time_Y.N)
validation_data$Reached.on.Time_Y.N <- factor(validation_data$Reached.on.Time_Y.N)
holdout_data$Reached.on.Time_Y.N <- factor(holdout_data$Reached.on.Time_Y.N)

# Train the classification tree using partition_data
tree_model <- rpart(Reached.on.Time_Y.N ~ ., data = partition_data, method = "class")

# plotting tree
prp(tree_model, type = 1, extra = 1, under = TRUE, split.font = 1, varlen = -10)

# Make predictions on the validation data
tree_predictions <- predict(tree_model, validation_data, type = "class")

# Evaluate the performance on the validation data
accuracy_validation <- sum(tree_predictions == validation_data$Reached.on.Time_Y.N, na.rm = TRUE) / length(tree_predictions)
cat("Accuracy on validation data:", accuracy_validation, "\n")

# Make predictions on the holdout data
tree_holdout_predictions <- predict(tree_model, holdout_data, type = "class")

# Evaluate the performance on the holdout data
accuracy_holdout <- sum(tree_holdout_predictions == holdout_data$Reached.on.Time_Y.N, na.rm = TRUE) / length(tree_holdout_predictions)
cat("Accuracy on holdout data:", accuracy_holdout, "\n")

  
```


```{r}
# Train the Naive Bayes classifier using partition_data
nb_model <- e1071::naiveBayes(Reached.on.Time_Y.N ~ ., data = partition_data)

# Make predictions on the validation data
nb_predictions <- predict(nb_model, validation_data)

# Evaluate the performance on the validation data
accuracy_validation <- mean(nb_predictions == validation_data$Reached.on.Time_Y.N, na.rm = TRUE)
cat("Accuracy on validation data:", accuracy_validation, "\n")

# Make predictions on the holdout data
nb_holdout_predictions <- predict(nb_model, holdout_data)

# Evaluate the performance on the holdout data
accuracy_holdout <- mean(nb_holdout_predictions == holdout_data$Reached.on.Time_Y.N, na.rm = TRUE)
cat("Accuracy on holdout data:", accuracy_holdout, "\n")
```


```{r}
# Train the logistic regression model using partition_data

partition_data <- subset(partition_data, select=-c(Warehouse_block_F, Mode_of_Shipment_Ship))
logistic_model <- glm(Reached.on.Time_Y.N ~ ., data = partition_data, family = binomial)
#trControl <- caret::trainControl(method="cv", number=5, allowParallel=TRUE)
#logit.reg <- caret::train(Reached.on.Time_Y.N ~ ., data=partition_data,trcontrol=trcontrol,
                       #fit logistic regression with a generalized linear model
                       #method="glm", family="binomial")


options(scipen = 999)
summary(logistic_model)


# Make predictions on the validation data
logistic_predictions <- predict(logistic_model, validation_data, type = "response")

# Convert predicted probabilities to binary predictions (0 or 1)
logistic_predictions_binary <- ifelse(logistic_predictions > 0.5, 1, 0)

# Evaluate the performance on the validation data
accuracy_validation <- mean(logistic_predictions_binary == validation_data$Reached.on.Time_Y.N, na.rm = TRUE)
cat("Accuracy on validation data:", accuracy_validation, "\n")

# Make predictions on the holdout data
logistic_holdout_predictions <- predict(logistic_model, holdout_data, type = "response")

# Convert predicted probabilities to binary predictions (0 or 1)
logistic_holdout_predictions_binary <- ifelse(logistic_holdout_predictions > 0.5, 1, 0)


# Evaluate the performance on the holdout data
accuracy_holdout <- mean(logistic_holdout_predictions_binary == holdout_data$Reached.on.Time_Y.N, na.rm = TRUE)
cat("Accuracy on holdout data:", accuracy_holdout, "\n")
```


```{r}
# confusion matrix before stepwise
# Confusion matrix for validation data 
confusion_matrix_validation <- confusionMatrix(factor(logistic_predictions_binary), factor(validation_data$Reached.on.Time_Y.N))
print(confusion_matrix_validation)

#Confusion matrix for holdout data
confusion_matrix_holdout <- confusionMatrix(factor(logistic_holdout_predictions_binary), factor(holdout_data$Reached.on.Time_Y.N))
print(confusion_matrix_holdout)
```


```{r}
#THIS IS AFTER STEPWISE VALIDATION AND HOLDOUT CONFUSION MATRIX 

final_model <- glm(Reached.on.Time_Y.N ~ Customer_care_calls + Customer_rating + Cost_of_the_Product + Prior_purchases + Product_importance + Discount_offered + Weight_in_gms + log_discount, data=partition_data, family = binomial)
summary(final_model)
prediction_validate <- predict(final_model, newdata = validation_data, type = "response")

# Convert predicted probabilities to predicted classes (0 or 1) based on a threshold (e.g., 0.5)
threshold <- 0.5
predicted_classes <- ifelse(prediction_validate > 0.5, 1, 0)
validation_data$predicted <- ifelse(predicted_classes >= threshold, 1, 0)

# Confusion matrix for validation data 
confusion_matrix_validation <- confusionMatrix(factor(predicted_classes), factor(validation_data$Reached.on.Time_Y.N))
print(confusion_matrix_validation)

# Confusion matrix for holdout data
holdout_validate <- predict(final_model, newdata = holdout_data, type = "response")

# Convert predicted probabilities to predicted classes (0 or 1) based on a threshold (e.g., 0.5)
threshold <- 0.5
holdout_classes <- ifelse(holdout_validate > 0.5, 1, 0)
holdout_data$predicted <- ifelse(holdout_validate >= threshold, 1, 0)


#Confusion matrix for holdout data
confusion_matrix_holdout <- confusionMatrix(factor(holdout_classes), factor(holdout_data$Reached.on.Time_Y.N))
print(confusion_matrix_holdout)



```
```{r}
# Fit the logistic regression model
final_model <- glm(Reached.on.Time_Y.N ~ Customer_care_calls + Customer_rating + Cost_of_the_Product + Prior_purchases + Product_importance + Discount_offered + Weight_in_gms + log_discount, data = partition_data, family = binomial)

# Summary of the final model
summary(final_model)

# Predict probabilities for the validation data
prediction_validate <- predict(final_model, newdata = validation_data, type = "response")

# Initialize vectors to store threshold and specificity values
threshold_values <- c()
specificity_values <- c()

# Loop through different threshold values
for (threshold in seq(0.1, 0.9, by = 0.1)) {
  # Convert predicted probabilities to predicted classes based on the threshold
  predicted_classes <- ifelse(prediction_validate > threshold, 1, 0)
  
  # Calculate confusion matrix for validation data
  confusion_matrix_validation <- confusionMatrix(factor(predicted_classes), factor(validation_data$Reached.on.Time_Y.N))
  
  # Calculate specificity
  specificity <- confusion_matrix_validation$byClass["Specificity"]
  
  # Store threshold and specificity values
  threshold_values <- c(threshold_values, threshold)
  specificity_values <- c(specificity_values, specificity)
}

# Combine threshold and specificity values into a data frame
threshold_specificity_df <- data.frame(Threshold = threshold_values, Specificity = specificity_values)

# Find the threshold value corresponding to the highest specificity
best_threshold <- threshold_specificity_df$Threshold[which.max(threshold_specificity_df$Specificity)]

# Print the best threshold
cat("Best Threshold (Highest Specificity):", best_threshold, "\n")

# Predict probabilities for the holdout data
holdout_validate <- predict(final_model, newdata = holdout_data, type = "response")

# Convert predicted probabilities to predicted classes based on the best threshold
holdout_classes <- ifelse(holdout_validate > best_threshold, 1, 0)

# Confusion matrix for holdout data using the best threshold
confusion_matrix_holdout <- confusionMatrix(factor(holdout_classes), factor(holdout_data$Reached.on.Time_Y.N))
print(confusion_matrix_holdout)

print(best_threshold)
```



```{r}
#THIS IS BEFORE STEPWISE NAIVE BAYES VALIDATION AND HOLDOUT CONFUSION MATRIX 


# Train the Naive Bayes classifier using partition_data
#data_combined$Discount_offered <- exp(data_combined$log_discount) - 1

nb_model <- e1071::naiveBayes(Reached.on.Time_Y.N ~ ., data = partition_data)

# Make predictions on the validation data
nb_predictions <- predict(nb_model, validation_data)

# Confusion matrix for validation data 
confusion_matrix_validation_nb <- confusionMatrix(nb_predictions, validation_data$Reached.on.Time_Y.N)
print(confusion_matrix_validation_nb)


#summary(nb_model)
# Specify a custom threshold for classification (e.g., 0.6)
prediction_holdout_log_model <- predict(nb_model, newdata = holdout_data, threshold = 0.5)

confusion_matrix_validation_nb <- confusionMatrix(prediction_holdout_log_model, holdout_data$Reached.on.Time_Y.N)
print(confusion_matrix_validation_nb)


```
```{r}
library(e1071)
library(caret)

# Train the Naive Bayes classifier using partition_data
nb_model <- naiveBayes(Reached.on.Time_Y.N ~ ., data = partition_data)

# Make predictions on the validation data
nb_predictions_val <- predict(nb_model, newdata = validation_data)

# Confusion matrix for validation data 
confusion_matrix_validation_nb <- confusionMatrix(nb_predictions_val, validation_data$Reached.on.Time_Y.N)
print("Validation Confusion Matrix:")
print(confusion_matrix_validation_nb)

# Initialize vectors to store results for different thresholds
threshold_values <- seq(0.1, 0.9, by = 0.1)  # Adjust the step size as needed
sensitivity_values <- numeric(length(threshold_values))
specificity_values <- numeric(length(threshold_values))
accuracy_values <- numeric(length(threshold_values))

# Loop through different threshold values
for (i in seq_along(threshold_values)) {
  # Make predictions on the holdout data with the current threshold
  prediction_holdout <- predict(nb_model, newdata = holdout_data, threshold = threshold_values[i])
  
  # Calculate confusion matrix for the holdout data
  cm_holdout <- confusionMatrix(prediction_holdout, holdout_data$Reached.on.Time_Y.N)
  
  # Store sensitivity, specificity, and accuracy for the current threshold
  sensitivity_values[i] <- cm_holdout$byClass["Sensitivity"]
  specificity_values[i] <- cm_holdout$byClass["Specificity"]
  accuracy_values[i] <- cm_holdout$overall["Accuracy"]
  
  # Print confusion matrix for the holdout data with the current threshold
  cat("\nThreshold:", threshold_values[i], "\n")
  print("Holdout Confusion Matrix:")
  print(cm_holdout)
}

# Find the best threshold based on specificity
best_threshold_spec <- threshold_values[which.max(specificity_values)]

# Print the best threshold based on specificity
cat("\nBest Threshold (Highest Specificity):", best_threshold_spec, "\n")

```


```{r}
#THIS IS  NAIVE BAYES AFTER STEPWISE CONFUSION MATRIX FOR VALIDATION AND HOLDOUT

# Train the Naive Bayes classifier using partition_data
nb_final_model <- e1071::naiveBayes(Reached.on.Time_Y.N ~ Customer_care_calls + Customer_rating + Cost_of_the_Product + Prior_purchases +  Product_importance + Discount_offered + Weight_in_gms + log_discount, data = partition_data)

# Make predictions on the validation data
nb_predictions <- predict(nb_final_model, validation_data)

confusion_matrix_validation_nb <- confusionMatrix(nb_predictions, validation_data$Reached.on.Time_Y.N)
print(confusion_matrix_validation_nb)


#summary(nb_model)
# Convert predicted probabilities to predicted classes (0 or 1) based on a threshold (e.g., 0.5)
prediction_holdout_log_model <- predict(nb_final_model, newdata = holdout_data)

confusion_matrix_validation_nb <- confusionMatrix(prediction_holdout_log_model, holdout_data$Reached.on.Time_Y.N)
print(confusion_matrix_validation_nb)


```
```{r}
library(e1071)
library(caret)

# Train the Naive Bayes classifier using partition_data
nb_final_model <- e1071::naiveBayes(Reached.on.Time_Y.N ~ Customer_care_calls + Customer_rating + Cost_of_the_Product + Prior_purchases +  Product_importance + Discount_offered + Weight_in_gms + log_discount, data = partition_data)

# Make predictions on the validation data
nb_predictions <- predict(nb_final_model, validation_data)

# Confusion matrix for validation data 
confusion_matrix_validation_nb <- confusionMatrix(nb_predictions_val, validation_data$Reached.on.Time_Y.N)
print("Validation Confusion Matrix:")
print(confusion_matrix_validation_nb)

# Initialize vectors to store results for different thresholds
threshold_values <- seq(0.1, 0.9, by = 0.1)  # Adjust the step size as needed
sensitivity_values <- numeric(length(threshold_values))
specificity_values <- numeric(length(threshold_values))
accuracy_values <- numeric(length(threshold_values))

# Loop through different threshold values
for (i in seq_along(threshold_values)) {
  # Make predictions on the holdout data with the current threshold
  prediction_holdout <- predict(nb_model, newdata = holdout_data, threshold = threshold_values[i])
  
  # Calculate confusion matrix for the holdout data
  cm_holdout <- confusionMatrix(prediction_holdout, holdout_data$Reached.on.Time_Y.N)
  
  # Store sensitivity, specificity, and accuracy for the current threshold
  sensitivity_values[i] <- cm_holdout$byClass["Sensitivity"]
  specificity_values[i] <- cm_holdout$byClass["Specificity"]
  accuracy_values[i] <- cm_holdout$overall["Accuracy"]
  
  # Print confusion matrix for the holdout data with the current threshold
  cat("\nThreshold:", threshold_values[i], "\n")
  print("Holdout Confusion Matrix:")
  print(cm_holdout)
}

# Find the best threshold based on specificity
best_threshold_spec <- threshold_values[which.max(specificity_values)]

# Print the best threshold based on specificity
cat("\nBest Threshold (Highest Specificity):", best_threshold_spec, "\n")

```



```{r}
#THIS IS CLASSIFICATION TREE BEFORE STEPWISE VALIDATION AND HOLDOUT CONFUSION MATRIX:

library(rpart)
partition_data$Reached.on.Time_Y.N <- factor(partition_data$Reached.on.Time_Y.N)
validation_data$Reached.on.Time_Y.N <- factor(validation_data$Reached.on.Time_Y.N)
holdout_data$Reached.on.Time_Y.N <- factor(holdout_data$Reached.on.Time_Y.N)

# Train the classification tree using partition_data
tree_model <- rpart(Reached.on.Time_Y.N ~ ., data = partition_data, method = "class")

# Make predictions on the validation data
tree_predictions_val <- predict(tree_model, validation_data, type = "class")

confusion_matrix_validation_class <- confusionMatrix(tree_predictions_val, validation_data$Reached.on.Time_Y.N)
print(confusion_matrix_validation_class)

tree_predictions_hold <- predict(tree_model, holdout_data, type = "class")

confusion_matrix_validation_class <- confusionMatrix(tree_predictions_hold, holdout_data$Reached.on.Time_Y.N)
print(confusion_matrix_validation_class)



```


```{r}
#THIS IS CLASSIFICATION TREE AFTER STEPWISE VALIDATION AND HOLDOUT CONFUSION MATRIX:

tree_model <- rpart(Reached.on.Time_Y.N ~ Customer_care_calls + Customer_rating + Cost_of_the_Product + Prior_purchases + Product_importance + Discount_offered + Weight_in_gms + log_discount, data = partition_data, method = "class")

# Make predictions on the validation data
tree_predictions_val <- predict(tree_model, validation_data, type = "class")

confusion_matrix_validation_class <- confusionMatrix(tree_predictions_val, validation_data$Reached.on.Time_Y.N)
print(confusion_matrix_validation_class)

tree_predictions_hold <- predict(tree_model, holdout_data, type = "class")

confusion_matrix_validation_class <- confusionMatrix(tree_predictions_hold, holdout_data$Reached.on.Time_Y.N)
print(confusion_matrix_validation_class)



```
```{r}
library(rpart)
library(caret)

# Fit the classification tree model using partition_data
tree_model <- rpart(Reached.on.Time_Y.N ~ Customer_care_calls + Customer_rating + Cost_of_the_Product + Prior_purchases + Product_importance + Discount_offered + Weight_in_gms + log_discount, data = partition_data, method = "class")

# Make predictions on the validation data
tree_predictions_val <- predict(tree_model, newdata = validation_data, type = "class")

# Ensure consistency in factor levels
levels(tree_predictions_val) <- levels(validation_data$Reached.on.Time_Y.N)

# Confusion matrix for validation data 
confusion_matrix_validation_tree <- confusionMatrix(tree_predictions_val, validation_data$Reached.on.Time_Y.N)
print("Validation Confusion Matrix:")
print(confusion_matrix_validation_tree)

# Initialize vectors to store results for different thresholds
threshold_values <- seq(0.1, 0.9, by = 0.1)  # Adjust the step size as needed
sensitivity_values <- numeric(length(threshold_values))
specificity_values <- numeric(length(threshold_values))
accuracy_values <- numeric(length(threshold_values))

# Loop through different threshold values
for (i in seq_along(threshold_values)) {
  # Make predictions on the holdout data with the current threshold
  tree_predictions_hold <- predict(tree_model, newdata = holdout_data, type = "prob")[, 2]
  prediction_holdout <- ifelse(tree_predictions_hold > threshold_values[i], levels(holdout_data$Reached.on.Time_Y.N)[2], levels(holdout_data$Reached.on.Time_Y.N)[1])
  
  # Ensure consistency in factor levels
  prediction_holdout <- factor(prediction_holdout, levels = levels(holdout_data$Reached.on.Time_Y.N))
  
  # Calculate confusion matrix for the holdout data
  cm_holdout <- confusionMatrix(prediction_holdout, holdout_data$Reached.on.Time_Y.N)
  
  # Store sensitivity, specificity, and accuracy for the current threshold
  sensitivity_values[i] <- cm_holdout$byClass["Sensitivity"]
  specificity_values[i] <- cm_holdout$byClass["Specificity"]
  accuracy_values[i] <- cm_holdout$overall["Accuracy"]
  
  # Print confusion matrix for the holdout data with the current threshold
  cat("\nThreshold:", threshold_values[i], "\n")
  print("Holdout Confusion Matrix:")
  print(cm_holdout)
}

# Find the best threshold based on specificity
best_threshold_spec <- threshold_values[which.max(specificity_values)]

# Print the best threshold based on specificity
cat("\nBest Threshold (Highest Specificity):", best_threshold_spec, "\n")

```


    



#### AIC and step wise
```{r}
data_combined$Reached.on.Time_Y.N <- factor(data_combined$Reached.on.Time_Y.N, levels = c(0, 1))
trControl <- caret::trainControl(method="none")
model <- caret::train(Reached.on.Time_Y.N ~ ., data=data_combined, trControl=trControl,
                      # select backward elmination
                      method="glmStepAIC", direction='backward')
coef(model$finalModel)

model <- caret::train(Reached.on.Time_Y.N ~ ., data=data_combined, trControl=trControl,
                      method="glmStepAIC", direction='forward')
coef(model$finalModel)
model <- caret::train(Reached.on.Time_Y.N ~ ., data=data_combined, trControl=trControl,
                      method="glmStepAIC", direction='both')
coef(model$finalModel)
```


##### Cross-validation
```{r}
set.seed(1)
library(caret)
# define 5-fold
trControl <- caret::trainControl(method="cv", number=5, allowParallel=TRUE)
model <- caret::train(Reached.on.Time_Y.N ~ ., data=data_combined,
                      method="glm",  # specify the model
                      trControl=trControl)
model
```

#### Exaustive
```{r}
library(fastDummies)
library(leaps)

search <- regsubsets(Reached.on.Time_Y.N ~ ., data=data_combined, nbest=1,
                     nvmax=ncol(data_combined), method="exhaustive")
sum <- summary(search)

# show models
sum$which
sum$rsq
sum$adjr2
sum$cp
optimal <- which.min(sum$cp)

# determine the variable names for the optimal model
X <- summary(search)$which[, -1]  # information about included predictors
xvars <- dimnames(X)[[2]] ## column names (all covariates except intercept)
xvars <- xvars[X[optimal,]]

optimal <- which.min(sum$cp)

# determine the variable names for the optimal model
X <- summary(search)$which[, -1]  # information about included predictors
xvars <- dimnames(X)[[2]] ## column names (all covariates except intercept)
xvars <- xvars[X[optimal,]]


xvars <- c("Customer_care_calls","Customer_rating","Cost_of_the_Product","Prior_purchases","Product_importance","Discount_offered","Weight_in_gms")

# rebuild model for best predictor set
set.seed(1)
trControl <- caret::trainControl(method="cv", number=5, allowParallel=TRUE)
model <- caret::train(Reached.on.Time_Y.N ~ ., data=data_combined[, c("Reached.on.Time_Y.N", xvars)],
                      method="glm",  # specify the model
                      trControl=trControl)
model

```


```{r}

```



```{r}
anova(logistic_model, test = 'Chisq')
```

<!-- Recommended Model: Decision Tree (or Random Forest)For this e-commerce shipping dataset, the Decision Tree (or an ensemble tree like Random Forest) is the best choice.Why Decision Tree is the Best Choice1. Superior Predictive AccuracyDecision Tree Holdout Accuracy: 68.36%Logistic Regression Holdout Accuracy: 65.82%Naïve Bayes Holdout Accuracy: 66.68%The tree model consistently outperforms linear and probabilistic baselines on unseen holdout test data.2. Non-Linear Feature InteractionsLinear models like Logistic Regression struggle with this dataset because key predictors—specifically Discount_offered and Weight_in_gms—do not have a strictly linear relationship with shipping delays.For instance:Discounts: Shipments with low to moderate discounts behave normally, but high discounts ($\ge 10\%$) trigger sharp, non-linear surges in shipping delays due to promotional volume spikes.Product Weights: Deliveries experience non-linear risk spikes in specific weight ranges (e.g., packages around 2,000–4,000g vs. light items under 2,000g).Decision trees handle these sharp step-function cutoffs naturally by creating conditional branches (e.g., if Discount > 10 AND Weight between 2000-4000g).3. Operational InterpretabilityUnlike "black-box" models, decision trees produce clear rule sets that logistics managers can easily turn into operational decisions (e.g., "Flag any promotional order over $200 with a discount above 10% for priority warehouse dispatch").Summary Recommendation MatrixMetric / NeedLogistic RegressionNaïve BayesDecision TreeAccuracy65.8%66.7%68.4%Handles Non-Linear Features❌ Poor⚠️ ModerateYES ExcellentActionable Business Rules⚠️ Coefficients only❌ Probability tablesYES Visual tree rules -->



```{r}
library(pROC)

# 1. Compute probabilities on your holdout data
prob_logistic <- predict(final_model, newdata = holdout_data, type = "response")
prob_tree     <- predict(tree_model, newdata = holdout_data, type = "prob")[, 2]

# 2. Calculate ROC curves
roc_glm  <- roc(holdout_data$Reached.on.Time_Y.N, prob_logistic)
roc_tree <- roc(holdout_data$Reached.on.Time_Y.N, prob_tree)

# 3. Save directly to PNG (without needing a 'roc_plot' object)
dir.create("images", showWarnings = FALSE)
png("images/roc_curve_comparison.png", width = 2100, height = 1800, res = 300, bg = "white")

# Plot standard ROC curve
plot(roc_tree, col = "#21908CFF", main = "ROC Curve Comparison (Holdout Set)", lwd = 2)
plot(roc_glm, col = "#440154FF", add = TRUE, lwd = 2)

legend("bottomright", 
       legend = c(paste("Decision Tree (AUC =", round(auc(roc_tree), 3), ")"),
                  paste("Logistic Reg  (AUC =", round(auc(roc_glm), 3), ")")),
       col = c("#21908CFF", "#440154FF"), 
       lwd = 2)

dev.off()
```