---
title: "How long does an account last? The survival worked example"
subtitle: "Topic-12 worked example (LAT-2369) — B2B account lifetime with censoring"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE, message = FALSE, warning = FALSE)
library(survival)
library(knitr)
```

## The business question

A B2B software company wants one number for its board: **how long does an account last?**
The finance team already has an answer — average the tenure of every account that has
churned. That number is wrong, and it is wrong in a predictable direction.

The reason is **censoring**. On the reporting date a large share of accounts are still
active. Their contracts have not ended, so their lifetimes are not known — but they are
not *missing*. An account that has been live for 900 days and has not churned tells you
something exact: **its lifetime is at least 900 days.** Averaging only the churned
accounts throws every one of those facts away, and every fact it throws away is good
news. The estimate can only come out short.

The second question is the actionable one: **does completing onboarding change how long
an account lasts?** Onboarding costs money, so it has to earn its place.

## The data

One row per account, as it would come out of a subscription billing system.

| column | meaning |
|---|---|
| `account_id` | the account |
| `segment` | SMB / Mid-Market / Enterprise |
| `seats` | licensed seats |
| `mrr_usd` | monthly recurring revenue, US dollars |
| `onboarding` | `completed` / `not_completed` — did the account finish guided onboarding |
| `started_at` | contract start date |
| `churned_at` | churn date, or empty if the account is **still active** |
| `tenure_days` | days from start to churn, or to the reporting date if still active |
| `churned` | 1 = churned (event observed), 0 = **censored** (still active) |

Accounts start on staggered dates, so each has a different length of observation. The
reporting date is fixed at 2026-06-30. Generation is seeded, so re-knitting reproduces
`data.csv` byte-for-byte.

```{r generate}
set.seed(20260824)

n <- 320
reporting_date <- as.Date("2026-06-30")

segment <- sample(c("SMB", "Mid-Market", "Enterprise"), n, replace = TRUE,
                  prob = c(0.55, 0.32, 0.13))
seats <- ifelse(segment == "SMB",         sample(3:25,    n, replace = TRUE),
         ifelse(segment == "Mid-Market",  sample(25:120,  n, replace = TRUE),
                                          sample(120:600, n, replace = TRUE)))
mrr_usd <- round(seats * runif(n, 28, 46), 2)

onboarding <- ifelse(runif(n) < 0.62, "completed", "not_completed")

# staggered entry: every account has a different observation window
follow_up_days <- sample(180:1500, n, replace = TRUE)
started_at     <- reporting_date - follow_up_days

# true lifetimes: Weibull, shorter for accounts that never finished onboarding
shape       <- 1.15
scale_days  <- ifelse(onboarding == "completed", 32 * 30.44, 14 * 30.44)
true_life   <- round(rweibull(n, shape = shape, scale = scale_days))

churned     <- as.integer(true_life <= follow_up_days)
tenure_days <- pmin(true_life, follow_up_days)
churned_at  <- as.Date(ifelse(churned == 1, as.character(started_at + tenure_days), NA))

accounts <- data.frame(
  account_id  = sprintf("ACC-%04d", seq_len(n)),
  segment     = segment,
  seats       = seats,
  mrr_usd     = mrr_usd,
  onboarding  = onboarding,
  started_at  = started_at,
  churned_at  = churned_at,
  tenure_days = tenure_days,
  churned     = churned,
  stringsAsFactors = FALSE
)
accounts <- accounts[order(accounts$account_id), ]
rownames(accounts) <- NULL
write.csv(accounts, "data.csv", row.names = FALSE)

n_events   <- sum(accounts$churned)
n_censored <- sum(accounts$churned == 0)
data.frame(accounts = n, churned = n_events, still_active = n_censored,
           censored_pct = round(100 * n_censored / n, 2))
```

```{r peek}
kable(head(accounts[, c("account_id","segment","seats","mrr_usd","onboarding",
                        "started_at","churned_at","tenure_days","churned")], 8),
      caption = "First eight rows as the billing export delivers them")
```

## The naive answer, and how far off it is

```{r naive}
churned_only <- accounts$tenure_days[accounts$churned == 1]
naive_mean   <- mean(churned_only)
naive_median <- median(churned_only)
c(naive_mean_days = round(naive_mean, 2), naive_median_days = round(naive_median, 2),
  naive_mean_months = round(naive_mean / 30.44, 2))
```

This is the finance team's number: it answers *"how long did the accounts that already
left last?"* — a question about the dead, not about the book of business.

## Kaplan-Meier: count every account for as long as it was actually observed

```{r km-overall}
fit_all <- survfit(Surv(tenure_days, churned) ~ 1, data = accounts)
km_med  <- summary(fit_all)$table
km_med[c("records", "events", "median", "0.95LCL", "0.95UCL")]
```

```{r km-anchors}
anchors <- summary(fit_all, times = c(365, 730, 1095))
data.frame(day = anchors$time,
           n_at_risk = anchors$n.risk,
           survival = round(anchors$surv, 4),
           lower = round(anchors$lower, 4),
           upper = round(anchors$upper, 4))
