Building an engagement score with RFM features, genetic-algorithm binning and least-squares weights

The business wanted to find its true fans and grow their number. Everyone agreed such customers existed; no single number identified them. Service usage, subscription count and recent activity each missed something: the heaviest data users were not the most loyal, and the longest-tenured customers were not always the most active. Counting usage across services rewards whichever service everyone uses daily and says little about loyalty. Before any model, the concept needed a definition the business would accept — and a way to say which services, how much and in what combination make a fan.

In this post we walk through the score we built. Engagement is defined by its consequences: an engaged customer pays more (ARPU) and leaves less (churn). Each service's use is described with RFM variables — recency, frequency and intensity. Optimal binning, with cut points searched by a genetic algorithm, turns each skewed variable into a few levels that best separate the target. Least-squares regression on the levelled variables gives a weight per service and facet, and from them a score and a rank for every customer. We close with the related rank-aggregation approach, used for target marketing where there was no target to fit.

Solution overview

The score is fitted on a training period and then applied as a lookup. Fitting defines the target, computes RFM variables per service, bins and encodes each variable, and regresses the target on the encodings. Scoring takes each customer's levels, each level's value, a weighted sum and a percentile — so every customer's score and rank can be explained level by level.

The numbered steps in the diagram:

1. Engagement is defined in relation to ARPU and churn: an engaged customer pays more and leaves less, combined into one target.
2. For each service, recency, frequency and intensity variables describe how each customer uses it; non-users get a level of their own.
3. Optimal binning cuts each variable into a few levels at the points that best separate the target, and each level is encoded by the target's mean in it.
4. A genetic algorithm searches cut-point combinations by selection, crossover and mutation over candidate thresholds.
5. Least-squares regression on the encoded levels gives a weight per service and facet; partial least squares handles collinear facets.
6. The weighted sum becomes a 0–100 score and a rank, and the weights show which services contribute most to customer value.

Technology stack

Step 1: Define engagement as a target

Most of the project's value came before any modelling: agreeing that engagement means higher revenue and lower churn, in combination. That turned a feeling into a target that every later choice — which variables, which cut points, which weights — could be tested against. A score that merely counted activity would have rewarded the service everyone uses daily and said nothing about loyalty.

y_i = z(ARPU_i) − λ · churn_i

z(·)      standardization over the training customers
churn_i   1 if the customer left within the outcome window, else 0
λ         how much a departure costs, in standard deviations of monthly revenue

λ is a business choice rather than a statistical one, and it changes the ranking: a larger λ favours customers who stay over customers who pay. The live model uses λ = 1.2 and labels churn over the following six months.

Why define the target first? Without it there is nothing to bin against and nothing to weight by, and every design choice becomes a matter of opinion. With it, each choice is a measurable improvement or not.

Step 2: Describe each service with RFM variables

Engagement shows up differently in each service, so usage is described per service with three facets: recency (days since last use), frequency (uses in a window) and intensity, the monetary-style facet — hours watched, data used or amount spent. With six services, as in the live model, that is eighteen variables per customer.

Non-users are not the low end of the users' scale: not using IPTV at all is a different state from watching it once a month. They get a level of their own (level 0), and the cut points are fitted on users only.

import pandas as pd

def rfm(events, customers, asof, window_days=90):
    """events: one row per use (customer, service, ts, amount).
    Returns one row per customer with <service>_R, _F, _M columns."""
    recent = events[events.ts > asof - pd.Timedelta(days=window_days)]
    g = recent.groupby(["customer", "service"])
    out = pd.DataFrame({
        "R": (asof - g.ts.max()).dt.days,    # recency: days since last use
        "F": g.size(),                        # frequency: uses in the window
        "M": g.amount.sum(),                  # intensity: hours, GB or spend
    }).unstack("service")
    out.columns = [f"{svc}_{facet}" for facet, svc in out.columns]
    return out.reindex(customers)             # non-users stay NaN -> level 0

Simplified. events is a generic usage log; the real services, windows and definitions are not shown.

Step 3: Cut usage into levels with optimal binning

Raw usage is skewed and its effect is non-linear: the tenth login in a month means less than the first, and some services have a threshold above which behaviour changes. Each RFM variable is therefore discretized into a few ordered levels, and optimal binning chooses the cut points that make those levels separate the target as well as possible.

The objective in the live model is the correlation ratio η² — the share of the target's variance explained by the level means — with a penalty for levels that hold too few customers to be stable.

level(x; c) = 1 + #{ k : c_k < x }          c = (c_1 < c_2 < c_3): four levels for users, level 0 for non-users

η²(c) = Σ_L n_L (ȳ_L − ȳ)² / Σ_i (y_i − ȳ)²  −  ρ · #{ L : n_L < ε·n }

c*    = argmax_c η²(c)          candidates: the 5%, 10%, …, 95% quantiles of x

import numpy as np

