Statistics Canada's Canadian Freight Analysis Framework records every commodity flow between provinces — 58,721 rows of shipments, tonnage, revenue and mode from 2011–2017. I engineered 13 corridor-level features, compressed them with PCA, and segmented 194 origin→destination corridors into five market archetypes with K-Means — then stress-tested the segmentation with bootstrap resampling, hierarchical and Gaussian-mixture cross-checks, permutation feature importance and a Mahalanobis anomaly audit. The result is a validated market map that feeds T-Baar's demand-model calibration, plus a live classifier you can steer yourself.
Each card is a cluster from the final k=5 model. Sizes, revenue shares and silhouettes are computed from the re-run pipeline; the top pills are each cluster's largest corridors by total revenue. Click a card to isolate that archetype in the explorer below.
Each dot is one of the 194 corridors, placed by t-SNE over the 5 retained principal components (or raw PC1×PC2 — toggle below). Size ∝ √revenue. Color by archetype, by anomaly score (Mahalanobis distance to own-cluster centroid), or by the corridor's 2011→2017 revenue-trend direction. Hover for details, click to pin a full profile.
Choosing k is where most clustering projects hand-wave. Here five independent criteria are computed for k = 2…9 on the PCA-reduced matrix: inertia (elbow), mean silhouette, Calinski–Harabasz, a spherical-Gaussian BIC for K-Means, an EM-fitted diagonal GMM BIC, and the gap statistic against 40 uniform reference sets. Silhouette and CH plateau at k≈4–6 while the elbow flattens past 5 — but the decisive evidence is stability: reclustering 120 bootstrap resamples and comparing each to the full-sample partition (Adjusted Rand Index) shows k=5 is the most reproducible choice, ARI 0.90.
A segmentation you can only get from one algorithm is fragile. Cutting a from-scratch Ward hierarchy at 5 clusters recovers the K-Means partition at ARI 0.75 — strong agreement given Ward's very different objective. A diagonal-covariance Gaussian mixture fitted by EM prefers similar mid-range k by BIC. And the re-run reproduces the originally archived labels at ARI 1.00 — the pipeline is fully deterministic and reproducible.
| Check | Method | Result | Reading |
|---|---|---|---|
| Reproducibility | re-run, seed 42, n_init 100 | ARI 1.00 vs archived | deterministic pipeline |
| Algorithm agreement | Ward linkage (Lance–Williams), cut k=5 | ARI 0.75 | same macro-structure |
| Model-based check | diagonal GMM via EM, BIC over k=2–8 | favors k≈5–7 | no strong over-split signal |
| Sampling robustness | 120 bootstrap resamples | ARI 0.90 [0.74, 1.00] | k=5 = stability peak |
| Separability | nearest-centroid classifier, 13-d scaled space | 100% training assignment recovery | clusters are convex & separable |
The heatmap shows each cluster's mean z-score on all 13 features — hover any cell for the
cluster's raw median. Two importance measures agree on the drivers: one-way ANOVA F (between/within variance)
and permutation importance (accuracy drop of the nearest-centroid classifier when a feature is shuffled, 30 reps).
Mode split (rl_share/tf_share) and the cross_border flag dominate — geography and
modal economics, not just size, define Canada's freight markets.
Segmentations get used blindly, so I audit the misfits. For each corridor I compute the
Mahalanobis distance to its own cluster's centroid (diagonal covariance, PC space) — corridors far from every
archetype are flagged rather than silently averaged in. The flags are sensible: territorial self-loops
(YT→YT, NT→NT) with unusual rate structures, SK→UM waste exports, and
BC→BC, which is trunk-like in volume but priced like a periphery corridor. In the T-Baar demand model
these get corridor-specific overrides instead of archetype-level parameters.
| # | Corridor | Archetype | Mahalanobis | Silhouette | Revenue | Dominant commodity | Why it's odd |
|---|
Total 2011–17 freight revenue for every origin (rows) → destination (columns) pair, log-scaled.
The bright diagonal is intra-provincial trucking; the off-diagonal hot cells — ON↔QC,
AB↔BC, and the US column — are exactly the corridors the Cross-Border Arteries and
Domestic Trunk archetypes capture. Hover any cell.
This is the fitted model running in your browser — the same scaler (log1p → z-score → ±4σ clip) and
the same 5 centroids in 13-d feature space. Load a real corridor as a starting point, then drag the sliders:
the nearest-centroid assignment, softmax confidence and distance profile update live. Try loading
ON→QC and dragging rail share up — you can watch a trunk route migrate into Rail-Heavy Bulk.
Parse comma-formatted numerics, drop incomplete rows, build the corridor key OrigProv→DestProv, and
drop "other international" pseudo-corridors (single-commodity data artifacts). 194 corridors survive.
Volume (shipments, revenue, value, tonne-km), economics (revenue per tonne-km, value/kg, revenue per shipment),
structure (median distance, shipment weight), mix (Shannon commodity entropy H = −Σ pᵢ log₂ pᵢ,
truck/rail/air shares), a cross-border flag, and a per-corridor OLS revenue-trend slope over the 7 years.
Five heavy-tailed volume features get log1p; all features are standardized, then clipped at 4σ so a
single mega-corridor (ON→ON) can't dominate distance computations without discarding it.
Five components retain 90.0% of variance. PC1 loads on scale (volume/revenue), PC2 on modal split, PC3 on cross-border reach — the loading heatmap in §04 shows the full rotation.
Elbow, silhouette, Calinski–Harabasz, K-Means BIC, GMM BIC and the gap statistic — then k=5 confirmed as the bootstrap-stability peak (ARI 0.90 over 120 resamples).
Ward-linkage cross-check (ARI 0.75), permutation feature importance, per-corridor silhouettes, and a Mahalanobis anomaly audit that flags 8 corridors for model-level overrides.
for k in range(2, 10):
km = KMeans(n_clusters=k, n_init=30, random_state=42).fit(X_pca)
sil = silhouette_score(X_pca, km.labels_) # peak/plateau at k≈5
ari = np.mean([adjusted_rand(km.labels_[idx], # bootstrap stability
KMeans(k).fit_predict(X_pca[idx]))
for idx in bootstrap_indices(B=120)])