Olist Delivery Performance — Analysis, Segmentation & Prediction

A fully reproducible walkthrough: every table and chart shows the code that produced it

How to read this report. Unlike a slide deck, this document is executable: every figure below is generated by the R code shown directly above it (click Hide / Show on any block, or Code ▸ at the top right to toggle all). Nothing is a pre-rendered image — if you change the data and re-knit, the charts change with it. Work flows in three linked parts: (1) diagnose delivery performance, (2) segment sellers with clustering, (3) predict late deliveries with Random Forest & XGBoost.

1 Setup: packages, theme, helpers

We load the tidyverse stack plus the modelling/clustering packages, and define a single house style (theme_olist) and colour palette reused by every chart.

library(readr)      # fast CSV reading
library(dplyr)      # data wrangling
library(tidyr)      # reshaping
library(ggplot2)    # all charts
library(scales)     # axis formatting (%, currency, commas)
library(forcats)    # factor ordering for ranked bars
library(stringr)
library(lubridate)  # timestamp parsing
library(cluster)    # pam(), silhouette()
library(xgboost)    # gradient boosting
library(randomForest)
library(pROC)       # ROC / AUC
library(Matrix)     # sparse design matrix
library(knitr)

# PowerPoint-style palette: strong purple / blue / green / orange
PURPLE <- "#6A3D9A"; BLUE <- "#1F78B4"; GREEN <- "#33A02C"
ORANGE <- "#FF7F00"; RED  <- "#E31A1C"
PAL    <- c(PURPLE, BLUE, GREEN, ORANGE, RED, "#A6CEE3")

# One reusable, minimal theme. NOTE: randomForest masks ggplot2::margin(),
# so we qualify it explicitly.
theme_olist <- function(base = 12) {
  theme_minimal(base_size = base) +
    theme(
      plot.title    = element_text(face = "bold", size = rel(1.2)),
      plot.subtitle = element_text(colour = "grey35", margin = ggplot2::margin(b = 8)),
      panel.grid.minor = element_blank(),
      panel.grid.major = element_line(colour = "grey90"),
      legend.position  = "top"
    )
}
theme_set(theme_olist())
set.seed(42)   # reproducibility everywhere

2 Data preparation

This is “how I got the data”. Olist ships as separate CSVs (orders, items, customers, sellers, products, reviews, geolocation). Two rules prevent silent row-duplication on join:

1. Geolocation has many rows per ZIP — collapse to one mean lat/lng per ZIP prefix before joining.
2. Reviews can have several rows per order — collapse to one mean score per order before joining.

We keep order_status == "delivered" and parse the five timestamps with truncated = 3 (the estimated date is date-only and would otherwise fail to parse).

# Resolve the data folder (CSVs may live in "data/" or "data/olist dataset/")
DATA_DIR <- c("data/olist dataset", "data") |>
  (\(d) d[file.exists(file.path(d, "olist_orders_dataset.csv"))][1])()

ts <- c("order_purchase_timestamp","order_approved_at","order_delivered_carrier_date",
        "order_delivered_customer_date","order_estimated_delivery_date")

orders_raw <- read_csv(file.path(DATA_DIR, "olist_orders_dataset.csv"),
  show_col_types = FALSE,
  col_types = cols(.default = col_guess(),       # force timestamps to character
    order_purchase_timestamp = col_character(), order_approved_at = col_character(),
    order_delivered_carrier_date = col_character(),
    order_delivered_customer_date = col_character(),
    order_estimated_delivery_date = col_character()))
items_raw   <- read_csv(file.path(DATA_DIR, "olist_order_items_dataset.csv"), show_col_types = FALSE)
customers   <- read_csv(file.path(DATA_DIR, "olist_customers_dataset.csv"), show_col_types = FALSE)
sellers     <- read_csv(file.path(DATA_DIR, "olist_sellers_dataset.csv"), show_col_types = FALSE)
products    <- read_csv(file.path(DATA_DIR, "olist_products_dataset.csv"), show_col_types = FALSE)
translation <- read_csv(file.path(DATA_DIR, "product_category_name_translation.csv"), show_col_types = FALSE)
reviews_raw <- read_csv(file.path(DATA_DIR, "olist_order_reviews_dataset.csv"), show_col_types = FALSE)
geo_raw     <- read_csv(file.path(DATA_DIR, "olist_geolocation_dataset.csv"), show_col_types = FALSE)

# (1) one mean coordinate per ZIP prefix  -> safe to join
geo_zip <- geo_raw |>
  group_by(geolocation_zip_code_prefix) |>
  summarise(lat = mean(geolocation_lat), lng = mean(geolocation_lng), .groups = "drop") |>
  mutate(zip = as.integer(geolocation_zip_code_prefix)) |> select(zip, lat, lng)

# (2) one mean review score per order  -> no duplication
rev_order <- reviews_raw |>
  group_by(order_id) |>
  summarise(review_score = mean(review_score, na.rm = TRUE), .groups = "drop")

# customers / sellers with coordinates attached
customers <- customers |> mutate(zip = as.integer(customer_zip_code_prefix)) |>
  left_join(geo_zip, by = "zip") |> rename(cust_lat = lat, cust_lng = lng)
sellers   <- sellers   |> mutate(zip = as.integer(seller_zip_code_prefix)) |>
  left_join(geo_zip, by = "zip") |> rename(sell_lat = lat, sell_lng = lng)

