1 Purpose and method

This report describes delivery-performance patterns among delivered Olist orders. Each observation supplied to k-means is one order. The clustering features are delivery duration, performance against the promised date, and the order’s mean review score.

The analysis is descriptive. Cluster membership identifies groups of similar orders in these variables; it does not establish that delivery performance caused a particular review score.

2 Data preparation

2.1 Import the three source files

# Read only the order columns needed for this analysis.
orders <- read_csv(
  "olist_orders_dataset.csv",
  col_types = cols_only(
    order_id = col_character(),
    customer_id = col_character(),
    order_status = col_character(),
    order_purchase_timestamp = col_datetime(),
    order_delivered_customer_date = col_datetime(),
    order_estimated_delivery_date = col_datetime()
  )
)

# Read only the review identifier and score, keeping review_score numeric.
reviews <- read_csv(
  "olist_order_reviews_dataset.csv",
  col_types = cols_only(
    order_id = col_character(),
    review_score = col_double()
  )
)

# Read customer location fields, with identifiers retained as character data.
customers <- read_csv(
  "olist_customers_dataset.csv",
  col_types = cols_only(
    customer_id = col_character(),
    customer_state = col_character(),
    customer_city = col_character()
  )
)

2.2 Prevent duplicate-review row inflation

Some orders have more than one review record. Joining those records directly to orders would duplicate orders and incorrectly give them extra weight. Therefore, reviews are first collapsed to exactly one row per order_id.

# Group all review records belonging to the same order.
reviews_by_order <- reviews %>%
  group_by(order_id) %>%
  # Calculate one mean score per order and ignore missing review values.
  summarise(
    review_score = mean(review_score, na.rm = TRUE),
    .groups = "drop"
  ) %>%
  # Convert NaN to NA when an order has no usable score after aggregation.
  mutate(review_score = if_else(is.nan(review_score), NA_real_, review_score))

# Confirm that the collapsed table has no duplicated order identifiers.
stopifnot(!anyDuplicated(reviews_by_order$order_id))

2.3 Build the order-level analysis table

The required three-table chain is orders → reviews → customers. In code, joining three tables requires two join operations. Delivered orders remain one row each because both joined lookup tables now have unique keys.

# Start with delivered orders because the requested population is completed deliveries.
analysis_orders <- orders %>%
  filter(order_status == "delivered") %>%
  # First join: attach the one-row-per-order mean review score.
  left_join(reviews_by_order, by = "order_id") %>%
  # Second join: attach the customer's state and city.
  left_join(customers, by = "customer_id") %>%
  # Engineer elapsed delivery time in days, allowing fractional days.
  mutate(
    delivery_days = as.numeric(
      difftime(
        order_delivered_customer_date,
        order_purchase_timestamp,
        units = "days"
      )
    ),
    # Positive delay_days means late; negative delay_days means early.
    delay_days = as.numeric(
      difftime(
        order_delivered_customer_date,
        order_estimated_delivery_date,
        units = "days"
      )
    )
  ) %>%
  # Apply the specifically required missing-value filter for delivery measures.
  filter(!is.na(delivery_days), !is.na(delay_days))

# k-means cannot accept missing values, so remove orders without a usable review.
cluster_orders <- analysis_orders %>%
  filter(!is.na(review_score))

# Verify that each order remains a single clustering observation after the joins.
stopifnot(!anyDuplicated(cluster_orders$order_id))

# Record the number of usable order-level observations for the narrative.
n_cluster_orders <- nrow(cluster_orders)

After filtering, the clustering dataset contains 96,470 delivered orders, with each order appearing once.

3 Feature scaling

# Select the three required numeric features and convert them to a plain matrix.
cluster_features <- cluster_orders %>%
  select(delivery_days, delay_days, review_score) %>%
  as.matrix()

# Standardise every feature to mean 0 and standard deviation 1.
# This is essential because k-means uses Euclidean distances: without scaling,
# variables measured in days could dominate review_score merely due to units.
cluster_features_scaled <- scale(cluster_features)

