Chapter 3 Project 1: Customer Churn Prediction & Retention Automation

Customer churn — the loss of customers over time — is one of the most expensive problems in business. Acquiring a new customer costs 5-25x more than retaining an existing one. In this chapter, we build a complete churn prediction system: from database-driven feature engineering to a real-time scoring API that triggers automated retention workflows.

3.1 Business Context

Our goal is to predict which customers will churn within the next 30 days and automatically trigger retention interventions. The system must:

  • Refresh features daily from production databases
  • Handle severe class imbalance (typically <10% churn rate)
  • Provide explainable predictions for customer success teams
  • Integrate with retention workflows via API

3.2 Architecture Overview

PostgreSQL ──► Feature Engineering ──► Model Training ──► Plumber API
    │                │                      │                  │
    │                ▼                      ▼                  ▼
    │         dbplyr lazy eval      XGBoost + SMOTE     Retention trigger
    │                                                           │
    └───────────────────────────────────────────────────────────┘
                              Feedback loop

3.3 Data Engineering with dbplyr

We start by building a feature-rich customer snapshot directly from the database. Using dbplyr, we push computation to PostgreSQL and only bring aggregated results into R.

library(tidyverse)
library(tidymodels)
library(DBI)
library(RPostgres)
library(dbplyr)

# Connect to PostgreSQL (see Chapter 2 for connection pooling)
con <- dbConnect(
  RPostgres::Postgres(),
  dbname = "analytics",
  host = "prod-db.company.com",
  port = 5432,
  user = Sys.getenv("DB_USER"),
  password = Sys.getenv("DB_PASS")
)

# Build feature-rich customer snapshot using dbplyr (lazy eval on DB)
customer_features <- tbl(con, "customers") %>%
  left_join(
    tbl(con, "user_activity") %>%
      filter(session_date >= Sys.Date() - 30) %>%
      group_by(customer_id) %>%
      summarise(
        active_days_30d = n_distinct(session_date),
        total_events_30d = sum(events, na.rm = TRUE),
        avg_session_duration = mean(session_duration_sec, na.rm = TRUE),
        last_active_date = max(session_date, na.rm = TRUE),
        .groups = "drop"
      ),
    by = "customer_id"
  ) %>%
  left_join(
    tbl(con, "support_tickets") %>%
      filter(created_at >= Sys.Date() - 30) %>%
      group_by(customer_id) %>%
      summarise(
        ticket_count_30d = n(),
        avg_satisfaction = mean(satisfaction_rating, na.rm = TRUE),
        critical_tickets = sum(priority == "critical", na.rm = TRUE),
        .groups = "drop"
      ),
    by = "customer_id"
  ) %>%
  mutate(
    days_since_last_active = as.numeric(Sys.Date() - last_active_date),
    avg_satisfaction = coalesce(avg_satisfaction, 3.0),
    ticket_count_30d = coalesce(ticket_count_30d, 0L)
  ) %>%
  collect()

Note: The coalesce() function replaces NA values with defaults. For satisfaction, we assume neutral (3.0) when no tickets exist. For ticket counts, we assume zero.

3.4 Feature Engineering

Feature engineering transforms raw data into signals that models can learn from. We design features that capture engagement decay, behavioral patterns, and RFM-style scoring.

build_churn_features <- function(df) {
  df %>%
    mutate(
      # Engagement velocity: recent vs. expected activity
      # If events_7d is much lower than the 30-day average, engagement is decaying
      engagement_velocity = (events_7d - events_30d / 4.29) /
        (events_30d / 4.29 + 1e-6),

      # Normalized recency: how long since last login relative to tenure
      days_since_login_ratio = days_since_last_active / tenure_days,

      # Behavioral flags
      power_user_flag = as.integer(active_days_30d >= 20),
      at_risk_support = as.integer(critical_tickets >= 2),

      # RFM-style score (Recency, Frequency, Monetary)
      r_score = ntile(desc(days_since_last_active), 5),
      f_score = ntile(total_events_30d, 5),
      m_score = ntile(mrr, 5),
      rfm_score = r_score * 100 + f_score * 10 + m_score,

      # Target variable
      churned = as.factor(churned_within_30d)
    )
}

