When the experiment defines the event

In the algae experiment that motivates this tutorial, three treatment groups (Control, Low Pollution, High Pollution) were monitored over 72 hours. The “event” can be defined in two ways:

1. End-of-experiment censoring — every culture is censored at hour 72, because the experiment ended before all cultures could show a clear biological endpoint. The Kaplan–Meier curves in this design show relative patterns of continued growth.

2. Threshold crossing — define a biologically meaningful absorbance level (e.g., OD ≥ 0.6). Cultures that reach it are “events”; those that do not are censored at their last measurement. This approach requires a defensible threshold.

Both designs are valid. This tutorial implements the threshold approach, which produces the most interpretable Kaplan–Meier curves and Cox model.

The survival function S(t)

The survival function S(t) gives the probability that the event has not yet occurred by time t:

\[S(t) = P(T > t)\]

where T is the random variable representing time-to-event. S(0) = 1 (no events at the start) and S(t) is non-increasing over time.

Censoring

An observation is right-censored at time c if the event had not occurred by the time the study ended or the subject left observation. The censored observation tells us only that T > c — the event would have happened eventually, but we did not observe it.

In R’s Surv() function, status = 1 indicates an observed event; status = 0 indicates censoring.

The hazard function h(t)

The hazard h(t) is the instantaneous rate at which events occur among those who have not yet had the event. A Cox proportional hazards model estimates how covariates (e.g., treatment group) scale the baseline hazard:

\[h(t \mid x) = h_0(t) \cdot \exp(\beta_1 x_1 + \beta_2 x_2 + \ldots)\]

The hazard ratio HR = exp(β) compares event rates between groups. HR > 1 means the treatment group reaches the threshold faster; HR < 1 means it reaches it more slowly.

Required packages

library(survival)   # Surv(), survfit(), coxph(), cox.zph()
library(survminer)  # ggsurvplot(), pairwise_survdiff()
library(dplyr)      # data wrangling
library(tidyr)
library(broom)      # tidy() for Cox output
library(ggplot2)

Simulated algae growth data

The data below simulate a 72-hour algae growth experiment with 10 biological replicates per treatment. Growth follows a logistic curve with treatment-specific rates and noise.

set.seed(42)

# Simulate one replicate
sim_culture <- function(sample_id, treatment, times = seq(0, 72, by = 8)) {
  params <- list(
    A = list(K = 1.2, r = 0.08, n = 0.06),
    B = list(K = 1.0, r = 0.06, n = 0.09),
    C = list(K = 0.7, r = 0.04, n = 0.12)
  )[[treatment]]
  
  abs_val <- params$K / (1 + exp(-params$r * (times - 36))) +
             rnorm(length(times), 0, params$n)
  abs_val <- pmax(abs_val, 0)
  
  data.frame(sample = sample_id, treatment = treatment,
             time = times, absorbance = abs_val)
}

treatments <- c("A", "B", "C")
df_raw <- do.call(rbind, lapply(treatments, function(trt) {
  do.call(rbind, lapply(1:10, function(i) {
    sim_culture(paste0(trt, "_", i), trt)
  }))
}))
df_raw$treatment <- factor(df_raw$treatment, levels = c("A", "B", "C"))

Define the event: time to reach threshold

Choose a biologically meaningful absorbance threshold. Here OD ≥ 0.6 represents the point at which growth is considered well established. For each replicate, find the first time that threshold is reached. Replicates that never reach it are right-censored at their final measurement time.

threshold <- 0.6

time_to_event <- df_raw |>
  group_by(sample, treatment) |>
  arrange(time) |>
  summarise(
    event_time = {
      hit <- which(absorbance >= threshold)
      if (length(hit) == 0) NA_real_ else time[min(hit)]
    },
    last_time = max(time),
    .groups = "drop"
  ) |>
  mutate(
    status = if_else(is.na(event_time), 0L, 1L),  # 1 = event, 0 = censored
    time   = if_else(is.na(event_time), last_time, event_time)
  )

# Event count by treatment
time_to_event |> count(treatment, status)

## # A tibble: 4 × 3
##   treatment status     n
##   <fct>      <int> <int>
## 1 A              1    10
## 2 B              1    10
## 3 C              0     5
## 4 C              1     5

Fitting the model

The Kaplan–Meier estimator is a non-parametric method that estimates S(t) directly from the data without assuming any underlying distribution. It updates at each observed event time using:

\[\hat{S}(t) = \prod_{t_i \leq t} \left(1 - \frac{d_i}{n_i}\right)\]

where \(d_i\) is the number of events and \(n_i\) is the number still at risk at time \(t_i\).

km_fit <- survfit(Surv(time, status) ~ treatment, data = time_to_event)
summary(km_fit)

## Call: survfit(formula = Surv(time, status) ~ treatment, data = time_to_event)
## 
##                 treatment=A 
##  time n.risk n.event survival std.err lower 95% CI upper 95% CI
##    32     10       1      0.9  0.0949       0.7320        1.000
##    40      9       8      0.1  0.0949       0.0156        0.642
##    48      1       1      0.0     NaN           NA           NA
## 
##                 treatment=B 
##  time n.risk n.event survival std.err lower 95% CI upper 95% CI
##    40     10       3      0.7  0.1449       0.4665        1.000
##    48      7       6      0.1  0.0949       0.0156        0.642
##    56      1       1      0.0     NaN           NA           NA
## 
##                 treatment=C 
##  time n.risk n.event survival std.err lower 95% CI upper 95% CI
##    40     10       1      0.9  0.0949        0.732        1.000
##    48      9       1      0.8  0.1265        0.587        1.000
##    56      8       1      0.7  0.1449        0.467        1.000
##    64      7       1      0.6  0.1549        0.362        0.995
##    72      6       1      0.5  0.1581        0.269        0.929