# Confirm that scaling did not create non-finite values.
stopifnot(all(is.finite(cluster_features_scaled)))

4 Selecting the number of clusters

Both the elbow and average silhouette methods are inspected. Silhouette calculation requires pairwise distances, whose memory use grows approximately with the square of the number of observations. A reproducible sample of up to 3,000 orders is therefore used only to diagnose k. The final k-means model is fitted to every eligible order.

# Limit the diagnostic sample to a memory-safe size for silhouette distances.
diagnostic_n <- min(3000L, nrow(cluster_features_scaled))

# Reset the seed so the diagnostic sample is identical whenever the report runs.
set.seed(1)

# Sample row positions without replacement for the k diagnostics.
diagnostic_rows <- sample(
  seq_len(nrow(cluster_features_scaled)),
  size = diagnostic_n,
  replace = FALSE
)

# Extract the same scaled features for the reproducible diagnostic sample.
k_diagnostic_data <- cluster_features_scaled[diagnostic_rows, , drop = FALSE]

# Draw the within-cluster sum-of-squares curve used for the elbow assessment.
fviz_nbclust(
  k_diagnostic_data,
  kmeans,
  method = "wss",
  k.max = 8,
  nstart = 25
) +
  labs(
    title = "Elbow method for selecting k",
    x = "Number of clusters (k)",
    y = "Total within-cluster sum of squares"
  ) +
  theme_minimal()

# Draw the mean silhouette curve, where larger values indicate clearer separation.
fviz_nbclust(
  k_diagnostic_data,
  kmeans,
  method = "silhouette",
  k.max = 8,
  nstart = 25
) +
  labs(
    title = "Silhouette method for selecting k",
    x = "Number of clusters (k)",
    y = "Average silhouette width"
  ) +
  theme_minimal()

The elbow curve shows diminishing improvement after a small number of clusters, while the silhouette diagnostic favours a relatively simple solution. This report chooses k = 3 as a transparent compromise between statistical separation and business interpretability: it permits an early/strong group, a typical group, and a late/weak group without creating many small segments that are difficult to communicate. This remains an analytical judgement rather than an objectively true number of clusters.

5 Fit the order-level k-means model

# Store the chosen number of clusters explicitly for auditability.
chosen_k <- 3L

# Reset the seed immediately before fitting the final model.
set.seed(1)

# Fit k-means to every eligible order, using 25 random starts for stability.
kmeans_fit <- kmeans(
  cluster_features_scaled,
  centers = chosen_k,
  nstart = 25
)

# Attach the integer cluster label to its corresponding order-level observation.
clustered_orders <- cluster_orders %>%
  mutate(cluster = factor(kmeans_fit$cluster))

6 Cluster profiles and interpretive names

# Summarise the original, unscaled variables so profiles retain business meaning.
cluster_profile_raw <- clustered_orders %>%
  group_by(cluster) %>%
  summarise(
    n = n(),
    mean_delivery_days = mean(delivery_days),
    mean_delay_days = mean(delay_days),
    mean_review_score = mean(review_score),
    .groups = "drop"
  )

# Identify the relatively worst cluster: latest first, then lowest review score.
worst_cluster_id <- cluster_profile_raw %>%
  arrange(desc(mean_delay_days), mean_review_score) %>%
  slice(1) %>%
  pull(cluster) %>%
  as.character()

# Identify the relatively strongest cluster: earliest first, then highest score.
best_cluster_id <- cluster_profile_raw %>%
  arrange(mean_delay_days, desc(mean_review_score)) %>%
  slice(1) %>%
  pull(cluster) %>%
  as.character()

# Give each cluster a concise relative interpretation based on its profile.
cluster_names <- cluster_profile_raw %>%
  mutate(
    cluster_name = case_when(
      as.character(cluster) == worst_cluster_id ~
        "slower, later, less satisfied",
      as.character(cluster) == best_cluster_id ~
        "faster, earlier, more satisfied",
      TRUE ~ "typical timing, mixed satisfaction"
    ),
    cluster_label = paste0("Cluster ", cluster, ": ", cluster_name)
  ) %>%
  select(cluster, cluster_name, cluster_label)

