Understanding Freight Cost Drivers in Float Glass Export Operations: An Exploratory and Inferential Analytics Study

Author

Hammed Lawal

Published

May 20, 2026

1 Executive Summary

This case study investigates the cost structure of freight operations for float glass exports from Nigerian seaports. The data covers 100 shipments booked between January and August 2024 from three port origins — Lekki, Apapa, and Tincan — dispatched to destinations across Sub-Saharan Africa, North Africa and the Mediterranean, the Middle East, and the Americas. The central business problem is: what drives freight cost, and how can operations management use this understanding to price shipments, benchmark routes, and optimise tonnage loading decisions?

Using five analytical techniques — Exploratory Data Analysis, Data Visualisation, Hypothesis Testing, Correlation Analysis, and Linear Regression — this study finds that tonnage is the single strongest predictor of total freight cost (r = 0.96), and that destination region is a highly significant cost driver even after controlling for cargo volume. Shipments to North Africa and the Mediterranean attract rates roughly three times higher per tonne than those to Sub-Saharan Africa, while the Americas cluster at an intermediate level. A two-variable regression model (tonnage + region) explains 97% of freight cost variance, providing a practical pricing formula for the commercial team.

2 Professional Disclosure

Job Title: Freight & Logistics Coordinator — Export Operations

Organisation Domain: The organisation is a freight forwarding and logistics company operating within Nigeria’s export supply chain. It manages the end-to-end shipment of manufactured goods — primarily flat and float glass — from Nigerian seaports to international buyers across Africa, the Mediterranean basin, the Middle East, and the Americas. Day-to-day work involves booking ocean freight, coordinating with shipping lines, tracking vessel schedules, and generating shipping instructions and commercial documentation.

Relevance of the Five Techniques to Daily Work:

Exploratory Data Analysis (EDA): Before any pricing or routing decision, I need to understand the spread and central tendency of freight rates in our booking records. EDA allows me to quickly identify outliers — shipments priced unusually high or low — and flag data-entry errors in the booking system before they propagate into financial reports.
Data Visualisation: Management and clients require visual summaries of freight trends across routes and periods. Visualisations allow me to communicate monthly freight spend, cost per tonne by destination, and seasonal patterns in a format accessible to non-technical stakeholders such as sales managers and CFO.
Hypothesis Testing: A recurring operational question is whether different loading ports (Lekki vs Apapa) yield meaningfully different freight costs. Formal hypothesis testing replaces anecdotal reasoning with statistically defensible answers, informing decisions about which port to assign future bookings.
Correlation Analysis: Understanding how freight cost relates to tonnage, dwell time at port, and transit days allows the team to build an intuition for which variables to monitor. For example, if transit days and freight are highly correlated, route selection becomes a cost lever, not merely a scheduling preference.
Linear Regression: A regression model formalises the relationship between freight cost and its predictors, giving the commercial team a pricing equation. When quoting a new shipment, the model provides an evidence-based estimate rather than relying solely on the operator’s experience or manual lookup tables.

3 Data Collection & Sampling

3.1 Source and Collection Method

The dataset was extracted from the company’s internal freight booking management system (FMS), which logs every shipment booking from initiation to vessel departure. Each record is created when a booking confirmation is received from the shipping line and is updated as the shipment progresses. The export was performed as a CSV extract covering all float glass bookings within the 2024 calendar year up to the date of extraction (late July 2024).

3.2 Sampling Frame and Period

Population: All ocean freight bookings involving float glass commodity from Nigerian ports managed by the organisation.
Sample frame: Bookings with a confirmed ship date between 1 January 2024 and 31 July 2024.
Sample size: 100 booking records (the complete population for the period — this is a census, not a random sample).
Variables captured: Booking reference, ship date, commodity, port of loading (POL), port of discharge (POD), freight amount (USD), cargo tonnage, vessel dwell time at port of loading (stay), and ocean transit time.

4 Data Description

4.1 Variable Dictionary

Variable	Original Name	Type	Description
`booking_no`	BOOKING NO	Character	Unique freight booking reference
`ship_date`	SHIP DATE	Date	Date cargo departed port of loading
`pol`	POL	Categorical	Port of loading (LEKKI, APAPA, TINCAN)
`pod`	POD	Categorical	Port of discharge (25 unique destinations)
`region`	— (derived)	Categorical	Geographic region of destination
`freight_usd`	FREIGHT USD	Numeric	Total ocean freight in US Dollars
`tonnage`	TONNAGE	Numeric	Cargo weight in metric tonnes
`stay_days`	STAY	Numeric	Vessel dwell time at POL (days)
`transit_days`	TRANSIT TIME	Numeric	Ocean transit time POL → POD (days)
`freight_rate_per_ton`	— (derived)	Numeric	Freight USD ÷ Tonnage (USD/MT)

Note: COMMODITY (all records = “FLOAT GLASS”) and S/N (row index) were excluded from analysis as they carry no analytical information.

4.2 Data Loading and Cleaning

Code

# ── Libraries ──────────────────────────────────────────────────────────────────
library(tidyverse)
library(lubridate)
library(janitor)
library(knitr)
library(kableExtra)
library(scales)
library(corrplot)
library(broom)
library(car)
library(rstatix)
library(ggthemes)
library(patchwork)
library(nortest)

# ── Load raw data ──────────────────────────────────────────────────────────────
df_raw <- read_csv("data.csv", show_col_types = FALSE)

# ── Clean column names ─────────────────────────────────────────────────────────
df_raw <- df_raw %>% clean_names()

# ── Step 1: Strip invisible Unicode characters from all character columns ──────
# Some dates contain Unicode left-to-right marks (U+202A, U+202C) and
# other non-printable bytes — a known artefact of copy-pasting from certain
# logistics portals into Excel before CSV export.
# iconv converts each string to UTF-8, replacing unmappable bytes with "".
# Then we strip any remaining non-ASCII (codes > 127) and trim whitespace.
df_raw <- df_raw %>%
  mutate(across(
    where(is.character),
    ~ iconv(., from = "UTF-8", to = "ASCII", sub = "") %>%
      str_trim()
  ))

# ── Step 2: Parse ship_date (multiple date formats present) ───────────────────
df_raw <- df_raw %>%
  mutate(ship_date = parse_date_time(ship_date,
                                     orders = c("m/d/Y", "d/m/Y", "m/d/y"),
                                     quiet  = TRUE))

# Correct one clear year-entry error: row 13 records "2/25/25" — contextually
# this must be 25 Feb 2024 given the surrounding entries.
df_raw <- df_raw %>%
  mutate(ship_date = if_else(
    !is.na(ship_date) & year(ship_date) == 2025,
    ship_date - years(1),
    ship_date
  ))

# ── Step 3: Parse stay_days and transit_days (stored as "6DAYS", "15DAY") ─────
df_raw <- df_raw %>%
  mutate(
    stay_days    = as.numeric(gsub("[A-Za-z]", "", stay)),
    transit_days = as.numeric(gsub("[A-Za-z]", "", transit_time))
  )

# ── Step 4: Parse freight_usd (commas in numeric strings) ─────────────────────
df_raw <- df_raw %>%
  mutate(freight_usd = as.numeric(gsub(",", "", freight_usd)))

# ── Step 5: Standardise POD spelling variants ─────────────────────────────────
# "BUEVENATURA" and "BUENAVENTURA" refer to the same Colombian port.
# "Iquique,Chile" and "Iquique" refer to the same Chilean port.
df_raw <- df_raw %>%
  mutate(pod = case_when(
    str_detect(pod, regex("BUEV|BUEN", ignore_case = TRUE)) ~ "BUENAVENTURA",
    str_detect(pod, regex("Iquique", ignore_case = TRUE))   ~ "IQUIQUE",
    str_detect(pod, regex("Guayaquil", ignore_case = TRUE)) ~ "GUAYAQUIL",
    TRUE ~ toupper(str_trim(pod))
  ))

# ── Step 6: Assign destination regions ────────────────────────────────────────
df_raw <- df_raw %>%
  mutate(region = case_when(
    pod %in% c("ABIDJAN", "MONROVIA", "POINTE NOIRE",
               "DAR ES SALAM")                                    ~ "Sub-Saharan Africa",
    pod %in% c("SKIKDA", "BEJAIA", "ALGER", "DJEN DJEN",
               "DAMIETTA", "EL KHOMS", "MISURATA")               ~ "N. Africa / Mediterranean",
    pod %in% c("HAMAD")                                           ~ "Middle East",
    pod %in% c("CALLAO", "CAUCEDO", "GUAYAQUIL", "CARTAGENA",
               "BUENAVENTURA", "BARRANQUILA", "PUERTO CORTES",
               "ARICA", "IQUIQUE", "PUERTO CABELLO")             ~ "Americas",
    TRUE ~ "Other"
  ))

# ── Step 7: Derived variables ─────────────────────────────────────────────────
df_raw <- df_raw %>%
  mutate(
    freight_rate_per_ton = round(freight_usd / tonnage, 2),
    ship_month           = floor_date(ship_date, "month"),
    ship_month_label     = format(ship_date, "%b %Y")
  )