Plot with risk table

ggsurvplot(
  km_fit,
  data        = time_to_event,
  pval        = TRUE,
  conf.int    = TRUE,
  risk.table  = TRUE,
  palette     = c("#4a7c5e", "#c8832a", "#8b3a2a"),
  ggtheme     = theme_minimal(base_size = 12),
  title       = "Kaplan–Meier: time to reach absorbance threshold",
  xlab        = "Time (hours)",
  ylab        = "Probability of not yet reaching threshold",
  legend.title = "Treatment",
  risk.table.height = 0.25
)

Kaplan–Meier survival curves showing the probability that each culture has not yet reached OD = 0.6. Shading = 95% CI. The risk table below shows how many replicates remain uncensored at each time point.

Interpreting the curves. S(t) is the probability that a replicate has not yet reached OD 0.6 by time t. A curve that drops steeply early means the treatment reached the threshold quickly; a flat curve means most replicates were still below threshold (censored). Treatment A should show the steepest drop; treatment C the flattest, because pollution suppresses growth.

Pairwise comparisons (adjusted)

If the omnibus test is significant, follow up with pairwise comparisons. The Benjamini–Hochberg (BH) method controls the false discovery rate.

pairwise_res <- pairwise_survdiff(
  Surv(time, status) ~ treatment,
  data = time_to_event,
  p.adjust.method = "BH"
)
pairwise_res

## 
##  Pairwise comparisons using Log-Rank test 
## 
## data:  time_to_event and treatment 
## 
##   A       B      
## B 0.00657 -      
## C 0.00027 0.00148
## 
## P value adjustment method: BH

Extracting and plotting hazard ratios

hr <- tidy(cox1, exponentiate = TRUE, conf.int = TRUE)
print(hr)

## # A tibble: 2 × 7
##   term       estimate std.error statistic   p.value conf.low conf.high
##   <chr>         <dbl>     <dbl>     <dbl>     <dbl>    <dbl>     <dbl>
## 1 treatmentB   0.305      0.490     -2.42 0.0154     0.117       0.797
## 2 treatmentC   0.0359     0.774     -4.30 0.0000171  0.00787     0.163

ggplot(hr, aes(x = term, y = estimate)) +
  geom_point(size = 3, color = "#4a7c5e") +
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high),
                width = 0.15, color = "#4a7c5e") +
  geom_hline(yintercept = 1, linetype = "dashed", color = "#c8832a") +
  coord_flip() +
  labs(y = "Hazard Ratio (HR)", x = "",
       title = "Cox model: HR and 95% CI (reference = Treatment A)") +
  theme_minimal(base_size = 13)

Cox model hazard ratios (HR) with 95% confidence intervals. HR > 1 means the treatment reaches the threshold faster than treatment A; HR < 1 means slower.

Interpreting hazard ratios. An HR of 2 for treatment C means that at any point in time, a culture in treatment C is twice as likely to reach the threshold as one in treatment A, conditional on not yet having done so. HR < 1 would indicate slower progression.

Parametric alternative: Weibull model

weibull_mod <- survreg(Surv(time, status) ~ treatment,
                       data = time_to_event, dist = "weibull")
summary(weibull_mod)

## 
## Call:
## survreg(formula = Surv(time, status) ~ treatment, data = time_to_event, 
##     dist = "weibull")
##               Value Std. Error      z      p
## (Intercept)  3.7195     0.0366 101.69 <2e-16
## treatmentB   0.1569     0.0511   3.07 0.0021
## treatmentC   0.5872     0.0638   9.20 <2e-16
## Log(scale)  -2.1694     0.1596 -13.59 <2e-16
## 
## Scale= 0.114 
## 
## Weibull distribution
## Loglik(model)= -86.1   Loglik(intercept only)= -109.4
##  Chisq= 46.64 on 2 degrees of freedom, p= 7.5e-11 
## Number of Newton-Raphson Iterations: 7 
## n= 30

Complete time-to-event table

knitr::kable(
  time_to_event[order(time_to_event$treatment, time_to_event$time), ],
  col.names = c("Sample", "Treatment", "Event time", "Last time", "Status", "Time used"),
  caption   = "Time-to-event data. Status: 1 = threshold reached, 0 = censored.",
  digits    = 1
)

What’s next?

• Logistic Regression — model binary outcomes with a generalized linear model when there is no time dimension
• Mixed-Effects Models — extend survival frameworks to nested or repeated-measures designs

library(survival)
data(lung)
km <- survfit(Surv(time, status) ~ sex, data = lung)
plot(km, col = c("#1e3a2f","#c8832a"), lwd = 2,
     xlab = "Days", ylab = "Survival probability",
     main = "Lung cancer survival by sex")
legend("topright", legend = c("Male","Female"),
       col = c("#1e3a2f","#c8832a"), lwd = 2)
survdiff(Surv(time, status) ~ sex, data = lung)