# Attach the business-readable name to every order.
clustered_orders <- clustered_orders %>%
  left_join(cluster_names, by = "cluster")

# Produce the requested final profile table, rounded only for presentation.
cluster_profile <- cluster_profile_raw %>%
  left_join(cluster_names, by = "cluster") %>%
  select(
    cluster,
    cluster_name,
    n,
    mean_delivery_days,
    mean_delay_days,
    mean_review_score
  ) %>%
  arrange(cluster) %>%
  mutate(
    mean_delivery_days = round(mean_delivery_days, 2),
    mean_delay_days = round(mean_delay_days, 2),
    mean_review_score = round(mean_review_score, 2)
  )

# Display the requested cluster summary in a compact HTML table.
knitr::kable(
  cluster_profile,
  caption = "Order-level cluster profiles in original measurement units",
  col.names = c(
    "Cluster",
    "Interpretive name",
    "Orders (n)",
    "Mean delivery days",
    "Mean delay days",
    "Mean review score"
  )
)

Order-level cluster profiles in original measurement units
Cluster	Interpretive name	Orders (n)	Mean delivery days	Mean delay days	Mean review score
1	typical timing, mixed satisfaction	15760	11.94	-12.25	2.07
2	faster, earlier, more satisfied	74719	10.69	-12.75	4.75
3	slower, later, less satisfied	5991	37.44	11.21	1.98

The labels are relative descriptions of this dataset:

Cluster 2: faster, earlier, more satisfied has the strongest combination of earlier delivery and higher satisfaction.
Cluster 1: typical timing, mixed satisfaction represents the intermediate profile.
Cluster 3: slower, later, less satisfied has the weakest combination of timing and satisfaction.

6.1 Chart 1: cluster profile

# Plot cluster means so the chart represents the profile rather than dense raw points.
ggplot(
  cluster_profile_raw %>% left_join(cluster_names, by = "cluster"),
  aes(
    x = mean_delay_days,
    y = mean_review_score,
    colour = cluster_label,
    size = n
  )
) +
  # Add a vertical reference line at zero: right of zero indicates average lateness.
  geom_vline(xintercept = 0, linetype = "dashed", colour = "grey55") +
  # Draw one point per cluster, sized by its number of orders.
  geom_point(alpha = 0.9) +
  # Label points directly to make the chart understandable without colour alone.
  geom_text(
    aes(label = paste0("Cluster ", cluster)),
    vjust = -1.2,
    show.legend = FALSE
  ) +
  scale_colour_manual(values = cluster_colours[seq_len(chosen_k)]) +
  scale_size_continuous(range = c(6, 13)) +
  labs(
    title = "Cluster profile: delivery timing and customer reviews",
    subtitle = "Points are cluster means; point size represents order count",
    x = "Mean delay versus estimated date (days)",
    y = "Mean review score (1–5)",
    colour = "Order cluster",
    size = "Orders"
  ) +
  theme_minimal()

6.2 Chart 2: distribution of performance against estimate

# Restrict extreme tails in the display so the central distribution remains readable.
delay_plot_limits <- quantile(
  clustered_orders$delay_days,
  probs = c(0.01, 0.99),
  na.rm = TRUE
)

# Compare the delay distribution across the named order clusters.
ggplot(
  clustered_orders,
  aes(x = cluster_label, y = delay_days, fill = cluster_label)
) +
  # Use boxplots because delay_days is skewed and contains meaningful outliers.
  geom_boxplot(outlier.alpha = 0.08, width = 0.65) +
  # Mark the on-time boundary at zero days.
  geom_hline(yintercept = 0, linetype = "dashed", colour = "grey35") +
  # Zoom without discarding data used to calculate the boxplots.
  coord_cartesian(ylim = delay_plot_limits) +
  scale_fill_manual(values = cluster_colours[seq_len(chosen_k)]) +
  labs(
    title = "Delivery performance differs markedly across clusters",
    subtitle = "Negative values are early; positive values are late",
    x = "Order cluster",
    y = "Days delivered relative to estimate",
    fill = "Order cluster"
  ) +
  theme_minimal() +
  theme(
    legend.position = "none",
    axis.text.x = element_text(angle = 20, hjust = 1)
  )

