RFM with K=4: Four Clusters, But Only Two Actually Say Something

Chapter 4 trained the models and compared metrics. Now I open up the actual black box: SHAP for the delay and review scenarios, and the final decision on how many customer groups exist in the RFM data. Full executed notebook, public on Colab. This notebook also starts from scratch (reloads the 9 CSVs, rebuilds the master dataframe) and rebuilds Chapter 3's features:

geo_media = geolocation.groupby('geolocation_zip_code_prefix')[['geolocation_lat', 'geolocation_lng']].mean().reset_index()

def haversine(lat1, lon1, lat2, lon2):
    R = 6371
    phi1, phi2 = np.radians(lat1), np.radians(lat2)
    dphi = np.radians(lat2 - lat1)
    dlambda = np.radians(lon2 - lon1)
    a = np.sin(dphi / 2) ** 2 + np.cos(phi1) * np.cos(phi2) * np.sin(dlambda / 2) ** 2
    return 2 * R * np.arcsin(np.sqrt(a))

m = mestre.copy()
m['customer_zip_code_prefix'] = m['customer_zip_code_prefix'].astype('int64')
m['seller_zip_code_prefix'] = m['seller_zip_code_prefix'].astype('int64')
geo_media['geolocation_zip_code_prefix'] = geo_media['geolocation_zip_code_prefix'].astype('int64')
m = m.merge(geo_media.rename(columns={'geolocation_zip_code_prefix': 'customer_zip_code_prefix', 'geolocation_lat': 'cust_lat', 'geolocation_lng': 'cust_lng'}), on='customer_zip_code_prefix', how='left')
m = m.merge(geo_media.rename(columns={'geolocation_zip_code_prefix': 'seller_zip_code_prefix', 'geolocation_lat': 'sell_lat', 'geolocation_lng': 'sell_lng'}), on='seller_zip_code_prefix', how='left')
m['distance_km'] = haversine(m['cust_lat'], m['cust_lng'], m['sell_lat'], m['sell_lng'])
m['approval_hours'] = (m['order_approved_at'] - m['order_purchase_timestamp']).dt.total_seconds() / 3600
m['month'] = m['order_purchase_timestamp'].dt.month
m['atraso_dias'] = (m['order_delivered_customer_date'] - m['order_estimated_delivery_date']).dt.days
m['delivery_days'] = (m['order_delivered_customer_date'] - m['order_purchase_timestamp']).dt.days

Same code as always: distance_km (Haversine, with Chapter 3's int64 fix), approval_hours, month, atraso_dias, and delivery_days, all already covered carefully back there.

SHAP: why the model thinks it'll be late

feature_importances_, which I already used in Chapter 3, is a global average computed during training, it measures how much each feature helped split the classes across the whole tree. SHAP (SHapley Additive exPlanations) is more rigorous: for each individual prediction, it computes how much each feature pushed that specific prediction up or down, starting from an "average prediction" baseline. The sum of those contributions explains the exact final value. Doing that math for every order in the test set and averaging the absolute value gets me a more reliable importance than Chapter 3's.

I used XGBoost here (not Random Forest), because it was the only one that actually "bets" on the positive class in Chapter 4, it makes more sense to open the black box of a model that decides something rather than one that never crosses the threshold.

delay_df = (
    m.drop_duplicates(subset='order_id')
    [['order_id', 'distance_km', 'product_weight_g', 'price', 'freight_value', 'approval_hours', 'month', 'atraso_dias']]
    .dropna()
)
delay_df['atrasado'] = (delay_df['atraso_dias'] > 0).astype(int)

features_1 = ['distance_km', 'product_weight_g', 'price', 'freight_value', 'approval_hours', 'month']
X1, y1 = delay_df[features_1], delay_df['atrasado']
X1tr, X1te, y1tr, y1te = train_test_split(X1, y1, test_size=0.2, random_state=42, stratify=y1)

modelo_delay = xgb.XGBClassifier(n_estimators=100, max_depth=6, random_state=42, eval_metric='logloss')
modelo_delay.fit(X1tr, y1tr)

explainer_delay = shap.TreeExplainer(modelo_delay)
shap_values_delay = explainer_delay.shap_values(X1te)
importancia_delay = np.abs(shap_values_delay).mean(axis=0)

I retrain XGBoost with the same hyperparameters from Chapter 4 and call shap.TreeExplainer(modelo_delay), an explainer built specifically for tree-based models that computes each feature's exact contribution without needing sampling-based approximation. explainer_delay.shap_values(X1te) returns a matrix the size of the whole test set (one row per order, one column per feature), each cell is how much that specific feature pushed that specific prediction up or down. np.abs(...).mean(axis=0) takes the absolute value (direction doesn't matter, just how much it pushed) and averages by column, that's the number that becomes the chart below.

Loading real data...

Same order as Chapter 3: month leads by a wide margin, distance_km in second place. Two different methods (tree importance and SHAP) agreeing on the same order is a good sign, the finding isn't an artifact of one specific way of measuring, it's robust to the method.