# ── Step 8: Final analytical dataset ──────────────────────────────────────────
df <- df_raw %>%
  select(booking_no, ship_date, ship_month, ship_month_label,
         pol, pod, region,
         freight_usd, tonnage, stay_days, transit_days,
         freight_rate_per_ton) %>%
  filter(!is.na(freight_usd))    # drop any rows where freight failed to parse

cat("Final dataset dimensions:", nrow(df), "rows ×", ncol(df), "columns\n")

Final dataset dimensions: 100 rows × 12 columns

4.3 Summary Statistics

Code

# Numeric summary
num_summary <- df %>%
  select(freight_usd, tonnage, stay_days, transit_days, freight_rate_per_ton) %>%
  pivot_longer(everything(), names_to = "variable", values_to = "value") %>%
  group_by(variable) %>%
  summarise(
    n       = sum(!is.na(value)),
    mean    = round(mean(value, na.rm = TRUE), 2),
    median  = round(median(value, na.rm = TRUE), 2),
    sd      = round(sd(value, na.rm = TRUE), 2),
    min     = round(min(value, na.rm = TRUE), 2),
    max     = round(max(value, na.rm = TRUE), 2),
    missing = sum(is.na(value))
  )

num_summary %>%
  kbl(caption = "Summary Statistics — Numeric Variables",
      col.names = c("Variable", "N", "Mean", "Median", "SD", "Min", "Max", "Missing")) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = FALSE)

Summary Statistics — Numeric Variables
Variable	N	Mean	Median	SD	Min	Max
freight_rate_per_ton	100	39.92	37.74	12.57	6.6	86.79
freight_usd	100	4471.48	2672.00	6070.63	500.0	35437.30
stay_days	100	11.24	10.00	3.65	5.0	25.00
tonnage	100	105.95	79.50	99.09	26.5	530.00
transit_days	100	42.31	45.00	19.23	3.0	67.00

Code

import pandas as pd
import numpy as np

df_py = pd.read_csv("data.csv", encoding="latin-1")

# --- clean column names ---
df_py.columns = (df_py.columns
                 .str.strip()
                 .str.lower()
                 .str.replace(r'\s+', '_', regex=True)
                 .str.replace(r'[^\w]', '_', regex=True))

# --- strip non-ASCII from string columns ---
for col in df_py.select_dtypes(include='object').columns:
    df_py[col] = df_py[col].astype(str).str.replace(r'[^\x20-\x7E]', '', regex=True).str.strip()

# --- parse freight_usd ---
df_py['freight_usd'] = pd.to_numeric(df_py['freight_usd'].str.replace(',', '', regex=False), errors='coerce')

# --- parse tonnage ---
df_py['tonnage'] = pd.to_numeric(df_py['tonnage'], errors='coerce')

# --- parse stay_days and transit_days ---
df_py['stay_days']    = pd.to_numeric(df_py['stay'].str.replace(r'[A-Za-z]', '', regex=True), errors='coerce')
df_py['transit_days'] = pd.to_numeric(df_py['transit_time'].str.replace(r'[A-Za-z]', '', regex=True), errors='coerce')

# --- derived ---
df_py['freight_rate_per_ton'] = df_py['freight_usd'] / df_py['tonnage']

num_cols = ['freight_usd', 'tonnage', 'stay_days', 'transit_days', 'freight_rate_per_ton']
summary = df_py[num_cols].describe().T
summary['missing'] = df_py[num_cols].isna().sum()
print(summary.round(2).to_string())

                      count     mean      std    min      25%      50%      75%       max  missing
freight_usd           100.0  4471.48  6070.63  500.0  1387.50  2672.00  6000.00  35437.30        0
tonnage               100.0   105.95    99.09   26.5    53.00    79.50   132.50    530.00        0
stay_days             100.0    11.24     3.65    5.0    10.00    10.00    15.00     25.00        0
transit_days          100.0    42.31    19.23    3.0    30.00    45.00    60.00     67.00        0
freight_rate_per_ton  100.0    39.92    12.57    6.6    34.53    37.74    49.33     86.79        0

5 Technique 1: Exploratory Data Analysis (EDA)

5.1 Theoretical Background

Exploratory Data Analysis, formalised by Tukey (1977), is the practice of summarising and visualising a dataset’s core structure before any formal modelling. It encompasses four pillars: summary statistics (central tendency, dispersion), distribution analysis (skewness, modality), missing-value and outlier detection, and bivariate relationships. EDA forces the analyst to question data quality and to develop intuitions about which transformations or model structures are appropriate (James et al., 2021, Ch. 4).

5.2 Business Justification

In freight operations, undetected data errors — miskeyed amounts, wrong date formats, duplicate bookings — directly distort financial reports and route-profitability analyses. Before presenting any cost figures to management, the operations team must verify that the underlying booking data is complete, internally consistent, and free of anomalous entries. EDA is the first line of defence against garbage-in/garbage-out reporting.

5.3 Analysis

5.3.1 Missing Values and Data Quality Issues

Code

# ── Missing value map ──────────────────────────────────────────────────────────
miss_df <- df %>%
  summarise(across(everything(), ~ sum(is.na(.)))) %>%
  pivot_longer(everything(), names_to = "variable", values_to = "n_missing") %>%
  mutate(pct_missing = round(n_missing / nrow(df) * 100, 1))

miss_df %>%
  filter(n_missing > 0) %>%
  kbl(caption = "Variables with Missing Values",
      col.names = c("Variable", "N Missing", "% Missing")) %>%
  kable_styling(bootstrap_options = c("striped", "condensed"), full_width = FALSE)

Variables with Missing Values
Variable	N Missing	% Missing

Code

if (nrow(filter(miss_df, n_missing > 0)) == 0) {
  cat("No missing values remain in the analytical dataset after cleaning.\n")
}

No missing values remain in the analytical dataset after cleaning.

Code

# ── Data quality issues identified ────────────────────────────────────────────

# Issue 1: Date format inconsistency
# Raw data contains at least four date formats (d/m/Y, m/d/Y, m/d/y, and
# dates with embedded Unicode control characters). After cleaning, one record
# had a year of 2025 — almost certainly a keying error for 2024.

n_date_problems <- df_raw %>% filter(is.na(ship_date)) %>% nrow()
cat("Issue 1 — Malformed date strings: ", n_date_problems, "record(s) could not be parsed.\n")

Issue 1 — Malformed date strings:  0 record(s) could not be parsed.

Code

# Issue 2: STAY and TRANSIT TIME as text strings
cat("Issue 2 — Numeric fields stored as text: STAY and TRANSIT TIME required ",
    "regex extraction before use.\n")

Issue 2 — Numeric fields stored as text: STAY and TRANSIT TIME required  regex extraction before use.

Code

# Issue 3: FREIGHT_USD outliers (very large DJEN DJEN shipments)
q_upper <- quantile(df$freight_usd, 0.75) + 1.5 * IQR(df$freight_usd)
outliers <- df %>% filter(freight_usd > q_upper)
cat("Issue 3 — Upper-fence outliers in FREIGHT_USD (IQR rule): ", nrow(outliers),
    " record(s). Max = $", format(max(df$freight_usd), big.mark = ","), "\n")

Issue 3 — Upper-fence outliers in FREIGHT_USD (IQR rule):  3  record(s). Max = $ 35,437.3

Code

# Issue 4: POD spelling variants
cat("Issue 4 — Spelling variants in POD: 'BUEVENATURA'/'BUENAVENTURA' and ",
    "'Iquique,Chile'/'Iquique' were standardised.\n")

Issue 4 — Spelling variants in POD: 'BUEVENATURA'/'BUENAVENTURA' and  'Iquique,Chile'/'Iquique' were standardised.

5.3.2 Distributions of Key Numeric Variables

Code

p_freight <- ggplot(df, aes(x = freight_usd)) +
  geom_histogram(binwidth = 2000, fill = "#2C7BB6", colour = "white", alpha = 0.85) +
  geom_vline(xintercept = mean(df$freight_usd), colour = "#D7191C",
             linetype = "dashed", linewidth = 0.8) +
  geom_vline(xintercept = median(df$freight_usd), colour = "#1A9641",
             linetype = "dotted", linewidth = 0.8) +
  annotate("text", x = mean(df$freight_usd) + 1200, y = 18,
           label = paste0("Mean\n$", round(mean(df$freight_usd)/1000,1), "k"),
           colour = "#D7191C", size = 3) +
  annotate("text", x = median(df$freight_usd) - 1800, y = 18,
           label = paste0("Median\n$", round(median(df$freight_usd)/1000,1), "k"),
           colour = "#1A9641", size = 3) +
  scale_x_continuous(labels = label_dollar(scale = 1e-3, suffix = "k")) +
  labs(title = "Distribution of Freight USD",
       x = "Freight (USD)", y = "Count") +
  theme_minimal(base_size = 12)

p_tonnage <- ggplot(df, aes(x = tonnage)) +
  geom_histogram(binwidth = 25, fill = "#FC8D59", colour = "white", alpha = 0.85) +
  labs(title = "Distribution of Tonnage (MT)",
       x = "Tonnage (metric tonnes)", y = "Count") +
  theme_minimal(base_size = 12)

p_stay <- ggplot(df, aes(x = stay_days)) +
  geom_histogram(binwidth = 2, fill = "#91CF60", colour = "white", alpha = 0.85) +
  labs(title = "Dwell Time at POL (days)",
       x = "Stay (days)", y = "Count") +
  theme_minimal(base_size = 12)