7 Geographic breakdown

States are compared using the share of their orders in each cluster. To reduce unstable rankings based on small samples, the ranking includes only states with at least 200 eligible delivered orders.

# Set and document the minimum number of orders needed for state comparisons.
minimum_state_orders <- 200L

# Count eligible clustered orders for each state and cluster combination.
state_cluster_counts <- clustered_orders %>%
  filter(!is.na(customer_state)) %>%
  count(customer_state, cluster, cluster_label, name = "cluster_orders")

# Calculate the total eligible order count within each state.
state_totals <- state_cluster_counts %>%
  group_by(customer_state) %>%
  summarise(state_orders = sum(cluster_orders), .groups = "drop")

# Convert counts to within-state shares and retain sufficiently large states.
state_cluster_mix <- state_cluster_counts %>%
  left_join(state_totals, by = "customer_state") %>%
  mutate(cluster_share = cluster_orders / state_orders) %>%
  filter(state_orders >= minimum_state_orders)

# Rank eligible states by the share assigned to the relatively worst cluster.
state_ranking <- state_cluster_mix %>%
  filter(as.character(cluster) == worst_cluster_id) %>%
  arrange(desc(cluster_share))

# Preserve that ranking as the state order used in the chart.
state_levels <- state_ranking$customer_state

# Convert state to an ordered factor so the worst-cluster ranking is visible.
state_cluster_mix <- state_cluster_mix %>%
  mutate(customer_state = factor(customer_state, levels = rev(state_levels)))

# Display the five states with the largest share in the relatively worst cluster.
knitr::kable(
  state_ranking %>%
    select(customer_state, state_orders, cluster_share) %>%
    slice_head(n = 5) %>%
    mutate(cluster_share = round(100 * cluster_share, 1)),
  caption = paste0(
    "States with the highest share of the relatively worst cluster (n ≥ ",
    minimum_state_orders,
    " orders)"
  ),
  col.names = c("State", "Eligible orders", "Worst-cluster share (%)")
)

States with the highest share of the relatively worst cluster (n ≥ 200 orders)
State	Eligible orders	Worst-cluster share (%)
AL	397	26.4
MA	717	18.8
SE	335	18.8
PA	946	18.0
CE	1279	15.9

7.1 Chart 3: cluster mix by customer state

# Plot within-state proportions so states of different sizes can be compared fairly.
ggplot(
  state_cluster_mix,
  aes(x = customer_state, y = cluster_share, fill = cluster_label)
) +
  # Stacked bars sum to 100% because cluster_share is calculated within each state.
  geom_col(width = 0.75) +
  coord_flip() +
  # Format proportions as percentages without loading an additional package.
  scale_y_continuous(labels = function(x) paste0(round(100 * x), "%")) +
  scale_fill_manual(values = cluster_colours[seq_len(chosen_k)]) +
  labs(
    title = "Order-cluster mix by customer state",
    subtitle = paste0(
      "States with at least ",
      minimum_state_orders,
      " eligible delivered orders; ranked by worst-cluster share"
    ),
    x = "Customer state",
    y = "Share of the state's eligible orders",
    fill = "Order cluster"
  ) +
  theme_minimal() +
  theme(legend.position = "bottom")

8 Business insights

