Survival analysis sounds like it's only for medical doctors tracking patient outcomes, but in the software world, we use these same patterns for "churn analysis"—basically, how long before a user cancels their subscription. The core challenge is "censoring." If a user is still active when you run your report, you don't know when they'll leave; you just know they haven't left yet. You can't just ignore them, or you'll bias your results. That's where the survival package comes in.

Setting up our patient recovery data

Let's build a small simulation. Imagine we're tracking how many days it takes for patients to recover from a specific surgery. We have two groups: those who received a standard treatment and those who received a new experimental drug. Some patients recovered, but others left the study early or the study ended before they recovered (these are our "censored" observations).

library(survival)

# Creating a small, messy dataset
set.seed(42)
recovery_data <- data.frame(
  patient_id = 1:100,
  days_to_event = sample(1:60, 100, replace = TRUE),
  status = sample(c(1, 0), 100, replace = TRUE, prob = c(0.7, 0.3)), # 1 = recovered, 0 = censored
  group = sample(c("Standard", "Experimental"), 100, replace = TRUE)
)

head(recovery_data)

Wrapping time and status in a Surv object

If you try to run a standard regression on days_to_event, you're lying to your model. You're treating the censored patients as if they recovered on day X, when in reality, they just stopped being tracked on day X. To fix this, we use the Surv() function. It creates a special object that tells R: "Here is the time, and here is whether the event actually happened."

# This is the magic sauce of the survival package
surv_obj <- Surv(time = recovery_data$days_to_event, event = recovery_data$status)
# If you print surv_obj, you'll see it looks like "60+" for censored cases.

Visualizing the Kaplan-Meier curve

I usually start with a Kaplan-Meier plot. It's a step function that shows the probability of "survival" (or in our case, not yet recovering) over time. It's the most honest way to see the difference between two groups.

# Fit the survival curve
km_fit <- survfit(surv_obj ~ group, data = recovery_data)

# Plotting it
plot(km_fit, col = c("red", "blue"), lty = 1:2, 
     xlab = "Days", ylab = "Probability of NOT yet recovered", 
     main = "Recovery Rates by Treatment Group")
legend("topright", legend = c("Experimental", "Standard"), col = c("red", "blue"), lty = 1:2)

The mistake: Trying to use linear regression

When I first started with survival data, I tried to use lm() to see if the group affected the recovery time. I thought, "It's just a number of days, right?" I wrote lm(days_to_event ~ group, data = recovery_data). The problem? The model treated a patient who was censored at day 10 (didn't recover) the same as a patient who recovered on day 10. This completely skewed the average recovery time downward.

You can't use lm() because survival data is rarely normally distributed and, more importantly, because linear regression can't handle the "plus" sign (the censoring) in our Surv object. We need a model that calculates the hazard—the instantaneous risk of the event happening.

Modeling the hazards with Cox Proportional Hazards

The Cox model is the industry standard here. Instead of predicting the exact day of recovery, it predicts the Hazard Ratio. If the ratio is 2.0, the Experimental group is recovering twice as fast as the Standard group at any given time.

# The coxph function takes the Surv object as the response
cox_model <- coxph(Surv(days_to_event, status) ~ group, data = recovery_data)

summary(cox_model)

When you look at the summary, focus on exp(coef). That's your Hazard Ratio. If it's greater than 1, the "event" (recovery) is happening faster in that group. If it's less than 1, it's slower. It's a much more robust way to handle time-to-event data than trying to force it into a linear model.

📋 Practical Task

Analyzing Churn Hazard Ratios for Software Tiers

You have been given a dataset of 200 users of a SaaS product. You need to determine if users on the "Enterprise" plan churn at a different rate than those on the "Basic" plan.

Your task:

1. Create a dataframe with 200 observations containing: tenure_days (random integers between 1 and 500), churned (binary 0 or 1), and plan ("Basic" or "Enterprise").
2. Using the survival package, create a Cox Proportional Hazards model to see if the plan significantly impacts the hazard of churning.
3. Extract the Hazard Ratio (exp(coef)) for the Enterprise plan.
4. Print a sentence stating whether Enterprise users are more or less likely to churn based on that ratio.