Executive Summary
Time-to-event verdict across 600 subjects.
Median churn time is 11.16 months for Plan A and 9.25 months for Plan B. The log-rank test confirms these curves differ significantly (p = 9.72e-05). Plan B customers churn 1.5× faster than Plan A customers, a substantial and statistically robust difference. Across both plans, 30.3% of the 600 customers remained active at last observation, meaning the median churn time of 10.16 months overall is a reliable estimate unbiased by censoring.
Analysis Overview
Kaplan-Meier time-to-event analysis of 600 subjects across 2 groups.
This analysis tracks 600 customers across two subscription plans to determine when they churn. Of these, 418 customers churned and 182 (30.3%) were still active at last observation. Censoring—the presence of active customers whose churn time remains unknown—requires specialized statistical handling. Naive averages would artificially shorten duration estimates by ignoring still-active accounts. Kaplan-Meier estimation counts censored customers for exactly as long as they remained observed, producing honest survival curves. Log-rank testing and Cox regression then compare churn risk between Plan A and Plan B, expressing relative hazard as a ratio.
Data Quality
Row exclusions, event-flag interpretation, and censoring rate.
All 600 rows loaded were usable; no records were excluded. Churn was defined as Churned = "1"; all other values were treated as censored (customer still active). The censoring rate of 30.3% means 182 customers had not churned when observation ended. This proportion is correctly incorporated into Kaplan-Meier estimates rather than discarded, ensuring that durations reflect the actual time each customer was observed and at risk of churning. Data quality was clean with no missing or invalid duration values, allowing the full cohort to contribute to survival estimates.
Survival Curves
Percentage still event-free over Months Active, by Customer Plan.
The Kaplan-Meier curves show the percentage of customers still active (not churned) over time. Plan A maintains higher retention: at month 0, both plans start at 100% active. Plan B's curve drops faster, crossing the 50% survival line at 9.25 months compared to Plan A's 11.16 months. The curves separate early and remain apart throughout the observation period, indicating a persistent and genuine difference in churn risk between the two plans. Censored customers (still active) leave the curve without pulling it downward, preserving the integrity of survival estimates and accurately reflecting the underlying churn dynamics.
Median Time to Event
Per-Customer Plan medians with 95% confidence intervals.
| Group | N | Events | Censored PCT | Median Time | CI Low | CI High |
|---|---|---|---|---|---|---|
| Plan A | 300 | 208 | 30.7 | 11.16 | 9.68 | 15.21 |
| Plan B | 300 | 210 | 30 | 9.25 | 8.51 | 10.78 |
Plan A (n=300) recorded 208 churn events with a median churn time of 11.16 months (95% CI: 9.68–15.21). Plan B (n=300) recorded 210 churn events with a median churn time of 9.25 months (95% CI: 8.51–10.78). Both plans had similar censoring rates (Plan A: 30.7%; Plan B: 30%). The median—the point at which half the customers have churned—is more robust than the mean for churn data and directly answers "how long until half the customers leave?" Plan B's narrower confidence interval (9.25 ± ~1.14 months) reflects tighter clustering of churn times compared to Plan A's wider interval (11.16 ± ~2.77 months).
Statistical Tests
Log-rank test and Cox hazard ratios with confidence intervals.
| Test | Statistic | P Value | Interpretation |
|---|---|---|---|
| Log-rank test (do the curves differ?) | 15.19 | 0.0001 | The survival curves differ significantly across Customer Plan groups (p = 9.72e-05). |
| Cox hazard ratio: Plan B vs Plan A | 1.48 | 0.0001 | Plan B reaches the event 1.5× faster than Plan A (95% CI 1.21–1.80). |
The log-rank test statistic of 15.19 (p = 9.72e-05) provides strong evidence that Plan A and Plan B have genuinely different churn curves. The Cox hazard ratio for Plan B versus Plan A is 1.48 (95% CI: 1.21–1.80, p = 0.0001), meaning Plan B customers churn 1.5× faster at any given moment. This hazard ratio is both statistically significant and clinically meaningful: a 48% elevation in churn risk is substantial enough to warrant strategic attention. The tight confidence interval (1.21–1.80) indicates precise estimation, ruling out chance variation as an explanation for the observed difference.
When Events Happen
Distribution of observed event times (censored subjects excluded).
Among the 418 observed churn events, the distribution is front-loaded: the first quartile occurs at 3.05 months, the median at 7.42 months, and the third quartile at 14 months. This skew toward early churn means most customers who churn do so within the first few months of the observation window. The data range spans from 0.11 months to 61.34 months, showing that while most churn is concentrated early, a tail of long-surviving customers extends substantially. This early-risk pattern is typical of subscription models and underscores the importance of retention efforts in the initial months after signup.
Methodology
Statistical methodology and diagnostics for Survival Analysis — Time to Event
Statistical Method
Standard-library analysis: how long until customers churn, machines fail, or patients relapse — with censoring handled properly. Kaplan-Meier survival curves per group, median time-to-event with confidence intervals, a log-rank test of whether the groups truly differ, and Cox hazard ratios that say how much faster one group reaches the event than another. Works on any duration data: map a time column, an event flag, and a group.
- Durations are positive numbers on a consistent scale (days, months, cycles)
- Censoring is non-informative — subjects still active are not systematically different
- Cox hazard ratios assume proportional hazards (a roughly constant relative risk over time)
- Groups with fewer than 5 subjects are excluded; group levels beyond the 6 largest are lumped into Other
- Median survival is 'not reached' when more than half a group is still event-free — the data cannot say more
- The log-rank test says the curves differ, not why — confounding is not addressed
- Rows with missing, zero, or negative durations are dropped and reported
Analysis Code
Complete R source code for this analysis
Survival Analysis — Time to Event
Estimates how long subjects last until an event (churn, failure, relapse) with censoring handled properly: Kaplan-Meier curves per group, median time-to-event with confidence intervals, a log-rank test of whether the groups differ, and Cox hazard ratios versus the largest group.
Why This Method?
Duration data almost always contains open cases — customers still active, machines still running. Naive averages of observed durations treat those as finished and systematically mislead. Kaplan-Meier estimation, the log-rank test, and Cox regression are the standard toolkit that counts censored subjects correctly.
What This Analysis Covers
- Kaplan-Meier survival curves per group
- Median time-to-event per group with 95% confidence intervals
- Log-rank test of whether the curves differ
- Cox hazard ratios versus the reference (largest) group
- Distribution of observed event times
Standard Library
Platform standard-library module (LAT-1441): runs on ANY dataset via the semantic mapping {time, event, group}. All narrative is derived from the user's own column names and computed values.
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))
suppressPackageStartupMessages(library(survival))Core Analysis Pipeline
Step 1: Discover mapped columns (humanized names for ALL prose)
initial_rows <- nrow(df)
for (k in c("time", "event", "group")) {
if (!k %in% names(df)) {
stop(sprintf("column_mapping must map the '%s' column.", k))
}
}
h_time <- humanize_semantic("time", col_map)
h_event <- humanize_semantic("event", col_map)
h_group <- humanize_semantic("group", col_map)Step 2: Coerce time numeric (95% rule); drop invalid (<=0 / NA) durations
tv <- df$time
if (!is.numeric(tv)) {
conv <- suppressWarnings(as.numeric(as.character(tv)))
n_orig <- sum(!is.na(tv) & as.character(tv) != "")
if (n_orig > 0 && sum(!is.na(conv)) >= 0.95 * n_orig) {
tv <- conv
} else {
stop(sprintf("The duration column(%s) could not be read as numeric.", h_time))
}
}
bad_time <- is.na(tv) | tv <= 0
dropped_time_rows <- sum(bad_time)Step 3: Binarize the event flag (robust two-state logic)
eb <- .binarize_event(df$event, h_event)
ev <- eb$ev
event_level <- eb$level
bad_event <- is.na(ev)
dropped_event_rows <- sum(bad_event & !bad_time)
keep <- !bad_time & !bad_event
d <- data.frame(
time = tv[keep],
event = as.integer(ev[keep]),
group = trimws(as.character(df$group[keep])),
stringsAsFactors = FALSE
)
d$group[is.na(d$group) | d$group == ""] <- "Missing"Step 4: Group hygiene — lump to top 6 levels, drop n<5 groups
tab <- sort(table(d$group), decreasing = TRUE)
lumped_levels <- character(0)
if (length(tab) > 6) {
lumped_levels <- names(tab)[-(1:6)]
d$group[d$group %in% lumped_levels] <- "Other"
tab <- sort(table(d$group), decreasing = TRUE)
}
small <- names(tab)[tab < 5]
dropped_small_groups <- small
if (length(small) > 0) {
d <- d[!d$group %in% small, , drop = FALSE]
tab <- sort(table(d$group), decreasing = TRUE)
}
if (nrow(d) == 0 || length(tab) == 0) {
stop(sprintf("No usable groups remained in %s after cleaning(each group needs at least 5 rows).", h_group))
}
final_rows <- nrow(d)
rows_removed <- initial_rows - final_rows
if (final_rows < 30) {
stop(sprintf("Only %d usable rows across %s, %s and %s — survival analysis needs at least 30.",
final_rows, h_time, h_event, h_group))
}
n_events <- sum(d$event)
n_censored <- final_rows - n_events
censored_pct <- 100 * n_censored / final_rows
if (n_events == 0) {
stop(sprintf("No events observed in %s(every row is censored) — cannot estimate time-to-event.", h_event))
}
group_levels <- names(tab) # largest first
single_group <- length(group_levels) < 2
ref_group <- group_levels[1]
d$group <- factor(d$group, levels = group_levels)Step 5: Kaplan-Meier fit (per group, or overall in single-group mode)
fit <- if (single_group) {
survfit(Surv(time, event) ~ 1, data = d, conf.int = 0.95)
} else {
survfit(Surv(time, event) ~ group, data = d, conf.int = 0.95)
}
s <- summary(fit)
strata_lab <- if (is.null(s$strata)) {
rep(if (single_group) group_levels[1] else "All", length(s$time))
} else {
sub("^group=", "", as.character(s$strata))
}KM curve dataset: downsample to <=100 points per group, prepend (0, 100%)
km_parts <- list()
for (g in unique(strata_lab)) {
idx <- which(strata_lab == g)
if (length(idx) > 100) {
idx <- idx[unique(round(seq(1, length(idx), length.out = 100)))]
}
km_parts[[g]] <- data.frame(
time_point = c(0, s$time[idx]),
survival_pct = round(100 * c(1, s$surv[idx]), 2),
group_label = g,
ci_low = round(100 * c(1, s$lower[idx]), 2),
ci_high = round(100 * c(1, s$upper[idx]), 2),
stringsAsFactors = FALSE
)
}
km_df <- do.call(rbind, km_parts)
rownames(km_df) <- NULLStep 6: Median survival per group (NA-safe: "not reached")
st <- summary(fit)$table
if (is.null(dim(st))) st <- matrix(st, nrow = 1, dimnames = list(group_levels[1], names(st)))
med_groups <- sub("^group=", "", rownames(st))
median_num <- setNames(as.numeric(st[, "median"]), med_groups)
median_df <- data.frame(
group = med_groups,
n = as.integer(st[, "records"]),
events = as.integer(st[, "events"]),
censored_pct = round(100 * (st[, "records"] - st[, "events"]) / st[, "records"], 1),
median_time = .fmt_time(st[, "median"]),
ci_low = .fmt_time(st[, "0.95LCL"], na_label = "—"),
ci_high = .fmt_time(st[, "0.95UCL"], na_label = "—"),
stringsAsFactors = FALSE
)
rownames(median_df) <- NULL
fit_all <- survfit(Surv(time, event) ~ 1, data = d)
overall_median <- as.numeric(summary(fit_all)$table["median"])Step 7: Log-rank test + Cox hazard ratios vs the largest group
logrank_p <- NA_real_; logrank_chisq <- NA_real_
hr_df <- NULL
test_rows <- list()
if (!single_group) {
sd <- tryCatch(survdiff(Surv(time, event) ~ group, data = d), error = function(e) NULL)
if (!is.null(sd)) {
logrank_chisq <- as.numeric(sd$chisq)
logrank_p <- pchisq(logrank_chisq, df = length(sd$n) - 1, lower.tail = FALSE)
}
lr_interp <- if (is.na(logrank_p)) {
sprintf("The log-rank test could not be computed for %s.", h_group)
} else if (logrank_p < 0.05) {
sprintf("The survival curves differ significantly across %s groups(p = %s).",
h_group, signif(logrank_p, 3))
} else {
sprintf("No significant difference between the %s survival curves(p = %s).",
h_group, signif(logrank_p, 3))
}
test_rows[[1]] <- data.frame(
test = "Log-rank test(do the curves differ?)",
statistic = ifelse(is.na(logrank_chisq), NA_real_, round(logrank_chisq, 2)),
p_value = ifelse(is.na(logrank_p), NA_real_, signif(logrank_p, 3)),
interpretation = lr_interp,
stringsAsFactors = FALSE
)
cx <- tryCatch(
suppressWarnings(coxph(Surv(time, event) ~ group, data = d)),
error = function(e) NULL
)
if (!is.null(cx)) {
sm <- summary(cx)
ci <- sm$conf.int; cf <- sm$coefficients
if (is.null(dim(ci))) ci <- matrix(ci, nrow = 1, dimnames = list(rownames(sm$conf.int), colnames(sm$conf.int)))
g_names <- sub("^group", "", rownames(ci))
hr_df <- data.frame(
group = g_names,
hr = as.numeric(ci[, "exp(coef)"]),
ci_low = as.numeric(ci[, "lower .95"]),
ci_high = as.numeric(ci[, "upper .95"]),
p = as.numeric(cf[, "Pr(>|z|)"]),
stringsAsFactors = FALSE
)
for (i in seq_len(nrow(hr_df))) {
hr <- hr_df$hr[i]; g <- hr_df$group[i]
interp <- if (is.na(hr)) {
sprintf("Hazard ratio for %s could not be estimated.", g)
} else if (hr < 0.01 || hr > 100) {
sprintf("Risk for %s could not be reliably estimated(no or nearly no events in the group).", g)
} else if (hr >= 1) {
sprintf("%s reaches the event %.1f× faster than %s(95%% CI %.2f–%.2f).",
g, hr, ref_group, hr_df$ci_low[i], hr_df$ci_high[i])
} else {
sprintf("%s reaches the event %.1f× slower than %s(HR %.2f, 95%% CI %.2f–%.2f).",
g, 1 / hr, ref_group, hr, hr_df$ci_low[i], hr_df$ci_high[i])
}
test_rows[[length(test_rows) + 1]] <- data.frame(
test = sprintf("Cox hazard ratio: %s vs %s", g, ref_group),
statistic = ifelse(is.na(hr), NA_real_, round(hr, 2)),
p_value = ifelse(is.na(hr_df$p[i]), NA_real_, signif(hr_df$p[i], 3)),
interpretation = interp,
stringsAsFactors = FALSE
)
}
} else {
test_rows[[length(test_rows) + 1]] <- data.frame(
test = "Cox proportional hazards",
statistic = NA_real_, p_value = NA_real_,
interpretation = sprintf("The Cox model could not be fitted for %s.", h_group),
stringsAsFactors = FALSE
)
}
} else {
test_rows[[1]] <- data.frame(
test = "Log-rank test",
statistic = NA_real_, p_value = NA_real_,
interpretation = sprintf(
"Only one %s group was found — the analysis reports a single overall survival curve; no between-group test applies.",
h_group),
stringsAsFactors = FALSE
)
}
tests_df <- do.call(rbind, test_rows)
rownames(tests_df) <- NULLBiggest hazard ratio in plain English (NA-safe: never which.max over all-NA)
top_hr_text <- ""
if (!is.null(hr_df) && nrow(hr_df) > 0) {
mag <- pmax(hr_df$hr, 1 / hr_df$hr)
ok <- which(!is.na(mag) & hr_df$hr > 0.01 & hr_df$hr < 100)
if (length(ok) > 0) {
i <- ok[which.max(mag[ok])]
top_hr_text <- if (hr_df$hr[i] >= 1) {
sprintf("%s reaches the event %.1f× faster than %s", hr_df$group[i], hr_df$hr[i], ref_group)
} else {
sprintf("%s reaches the event %.1f× slower than %s", hr_df$group[i], 1 / hr_df$hr[i], ref_group)
}
}
}Step 8: Event-time distribution (events only, <=2000 rows)
et <- d$time[d$event == 1]
set.seed(42)
if (length(et) > 2000) et <- sample(et, 2000)
event_times_df <- data.frame(event_time = round(et, 4), stringsAsFactors = FALSE)
metrics <- list(
`Subjects` = final_rows,
`Events` = as.integer(n_events),
`Censored %` = round(censored_pct, 1),
`Groups` = length(group_levels),
`Median Time` = .fmt_time(overall_median),
`Log-rank p` = if (is.na(logrank_p)) "—" else signif(logrank_p, 3)
)
med_phrase <- if (single_group) {
sprintf("median %s to event %s", h_time, .fmt_time(overall_median))
} else {
paste(sapply(seq_len(nrow(median_df)), function(i) {
sprintf("%s: %s", median_df$group[i], median_df$median_time[i])
}), collapse = "; ")
}
lr_phrase <- if (single_group) {
"single overall curve(one group)"
} else if (is.na(logrank_p)) {
"log-rank unavailable"
} else if (logrank_p < 0.05) {
sprintf("curves differ significantly(log-rank p = %s)", signif(logrank_p, 3))
} else {
sprintf("no significant curve difference(log-rank p = %s)", signif(logrank_p, 3))
}
json_output <- list(
answer = paste0(
"Kaplan-Meier survival analysis of ", format(final_rows, big.mark = ","),
" subjects(", n_events, " events, ", round(censored_pct, 1), "% censored)",
if (single_group) "" else paste0(" across ", length(group_levels), " ", h_group, " groups"),
". Median time to event — ", med_phrase, ". ", lr_phrase,
if (nchar(top_hr_text) > 0) paste0("; ", top_hr_text, ".") else "."
),
cards = lapply(
c("tldr", "overview", "preprocessing", "km_curves",
"median_survival", "test_results", "event_distribution"),
function(cid) list(id = cid, metrics = metrics)
)
)
list(
initial_rows = initial_rows, final_rows = final_rows,
rows_removed = rows_removed,
h_time = h_time, h_event = h_event, h_group = h_group,
event_level = event_level,
n_events = n_events, n_censored = n_censored, censored_pct = censored_pct,
dropped_time_rows = dropped_time_rows, dropped_event_rows = dropped_event_rows,
dropped_small_groups = dropped_small_groups, lumped_levels = lumped_levels,
single_group = single_group, group_levels = group_levels, ref_group = ref_group,
km_df = km_df, median_df = median_df, median_num = median_num,
tests_df = tests_df, logrank_p = logrank_p, logrank_chisq = logrank_chisq,
hr_df = hr_df, top_hr_text = top_hr_text,
event_times_df = event_times_df, overall_median = overall_median,
metrics = metrics, json_output = json_output
)
}Fastest-declining group: lowest median, NA-safe ("not reached" = best survival)
fastest_note <- ""
if (!shared$single_group) {
mn <- shared$median_num
ok <- which(!is.na(mn))
if (length(ok) > 0) {
worst <- names(mn)[ok][which.min(mn[ok])]
fastest_note <- sprintf(
" %s declines fastest(median %s %s), so it reaches the event soonest.",
worst, shared$h_time, .fmt_time(min(mn[ok])))
}
}
list(
title = "Survival Curves",
description = paste0("Percentage still event-free over ", shared$h_time,
if (shared$single_group) "" else paste0(", by ", shared$h_group), "."),
text = paste0(
"Each curve tracks the percentage of subjects still event-free as ",
shared$h_time, " increases, starting at 100% and stepping down at each ",
"observed event. Where a curve crosses the 50% line is that group's ",
"median time to event.", fastest_note,
if (shared$single_group) "" else
" Curves that separate early and stay apart indicate a real, persistent difference in risk between the groups.",
" Censored subjects leave the curve without pulling it down — that is what makes these estimates honest."
),
chart_labels = list(
time_point = shared$h_time,
survival_pct = "% still active"
),
data = list(km_curves = shared$km_df)
)
}
# Card: median_survival (table)
card_median_survival <- function(shared, df, params) {
nr_groups <- shared$median_df$group[shared$median_df$median_time == "not reached"]
nr_note <- if (length(nr_groups) > 0) {
paste0(" For ", paste(nr_groups, collapse = ", "),
", the median is \"not reached\" — more than half of those subjects were ",
"still event-free at last observation, which itself signals strong retention.")
} else ""
list(
title = "Median Time to Event",
description = paste0("Per-", shared$h_group, " medians with 95% confidence intervals."),
text = paste0(
"The table shows each group's subjects, observed events, censoring rate, ",
"and median ", shared$h_time, " to event with a 95% confidence interval. ",
"The median is the point where half the group has had the event — a more ",
"robust summary than the mean for duration data.", nr_note
),
data = list(median_survival = shared$median_df)
)
}
# Card: test_results (table)
card_test_results <- function(shared, df, params) {
intro <- if (shared$single_group) {
paste0("Only one ", shared$h_group, " group is present, so no between-group ",
"test applies; the row below records that explicitly.")
} else {
paste0(
"The log-rank test asks whether the survival curves differ anywhere along ",
"their length; the Cox rows translate each group's risk into a hazard ratio ",
"versus ", shared$ref_group, " (the largest group). A hazard ratio of 1.5 ",
"means that group reaches the event 1.5× faster at any given moment. ",
if (!is.na(shared$logrank_p) && shared$logrank_p < 0.05)
paste0("Here the curves differ significantly(p = ", signif(shared$logrank_p, 3), ").")
else if (!is.na(shared$logrank_p))
paste0("Here the difference is not statistically significant(p = ", signif(shared$logrank_p, 3), ").")
else ""
)
}
list(
title = "Statistical Tests",
description = "Log-rank test and Cox hazard ratios with confidence intervals.",
text = intro,
data = list(test_results = shared$tests_df)
)
}
# Card: event_distribution (histogram)
card_event_distribution <- function(shared, df, params) {
et <- shared$event_times_df$event_time
q <- quantile(et, c(0.25, 0.5, 0.75), na.rm = TRUE)
concentration <- {
spread <- (q[3] - q[1]) / max(q[2], .Machine$double.eps)
if (spread < 0.8) "events concentrate tightly around the median"
else if (q[2] < mean(range(et)) * 0.6) "events skew early — most risk is front-loaded"
else "events spread across the observed timeline"
}
list(
title = "When Events Happen",
description = paste0("Distribution of observed event times(censored subjects excluded)."),
text = paste0(
"Among the ", format(length(et), big.mark = ","), " observed events, half ",
"happened within ", .fmt_time(q[2]), " ", shared$h_time,
" (quartiles: ", .fmt_time(q[1]), " / ", .fmt_time(q[2]), " / ", .fmt_time(q[3]),
"). In this data, ", concentration, ". Note this histogram shows only ",
"subjects whose event was observed — the ", round(shared$censored_pct, 1),
"% censored are visible in the survival curves, not here."
),
chart_labels = list(event_time = shared$h_time),
data = list(event_distribution = shared$event_times_df)
)
}