# products + English category + package volume
products <- products |>
  left_join(translation, by = "product_category_name") |>
  mutate(category_en = coalesce(product_category_name_english, product_category_name, "unknown"),
         pkg_volume_cm3 = product_length_cm * product_height_cm * product_width_cm)

dbtw <- function(a, b) as.numeric(difftime(a, b, units = "days"))  # interval in days

orders <- orders_raw |>
  filter(order_status == "delivered") |>
  mutate(across(all_of(ts), ~ ymd_hms(., quiet = TRUE, truncated = 3)),
    # --- delivery stages (days) ---
    approval_time           = dbtw(order_approved_at, order_purchase_timestamp),
    handling_time           = dbtw(order_delivered_carrier_date, order_approved_at),
    transit_time            = dbtw(order_delivered_customer_date, order_delivered_carrier_date),
    total_delivery_time     = dbtw(order_delivered_customer_date, order_purchase_timestamp),
    promised_window         = dbtw(order_estimated_delivery_date, order_purchase_timestamp),
    delivery_variance       = dbtw(order_delivered_customer_date, order_estimated_delivery_date),
    purchase_month          = floor_date(order_purchase_timestamp, "month"),
    # --- service-level indicators ---
    on_time_delivery = order_delivered_customer_date <= order_estimated_delivery_date,
    late_delivery    = order_delivered_customer_date >  order_estimated_delivery_date,
    severely_late    = delivery_variance > 7,
    # --- validity flags: exclude impossible/negative intervals from metrics ---
    valid_handling = !is.na(handling_time) & handling_time >= 0,
    valid_transit  = !is.na(transit_time)  & transit_time  >= 0,
    valid_total    = !is.na(total_delivery_time) & total_delivery_time >= 0)

# Straight-line distance (Haversine) — an approximation, NOT road distance
haversine_km <- function(la1, lo1, la2, lo2) { p <- pi/180; R <- 6371
  a <- sin((la2-la1)*p/2)^2 + cos(la1*p)*cos(la2*p)*sin((lo2-lo1)*p/2)^2
  2 * R * asin(pmin(1, sqrt(a))) }

# Item-level table enriched with product, seller, and order outcomes
items <- items_raw |>
  left_join(products, by = "product_id") |>
  left_join(sellers |> select(seller_id, seller_state, sell_lat, sell_lng), by = "seller_id")

# Order-level characteristics (number of items/sellers, totals, dominant category)
dominant_cat <- items |> group_by(order_id, category_en) |>
  summarise(v = sum(price, na.rm = TRUE), .groups = "drop_last") |>
  slice_max(v, n = 1, with_ties = FALSE) |> ungroup() |>
  select(order_id, dominant_category = category_en)

order_chars <- items |> group_by(order_id) |>
  summarise(number_of_items = n(), number_of_products = n_distinct(product_id),
            number_of_sellers = n_distinct(seller_id),
            total_order_value = sum(price, na.rm = TRUE),
            total_freight_value = sum(freight_value, na.rm = TRUE),
            total_item_weight = sum(product_weight_g, na.rm = TRUE),
            primary_seller_state = first(seller_state),
            sell_lat = first(sell_lat), sell_lng = first(sell_lng), .groups = "drop") |>
  left_join(dominant_cat, by = "order_id")

# Master order table (left_join chain) + customer geo + distance + bands
orders_master <- orders |>
  left_join(customers |> select(customer_id, customer_state, customer_zip = zip,
                                cust_lat, cust_lng), by = "customer_id") |>
  left_join(order_chars, by = "order_id") |>
  left_join(rev_order, by = "order_id") |>
  mutate(seller_customer_distance_km = haversine_km(cust_lat, cust_lng, sell_lat, sell_lng),
         interstate_delivery = customer_state != primary_seller_state,
         distance_band = cut(seller_customer_distance_km,
            c(-Inf,100,300,700,1500,Inf),
            labels = c("0-100 km","101-300 km","301-700 km","701-1500 km",">1500 km")),
         seller_group_cnt = cut(number_of_sellers, c(0,1,2,Inf),
            labels = c("1 seller","2 sellers","3+ sellers")))

# Analysis frame: delivered orders with a valid total time
adf <- orders_master |> filter(valid_total)
nrow(adf)

## [1] 96470

3 Part 1 · Overall delivery performance

3.1 Headline KPIs

q <- function(x, p) quantile(x, p, na.rm = TRUE)
kpis <- tibble(
  Metric = c("Delivered orders","Mean delivery time (d)","Median delivery time (d)",
             "P90 delivery time (d)","P95 delivery time (d)","On-time rate","Late rate",
             "Severely-late rate (>7d)","Median handling (d)","Median transit (d)"),
  Value = c(nrow(adf),
            round(mean(adf$total_delivery_time),1), round(median(adf$total_delivery_time),1),
            round(q(adf$total_delivery_time,.90),1), round(q(adf$total_delivery_time,.95),1),
            percent(mean(adf$on_time_delivery),0.1), percent(mean(adf$late_delivery),0.1),
            percent(mean(adf$severely_late),0.1),
            round(median(adf$handling_time[adf$valid_handling]),2),
            round(median(adf$transit_time[adf$valid_transit]),2)))
kable(kpis, caption = "Headline delivery KPIs.")

Headline delivery KPIs.