```

```{r gap}
km_median_days <- unname(km_med["median"])
gap_pct <- 100 * (km_median_days - naive_median) / km_median_days
c(naive_median_days = naive_median, km_median_days = km_median_days,
  naive_understates_pct = round(gap_pct, 2))
```

## Restricted mean survival over a fixed two-year horizon

The median answers "half the book outlasts this." If you need an *average* instead, the
honest one is the **restricted mean** — the area under the survival curve out to a stated
horizon. The horizon is part of the number; quoting it without the horizon is meaningless.

```{r rmst}
rmst <- function(fit, horizon) {
  t <- c(0, fit$time); s <- c(1, fit$surv)
  keep <- t <= horizon
  t <- c(t[keep], horizon); s <- c(s[keep], tail(s[keep], 1))
  sum(diff(t) * head(s, -1))          # left-continuous step integral
}
rmst_all <- rmst(fit_all, 730)
c(rmst_730_days = round(rmst_all, 2), rmst_730_months = round(rmst_all / 30.44, 2))
```

## Does onboarding change account lifetime?

```{r km-group}
fit_grp <- survfit(Surv(tenure_days, churned) ~ onboarding, data = accounts)
grp_tab <- summary(fit_grp)$table
kable(round(grp_tab[, c("records","events","median","0.95LCL","0.95UCL")], 2),
      caption = "Median account lifetime in days, by onboarding status")
```

```{r logrank}
lr <- survdiff(Surv(tenure_days, churned) ~ onboarding, data = accounts)
lr
lr_p <- 1 - pchisq(lr$chisq, df = length(lr$n) - 1)
signif(lr_p, 6)
```

```{r cox}
cx <- coxph(Surv(tenure_days, churned) ~ onboarding, data = accounts)
summary(cx)
hr    <- unname(exp(coef(cx)))
hr_ci <- unname(exp(confint(cx)))
c(hazard_ratio_notcompleted_vs_completed = round(hr, 4),
  ci_lo = round(hr_ci[1], 4), ci_hi = round(hr_ci[2], 4))
```

```{r rmst-group}
fits <- survfit(Surv(tenure_days, churned) ~ onboarding, data = accounts)
strata_names <- names(fits$strata)
idx <- rep(seq_along(fits$strata), fits$strata)
rmst_by <- sapply(seq_along(strata_names), function(k) {
  sub <- list(time = fits$time[idx == k], surv = fits$surv[idx == k])
  rmst(sub, 730)
})
names(rmst_by) <- strata_names
round(rmst_by, 2)
```

## The charts a practitioner reads

```{r km-plot, fig.width=9, fig.height=5}
plot(fit_grp, col = c("#F97316", "#5fa9dd"), lwd = 2.4, mark.time = TRUE,
     xlab = "days since contract start", ylab = "share of accounts still active",
     main = "Kaplan-Meier: account survival by onboarding status",
     xlim = c(0, 1200), yaxs = "i", ylim = c(0, 1.02))
abline(h = 0.5, lty = 3, col = "#888888")
legend("topright", legend = c("onboarding completed", "onboarding not completed"),
       col = c("#F97316", "#5fa9dd"), lwd = 2.4, bty = "n")
mtext("vertical ticks = censored accounts (still active at the reporting date)",
      side = 3, line = 0.2, cex = 0.8, col = "#666666")
```

```{r gap-plot, fig.width=7, fig.height=4}
vals <- c(`churned accounts only\n(the naive answer)` = naive_median,
          `Kaplan-Meier median\n(every account counted)` = km_median_days)
bp <- barplot(vals, col = c("#5fa9dd", "#F97316"), ylim = c(0, max(vals) * 1.25),
              ylab = "median account lifetime (days)",
              main = "Dropping the still-active accounts shortens the answer")
text(bp, vals + max(vals) * 0.05, sprintf("%.0f days", vals), font = 2)
```

## What the lesson teaches

1. **Censored is not missing.** An account still active at the reporting date carries a
   hard fact — *at least this long*. Methods that only count finished lifetimes answer a
   question about the accounts you already lost.
2. **Report the median off the curve**, not the average of the churned. If you need an
   average, use the restricted mean **and say the horizon**.
3. **Compare groups with the log-rank test**, which uses every account and every censoring
   time, not a single summary point.
4. **State the effect as a hazard ratio with its interval.** An interval that clears 1 is
   what makes the difference reportable.
5. **The tool needs row-level data**: one row per account, a duration, an event flag, and
   a group. Totals cannot be un-summed back into a survival curve.