def levels(x, cuts):
    """Levels 1 … len(cuts)+1 for users; level 0 is reserved for non-users."""
    lv = np.searchsorted(np.sort(cuts), x, side="left") + 1
    return np.where(np.isnan(x), 0, lv)

def eta_squared(y, lv, min_share=0.05, penalty=0.05):
    """Correlation ratio: share of the target's variance explained by level means,
    minus a penalty for each level too thin to be stable."""
    total = ((y - y.mean()) ** 2).sum()
    between, thin = 0.0, 0
    for L in np.unique(lv):
        yl = y[lv == L]
        between += len(yl) * (yl.mean() - y.mean()) ** 2
        thin += len(yl) < min_share * len(y)
    return between / total - penalty * thin

def level_means(lv, y):
    """Encode each level by the target's mean in it (training customers)."""
    return {L: y[lv == L].mean() for L in np.unique(lv)}

Simplified objective and encoding, as in the live model: three cuts give four user levels, and levels with under 5% of customers are penalized.

Each level is then encoded by the target's mean among the training customers in it. The encoded value is what the level predicts, which is what keeps the final score readable: a level's contribution is a number on the target's own scale.

Step 4: Search the cut points with a genetic algorithm

Cut points interact — moving one changes which customers the others separate — so they are searched jointly. A genetic algorithm treats a set of cut points as a genome and evolves a population of them: score each genome by the objective, keep the best few unchanged (elitism), breed children from well-ranked parents by crossover, and mutate a child by nudging one cut to a neighbouring candidate. After a few dozen generations the population converges on a good combination without enumerating them all.

import numpy as np

def ga_cuts(x, y, n_cuts=3, pop=24, gens=30, elite=4, p_mut=0.6, seed=0):
    """Search cut points among the 5%…95% quantiles of x to maximize η²(y)."""
    rng = np.random.default_rng(seed)
    cand = np.quantile(x, np.linspace(0.05, 0.95, 19))
    fit = lambda g: eta_squared(y, levels(x, cand[g]))
    P = [np.sort(rng.choice(len(cand), n_cuts, replace=False)) for _ in range(pop)]
    for _ in range(gens):
        P.sort(key=fit, reverse=True)                        # rank by fitness
        nxt = P[:elite]                                      # elitism
        while len(nxt) < pop:
            a, b = (P[min(rng.integers(8), rng.integers(8))] for _ in range(2))
            child = np.where(rng.random(n_cuts) < 0.5, a, b)        # uniform crossover
            if rng.random() < p_mut:                                 # mutation
                k = rng.integers(n_cuts)
                child[k] = np.clip(child[k] + round(rng.normal(0, 2)), 0, len(cand) - 1)
            child = np.unique(child)                         # sorted, no duplicates
            if len(child) == n_cuts:
                nxt.append(child)
        P = nxt
    best = max(P, key=fit)
    return cand[best], fit(best)

Simplified GA, as in the live model: 24 genomes, 30 generations, elitism of 4, parents drawn with a bias towards the top eight.

Why a genetic algorithm? For three cuts among nineteen candidates, exhaustive search is still affordable (969 combinations). The GA earns its keep when the candidate grid is finer, the number of levels grows or many variables are binned together, and it accepts any objective, including penalties that break the structure exact methods rely on. Where effects are smooth, its cuts land close to equal-count bins — the live model lets you compare the two.

Step 5: Weight services and facets with least squares

With every variable encoded as its level's target mean, the score is a linear combination of the encodings. Regressing the target on them gives one weight per service and facet; the weight says how much that facet adds once the others are known.

Least squares is the base method. RFM facets of the same service are strongly correlated, so partial least squares, which regresses on a few components chosen to maximize covariance with the target, is the alternative when facets crowd each other out. The live model uses ridge least squares, which serves the same purpose with a small penalty.

score_i = β_0 + Σ_(service s) Σ_(facet f ∈ {R, F, M})  β_sf · v_sf( level_sf,i )

v_sf(L) = mean of y over training customers at level L of variable (s, f)
β       = argmin ‖y − Vβ‖² + α‖β‖²          α = 0 for plain least squares

import numpy as np
from sklearn.linear_model import Ridge
from sklearn.cross_decomposition import PLSRegression

def encode(levels_by_var, means_by_var):
    """One column per (service, facet): the target mean of the customer's level."""
    cols = [np.array([means_by_var[v].get(L, 0.0) for L in lv])
            for v, lv in levels_by_var.items()]
    return np.column_stack(cols)

V_train = encode(train_levels, means)                 # customers × (services × R, F, M)
ridge = Ridge(alpha=2.0).fit(V_train, y_train)        # intercept is not penalized
pls = PLSRegression(n_components=4).fit(V_train, y_train)   # when facets crowd each other

raw_train = ridge.predict(V_train)
raw_test = ridge.predict(encode(test_levels, means))

Simplified. Ridge does not penalize the intercept; PLSRegression is shown as the alternative for collinear facets.

Step 6: Turn weights into a score, a rank and service priorities