p_transit <- ggplot(df, aes(x = transit_days)) +
  geom_histogram(binwidth = 5, fill = "#762A83", colour = "white", alpha = 0.85) +
  labs(title = "Ocean Transit Time (days)",
       x = "Transit (days)", y = "Count") +
  theme_minimal(base_size = 12)

(p_freight | p_tonnage) / (p_stay | p_transit) +
  plot_annotation(
    title = "Figure 1: Distributions of Key Numeric Variables",
    theme = theme(plot.title = element_text(size = 14, face = "bold"))
  )

Distribution of key numeric variables. FREIGHT_USD shows strong positive skew driven by large DJEN DJEN shipments.

5.3.3 Outlier Detection — Box Plots

Code

df_long_num <- df %>%
  select(freight_usd, tonnage, stay_days, transit_days) %>%
  pivot_longer(everything(), names_to = "variable", values_to = "value") %>%
  mutate(variable = case_match(variable,
    "freight_usd"  ~ "Freight (USD)",
    "tonnage"      ~ "Tonnage (MT)",
    "stay_days"    ~ "Stay (days)",
    "transit_days" ~ "Transit (days)",
    .default = variable
  ))

ggplot(df_long_num, aes(x = variable, y = value, fill = variable)) +
  geom_boxplot(outlier.colour = "#D7191C", outlier.shape = 16,
               outlier.size = 2.5, alpha = 0.7) +
  scale_fill_brewer(palette = "Set2") +
  facet_wrap(~variable, scales = "free", nrow = 1) +
  labs(title = "Figure 2: Box Plots — Outlier Detection",
       x = NULL, y = "Value") +
  theme_minimal(base_size = 12) +
  theme(legend.position = "none",
        strip.text = element_text(face = "bold"))

Box plots reveal that FREIGHT_USD and TONNAGE contain notable upper outliers (the three large DJEN DJEN shipments at 530 MT each).

5.3.4 Categorical Variable Frequencies

Code

p_pol <- df %>%
  count(pol) %>%
  ggplot(aes(x = reorder(pol, n), y = n, fill = pol)) +
  geom_col(alpha = 0.85) +
  geom_text(aes(label = n), hjust = -0.2, size = 3.5) +
  coord_flip() +
  scale_fill_manual(values = c("#2C7BB6", "#FC8D59", "#1A9641")) +
  labs(title = "Bookings by Port of Loading", x = NULL, y = "Count") +
  theme_minimal(base_size = 11) +
  theme(legend.position = "none") +
  ylim(0, max(count(df, pol)$n) * 1.15)

p_region <- df %>%
  count(region) %>%
  ggplot(aes(x = reorder(region, n), y = n, fill = region)) +
  geom_col(alpha = 0.85) +
  geom_text(aes(label = n), hjust = -0.2, size = 3.5) +
  coord_flip() +
  scale_fill_brewer(palette = "Set1") +
  labs(title = "Bookings by Destination Region", x = NULL, y = "Count") +
  theme_minimal(base_size = 11) +
  theme(legend.position = "none") +
  ylim(0, max(count(df, region)$n) * 1.15)

p_pol | p_region

5.4 Interpretation

For a non-technical manager: The freight cost data are right-skewed — most shipments cost between $1,000 and $10,000, but three very large shipments to DJEN DJEN, Algeria (each carrying 530 metric tonnes) pushed the average cost to $4,471 versus a median of $2,672. This gap between mean and median is a red flag that averages alone should not be used in route-profitability reporting without flagging the outsized DJEN DJEN effect.

Four data quality issues were identified and resolved before analysis: (1) inconsistent date formats with embedded invisible characters, (2) numeric fields encoded as text, (3) a year-entry error (2025 for a date that was clearly 2024), and (4) spelling variants for the same port. These are common artefacts of manual data entry in logistics systems and underscore the need for input validation at the booking stage.

Lekki dominates as the loading port (76 of 100 bookings), with Apapa accounting for 22 and Tincan for 2. The Americas and North Africa/Mediterranean account for the bulk of destinations by shipment count.

6 Technique 2: Data Visualisation

6.1 Theoretical Background

Data visualisation translates quantitative patterns into perceptual signals — position, length, colour, shape — that the human visual system processes preattentively (Ware, 2012). The Grammar of Graphics framework (Wilkinson, 2005; Wickham, 2016) organises visual communication into layers: data, aesthetic mappings, geometric objects, scales, and coordinate systems. Effective chart selection requires matching the geometric object to the data type and the question: scatter plots for continuous-continuous relationships, bar charts for categorical comparisons, line plots for temporal trends (Evergreen, 2017, as cited in the module textbook Ch. 5).

6.2 Business Justification

The operations team currently reports freight costs as a single monthly total. This obscures which routes are most expensive, whether costs are trending upward, and whether different loading ports achieve different results. A visualisation narrative transforms the booking ledger into a story about route profitability and cost drivers — information that directly supports pricing negotiations with shipping lines at the next contract review.

6.3 Visualisation Narrative: Five Plots That Tell One Story

6.3.1 Plot 1 — Freight Cost by Destination Region

Code

region_summary <- df %>%
  group_by(region) %>%
  summarise(
    n_shipments  = n(),
    total_freight = sum(freight_usd),
    avg_freight   = mean(freight_usd),
    avg_rate_pt   = mean(freight_rate_per_ton)
  ) %>%
  arrange(desc(avg_freight))

p_region_avg <- ggplot(region_summary,
       aes(x = reorder(region, avg_freight), y = avg_freight, fill = region)) +
  geom_col(alpha = 0.88) +
  geom_text(aes(label = paste0("$", format(round(avg_freight), big.mark = ","))),
            hjust = -0.1, size = 3.5, fontface = "bold") +
  coord_flip() +
  scale_y_continuous(labels = label_dollar(scale = 1e-3, suffix = "k"),
                     expand = expansion(mult = c(0, 0.18))) +
  scale_fill_brewer(palette = "Set1") +
  labs(title = "Average Freight Cost by Destination Region",
       subtitle = "N. Africa / Med carries the highest average freight bill per booking",
       x = NULL, y = "Average Freight (USD)") +
  theme_minimal(base_size = 12) +
  theme(legend.position = "none")

p_region_rate <- ggplot(region_summary,
       aes(x = reorder(region, avg_rate_pt), y = avg_rate_pt, fill = region)) +
  geom_col(alpha = 0.88) +
  geom_text(aes(label = paste0("$", round(avg_rate_pt, 1), "/MT")),
            hjust = -0.1, size = 3.5, fontface = "bold") +
  coord_flip() +
  scale_y_continuous(expand = expansion(mult = c(0, 0.2))) +
  scale_fill_brewer(palette = "Set1") +
  labs(title = "Average Freight Rate per Metric Tonne",
       subtitle = "N. Africa / Med commands ~3× the Sub-Saharan Africa rate",
       x = NULL, y = "Freight Rate (USD / MT)") +
  theme_minimal(base_size = 12) +
  theme(legend.position = "none")

p_region_avg / p_region_rate +
  plot_annotation(
    title = "Figure 3: Freight Economics by Destination Region",
    theme = theme(plot.title = element_text(size = 14, face = "bold"))
  )

Figure 3: Total and per-shipment freight cost vary significantly by destination region.

6.3.2 Plot 2 — Tonnage vs Freight (The Core Relationship)

Code

ggplot(df, aes(x = tonnage, y = freight_usd, colour = region, shape = pol)) +
  geom_point(size = 3, alpha = 0.80) +
  geom_smooth(aes(group = region), method = "lm", se = FALSE,
              linewidth = 0.9, linetype = "dashed") +
  scale_colour_brewer(palette = "Set1", name = "Destination Region") +
  scale_shape_manual(values = c(16, 17, 15), name = "Port of Loading") +
  scale_x_continuous(breaks = seq(0, 600, 100)) +
  scale_y_continuous(labels = label_dollar(scale = 1e-3, suffix = "k")) +
  labs(
    title = "Figure 4: Tonnage vs Freight Cost, Coloured by Region",
    subtitle = "Within each region, freight scales linearly with tonnage — but the rate per MT varies dramatically",
    x = "Tonnage (metric tonnes)",
    y = "Freight (USD)"
  ) +
  theme_minimal(base_size = 12) +
  theme(legend.position = "right")

Figure 4: Freight cost is almost perfectly linearly proportional to tonnage within each destination region, but the slope (rate per MT) differs sharply by region.

6.3.3 Plot 3 — Monthly Freight Trend

Code

monthly_trend <- df %>%
  group_by(ship_month) %>%
  summarise(
    total_freight = sum(freight_usd),
    n_shipments   = n(),
    avg_freight   = mean(freight_usd)
  )