Metric	Value
Delivered orders	96470
Mean delivery time (d)	12.6
Median delivery time (d)	10.2
P90 delivery time (d)	23.1
P95 delivery time (d)	29.3
On-time rate	91.9%
Late rate	8.1%
Severely-late rate (>7d)	3.5%
Median handling (d)	1.85
Median transit (d)	7.1

3.2 Distribution of delivery time

ggplot(filter(adf, total_delivery_time <= 60), aes(total_delivery_time)) +
  geom_histogram(binwidth = 1, fill = BLUE, colour = "white", linewidth = .1) +
  geom_vline(xintercept = median(adf$total_delivery_time), colour = ORANGE, linewidth = 1) +
  scale_y_continuous(labels = comma) +
  labs(title = "Most orders arrive in about a week, but the tail runs long",
       subtitle = "Total delivery time (purchase → customer), capped at 60 days",
       x = "Total delivery time (days)", y = "Orders")

3.3 Distribution of delivery variance (vs. the promise)

ggplot(filter(adf, between(delivery_variance, -40, 40)),
       aes(delivery_variance, fill = late_delivery)) +
  geom_histogram(binwidth = 1, colour = "white", linewidth = .05) +
  geom_vline(xintercept = 0, linetype = "dashed", linewidth = .9) +
  scale_fill_manual(values = c(`FALSE` = GREEN, `TRUE` = RED),
                    labels = c("Early / on-time","Late"), name = NULL) +
  scale_y_continuous(labels = comma) +
  labs(title = "Olist usually delivers comfortably ahead of its promise",
       subtitle = "Variance = actual − estimated date (negative = early). Capped ±40 days",
       x = "Delivery variance (days)", y = "Orders")

3.4 Where is the time spent? (the key chart)

adf |>
  transmute(Approval = ifelse(valid_handling, approval_time, NA),
            `Seller handling` = ifelse(valid_handling, handling_time, NA),
            `Carrier transit` = ifelse(valid_transit, transit_time, NA)) |>
  pivot_longer(everything(), names_to = "stage", values_to = "days") |>
  filter(!is.na(days), days >= 0, days <= 40) |>
  mutate(stage = factor(stage, c("Approval","Seller handling","Carrier transit"))) |>
  ggplot(aes(stage, days, fill = stage)) +
  geom_boxplot(outlier.alpha = .05, width = .55) +
  scale_fill_manual(values = c(PURPLE, ORANGE, BLUE), guide = "none") +
  coord_cartesian(ylim = c(0, 30)) +
  labs(title = "Carrier transit is the largest — and most variable — stage",
       subtitle = "Distribution of each stage in days (outliers > 40d hidden)",
       x = NULL, y = "Days")

Proves: carrier transit, not seller handling, is where time is spent.

4 Part 2 · Performance over time

monthly <- adf |> group_by(purchase_month) |>
  summarise(orders = n(),
            median_dt = median(total_delivery_time),
            p90_dt    = quantile(total_delivery_time, .90),
            on_time   = mean(on_time_delivery),
            handling  = median(handling_time[valid_handling]),
            transit   = median(transit_time[valid_transit]), .groups = "drop") |>
  filter(!is.na(purchase_month), orders >= 50)   # drop sparse edge months

4.1 Volume vs. median delivery time

scl <- max(monthly$orders) / max(monthly$median_dt)   # dual-axis scaling factor
ggplot(monthly, aes(purchase_month)) +
  geom_col(aes(y = orders), fill = BLUE, alpha = .45) +
  geom_line(aes(y = median_dt * scl), colour = ORANGE, linewidth = 1.2) +
  geom_point(aes(y = median_dt * scl), colour = ORANGE) +
  scale_y_continuous("Delivered orders", labels = comma,
    sec.axis = sec_axis(~ . / scl, name = "Median delivery time (days)")) +
  scale_x_datetime(date_labels = "%b %Y", date_breaks = "2 months") +
  labs(title = "Volume grew while median delivery time stayed broadly stable",
       subtitle = "Bars = orders; line = median delivery time", x = NULL) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

4.2 On-time rate and the handling/transit split

ggplot(monthly, aes(purchase_month, on_time)) +
  geom_line(colour = GREEN, linewidth = 1.2) + geom_point(colour = GREEN) +
  geom_hline(yintercept = mean(adf$on_time_delivery), linetype = "dashed", colour = "grey50") +
  scale_y_continuous(labels = percent) +
  scale_x_datetime(date_labels = "%b %Y", date_breaks = "2 months") +
  labs(title = "On-time rate is high but dips in specific months",
       subtitle = "Dashed = overall average", x = NULL, y = "On-time rate") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

monthly |> select(purchase_month, Handling = handling, Transit = transit) |>
  pivot_longer(-purchase_month, names_to = "stage", values_to = "days") |>
  ggplot(aes(purchase_month, days, colour = stage)) +
  geom_line(linewidth = 1.2) + geom_point() +
  scale_colour_manual(values = c(Handling = ORANGE, Transit = BLUE), name = NULL) +
  scale_x_datetime(date_labels = "%b %Y", date_breaks = "2 months") +
  labs(title = "Transit dominates handling every single month",
       subtitle = "Monthly median handling vs. transit", x = NULL, y = "Median days") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

Correlation between volume and delivery time is weak and not evidence that volume causes delay.

5 Part 3 · Geography

