Executive Summary
Predicted customer value over the next 365 days
The model projects 17,155 in total value over the next 365 days—an average of 14.3 per customer, but a median of only 0.33. The top 10% account for 83.5% of all projected value. Validation on held-out data (the final 219 days) was good: the model predicted 733.6 repeat purchases against 648 actually observed, correlating 0.67 with individual customer behaviour. These are forward-looking projections under the assumption that past purchasing patterns continue and do not account for costs, margin, or discounting.
Analysis Overview
BG/NBD + Gamma-Gamma lifetime value across 1,200 customers and 6,259 transactions.
Two probabilistic models work together to forecast customer value: BG/NBD separates quiet customers (who buy rarely) from churned ones by fitting each customer's latent purchase rate and dropout probability; Gamma-Gamma estimates what each transaction is worth by blending individual and population averages. Fitted on 6,259 transactions across 1,200 customers over 729 days, the models produce a 365-day forward projection. This is a forecast assuming past behaviour persists; it excludes costs, margin, and discounting. The repeat-customer base (1,122 customers, 94% of the total) provides the only signal for churn detection.
Data Quality
Row accounting, date parsing, and how purchase occasions were formed.
6,271 rows loaded; 6,259 usable after dropping 12 unparseable dates. The data spans 2024-01-01 to 2025-12-30 (729 days) and covers 1,200 customers. Same-day purchases by the same customer were merged into single purchase occasions (212 transactions consolidated). Customer age (T) is measured from each customer's first purchase to 2025-12-30; recency is the date of their most recent purchase. These two quantities feed the BG/NBD model. Customers appearing only late in the window have minimal history, so their individual forecasts rely more heavily on the population estimate.
How Predicted Value Is Distributed
Projected 365-day value across every customer.
Customer value is heavily right-skewed: the mean (14.3) sits far above the median (0.33), only 15.3% of customers are worth more than the average, and 88.6% are projected at effectively zero because the model judges them very likely already churned. This concentration is the practical point: setting acquisition or retention budgets against the mean CLV will systematically overpay for the tail and underpay for the head. The distribution shows that almost no customer is worth the average, making segment-based strategies essential for efficient spend allocation.
Highest-Value Customers
Top 20 by predicted 365-day value, with P(alive).
| Customer | Predicted Clv | P Alive | Expected Transactions | Expected Avg Value | Past Purchases |
|---|---|---|---|---|---|
| CUST-01071 | 1733 | 0.934 | 27.19 | 63.72 | 59 |
| CUST-00659 | 881.9 | 0.943 | 14.99 | 58.83 | 28 |
| CUST-00258 | 611.3 | 0.85 | 18.82 | 32.48 | 44 |
| CUST-00475 | 534.3 | 0.899 | 9.5 | 56.26 | 23 |
| CUST-00978 | 445.3 | 0.726 | 5.86 | 76.02 | 18 |
| CUST-00637 | 418.6 | 0.931 | 7.28 | 57.48 | 18 |
| CUST-00573 | 410.9 | 0.938 | 13.22 | 31.07 | 32 |
| CUST-00309 | 377 | 0.253 | 6.72 | 56.12 | 48 |
| CUST-00749 | 366.6 | 0.805 | 6.53 | 56.17 | 16 |
| CUST-00758 | 306.1 | 0.933 | 14.05 | 21.79 | 31 |
| CUST-00887 | 304.4 | 0.978 | 22.81 | 13.35 | 60 |
| CUST-00189 | 299.1 | 0.96 | 18.33 | 16.31 | 42 |
| CUST-00151 | 282.4 | 0.757 | 3.84 | 73.55 | 10 |
| CUST-00019 | 249.9 | 0.912 | 9.12 | 27.41 | 20 |
| CUST-00009 | 234.6 | 0.819 | 3.14 | 74.61 | 9 |
| CUST-00043 | 220.9 | 0.926 | 6.72 | 32.89 | 18 |
| CUST-00747 | 219.9 | 0.88 | 6.46 | 34.05 | 16 |
| CUST-00650 | 188.8 | 0.709 | 9.88 | 19.12 | 25 |
| CUST-01173 | 179 | 0.739 | 2.15 | 83.26 | 5 |
| CUST-00328 | 169.9 | 0.941 | 13.06 | 13.01 | 30 |
The top 20 customers account for 49.2% of the 17,155 total projected value. Nineteen of the 20 have P(alive) above 0.5, indicating the model judges them likely still active. CUST-01071 leads at 1,732.71, driven by 27.19 expected transactions at an average value of 63.72 per transaction. The notable exception is CUST-00309 (projected value 376.97) with P(alive) of only 0.253—a reactivation case: high value if they return, but the model doubts they will on their own. Each customer's projection is split between expected transaction count and expected per-transaction value, showing which driver dominates their forecast.
Does the Forecast Hold Up?
Predicted vs actual repeat purchases over 219 held-out days.
The model was refit on the first 510 days of history, then asked to predict repeat purchases over the next 219 days it had never seen. Correlation with actual behaviour is 0.668, with a mean absolute error of 0.577 purchases per customer. In aggregate, the model predicted 733.6 purchases against 648 actually observed—a 13% over-prediction. This rates as good validation. Vertical scatter is expected because individual purchase timing is random even when the underlying rate is correct; the aggregate total and correlation are the valid measures of forecast quality, not individual point accuracy.
Fitted Model Parameters
Maximum-likelihood estimates and fit diagnostics for both models.
| Parameter | Value | Interpretation |
|---|---|---|
| r (purchase-rate shape) | 0.7754 | Spread of the underlying purchase rate across customers |
| alpha (purchase-rate scale) | 13.67 | Purchase rate is measured in transactions per day of Transaction Date |
| a (dropout shape) | 1.33 | Shape of the churn-probability distribution |
| b (dropout shape) | 2.441 | Second shape of the churn-probability distribution |
| mean purchases per day | 0.0567 | Model-implied average rate: one purchase every 18 days |
| mean dropout probability | 0.3526 | Model-implied chance a customer stops after any given purchase |
| p (value shape) | 5.195 | Spread of transaction values within a customer |
| q (value shape) | 4.067 | Spread of average transaction value between customers (q>1 required) |
| gamma (value scale) | 17.93 | Scale of the value distribution, in units of Transaction Amount |
| population mean transaction value | 30.38 | What a randomly chosen customer's typical Transaction Amount is worth |
| BG/NBD log-likelihood | -2.25e+04 | Maximised log-likelihood of the transaction/churn fit |
| Gamma-Gamma log-likelihood | -4645 | Maximised log-likelihood of the monetary fit |
| frequency vs value correlation | 0.075 | Near zero, as the Gamma-Gamma model assumes |
| customers fitted | 1200 | Distinct values of Customer ID used in the fit |
| repeat customers | 1122 | Customers with more than one purchase — the only ones informing churn |
BG/NBD converged to r = 0.7754 and alpha = 13.6709 (purchase-rate parameters); a = 1.3296 and b = 2.4408 (dropout parameters). These imply a mean purchase rate of 0.0567 per day (roughly one purchase every 18 days for an active customer) and a 35.3% dropout probability after any given purchase. Gamma-Gamma fit with p = 5.1954, q = 4.0674, and gamma = 17.9345, yielding a population mean transaction value of 30.38. Repeat frequency and average transaction value correlate at r = 0.075, validating the independence assumption. The BG/NBD log-likelihood is −22,497.1 and Gamma-Gamma is −4,645.2; both models converged successfully on 1,200 customers, of which 1,122 had repeat purchases.
Methodology
Statistical methodology and diagnostics for Predictive Customer Lifetime Value (BG/NBD + Gamma-Gamma)
Statistical Method
Standard-library analysis: forecast what your customers are worth from here, straight from a raw transactions table (one row per order). The industry-standard buy-till-you-die model — BG/NBD for how often each customer buys and how likely they are to have quietly churned, Gamma-Gamma for what a transaction of theirs is worth — gives every customer an expected number of future purchases, a probability they are still active, an expected transaction value, and a predicted lifetime value over a stated horizon. Every run is back-tested: part of your own history is held out, the model is refit without it, and predicted purchases are compared against what actually happened, so you can see whether the forecast is worth trusting before you act on it.
- One row per transaction (transaction-level data, not pre-aggregated per customer)
- Purchases while a customer is active follow a Poisson process at that customer's own rate; rates vary across customers as a Gamma distribution
- A customer may churn permanently after any repeat purchase, with a probability that varies across customers as a Beta distribution
- Transaction values are Gamma-distributed within a customer and independent of that customer's purchase frequency (this assumption is measured and reported, not assumed silently)
- The last date in the data is the end of observation; past purchasing behaviour continues into the horizon
- This is a model-based projection, not a guarantee — it assumes the future resembles the observed past, and a business change (pricing, product, market) invalidates it
- Value is gross transaction value: costs, margin, and discounting are NOT applied unless the amount column you map already reflects margin
- The horizon is not open-ended lifetime value — it is the stated number of days, chosen from your data's own span and capped at one year
- The model is non-contractual: it infers churn from silence, so it cannot see a cancellation you already know about
Analysis Code
Complete R source code for this analysis
Predictive Customer Lifetime Value — BG/NBD + Gamma-Gamma
Fits the standard buy-till-you-die model to a raw transaction log and projects each customer forward: how many more purchases they are likely to make, how likely they are to still be active at all, what a typical transaction of theirs is worth, and therefore what they are worth over a stated horizon.
Why This Method?
RFM tells you what a customer HAS done — it is a backward-looking segmentation of observed behaviour. BG/NBD + Gamma-Gamma is the forward-looking counterpart: it estimates the latent purchase rate and the latent churn probability behind that behaviour, so a customer who is quiet because they buy twice a year is separated from one who is quiet because they are gone. The two are complements — segment with RFM, forecast with this.
What This Analysis Covers
- BG/NBD fit by maximum likelihood (transaction + dropout process)
- Gamma-Gamma fit by maximum likelihood (monetary value)
- Per-customer expected transactions, P(alive), expected value, and CLV
- A calibration/holdout back-test: predicted vs actual repeat purchases
- An explicit check of the frequency/monetary independence assumption
Standard Library
Platform standard-library module (LAT-1545): runs on ANY transaction log via the semantic mapping {customer, date, amount}. All narrative is derived from the user's own column names and computed values. The model is fit here by hand (optim on the exact log-likelihoods) — no BTYD package.
suppressPackageStartupMessages(library(DT))
suppressPackageStartupMessages(library(htmlwidgets))
suppressPackageStartupMessages(library(arrow))
suppressPackageStartupMessages(library(knitr))
suppressPackageStartupMessages(library(rmarkdown))
suppressPackageStartupMessages(library(dplyr))
suppressPackageStartupMessages(library(tidyr))
suppressPackageStartupMessages(library(ggplot2))
suppressPackageStartupMessages(library(stringr))
suppressPackageStartupMessages(library(lubridate))
suppressPackageStartupMessages(library(broom))
suppressPackageStartupMessages(library(Matrix))
suppressPackageStartupMessages(library(cluster))
suppressPackageStartupMessages(library(data.table))Core Analysis Pipeline
Multi-format date parser
Tries each candidate format across the whole vector, keeps the format that parses the most values, then back-fills stragglers from the remaining formats. Returns a Date vector (NA where unparseable).
parse_dates_multi <- function(x) {
if (inherits(x, "Date")) return(x)
if (inherits(x, "POSIXct") || inherits(x, "POSIXlt")) return(as.Date(x))
s <- trimws(as.character(x))
s[is.na(s) | s == ""] <- NA_character_
fmts <- c("%Y-%m-%d", "%Y/%m/%d", "%m/%d/%Y", "%d/%m/%Y", "%m-%d-%Y",
"%d-%m-%Y", "%d.%m.%Y", "%b %d, %Y", "%B %d, %Y",
"%d %b %Y", "%d %B %Y")
parsed <- lapply(fmts, function(f) suppressWarnings(as.Date(s, format = f)))
counts <- vapply(parsed, function(p) sum(!is.na(p)), integer(1))
best <- parsed[[which.max(counts)]] # counts is integer, never all-NA
for (p in parsed[order(-counts)]) {
fill <- is.na(best) & !is.na(p)
if (any(fill)) best[fill] <- p[fill]
}
best
}log(exp(u) + exp(v)) without overflow
v is allowed to be -Inf (the x == 0 branch of the BG/NBD likelihood).
log_sum_exp2 <- function(u, v) {
m <- pmax(u, v)
out <- rep(-Inf, length(m))
ok <- is.finite(m)
out[ok] <- m[ok] + log(exp(u[ok] - m[ok]) + exp(v[ok] - m[ok]))
out
}Gaussian hypergeometric 2F1 by its defining power series
2F1(a,b;c;z) = sum_n (a)_n (b)_n / (c)_n * z^n / n!. Convergent for |z| < 1, which always holds here because z = t / (alpha + T + t) with all terms positive. Vectorised over z; iterates until every element stops changing. Returns NA for elements that fail to settle so callers can refuse rather than silently propagate a truncated series.
h2f1_series <- function(a, b, c, z, max_iter = 5000) {
n <- length(z)
term <- rep(1, n)
total <- rep(1, n)
settled <- rep(FALSE, n)
for (j in seq_len(max_iter)) {
prev <- total
term <- term * (a + j - 1) * (b + j - 1) / ((c + j - 1) * j) * z
total <- total + term
settled <- (total == prev)
if (all(settled)) break
}
total[!settled] <- NA_real_
total
}BG/NBD negative log-likelihood (Fader, Hardie & Lee 2005)
Parameters are carried on the log scale so optim searches an unconstrained space and every parameter stays strictly positive. x = number of REPEAT transactions t_x = age at the last transaction (recency, same units as T) T = total observed age (first transaction -> end of observation)
bgnbd_nll <- function(par_log, x, t_x, T_obs) {
pr <- exp(par_log)
if (any(!is.finite(pr)) || any(pr <= 0) || any(pr > 1e6)) return(1e12)
r <- pr[1]; alpha <- pr[2]; a <- pr[3]; b <- pr[4]
ln_a1 <- lgamma(r + x) - lgamma(r) + r * log(alpha)
ln_a2 <- lgamma(a + b) + lgamma(b + x) - lgamma(b) - lgamma(a + b + x)
ln_a3 <- -(r + x) * log(alpha + T_obs)
ln_a4 <- rep(-Inf, length(x))
pos <- x > 0
if (any(pos)) {
ln_a4[pos] <- log(a) - log(b + x[pos] - 1) -
(r + x[pos]) * log(alpha + t_x[pos])
}
ll <- ln_a1 + ln_a2 + log_sum_exp2(ln_a3, ln_a4)
if (any(!is.finite(ll))) return(1e12)
-sum(ll)
}Fit BG/NBD by maximum likelihood
Multi-start Nelder-Mead (the surface is flat and ridged, so a single start is not trustworthy) refined with BFGS. Returns ok = FALSE with a reason whenever the optimiser fails or the estimate lands somewhere the model cannot be believed at — the caller must refuse rather than degrade.
fit_bgnbd <- function(x, t_x, T_obs) {
mean_rate <- mean(x / pmax(T_obs, 1e-6))
alpha_hint <- if (is.finite(mean_rate) && mean_rate > 0) 1 / mean_rate else 10
starts <- list(
log(c(1, 1, 1, 1)),
log(c(0.5, max(alpha_hint, 1e-3), 1.2, 2.5)),
log(c(1, max(alpha_hint / 4, 1e-3), 2, 4)),
log(c(2, 10, 1.5, 1.5)),
log(c(0.25, 50, 0.8, 3))
)
best <- NULL
for (s in starts) {
fit <- tryCatch(
optim(s, bgnbd_nll, x = x, t_x = t_x, T_obs = T_obs,
method = "Nelder-Mead",
control = list(maxit = 20000, reltol = 1e-12)),
error = function(e) NULL)
if (is.null(fit) || !is.finite(fit$value)) next
if (is.null(best) || fit$value < best$value) best <- fit
}
if (is.null(best)) {
return(list(ok = FALSE,
reason = "the optimiser could not evaluate the BG/NBD likelihood at any starting point"))
}
refined <- tryCatch(
optim(best$par, bgnbd_nll, x = x, t_x = t_x, T_obs = T_obs,
method = "BFGS", control = list(maxit = 5000, reltol = 1e-12)),
error = function(e) NULL)
if (!is.null(refined) && is.finite(refined$value) && refined$value <= best$value) {
best <- refined
}
pr <- exp(best$par)
names(pr) <- c("r", "alpha", "a", "b")
if (best$convergence != 0) {
return(list(ok = FALSE, params = pr, convergence = best$convergence,
reason = sprintf("optim did not converge(code %d)", best$convergence)))
}
if (any(!is.finite(pr)) || any(pr < 1e-4) || any(pr > 1e4)) {
return(list(ok = FALSE, params = pr, convergence = best$convergence,
reason = paste0("the fitted parameters ran to an implausible boundary(",
paste(sprintf("%s=%.4g", names(pr), pr), collapse = ", "), ")")))
}
list(ok = TRUE, params = pr, loglik = -best$value,
convergence = best$convergence, reason = "")
}P(customer is still active | x, t_x, T)
bgnbd_palive <- function(pm, x, t_x, T_obs) {
r <- pm[["r"]]; alpha <- pm[["alpha"]]; a <- pm[["a"]]; b <- pm[["b"]]
extra <- rep(0, length(x))
pos <- x > 0
if (any(pos)) {
extra[pos] <- (a / (b + x[pos] - 1)) *
((alpha + T_obs[pos]) / (alpha + t_x[pos]))^(r + x[pos])
}
1 / (1 + extra)
}E[transactions in the next horizon | x, t_x, T] (FHL 2005, eq. 10)
Returns NA where the hypergeometric series or the (a-1) denominator make the closed form untrustworthy; the caller refuses on any NA.
bgnbd_expected_transactions <- function(pm, x, t_x, T_obs, horizon) {
r <- pm[["r"]]; alpha <- pm[["alpha"]]; a <- pm[["a"]]; b <- pm[["b"]]
if (abs(a - 1) < 1e-6) return(rep(NA_real_, length(x)))
cc <- a + b + x - 1
if (any(cc <= 0)) return(rep(NA_real_, length(x)))
z <- horizon / (alpha + T_obs + horizon)
hyp <- h2f1_series(r + x, b + x, cc, z)
num <- (cc / (a - 1)) *
(1 - ((alpha + T_obs) / (alpha + T_obs + horizon))^(r + x) * hyp)
den <- 1 / bgnbd_palive(pm, x, t_x, T_obs) # == 1 + delta * (...)
out <- num / den
out[!is.finite(out)] <- NA_real_
out
}Gamma-Gamma negative log-likelihood (Fader & Hardie 2013)
Fit on customers with at least one repeat transaction, using their mean repeat transaction value mbar and repeat count x.
gg_nll <- function(par_log, x, mbar) {
pr <- exp(par_log)
if (any(!is.finite(pr)) || any(pr <= 0) || any(pr > 1e6)) return(1e12)
p <- pr[1]; q <- pr[2]; g <- pr[3]
ll <- lgamma(p * x + q) - lgamma(p * x) - lgamma(q) +
q * log(g) + (p * x - 1) * log(mbar) + p * x * log(x) -
(p * x + q) * log(g + mbar * x)
if (any(!is.finite(ll))) return(1e12)
-sum(ll)
}Fit Gamma-Gamma by maximum likelihood
q > 1 is required for the expected value to be finite; anything at or below it is an outright refusal, not a caveat.
fit_gg <- function(x, mbar) {
m0 <- mean(mbar)
starts <- list(
log(c(1, 2, max(m0, 1e-3))),
log(c(2, 3, max(2 * m0, 1e-3))),
log(c(6, 4, max(m0 / 2, 1e-3))),
log(c(0.5, 1.5, max(m0, 1e-3)))
)
best <- NULL
for (s in starts) {
fit <- tryCatch(
optim(s, gg_nll, x = x, mbar = mbar, method = "Nelder-Mead",
control = list(maxit = 20000, reltol = 1e-12)),
error = function(e) NULL)
if (is.null(fit) || !is.finite(fit$value)) next
if (is.null(best) || fit$value < best$value) best <- fit
}
if (is.null(best)) {
return(list(ok = FALSE,
reason = "the optimiser could not evaluate the Gamma-Gamma likelihood at any starting point"))
}
refined <- tryCatch(
optim(best$par, gg_nll, x = x, mbar = mbar, method = "BFGS",
control = list(maxit = 5000, reltol = 1e-12)),
error = function(e) NULL)
if (!is.null(refined) && is.finite(refined$value) && refined$value <= best$value) {
best <- refined
}
pr <- exp(best$par)
names(pr) <- c("p", "q", "gamma")
if (best$convergence != 0) {
return(list(ok = FALSE, params = pr, convergence = best$convergence,
reason = sprintf("optim did not converge(code %d)", best$convergence)))
}
if (any(!is.finite(pr)) || any(pr < 1e-4) || any(pr > 1e4)) {
return(list(ok = FALSE, params = pr, convergence = best$convergence,
reason = paste0("the fitted parameters ran to an implausible boundary(",
paste(sprintf("%s=%.4g", names(pr), pr), collapse = ", "), ")")))
}
if (pr[["q"]] <= 1.0001) {
return(list(ok = FALSE, params = pr, convergence = best$convergence,
reason = sprintf(paste0("the monetary shape parameter q estimated at %.3f; ",
"at q <= 1 the model's expected transaction value is ",
"infinite, so no value figure can be believed"),
pr[["q"]])))
}
list(ok = TRUE, params = pr, loglik = -best$value,
convergence = best$convergence, reason = "")
}E[transaction value | p, q, gamma, mbar, x]
A precision-weighted blend of the customer's own observed mean and the population mean; customers with no repeat history fall back entirely to the population mean.
gg_expected_value <- function(pm, x, mbar) {
p <- pm[["p"]]; q <- pm[["q"]]; g <- pm[["gamma"]]
pop_mean <- g * p / (q - 1)
out <- rep(pop_mean, length(x))
pos <- x > 0 & is.finite(mbar) & mbar > 0
if (any(pos)) {
den <- p * x[pos] + q - 1
out[pos] <- (g * p + mbar[pos] * p * x[pos]) / den
}
out[!is.finite(out)] <- NA_real_
out
}Build the BTYD per-customer summary (x, t_x, T, mean repeat value)
Multiple transactions on the SAME DAY are treated as one transaction, the standard BTYD convention: the models describe purchase OCCASIONS.
btyd_summary <- function(tx, obs_end) {
daily <- tx %>%
group_by(customer, txn_date) %>%
summarise(amount = sum(amount), .groups = "drop")
daily %>%
arrange(customer, txn_date) %>%
group_by(customer) %>%
summarise(
first_date = min(txn_date),
last_date = max(txn_date),
n_txn = n(),
x = n() - 1,
t_x = as.numeric(max(txn_date) - min(txn_date)),
T_obs = as.numeric(obs_end - min(txn_date)),
mean_repeat_value = if (n() > 1) {
v <- amount[txn_date > min(txn_date)]
v <- v[is.finite(v) & v > 0]
if (length(v) > 0) mean(v) else NA_real_
} else NA_real_,
total_value = sum(amount),
.groups = "drop"
) %>%
as.data.frame(stringsAsFactors = FALSE)
}
compute_shared <- function(df, params, col_map = list()) {
# === SHARED EXPORTS ===
# initial_rows/final_rows/rows_removed $ row accounting
# cust_name/date_name/amount_name $ humanized user column names
# n_bad_dates/n_missing_id/n_nonpositive/n_same_day_merged $ prep counts
# date_min/date_max/span_days/horizon_days
# n_customers/n_repeaters $ integer
# bg $ list — BG/NBD fit (ok, params r/alpha/a/b, loglik)
# gg $ list — Gamma-Gamma fit (ok, params p/q/gamma, loglik)
# cust_df $ data.frame per customer incl. expected_transactions, p_alive,
# expected_avg_value, predicted_clv
# total_clv/avg_clv/median_clv/top10_share/total_expected_txn
# fm_cor / fm_cor_p / independence_violated
# hold $ list — calibration/holdout back-test (available, cor, mae,
# predicted_total, actual_total, ratio, verdict, n)
# clv_hist_df / top_customers_df / holdout_df / model_params_df
# metrics / json_output
# === /SHARED EXPORTS ===Step 1: Required semantic columns + humanized names
initial_rows <- nrow(df)
cust_name <- humanize_semantic("customer", col_map)[1]
date_name <- humanize_semantic("date", col_map)[1]
amount_name <- humanize_semantic("amount", col_map)[1]
for (need in c("customer", "date", "amount")) {
if (!need %in% names(df)) {
stop(sprintf("Required column '%s' is not mapped.",
humanize_semantic(need, col_map)[1]))
}
}Step 2: Customer id — character, drop blank/missing
cid <- trimws(as.character(df$customer))
keep_id <- !is.na(cid) & cid != ""
n_missing_id <- sum(!keep_id)
df <- df[keep_id, , drop = FALSE]
cid <- cid[keep_id]Step 3: Dates — multi-format parse; refuse if >5% unparseable
dts <- parse_dates_multi(df$date)
n_bad_dates <- sum(is.na(dts))
if (nrow(df) > 0 && n_bad_dates / nrow(df) > 0.05) {
stop(sprintf(
paste0("Could not read %s of %s values in '%s' as dates (%.1f%%). More ",
"than 5%% of the transaction dates are unparseable, so customer ",
"age and recency cannot be measured — please check that column's ",
"date format."),
format(n_bad_dates, big.mark = ","), format(nrow(df), big.mark = ","),
date_name, 100 * n_bad_dates / nrow(df)))
}
keep_dt <- !is.na(dts)
df <- df[keep_dt, , drop = FALSE]; cid <- cid[keep_dt]; dts <- dts[keep_dt]Step 4: Amount — 95% numeric coercion rule; non-positive amounts
(refunds, zero-value rows) cannot be used by the Gamma-Gamma model.
val <- df$amount
if (!is.numeric(val)) {
conv <- suppressWarnings(as.numeric(as.character(val)))
n_orig <- sum(!is.na(val) & trimws(as.character(val)) != "")
if (n_orig == 0 || sum(!is.na(conv)) < 0.95 * n_orig) {
stop(sprintf(paste0("Column '%s' is not numeric enough to use as a ",
"transaction amount — fewer than 95%% of its values ",
"could be read as numbers."), amount_name))
}
val <- conv
}
keep_val <- !is.na(val)
cid <- cid[keep_val]; dts <- dts[keep_val]; val <- val[keep_val]
n_nonpositive <- sum(val <= 0)
if (length(val) > 0 && n_nonpositive / length(val) > 0.5) {
stop(sprintf(
paste0("%s of %s values in '%s' are zero or negative (%.1f%%). The ",
"Gamma-Gamma monetary model is defined only on positive ",
"transaction values, so with the majority non-positive no ",
"lifetime value can be estimated. Supply gross transaction ",
"amounts rather than a net/refund ledger."),
format(n_nonpositive, big.mark = ","), format(length(val), big.mark = ","),
amount_name, 100 * n_nonpositive / max(1, length(val))))
}
final_rows <- length(val)
rows_removed <- initial_rows - final_rows
tx <- data.frame(customer = cid, txn_date = dts, amount = val,
stringsAsFactors = FALSE)
if (final_rows < 1) {
stop(sprintf("No usable transactions remain after cleaning '%s', '%s', and '%s'.",
cust_name, date_name, amount_name))
}Step 5: Observation window + the BTYD per-customer summary
date_min <- min(tx$txn_date); date_max <- max(tx$txn_date)
span_days <- as.numeric(date_max - date_min)
obs_end <- date_max
cust_df <- btyd_summary(tx, obs_end)
n_same_day_merged <- final_rows - sum(cust_df$n_txn)
n_customers <- nrow(cust_df)
n_repeaters <- sum(cust_df$x > 0)Step 6: Identifiability guards — refuse loudly, never degrade
if (n_customers < 50) {
stop(sprintf(
paste0("Only %d distinct customers in '%s'. Predictive lifetime value ",
"needs at least 50 customers before the purchase-rate and churn ",
"distributions mean anything — with fewer, the fitted parameters ",
"are noise."),
n_customers, cust_name))
}
if (n_repeaters < 20) {
stop(sprintf(
paste0("Only %d of %d customers in '%s' made more than one purchase. ",
"The BG/NBD model learns the churn process from REPEAT ",
"behaviour, so with fewer than 20 repeat customers it is ",
"unidentifiable — no honest lifetime-value projection is ",
"possible from this data. This usually means the export covers ",
"too short a window, or '%s' is not a stable customer identifier."),
n_repeaters, n_customers, cust_name, cust_name))
}
if (span_days < 60) {
stop(sprintf(
paste0("The transactions in '%s' span only %.0f days (%s to %s). A ",
"lifetime-value model needs at least 60 days so that part of the ",
"history can be held back to validate the forecast; without that ",
"back-test the projection would be untestable."),
date_name, span_days, format(date_min), format(date_max)))
}Step 7: Fit BG/NBD on the full observation window
bg <- fit_bgnbd(cust_df$x, cust_df$t_x, cust_df$T_obs)
if (!isTRUE(bg$ok)) {
stop(sprintf(
paste0("The BG/NBD purchase model could not be fitted to this data: %s. ",
"Rather than present lifetime-value figures the model does not ",
"support, the analysis stops here. This usually means the repeat ",
"purchase pattern in '%s' is too sparse or too irregular for the ",
"buy-till-you-die assumptions."),
bg$reason, cust_name))
}Step 8: Fit Gamma-Gamma on repeat customers with positive value
gg_rows <- cust_df$x > 0 & is.finite(cust_df$mean_repeat_value) &
cust_df$mean_repeat_value > 0
if (sum(gg_rows) < 20) {
stop(sprintf(
paste0("Only %d customers have a positive-value repeat purchase in ",
"'%s'. The Gamma-Gamma monetary model needs at least 20 to ",
"estimate transaction value, so no lifetime value can be ",
"reported."),
sum(gg_rows), amount_name))
}
gg <- fit_gg(cust_df$x[gg_rows], cust_df$mean_repeat_value[gg_rows])
if (!isTRUE(gg$ok)) {
stop(sprintf(
paste0("The Gamma-Gamma monetary model could not be fitted to '%s': %s. ",
"Expected transaction value is the money half of lifetime value, ",
"so the analysis stops rather than report a currency figure it ",
"cannot stand behind."),
amount_name, gg$reason))
}Step 9: Independence check — the Gamma-Gamma assumption, measured
fm_cor <- NA_real_; fm_cor_p <- NA_real_
if (sum(gg_rows) >= 3) {
ct <- tryCatch(cor.test(cust_df$x[gg_rows],
cust_df$mean_repeat_value[gg_rows]),
error = function(e) NULL)
if (!is.null(ct)) { fm_cor <- unname(ct$estimate); fm_cor_p <- ct$p.value }
}
independence_violated <- is.finite(fm_cor) && abs(fm_cor) >= 0.1 &&
is.finite(fm_cor_p) && fm_cor_p < 0.05Step 10: Horizon — the data's own span, capped at one year
horizon_days <- max(30, min(365, round(span_days)))Step 11: Per-customer projection
exp_txn <- bgnbd_expected_transactions(bg$params, cust_df$x, cust_df$t_x,
cust_df$T_obs, horizon_days)
if (any(!is.finite(exp_txn)) || any(exp_txn < 0, na.rm = TRUE)) {
stop(sprintf(
paste0("The fitted BG/NBD model produced impossible transaction ",
"forecasts(non-finite or negative) for %d of %d customers, ",
"which means the parameter estimate(r=%.4g, alpha=%.4g, a=%.4g, ",
"b=%.4g) is not usable for projection. No lifetime-value figures ",
"are reported."),
sum(!is.finite(exp_txn) | exp_txn < 0), n_customers,
bg$params[["r"]], bg$params[["alpha"]], bg$params[["a"]], bg$params[["b"]]))
}
p_alive <- bgnbd_palive(bg$params, cust_df$x, cust_df$t_x, cust_df$T_obs)
exp_val <- gg_expected_value(gg$params, cust_df$x, cust_df$mean_repeat_value)
if (any(!is.finite(exp_val))) {
stop("The fitted Gamma-Gamma model produced non-finite expected transaction values; no lifetime-value figures are reported.")
}
cust_df$expected_transactions <- exp_txn
cust_df$p_alive <- p_alive
cust_df$expected_avg_value <- exp_val
cust_df$predicted_clv <- exp_txn * exp_val
total_clv <- sum(cust_df$predicted_clv)
avg_clv <- mean(cust_df$predicted_clv)
median_clv <- median(cust_df$predicted_clv)
total_expected_txn <- sum(cust_df$expected_transactions)
ord <- order(-cust_df$predicted_clv)
n_top10 <- max(1, round(0.10 * n_customers))
top10_share <- 100 * sum(cust_df$predicted_clv[ord][seq_len(n_top10)]) /
max(total_clv, .Machine$double.eps)
pop_mean_value <- gg$params[["gamma"]] * gg$params[["p"]] /
(gg$params[["q"]] - 1)
observed_total <- sum(tx$amount)Step 12: Calibration / holdout back-test — the honesty backbone
Split the observation window at 70% of its span, refit BG/NBD on the calibration half only, and compare each customer's PREDICTED repeat transactions in the holdout window against what they actually did.
hold <- list(available = FALSE, reason = "")
cal_end <- date_min + as.difftime(round(0.7 * span_days), units = "days")
holdout_days <- as.numeric(date_max - cal_end)
cal_tx <- tx[tx$txn_date <= cal_end, , drop = FALSE]
if (holdout_days < 14) {
hold$reason <- sprintf("the holdout window would be only %.0f days long", holdout_days)
} else if (nrow(cal_tx) < 1) {
hold$reason <- "no transactions fall inside the calibration window"
} else {
cal_sum <- btyd_summary(cal_tx, cal_end)
cal_sum <- cal_sum[cal_sum$T_obs > 0, , drop = FALSE]
if (sum(cal_sum$x > 0) < 20) {
hold$reason <- sprintf(
"only %d customers made a repeat purchase inside the calibration window",
sum(cal_sum$x > 0))
} else {
bg_cal <- fit_bgnbd(cal_sum$x, cal_sum$t_x, cal_sum$T_obs)
if (!isTRUE(bg_cal$ok)) {
hold$reason <- paste0("the calibration-period model did not fit(",
bg_cal$reason, ")")
} else {
pred <- bgnbd_expected_transactions(bg_cal$params, cal_sum$x,
cal_sum$t_x, cal_sum$T_obs,
holdout_days)
hold_tx <- tx[tx$txn_date > cal_end, , drop = FALSE]
actual_tbl <- hold_tx %>%
group_by(customer, txn_date) %>%
summarise(.groups = "drop") %>%
group_by(customer) %>%
summarise(actual = n(), .groups = "drop") %>%
as.data.frame(stringsAsFactors = FALSE)
actual <- actual_tbl$actual[match(cal_sum$customer, actual_tbl$customer)]
actual[is.na(actual)] <- 0
ok <- is.finite(pred) & is.finite(actual)
if (sum(ok) < 20) {
hold$reason <- "the calibration model could not produce forecasts for enough customers"
} else {
pr <- pred[ok]; ac <- actual[ok]
ct <- tryCatch(cor.test(pr, ac), error = function(e) NULL)
hcor <- if (!is.null(ct)) unname(ct$estimate) else NA_real_
mae <- mean(abs(pr - ac))
ptot <- sum(pr); atot <- sum(ac)
ratio <- if (atot > 0) ptot / atot else NA_real_
verdict <- if (!is.finite(hcor) || !is.finite(ratio)) "inconclusive"
else if (hcor >= 0.5 && ratio >= 0.8 && ratio <= 1.25) "good"
else if (hcor >= 0.3 && ratio >= 0.6 && ratio <= 1.6) "fair"
else "poor"
set.seed(42)
sidx <- if (length(pr) > 2000) sample(length(pr), 2000) else seq_along(pr)
hold <- list(
available = TRUE, reason = "",
n = length(pr), cor = hcor, mae = mae,
predicted_total = ptot, actual_total = atot, ratio = ratio,
verdict = verdict, holdout_days = holdout_days,
cal_end = cal_end, cal_params = bg_cal$params,
df = data.frame(predicted_transactions = round(pr[sidx], 4),
actual_transactions = as.numeric(ac[sidx]),
stringsAsFactors = FALSE)
)
}
}
}
}