ggplot(monthly_trend, aes(x = ship_month)) +
  geom_col(aes(y = total_freight), fill = "#2C7BB6", alpha = 0.75, width = 20) +
  geom_line(aes(y = avg_freight * 4), colour = "#D7191C",
            linewidth = 1.1, linetype = "solid") +
  geom_point(aes(y = avg_freight * 4), colour = "#D7191C", size = 3) +
  geom_text(aes(y = total_freight + 2000, label = paste0("n=", n_shipments)),
            size = 3, colour = "grey30") +
  scale_y_continuous(
    name    = "Total Monthly Freight (USD)",
    labels  = label_dollar(scale = 1e-3, suffix = "k"),
    sec.axis = sec_axis(~ . / 4,
                        name   = "Avg. Freight per Booking (USD)",
                        labels = label_dollar(scale = 1e-3, suffix = "k"))
  ) +
  scale_x_date(date_labels = "%b %Y", date_breaks = "1 month") +
  labs(
    title    = "Figure 5: Monthly Freight Spend (Jan–Aug 2024)",
    subtitle = "Bars = total spend; red line = average freight per booking. June was the costliest month.",
    x        = "Month"
  ) +
  theme_minimal(base_size = 12) +
  theme(axis.text.x = element_text(angle = 30, hjust = 1))

Figure 5: Total monthly freight spend fluctuates between January and July 2024.

6.3.4 Plot 4 — Freight Rate per Tonne by Destination (Top 10 Ports)

Code

top_pods <- df %>%
  group_by(pod) %>%
  summarise(avg_rate = mean(freight_rate_per_ton), n = n()) %>%
  filter(n >= 2) %>%                         # exclude one-off ports
  arrange(desc(avg_rate)) %>%
  slice_head(n = 12)

ggplot(top_pods, aes(x = reorder(pod, avg_rate), y = avg_rate, fill = avg_rate)) +
  geom_col(alpha = 0.88) +
  geom_text(aes(label = paste0("$", round(avg_rate, 0), "/MT")),
            hjust = -0.1, size = 3.2) +
  coord_flip() +
  scale_fill_gradient(low = "#91CF60", high = "#D7191C",
                      name = "USD / MT") +
  scale_y_continuous(expand = expansion(mult = c(0, 0.18))) +
  labs(
    title    = "Figure 6: Freight Rate per Metric Tonne by Port of Discharge",
    subtitle = "Ports with ≥ 2 bookings shown. DJEN DJEN is the most expensive route; Abidjan the cheapest.",
    x        = NULL,
    y        = "Average Freight Rate (USD / MT)"
  ) +
  theme_minimal(base_size = 12)

Figure 6: DJEN DJEN commands the highest per-tonne rate, more than double SKIKDA. West African ports are the cheapest to serve.

6.3.5 Plot 5 — Loading Port Performance Comparison

Code

pol_summary <- df %>%
  group_by(pol) %>%
  summarise(
    n_shipments    = n(),
    avg_tonnage    = mean(tonnage),
    avg_freight    = mean(freight_usd),
    avg_rate_pt    = mean(freight_rate_per_ton),
    total_freight  = sum(freight_usd)
  )

p_pol_tonnage <- ggplot(pol_summary,
  aes(x = pol, y = avg_tonnage, fill = pol)) +
  geom_col(alpha = 0.85, width = 0.5) +
  geom_text(aes(label = paste0(round(avg_tonnage, 0), " MT")),
            vjust = -0.5, size = 4, fontface = "bold") +
  scale_fill_manual(values = c("#2C7BB6", "#FC8D59", "#1A9641")) +
  labs(title = "Avg. Tonnage per Booking",
       x = "Port of Loading", y = "MT") +
  theme_minimal(base_size = 11) +
  theme(legend.position = "none")

p_pol_rate <- ggplot(pol_summary,
  aes(x = pol, y = avg_rate_pt, fill = pol)) +
  geom_col(alpha = 0.85, width = 0.5) +
  geom_text(aes(label = paste0("$", round(avg_rate_pt, 0), "/MT")),
            vjust = -0.5, size = 4, fontface = "bold") +
  scale_fill_manual(values = c("#2C7BB6", "#FC8D59", "#1A9641")) +
  labs(title = "Avg. Freight Rate (USD/MT)",
       x = "Port of Loading", y = "USD / MT") +
  theme_minimal(base_size = 11) +
  theme(legend.position = "none")

p_pol_tonnage | p_pol_rate

Figure 7: Apapa handles disproportionately larger cargo loads than Lekki, resulting in higher total freight despite fewer bookings.

6.4 Interpretation

For a non-technical manager: Five key insights emerge from the visualisation narrative:

North Africa / Mediterranean routes are the most expensive per tonne — up to three times the cost of serving West African ports like Abidjan. When quoting new business to Algerian buyers, the freight component alone makes the deal structurally different.
Freight cost is almost entirely determined by tonnage and destination — the scatter plot (Figure 4) shows near-perfect straight lines within each region. This means pricing is predictable and should not require ad-hoc negotiation.
June 2024 was the costliest month, driven by a cluster of large bookings to Algeria and South America. Understanding volume seasonality helps negotiate forward freight agreements with shipping lines.
DJEN DJEN (Algeria) charges the highest per-tonne rate ($67 /MT) — more than double SKIKDA. If alternative Algerian port calls exist, benchmarking their cost is a priority.
Apapa handles larger average cargo sizes than Lekki (123 MT vs 102 MT), which inflates its total freight spend. However, the per-tonne rate is comparable between the two ports.

7 Technique 3: Hypothesis Testing

7.1 Theoretical Background

Hypothesis testing is a formal framework for deciding, under uncertainty, whether an observed difference between groups is likely to reflect a real population-level effect or could plausibly be attributed to sampling noise. The procedure specifies a null hypothesis (H₀) — the default position of no effect — and an alternative (H₁). A test statistic is computed and compared to a reference distribution. The resulting p-value measures the probability of observing a result at least as extreme as the one obtained, given that H₀ is true. A p-value below the pre-specified significance level α (here, α = 0.05) leads to rejection of H₀ (Spiegel & Stephens, 2018, as cited in Ch. 6). Effect sizes (Cohen’s d, η²) quantify the practical magnitude of any detected difference, independent of sample size (Ellis, 2010).

7.2 Business Justification

The operations team regularly debates whether routing cargo through Lekki versus Apapa affects freight cost. A hypothesis test turns this debate into a defensible, data-driven decision. Similarly, management needs to know whether destination region significantly explains freight variability — if it does, region should be a first-order variable in cost budgets; if it does not, tonnage alone may suffice.

7.3 Test 1: Do Freight Costs Differ Between Lekki and Apapa?

H₀: The mean freight cost for bookings from Lekki equals the mean for Apapa (μ_Lekki = μ_Apapa). H₁: The mean freight costs are different (μ_Lekki ≠ μ_Apapa). Test: Independent two-sample t-test (or Wilcoxon rank-sum test if normality assumption is violated).

Code

# Restrict to LEKKI and APAPA (TINCAN has only 3 observations)
df_ttest <- df %>% filter(pol %in% c("LEKKI", "APAPA"))

# Step 1: Check normality within each group (Shapiro-Wilk)
sw_results <- df_ttest %>%
  group_by(pol) %>%
  summarise(
    n       = n(),
    mean    = round(mean(freight_usd), 0),
    sd      = round(sd(freight_usd), 0),
    sw_stat = round(shapiro.test(freight_usd)$statistic, 4),
    sw_p    = round(shapiro.test(freight_usd)$p.value, 4)
  )

sw_results %>%
  kbl(caption = "Shapiro-Wilk Normality Test by Port of Loading",
      col.names = c("POL", "N", "Mean Freight", "SD", "W Statistic", "p-value")) %>%
  kable_styling(bootstrap_options = c("striped", "condensed"), full_width = FALSE)

Shapiro-Wilk Normality Test by Port of Loading
POL	N	Mean Freight	SD	W Statistic	p-value
APAPA	22	5712	9740	0.4615	0
LEKKI	76	4190	4616	0.6092	0

Code

# Both groups have p < 0.05 in Shapiro-Wilk → non-normal → use Wilcoxon test
wilcox_result <- wilcox.test(freight_usd ~ pol, data = df_ttest,
                              exact = FALSE, conf.int = TRUE)

# Effect size: rank-biserial correlation
eff_size <- df_ttest %>%
  wilcox_effsize(freight_usd ~ pol)

cat("Wilcoxon Rank-Sum Test: Lekki vs Apapa Freight\n")

Wilcoxon Rank-Sum Test: Lekki vs Apapa Freight

Code

cat("----------------------------------------------\n")

----------------------------------------------

Code

cat("W statistic : ", wilcox_result$statistic, "\n")

W statistic :  811.5

Code

cat("p-value     : ", round(wilcox_result$p.value, 4), "\n")

p-value     :  0.838

Code

cat("95% CI diff : [", round(wilcox_result$conf.int[1], 0),
    ",", round(wilcox_result$conf.int[2], 0), "]\n")

95% CI diff : [ -1100 , 700 ]

Code

cat("Effect size (r) : ", round(eff_size$effsize, 3),
    "(", eff_size$magnitude, ")\n")

Effect size (r) :  0.021 ( 1 )

Code

ggplot(df_ttest, aes(x = pol, y = freight_usd, fill = pol)) +
  geom_violin(alpha = 0.4, trim = FALSE) +
  geom_boxplot(width = 0.25, alpha = 0.7, outlier.colour = "red") +
  stat_summary(fun = mean, geom = "point", shape = 23,
               size = 4, fill = "white", colour = "black") +
  scale_fill_manual(values = c("#2C7BB6", "#FC8D59")) +
  scale_y_continuous(labels = label_dollar(scale = 1e-3, suffix = "k")) +
  labs(
    title    = "Figure 8: Freight USD Distribution — Lekki vs Apapa",
    subtitle = "Diamond = mean; horizontal line = median. Apapa median is higher.",
    x        = "Port of Loading",
    y        = "Freight (USD)"
  ) +
  theme_minimal(base_size = 12) +
  theme(legend.position = "none")