Timing and satisfaction move together descriptively. The relatively late cluster also has the lowest mean review score, whereas the relatively early cluster has the highest. This is an association within delivered orders and should not be interpreted as proof that lateness caused the lower scores.
Delivery outcomes are heterogeneous rather than uniform. The cluster profiles separate orders with materially different delivery durations and performance against the estimate. Operational teams could use these segments as a starting point for investigating carriers, routes, product types, or fulfilment processes, but the clustering alone does not identify the underlying mechanism.
The share of the relatively worst cluster varies by state. Among states with at least 200 eligible orders, the geographic mix is not identical. This may help prioritise further diagnostic work, but it should not be read as evidence that a state itself causes poorer service.

9 Supervised learning: predicting delivery lateness

9.1 Why supervised learning answers a different question

The k-means analysis asks, “What kinds of delivered orders appear in the data?” It is unsupervised: no outcome label directs the clustering, and the result describes groups with similar delivery and review characteristics.

The random forest instead asks, “Using only information available when an order is placed, can the model predict whether that order will arrive late?” This is supervised learning because the historical late/on-time outcome trains the model. It also ranks predictors by their usefulness for prediction, which can guide further operational investigation.

These are complementary rather than interchangeable questions. Clustering describes realised delivery experiences, whereas supervised learning tests advance prediction and ranks potential warning signals.

9.2 Create a leakage-safe modelling frame

# Build the supervised dataset from the existing delivered-order table.
rf_data <- analysis_orders %>%
  mutate(
    # Define the required binary outcome from realised delivery performance.
    is_late = factor(
      ifelse(delay_days > 0, "late", "on_time"),
      levels = c("on_time", "late")
    ),
    # Measure the delivery time promised when the customer placed the order.
    promised_days = as.numeric(
      difftime(
        order_estimated_delivery_date,
        order_purchase_timestamp,
        units = "days"
      )
    ),
    # Represent seasonal ordering patterns using the purchase month.
    purchase_month = factor(format(order_purchase_timestamp, "%m")),
    # Represent weekly ordering patterns using the purchase weekday.
    purchase_wday = factor(weekdays(order_purchase_timestamp)),
    # Treat state as a categorical location predictor rather than a number.
    customer_state = factor(customer_state)
  ) %>%
  # Remove rows that cannot supply a target or one of the required predictors.
  filter(
    !is.na(is_late),
    !is.na(promised_days),
    !is.na(customer_state)
  ) %>%
  # Keep only the target and information available at order time.
  select(
    is_late,
    customer_state,
    promised_days,
    purchase_month,
    purchase_wday
  )

# Post-delivery variables are deliberately excluded to prevent target leakage.
# In particular, delay_days defines the target, while delivery_days and
# review_score are known only at or after delivery and cannot be valid predictors.

# randomForest permits factor predictors with at most 53 levels.
# customer_state has about 27 levels, so it is safely within that limit.
stopifnot(nlevels(rf_data$customer_state) <= 53)

# Confirm that the modelling frame contains no missing values.
stopifnot(!anyNA(rf_data))

# Calculate the observed class balance for interpretation.
rf_class_balance <- rf_data %>%
  count(is_late, name = "orders") %>%
  mutate(share = orders / sum(orders))

# Display the target distribution before model fitting.
knitr::kable(
  rf_class_balance %>%
    mutate(share = round(100 * share, 1)),
  caption = "Class distribution in the random-forest modelling frame",
  col.names = c("Delivery outcome", "Orders", "Share (%)")
)

Class distribution in the random-forest modelling frame
Delivery outcome	Orders	Share (%)
on_time	88644	91.9
late	7826	8.1

Most delivered orders are early or on time, making late delivery the minority class. Consequently, overall accuracy is not sufficient: a model could predict nearly everything as on time, achieve apparently strong accuracy, and still fail at the business task of identifying late orders. The evaluation therefore reports recall for the "late" class explicitly.

A leaked model using delay_days, delivery_days, or review_score would effectively be given information from or after the event it is supposed to predict. Its accuracy would be meaningless for advance operational use because those values do not exist when a new order is placed.

9.3 Train/test split and model fitting

# Set the seed immediately before sampling to reproduce the same split.
set.seed(1)