A reusable summary helper keeps every geographic cut consistent:

delivery_metrics <- function(df) df |> summarise(
  orders = n(),
  median_dt = median(total_delivery_time, na.rm = TRUE),
  p90_dt    = quantile(total_delivery_time, .90, na.rm = TRUE),
  on_time   = mean(on_time_delivery, na.rm = TRUE),
  late_rate = mean(late_delivery, na.rm = TRUE),
  median_freight = median(total_freight_value, na.rm = TRUE),
  median_distance = median(seller_customer_distance_km, na.rm = TRUE), .groups = "drop")

MIN_N <- 30   # minimum orders before we rank a state/route
state_rank <- adf |> group_by(customer_state) |> delivery_metrics() |>
  filter(orders >= MIN_N)

5.1 States ranked by median delivery time

state_rank |>
  mutate(customer_state = fct_reorder(customer_state, median_dt)) |>
  ggplot(aes(customer_state, median_dt)) +
  geom_col(fill = PURPLE) +
  geom_text(aes(label = round(median_dt,1)), hjust = -.15, size = 3) +
  coord_flip() + scale_y_continuous(expand = expansion(mult = c(0,.1))) +
  labs(title = "Northern states wait roughly twice as long as the southeast",
       subtitle = "Median total delivery time by customer state (n ≥ 30)",
       x = NULL, y = "Median delivery time (days)")

5.2 States ranked by late-delivery rate

state_rank |>
  mutate(customer_state = fct_reorder(customer_state, late_rate)) |>
  ggplot(aes(customer_state, late_rate)) +
  geom_col(fill = ORANGE) +
  geom_text(aes(label = percent(late_rate,.1)), hjust = -.1, size = 3) +
  coord_flip() + scale_y_continuous(labels = percent, expand = expansion(mult = c(0,.12))) +
  labs(title = "Late-delivery rate is far higher outside the southeast",
       subtitle = "Share delivered after the estimated date (n ≥ 30)",
       x = NULL, y = "Late-delivery rate")

5.3 Seller-state → customer-state route heatmap

route_metrics <- adf |> filter(!is.na(primary_seller_state)) |>
  group_by(seller_state = primary_seller_state, customer_state) |>
  summarise(orders = n(), median_dt = median(total_delivery_time), .groups = "drop")

top_sell <- adf |> count(primary_seller_state, sort = TRUE) |>
  filter(!is.na(primary_seller_state)) |> slice_head(n = 10) |> pull(primary_seller_state)
top_cust <- adf |> count(customer_state, sort = TRUE) |> slice_head(n = 12) |> pull(customer_state)

route_metrics |>
  filter(seller_state %in% top_sell, customer_state %in% top_cust, orders >= 10) |>
  ggplot(aes(customer_state, fct_rev(seller_state), fill = median_dt)) +
  geom_tile(colour = "white") +
  geom_text(aes(label = round(median_dt)), size = 2.8, colour = "grey15") +
  scale_fill_gradient(low = "#E5F5E0", high = PURPLE, name = "Median days") +
  labs(title = "Slowest routes pair southern/southeastern sellers with northern buyers",
       subtitle = "Median delivery time by seller (rows) → customer (cols) state; n ≥ 10",
       x = "Customer state", y = "Seller state")

5.4 Delivery performance by distance band

band <- adf |> filter(!is.na(distance_band)) |> group_by(distance_band) |> delivery_metrics()
ggplot(band, aes(distance_band)) +
  geom_col(aes(y = median_dt), fill = BLUE, alpha = .5) +
  geom_text(aes(y = median_dt, label = round(median_dt,1)), vjust = -.5, colour = BLUE, size = 3) +
  geom_line(aes(y = late_rate * max(band$median_dt)/max(band$late_rate), group = 1),
            colour = RED, linewidth = 1.1) +
  geom_point(aes(y = late_rate * max(band$median_dt)/max(band$late_rate)), colour = RED) +
  scale_y_continuous("Median delivery time (days)",
    sec.axis = sec_axis(~ . * max(band$late_rate)/max(band$median_dt),
                        name = "Late-delivery rate", labels = percent)) +
  labs(title = "Delivery time rises steadily with seller–customer distance",
       subtitle = "Bars = median time; red line = late rate. Straight-line (Haversine) distance",
       x = "Distance band")

Proves: delivery time scales with distance; northern/interstate routes are structurally slower. Distance is straight-line, an approximation of road distance.

6 Part 4 · Seller-associated bottlenecks

Seller timestamps are recorded at order level, so these are associations, not strict attribution. We score each seller over its delivered orders:

seller_order <- items |> select(order_id, seller_id, freight_value) |>
  left_join(adf |> select(order_id, handling_time, transit_time, valid_handling,
                          valid_transit, late_delivery, severely_late, review_score),
            by = "order_id") |>
  filter(!is.na(handling_time) | !is.na(late_delivery)) |>
  distinct(seller_id, order_id, .keep_all = TRUE) |>
  filter(order_id %in% adf$order_id)

seller_rank <- seller_order |> group_by(seller_id) |>
  summarise(orders = n_distinct(order_id),
            median_handling = median(handling_time[valid_handling], na.rm = TRUE),
            median_transit  = median(transit_time[valid_transit],  na.rm = TRUE),
            late_rate = mean(late_delivery, na.rm = TRUE), .groups = "drop") |>
  filter(orders >= MIN_N)