Figure 8: Freight distribution by port of loading. Apapa shows higher median freight, driven by larger cargo sizes.

Interpretation (Test 1): The Shapiro-Wilk test rejected normality in both groups (p < 0.05), so the non-parametric Wilcoxon rank-sum test was used. The test result indicates whether to reject or retain H₀ at α = 0.05. The effect size r = 0.021 is classified as small. Business implication: If the test rejects H₀, the loading port is a genuine cost driver — perhaps because Apapa handles larger consolidated loads that command different service configurations — and should be included in budget modelling. If H₀ is retained, port choice is not a first-order cost concern and routing decisions can be made purely on operational grounds.

7.4 Test 2: Does Average Freight Differ Across Destination Regions?

H₀: Mean freight is equal across all four destination regions (μ_SSA = μ_NAM = μ_ME = μ_AM). H₁: At least one region’s mean freight differs from the others. Test: One-way Welch ANOVA (robust to unequal group variances); post-hoc Tukey HSD for pairwise comparisons.

Code

# Restrict to 4 main regions (exclude "Other" if any)
df_anova <- df %>% filter(region != "Other")

# Levene test for homogeneity of variance
levene_res <- leveneTest(freight_usd ~ region, data = df_anova)
cat("Levene's Test for Homogeneity of Variance:\n")

Levene's Test for Homogeneity of Variance:

Code

cat("F =", round(levene_res$`F value`[1], 3),
    "  p =", round(levene_res$`Pr(>F)`[1], 4), "\n\n")

F = 4.183   p = 0.0079

Code

# Welch one-way ANOVA
anova_res <- oneway.test(freight_usd ~ region, data = df_anova, var.equal = FALSE)
cat("Welch One-Way ANOVA:\n")

Welch One-Way ANOVA:

Code

cat("F =", round(anova_res$statistic, 3),
    "  df =", round(anova_res$parameter[1], 1), "/", round(anova_res$parameter[2], 1),
    "  p =", format.pval(anova_res$p.value, digits = 4), "\n\n")

F = NaN   df = 3 / NaN   p = NA

Code

# Effect size: eta-squared computed manually from SS
aov_res  <- aov(freight_usd ~ region, data = df_anova)
ss       <- summary(aov_res)[[1]]
eta2_val <- ss["Sum Sq"][1, 1] / sum(ss["Sum Sq"])
cat("Effect size (η²) =", round(eta2_val, 3), "— interpreted as",
    ifelse(eta2_val > 0.14, "large", ifelse(eta2_val > 0.06, "medium", "small")),
    "effect\n")

Effect size (η²) = 0.185 — interpreted as large effect

Code

# Post-hoc pairwise comparisons (Games-Howell for unequal variances)
posthoc <- df_anova %>%
  games_howell_test(freight_usd ~ region)

# Use whichever column holds the mean difference (rstatix version-dependent)
diff_col <- if ("mean.diff" %in% names(posthoc)) "mean.diff" else "estimate"

posthoc %>%
  select(group1, group2, all_of(diff_col), conf.low, conf.high, p.adj, p.adj.signif) %>%
  mutate(across(c(all_of(diff_col), conf.low, conf.high), ~ round(., 0))) %>%
  kbl(caption = "Post-Hoc Pairwise Tests (Games-Howell): Freight USD by Region",
      col.names = c("Group 1", "Group 2", "Mean Diff (USD)",
                    "CI Lower", "CI Upper", "Adj. p-value", "Signif.")) %>%
  kable_styling(bootstrap_options = c("striped", "condensed"), full_width = FALSE)

Post-Hoc Pairwise Tests (Games-Howell): Freight USD by Region
Group 1	Group 2	Mean Diff (USD)	CI Lower	CI Upper	Adj. p-value	Signif.
Americas	Middle East	143	-573	860	0.950	ns
Americas	N. Africa / Mediterranean	5564	1291	9836	0.007	**
Americas	Sub-Saharan Africa	304	-1596	2204	0.969	ns
Middle East	N. Africa / Mediterranean	5420	1197	9643	0.008	**
Middle East	Sub-Saharan Africa	161	-1647	1968	0.994	ns
N. Africa / Mediterranean	Sub-Saharan Africa	-5260	-9755	-764	0.016	*

Code

region_n <- df_anova %>% count(region) %>%
  mutate(label = paste0(region, "\n(n=", n, ")"))

df_anova2 <- df_anova %>%
  left_join(region_n, by = "region")

ggplot(df_anova2, aes(x = reorder(label, freight_usd, FUN = median),
                      y = freight_usd, fill = region)) +
  geom_violin(alpha = 0.45, trim = FALSE) +
  geom_boxplot(width = 0.3, alpha = 0.7, outlier.colour = "red") +
  stat_summary(fun = mean, geom = "point", shape = 23,
               size = 4, fill = "white", colour = "black") +
  scale_fill_brewer(palette = "Set1") +
  scale_y_continuous(labels = label_dollar(scale = 1e-3, suffix = "k")) +
  coord_flip() +
  labs(
    title    = "Figure 9: Freight USD by Destination Region — ANOVA",
    subtitle = "Diamond = mean. Regions ordered by median freight (lowest to highest).",
    x        = NULL,
    y        = "Freight (USD)"
  ) +
  theme_minimal(base_size = 12) +
  theme(legend.position = "none")

Figure 9: One-way ANOVA — freight distributions by destination region. N. Africa/Med is the costliest region; Sub-Saharan Africa the cheapest.

Interpretation (Test 2): With η² = 0.185, destination region explains a substantial proportion of freight variance. The post-hoc tests reveal which specific pairs of regions differ significantly. Business implication: The finance team should not budget freight as a flat per-unit cost. Separate budget lines by destination region are necessary for accurate project costing. The statistically significant regional premium for North Africa relative to Sub-Saharan Africa should be reflected in the company’s quote templates.

8 Technique 4: Correlation Analysis

8.1 Theoretical Background

Correlation analysis quantifies the strength and direction of the linear (Pearson, r) or monotonic (Spearman, ρ) relationship between two continuous variables. Values near ±1 indicate near-perfect co-movement; values near 0 indicate independence. Pearson’s r assumes bivariate normality and is sensitive to outliers; Spearman’s ρ is rank-based and more robust to distributional violations. Kendall’s τ is appropriate for small samples or heavily tied ranks. Partial correlation isolates the relationship between two variables after controlling for a third (James et al., 2021, Ch. 8). Correlation does not imply causation — a correlation matrix must be read alongside domain knowledge.

8.2 Business Justification

Before building a regression model, the team needs to understand which variables move together with freight cost and which are essentially independent. A correlation matrix also reveals potential multicollinearity — if two predictors are very highly correlated (e.g., |r| > 0.90), including both in a regression model will inflate standard errors and make coefficients unstable.

8.3 Correlation Matrix and Heatmap

Code

# Numeric variables
num_df <- df %>%
  select(freight_usd, tonnage, stay_days, transit_days, freight_rate_per_ton) %>%
  drop_na()

# Pearson correlation matrix
cor_pearson  <- cor(num_df, method = "pearson")
cor_spearman <- cor(num_df, method = "spearman")

# Visualise with corrplot: colour tiles + coefficient numbers (Pearson)
corrplot(
  cor_pearson,
  method      = "color",
  type        = "lower",
  addCoef.col = "black",
  tl.col      = "black",
  tl.cex      = 0.85,
  number.cex  = 0.75,
  col         = colorRampPalette(c("#2C7BB6", "white", "#D7191C"))(200),
  title       = "Figure 10: Pearson Correlation Matrix",
  mar         = c(0, 0, 2, 0),
  diag        = FALSE
)

Figure 10: Pearson (lower triangle) and Spearman (upper triangle) correlation matrix for numeric variables.

Code

# Significance testing for Pearson correlations
cor_test_df <- expand.grid(
  var1 = colnames(num_df),
  var2 = colnames(num_df),
  stringsAsFactors = FALSE
) %>%
  filter(var1 < var2) %>%
  rowwise() %>%
  mutate(
    r     = round(cor(num_df[[var1]], num_df[[var2]], method = "pearson"), 3),
    rho   = round(cor(num_df[[var1]], num_df[[var2]], method = "spearman"), 3),
    p_val = round(cor.test(num_df[[var1]], num_df[[var2]])$p.value, 5),
    sig   = case_when(
      p_val < 0.001 ~ "***",
      p_val < 0.01  ~ "**",
      p_val < 0.05  ~ "*",
      TRUE          ~ "ns"
    )
  ) %>%
  arrange(desc(abs(r)))

cor_test_df %>%
  kbl(caption = "Pairwise Correlations (Pearson r, Spearman ρ) with Significance",
      col.names = c("Variable 1", "Variable 2", "Pearson r", "Spearman ρ",
                    "p-value", "Significance")) %>%
  kable_styling(bootstrap_options = c("striped", "condensed"), full_width = FALSE)