# Randomly select 70% of rows for training without replacement.
train_rows <- sample(
  seq_len(nrow(rf_data)),
  size = floor(0.70 * nrow(rf_data)),
  replace = FALSE
)

# Fit the model only on the training observations.
rf_train <- rf_data[train_rows, , drop = FALSE]

# Reserve the remaining 30% as unseen test observations.
rf_test <- rf_data[-train_rows, , drop = FALSE]

# Verify that neither outcome class was lost during the random split.
stopifnot(all(levels(rf_data$is_late) %in% unique(rf_train$is_late)))
stopifnot(all(levels(rf_data$is_late) %in% unique(rf_test$is_late)))

# Reset the seed so the forest itself is reproducible independently of the split.
set.seed(1)

# Train 200 classification trees using only the training set.
rf_model <- randomForest(
  is_late ~ .,
  data = rf_train,
  ntree = 200,
  importance = TRUE
)

9.4 Test-set evaluation

# Generate class predictions for observations not used to train the model.
rf_predictions <- predict(
  rf_model,
  newdata = rf_test,
  type = "class"
)

# Cross-tabulate predicted classes in rows against actual classes in columns.
rf_confusion <- table(
  Predicted = rf_predictions,
  Actual = rf_test$is_late
)

# Display the full confusion matrix so both error directions remain visible.
rf_confusion

##          Actual
## Predicted on_time  late
##   on_time   26337  2247
##   late        208   149

# Calculate the proportion of all test observations classified correctly.
rf_accuracy <- sum(diag(rf_confusion)) / sum(rf_confusion)

# Calculate minority-class recall: correctly predicted late orders / all late orders.
rf_late_recall <- rf_confusion["late", "late"] /
  sum(rf_confusion[, "late"])

# Combine the two headline metrics into a compact audit table.
rf_metrics <- data.frame(
  metric = c("Overall test accuracy", "Recall for late deliveries"),
  value = c(rf_accuracy, rf_late_recall)
)

# Present metrics as percentages while retaining numeric calculations above.
knitr::kable(
  rf_metrics %>%
    mutate(value = paste0(round(100 * value, 1), "%")),
  caption = "Random-forest performance on the held-out 30% test set",
  col.names = c("Metric", "Result")
)

Random-forest performance on the held-out 30% test set
Metric	Result
Overall test accuracy	91.5%
Recall for late deliveries	6.2%

The model achieves 91.5% overall test accuracy, but its recall for late deliveries is 6.2%. The recall figure is the more demanding operational measure: it states the share of genuinely late test orders that the model successfully flags. The difference between accuracy and late-class recall demonstrates why the imbalanced target cannot be judged by accuracy alone.

9.5 Predictor importance and business interpretation

# Draw the random forest's standard permutation and node-purity importance plots.
varImpPlot(
  rf_model,
  main = "Random-forest predictor importance for delivery lateness"
)

# Extract numeric importance values for a transparent ranking in the prose.
rf_importance <- importance(rf_model)

# Use mean decrease in accuracy because it measures predictive loss after shuffling.
rf_importance_ranking <- data.frame(
  predictor = rownames(rf_importance),
  mean_decrease_accuracy = rf_importance[, "MeanDecreaseAccuracy"],
  row.names = NULL
) %>%
  arrange(desc(mean_decrease_accuracy))

# Store the two leading predictors for reproducible inline interpretation.
most_important_predictor <- rf_importance_ranking$predictor[1]
second_important_predictor <- rf_importance_ranking$predictor[2]

# Display the ranking used for the written interpretation.
knitr::kable(
  rf_importance_ranking %>%
    mutate(mean_decrease_accuracy = round(mean_decrease_accuracy, 2)),
  caption = "Predictor importance ranked by mean decrease in accuracy",
  col.names = c("Predictor", "Mean decrease in accuracy")
)

Predictor importance ranked by mean decrease in accuracy
Predictor	Mean decrease in accuracy
customer_state	93.58
promised_days	71.70
purchase_month	47.83
purchase_wday	21.12