SHAP: what weighs on the score

Same logic for the XGBoost model predicting review_score:

s2 = (
    m.drop_duplicates(subset='order_id')
    [['order_id', 'atraso_dias', 'delivery_days', 'price', 'freight_value', 'payment_installments', 'review_score']]
    .dropna()
)
features_2 = ['atraso_dias', 'delivery_days', 'price', 'freight_value', 'payment_installments']
X2, y2 = s2[features_2], s2['review_score']
X2tr, X2te, y2tr, y2te = train_test_split(X2, y2, test_size=0.2, random_state=42)

modelo_review = xgb.XGBRegressor(n_estimators=100, max_depth=6, random_state=42)
modelo_review.fit(X2tr, y2tr)

explainer_review = shap.TreeExplainer(modelo_review)
shap_values_review = explainer_review.shap_values(X2te)
importancia_review = np.abs(shap_values_review).mean(axis=0)

Same recipe as the previous section (train XGBoost, TreeExplainer, shap_values, average absolute value), only what comes out the other end changes because the target here is continuous (review_score), not a class: it's a regressor, not a classifier.

Loading real data...

Same order as Chapter 3's mutual information again: delay and delivery time dominate, everything else is noise by comparison. Notice something: XGBoost showed up in both scenarios I chose to open with SHAP, but that doesn't mean it's always the better model. In the threshold sweep I already prepared for Chapter 0, at the optimal point (0.14, not the 0.5 default Chapter 4 used) Random Forest actually beats XGBoost on F1. Hold onto that, Chapter 0 lets you play with that threshold live.

RFM: the final segmentation

Chapter 4 swept K from 2 to 8 by silhouette. K=2 won by a wide margin (0.743), but that was just the trivial cut Chapter 3 had already found without any clustering: bought more than once or not. Excluding K=2, the next best by silhouette was K=4 (0.49, higher than 5, 6, 7 and 8), so that's what I use for the real segmentation, across Chapter 3's 95,560 unique customers.

data_max = orders['order_purchase_timestamp'].max()
orders_validos = orders[orders['order_status'] != 'canceled']
oc = orders_validos.merge(customers, on='customer_id', how='left')
pagamento_por_pedido = payments.groupby('order_id')['payment_value'].sum().reset_index()
oc = oc.merge(pagamento_por_pedido, on='order_id', how='left')

rfm = oc.groupby('customer_unique_id').agg(
    recencia_dias=('order_purchase_timestamp', lambda x: (data_max - x.max()).days),
    frequencia=('order_id', 'nunique'),
    valor_monetario=('payment_value', 'sum'),
).reset_index().dropna()

scaler = StandardScaler()
rfm_escalado = scaler.fit_transform(rfm[['recencia_dias', 'frequencia', 'valor_monetario']])

kmeans_final = KMeans(n_clusters=4, random_state=42, n_init=10)
rfm['cluster'] = kmeans_final.fit_predict(rfm_escalado)

perfil = rfm.groupby('cluster').agg(
    recencia_media=('recencia_dias', 'mean'),
    frequencia_media=('frequencia', 'mean'),
    valor_medio=('valor_monetario', 'mean'),
    n_clientes=('customer_unique_id', 'count'),
).reset_index()

I rebuild RFM identically to Chapter 4, but now with K fixed at 4 instead of sweeping 2 through 8: kmeans_final.fit_predict(rfm_escalado) trains and already returns each customer's cluster in one pass, becoming a new column (rfm['cluster']). The .groupby('cluster').agg(...) that follows is what builds the profile table below, the average of each RFM variable within each cluster plus the customer count.

Naming each one by what the numbers actually show: cluster 0 is the recent entry-level customer (bought once, recently, low ticket), cluster 2 is the dormant customer (bought once, almost a year and two months ago, same low ticket as cluster 0), cluster 1 is the occasional big spender (ticket nearly 9 times higher than the previous two, but almost never comes back), and cluster 3 is the genuine repeat customer, the only one with average frequency above 1 (2.114 orders), matching the 3.06% of repeat buyers Chapter 3 already found without any clustering at all.

Here's the chapter's real finding: clusters 0 and 2 together are 94.3% of the entire base, and the only thing separating them is recency, average ticket and frequency are practically identical between them (R$ 134.58 vs. R$ 133.72, both at 1.000 frequency). K=4 is statistically the best non-trivial choice, but in practice it only delivers two segments that actually matter to the business (the occasional big spender and the repeat customer), the rest is a fancier way of saying "bought recently" versus "bought a long time ago". A statistically optimal K isn't synonymous with a rich segmentation, you have to look at each cluster's centroid and ask whether the difference actually matters for some business decision.

Preparing the interactive showcase

Three more exports that don't produce a chart here in the chapter, but feed Chapter 0's interactive widgets: an RFM sample with each customer's cluster, the threshold sweep on the two models that actually decide something, and the region → delivery → score flow.