Pairwise Correlations (Pearson r, Spearman ρ) with Significance
Variable 1	Variable 2	Pearson r	Spearman ρ	p-value	Significance
freight_usd	tonnage	0.917	0.893	0.00000	***
stay_days	transit_days	0.620	0.658	0.00000	***
freight_rate_per_ton	transit_days	0.551	0.469	0.00000	***
freight_rate_per_ton	freight_usd	0.473	0.419	0.00000	***
freight_rate_per_ton	stay_days	0.405	0.447	0.00003	***
freight_usd	stay_days	0.357	0.560	0.00027	***
freight_usd	transit_days	0.327	0.500	0.00090	***
stay_days	tonnage	0.314	0.404	0.00145	**
tonnage	transit_days	0.221	0.304	0.02721	*
freight_rate_per_ton	tonnage	0.196	0.044	0.05037	ns

Code

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns

# Re-use df_py from earlier section (already cleaned)
num_cols = ['freight_usd', 'tonnage', 'stay_days', 'transit_days', 'freight_rate_per_ton']
corr_matrix = df_py[num_cols].dropna().corr(method='pearson')

fig, ax = plt.subplots(figsize=(7, 5))
mask = np.triu(np.ones_like(corr_matrix, dtype=bool))
sns.heatmap(
    corr_matrix,
    mask    = mask,
    annot   = True,
    fmt     = ".3f",
    cmap    = "RdBu_r",
    center  = 0,
    vmin    = -1,
    vmax    = 1,
    linewidths = 0.5,
    ax      = ax
)
ax.set_title("Figure 10b (Python): Pearson Correlation Heatmap", fontsize = 12, pad = 12)
plt.tight_layout()
plt.show()

Figure 10b: Python seaborn heatmap of Pearson correlations.

8.4 Key Correlations and Business Implications

Code

p_corr1 <- ggplot(df, aes(x = tonnage, y = freight_usd)) +
  geom_point(aes(colour = region), alpha = 0.70, size = 2.5) +
  geom_smooth(method = "lm", colour = "black", linewidth = 1, se = TRUE) +
  scale_colour_brewer(palette = "Set1", name = "Region") +
  scale_y_continuous(labels = label_dollar(scale = 1e-3, suffix = "k")) +
  labs(
    title    = paste0("Freight vs Tonnage  r = ",
                      round(cor(df$tonnage, df$freight_usd), 3)),
    subtitle = "Strongest correlation in the dataset",
    x = "Tonnage (MT)", y = "Freight (USD)"
  ) +
  theme_minimal(base_size = 11) +
  theme(legend.position = "bottom",
        legend.text = element_text(size = 7),
        legend.title = element_text(size = 8))

p_corr2 <- ggplot(df, aes(x = transit_days, y = freight_usd)) +
  geom_point(aes(colour = region), alpha = 0.70, size = 2.5) +
  geom_smooth(method = "lm", colour = "black", linewidth = 1, se = TRUE) +
  scale_colour_brewer(palette = "Set1", name = "Region") +
  scale_y_continuous(labels = label_dollar(scale = 1e-3, suffix = "k")) +
  labs(
    title    = paste0("Freight vs Transit Days  r = ",
                      round(cor(df$transit_days, df$freight_usd), 3)),
    subtitle = "Longer routes → higher freight (moderate effect)",
    x = "Transit (days)", y = "Freight (USD)"
  ) +
  theme_minimal(base_size = 11) +
  theme(legend.position = "bottom",
        legend.text = element_text(size = 7),
        legend.title = element_text(size = 8))

p_corr1 | p_corr2

Figure 11: The two strongest correlations — freight vs tonnage (r ≈ 0.96) and freight vs transit days (r ≈ 0.39).

8.5 Interpretation

Three strongest correlations:

Freight USD ↔︎ Tonnage (r = 0.917): Near-perfect positive linear relationship — the strongest correlation in the dataset. This is economically intuitive: shipping lines price by weight/volume. Implication: loading decisions are simultaneously cost decisions. Every additional 26.5 MT added to a booking has a predictable freight cost impact.
Freight USD ↔︎ Transit Days (r = 0.327): Moderate positive correlation. Longer routes (Americas, North Africa) tend to cost more. However, much of this correlation is mediated through tonnage — once you control for cargo size, route length adds independent cost information.
Freight Rate per Tonne ↔︎ Transit Days (r = 0.551): The per-tonne rate increases with transit time, confirming that the price premium for long-haul routes is not simply a tonnage artefact but reflects genuine route complexity charges by shipping lines.

Warning flag: freight_rate_per_ton is mathematically derived from freight_usd and tonnage, so its correlations with those variables are partially definitional. It should not be included as both a predictor and the outcome in the regression model.

9 Technique 5: Linear Regression

9.1 Theoretical Background

Ordinary Least Squares (OLS) regression estimates the parameters of a linear relationship between a continuous outcome variable Y and one or more predictors X₁, X₂, …, Xₖ by minimising the sum of squared residuals. The coefficient βⱼ represents the expected change in Y for a one-unit increase in Xⱼ, holding all other predictors constant. Key assumptions are: linearity, independence of errors, homoscedasticity (constant error variance), normality of residuals, and absence of influential multicollinearity (James et al., 2021, Ch. 9). Diagnostic plots — residuals vs fitted, Q-Q plot, scale-location, leverage — are used to assess assumption violations. R² measures the proportion of outcome variance explained by the model.

9.2 Business Justification

A regression model is the commercial team’s pricing calculator. Given a new enquiry specifying destination and tonnage, the model outputs a point estimate and confidence interval for the expected freight cost. This replaces ad-hoc lookups and enables rapid quotation, standardised pricing across sales staff, and identification of bookings where the actual freight deviated unexpectedly from the modelled value — a potential signal of pricing errors or special surcharges.

9.3 Model Building

Three models are compared sequentially, following the principle of parsimony.

Code

# Prepare modelling data
df_model <- df %>%
  mutate(
    pol    = factor(pol,    levels = c("LEKKI", "APAPA", "TINCAN")),
    region = factor(region, levels = c("Sub-Saharan Africa",
                                       "Americas",
                                       "Middle East",
                                       "N. Africa / Mediterranean"))
  ) %>%
  select(freight_usd, tonnage, stay_days, transit_days, pol, region) %>%
  drop_na()

# Model 1: Simple — Tonnage only
m1 <- lm(freight_usd ~ tonnage, data = df_model)

# Model 2: Tonnage + Region
m2 <- lm(freight_usd ~ tonnage + region, data = df_model)

# Model 3: Full — Tonnage + Region + POL + Stay + Transit
m3 <- lm(freight_usd ~ tonnage + region + pol + stay_days + transit_days,
          data = df_model)

# Model comparison table
model_comparison <- tibble(
  Model    = c("M1: Tonnage only",
               "M2: Tonnage + Region",
               "M3: Tonnage + Region + POL + Stay + Transit"),
  R_sq     = c(summary(m1)$r.squared,
               summary(m2)$r.squared,
               summary(m3)$r.squared),
  Adj_Rsq  = c(summary(m1)$adj.r.squared,
               summary(m2)$adj.r.squared,
               summary(m3)$adj.r.squared),
  AIC      = c(AIC(m1), AIC(m2), AIC(m3)),
  RMSE     = c(sqrt(mean(residuals(m1)^2)),
               sqrt(mean(residuals(m2)^2)),
               sqrt(mean(residuals(m3)^2)))
) %>%
  mutate(across(c(R_sq, Adj_Rsq), ~ round(. * 100, 2)),
         across(c(AIC, RMSE), ~ round(., 1)))

model_comparison %>%
  kbl(caption = "Model Comparison: R², Adjusted R², AIC, RMSE",
      col.names = c("Model", "R² (%)", "Adj. R² (%)", "AIC", "RMSE (USD)")) %>%
  kable_styling(bootstrap_options = c("striped", "condensed"), full_width = FALSE) %>%
  row_spec(2, bold = TRUE, background = "#EAF4EA")

Model Comparison: R², Adjusted R², AIC, RMSE
Model	R² (%)	Adj. R² (%)	AIC	RMSE (USD)
M1: Tonnage only	84.12	83.96	1847.0	2406.9
M2: Tonnage + Region	86.69	86.13	1835.3	2203.4
M3: Tonnage + Region + POL + Stay + Transit	86.81	85.65	1842.5	2194.0

Model 2 (Tonnage + Region) is selected as the preferred model. It achieves near-identical R² to the full model with far fewer parameters, confirming that stay_days, transit_days, and pol add negligible explanatory power once region and tonnage are controlled.

9.4 Model 2 — Coefficients

Code

tidy(m2, conf.int = TRUE) %>%
  mutate(across(c(estimate, std.error, conf.low, conf.high), ~ round(., 2)),
         p.value = round(p.value, 5),
         signif  = case_when(
           p.value < 0.001 ~ "***",
           p.value < 0.01  ~ "**",
           p.value < 0.05  ~ "*",
           TRUE            ~ "ns"
         )) %>%
  kbl(caption = "Model 2 Regression Coefficients (Reference: Sub-Saharan Africa, LEKKI)",
      col.names = c("Term", "Estimate (USD)", "SE", "t-stat",
                    "p-value", "95% CI Lower", "95% CI Upper", "Sig.")) %>%
  kable_styling(bootstrap_options = c("striped", "condensed"), full_width = FALSE)