The importance ranking places customer_state first and promised_days second. Importance means that predictions become less accurate when a variable’s information is disrupted; it does not mean that the variable causes lateness.

If customer_state ranks highly, the result is consistent with structural logistics differences associated with customer location, including distance from the São Paulo seller hub and carrier-network coverage. This would corroborate the clustering result in which the “slower, later, less satisfied” cluster is more concentrated in particular states: two different methods would then triangulate on geography as a useful delivery-performance signal.

If promised_days ranks highly, it suggests that the length of the quoted delivery window is associated with late-delivery risk. This could reflect systematically optimistic estimates on more difficult routes or products, but the model cannot establish that the estimate itself causes lateness. Month and weekday importance would similarly indicate predictive seasonality or operational timing patterns, not causal calendar effects.

10 Combined conclusion

The two learning paradigms produce one integrated argument. K-means identifies distinct realised delivery experiences, including a relatively slower, later, and less-satisfied segment whose geographic mix varies across states. The random forest then asks whether lateness can be anticipated using only order-time information and shows which available signals contribute most to that prediction. Together, the methods suggest that delivery risk is structured rather than random, while still requiring operational investigation before any causal or policy conclusion is drawn.

11 Limitations

Survivorship bias: the analysis filters to delivered orders only. Undelivered, unavailable, or cancelled orders—arguably including some of the worst delivery experiences—are excluded.
k is assumed, not objectively true: elbow and silhouette diagnostics inform the decision, but k-means does not reveal a single indisputable number of naturally existing customer groups.
Ordinal reviews are treated as continuous: review_score runs from 1 to 5 and is ordinal, yet means, scaling, and Euclidean distance treat the gaps between adjacent scores as equal.
Geography is associational: state patterns may reflect logistics distance from São Paulo, seller location, carrier networks, infrastructure, or product mix rather than regional service quality. The analysis does not establish causation.
Diagnostic sampling: k-selection plots use a reproducible subsample to keep silhouette calculations computationally feasible. A different sample could alter the visual diagnostics slightly.
k-means assumptions: clusters are encouraged to be roughly spherical in scaled Euclidean space and may be influenced by skewness or outliers.
Target leakage handling: post-delivery variables were excluded from the random forest, which protects the validity of advance prediction but also limits the model to the relatively small set of order-time predictors available in these files.
Class imbalance constrains practical use: late deliveries are the minority class, so a model with high accuracy may still have insufficient late-class recall for operational deployment.
Single train/test split: reported performance depends partly on one reproducible 70/30 split. Cross-validation or repeated resampling would give a more stable estimate of performance.
Importance is not causation: random-forest importance ranks predictive contribution within this model; it does not prove that changing a predictor would change delivery outcomes.

12 AI AUDIT NOTE

Codex generated the initial RMarkdown scaffold, including the data-joining, feature-engineering, clustering, visualisation, and reporting structure, and later added the random-forest prediction and evaluation scaffold. Before submission, the student checked or modified the following:

confirmed that the three CSV filenames and columns match the local files;
ran the complete document from a clean R session and confirmed that it knits;
inspected the elbow and silhouette plots and reviewed the justification for choosing k = 3;
checked that duplicate reviews are collapsed before joining and that one delivered order remains one clustering observation;
reviewed the cluster names, state threshold, charts, insights, and limitations for accurate interpretation; and
verified that the random-forest predictors are available at order time, checked the 70/30 split, confusion matrix, accuracy, late-class recall, and variable-importance interpretation;
reviewed the class-imbalance and target-leakage discussion and confirmed that importance is interpreted as predictive association rather than causation; and
confirmed that the student can explain each line, transformation, modelling choice, clustering step, random-forest step, chart, metric, and limitation in their own words.

Delivery Performance in the Brazilian Olist Dataset

Order-level k-means clustering and random-forest prediction of delivery lateness

Student analysis

23 June 2026