Individual Assignment JE
Jonathan Edel
2025-09-22
Graded Assignment
In this assignment, we are given a typical e-commerce dataset from an appearal seller, and tasked with predicting whether an item will be returned. We will attempt to train a predictive model to best predict the likelihood of an order being returned, based on the information available about the order.
Data cleanup
We will start with data cleanup. It is generally accepted that 50=80% of the time spent on data science projects is devoted to data cleanup.
We will go variable by variable in this exercise. Alongside this markup, I have also submitted the full data cleanup script I am using. Please note you will need to change the WD where relevant.
After cursory initial investigation, we start with the basic step of merging the user and orders datasets. This is possible, as they both contain the same unique order identifier (order_item_id). We also load all packages we will use in data cleanup.
#Load the necessary packages
library(lubridate) # date management## Warning: package 'lubridate' was built under R version 4.4.3##
## Attaching package: 'lubridate'## The following objects are masked from 'package:base':
##
## date, intersect, setdiff, unionlibrary(ggplot2) # charts## Warning: package 'ggplot2' was built under R version 4.4.3library(dplyr) # data handling## Warning: package 'dplyr' was built under R version 4.4.3##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(glmnet) # implementation of the methods## Warning: package 'glmnet' was built under R version 4.4.3## Loading required package: Matrix## Loaded glmnet 4.1-10library(caret) # machine-learning workflow## Warning: package 'caret' was built under R version 4.4.3## Loading required package: lattice# set working directory:
setwd("C:/Users/jedel/OneDrive/Desktop/Data for ML class/Individual Assignment")
#Load data
users <- read.csv("users.csv", sep = ",")
orders <- read.csv("orders.csv", sep = ",")
# Merge data into a single source
orders_full <- merge(orders, users, by = "user_id")In addition to looking through the table, we will also compute simple summary statistics. This is an early step when examining a new dataset. Not all details of preparatory data examination for this project are included in this markup.
summary(orders_full)## user_id order_item_id order_date delivery_date
## Min. : 9 Min. : 1 Length:75000 Length:75000
## 1st Qu.:14674 1st Qu.:18751 Class :character Class :character
## Median :30479 Median :37501 Mode :character Mode :character
## Mean :25739 Mean :37501
## 3rd Qu.:36477 3rd Qu.:56250
## Max. :43801 Max. :75000
## item_id item_size item_color brand_id
## Min. : 1.0 Length:75000 Length:75000 Min. : 1.00
## 1st Qu.: 175.0 Class :character Class :character 1st Qu.: 5.00
## Median : 479.0 Mode :character Mode :character Median : 24.00
## Mean : 784.7 Mean : 28.83
## 3rd Qu.:1516.0 3rd Qu.: 43.00
## Max. :2224.0 Max. :137.00
## item_price return user_title user_dob
## Min. : 0.00 Min. :0.0000 Length:75000 Length:75000
## 1st Qu.: 29.90 1st Qu.:0.0000 Class :character Class :character
## Median : 49.90 Median :0.0000 Mode :character Mode :character
## Mean : 60.67 Mean :0.4501
## 3rd Qu.: 74.90 3rd Qu.:1.0000
## Max. :999.00 Max. :1.0000
## user_state user_reg_date
## Length:75000 Length:75000
## Class :character Class :character
## Mode :character Mode :character
##
##
## Next, we will convert returns to a factor, and set all dates to ymd format. This is directly from the assignment. We will calculate delivery time once we get to order date and delivery date.
# Convert all dates into proper proper date format
orders_full$order_date <- ymd(orders_full$order_date)
orders_full$delivery_date <- ymd(orders_full$delivery_date)
orders_full$user_dob <- ymd(orders_full$user_dob)
orders_full$user_reg_date <- ymd(orders_full$user_reg_date)
#convert returns to "Yes" or "No"
orders_full$return <- factor(orders_full$return,
levels = c(0, 1),
labels = c("No", "Yes"))Now, we will start looking at individual variables. These are in no particular order. Let’s start with item price.
Price
We first pull summary statistics, and view as a table.
# Item price
summary(orders_full$item_price)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00 29.90 49.90 60.67 74.90 999.00Prices seem to fall within a reasonable range for apearal, with some extremes at both ends.
There are also two outliers with a price of 999. Appear can be expensive, so we will assume these are correct.
We can also see that 324 observations where the price is 0. All of these observations are unsized. There seem to be a mix of returned and unreturned items. This could indicate missing data, a promotion, or a return. Because these observations account for only .00432% of the dataset, we will exclude them:
orders_full <- orders_full %>%
filter(item_price > 0)Brand ID
Next, let’s look at brand. We have no information on what these brands are beyond a numeric ID. First, we’ll pull some summary information - I’d like to know how many unique brands there are, and if any of them account for a particularly large portion of sales.
# Check brand ID counts and how much of the data they make up
brand_counts <- orders_full %>%
count(brand_id, sort = TRUE)
n_unique_brands <- nrow(brand_counts)
n_unique_brands## [1] 129head(brand_counts, 25) # preview top 10 brands## brand_id n
## 1 3 8399
## 2 1 6546
## 3 5 4388
## 4 37 3168
## 5 11 2812
## 6 20 2205
## 7 17 1886
## 8 6 1875
## 9 38 1678
## 10 14 1631
## 11 40 1473
## 12 46 1466
## 13 43 1448
## 14 28 1403
## 15 49 1370
## 16 42 1360
## 17 29 1229
## 18 25 1082
## 19 33 982
## 20 21 976
## 21 19 962
## 22 30 962
## 23 47 931
## 24 48 900
## 25 53 884Here we can see that only three brands (IDs 3, 1, and 5) account for >5% of observations (computations were done outside of R, but can be easily reproduced). For simplicity and time, we will keep these, and merge all others into an “other” category.
We will also need to how brand data is stored. Currently, brand is stored as an integer. Because some models may multiply this by a coefficient, we will change brands to a string.
# Create buckets and set to string
orders_full <- orders_full %>%
mutate(
brand_id = as.character(brand_id), # ensure string type
brand_id = if_else(brand_id %in% c("1", "3", "5"),
brand_id, "Other")
)We can now see counts per brand over 2/3 of observations are “other”. That’s not perfect, but with our limited information, its what we will go for.
# See a table of the above
table(orders_full$brand_id)##
## 1 3 5 Other
## 6546 8399 4388 55339User State
Next, we will look at which Bundesland a given user is from.
# Show observation counts per state
unique(orders_full$user_state)## [1] "Rhineland-Palatinate" "Lower Saxony"
## [3] "Berlin" "Saxony"
## [5] "Hesse" "Schleswig-Holstein"
## [7] "North Rhine-Westphalia" "Baden-Wuerttemberg"
## [9] "Bavaria" "Mecklenburg-Western Pomerania"
## [11] "Saarland" "Bremen"
## [13] "Hamburg" "Brandenburg"
## [15] "Saxony-Anhalt" "Thuringia"table(orders_full$user_state)##
## Baden-Wuerttemberg Bavaria
## 9162 10443
## Berlin Brandenburg
## 2948 1676
## Bremen Hamburg
## 923 1934
## Hesse Lower Saxony
## 5506 11251
## Mecklenburg-Western Pomerania North Rhine-Westphalia
## 1159 16930
## Rhineland-Palatinate Saarland
## 3621 656
## Saxony Saxony-Anhalt
## 2335 868
## Schleswig-Holstein Thuringia
## 4164 1096This seems reasonable. All 16 federal subjects are represented, and distribution seems believable. Bremen having low total orders while Lower Saxony has high total orders seem odd, but we don’t know much about this company beyond this dataset. No cleanup action will be taken here.
Delivery Date/Time
Next, we will look at delivery date, order date, and the calculated variable delivery time.
Authors note: This was the first section I prepared, and served as my Feuertaufe for R. Please accept my apologies for any inefficiencies and redundancies in this section.
We will start by calculating delivery time, and creating a new variable.
#Calcualte delivery time
orders_full$delivery_time <- orders_full$delivery_date - orders_full$order_dateNext, we will ID possible issues by counting NA values. Note that we also did this for all observations in the initial data checks.
# Count NA values in observations
colSums(is.na(orders_full)) ## user_id order_item_id order_date delivery_date item_id
## 0 0 0 7900 0
## item_size item_color brand_id item_price return
## 0 0 0 0 0
## user_title user_dob user_state user_reg_date delivery_time
## 0 6351 0 0 7900We see ~6000 NA values in user_dob, and ~8000 in both delivery date and time.
On further investigation, we can see that NA values in delivery_date are only present when an order has been returned. This could be due to deletion of delivery data upon return, or a similar operational issue. We would like to avoid unnecessarily discarding these observations if possible.
# Conditional NA values
sum(is.na(orders_full$delivery_date[orders_full$return == "Yes"]))## [1] 0sum(is.na(orders_full$delivery_date[orders_full$return == "No"]))## [1] 7900In addition, we see that there are illogical delivery times in the 1990s. We will remove them, as these account for ~1% of the observations. Imputation could lead to inaccurate predictions as it may distort relationships between delivery date/time and other predictors.
# Summarize delivery date entries
summary(orders_full$delivery_date)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## "1994-12-31" "2016-07-04" "2016-07-26" "2016-04-13" "2016-08-11" "2017-01-23"
## NA's
## "7900"# Remove observations with impossible delivery dates
orders_full <- orders_full %>%
filter(delivery_date >= as.Date("2016-01-01") | is.na(delivery_date))We can create a histogram demonstrating that delivery dates are within a believable time-frame. Remember, all order dates are in June-August 2016.
# Create a year-month variable
orders_full$year_month <- floor_date(orders_full$delivery_date, unit = "month")
# update data for histogram to exclude the removed observations
monthly_counts_delivery_date <- orders_full %>%
group_by(year_month) %>%
summarize(count = n(), .groups = "drop")
#Plot histogram
ggplot(monthly_counts_delivery_date, aes(x = year_month, y = count)) +
geom_col(fill = "steelblue") +
scale_x_date(date_labels = "%Y-%m", date_breaks = "1 month") +
labs(title = "Number of Observations by Year-Month of Delivery",
x = "Year-Month",
y = "Number of Orders") +
theme_minimal() +
theme(axis.text.x = element_text(angle = 45, hjust = 1))## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_col()`).
We now create a similar visualization to better understand delivery time. For readability, we will only chart observations where delivery >100 days.
#
ggplot(orders_full, aes(x = delivery_time)) +
geom_histogram(binwidth = 1, fill = "steelblue", color = "black") +
labs(title = "Histogram of Delivery Time",
x = "Delivery Time (days)",
y = "Number of Observations") +
coord_cartesian(xlim = c(0, 50)) + # limit x-axis to 0-100
theme_minimal()## Don't know how to automatically pick scale for object of type <difftime>.
## Defaulting to continuous.## Warning: Removed 7900 rows containing non-finite outside the scale range
## (`stat_bin()`).
We can see that there the majority of orders are fulfilled within a week. Because high but non-outlier delivery times may disproportionately impact a model, we will bucket these. Because ~10% of our dataset has a value of NA for delivery time, we will also include an “unknown” bucket, as this category may have explanatory power. We choose roughly even buckets based on total delivery time.
# Create buckets
orders_full <- orders_full %>%
mutate(
delivery_days = as.numeric(delivery_time, units = "days"), # convert difftime → numeric days
delivery_time_bucket = case_when(
is.na(delivery_days) ~ "unknown",
delivery_days >= 0 & delivery_days <= 1 ~ "0-1 days",
delivery_days >= 1 & delivery_days <= 2 ~ "0-2 days",
delivery_days >= 2 & delivery_days <= 3 ~ "2-3 days",
delivery_days >= 3 & delivery_days <= 7 ~ "3-7 days",
delivery_days > 7 ~ ">7 days"
)
)
# Display table
table(orders_full$delivery_time_bucket)##
## >7 days 0-1 days 0-2 days 2-3 days 3-7 days unknown
## 12536 11570 13226 12161 16395 7900Order Date
We see that order date can be one of three months:
orders_month_count <- orders_full %>%
mutate(order_month = floor_date(order_date, "month")) # sets all dates to first day of month
# Plot histogram of orders per month
ggplot(orders_month_count, aes(x = order_month)) +
geom_histogram(stat = "count", fill = "steelblue", color = "black") +
labs(title = "Number of Orders per Month",
x = "Month",
y = "Number of Orders") +
scale_x_date(date_labels = "%b %Y", date_breaks = "1 month") +
theme_minimal() +
theme(axis.text.x = element_text(angle = 45, hjust = 1))## Warning in geom_histogram(stat = "count", fill = "steelblue", color = "black"):
## Ignoring unknown parameters: `binwidth`, `bins`, and `pad`
Seasonality can be significant, but this data is noisy. We will instead use three buckets for each month.
orders_full <- orders_full %>%
mutate(order_date_bucket = case_when(
format(order_date, "%Y-%m") == "2016-06" ~ "June",
format(order_date, "%Y-%m") == "2016-07" ~ "July",
format(order_date, "%Y-%m") == "2016-08" ~ "August",
TRUE ~ "Other"
))
# table
table(orders_full$order_date_bucket)##
## August July June
## 27855 25504 20429Account Age
We will look at an additional calculated variable: account age is the difference between order date and user registration date. First we calculate this value:
# Add account age
orders_full <- orders_full %>%
mutate(
# account age in days
account_age_days = as.numeric(order_date - user_reg_date)
)Numerous observations have a value of -1. This may be because the account was created for a new user. While ‘-1’ is not an intuitive value for for an account’s age, I believe this is realistic. Account creation systems may lag order placement and tracking systems for operational reasons.
As with delivery time, we will create buckets. An age of ‘-1’ will be represented as ‘New’:
# Bucket account age
orders_full <- orders_full %>%
mutate(
account_age_bucket = case_when(
account_age_days == -1 ~ "New",
account_age_days >= 0 & account_age_days <= 365 ~ "0-1 years",
account_age_days > 365 ~ "> 1 year",
TRUE ~ NA_character_ # catch-all for unexpected cases
)
)
table(orders_full$account_age_bucket)##
## > 1 year 0-1 years New
## 32158 27696 13934We see here a reasonably even split between buckets, with the exception ‘New’, which has fewer observations.
Dates of Birth
To get us started, we will use a histogram to visualize user birth dates:
#create counts for each unique string entry
counts_dob <- orders_full %>%
count(user_dob) %>%
arrange(desc(n))
#create chart
ggplot(counts_dob, aes(x = user_dob, y = n)) +
geom_col(fill = "steelblue", colour = "black", width = 0.8) +
labs(title = "Frequency of DOB",
x = "DOB",
y = "Count") +
theme_minimal() +
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1))## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_col()`).
Looks like we have some outliers, and some unrealistic birth dates. This may be because default dates are used when users are creating some accounts. There may be varying default dates on different platforms. There are also a significant number of NA values.
Let’s take a look at some of the outliers:
# Count DOB frequencies
dob_counts <- as.data.frame(table(orders_full$user_dob))
# Sort in decreasing order of frequency
dob_counts <- dob_counts[order(-dob_counts$Freq), ]
# Rename columns for clarity
colnames(dob_counts) <- c("user_dob", "count")
# Show the result
head(dob_counts)## user_dob count
## 1 1900-11-21 712
## 576 1949-11-21 344
## 4259 1965-06-10 182
## 7581 1980-12-11 107
## 869 1951-12-30 97
## 1278 1954-04-07 88It looks like some of the outliers occur pretty frequently. To prevent these from interfering with our predictive models, let’s convert them to NA. Specifically, we will convert all dates before 1930 and after 200 to NA. We also convert all values which occur in over 75 observations to NA - this is to remove default dates which could be realistic birth dates. We convert them to NA as we will bucket these values as ‘unknown’ at a later step.
Note on the histogram - we will convert DoB back to ymd at a later step. This histogram is only intended to visualize distribution.
# Convert DOBs that occur > 100 times to NA
orders_full <- orders_full %>%
add_count(user_dob, name = "dob_freq") %>%
mutate(user_dob = if_else(dob_freq > 75, NA_Date_, user_dob)) %>%
select(-dob_freq)
# Convert DOBs before 1930 or after 2000 to NA.
orders_full$user_dob <- ifelse(
year(orders_full$user_dob) < 1930 | year(orders_full$user_dob) > 2000,
NA,
orders_full$user_dob
)
#Re-plot the histogram to see practical results
#update counts for each unique string entry
counts_dob <- orders_full %>%
count(user_dob) %>%
arrange(desc(n))
#create chart
ggplot(counts_dob, aes(x = user_dob, y = n)) +
geom_col(fill = "steelblue", colour = "black", width = 0.8) +
labs(title = "Frequency of DOB",
x = "DOB",
y = "Count") +
theme_minimal() +
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1))## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_col()`).
While there are still some outliers here, this looks like a much more realistic distribution. As before, we will create buckets. Because we don’t know much about user behavior outside this dataset, we will create 5 buckets: 4 quartiles representing the highest 25% of ages, the next highest 25%, etc, and unknown:
# 1. Compute quartiles of user_dob (ignoring NAs)
dob_quartiles <- quantile(orders_full$user_dob, probs = c(0.25, 0.5, 0.75), na.rm = TRUE)
# 2. Create DOB bins based on quartiles, plus "Unknown"
orders_full <- orders_full %>%
mutate(dob_bin = case_when(
is.na(user_dob) ~ "Unknown",
user_dob <= dob_quartiles[1] ~ "Q1", # oldest
user_dob <= dob_quartiles[2] ~ "Q2",
user_dob <= dob_quartiles[3] ~ "Q3",
user_dob > dob_quartiles[3] ~ "Q4"
))
# Make it a factor in quartile order
orders_full$dob_bin <- factor(orders_full$dob_bin,
levels = c("Q1", "Q2", "Q3", "Q4", "Unknown"))
# Convert dates back to ymd
orders_full$user_dob <- as.Date(orders_full$user_dob, origin = "1970-01-01")User Title
We take a quick look at the data on user title:
# Summarize possible values and counts for user_title
unique(orders_full$user_title)## [1] "Mrs" "Company" "Mr" "Family" "not reported"table(orders_full$user_title)##
## Company Family Mr Mrs not reported
## 90 297 2755 70549 97The outputs here seem reasonable. While some bins are small, they appear distint, and intuitively may have predictive power. We move on without taking action.
Item Color
Let’s move on to item color. It is odd that this online retailer does not have some identifier for item category (pants, shirts, jackets, etc), beyond just color, but we will work with the data we have available.
# Summarize possible values and counts for item_color
unique(orders_full$item_color)## [1] "black" "ancient" "grey" "turquoise"
## [5] "green" "ash" "red" "denim"
## [9] "petrol" "brown" "ocher" "berry"
## [13] "aquamarine" "mahagoni" "blue" "white"
## [17] "khaki" "silver" "purple" "dark oliv"
## [21] "pink" "anthracite" "aubergine" "orange"
## [25] "navy" "mint" "olive" "curry"
## [29] "stained" "mocca" "magenta" "darkblue"
## [33] "basalt" "?" "nature" "pallid"
## [37] "bordeaux" "fuchsia" "dark denim" "blau"
## [41] "azure" "hibiscus" "beige" "striped"
## [45] "yellow" "curled" "ecru" "dark garnet"
## [49] "dark navy" "floral" "gold" "aqua"
## [53] "champagner" "jade" "coral" "cognac"
## [57] "dark grey" "ebony" "baltic blue" "cobalt blue"
## [61] "antique pink" "creme" "habana" "mango"
## [65] "currant purple" "kanel" "almond" "terracotta"
## [69] "ivory" "brwon" "apricot" "caramel"
## [73] "aviator" "avocado" "amethyst" "copper coin"
## [77] "opal"table(orders_full$item_color)##
## ? almond amethyst ancient anthracite
## 25 13 5 121 1723
## antique pink apricot aqua aquamarine ash
## 10 10 35 254 1191
## aubergine aviator avocado azure baltic blue
## 699 8 4 397 19
## basalt beige berry black blau
## 189 609 1465 13297 178
## blue bordeaux brown brwon caramel
## 8228 481 6328 10 9
## champagner cobalt blue cognac copper coin coral
## 42 25 58 6 133
## creme curled currant purple curry dark denim
## 3 11 6 244 561
## dark garnet dark grey dark navy dark oliv darkblue
## 11 32 21 30 59
## denim ebony ecru floral fuchsia
## 1543 11 230 57 92
## gold green grey habana hibiscus
## 37 3814 6623 14 157
## ivory jade kanel khaki magenta
## 37 46 31 338 537
## mahagoni mango mint mocca nature
## 92 62 125 1366 649
## navy ocher olive opal orange
## 236 1859 1616 1 723
## pallid petrol pink purple red
## 229 2495 788 2936 5406
## silver stained striped terracotta turquoise
## 102 1227 83 7 727
## white yellow
## 2693 249This seems noisy. We will manually combine some obvious pairs, based on intuition. This was done manually, by using intuition to check and modify groups suggested by an LLM.
It is worth noting that an approach like factor analysis may be applied to similar problems where there is too much noise/too many variables. Because colors are relatively intuitive, we took a simple approach in this case.
orders_full <- orders_full %>%
mutate(item_color_simple = case_when(
item_color %in% c("black", "ebony", "anthracite", "basalt") ~ "black",
item_color %in% c("grey", "ash", "dark grey") ~ "grey",
item_color %in% c("white", "ecru", "creme", "ivory") ~ "white",
item_color %in% c("red", "dark garnet", "bordeaux", "terracotta") ~ "red",
item_color %in% c("blue", "navy", "darkblue", "blau", "denim", "dark denim", "cobalt blue", "baltic blue", "azure") ~ "blue",
item_color %in% c("green", "mint", "olive", "avocado", "petrol", "turquoise", "aqua", "aquamarine") ~ "green",
item_color %in% c("yellow", "gold", "mango", "apricot", "curry", "champagner", "ocher", "orange") ~ "yellow",
item_color %in% c("pink", "magenta", "hibiscus", "antique pink", "fuchsia", "purple", "aubergine", "currant purple") ~ "magenta",
item_color %in% c("brown", "brwon", "mocca", "cognac", "caramel", "mahagoni", "terracotta") ~ "brown",
TRUE ~ "other"
))Item Size
Item size presents us with a challenge. There is no unified system for sizing, and sizes are represented numerically, semi-numerically (46+, etc), and non-numerically. Within the numerical sizes, there is no obvious sizing schema.
Let’s take a quick look:
# Summarize possible values and counts for item_size
unique(orders_full$item_size)## [1] "40" "37" "41" "l" "m" "42" "xl"
## [8] "unsized" "33" "44" "11" "39" "152" "xxl"
## [15] "38" "41+" "43" "28" "34" "1" "13"
## [22] "30" "32" "2" "35" "3" "29" "31"
## [29] "27" "10" "36+" "36" "39+" "26" "38+"
## [36] "s" "18" "7" "7+" "42+" "24" "21"
## [43] "22" "46" "3332" "8" "9" "5+" "5"
## [50] "4+" "8+" "40+" "4" "48" "6+" "10+"
## [57] "37+" "104" "6" "50" "45" "23" "12"
## [64] "116" "128" "14" "47" "xs" "9+" "164"
## [71] "25" "140" "20" "3634" "xxxl" "4032" "45+"
## [78] "19" "3632" "44+" "80" "46+" "54" "2+"
## [85] "100" "56" "3+" "52" "3832" "3132" "3432"
## [92] "58" "11+" "43+" "176" "90" "49" "4034"table(orders_full$item_size)##
## 1 10 10+ 100 104 11 11+ 116 12 128
## 28 192 136 11 44 137 15 70 56 87
## 13 14 140 152 164 176 18 19 2 2+
## 28 6 92 77 21 6 73 78 6 2
## 20 21 22 23 24 25 26 27 28 29
## 84 111 99 77 86 50 79 82 62 117
## 3 3+ 30 31 3132 32 33 3332 34 3432
## 23 10 101 143 1 121 106 24 132 8
## 35 36 36+ 3632 3634 37 37+ 38 38+ 3832
## 176 1039 18 15 1 1719 55 4217 110 7
## 39 39+ 4 4+ 40 40+ 4032 4034 41 41+
## 4926 151 118 76 5591 122 9 1 3642 93
## 42 42+ 43 43+ 44 44+ 45 45+ 46 46+
## 3666 92 1200 8 1493 8 165 7 934 7
## 47 48 49 5 5+ 50 52 54 56 58
## 31 663 1 190 196 116 18 18 15 2
## 6 6+ 7 7+ 8 8+ 80 9 9+ 90
## 342 326 308 264 382 184 2 297 149 7
## l m s unsized xl xs xxl xxxl
## 9221 7502 2841 3257 8289 28 6723 69# Since the table command is difficult to intuitively read, we create a separate dataset to interpret frequencies of sizes
item_size_freq <- orders_full %>%
group_by(item_size) %>%
summarize(count = n(), .groups = "drop") %>%
arrange(desc(count))
item_size_freq$percent <- item_size_freq$count / sum(item_size_freq$count)
item_size_freq$percent <- paste0(round(item_size_freq$percent * 100, 1), "%")
head(item_size_freq)## # A tibble: 6 × 3
## item_size count percent
## <chr> <int> <chr>
## 1 l 9221 12.5%
## 2 xl 8289 11.2%
## 3 m 7502 10.2%
## 4 xxl 6723 9.1%
## 5 40 5591 7.6%
## 6 39 4926 6.7%From this, we can identify several patterns:
Numerous observations fit into expected “s, m, l, xl…” categories
A meaningful subset of sizes have a “+” symbol next to them - we will make an assumption and consider all of these “plus sizes”
The majority of non-‘+’ numeric sizes appear to have a numeric value <50.
While they may follow different sizing schemas across observations, the numeric values may nonetheless have explanatory power.
We will take several actions here:
We will merge the letter sizes into: S, M, L, and XL, and XXL
We remove all 536 observations in which numeric size is >50
We merge all remaining numeric sizes into quartiles.
Plus non-numeric, and we will merge these into a separate ‘plus’ category
# Merge the letter sizes into: S, M, L, and XL, and XXL:
orders_full <- orders_full %>%
mutate(
item_size = dplyr::recode(item_size,
"XS" = "S",
"XXXL" = "XXL")
)%>%
# remove observations in which numeric size is >50
mutate(size_numeric = suppressWarnings(as.numeric(item_size))) %>%
filter(is.na(size_numeric) | size_numeric <= 50)
# Merge numeric into quartiles
orders_full <- orders_full %>%
# Standardize letter sizes and trim whitespace
mutate(item_size = toupper(trimws(item_size))) %>%
# Compute quartiles for numeric sizes only
mutate(size_quartiles = ifelse(!is.na(size_numeric),
ntile(size_numeric, 4),
NA)) %>%
# Create unified clean size column
mutate(size_category = case_when(
item_size %in% c("S","M","L","XL","XXL") ~ item_size, # keep standard letters
!is.na(size_numeric) & size_quartiles == 1 ~ "numeric Q1", # numeric quartiles
!is.na(size_numeric) & size_quartiles == 2 ~ "numeric Q2",
!is.na(size_numeric) & size_quartiles == 3 ~ "numeric Q3",
!is.na(size_numeric) & size_quartiles == 4 ~ "numeric Q4",
TRUE ~ "plus" # everything else → plus
)) %>%
# convert to factor with sensible order
mutate(size_category = factor(size_category,
levels = c("S","M","L","XL","XXL",
"numeric Q1","numeric Q2","numeric Q3","numeric Q4",
"plus"))) %>%
select(-size_numeric, -size_quartiles) # remove helper columns
# Let's take a look at the new Size category variable
item_size_category_freq <- orders_full %>%
group_by(size_category) %>%
summarize(count = n(), .groups = "drop") %>%
arrange(desc(count))
item_size_category_freq$percent <- item_size_category_freq$count / sum(item_size_category_freq$count)
item_size_category_freq$percent <- paste0(round(item_size_category_freq$percent * 100, 1), "%")
table(orders_full$size_category) # optional to view above##
## S M L XL XXL numeric Q1 numeric Q2
## 2841 7502 9221 8289 6723 8324 8323
## numeric Q3 numeric Q4 plus
## 8323 8323 5383Last steps
Now that we have completed data cleaning, we will save a CSV of our cleaned data locally. In case you encounter errors if following along with this markup, I have attached the unmodified R Script I use to create orders_clean.csv. Code included as text to enable R markup knitting.
#write.csv(orders_full, "orders_clean.csv", row.names = FALSE)Classification models
Now that our data is cleaned, we will try to fit a number of classification models, with the goal of predicting whether an order will be returned. Results are mixed - the lowest error rate we are able to achieve is: 0.3582545.
Logistic regression
We will start with a logistic regression. Before we get started, let’s double check that we have the right packages:
library(lubridate) # date management
library(ggplot2) # charts
library(dplyr) # data handling
library(glmnet) # implementation of the methods
library(caret) # machine-learning workflow
library(car) # multicolinearity check## Warning: package 'car' was built under R version 4.4.3## Loading required package: carData## Warning: package 'carData' was built under R version 4.4.3##
## Attaching package: 'car'## The following object is masked from 'package:dplyr':
##
## recodeNext, we will set directory, and load our clean data that we saved in the last step:
# set working directory:
setwd("C:/Users/jedel/OneDrive/Desktop/Data for ML class/Individual Assignment")
#Load data
orders_clean <- read.csv("orders_clean.csv", sep = ",")Next, we will prepare our training and test data. Seed is set to 8888 for reporduceability. We will do 70/30 split:
# Split the sample into training and test set. .7/1 is used as training data.
# Set seed
set.seed(8888) # for reproducibility. I use seed "8888"
# Determine sample size
train_size <- floor(0.7 * nrow(orders_clean))
# Generate random row indices for training set
train_indices <- sample(seq_len(nrow(orders_clean)), size = train_size)
# Split the dataset
orders_train <- orders_clean[train_indices, ]
orders_test <- orders_clean[-train_indices, ]Next, we start fitting our logistic regression:
# First, Double-check "return" is coded as a factor
orders_train$return <- factor(orders_train$return, levels = c("No", "Yes"))
# Fit logistic regression
base_model <- glm(return ~
size_category + # cleaned size
item_color_simple + # simplified color
brand_id + # brand, note that there are now 4 possible options
item_price + # price (cleaned)
user_title + # title
user_state + # state
dob_bin + # bucket of age ranges + unknown
delivery_time_bucket + # bucket of delivery times + unknown
order_date_bucket + # order date bucket (month)
account_age_bucket, # account age bucket
data = orders_train,
family = binomial)
summary(base_model)##
## Call:
## glm(formula = return ~ size_category + item_color_simple + brand_id +
## item_price + user_title + user_state + dob_bin + delivery_time_bucket +
## order_date_bucket + account_age_bucket, family = binomial,
## data = orders_train)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.297e-01 2.921e-01 0.444 0.657041
## size_categoryM 3.713e-02 3.936e-02 0.943 0.345460
## size_categorynumeric Q1 -6.600e-02 4.054e-02 -1.628 0.103493
## size_categorynumeric Q2 -1.147e-01 4.148e-02 -2.765 0.005700
## size_categorynumeric Q3 -6.141e-02 4.126e-02 -1.488 0.136641
## size_categorynumeric Q4 3.473e-02 4.068e-02 0.854 0.393278
## size_categoryplus -3.964e-01 4.574e-02 -8.666 < 2e-16
## size_categoryS -1.034e-02 5.461e-02 -0.189 0.849821
## size_categoryXL 2.534e-02 3.891e-02 0.651 0.514926
## size_categoryXXL -5.643e-02 4.130e-02 -1.366 0.171798
## item_color_simpleblue 1.819e-01 3.266e-02 5.568 2.57e-08
## item_color_simplebrown 1.188e-01 3.578e-02 3.321 0.000897
## item_color_simplegreen 1.397e-01 3.476e-02 4.019 5.84e-05
## item_color_simplegrey 1.547e-01 3.629e-02 4.264 2.01e-05
## item_color_simplemagenta 8.146e-02 4.199e-02 1.940 0.052371
## item_color_simpleother 3.777e-02 4.216e-02 0.896 0.370306
## item_color_simplered 1.158e-01 4.023e-02 2.878 0.004007
## item_color_simplewhite 1.902e-01 5.333e-02 3.565 0.000363
## item_color_simpleyellow 1.110e-01 5.141e-02 2.159 0.030877
## brand_id3 -2.149e-01 5.172e-02 -4.156 3.24e-05
## brand_id5 -6.065e-02 5.733e-02 -1.058 0.290172
## brand_idOther 1.517e-01 3.672e-02 4.133 3.59e-05
## item_price 6.471e-03 2.714e-04 23.847 < 2e-16
## user_titleFamily -1.034e+00 3.270e-01 -3.161 0.001571
## user_titleMr -1.069e+00 2.874e-01 -3.718 0.000201
## user_titleMrs -9.065e-01 2.832e-01 -3.201 0.001368
## user_titlenot reported -1.118e+00 3.924e-01 -2.849 0.004388
## user_stateBavaria 1.371e-02 3.721e-02 0.368 0.712612
## user_stateBerlin -1.570e-02 5.537e-02 -0.283 0.776806
## user_stateBrandenburg -6.652e-03 6.848e-02 -0.097 0.922625
## user_stateBremen -2.365e-01 8.937e-02 -2.647 0.008126
## user_stateHamburg -5.719e-02 6.547e-02 -0.874 0.382337
## user_stateHesse -3.211e-03 4.408e-02 -0.073 0.941927
## user_stateLower Saxony 1.535e-02 3.669e-02 0.419 0.675577
## user_stateMecklenburg-Western Pomerania 1.201e-01 8.211e-02 1.463 0.143416
## user_stateNorth Rhine-Westphalia -5.863e-02 3.374e-02 -1.738 0.082248
## user_stateRhineland-Palatinate -5.424e-02 5.106e-02 -1.062 0.288103
## user_stateSaarland 1.677e-01 1.046e-01 1.604 0.108726
## user_stateSaxony 4.150e-03 5.979e-02 0.069 0.944660
## user_stateSaxony-Anhalt -6.938e-03 9.046e-02 -0.077 0.938869
## user_stateSchleswig-Holstein 4.117e-02 4.859e-02 0.847 0.396821
## user_stateThuringia -1.707e-01 8.297e-02 -2.057 0.039681
## dob_binQ2 8.375e-02 2.856e-02 2.933 0.003359
## dob_binQ3 1.679e-01 2.852e-02 5.885 3.97e-09
## dob_binQ4 3.028e-01 2.912e-02 10.399 < 2e-16
## dob_binUnknown 9.276e-02 3.535e-02 2.624 0.008699
## delivery_time_bucket0-1 days 2.492e-01 3.214e-02 7.753 8.95e-15
## delivery_time_bucket0-2 days 2.039e-01 3.075e-02 6.630 3.36e-11
## delivery_time_bucket2-3 days 3.303e-01 3.149e-02 10.490 < 2e-16
## delivery_time_bucket3-7 days 3.030e-01 2.953e-02 10.264 < 2e-16
## delivery_time_bucketunknown -1.741e+01 5.254e+01 -0.331 0.740344
## order_date_bucketJuly 4.247e-02 2.317e-02 1.833 0.066771
## order_date_bucketJune 4.160e-02 2.603e-02 1.598 0.110046
## account_age_bucket0-1 years -1.743e-01 2.135e-02 -8.161 3.33e-16
## account_age_bucketNew -7.269e-02 2.745e-02 -2.648 0.008094
##
## (Intercept)
## size_categoryM
## size_categorynumeric Q1
## size_categorynumeric Q2 **
## size_categorynumeric Q3
## size_categorynumeric Q4
## size_categoryplus ***
## size_categoryS
## size_categoryXL
## size_categoryXXL
## item_color_simpleblue ***
## item_color_simplebrown ***
## item_color_simplegreen ***
## item_color_simplegrey ***
## item_color_simplemagenta .
## item_color_simpleother
## item_color_simplered **
## item_color_simplewhite ***
## item_color_simpleyellow *
## brand_id3 ***
## brand_id5
## brand_idOther ***
## item_price ***
## user_titleFamily **
## user_titleMr ***
## user_titleMrs **
## user_titlenot reported **
## user_stateBavaria
## user_stateBerlin
## user_stateBrandenburg
## user_stateBremen **
## user_stateHamburg
## user_stateHesse
## user_stateLower Saxony
## user_stateMecklenburg-Western Pomerania
## user_stateNorth Rhine-Westphalia .
## user_stateRhineland-Palatinate
## user_stateSaarland
## user_stateSaxony
## user_stateSaxony-Anhalt
## user_stateSchleswig-Holstein
## user_stateThuringia *
## dob_binQ2 **
## dob_binQ3 ***
## dob_binQ4 ***
## dob_binUnknown **
## delivery_time_bucket0-1 days ***
## delivery_time_bucket0-2 days ***
## delivery_time_bucket2-3 days ***
## delivery_time_bucket3-7 days ***
## delivery_time_bucketunknown
## order_date_bucketJuly .
## order_date_bucketJune
## account_age_bucket0-1 years ***
## account_age_bucketNew **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 70675 on 51275 degrees of freedom
## Residual deviance: 61978 on 51221 degrees of freedom
## AIC: 62088
##
## Number of Fisher Scoring iterations: 16The summary data gives us insight into the predictive power of different variables, as well an overview of the regression function. We can see that order_date_bucket and user_state are not very predictive, and that only some states and size categories are predictive.
Before we move on to our next steps, let’s also check multicolinearity:
# Check multidisciplinary. 1=non, 5=moderate, 10=high
vif(base_model)## GVIF Df GVIF^(1/(2*Df))
## size_category 1.586985 9 1.025990
## item_color_simple 1.188109 9 1.009622
## brand_id 1.661080 3 1.088258
## item_price 1.442735 1 1.201139
## user_title 1.052868 4 1.006461
## user_state 1.076232 15 1.002452
## dob_bin 1.088840 4 1.010696
## delivery_time_bucket 1.098420 5 1.009432
## order_date_bucket 1.215751 2 1.050053
## account_age_bucket 1.107194 2 1.025784Everything looks good on multicolinearity - let’s test our model’s accuracy, and see if we can find any improvements.
#Run logistic regression on test data
# Make sure "return" is a factor with the same levels as in training
orders_test$return <- factor(orders_test$return, levels = levels(orders_train$return))
# Make predictions (predicted probabilities that "return" = "Yes")
# Predicted probabilities
pred_probs <- predict(base_model, newdata = orders_test, type = "response")
#Convert probabilities to class preditctions (.5 is cuttoff)
pred_class <- ifelse(pred_probs > 0.5, "Yes", "No")
pred_class <- factor(pred_class, levels = levels(orders_test$return))
#Compare the error rate
# Compare predicted vs actual
error_rate <- mean(pred_class != orders_test$return)
error_rate## [1] 0.3814161Our base model is getting an error rate of 0.3813161. Let’s try to improve this. Based on the summary and error rate, we will test a model with fewer variables, removing order date, size, and state:
# Fit new logistic regression
reduced_model <- glm(return ~
item_color_simple + # simplified color
brand_id + # brand, note that there are now 4 possible options
item_price + # price (cleaned)
user_title + # title
dob_bin + # bucket of age ranges + unknown
delivery_time_bucket + # bucket of delivery times + unknown
account_age_bucket, # account age bucket
data = orders_train,
family = binomial)
summary(reduced_model)##
## Call:
## glm(formula = return ~ item_color_simple + brand_id + item_price +
## user_title + dob_bin + delivery_time_bucket + account_age_bucket,
## family = binomial, data = orders_train)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 7.761e-02 2.876e-01 0.270 0.787296
## item_color_simpleblue 2.038e-01 3.253e-02 6.264 3.75e-10 ***
## item_color_simplebrown 1.143e-01 3.567e-02 3.204 0.001354 **
## item_color_simplegreen 1.328e-01 3.452e-02 3.846 0.000120 ***
## item_color_simplegrey 1.526e-01 3.611e-02 4.226 2.38e-05 ***
## item_color_simplemagenta 8.743e-02 4.188e-02 2.088 0.036827 *
## item_color_simpleother 1.256e-02 4.179e-02 0.301 0.763793
## item_color_simplered 1.183e-01 4.005e-02 2.953 0.003151 **
## item_color_simplewhite 2.244e-01 5.290e-02 4.242 2.22e-05 ***
## item_color_simpleyellow 1.255e-01 5.112e-02 2.455 0.014069 *
## brand_id3 -1.303e-01 4.678e-02 -2.785 0.005352 **
## brand_id5 2.367e-02 5.341e-02 0.443 0.657721
## brand_idOther 1.812e-01 3.530e-02 5.134 2.84e-07 ***
## item_price 6.463e-03 2.534e-04 25.511 < 2e-16 ***
## user_titleFamily -1.112e+00 3.256e-01 -3.415 0.000637 ***
## user_titleMr -1.103e+00 2.868e-01 -3.847 0.000119 ***
## user_titleMrs -9.335e-01 2.825e-01 -3.305 0.000951 ***
## user_titlenot reported -1.137e+00 3.906e-01 -2.910 0.003617 **
## dob_binQ2 8.614e-02 2.846e-02 3.027 0.002469 **
## dob_binQ3 1.706e-01 2.841e-02 6.006 1.90e-09 ***
## dob_binQ4 3.005e-01 2.877e-02 10.444 < 2e-16 ***
## dob_binUnknown 9.483e-02 3.523e-02 2.691 0.007114 **
## delivery_time_bucket0-1 days 2.515e-01 3.199e-02 7.863 3.75e-15 ***
## delivery_time_bucket0-2 days 2.032e-01 3.061e-02 6.639 3.17e-11 ***
## delivery_time_bucket2-3 days 3.275e-01 3.127e-02 10.472 < 2e-16 ***
## delivery_time_bucket3-7 days 3.045e-01 2.930e-02 10.394 < 2e-16 ***
## delivery_time_bucketunknown -1.743e+01 5.266e+01 -0.331 0.740567
## account_age_bucket0-1 years -1.721e-01 2.125e-02 -8.100 5.50e-16 ***
## account_age_bucketNew -7.175e-02 2.734e-02 -2.625 0.008668 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 70675 on 51275 degrees of freedom
## Residual deviance: 62136 on 51247 degrees of freedom
## AIC: 62194
##
## Number of Fisher Scoring iterations: 16# Make sure "return" is a factor with the same levels as in training
orders_test$return <- factor(orders_test$return, levels = levels(orders_train$return))
# Make predictions (predicted probabilities that "return" = "Yes")
# Predicted probabilities
pred_probs <- predict(reduced_model, newdata = orders_test, type = "response")
#Convert probabilities to class preditctions (.5 is cuttoff)
pred_class <- ifelse(pred_probs > 0.5, "Yes", "No")
pred_class <- factor(pred_class, levels = levels(orders_test$return))
#Compare the error rate
# Compare predicted vs actual
error_rate <- mean(pred_class != orders_test$return)
error_rate## [1] 0.3833728We see that this has a marginally higher error rate than before. We will try adding size back:
refined_model <- glm(return ~
size_category + # cleaned size
item_color_simple + # simplified color
brand_id + # brand, note that there are now 4 possible options
item_price + # price (cleaned)
user_title + # title
dob_bin + # bucket of age ranges + unknown
delivery_time_bucket + # bucket of delivery times + unknown
account_age_bucket, # account age bucket
data = orders_train,
family = binomial)
summary(refined_model)##
## Call:
## glm(formula = return ~ size_category + item_color_simple + brand_id +
## item_price + user_title + dob_bin + delivery_time_bucket +
## account_age_bucket, family = binomial, data = orders_train)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.600e-01 2.894e-01 0.553 0.580460
## size_categoryM 3.790e-02 3.932e-02 0.964 0.335078
## size_categorynumeric Q1 -5.850e-02 4.043e-02 -1.447 0.147881
## size_categorynumeric Q2 -1.094e-01 4.134e-02 -2.646 0.008137 **
## size_categorynumeric Q3 -5.554e-02 4.111e-02 -1.351 0.176735
## size_categorynumeric Q4 3.633e-02 4.064e-02 0.894 0.371301
## size_categoryplus -3.983e-01 4.564e-02 -8.727 < 2e-16 ***
## size_categoryS -9.631e-03 5.454e-02 -0.177 0.859849
## size_categoryXL 2.682e-02 3.884e-02 0.690 0.489944
## size_categoryXXL -5.701e-02 4.124e-02 -1.383 0.166782
## item_color_simpleblue 1.824e-01 3.263e-02 5.591 2.25e-08 ***
## item_color_simplebrown 1.206e-01 3.576e-02 3.372 0.000745 ***
## item_color_simplegreen 1.362e-01 3.470e-02 3.925 8.69e-05 ***
## item_color_simplegrey 1.588e-01 3.621e-02 4.386 1.16e-05 ***
## item_color_simplemagenta 8.057e-02 4.196e-02 1.920 0.054822 .
## item_color_simpleother 3.972e-02 4.212e-02 0.943 0.345649
## item_color_simplered 1.178e-01 4.020e-02 2.930 0.003387 **
## item_color_simplewhite 1.972e-01 5.315e-02 3.711 0.000207 ***
## item_color_simpleyellow 1.163e-01 5.133e-02 2.265 0.023501 *
## brand_id3 -2.292e-01 5.125e-02 -4.472 7.73e-06 ***
## brand_id5 -5.943e-02 5.720e-02 -1.039 0.298827
## brand_idOther 1.420e-01 3.650e-02 3.890 0.000100 ***
## item_price 6.299e-03 2.569e-04 24.517 < 2e-16 ***
## user_titleFamily -1.072e+00 3.260e-01 -3.290 0.001001 **
## user_titleMr -1.075e+00 2.869e-01 -3.745 0.000180 ***
## user_titleMrs -9.081e-01 2.827e-01 -3.212 0.001316 **
## user_titlenot reported -1.087e+00 3.917e-01 -2.776 0.005510 **
## dob_binQ2 8.567e-02 2.850e-02 3.006 0.002649 **
## dob_binQ3 1.707e-01 2.846e-02 5.998 2.00e-09 ***
## dob_binQ4 3.050e-01 2.891e-02 10.548 < 2e-16 ***
## dob_binUnknown 9.302e-02 3.529e-02 2.636 0.008400 **
## delivery_time_bucket0-1 days 2.480e-01 3.206e-02 7.736 1.03e-14 ***
## delivery_time_bucket0-2 days 2.004e-01 3.067e-02 6.534 6.40e-11 ***
## delivery_time_bucket2-3 days 3.265e-01 3.134e-02 10.416 < 2e-16 ***
## delivery_time_bucket3-7 days 3.006e-01 2.936e-02 10.239 < 2e-16 ***
## delivery_time_bucketunknown -1.740e+01 5.256e+01 -0.331 0.740599
## account_age_bucket0-1 years -1.747e-01 2.129e-02 -8.205 2.30e-16 ***
## account_age_bucketNew -7.424e-02 2.741e-02 -2.708 0.006762 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 70675 on 51275 degrees of freedom
## Residual deviance: 62011 on 51238 degrees of freedom
## AIC: 62087
##
## Number of Fisher Scoring iterations: 16#Run logistic regression on test data
# Make sure "return" is a factor with the same levels as in training
orders_test$return <- factor(orders_test$return, levels = levels(orders_train$return))
# Make predictions (predicted probabilities that "return" = "Yes")
# Predicted probabilities
pred_probs <- predict(refined_model, newdata = orders_test, type = "response")
#Convert probabilities to class preditctions (.5 is cuttoff)
pred_class <- ifelse(pred_probs > 0.5, "Yes", "No")
pred_class <- factor(pred_class, levels = levels(orders_test$return))
#Compare the error rate
# Compare predicted vs actual
error_rate <- mean(pred_class != orders_test$return)
error_rate## [1] 0.3802785We see that this is marginally higher success rate than the base model. We can visualize this in a confusion matrix:
#Visualize above via a confusion matrix
cm <- confusionMatrix(pred_class, orders_test$return)
cm_df <- as.data.frame(cm$table)
colnames(cm_df) <- c("Predicted", "Actual", "Count")
cm_df$Type <- with(cm_df, ifelse(
(Predicted == "Yes" & Actual == "Yes") | (Predicted == "No" & Actual == "No"),
"Correct", "Incorrect"
))
ggplot(cm_df, aes(x = Actual, y = Predicted, fill = Type)) +
geom_tile(color = "white", size = 0.5) +
geom_text(aes(label = Count), size = 6) +
scale_fill_manual(values = c("Correct" = "lightgreen", "Incorrect" = "lightcoral")) +
labs(title = "Confusion Matrix",
x = "Actual Class",
y = "Predicted Class") +
theme_minimal(base_size = 14) +
theme(legend.position = "none")## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
We’ve tried a number of of optimization and have only reached ~62% accuracy. Let’s try some alternative models.
Classification Tree
Next, we will try running a classification tree. This tree splits data on predictors to predict whether and order will be returned or not.
As always, we’ll start off by loading our packages and data:
# Packages
library(rpart) # for building decision trees## Warning: package 'rpart' was built under R version 4.4.3library(rpart.plot)# for visualizing trees## Warning: package 'rpart.plot' was built under R version 4.4.3library(lubridate) # date management
library(ggplot2) # charts
library(dplyr) # data handling
library(glmnet) # implementation of the methods
library(caret) # machine-learning workflow
library(car) # multicolinearity check
# set working directory:
setwd("C:/Users/jedel/OneDrive/Desktop/Data for ML class/Individual Assignment")
#Load data
orders_clean <- read.csv("orders_clean.csv", sep = ",")Next, we’ll split our data into training and test sets:
# Split the sample into training and test set. .7/1 is used as training data.
# Set seed
set.seed(8888) # for reproducibility. I use seed "8888"
# Determine sample size
train_size <- floor(0.7 * nrow(orders_clean))
# Generate random row indices for training set
train_indices <- sample(seq_len(nrow(orders_clean)), size = train_size)
# Split the dataset
orders_train <- orders_clean[train_indices, ]
orders_test <- orders_clean[-train_indices, ]Next, we’ll user user our data to build our baseline tree
# Build a classification tree
tree_model <- rpart(
return ~
size_category + # cleaned size
item_color_simple + # simplified color
brand_id + # brand, note that there are now 4 possible options
item_price + # price (cleaned)
user_title + # title
dob_bin + # bucket of age ranges + unknown
delivery_time_bucket + # bucket of delivery times + unknown
account_age_bucket, # account age bucket
data = orders_train,
method = "class", # classification tree
control = rpart.control(cp = 0.01, minsplit = 20) # tuning parameters
)
# Visualize the tree
rpart.plot(tree_model, type = 2, extra = 104, fallen.leaves = TRUE)
I have included a visualization of our first tree above. Before testing our tree on the test data, let’s take a look at variable importance: It looks like item color does not provide much explanatory power.
#Variable importance
tree_model$variable.importance## delivery_time_bucket item_price brand_id
## 2583.266571 593.999012 105.171679
## size_category item_color_simple
## 10.177667 2.008834Let’s take a look at how the model performs:
# Predict class labels on the training set
pred_class <- predict(tree_model, newdata = orders_train, type = "class")
# Actual values
actual <- orders_train$return
# Create a confusion table
conf_matrix <- table(Predicted = pred_class, Actual = actual)
conf_matrix## Actual
## Predicted No Yes
## No 13079 4612
## Yes 14845 18740error_rate <- sum(pred_class != actual) / length(actual)
error_rate## [1] 0.3794563Error rate is not meaningfully below our results from logistic regression. Let’s try optimizing some the complexity parameter, which controls how many branches and hence how precise (to the training data) the tree can be. Below, we test a grid of possible values to find an optimal complexity parameter:
#optimize model
# Define parameter grid
grid <- expand.grid(cp = seq(0.0001, 0.05, length.out = 20))
grid## cp
## 1 0.000100000
## 2 0.002726316
## 3 0.005352632
## 4 0.007978947
## 5 0.010605263
## 6 0.013231579
## 7 0.015857895
## 8 0.018484211
## 9 0.021110526
## 10 0.023736842
## 11 0.026363158
## 12 0.028989474
## 13 0.031615789
## 14 0.034242105
## 15 0.036868421
## 16 0.039494737
## 17 0.042121053
## 18 0.044747368
## 19 0.047373684
## 20 0.050000000# 5-fold cross validation
train_control <- trainControl(method = "cv", number = 5)
# Train a tuned tree
tree_tuned <- train(
return ~ size_category + item_color_simple + brand_id + item_price +
user_title + user_state + dob_bin + delivery_time_bucket +
order_date_bucket + account_age_bucket,
data = orders_train,
method = "rpart",
trControl = train_control,
tuneGrid = grid
)
#Show best CP. This step takes some time
tree_tuned$bestTune # shows best cp## cp
## 2 0.002726316# We see an optimal cp of 0.002726316We can see that the optimal CP is ~0.0027. Let’s plug that into our model and see what kind of error rate we get:
# We now use this CP for our optimized tree
final_tree <- rpart(
return ~ size_category + item_color_simple + brand_id + item_price +
user_title + user_state + dob_bin + delivery_time_bucket +
order_date_bucket + account_age_bucket,
data = orders_train,
method = "class",
control = rpart.control(cp = tree_tuned$bestTune$cp)
)
# Visualize the tree
rpart.plot(final_tree, type = 2, extra = 104, fallen.leaves = TRUE)
#Variable importance
final_tree$variable.importance## delivery_time_bucket item_price brand_id
## 2583.266571 593.999012 130.967684
## size_category order_date_bucket item_color_simple
## 10.177667 7.634840 2.008834# Predict class labels on the training set
pred_class <- predict(final_tree, newdata = orders_train, type = "class")
# Actual values
actual <- orders_train$return
# Create a confusion table
conf_matrix <- table(Predicted = pred_class, Actual = actual)
conf_matrix## Actual
## Predicted No Yes
## No 14651 5982
## Yes 13273 17370#Show error rate
error_rate <- sum(pred_class != actual) / length(actual)
error_rate## [1] 0.3755168We’ve reached our highest accuracy yet, but there still isn’t a meaningful improvement from logistic regression. Let’s see if we can do better.
Random Forests
Random forests models create a series of classification trees. These trees are artificially limited in the number of predictors they can ‘choose’ from, meaning that predictive value can be drawn from less dominant predictors. For example, in our single classification tree, delivery time and item price had outsize importance.
As always we start with packages and data:
library(rpart) # for building decision trees
library(rpart.plot) # for visualizing trees
library(lubridate) # date management
library(ggplot2) # charts
library(dplyr) # data handling
library(glmnet) # implementation of the methods
library(caret) # machine-learning workflow
library(car) # multicolinearity check
library(randomForest) # Makes random forests## Warning: package 'randomForest' was built under R version 4.4.3## randomForest 4.7-1.2## Type rfNews() to see new features/changes/bug fixes.##
## Attaching package: 'randomForest'## The following object is masked from 'package:dplyr':
##
## combine## The following object is masked from 'package:ggplot2':
##
## margin# set working directory:
setwd("C:/Users/jedel/OneDrive/Desktop/Data for ML class/Individual Assignment")
#Load data
orders_clean <- read.csv("orders_clean.csv", sep = ",")Next, we set our seed and split our data into test and training:
# Set seed
set.seed(8888) # for reproducibility. I use seed "8888"
# Determine sample size
train_size <- floor(0.7 * nrow(orders_clean))
# Generate random row indices for training set
train_indices <- sample(seq_len(nrow(orders_clean)), size = train_size)
# Split the dataset
orders_train <- orders_clean[train_indices, ]
orders_test <- orders_clean[-train_indices, ]
# Ensure return is a factor
orders_train$return <- as.factor(orders_train$return)
orders_test$return <- as.factor(orders_test$return)Next, we’ll train our baseline model. I have an off-the-shelf non-gaming laptop, so I will be limited by computing power. I limit initial tree count to 500:
# Fit base model
rf_model <- randomForest(
return ~ size_category + item_color_simple + brand_id + item_price +
user_title + user_state + dob_bin + delivery_time_bucket +
order_date_bucket + account_age_bucket,
data = orders_train,
ntree = 500, # number of trees
mtry = 3, # number of variables considered at each split
importance = TRUE # to check variable importance later
)
# Predict on test set
rf_preds <- predict(rf_model, orders_test, type = "class")
# Error rate
rf_error_rate <- mean(rf_preds != orders_test$return)
rf_error_rate## [1] 0.3633964Before optimizing the random forest model, we will see if bagging has an impact. Bagging sets the number of predictors a available to a given tree to equal the total number of predictors. In our case, we are going fro 3 to 10 predictors.
Bagging follows a similar set of operations:
# Bagging (all predictors at each split)
rf_bag <- randomForest(
return ~ size_category + item_color_simple + brand_id + item_price +
user_title + user_state + dob_bin + delivery_time_bucket +
order_date_bucket + account_age_bucket,
data = orders_train,
ntree = 500,
mtry = 10, # use all predictors at each split
importance = TRUE
)
# Check on test data
bag_preds <- predict(rf_bag, orders_test, type = "class")
bag_error_rate <- mean(bag_preds != orders_test$return)
bag_error_rate## [1] 0.374636This is a little bit worse than our previous model, so we will proceed with optimizing the original model. We will test several variables to first attempt to optimize the mtry and ntree parameters, then see if we can reduce the number of variables used.
Let’s start by increasing the prdictors a tree can split into to 4:
# Build random forest with 4 variables considered per split
rf_4split <- randomForest(
return ~ size_category + item_color_simple + brand_id + item_price +
user_title + user_state + dob_bin + delivery_time_bucket +
order_date_bucket + account_age_bucket,
data = orders_train,
ntree = 500, # number of trees (default is 500)
mtry = 4, # number of variables considered at each split
importance = TRUE # to check variable importance later
)
# Predict on test set
rf_preds <- predict(rf_4split, orders_test, type = "class")
# Error rate
rf_error_rate <- mean(rf_preds != orders_test$return)
rf_error_rate## [1] 0.3684929That’s a bit worse than our base model. Let’s try 2 instead of 4:
# Build random forest with 2 variables considered per split
rf_2split <- randomForest(
return ~ size_category + item_color_simple + brand_id + item_price +
user_title + user_state + dob_bin + delivery_time_bucket +
order_date_bucket + account_age_bucket,
data = orders_train,
ntree = 500, # number of trees (default is 500)
mtry = 2, # number of variables considered at each split
importance = TRUE # to check variable importance later
)
# Predict on test set
rf_preds <- predict(rf_2split, orders_test, type = "class")
# Error rate
rf_error_rate <- mean(rf_preds != orders_test$return)
rf_error_rate## [1] 0.358755This is our best yet! Now let’s test out a variation where we grow 1000 rather than 500 trees:
# Build random forest with mtry = 2 and ntree = 1000
rf_2split1000 <- randomForest(
return ~ size_category + item_color_simple + brand_id + item_price +
user_title + user_state + dob_bin + delivery_time_bucket +
order_date_bucket + account_age_bucket,
data = orders_train,
ntree = 1000, # number of trees (default is 500)
mtry = 2, # number of variables considered at each split
importance = TRUE # to check variable importance later
)
# Predict on test set
rf_preds <- predict(rf_2split1000, orders_test, type = "class")
# Error rate
rf_error_rate <- mean(rf_preds != orders_test$return)
rf_error_rate## [1] 0.3604842We see that this is slightly worse, so we’ll stay at 500. This is inconsistent, and I experienced different results in testing with a marginally different dataset. In the case of this project, we’ll keep it at 500, and my CPUs will be relieved.
Now let’s see if we can optimize any of the predictors. We start by looking at variable importance. Higher mean indicates predictive importance, and higher gini indicates relative importance of a variable for splitting/partitioning data:
# Since it's the best we've gotten so far, let's look at variable importance for the rf_2split1000 model
varImpPlot(rf_2split)
Looking at the relative impact of all variables across both accuracy and gini, we see that user_title is relatively unimportant. Let’s try excluding it:
# Build random forest with mtry = 2 and without user_title
rf_2splitnotitle <- randomForest(
return ~ size_category + item_color_simple + brand_id + item_price + user_state + dob_bin + delivery_time_bucket +
order_date_bucket + account_age_bucket,
data = orders_train,
ntree = 500, # number of trees (default is 500)
mtry = 2, # number of variables considered at each split
importance = TRUE # to check variable importance later
)
# Predict on test set
rf_preds <- predict(rf_2splitnotitle, orders_test, type = "class")
# Error rate
rf_error_rate <- mean(rf_preds != orders_test$return)
rf_error_rate## [1] 0.3586185This is our best yet! Before optimizing further, let’s try a neural network.
Neural network
As a final step, we attempted to fit a neural network. Due to time limitations, we were not able to explore this approach fully, but were able to fit a limited model with 5 hidden layers. In 50 iterations, we were able to reduce loss from 41412 to 30906. I am personally interested in exploring further, but this is outside the scope of this project.
Based on the above, a random forest is the best model we were able to identify over the course of this project.