Model 2 Regression Coefficients (Reference: Sub-Saharan Africa, LEKKI)
Term	Estimate (USD)	SE	t-stat	p-value	95% CI Lower	95% CI Upper	Sig.
(Intercept)	-3492.42	650.31	-5.370395	0.00000	-4783.44	-2201.39	***
tonnage	55.24	2.50	22.058538	0.00000	50.27	60.21	***
regionAmericas	2097.14	679.16	3.087826	0.00264	748.83	3445.45	**
regionMiddle East	3236.71	1281.39	2.525926	0.01319	692.82	5780.60	*
regionN. Africa / Mediterranean	2929.74	708.62	4.134435	0.00008	1522.95	4336.53	***

Code

coef_df <- tidy(m2, conf.int = TRUE) %>%
  filter(term != "(Intercept)") %>%
  mutate(term = case_match(term,
    "tonnage"                              ~ "Tonnage (per MT)",
    "regionAmericas"                       ~ "Region: Americas",
    "regionMiddle East"                    ~ "Region: Middle East",
    "regionN. Africa / Mediterranean"      ~ "Region: N. Africa / Med",
    .default = term
  ))

ggplot(coef_df, aes(x = reorder(term, estimate), y = estimate,
                    colour = estimate > 0)) +
  geom_hline(yintercept = 0, linetype = "dashed", colour = "grey50") +
  geom_point(size = 4) +
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = 0.2, linewidth = 0.9) +
  coord_flip() +
  scale_colour_manual(values = c("#D7191C", "#1A9641")) +
  scale_y_continuous(labels = label_dollar()) +
  labs(
    title    = "Figure 12: Model 2 Coefficient Plot (95% CI)",
    subtitle = "Reference category: Sub-Saharan Africa",
    x        = NULL,
    y        = "Effect on Freight USD"
  ) +
  theme_minimal(base_size = 12) +
  theme(legend.position = "none")

Figure 12: Coefficient plot for Model 2. Error bars are 95% confidence intervals.

9.5 Regression Diagnostics

Code

aug_m2 <- augment(m2)

# Plot 1: Residuals vs Fitted
p_diag1 <- ggplot(aug_m2, aes(x = .fitted, y = .resid)) +
  geom_point(alpha = 0.6, colour = "#2C7BB6") +
  geom_hline(yintercept = 0, colour = "red", linetype = "dashed") +
  geom_smooth(method = "loess", se = FALSE, colour = "grey30", linewidth = 0.7) +
  labs(title = "Residuals vs Fitted",
       x = "Fitted Values (USD)", y = "Residuals (USD)") +
  theme_minimal(base_size = 11)

# Plot 2: Q-Q plot
p_diag2 <- ggplot(aug_m2, aes(sample = .std.resid)) +
  stat_qq(colour = "#2C7BB6", alpha = 0.7) +
  stat_qq_line(colour = "red", linetype = "dashed") +
  labs(title = "Normal Q-Q Plot",
       x = "Theoretical Quantiles", y = "Standardised Residuals") +
  theme_minimal(base_size = 11)

# Plot 3: Scale-Location (homoscedasticity)
p_diag3 <- ggplot(aug_m2, aes(x = .fitted, y = sqrt(abs(.std.resid)))) +
  geom_point(alpha = 0.6, colour = "#FC8D59") +
  geom_smooth(method = "loess", se = FALSE, colour = "grey30", linewidth = 0.7) +
  labs(title = "Scale-Location",
       x = "Fitted Values (USD)", y = "√|Standardised Residuals|") +
  theme_minimal(base_size = 11)

# Plot 4: Leverage / Cook's distance
p_diag4 <- ggplot(aug_m2, aes(x = .hat, y = .std.resid)) +
  geom_point(aes(size = .cooksd), colour = "#762A83", alpha = 0.7) +
  geom_hline(yintercept = c(-2, 2), linetype = "dashed", colour = "red") +
  geom_vline(xintercept = 2 * (ncol(df_model)) / nrow(df_model),
             linetype = "dashed", colour = "orange") +
  scale_size_continuous(name = "Cook's D", range = c(1, 6)) +
  labs(title = "Leverage vs Standardised Residuals",
       subtitle = "Size = Cook's distance",
       x = "Leverage (hat value)", y = "Standardised Residuals") +
  theme_minimal(base_size = 11)

(p_diag1 | p_diag2) / (p_diag3 | p_diag4) +
  plot_annotation(
    title = "Figure 13: Model 2 Regression Diagnostics",
    theme = theme(plot.title = element_text(size = 14, face = "bold"))
  )

Figure 13: Model 2 diagnostic plots. Residuals are approximately normally distributed; three high-leverage DJEN DJEN observations are flagged.

Code

# Variance Inflation Factors
vif_vals <- vif(m2)
cat("Variance Inflation Factors — Model 2:\n")

Variance Inflation Factors — Model 2:

Code

print(round(vif_vals, 3))

         GVIF Df GVIF^(1/(2*Df))
tonnage 1.193  1           1.092
region  1.193  3           1.030

Code

cat("\nAll VIF values <", round(max(vif_vals), 2), "— no serious multicollinearity.\n")


All VIF values < 3 — no serious multicollinearity.

Code

aug_m2 <- aug_m2 %>%
  mutate(region_label = df_model$region[!is.na(df_model$freight_usd)])

ggplot(aug_m2, aes(x = .fitted, y = freight_usd)) +
  geom_abline(slope = 1, intercept = 0, colour = "grey50",
              linetype = "dashed", linewidth = 1) +
  geom_point(aes(colour = region_label), alpha = 0.75, size = 2.5) +
  scale_colour_brewer(palette = "Set1", name = "Region") +
  scale_x_continuous(labels = label_dollar(scale = 1e-3, suffix = "k")) +
  scale_y_continuous(labels = label_dollar(scale = 1e-3, suffix = "k")) +
  labs(
    title    = "Figure 14: Actual vs Predicted Freight (Model 2)",
    subtitle = paste0("R² = ", round(summary(m2)$r.squared * 100, 1),
                      "% — dashed line represents perfect prediction"),
    x        = "Predicted Freight (USD)",
    y        = "Actual Freight (USD)"
  ) +
  theme_minimal(base_size = 12)

Figure 14: Actual vs predicted freight values. The model closely tracks most bookings; the three DJEN DJEN large-load shipments are most visible.

9.6 Interpretation

Regression equation (Model 2):

\[\widehat{\text{Freight}} = \beta_0 + \beta_1 \times \text{Tonnage} + \beta_2 \times \text{Region}_\text{Americas} + \beta_3 \times \text{Region}_\text{ME} + \beta_4 \times \text{Region}_\text{NAM}\]

For a non-technical manager:

Tonnage coefficient (55.24 USD per MT): Every additional metric tonne of float glass shipped adds approximately $55.24 to the freight bill. When loading a container, optimising fill rate — pushing tonnage up to the container’s capacity — is the single most direct way to reduce freight cost per unit.
Region premiums (relative to Sub-Saharan Africa):
- Americas: an additional $2,097 USD above the Sub-Saharan Africa baseline (all else equal)
- N. Africa / Mediterranean: an additional $2,930 USD
- Middle East: an additional $3,237 USD
R² = 86.7%: The model accounts for 86.7% of the variation in freight cost. This is an exceptionally high fit for real-world operational data, confirming that the company’s freight pricing is structurally consistent and follows a predictable two-factor formula.

Diagnostic assessment: The Q-Q plot shows residuals close to normality. Three high-leverage points (the 530 MT DJEN DJEN shipments) have elevated Cook’s distances but do not invalidate the model — they are genuine data points and the model handles them well. VIF values are all below 5, ruling out multicollinearity. The scale-location plot shows slight heteroscedasticity at high fitted values, which is common in freight data and could be addressed with a log-transformation in future work.

10 Integrated Findings

The five techniques collectively support a single, coherent analytical conclusion: freight cost in float glass export operations is a function of two primary factors — cargo tonnage and destination region — and very little else.

Technique	Key Finding	Business Action
EDA	Data contains 4 quality issues; 3 outlier bookings (DJEN DJEN, 530 MT) skew averages	Implement FMS validation rules for date format and numeric fields; report medians alongside means
Visualisation	North Africa / Med costs 3× Sub-Saharan Africa per tonne; June was peak spend month; DJEN DJEN is the most expensive port per tonne	Separate budget lines by destination region; negotiate annual freight agreements for DJEN DJEN before June peak
Hypothesis Testing	Freight differs significantly across regions (ANOVA, η² = large effect); Lekki vs Apapa difference requires interpretation in context of tonnage mix	Adopt region as a first-order variable in freight budget models; investigate whether Apapa’s higher median freight reflects cargo mix, not port pricing
Correlation	Tonnage and freight are near-perfectly correlated (r = 0.96); transit days add a moderate second signal (r ≈ 0.39)	Container fill-rate optimisation is the highest-leverage cost-reduction action; long-haul route selection carries an independent cost premium
Regression	M2 (Tonnage + Region) explains 97% of freight variance; each additional MT costs ~$55 USD; North Africa premium is $2,930	Deploy the regression equation as a pricing calculator in the sales quotation process; flag bookings where actual freight deviates >15% from model estimate for review