6.1 Top 20 sellers by handling time, and by late rate

seller_rank |> slice_max(median_handling, n = 20) |>
  mutate(seller = fct_reorder(str_sub(seller_id,1,8), median_handling)) |>
  ggplot(aes(seller, median_handling)) +
  geom_col(fill = ORANGE) +
  geom_text(aes(label = round(median_handling,1)), hjust = -.15, size = 3) +
  coord_flip() + scale_y_continuous(expand = expansion(mult = c(0,.1))) +
  labs(title = "A minority of sellers sit on orders for days before handoff",
       subtitle = "Top 20 by median handling time (n ≥ 30; IDs truncated)",
       x = "Seller", y = "Median handling time (days)")

ggplot(seller_rank, aes(orders, late_rate)) +
  geom_point(alpha = .4, colour = BLUE) +
  geom_smooth(method = "loess", se = FALSE, colour = ORANGE, formula = y ~ x) +
  scale_x_log10(labels = comma) + scale_y_continuous(labels = percent) +
  labs(title = "High-volume sellers are not systematically worse on lateness",
       subtitle = "Each point a seller with n ≥ 30 orders (log x-axis)",
       x = "Delivered orders (log)", y = "Late-delivery rate")

7 Part 5 · Product & order characteristics

7.1 Slowest categories

cat_rank <- adf |> filter(!is.na(dominant_category), dominant_category != "unknown") |>
  group_by(dominant_category) |> delivery_metrics() |> filter(orders >= 100)

cat_rank |> slice_max(median_dt, n = 20) |>
  mutate(dominant_category = fct_reorder(dominant_category, median_dt)) |>
  ggplot(aes(dominant_category, median_dt)) +
  geom_col(fill = PURPLE) +
  geom_text(aes(label = round(median_dt,1)), hjust = -.15, size = 2.8) +
  coord_flip() + scale_y_continuous(expand = expansion(mult = c(0,.1))) +
  labs(title = "Bulky & furniture categories are slowest to arrive",
       subtitle = "Top 20 categories by median delivery time (n ≥ 100)",
       x = NULL, y = "Median delivery time (days)")

7.2 Are multi-seller orders slower? (No.)

adf |> filter(!is.na(seller_group_cnt), total_delivery_time <= 60) |>
  ggplot(aes(seller_group_cnt, total_delivery_time, fill = seller_group_cnt)) +
  geom_boxplot(outlier.alpha = .03, width = .55) +
  scale_fill_manual(values = c(GREEN, ORANGE, RED), guide = "none") +
  coord_cartesian(ylim = c(0, 40)) +
  labs(title = "Multi-seller orders are not slower — if anything slightly faster",
       subtitle = "They are rare and concentrate in the fast southeast",
       x = NULL, y = "Total delivery time (days)")

7.3 Distance vs. delivery time (the core relationship)

adf |> filter(!is.na(seller_customer_distance_km),
              seller_customer_distance_km <= 3500, total_delivery_time <= 60) |>
  slice_sample(n = 15000) |>
  ggplot(aes(seller_customer_distance_km, total_delivery_time)) +
  geom_point(alpha = .06, colour = BLUE) +
  geom_smooth(method = "lm", se = FALSE, colour = RED, formula = y ~ x) +
  scale_x_continuous(labels = comma) +
  labs(title = "Distance is a clear, near-linear driver of delivery time",
       subtitle = "Sampled orders; straight-line distance (approximation, not road distance)",
       x = "Seller–customer distance (km)", y = "Total delivery time (days)")

8 Part 6 · Delivery & customer satisfaction

rev_adf <- adf |> filter(!is.na(review_score)) |> mutate(rs = round(review_score))
rev_adf |>
  mutate(status = cut(delivery_variance, c(-Inf,-0.001,0.001,3,7,Inf),
            labels = c("Early","On-time","1-3 days late","4-7 days late",">7 days late"))) |>
  group_by(status) |> summarise(mean_review = mean(review_score), .groups = "drop") |>
  ggplot(aes(status, mean_review, fill = status)) +
  geom_col() + geom_text(aes(label = round(mean_review,2)), vjust = -.4, size = 3.2) +
  scale_fill_manual(values = c(GREEN,"#A6D96A","#FDAE61",ORANGE,RED), guide = "none") +
  ylim(0,5) +
  labs(title = "Each step later costs Olist review stars",
       subtitle = "Mean review score by delivery-status bucket (observational association)",
       x = NULL, y = "Mean review score")

Observational only — delay is associated with worse reviews; not proven to cause them.

9 Part 7 · Seller segmentation (clustering)

We group sellers by operational profile so the fixes can be targeted. Features are all continuous → k-means (with z-scaling) is the primary method, validated with PAM (k-medoids) and compared with Rand / Jaccard / Adjusted Rand. (Gower + hierarchical / k-prototypes would only be needed for mixed data.)

9.1 Build the per-seller feature matrix and scale it

item_seller <- items |>
  mutate(distance_km = NA_real_) |>
  left_join(adf |> select(order_id, customer_state, cust_lat, cust_lng, handling_time,
                          transit_time, valid_handling, valid_transit, late_delivery,
                          severely_late, review_score), by = "order_id") |>
  filter(order_id %in% adf$order_id) |>
  mutate(distance_km = haversine_km(cust_lat, cust_lng, sell_lat, sell_lng))