The raw linear score is mapped to 0–100 as its percentile among the training customers, so a score of 82 means more engaged than 82% of the reference base, and the rank follows directly.

Service weights summarize which services make a fan. Raw coefficients are not comparable across variables with different spreads, so each variable's contribution is measured as |β| times the standard deviation of its encoded column, summed per service and normalized to shares.

import numpy as np
from sklearn.metrics import roc_auc_score

def to_score(raw, raw_train):
    """0–100: percentile of a raw score among the training customers."""
    return np.round(100 * np.searchsorted(np.sort(raw_train), raw) / len(raw_train))

def service_weights(beta, V, var_service):
    """Share of the score's spread carried by each service: Σ |β_k| · sd(V_k)."""
    spread = np.abs(beta) * V.std(axis=0)
    shares = {}
    for k, svc in enumerate(var_service):
        shares[svc] = shares.get(svc, 0.0) + spread[k]
    total = sum(shares.values())
    return {svc: v / total for svc, v in shares.items()}

def deciles(score, churned, arpu):
    order = np.argsort(-score)                            # top decile first
    return [(churned[i].mean(), arpu[i].mean()) for i in np.array_split(order, 10)]

score = to_score(raw_test, raw_train)
weights = service_weights(ridge.coef_, V_train, var_service)
auc = roc_auc_score(churned_test, -score)                 # higher score, less churn

Simplified, as in the live model. roc_auc_score on the negated score measures how well a higher score picks out customers who stay.

On held-out customers, the live model compares the score with the obvious baseline — ranking by raw usage count — using churn and ARPU by decile and the AUC for churn. Counting uses rewards mobile data, which in the generated base follows plan and price sensitivity rather than loyalty; the score weights it only by what it actually predicts.

Step 7: Aggregate ranks when there is no target to fit

A related target-marketing project for LG Hausys faced a different version of the same question. Interest indices for marriage, moving and interior remodelling had to be built from several sources — related web and app activity, related product purchases — and each source produced its own ranking of customers. The sources disagreed, had no common scale, and there was no outcome to calibrate against.

There, rank aggregation combined the rankings directly. Methods such as the Borda count and hierarchical rank aggregation need only each source's order, not its scale. In the Borda count each customer collects points for their position in every source, and the totals give the combined order.

Borda(i) = Σ_(source k) ( N − rank_k(i) )          rank 1 = best,  N customers

import numpy as np

def borda(ranks):
    """ranks: sources × customers, 0 = best; NaN where a source has no signal.
    Each customer collects (N − 1 − rank) points from every source."""
    R = np.asarray(ranks, dtype=float)
    N = R.shape[1]
    R = np.where(np.isnan(R), N - 1, R)                   # missing = worst rank
    points = ((N - 1) - R).sum(axis=0)
    return np.argsort(-points, kind="stable")             # aggregated order, best first

Simplified Borda count; customers missing from a source get its worst rank.

Scores or ranks? Rank aggregation is the right tool when sources cannot be put on one scale and no target exists. When a target exists, as with engagement, calibrating levels and weights against it says more — including how much each source matters.

Results

The project converted an intuitive concept of fans into an operational Engagement Score and Rank. It supported customer segmentation, service-usage campaigns and retention strategies, and the prioritisation of services that contribute to customer value.

Because every customer's score decomposes into levels and weights, a campaign can aim at a specific move — from a lower to a higher level of a service that carries weight, or from not using a service to using it — rather than at "more usage" in general.

Lessons learned

• The target does the defining. Agreeing that engagement means higher revenue and lower churn made every later choice — variables, cut points, weights — testable.
• Levels make a score explainable. Binning costs a little precision and buys a plain meaning for each level and a contribution the business can read.
• Search hard only where thresholds matter. A genetic search finds thresholds the data knows and the analyst does not; where effects are smooth it lands near equal-count bins at little cost.
• Report weights on a common scale. |β| times the spread of each encoded variable makes service weights comparable and turns the model into a statement about which services to grow.

Conclusion

Defining engagement against ARPU and churn, describing usage with per-service RFM variables, cutting them into levels with genetic-algorithm-searched optimal binning and weighting the levels by least squares turned an intuitive idea of fans into a score that ranks every customer and says which services contribute to customer value.

The pattern — fix the target first, discretize against it, weight with a transparent linear model — fits any composite index that has to be both predictive and explainable, such as loyalty, health or risk scores. Where no target exists, rank aggregation is the fallback.

Limitations

• The score inherits its definition: a different weighting of revenue against churn produces a different ranking, and that choice belongs to the business.
• Levels fitted on one period drift as services and plans change; cut points and weights need periodic refitting.
• In the live model the engagement structure is planted and smooth, which flatters both the usage count and the score, and the gain from optimized binning there is modest; its weights come from ridge least squares on level means.

About the demo and confidentiality

Customers, services, usage, revenue and churn in the embedded model are generated. No usage, billing or churn data, service list, cut point or weight from the real project appears in this post. Code is simplified and written for illustration; the libraries it uses are not a description of the production stack.