df <- customer_features %>% build_churn_features()

Tip: The 1e-6 epsilon prevents division by zero when a customer has no 30-day activity. This is a common defensive programming pattern.

3.5 Train/Test Split with Stratification

Stratified sampling ensures both training and test sets maintain the same churn rate as the full dataset. This is critical for imbalanced classification.

set.seed(42)
split <- initial_split(df, prop = 0.8, strata = churned)
train <- training(split)
test <- testing(split)

# 5-fold stratified cross-validation
cv_folds <- vfold_cv(train, v = 5, strata = churned)

3.6 Recipe: The Feature Pipeline

recipes from tidymodels defines a reproducible preprocessing pipeline. We handle missing values, normalize features, and apply SMOTE for class imbalance.

churn_recipe <- recipe(churned ~ ., data = train) %>%
  # Remove identifiers and intermediate scores
  step_rm(customer_id, last_active_date, r_score, f_score, m_score) %>%

  # Impute missing values
  step_impute_median(all_numeric_predictors()) %>%
  step_impute_mode(all_nominal_predictors()) %>%

  # Normalize numeric features (zero mean, unit variance)
  step_normalize(all_numeric_predictors()) %>%

  # SMOTE for class imbalance: oversample minority to 50% of majority
  step_smote(churned, over_ratio = 0.5) %>%

  # Remove zero-variance predictors
  step_zv(all_predictors())

Warning: SMOTE should only be applied to the training data, never to validation or test sets. The recipe handles this correctly when used within a workflow and resampling.

3.7 XGBoost Model with Hyperparameter Tuning

XGBoost is our algorithm of choice for churn prediction due to its ability to handle mixed data types, missing values, and non-linear relationships.

library(xgboost)
library(tune)
library(vip)

xgb_spec <- boost_tree(
  trees = tune(),        # Number of trees
  tree_depth = tune(),   # Max depth per tree
  min_n = tune(),        # Minimum observations in terminal nodes
  loss_reduction = tune(), # Minimum loss reduction for split
  sample_size = tune(),  # Row subsampling proportion
  mtry = tune(),         # Column subsampling proportion
  learn_rate = tune()    # Learning rate (shrinkage)
) %>%
  set_engine("xgboost", scale_pos_weight = 7, eval_metric = "aucpr") %>%
  set_mode("classification")

xgb_workflow <- workflow() %>%
  add_recipe(churn_recipe) %>%
  add_model(xgb_spec)

The scale_pos_weight = 7 parameter tells XGBoost that positive cases (churners) are 7x rarer than negative cases, adjusting the learning objective accordingly. We optimize for aucpr (Area Under the Precision-Recall Curve) because accuracy is misleading with imbalanced data.

3.7.2 Parallel Tuning

library(doParallel)
cl <- makePSOCKcluster(parallel::detectCores() - 1)
registerDoParallel(cl)

xgb_results <- tune_grid(
  xgb_workflow,
  resamples = cv_folds,
  grid = xgb_grid,
  metrics = metric_set(roc_auc, pr_auc, f_meas),
  control = control_grid(save_pred = TRUE, parallel_over = "everything")
)

stopCluster(cl)

# Select best hyperparameters by PR-AUC
best_params <- select_best(xgb_results, metric = "pr_auc")
final_xgb <- finalize_workflow(xgb_workflow, best_params)
final_fit <- fit(final_xgb, data = train)

3.8 Model Evaluation

For churn prediction, we care about precision (of flagged customers, how many actually churn?) and recall (of all churners, how many did we catch?).

test_pred <- predict(final_fit, test, type = "prob") %>%
  bind_cols(predict(final_fit, test)) %>%
  bind_cols(test)