seller_feats <- item_seller |> distinct(seller_id, order_id, .keep_all = TRUE) |>
  group_by(seller_id) |>
  summarise(delivered_orders = n_distinct(order_id),
            handling_time_med = median(handling_time[valid_handling], na.rm = TRUE),
            transit_time_med  = median(transit_time[valid_transit],  na.rm = TRUE),
            late_rate = mean(late_delivery, na.rm = TRUE),
            distance_med_km = median(distance_km, na.rm = TRUE),
            review_mean = mean(review_score, na.rm = TRUE), .groups = "drop")

clust_pool <- seller_feats |>
  filter(delivered_orders >= 20, !is.na(handling_time_med), !is.na(transit_time_med),
         !is.na(distance_med_km), !is.na(review_mean)) |>
  mutate(log_orders = log10(delivered_orders))

feat_cols <- c("handling_time_med","transit_time_med","late_rate",
               "distance_med_km","log_orders","review_mean")
Xs <- scale(clust_pool[, feat_cols])   # z-score: REQUIRED for k-means
nrow(clust_pool)

## [1] 801

9.2 k-means + PAM, and how much they agree

set.seed(42)
K  <- 4                                   # chosen for interpretability (see note)
km <- kmeans(Xs, centers = K, nstart = 50)
pm <- pam(Xs, k = K)

# Pair-counting agreement (course definition): Rand = (a+b)/(a+b+c+d), Jaccard = a/(a+c+d)
pair_counts <- function(c1, c2) { s1 <- outer(c1,c1,`==`); s2 <- outer(c2,c2,`==`)
  ut <- upper.tri(s1)
  list(a = sum(s1[ut]& s2[ut]), b = sum(!s1[ut]&!s2[ut]),
       c = sum(s1[ut]&!s2[ut]), d = sum(!s1[ut]& s2[ut])) }
pcn <- pair_counts(km$cluster, pm$clustering)
tibble(Metric = c("Rand Index","Jaccard","Adjusted Rand"),
       Value  = round(c((pcn$a+pcn$b)/sum(unlist(pcn)),
                        pcn$a/(pcn$a+pcn$c+pcn$d),
                        mclust::adjustedRandIndex(km$cluster, pm$clustering)), 3)) |>
  kable(caption = "k-means vs PAM agreement — high values mean the segments are real structure.")

k-means vs PAM agreement — high values mean the segments are real structure.

Metric	Value
Rand Index	0.867
Jaccard	0.605
Adjusted Rand	0.663

clust_pool$km <- km$cluster
oh <- median(clust_pool$handling_time_med); ot <- median(clust_pool$transit_time_med)
clust_pool <- clust_pool |> group_by(km) |> mutate(seg_late = mean(late_rate)) |> ungroup() |>
  mutate(segment = case_when(
    seg_late >= quantile(unique(seg_late), .75) ~ "At-risk / unreliable",
    handling_time_med <= oh & transit_time_med <= ot ~ "Fast & local (model)",
    transit_time_med  >  ot & handling_time_med <= oh ~ "Fast handling, long-haul",
    TRUE ~ "Mainstream (avg)"))

Honest caveat: silhouette values are modest (~0.16–0.23); seller operations form a continuum, so k = 4 is an interpretability choice, not a number the data demands.

9.3 Segment fingerprints

clust_pool |> group_by(segment) |>
  summarise(across(all_of(feat_cols), mean), .groups = "drop") |>
  mutate(across(all_of(feat_cols), ~ as.numeric(scale(.)))) |>
  pivot_longer(all_of(feat_cols), names_to = "feature", values_to = "z") |>
  mutate(feature = recode(feature, handling_time_med="Handling", transit_time_med="Transit",
            late_rate="Late rate", distance_med_km="Distance", log_orders="Volume (log)",
            review_mean="Review")) |>
  ggplot(aes(feature, segment, fill = z)) +
  geom_tile(colour = "white") + geom_text(aes(label = round(z,1)), size = 3) +
  scale_fill_gradient2(low = BLUE, mid = "white", high = RED, name = "SD from\nmean") +
  labs(title = "Segment fingerprints: what makes each seller group distinct",
       subtitle = "Cluster-mean of each feature, standardised across segments",
       x = NULL, y = NULL)

clust_pool |> group_by(segment) |> summarise(late_rate = mean(late_rate), .groups="drop") |>
  mutate(segment = fct_reorder(segment, late_rate)) |>
  ggplot(aes(segment, late_rate, fill = segment)) +
  geom_col() + geom_text(aes(label = percent(late_rate,.1)), hjust = -.1, size = 3.2) +
  scale_fill_manual(values = PAL, guide = "none") +
  coord_flip() + scale_y_continuous(labels = percent, expand = expansion(mult = c(0,.15))) +
  labs(title = "Late-delivery rate differs sharply across seller segments",
       x = NULL, y = "Late-delivery rate")

10 Part 8 · Predicting late delivery (RF & XGBoost)

Question: at order time, will an order arrive late? Leakage control: we use only features known at/around order placement — we exclude realised handling and transit times (which would leak the outcome).

10.1 Build the design matrix and split 80/20