amostra = rfm.sample(n=min(500, len(rfm)), random_state=42)[['recencia_dias', 'frequencia', 'valor_monetario', 'cluster']]

Just a sample of 500 customers, running the 3D scatter with nearly 95 thousand points would choke the browser without gaining anything in visual clarity, keeping the three original RFM variables plus each customer's cluster.

modelos_sweep = {
    'Random Forest': RandomForestClassifier(n_estimators=100, max_depth=8, random_state=42),
    'XGBoost': xgb.XGBClassifier(n_estimators=100, max_depth=6, random_state=42, eval_metric='logloss'),
}

sweep = []
for nome, modelo in modelos_sweep.items():
    modelo.fit(X1tr, y1tr)
    proba = modelo.predict_proba(X1te)[:, 1]
    for limiar in np.arange(0.0, 1.01, 0.02):
        pred = (proba >= limiar).astype(int)
        sweep.append({
            'model': nome,
            'threshold': round(float(limiar), 2),
            'precision': round(float(precision_score(y1te, pred, zero_division=0)), 4),
            'recall': round(float(recall_score(y1te, pred, zero_division=0)), 4),
            'f1': round(float(f1_score(y1te, pred, zero_division=0)), 4),
        })

I retrain Random Forest and XGBoost, the two that actually bet on a prediction, and sweep np.arange(0.0, 1.01, 0.02), 51 evenly spaced threshold values. For each one, (proba >= limiar).astype(int) recomputes the binary decision at that specific cut and recomputes the three metrics: that's the raw data behind the slider you drag in Chapter 0, each slider position is, literally, a different row of that sweep list.

regiao_por_uf = {
    'AC': 'Norte', 'AP': 'Norte', 'AM': 'Norte', 'PA': 'Norte', 'RO': 'Norte', 'RR': 'Norte', 'TO': 'Norte',
    'AL': 'Nordeste', 'BA': 'Nordeste', 'CE': 'Nordeste', 'MA': 'Nordeste', 'PB': 'Nordeste', 'PE': 'Nordeste', 'PI': 'Nordeste', 'RN': 'Nordeste', 'SE': 'Nordeste',
    'DF': 'Centro-Oeste', 'GO': 'Centro-Oeste', 'MT': 'Centro-Oeste', 'MS': 'Centro-Oeste',
    'ES': 'Sudeste', 'MG': 'Sudeste', 'RJ': 'Sudeste', 'SP': 'Sudeste',
    'PR': 'Sul', 'RS': 'Sul', 'SC': 'Sul',
}

sk = m.drop_duplicates(subset='order_id')[['order_id', 'customer_state', 'atraso_dias', 'review_score']].dropna()
sk['regiao'] = sk['customer_state'].map(regiao_por_uf)
sk['status_entrega'] = np.where(sk['atraso_dias'] > 0, 'Atrasado', 'No prazo')
sk['faixa_nota'] = pd.cut(sk['review_score'], bins=[0, 2, 3, 5], labels=['Nota baixa (1-2)', 'Nota média (3)', 'Nota alta (4-5)'])

fluxo1 = sk.groupby(['regiao', 'status_entrega']).size().reset_index(name='value')
fluxo2 = sk.groupby(['status_entrega', 'faixa_nota'], observed=True).size().reset_index(name='value')

And the sankey flow in two steps: first I map each state to one of Brazil's five regions (regiao_por_uf), then I build three categorical columns (region, whether it was late or not, score band) and count how many orders fall into each region→status and status→band combination. Two separate counts (fluxo1, fluxo2) that, together on the chart, form the whole region → delivery → score funnel you see in Chapter 0.

Wrapping up the project

The four scenarios, summed up by what they actually deliver: delay has an honest AUC (0.74) but no model ready for automated decisions without tuning the threshold first. Review score runs into a low ceiling because the only genuinely informative features (delay, delivery time) were already exhausted in Chapter 4, leaving little signal for any model to exploit. Freight is the most "solved" scenario of all, XGBoost with R² 0.64 is already usable for a real estimate. RFM has four clusters on paper, but only two actually differentiate anyone (together, just 5.6% of the base), the rest is recency noise dressed up as segmentation.

Limitations worth noting: the dataset is a frozen snapshot from 2016 to 2018, no continuous updates. There's no real churn data, a "dormant" customer might have simply stopped buying, or migrated to another marketplace, the dataset can't tell the two apart. SHAP only covers the two most narratively rich scenarios (delay and review), an editorial choice, not a technical limitation. No model went through hyperparameter search, the n_estimators/max_depth numbers are reasonable choices, not an optimum.

If you want to relive all of this without reading another line of prose, Chapter 0 is the showcase: the same threshold that only showed up as text here becomes a slider you drag by hand, RFM becomes a rotatable 3D scatter, and a handful more ways of looking at the same data that ran through these five chapters. And if you want one more technical angle, Chapter 6 closes the project with the one question the four original scenarios never touched: how many orders are coming next.