# ROC-AUC: Discrimination ability
test_pred %>% roc_auc(churned, .pred_1)

# Precision: Quality of positive predictions
test_pred %>% precision(churned, .pred_class)

# Recall: Coverage of actual churners
test_pred %>% recall(churned, .pred_class)

# Confusion matrix
test_pred %>% conf_mat(churned, .pred_class)

Interpreting Results: In churn prediction, false negatives (missing a churner) are typically more costly than false positives (unnecessary retention spend). Optimize recall if retention is cheap; optimize precision if retention is expensive.

3.9 Explainability with DALEX

Black-box models require explainability for production deployment. We use DALEX to understand global feature importance and individual prediction drivers.

library(DALEX)

explainer <- explain(
  final_fit,
  data = test %>% select(-churned),
  y = as.numeric(test$churned) - 1,
  label = "XGBoost Churn",
  verbose = FALSE
)

# Global variable importance
vi <- model_parts(explainer, N = 1000)
plot(vi)

# Individual prediction breakdown (SHAP-like)
bd <- predict_parts(explainer, new_observation = test[1, ], type = "break_down")
plot(bd)

3.10 Production API with Plumber

We expose the model as a REST API using plumber. The API fetches fresh features, scores the customer, and triggers retention workflows for high-risk customers.

library(plumber)

#* @apiTitle Churn Prediction API
#* @apiDescription Real-time churn scoring for customer retention

#* Score a single customer
#* @post /predict
#* @param customer_id:string Customer identifier
#* @serializer json
function(req, res, customer_id = NULL) {

  # Fetch fresh features from DB
  features <- tbl(con, "customers") %>%
    filter(customer_id == !!customer_id) %>%
    left_join(
      tbl(con, "user_activity") %>%
        filter(session_date >= Sys.Date() - 30) %>%
        group_by(customer_id) %>%
        summarise(
          active_days_30d = n_distinct(session_date),
          total_events_30d = sum(events, na.rm = TRUE),
          avg_session_duration = mean(session_duration_sec, na.rm = TRUE),
          last_active_date = max(session_date, na.rm = TRUE),
          .groups = "drop"
        ),
      by = "customer_id"
    ) %>%
    collect() %>%
    build_churn_features()

  # Predict
  pred <- predict(final_fit, features, type = "prob")
  class_pred <- predict(final_fit, features)

  # Trigger retention workflow if high risk
  if (pred$.pred_1 > 0.7) {
    trigger_retention_workflow(customer_id, pred$.pred_1)
  }

  list(
    customer_id = customer_id,
    churn_probability = round(pred$.pred_1, 4),
    risk_tier = case_when(
      pred$.pred_1 > 0.8 ~ "critical",
      pred$.pred_1 > 0.5 ~ "high",
      pred$.pred_1 > 0.2 ~ "medium",
      TRUE ~ "low"
    ),
    predicted_class = as.character(class_pred$.pred_class),
    top_drivers = get_top_shap_drivers(features, n = 3)
  )
}

Run the API:

pr("churn_api.R") %>% pr_run(port = 8000)

3.11 Monitoring & Retraining

Production models degrade over time. Implement:

  1. Performance monitoring: Track PR-AUC on a holdout set weekly
  2. Data drift detection: Monitor feature distributions with KS tests
  3. Automated retraining: Trigger retraining when PR-AUC drops below 0.65
  4. A/B testing: Compare new models against the production model

3.12 Exercises

  1. Implement the get_top_shap_drivers() function that returns the top 3 features driving a prediction.
  2. Add a batch scoring endpoint (/predict_batch) that accepts multiple customer IDs.
  3. Implement model versioning: save each trained model with a timestamp and load the latest by default.
  4. Write a data drift detection function that compares training vs. current feature distributions.

In the next chapter, we’ll build an interactive sales performance dashboard that gives executives real-time visibility into pipeline health.