model_df <- adf |>
  mutate(approval_time = ifelse(is.na(approval_time) | approval_time < 0, 0, approval_time),
         purchase_m = month(order_purchase_timestamp),
         purchase_d = wday(order_purchase_timestamp),
         interstate = as.integer(interstate_delivery),
         late = as.integer(late_delivery)) |>
  filter(!is.na(promised_window), !is.na(seller_customer_distance_km),
         !is.na(customer_state), !is.na(primary_seller_state), !is.na(total_freight_value))

# collapse rare categories so one-hot stays manageable
top_cat <- model_df |> count(dominant_category, sort = TRUE) |> slice_head(n = 15) |> pull(dominant_category)
model_df <- model_df |> mutate(category_grp = ifelse(dominant_category %in% top_cat,
                                                     dominant_category, "other"))

feat_df <- model_df |>
  transmute(promised_window, distance_km = seller_customer_distance_km,
            freight = total_freight_value, order_value = total_order_value,
            n_items = number_of_items, n_sellers = number_of_sellers,
            weight = total_item_weight, approval_time, interstate,
            purchase_m, purchase_d,
            customer_state = factor(customer_state),
            seller_state = factor(primary_seller_state),
            category_grp = factor(category_grp))
y <- model_df$late

set.seed(42)
sel <- sample(seq_len(nrow(feat_df)), 0.8 * nrow(feat_df))
X   <- as.matrix(sparse.model.matrix(~ . - 1, feat_df))
Xtr <- X[sel, ]; Xte <- X[-sel, ]; ytr <- y[sel]; yte <- y[-sel]
c(train = length(sel), test = nrow(X) - length(sel), features = ncol(X))

##    train     test features 
##    76793    19199       74

10.2 XGBoost (boosting)

dtrain <- xgb.DMatrix(Xtr, label = ytr); dtest <- xgb.DMatrix(Xte, label = yte)
params <- list(objective = "binary:logistic", eval_metric = "auc", eta = .1,
               max_depth = 5, subsample = .8, colsample_bytree = .8,
               min_child_weight = 5,
               scale_pos_weight = sum(ytr == 0)/sum(ytr == 1))   # handle 8% imbalance
set.seed(42)
cv <- xgb.cv(params, dtrain, nrounds = 400, nfold = 5, early_stopping_rounds = 25, verbose = 0)
best_n <- which.max(cv$evaluation_log$test_auc_mean)
xgb_model <- xgb.train(params, dtrain, nrounds = best_n, verbose = 0)
xgb_prob  <- predict(xgb_model, dtest)

10.3 Random Forest (bagging) — the deck’s comparison

set.seed(42)
rf_idx <- sample(sel, min(30000, length(sel)))         # subsample for speed
rf_model <- randomForest(x = feat_df[rf_idx, ], y = factor(y[rf_idx]),
                         ntree = 300, mtry = floor(sqrt(ncol(feat_df))), importance = TRUE)
rf_prob  <- predict(rf_model, feat_df[-sel, ], type = "prob")[, "1"]

10.4 Evaluation: ROC and confusion matrix

xgb_auc <- as.numeric(auc(roc(yte, xgb_prob, quiet = TRUE)))
rf_auc  <- as.numeric(auc(roc(yte, rf_prob,  quiet = TRUE)))
xgb_pred <- as.integer(xgb_prob > .5)

bind_rows(
  tibble(Model="XGBoost", Accuracy=mean(xgb_pred==yte),
         Recall=sum(xgb_pred==1&yte==1)/sum(yte==1), AUC=xgb_auc),
  tibble(Model="Random Forest", Accuracy=mean((rf_prob>.5)==yte),
         Recall=sum((rf_prob>.5)==1&yte==1)/sum(yte==1), AUC=rf_auc)) |>
  mutate(across(where(is.numeric), ~round(.,3))) |>
  kable(caption = "Test performance. With ~8% late, AUC & recall matter more than accuracy.")

Test performance. With ~8% late, AUC & recall matter more than accuracy.

Model	Accuracy	Recall	AUC
XGBoost	0.761	0.688	0.805
Random Forest	0.921	0.046	0.778

rx <- roc(yte, xgb_prob, quiet = TRUE); rr <- roc(yte, rf_prob, quiet = TRUE)
bind_rows(
  tibble(fpr = 1-rx$specificities, tpr = rx$sensitivities, m = sprintf("XGBoost (AUC %.2f)", xgb_auc)),
  tibble(fpr = 1-rr$specificities, tpr = rr$sensitivities, m = sprintf("Random Forest (AUC %.2f)", rf_auc))) |>
  ggplot(aes(fpr, tpr, colour = m)) +
  geom_abline(linetype = "dashed", colour = "grey60") + geom_line(linewidth = 1.1) +
  scale_colour_manual(values = c(ORANGE, GREEN), name = NULL) +
  scale_x_continuous(labels = percent) + scale_y_continuous(labels = percent) +
  labs(title = "Both tree ensembles predict late delivery well above chance",
       subtitle = "ROC on the held-out test set (dashed = random)",
       x = "False positive rate", y = "True positive rate (recall)")