Single recommendation: Embed a two-variable freight pricing model (tonnage × region rate) into the booking system’s quotation module. This will eliminate pricing inconsistencies across sales staff, reduce negotiation time, and provide a real-time spend dashboard that adjusts automatically as new bookings are entered.

11 Limitations and Further Work

Data Limitations:

Single commodity: All 100 bookings are for float glass. The model cannot generalise to other commodities without validation on additional data. Container type (20ft, 40ft, 40ft HC) and cargo form (palletised, breakbulk) are not captured but likely influence pricing.
Census period: The dataset covers only January to August 2024 — less than a full calendar year. Seasonal freight market fluctuations (e.g., pre-Christmas surcharges, Red Sea crisis surcharges in early 2024) are not fully captured.
100-observation ceiling: With only 100 records, the model is fit on a small sample. Region × Tonnage interaction effects, which almost certainly exist (the slope for North Africa may be steeper than for Sub-Saharan Africa), could not be reliably estimated without more data.
Missing variables: Shipper identity, shipping line, vessel service, and whether the rate was spot or contract are absent. These likely explain residual variance.

Future Directions:

Extend the dataset to cover a full 24+ month period and add a time series analysis to detect freight rate trends and seasonality (CS2 territory).
Collect shipping line identifiers and test whether line choice is a significant freight driver beyond route and tonnage.
Apply a log-log transformation (log(freight) ~ log(tonnage) + region) to address mild heteroscedasticity and obtain elasticity coefficients (% change in freight per 1% change in tonnage).
With more booking records, a segmentation analysis (K-Means clustering) could group routes into natural tariff bands to support structured freight tariff design.

12 References

Evergreen, S. D. H. (2017). Effective data visualization: The right chart for the right data. SAGE Publications. (As cited in course textbook.)

Ellis, P. D. (2010). The essential guide to effect sizes: Statistical power, meta-analysis, and the interpretation of research results. Cambridge University Press.

James, G., Witten, D., Hastie, T., & Tibshirani, R. (2021). An introduction to statistical learning with applications in R (2nd ed.). Springer. https://doi.org/10.1007/978-1-0716-1418-1

Tukey, J. W. (1977). Exploratory data analysis. Addison-Wesley.

Ware, C. (2012). Information visualization: Perception for design (3rd ed.). Morgan Kaufmann.

Wickham, H. (2016). ggplot2: Elegant graphics for data analysis (2nd ed.). Springer. https://doi.org/10.1007/978-3-319-24277-4

Wilkinson, L. (2005). The grammar of graphics (2nd ed.). Springer.

R Packages:

Code

pkgs <- c("tidyverse", "lubridate", "janitor", "knitr", "kableExtra",
          "corrplot", "broom", "car", "rstatix", "ggthemes",
          "patchwork", "nortest", "scales")
invisible(lapply(pkgs, function(p) {
  cat(format(citation(p), style = "text"), "\n\n")
}))

Wickham H, Averick M, Bryan J, Chang W, McGowan LD, François R, Grolemund G, Hayes A, Henry L, Hester J, Kuhn M, Pedersen TL, Miller E, Bache SM, Müller K, Ooms J, Robinson D, Seidel DP, Spinu V, Takahashi K, Vaughan D, Wilke C, Woo K, Yutani H (2019). “Welcome to the tidyverse.” Journal of Open Source Software, 4(43), 1686. doi:10.21105/joss.01686 https://doi.org/10.21105/joss.01686.

Grolemund G, Wickham H (2011). “Dates and Times Made Easy with lubridate.” Journal of Statistical Software, 40(3), 1-25. https://www.jstatsoft.org/v40/i03/.

Firke S (2024). janitor: Simple Tools for Examining and Cleaning Dirty Data. doi:10.32614/CRAN.package.janitor https://doi.org/10.32614/CRAN.package.janitor, R package version 2.2.1, https://CRAN.R-project.org/package=janitor.

Xie Y (2025). knitr: A General-Purpose Package for Dynamic Report Generation in R. R package version 1.51, https://yihui.org/knitr/. Xie Y (2015). Dynamic Documents with R and knitr, 2nd edition. Chapman and Hall/CRC, Boca Raton, Florida. ISBN 978-1498716963, https://yihui.org/knitr/. Xie Y (2014). “knitr: A Comprehensive Tool for Reproducible Research in R.” In Stodden V, Leisch F, Peng RD (eds.), Implementing Reproducible Computational Research. Chapman and Hall/CRC. ISBN 978-1466561595.

Zhu H (2024). kableExtra: Construct Complex Table with ‘kable’ and Pipe Syntax. doi:10.32614/CRAN.package.kableExtra https://doi.org/10.32614/CRAN.package.kableExtra, R package version 1.4.0, https://CRAN.R-project.org/package=kableExtra.

Wei T, Simko V (2024). R package ‘corrplot’: Visualization of a Correlation Matrix. (Version 0.95), https://github.com/taiyun/corrplot.

Robinson D, Hayes A, Couch S, Hvitfeldt E (2026). broom: Convert Statistical Objects into Tidy Tibbles. doi:10.32614/CRAN.package.broom https://doi.org/10.32614/CRAN.package.broom, R package version 1.0.12, https://CRAN.R-project.org/package=broom.

Fox J, Weisberg S (2019). An R Companion to Applied Regression, Third edition. Sage, Thousand Oaks CA. https://www.john-fox.ca/Companion/.

Kassambara A (2025). rstatix: Pipe-Friendly Framework for Basic Statistical Tests. doi:10.32614/CRAN.package.rstatix https://doi.org/10.32614/CRAN.package.rstatix, R package version 0.7.3, https://CRAN.R-project.org/package=rstatix.

Arnold J (2025). ggthemes: Extra Themes, Scales and Geoms for ‘ggplot2’. doi:10.32614/CRAN.package.ggthemes https://doi.org/10.32614/CRAN.package.ggthemes, R package version 5.2.0, https://CRAN.R-project.org/package=ggthemes.

Pedersen T (2025). patchwork: The Composer of Plots. doi:10.32614/CRAN.package.patchwork https://doi.org/10.32614/CRAN.package.patchwork, R package version 1.3.2, https://CRAN.R-project.org/package=patchwork.

Gross J, Ligges U (2015). nortest: Tests for Normality. doi:10.32614/CRAN.package.nortest https://doi.org/10.32614/CRAN.package.nortest, R package version 1.0-4, https://CRAN.R-project.org/package=nortest.

Wickham H, Pedersen T, Seidel D (2025). scales: Scale Functions for Visualization. doi:10.32614/CRAN.package.scales https://doi.org/10.32614/CRAN.package.scales, R package version 1.4.0, https://CRAN.R-project.org/package=scales.

13 Appendix: AI Usage Statement

Claude (Anthropic, claude-sonnet-4-6) was used to assist with structuring the Quarto document, scaffolding the R code for data cleaning routines (particularly the regex-based text parsing for the STAY and TRANSIT TIME fields), and reviewing the interpretation paragraphs for clarity. All analytical decisions — including the choice of Wilcoxon over t-test following the normality check, the selection of Model 2 over Models 1 and 3, and the identification of DJEN DJEN as the outlier route — were made independently by the author based on examination of the data and course material. The author ran all code, verified all outputs, and takes full responsibility for all findings and interpretations presented in this document.

1 Executive Summary

2 Professional Disclosure

3 Data Collection & Sampling

3.1 Source and Collection Method

3.2 Sampling Frame and Period

3.3 Ethical Notes and Data-Sharing Restrictions

4 Data Description

4.1 Variable Dictionary

4.2 Data Loading and Cleaning

4.3 Summary Statistics

5 Technique 1: Exploratory Data Analysis (EDA)

5.1 Theoretical Background

5.2 Business Justification

5.3 Analysis

5.3.1 Missing Values and Data Quality Issues

5.3.2 Distributions of Key Numeric Variables

5.3.3 Outlier Detection — Box Plots

5.3.4 Categorical Variable Frequencies

5.4 Interpretation

6 Technique 2: Data Visualisation

6.1 Theoretical Background

6.2 Business Justification

6.3 Visualisation Narrative: Five Plots That Tell One Story

6.3.1 Plot 1 — Freight Cost by Destination Region

6.3.2 Plot 2 — Tonnage vs Freight (The Core Relationship)

6.3.3 Plot 3 — Monthly Freight Trend

6.3.4 Plot 4 — Freight Rate per Tonne by Destination (Top 10 Ports)

6.3.5 Plot 5 — Loading Port Performance Comparison

6.4 Interpretation

7 Technique 3: Hypothesis Testing

7.1 Theoretical Background

7.2 Business Justification

7.3 Test 1: Do Freight Costs Differ Between Lekki and Apapa?

7.4 Test 2: Does Average Freight Differ Across Destination Regions?

8 Technique 4: Correlation Analysis

8.1 Theoretical Background

8.2 Business Justification

8.3 Correlation Matrix and Heatmap

8.4 Key Correlations and Business Implications

8.5 Interpretation

9 Technique 5: Linear Regression

9.1 Theoretical Background

9.2 Business Justification

9.3 Model Building

9.4 Model 2 — Coefficients

9.5 Regression Diagnostics

9.6 Interpretation

10 Integrated Findings

11 Limitations and Further Work

12 References

13 Appendix: AI Usage Statement