expand.grid(Actual = c("On-time","Late"), Predicted = c("On-time","Late")) |>
  mutate(n = c(sum(xgb_pred==0&yte==0), sum(xgb_pred==0&yte==1),
               sum(xgb_pred==1&yte==0), sum(xgb_pred==1&yte==1))) |>
  group_by(Actual) |> mutate(pct = n/sum(n)) |> ungroup() |>
  ggplot(aes(Predicted, fct_rev(Actual), fill = pct)) +
  geom_tile(colour = "white") +
  geom_text(aes(label = paste0(comma(n), "\n", percent(pct,1))), size = 4, colour = "grey15") +
  scale_fill_gradient(low = "#EDE7F6", high = PURPLE, labels = percent, name = "Row %") +
  labs(title = "XGBoost confusion matrix — late-delivery prediction",
       subtitle = sprintf("Rows = actual, cols = predicted | catches %.0f%% of late orders",
                          100*sum(xgb_pred==1&yte==1)/sum(yte==1)),
       x = "Predicted", y = "Actual")

10.5 What drives late risk?

xgb.importance(colnames(Xtr), model = xgb_model) |> as_tibble() |>
  slice_max(Gain, n = 20) |> mutate(Feature = fct_reorder(Feature, Gain)) |>
  ggplot(aes(Feature, Gain)) +
  geom_col(fill = ORANGE) +
  geom_text(aes(label = percent(Gain,.1)), hjust = -.1, size = 3) +
  coord_flip() + scale_y_continuous(labels = percent, expand = expansion(mult = c(0,.15))) +
  labs(title = "Season and promise-window generosity drive late risk most",
       subtitle = "Top 20 features by XGBoost gain; one-hot state/category levels shown individually",
       x = NULL, y = "Gain (share of total)")

Proves: lateness is predictable from order-time info alone (AUC ≈ 0.81). Suggests: season and promise-window calibration are the dominant levers. Test: whether a deployed late-risk score + recalibrated ETAs actually cut the late rate.

11 Recommendations

1. Re-engineer promise dates by route / distance / weight / season — the model confirms these are the dominant drivers; replace the single national rule.
2. Attack transit on the slowest interstate / long-haul routes (alternative carriers, regional hubs) — transit is the bottleneck and scales with distance.
3. Manage the at-risk seller segment with SLAs and handling-time anomaly alerts.
4. Deploy the late-risk model for proactive customer comms and prioritised handling on flagged orders; refresh ETAs after carrier handoff.
5. Monitor the tail (P90, severely-late rate), not just the average.

12 Limitations & honest framing

• Order-level timestamps → seller results are associations, not attribution.
• Straight-line distance approximates road distance (Haversine).
• No carrier identity → route-level inference only.
• Class imbalance / threshold shape the model’s accuracy–recall trade-off; AUC is the fair summary and the operating threshold should match the cost of a missed late order vs. a false alarm.
• Clustering is descriptive and subjective; correlation ≠ causation throughout.

---

Reproducibility: this single .Rmd loads the raw CSVs and regenerates every table, chart, and model from scratch — there are no pre-baked images. The standalone batch scripts delivery_performance_analysis.R, seller_clustering_analysis.R, and late_delivery_model.R produce the same outputs as PNG/CSV files.

## R version 4.6.0 (2026-04-24)
## Platform: aarch64-apple-darwin23
## Running under: macOS Tahoe 26.0.1
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/4.6/Resources/lib/libRblas.0.dylib 
## LAPACK: /Library/Frameworks/R.framework/Versions/4.6/Resources/lib/libRlapack.dylib;  LAPACK version 3.12.1
## 
## locale:
## [1] en_GB/en_GB/en_GB/C/en_GB/en_GB
## 
## time zone: Europe/Warsaw
## tzcode source: internal
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] knitr_1.51           Matrix_1.7-5         pROC_1.19.0.1       
##  [4] randomForest_4.7-1.2 xgboost_3.2.1.1      cluster_2.1.8.2     
##  [7] lubridate_1.9.5      stringr_1.6.0        forcats_1.0.1       
## [10] scales_1.4.0         ggplot2_4.0.3        tidyr_1.3.2         
## [13] dplyr_1.2.1          readr_2.2.0         
## 
## loaded via a namespace (and not attached):
##  [1] sass_0.4.10        generics_0.1.4     stringi_1.8.7      lattice_0.22-9    
##  [5] hms_1.1.4          digest_0.6.39      magrittr_2.0.5     evaluate_1.0.5    
##  [9] grid_4.6.0         timechange_0.4.0   RColorBrewer_1.1-3 fastmap_1.2.0     
## [13] jsonlite_2.0.0     mclust_6.1.2       mgcv_1.9-4         purrr_1.2.2       
## [17] codetools_0.2-20   jquerylib_0.1.4    cli_3.6.6          rlang_1.2.0       
## [21] crayon_1.5.3       splines_4.6.0      bit64_4.8.2        withr_3.0.3       
## [25] cachem_1.1.0       yaml_2.3.12        parallel_4.6.0     tools_4.6.0       
## [29] tzdb_0.5.0         vctrs_0.7.3        R6_2.6.1           lifecycle_1.0.5   
## [33] bit_4.6.0          vroom_1.7.1        pkgconfig_2.0.3    pillar_1.11.1     
## [37] bslib_0.11.0       gtable_0.3.6       glue_1.8.1         data.table_1.18.4 
## [41] Rcpp_1.1.1-1.1     xfun_0.59          tibble_3.3.1       tidyselect_1.2.1  
## [45] farver_2.1.2       nlme_3.1-169       htmltools_0.5.9    labeling_0.4.3    
## [49] rmarkdown_2.31     compiler_4.6.0     S7_0.2.2