V The V Lab

Survival Analysis in Depth

From Research Questions, Censoring, and Nonparametric Estimation to Complex Event Processes and Transparent Reporting

The V Lab Home Biostatistics Foundations Medical Statistical Tests Linear Regression Logistic Regression Multilevel Modelling Causal Inference Propensity Score Matching Survival Analysis AFT and Cox PH 中文版

AudiencePublic health, medical, epidemiology, and health data science learners

Study timeApproximately 180–240 minutes

PrerequisitesDescriptive statistics, confidence intervals, regression fundamentals, and basic R

About the tutorial data Every participant, event time, censoring time, competing event, and treatment change is simulated with a fixed random seed and contains no real personal health information. The main cohort deliberately includes a time-varying intervention effect to demonstrate crossing survival curves and the limitations of the log-rank test and a single HR. All values are for teaching only and are not real clinical or causal evidence.

How to use this tutorial

This tutorial follows time from a defined origin to the first study endpoint, but it does not reduce survival analysis to a Kaplan–Meier curve or a Cox model. Work through this sequence:

Question and estimand → time structure and censoring → nonparametric description → group comparison → regression and diagnosis → complex event processes → design and reporting

Read each section’s research question before running its R code and interpreting the output. Knowledge checks and exercises are folded by default. The analytical code relies only on base R and the survival package; knitr is used solely to display tables on the teaching page.

Learning objectives

After completing this tutorial, you should be able to:

• translate “how long until the event?” into an explicit population, time origin, event, and target estimand;
• distinguish right, left, and interval censoring from left truncation/delayed entry;
• explain survival, hazard, and cumulative hazard functions and their relationships;
• calculate and interpret Kaplan–Meier estimates, numbers at risk, and Nelson–Aalen estimates;
• explain what the log-rank test evaluates, when it is powerful, and why it cannot replace effect estimation;
• express absolute differences with fixed-time survival probabilities and restricted mean survival time (RMST);
• identify nonproportional hazards and choose time-varying effects or another suitable estimand;
• construct delayed-entry and time-varying-covariate data correctly;
• distinguish cause-specific hazards, cumulative incidence functions, and Aalen–Johansen state probabilities;
• identify the data structures needed for competing risks, recurrent events, and multistate processes;
• incorporate study design, loss to follow-up, missing data, and causal interpretation into an analysis plan;
• report absolute and relative quantities, uncertainty, diagnostic findings, and limitations together.

1 From the Question to the Target Estimand

1.1 “Time to event” is not a complete research question

Survival analysis considers both whether an event occurs and when it occurs. Before doing any calculations, write the following elements into the study protocol:

Element	Question to clarify	Definition in this tutorial’s main cohort
Target population	To whom should the results generalize?	Adults who meet the simulated study’s eligibility criteria
Time zero	When does a person first become at risk?	Date of randomized assignment to a management strategy
Time scale	Follow-up months, age, or calendar time?	Months since assignment
Event	Which reproducibly assessed endpoint?	First occurrence of the simulated study endpoint
Competing event	What would prevent the event of interest from occurring?	None in the main cohort; simulated separately below
End of observation	When does observation stop?	Event, loss to follow-up, or administrative cutoff
Unit of analysis	Whom or what does each row represent?	One participant per row in the main analysis

Time zero, eligibility assessment, and strategy assignment should be aligned whenever possible. If people must survive until treatment to be classified as treated, but follow-up begins at an earlier diagnosis date, the treated group receives a period during which survival is guaranteed. This creates immortal time bias.

1.2 Choose the estimand before choosing the model

The same research question can be answered on different scales:

Target estimand	Question answered	Typical expression
$S(12)$	What is the probability of remaining event-free at 12 months?	12-month survival probability and between-group difference
$F(12)=1-S(12)$	With no competing risk, what is the probability of an event by 12 months?	12-month cumulative event risk
Median event time	When have 50% of people experienced the event?	Months; may be “not reached”
RMST($\tau$)	How much event-free time is expected on average through $\tau$?	Mean event-free months through $\tau$ and the between-group difference
HR	How do conditional instantaneous event rates differ among people still at risk?	Cox PH or time-specific HR
TR	By what factor are conditional event-time quantiles multiplied?	AFT time ratio
Cumulative incidence function $F_k(t)$	With competing events, what is the real-world probability of the event of interest by $t$?	Cause-specific cumulative incidence probability

These estimands are not interchangeable versions of the same answer. An HR is not a risk ratio, a TR is not an HR, and an RMST difference is not uniquely determined by a single HR. The estimand that best supports the decision should be prespecified before examining the results.

Do not let a p value choose your research question Do not search across log-rank, Cox, AFT, and RMST analyses for the smallest p value and then call that scale the “primary result.” Doing so changes the question and amplifies selective reporting. The scientific question determines the primary estimand; other scales can serve as prespecified supplementary or sensitivity analyses.

Check your understanding: Is an HR the same as 12-month risk?

If a Cox model gives HR=0.70, can you write directly that “12-month event risk was reduced by 30%”?

Answer: No. The HR compares conditional instantaneous event rates among people still at risk; 12-month risk is a probability accumulated from time zero through 12 months. Calculate the 12-month risk or risk difference from the corresponding survival probabilities.

2 Time, Events, and Censoring

2.1 The complete event time is not always observed

Let the true event time be $T$ and the right-censoring time be $C$. We usually observe:

$$\widetilde T=\min(T,C), \qquad \delta=I(T\le C).$$

When $\delta=1$, the event is observed. When $\delta=0$, we know only that the true event time exceeds the observed time. A censored participant is not “event-free forever,” and the event time must not be recorded as infinity. Right-censoring methods retain the person’s event-free follow-up before censoring and remove that person from the risk set afterward.

2.1.1 Common observation structures

Structure	Information known	Example	Analysis reminder
Right censoring	$T>C$	No recurrence by the end of the study	Commonly Surv(time, event)
Left censoring	$T\le L$	Antibodies are already present at the first test	Not the same as delayed entry
Interval censoring	$L<T\le R$	Seroconversion occurs between two screening visits	Requires an interval-censoring likelihood
Left truncation/delayed entry	Only people satisfying $T>E$ are observed	Entry into a registry several months after diagnosis	The risk set begins at entry

Independent censoring does not require everyone to have the same censoring time. It means that, conditional on the information included in the analysis, the latent event time and censoring mechanism carry no additional information about one another. If people who deteriorate rapidly are more likely to be lost to follow-up and that deterioration is unrecorded, ordinary KM or Cox analyses may be biased.

2.2 Audit event coding and follow-up completeness

For a numeric status, Surv(time, event) usually uses 1=event and 0=censored. Check this explicitly before modeling instead of relying on software to infer the coding.

stopifnot(
  all(study_data$time_months > 0),
  all(study_data$event %in% c(0, 1)),
  !anyNA(study_data)
)

preview <- transform(
  head(study_data[, c(
    "participant_id", "time_months", "event",
    "treatment", "age", "severity", "biomarker"
  )]),
  event = ifelse(event == 1, "Event", "Right-censored"),
  treatment = treatment_labels[as.character(treatment)],
  severity = severity_labels[as.character(severity)]
)

knitr::kable(
  preview,
  col.names = c(
    "Participant", "Observed time (months)", "Observed endpoint", "Management strategy",
    "Age", "Baseline severity", "Standardized biomarker"
  ),
  caption = "First six rows of the simulated main cohort"
)

First six rows of the simulated main cohort

Participant	Observed time (months)	Observed endpoint	Management strategy	Age	Baseline severity	Standardized biomarker
P001	1.06	Event	Standard management	73	Mild	0.59
P002	4.34	Event	Standard management	48	Mild	-1.08
P003	10.17	Event	Intervention	37	Mild	-0.63
P004	19.39	Right-censored	Standard management	49	Moderate	2.11
P005	4.24	Event	Standard management	56	Mild	0.62
P006	22.77	Right-censored	Standard management	49	Severe	-1.39

followup_audit <- data.frame(
  Sample_size = nrow(study_data),
  Events = sum(study_data$event),
  Right_censored = sum(study_data$event == 0),
  Right_censored_percentage = pct(mean(study_data$event == 0)),
  Shortest_observation_months = min(study_data$time_months),
  Longest_observation_months = max(study_data$time_months),
  check.names = FALSE
)

knitr::kable(
  followup_audit,
  digits = 2,
  col.names = c(
    "Sample size", "Events", "Right-censored", "Right-censored (%)",
    "Shortest observation (months)", "Longest observation (months)"
  ),
  caption = "Follow-up audit for the main cohort"
)

Follow-up audit for the main cohort

Sample size	Events	Right-censored	Right-censored (%)	Shortest observation (months)	Longest observation (months)
720	575	145	20.1%	0.12	29.8

When describing follow-up, report at least the sample size, event count, censoring count, observation range, and numbers at risk at important time points. Reporting only “median follow-up” hides when censoring occurred and how many people support late estimates.

3 Survival, Hazard, and Cumulative Hazard

3.1 Three related but noninterchangeable functions

For a continuous event time $T$:

Survival function: $S(t)=P(T>t)$

Hazard function: $h(t)=\lim_{\Delta t\to 0}\dfrac{P(t\le T<t+\Delta t\mid T\ge t)}{\Delta t}$

Cumulative hazard: $H(t)=\int_0^t h(u)\,du$

They satisfy:

$$S(t)=\exp\{-H(t)\}, \qquad H(t)=-\log S(t).$$

• $S(t)$ is the probability of remaining event-free from time zero through $t$ and lies between 0 and 1;
• $h(t)$ is the instantaneous event rate immediately after $t$, conditional on still being at risk at $t$. It has units of “per unit time” and need not be less than 1;
• $H(t)$ accumulates the hazard over time. It is not a probability and may exceed 1.

A hazard is not an individual’s “probability of having the event right now” The hazard is a population-level conditional instantaneous rate. It depends on the composition of the people who remain in the risk set; as higher-risk people experience events earlier, the later risk set becomes selected. Thus, even though the HR has a precise mathematical definition, it should not be given a simplistic causal interpretation without reference to the risk set, time horizon, and study design.

4 Nonparametric Estimation: Let the Data Speak First

4.1 Kaplan–Meier estimation and numbers at risk

At each distinct event time $t_j$, let $d_j$ be the number of events and $n_j$ the number at risk just before that time. The Kaplan–Meier (KM) estimator is:

$$\widehat S_{KM}(t)= \prod_{t_j\le t}\left(1-\frac{d_j}{n_j}\right).$$

Events make the curve step downward. Censoring does not itself lower the curve, but it reduces the number at risk afterward. A flat segment of a KM curve does not establish “no risk”; it only means that no event was observed during that interval.

survival_outcome <- with(
  study_data,
  survival::Surv(time_months, event)
)

km_fit <- survival::survfit(
  survival_outcome ~ treatment,
  data = study_data,
  conf.type = "log-log"
)

plot(
  km_fit,
  col = c(palette_sa["orange"], palette_sa["teal"]),
  lwd = 2.4,
  mark.time = TRUE,
  conf.int = FALSE,
  xlab = "Months since assignment",
  ylab = "Estimated probability of remaining event-free",
  xlim = c(0, 30),
  ylim = c(0, 1),
  las = 1
)
abline(v = 12, lty = 3, col = palette_sa["gray"])
legend(
  "topright",
  legend = c(
    "Standard management", "Intervention",
    "Prespecified 12-month cutoff"
  ),
  col = c(
    palette_sa["orange"],
    palette_sa["teal"],
    palette_sa["gray"]
  ),
  lty = c(1, 1, 3),
  lwd = c(2.4, 2.4, 1.2),
  bty = "n"
)

Kaplan–Meier survival curves by management strategy. Short vertical marks indicate right censoring. The curves separate early, then converge and cross, suggesting that a single proportional effect may be inappropriate.

The graph must be read together with the numbers at risk. If only a few participants remain near the tail of a curve, one event can produce a large step and the confidence interval will widen.

km_selected <- summary(
  km_fit,
  times = c(6, 12, 18, 24)
)

km_selected_table <- data.frame(
  Management_strategy = treatment_labels[
    sub("treatment=", "", as.character(km_selected$strata))
  ],
  Time_months = km_selected$time,
  Number_at_risk = km_selected$n.risk,
  Survival_probability = km_selected$surv,
  Lower_95_CI = km_selected$lower,
  Upper_95_CI = km_selected$upper,
  check.names = FALSE
)

knitr::kable(
  km_selected_table,
  digits = 3,
  col.names = c(
    "Management strategy", "Time (months)", "Number at risk",
    "Survival probability", "Lower 95% CI", "Upper 95% CI"
  ),
  caption = paste(
    "KM survival probabilities, numbers at risk, and 95% confidence",
    "intervals at prespecified time points"
  )
)

KM survival probabilities, numbers at risk, and 95% confidence intervals at prespecified time points

Management strategy	Time (months)	Number at risk	Survival probability	Lower 95% CI	Upper 95% CI
Standard management	6	226	0.630	0.577	0.677
Standard management	12	157	0.437	0.386	0.488
Standard management	18	88	0.278	0.233	0.326
Standard management	24	32	0.216	0.172	0.263
Intervention	6	290	0.803	0.758	0.841
Intervention	12	231	0.640	0.588	0.687
Intervention	18	94	0.286	0.239	0.333
Intervention	24	22	0.155	0.115	0.201

km_risk_summary <- summary(
  km_fit,
  times = c(0, 6, 12, 18, 24)
)
km_risk_long <- data.frame(
  Management_strategy = treatment_labels[
    sub("treatment=", "", as.character(km_risk_summary$strata))
  ],
  Time_months = km_risk_summary$time,
  Number_at_risk = km_risk_summary$n.risk
)
km_risk_table <- xtabs(
  Number_at_risk ~ Management_strategy + Time_months,
  data = km_risk_long
)

knitr::kable(
  km_risk_table,
  caption = "Number still in the risk set just before each time point"
)

Number still in the risk set just before each time point

	0	6	12	18	24
Intervention	361	290	231	94	22
Standard management	359	226	157	88	32

km_median_raw <- summary(km_fit)$table
km_median_table <- data.frame(
  Management_strategy = treatment_labels[
    sub("treatment=", "", rownames(km_median_raw))
  ],
  Median_event_time_months = km_median_raw[, "median"],
  Lower_95_CI = km_median_raw[, "0.95LCL"],
  Upper_95_CI = km_median_raw[, "0.95UCL"],
  check.names = FALSE
)

knitr::kable(
  km_median_table,
  digits = 2,
  col.names = c(
    "Management strategy", "Median event time (months)",
    "Lower 95% CI", "Upper 95% CI"
  ),
  caption = "Unadjusted KM median event time by management strategy"
)

Unadjusted KM median event time by management strategy

	Management strategy	Median event time (months)	Lower 95% CI	Upper 95% CI
Standard	Standard management	10.3	8.42	11.8
Intervention	Intervention	14.0	13.16	15.2

If a curve remains above 0.50 throughout observation, report the median event time as “not reached”; do not substitute the final observed time. Curves should also not be compared using only their two medians, because doing so ignores the rest of follow-up.

4.2 Nelson–Aalen cumulative hazard estimation

The Nelson–Aalen estimator sums hazard increments at each event time:

$$\widehat H_{NA}(t)= \sum_{t_j\le t}\frac{d_j}{n_j}.$$

KM directly estimates the survival function, whereas Nelson–Aalen directly estimates the cumulative hazard. The exponential relationship between survival and cumulative hazard connects the two scales. When event increments are small, the resulting survival estimates are usually close, but they are not algebraically identical.

na_fit <- survival::survfit(
  survival_outcome ~ treatment,
  data = study_data,
  ctype = 1
)

plot(
  na_fit,
  fun = "cumhaz",
  col = c(palette_sa["orange"], palette_sa["teal"]),
  lwd = 2.4,
  conf.int = FALSE,
  xlab = "Months since assignment",
  ylab = "Nelson–Aalen cumulative hazard",
  xlim = c(0, 30),
  las = 1
)
abline(v = 12, lty = 3, col = palette_sa["gray"])
legend(
  "topleft",
  legend = c("Standard management", "Intervention"),
  col = c(palette_sa["orange"], palette_sa["teal"]),
  lwd = 2.4,
  bty = "n"
)

Nelson–Aalen cumulative hazard curves by management strategy. Their slopes reflect the pace of event accumulation; the intervention curve becomes steeper later, consistent with the rapid late decline in its KM curve.

na_selected <- summary(
  na_fit,
  times = c(6, 12, 18, 24)
)

na_table <- data.frame(
  Management_strategy = treatment_labels[
    sub("treatment=", "", as.character(na_selected$strata))
  ],
  Time_months = na_selected$time,
  Cumulative_hazard = na_selected$cumhaz,
  Standard_error = na_selected$std.chaz,
  check.names = FALSE
)

knitr::kable(
  na_table,
  digits = 3,
  col.names = c(
    "Management strategy", "Time (months)",
    "Cumulative hazard", "Standard error"
  ),
  caption = "Nelson–Aalen cumulative hazard at selected time points"
)

Nelson–Aalen cumulative hazard at selected time points

Management strategy	Time (months)	Cumulative hazard	Standard error
Standard management	6	0.462	0.040
Standard management	12	0.825	0.060
Standard management	18	1.275	0.085
Standard management	24	1.525	0.107
Intervention	6	0.219	0.026
Intervention	12	0.446	0.039
Intervention	18	1.249	0.084
Intervention	24	1.854	0.141

Check your understanding: Does cumulative hazard=1 mean a 100% event probability?

Answer: No. Cumulative hazard is not a probability. If $H(t)=1$, the corresponding survival probability is $S(t)=e^{-1}$, approximately 0.37; without competing risks, the cumulative event probability is approximately 0.63.

5 Between-Group Comparisons and Absolute Effects

5.1 What does the log-rank test evaluate?

At each event time, the standard log-rank test compares the observed event count in each group with the count expected under the null hypothesis, then standardizes the difference by its variance. For two groups, it can be summarized as:

$$Q=\frac{(O_1-E_1)^2}{\operatorname{Var}(O_1-E_1)},$$

Under the null hypothesis and the corresponding censoring assumptions, $Q$ is approximately chi-squared with 1 degree of freedom. The log-rank test compares the full survival processes but does not estimate an effect size. It is generally powerful for proportional-hazards-type differences, but early and late differences may cancel when curves cross.

logrank_fit <- survival::survdiff(
  survival::Surv(time_months, event) ~ treatment,
  data = study_data,
  rho = 0
)

logrank_df <- length(logrank_fit$n) - 1
logrank_p <- pchisq(
  logrank_fit$chisq,
  df = logrank_df,
  lower.tail = FALSE
)

logrank_table <- data.frame(
  Management_strategy = treatment_labels[
    sub("treatment=", "", names(logrank_fit$n))
  ],
  Sample_size = as.numeric(logrank_fit$n),
  Observed_events = as.numeric(logrank_fit$obs),
  Expected_events_under_null = as.numeric(logrank_fit$exp),
  check.names = FALSE
)

knitr::kable(
  logrank_table,
  digits = 1,
  col.names = c(
    "Management strategy", "Sample size", "Observed events",
    "Expected events under the null"
  ),
  caption = "Observed and expected event counts in the log-rank test"
)

Observed and expected event counts in the log-rank test

	Management strategy	Sample size	Observed events	Expected events under the null
Standard	Standard management	359	279	263
Intervention	Intervention	361	296	312

logrank_result <- data.frame(
  Chi_square_statistic = unname(logrank_fit$chisq),
  Degrees_of_freedom = logrank_df,
  p_value = format_p(logrank_p),
  check.names = FALSE
)

knitr::kable(
  logrank_result,
  digits = 3,
  col.names = c("Chi-square statistic", "Degrees of freedom", "p-value"),
  caption = "Standard log-rank test"
)

Standard log-rank test

Chi-square statistic	Degrees of freedom	p-value
1.78	1	0.182

In this example, the log-rank p-value is 0.182. This does not mean that the groups are the same at all times: the intervention curve is substantially higher at 12 months, whereas the difference has reversed by 24 months. The standard log-rank test combines these oppositely directed differences, so they may cancel.

5.1.1 Limits of the log-rank test

• The test result is not an HR, risk difference, or difference in survival time.
• A large p-value does not prove that the curves are identical; it may also reflect inadequate sample size or an effect that changes sign over time.
• A small p-value does not establish clinical importance or identify which period contributed most.
• Within each group, censoring should not carry unmodeled prognostic information.
• Weighted log-rank tests can emphasize early or late events, but weights should be prespecified rather than searched for significance.
• When covariate adjustment, delayed entry, clustering, or time-varying effects matter, use an appropriate model instead of relying only on an unadjusted test.

5.2 Fixed-time survival probabilities

Fixed-time results directly answer “What is the probability of remaining event-free at this time?” This example shows why time points must be prespecified rather than selected wherever the difference is largest.

fixed_time_table <- km_selected_table[
  km_selected_table$Time_months %in% c(12, 24),
]

knitr::kable(
  fixed_time_table,
  digits = 3,
  col.names = c(
    "Management strategy", "Time (months)", "Number at risk",
    "Survival probability", "Lower 95% CI", "Upper 95% CI"
  ),
  caption = "Unadjusted KM survival probabilities at 12 and 24 months"
)

Unadjusted KM survival probabilities at 12 and 24 months

	Management strategy	Time (months)	Number at risk	Survival probability	Lower 95% CI	Upper 95% CI
2	Standard management	12	157	0.437	0.386	0.488
4	Standard management	24	32	0.216	0.172	0.263
6	Intervention	12	231	0.640	0.588	0.687
8	Intervention	24	22	0.155	0.115	0.201

At 12 months, the KM survival probabilities are 43.7% for standard management and 64.0% for intervention. By 24 months, their ordering has reversed. This heterogeneity over time cannot be summarized by one “overall best” percentage.

5.3 Restricted mean survival time

Restricted mean survival time (RMST) is the area under the survival curve from 0 to a prespecified cutoff $\tau$:

$$\operatorname{RMST}(\tau)= \int_0^\tau S(t)\,dt.$$

It can be interpreted as the mean event-free time through $\tau$. An RMST difference is expressed in the original time units and does not require proportional hazards, but it depends on $\tau$. The cutoff should be prespecified from the clinical question and the range of follow-up supported in both groups.

rmst_tau <- 24
rmst_raw <- summary(km_fit, rmean = rmst_tau)$table
rmst_row_order <- match(
  paste0("treatment=", levels(study_data$treatment)),
  rownames(rmst_raw)
)
rmst_raw <- rmst_raw[rmst_row_order, , drop = FALSE]

rmst_group_table <- data.frame(
  Management_strategy = treatment_labels[levels(study_data$treatment)],
  RMST_months = rmst_raw[, "rmean"],
  Standard_error = rmst_raw[, "se(rmean)"],
  check.names = FALSE
)

rmst_difference <-
  rmst_group_table$RMST_months[2] - rmst_group_table$RMST_months[1]
rmst_difference_se <- sqrt(sum(rmst_group_table$Standard_error^2))
rmst_difference_ci <- rmst_difference +
  c(-1, 1) * qnorm(0.975) * rmst_difference_se
rmst_difference_p <- 2 * pnorm(
  -abs(rmst_difference / rmst_difference_se)
)

rmst_difference_table <- data.frame(
  Contrast = "Intervention - standard management",
  RMST_difference_at_24_months = rmst_difference,
  Lower_95_CI = rmst_difference_ci[1],
  Upper_95_CI = rmst_difference_ci[2],
  p_value = format_p(rmst_difference_p),
  check.names = FALSE
)

knitr::kable(
  rmst_group_table,
  digits = 2,
  col.names = c(
    "Management strategy", "24-month RMST (months)", "Standard error"
  ),
  caption = "Restricted mean event-free time through 24 months from the KM curves"
)

Restricted mean event-free time through 24 months from the KM curves

	Management strategy	24-month RMST (months)	Standard error
Standard	Standard management	11.6	0.45
Intervention	Intervention	13.7	0.38

knitr::kable(
  rmst_difference_table,
  digits = 3,
  col.names = c(
    "Contrast", "24-month RMST difference (months)",
    "Lower 95% CI", "Upper 95% CI", "p-value"
  ),
  caption = paste(
    "Between-group difference in 24-month RMST with large-sample",
    "approximate uncertainty"
  )
)

Between-group difference in 24-month RMST with large-sample approximate uncertainty

Contrast	24-month RMST difference (months)	Lower 95% CI	Upper 95% CI	p-value
Intervention - standard management	2.08	0.925	3.23	<0.001

Through 24 months, the intervention group’s estimated mean event-free time is 2.08 months longer than that of the standard-management group (95% CI, 0.92 to 3.23). This does not conflict with the log-rank p-value. The RMST difference is an effect estimate based on the area between curves, whereas log-rank is a weighted test of the full curves; with crossing curves, the two focus on different information.

RMST is not an “assumption-free” method RMST still depends on the event and censoring definitions, a common time horizon, and the censoring mechanism. Interpretation remains unreliable if almost no one remains at risk by the cutoff $\tau$, or if a cutoff is chosen only because it yields significance. Adjusted RMST additionally requires standardization, weighting, or a suitable survival model.

6 Overview of Regression Models and PH Diagnostics

6.1 Cox PH and AFT answer different questions

Model	Basic form	Primary effect measure	Core effect assumption
Cox PH	$h(t\mid X)=h_0(t)e^{X^T\beta}$	HR $=e^\beta$	The HR is constant over time
AFT	$\log T=X^T\beta+\sigma\varepsilon$	TR $=e^\beta$	The conditional event-time distribution is scaled by a fixed factor

The Cox model does not require a parametric form for the baseline hazard, but it still relies on proportional hazards, appropriate covariate functional forms, and conditionally independent censoring. Parametric AFT models directly provide time ratios and event-time quantiles, but require choosing an event-time distribution. For systematic derivations, parameterizations, diagnostics, and conversions, see the AFT and Cox PH module. This integrated module focuses on placing these models within a complete analysis workflow.

cox_fit <- survival::coxph(
  survival::Surv(time_months, event) ~
    treatment + age10 + severity + biomarker,
  data = study_data,
  ties = "efron",
  x = TRUE,
  y = TRUE
)

cox_summary <- summary(cox_fit)
cox_term_labels <- c(
  treatmentIntervention = "Intervention vs standard management",
  age10 = "Age (per 10 years)",
  severityModerate = "Moderate vs mild severity",
  severitySevere = "Severe vs mild severity",
  biomarker = "Biomarker (per 1 SD)"
)
cox_terms <- rownames(cox_summary$coefficients)
cox_table <- data.frame(
  Variable = unname(cox_term_labels[cox_terms]),
  HR = cox_summary$coefficients[, "exp(coef)"],
  `Lower 95% CI` = cox_summary$conf.int[, "lower .95"],
  `Upper 95% CI` = cox_summary$conf.int[, "upper .95"],
  `p-value` = vapply(
    cox_summary$coefficients[, "Pr(>|z|)"],
    format_p,
    character(1)
  ),
  check.names = FALSE
)

knitr::kable(
  cox_table,
  digits = 3,
  caption = "Adjusted Cox model assuming a constant HR"
)

Adjusted Cox model assuming a constant HR

	Variable	HR	Lower 95% CI	Upper 95% CI	p-value
treatmentIntervention	Intervention vs standard management	0.858	0.727	1.01	0.0679
age10	Age (per 10 years)	1.187	1.102	1.28	<0.001
severityModerate	Moderate vs mild severity	1.179	0.983	1.41	0.075
severitySevere	Severe vs mild severity	1.523	1.208	1.92	<0.001
biomarker	Biomarker (per 1 SD)	1.145	1.056	1.24	0.00106

The intervention HR in this table is a single summary over the entire follow-up period. Because the KM curves already suggest that the effect changes direction, PH must be assessed before interpreting this estimate. A constant HR should not be assumed simply because the software successfully returns one.

6.2 Schoenfeld residuals and time-varying effects

cox.zph() assesses whether scaled Schoenfeld residuals vary systematically over time. A small p-value suggests that the corresponding log HR may change over time. A large p-value indicates only that there is no strong evidence against PH; it does not prove that PH holds.

ph_check <- survival::cox.zph(cox_fit, transform = "km")
ph_table <- data.frame(
  Test = c(
    "Management strategy", "Age", "Baseline severity",
    "Biomarker", "Global"
  ),
  Chi_square_statistic = ph_check$table[, "chisq"],
  Degrees_of_freedom = ph_check$table[, "df"],
  p_value = vapply(ph_check$table[, "p"], format_p, character(1)),
  check.names = FALSE
)

knitr::kable(
  ph_table,
  digits = 3,
  col.names = c(
    "Test", "Chi-square statistic", "Degrees of freedom", "p-value"
  ),
  caption = "Proportional hazards tests based on scaled Schoenfeld residuals"
)

Proportional hazards tests based on scaled Schoenfeld residuals

	Test	Chi-square statistic	Degrees of freedom	p-value
treatment	Management strategy	45.555	1	<0.001
age10	Age	3.760	1	0.0525
severity	Baseline severity	0.434	2	0.805
biomarker	Biomarker	1.023	1	0.312
GLOBAL	Global	50.879	5	<0.001

plot(
  ph_check,
  var = 1,
  resid = TRUE,
  se = TRUE,
  col = palette_sa["teal"]
)
abline(h = 0, lty = 2, col = palette_sa["gray"])
title("Time-varying effect of management strategy")

Scaled Schoenfeld residual diagnostic for management strategy. The smooth curve departs substantially from a horizontal line, indicating that a single intervention log HR does not summarize the full follow-up period.

In this example, both the management-strategy test and the global test show clear departures from PH. The next model estimates early and late HRs using the prespecified 12-month mechanism. In a real analysis, the cutoff must be prespecified based on the protocol or a clinical mechanism; repeatedly searching the curves for the most significant cutoff is not valid.

# survSplit() requires Surv to appear directly as the function name
# on the left-hand side of the formula.
Surv <- survival::Surv
piecewise_data <- survival::survSplit(
  Surv(time_months, event) ~ .,
  data = study_data,
  cut = 12,
  start = "tstart",
  end = "time_months"
)
rm(Surv)
piecewise_data$after12 <- as.integer(piecewise_data$tstart >= 12)

piecewise_cox <- survival::coxph(
  survival::Surv(tstart, time_months, event) ~
    treatment_num + I(treatment_num * after12) +
    age10 + severity + biomarker + cluster(participant_id),
  data = piecewise_data,
  ties = "efron"
)

piecewise_names <- c(
  "treatment_num",
  "I(treatment_num * after12)"
)
piecewise_beta <- coef(piecewise_cox)[piecewise_names]
piecewise_vcov <- vcov(piecewise_cox)[piecewise_names, piecewise_names]
contrast_matrix <- rbind(early = c(1, 0), late = c(1, 1))
period_log_hr <- drop(contrast_matrix %*% piecewise_beta)
period_se <- sqrt(diag(
  contrast_matrix %*% piecewise_vcov %*% t(contrast_matrix)
))

period_hr_table <- data.frame(
  Period = c("0 to ≤12 months", ">12 months"),
  HR = exp(period_log_hr),
  `Lower 95% CI` = exp(period_log_hr - 1.96 * period_se),
  `Upper 95% CI` = exp(period_log_hr + 1.96 * period_se),
  check.names = FALSE
)

knitr::kable(
  period_hr_table,
  digits = 3,
  caption = paste(
    "Adjusted Cox model allowing the intervention effect",
    "to change at 12 months"
  )
)

Adjusted Cox model allowing the intervention effect to change at 12 months

	Period	HR	Lower 95% CI	Upper 95% CI
early	0 to ≤12 months	0.505	0.404	0.631
late	>12 months	1.826	1.393	2.394

The results show that the intervention is associated with a lower instantaneous event rate during the first 12 months, followed by a reversal in direction. These period-specific HRs represent the data-generating mechanism more faithfully than a single HR, but they should still be reported alongside fixed-time survival probabilities and RMST.

Check your understanding: Does a violation of PH automatically support an AFT model?

Answer: No. Failure of PH means only that a constant HR is inappropriate. An AFT model separately assumes that the conditional event-time distribution is scaled by a constant factor and also depends on the selected distribution. When survival curves cross, a simple AFT model should often be questioned as well.

7 Complex Time Structures

7.1 Delayed entry: Risk sets begin at entry

Registries often enroll individuals after time zero. Only people whose event time occurs after their entry time can be observed, a situation known as left truncation. When time is measured from a common origin, the correct outcome is Surv(entry, exit, event). An individual cannot enter the risk set before their entry time.

set.seed(20260816)
n_source <- 1500
lt_age <- round(pmin(pmax(rnorm(n_source, 62, 10), 30), 88))
lt_age10 <- (lt_age - 60) / 10
lt_treatment_num <- rbinom(n_source, 1, 0.50)
lt_treatment <- factor(
  lt_treatment_num,
  levels = 0:1,
  labels = c("Standard management", "Intervention")
)

# The intervention group enters later on average, illustrating the bias
# produced by an incorrect risk set.
entry_time <- runif(n_source, 0, 7) + 2.5 * lt_treatment_num
lt_event_time <- rexp(
  n_source,
  rate = 0.045 * exp(0.25 * lt_age10 - 0.35 * lt_treatment_num)
)
lt_censor_time <- runif(n_source, 18, 30)
exit_time <- pmin(lt_event_time, lt_censor_time)
observed_after_entry <- exit_time > entry_time

lt_data <- data.frame(
  id = seq_len(n_source),
  entry = entry_time,
  exit = exit_time,
  event = as.integer(lt_event_time <= lt_censor_time),
  treatment = lt_treatment,
  age10 = lt_age10
)[observed_after_entry, ]

correct_lt_cox <- survival::coxph(
  survival::Surv(entry, exit, event) ~ treatment + age10,
  data = lt_data
)
naive_lt_cox <- survival::coxph(
  survival::Surv(exit, event) ~ treatment + age10,
  data = lt_data
)

lt_comparison <- data.frame(
  Analysis = c(
    "Correct: entry–exit risk sets",
    "Incorrect: everyone assumed at risk from time 0"
  ),
  Intervention_HR = exp(c(
    coef(correct_lt_cox)["treatmentIntervention"],
    coef(naive_lt_cox)["treatmentIntervention"]
  )),
  check.names = FALSE
)

lt_risk_counts <- data.frame(
  Time_months = c(2, 5, 8, 12),
  Correct_risk_count = vapply(
    c(2, 5, 8, 12),
    function(t) sum(lt_data$entry < t & lt_data$exit >= t),
    numeric(1)
  ),
  Incorrect_risk_count = vapply(
    c(2, 5, 8, 12),
    function(t) sum(lt_data$exit >= t),
    numeric(1)
  ),
  check.names = FALSE
)

knitr::kable(
  lt_risk_counts,
  col.names = c(
    "Time (months)", "Correct risk count", "Incorrect risk count"
  ),
  caption = "Correct and incorrect risk counts under delayed entry"
)

Correct and incorrect risk counts under delayed entry

Time (months)	Correct risk count	Incorrect risk count
2	184	1229
5	621	1183
8	940	1076
12	921	921

knitr::kable(
  lt_comparison,
  digits = 3,
  col.names = c("Analysis", "Intervention HR"),
  caption = "Effect of delayed-entry handling on the adjusted intervention HR"
)

Effect of delayed-entry handling on the adjusted intervention HR

Analysis	Intervention HR
Correct: entry–exit risk sets	0.754
Incorrect: everyone assumed at risk from time 0	0.629

The incorrect analysis places people in the risk set before they have entered the study. Because entry times differ between groups in this example, the intervention appears excessively protective. Delayed-entry analyses also require that the relationship between the entry mechanism and the subsequent event process can be adequately handled conditional on the analysis variables.

7.2 Time-varying covariates: Use only information known at the time

When treatment begins during follow-up, coding “ever treated later” as a fixed baseline variable uses future information and incorrectly assigns the survival time required before treatment initiation to the treated group. Start–stop data split follow-up whenever treatment status changes, so each interval uses only the value known at the beginning of that interval.

set.seed(20260817)
n_tv <- 900
tv_base <- data.frame(
  id = seq_len(n_tv),
  age10 = (round(pmin(pmax(rnorm(n_tv, 60, 10), 30), 85)) - 60) / 10
)
will_start <- rbinom(n_tv, 1, 0.72)
tv_base$planned_start <- ifelse(
  will_start == 1,
  runif(n_tv, 3, 12),
  Inf
)
pre_rate <- 0.060 * exp(0.20 * tv_base$age10)
post_rate <- 0.032 * exp(0.20 * tv_base$age10)
pre_event <- rexp(n_tv, pre_rate)
post_event <- rexp(n_tv, post_rate)
tv_event_time <- ifelse(
  pre_event <= tv_base$planned_start,
  pre_event,
  tv_base$planned_start + post_event
)
tv_censor <- runif(n_tv, 15, 26)
tv_base$time <- pmin(tv_event_time, tv_censor)
tv_base$event <- as.integer(tv_event_time <= tv_censor)

tv_long <- survival::tmerge(
  data1 = tv_base,
  data2 = tv_base,
  id = id,
  tstart = 0,
  tstop = time,
  outcome = event(time, event)
)
tv_long <- survival::tmerge(
  data1 = tv_long,
  data2 = tv_base,
  id = id,
  on_treatment = tdc(planned_start)
)

tv_cox <- survival::coxph(
  survival::Surv(tstart, tstop, outcome) ~
    on_treatment + age10 + cluster(id),
  data = tv_long,
  ties = "efron"
)

tv_base$ever_treated <- as.integer(tv_base$planned_start < tv_base$time)
immortal_time_cox <- survival::coxph(
  survival::Surv(time, event) ~ ever_treated + age10,
  data = tv_base,
  ties = "efron"
)

tv_comparison <- data.frame(
  Analysis = c(
    "Correct: treatment as a time-varying covariate",
    "Incorrect: future treatment status used at baseline"
  ),
  HR = c(
    exp(coef(tv_cox)["on_treatment"]),
    exp(coef(immortal_time_cox)["ever_treated"])
  )
)

knitr::kable(
  head(tv_long[, c("id", "tstart", "tstop", "outcome", "on_treatment")], 8),
  digits = 2,
  caption = "Example start–stop structure for time-varying treatment"
)

Example start–stop structure for time-varying treatment

id	tstart	tstop	outcome	on_treatment
1	0.00	9.26	0	0
1	9.26	23.83	0	1
2	0.00	10.16	1	0
3	0.00	9.40	0	0
3	9.40	18.83	0	1
4	0.00	6.24	0	0
4	6.24	25.04	0	1
5	0.00	2.39	1	0

knitr::kable(
  tv_comparison,
  digits = 3,
  caption = paste(
    "Comparison of correct time updating with an",
    "immortal-time-biased analysis"
  )
)

Comparison of correct time updating with an immortal-time-biased analysis

	Analysis	HR
on_treatment	Correct: treatment as a time-varying covariate	0.554
ever_treated	Incorrect: future treatment status used at baseline	0.184

When one person contributes multiple rows, cluster(id) provides a cluster-robust variance estimate. It cannot correct unmeasured time-varying confounding. If current health status affects both subsequent treatment and the outcome, an ordinary time-varying Cox coefficient still does not automatically represent a causal treatment effect; methods such as marginal structural models may be needed.

7.3 Competing risks and Aalen–Johansen

If death prevents recurrence, death is a competing event for recurrence. The cumulative incidence function for the cause of interest is:

$$F_k(t)=P(T\le t,J=k)=\int_0^t S(u-)\,d\Lambda_k(u).$$

It depends on both the hazard for the cause of interest and the hazards for other causes. Computing 1-KM after treating competing deaths as ordinary censoring imagines that people who die could still recur later and will generally overestimate the real-world recurrence probability. In survival, Surv() uses a multistate representation for a factor status, and survfit() returns Aalen–Johansen state probabilities in pstate.

set.seed(20260818)
n_cr <- 700
cr_treatment_num <- rbinom(n_cr, 1, 0.50)
cr_treatment <- factor(
  cr_treatment_num,
  levels = 0:1,
  labels = c("Standard management", "Intervention")
)
recurrence_time <- rexp(
  n_cr,
  rate = 0.045 * ifelse(cr_treatment_num == 1, 0.68, 1)
)
death_time <- rexp(n_cr, rate = 0.026)
cr_censor_time <- runif(n_cr, 18, 30)
cr_time <- pmin(recurrence_time, death_time, cr_censor_time)
cr_code <- ifelse(
  recurrence_time <= death_time & recurrence_time <= cr_censor_time,
  1L,
  ifelse(death_time < recurrence_time & death_time <= cr_censor_time, 2L, 0L)
)

# The first factor level must denote censoring; the remaining levels
# represent distinct event states.
cr_status <- factor(
  cr_code,
  levels = 0:2,
  labels = c("Censor", "Recurrence", "Death")
)
cr_data <- data.frame(
  id = seq_len(n_cr),
  time = cr_time,
  status = cr_status,
  treatment = cr_treatment
)

aj_fit <- survival::survfit(
  survival::Surv(time, status) ~ treatment,
  data = cr_data
)
aj_summary <- summary(aj_fit, times = c(12, 24))
recurrence_column <- match("Recurrence", aj_summary$states)
stopifnot(!is.na(recurrence_column))

naive_recurrence_fit <- survival::survfit(
  survival::Surv(time, status == "Recurrence") ~ treatment,
  data = cr_data
)
naive_recurrence_summary <- summary(
  naive_recurrence_fit,
  times = c(12, 24)
)
stopifnot(
  identical(
    as.character(aj_summary$strata),
    as.character(naive_recurrence_summary$strata)
  ),
  identical(aj_summary$time, naive_recurrence_summary$time)
)

cr_table <- data.frame(
  Management_strategy = sub(
    "treatment=", "", as.character(aj_summary$strata)
  ),
  Time_months = aj_summary$time,
  Aalen_Johansen_recurrence_probability =
    aj_summary$pstate[, recurrence_column],
  Lower_95_CI = aj_summary$lower[, recurrence_column],
  Upper_95_CI = aj_summary$upper[, recurrence_column],
  Standard_error = aj_summary$std.err[, recurrence_column],
  Incorrect_1_minus_KM = 1 - naive_recurrence_summary$surv,
  check.names = FALSE
)

knitr::kable(
  cr_table,
  digits = 3,
  col.names = c(
    "Management strategy", "Time (months)",
    "Aalen–Johansen recurrence probability", "Lower 95% CI",
    "Upper 95% CI", "Standard error", "Incorrect 1-KM"
  ),
  caption = paste(
    "Aalen–Johansen recurrence probability versus 1-KM after treating",
    "death as censoring"
  )
)

Aalen–Johansen recurrence probability versus 1-KM after treating death as censoring

Management strategy	Time (months)	Aalen–Johansen recurrence probability	Lower 95% CI	Upper 95% CI	Standard error	Incorrect 1-KM
Standard management	12	0.371	0.323	0.425	0.026	0.425
Standard management	24	0.505	0.454	0.562	0.028	0.643
Intervention	12	0.259	0.217	0.309	0.023	0.300
Intervention	24	0.401	0.352	0.456	0.027	0.521

# Use numeric group indices to preserve the fitted objects' strata order
# rather than allowing split() to reorder groups by their labels.
aj_index <- split(
  seq_along(aj_fit$time),
  rep(seq_along(aj_fit$strata), as.numeric(aj_fit$strata))
)
naive_index <- split(
  seq_along(naive_recurrence_fit$time),
  rep(
    seq_along(naive_recurrence_fit$strata),
    as.numeric(naive_recurrence_fit$strata)
  )
)
group_colors <- c(palette_sa["orange"], palette_sa["teal"])
stopifnot(
  identical(names(aj_fit$strata), names(naive_recurrence_fit$strata)),
  length(aj_index) == length(group_colors)
)

plot(
  NA,
  xlim = c(0, 28),
  ylim = c(0, 0.75),
  xlab = "Follow-up (months)",
  ylab = "Cumulative incidence of recurrence",
  las = 1
)
for (i in seq_along(aj_index)) {
  idx <- aj_index[[i]]
  lines(
    c(0, aj_fit$time[idx]),
    c(0, aj_fit$pstate[idx, match("Recurrence", aj_fit$states)]),
    type = "s",
    col = group_colors[i],
    lwd = 2.4
  )
  naive_idx <- naive_index[[i]]
  lines(
    c(0, naive_recurrence_fit$time[naive_idx]),
    c(0, 1 - naive_recurrence_fit$surv[naive_idx]),
    type = "s",
    col = group_colors[i],
    lwd = 2.2,
    lty = 2
  )
}
legend(
  "topleft",
  legend = c(
    "Standard management: Aalen–Johansen",
    "Intervention: Aalen–Johansen",
    "Standard management: incorrect 1-KM",
    "Intervention: incorrect 1-KM"
  ),
  col = c(group_colors, group_colors),
  lty = c(1, 1, 2, 2),
  lwd = 2.2,
  bty = "n"
)

Cumulative incidence of recurrence in the presence of competing death. Solid lines show Aalen–Johansen estimates; dashed lines show 1-KM after treating death as ordinary censoring and are systematically higher.

Cause-specific Cox models answer, “How does the instantaneous rate for the cause of interest vary among people currently free of any event?” The cumulative incidence function answers, “What is the real-world probability of the event of interest by a given time when all competing processes coexist?” These are different scales, not interchangeable substitutes.

7.4 Overview of recurrent events and multistate processes

A single event does not always capture the process. Recurrent hospitalizations can be represented on a total-time or gap-time scale and require explicit treatment of event order, terminal events, and within-person dependence. “Disease-free → recurrence → death” is a multistate process. Aalen–Johansen can estimate state-occupation or transition probabilities, but regression effects require a clearly specified model for each transition.

# Andersen–Gill form: multiple rows per person with a cluster-robust
# variance by id.
recurrent_fit <- survival::coxph(
  survival::Surv(start, stop, event) ~ exposure + cluster(id),
  data = recurrent_long
)

# Multistate Surv: status is a factor whose first level denotes censoring
# or no transition. istate specifies the initial state; genuine
# multitransition processes usually use start–stop data.
multistate_fit <- survival::survfit(
  survival::Surv(start, stop, status) ~ group,
  data = multistate_long,
  id = id,
  istate = initial_state
)

This code skeleton is not a universal prescription. Andersen–Gill, event-order models, shared frailty models, and multistate models use different risk sets and estimands. First draw the allowed transition diagram, then determine which time interval each row represents.

8 Study Design, Missing Data, and Causal Interpretation

8.1 Survival models cannot repair design flaws

Design problem	Possible consequence	What to do before analysis
Misalignment of time zero and assignment	Immortal time or selection bias	Align eligibility, assignment, and the start of follow-up
Outcome-assessment frequency differs by group	Differential detection opportunity	Standardize follow-up or model the observation process
Loss to follow-up relates to unrecorded illness	Informative censoring	Collect reasons and use weighting or sensitivity analyses
Missing baseline covariates	Complete-case selection bias and lower precision	Describe patterns and plan multiple imputation
Treatment changes during follow-up	Exposure misclassification and time-varying confounding	Use time-updated data and define the causal strategy
Competing events are not distinguished	Ambiguous estimand	Distinguish cause-specific hazards from cumulative incidence

A multiple-imputation model should include the event indicator, follow-up information, important auxiliary variables, and variables related to missingness. Do not impute an unknown event time after censoring as if it were an ordinary missing continuous value. Event count, rather than total sample size, often places the tighter limit on model complexity; nonlinearities, interactions, and time-varying effects all require additional information.

8.2 An adjusted association is not automatically a causal effect

Causal interpretation also requires a well-defined intervention, comparison strategy, target population, grace period, treatment-switching rule, handling of loss to follow-up, and assumptions of consistency, exchangeability, and positivity. Baseline confounders should be chosen from temporal ordering and causal structure, not screened by univariable p-values; posttreatment variables may be mediators or colliders. See the Causal Inference module for a fuller framework.

Survival prediction also requires a fairness review A higher predicted risk may reflect differences in diagnostic opportunity, access to care, or structural inequity. Before deployment, compare censoring, calibration, and errors across relevant groups; examine variable meaning and potential harms; and explain how predictions will affect resource allocation.

9 Complete Applied Example

9.1 Integrating the analysis plan and results

Research question: In a simulated randomized cohort, how does the intervention affect the time distribution of the first study endpoint over 24 months?

The prespecified primary absolute estimand is the 24-month RMST difference. Supplementary results include KM survival probabilities at 12 and 24 months, the standard log-rank test, and adjusted HRs that allow the effect to change after 12 months. Censoring is independent by construction, and every result remains a teaching simulation only.

standard_s12 <- km_selected_table$Survival_probability[
  km_selected_table$Management_strategy == "Standard management" &
    km_selected_table$Time_months == 12
]
intervention_s12 <- km_selected_table$Survival_probability[
  km_selected_table$Management_strategy == "Intervention" &
    km_selected_table$Time_months == 12
]
standard_s24 <- km_selected_table$Survival_probability[
  km_selected_table$Management_strategy == "Standard management" &
    km_selected_table$Time_months == 24
]
intervention_s24 <- km_selected_table$Survival_probability[
  km_selected_table$Management_strategy == "Intervention" &
    km_selected_table$Time_months == 24
]

case_results <- data.frame(
  Result = c(
    "Sample and observed endpoints",
    "12-month survival probability",
    "24-month survival probability",
    "24-month RMST difference",
    "Log-rank test",
    "Adjusted HR, 0 to ≤12 months",
    "Adjusted HR, >12 months"
  ),
  Estimate_or_test = c(
    sprintf("n=%d; events=%d; right-censored=%d", nrow(study_data),
            sum(study_data$event), sum(study_data$event == 0)),
    sprintf("Standard management %.1f%%; intervention %.1f%%",
            100 * standard_s12, 100 * intervention_s12),
    sprintf("Standard management %.1f%%; intervention %.1f%%",
            100 * standard_s24, 100 * intervention_s24),
    sprintf("%.2f months (95%% CI %.2f to %.2f)",
            rmst_difference, rmst_difference_ci[1], rmst_difference_ci[2]),
    sprintf("Chi-square=%.2f; p=%s", logrank_fit$chisq, format_p(logrank_p)),
    sprintf("%.2f (95%% CI %.2f to %.2f)",
            period_hr_table$HR[1], period_hr_table$`Lower 95% CI`[1],
            period_hr_table$`Upper 95% CI`[1]),
    sprintf("%.2f (95%% CI %.2f to %.2f)",
            period_hr_table$HR[2], period_hr_table$`Lower 95% CI`[2],
            period_hr_table$`Upper 95% CI`[2])
  ),
  check.names = FALSE
)

knitr::kable(
  case_results,
  col.names = c("Result", "Estimate or test"),
  caption = "Prespecified analysis results for the main cohort"
)

Prespecified analysis results for the main cohort

Result	Estimate or test
Sample and observed endpoints	n=720; events=575; right-censored=145
12-month survival probability	Standard management 43.7%; intervention 64.0%
24-month survival probability	Standard management 21.6%; intervention 15.5%
24-month RMST difference	2.08 months (95% CI 0.92 to 3.23)
Log-rank test	Chi-square=1.78; p=0.182
Adjusted HR, 0 to ≤12 months	0.50 (95% CI 0.40 to 0.63)
Adjusted HR, >12 months	1.83 (95% CI 1.39 to 2.39)

Among 720 simulated participants, 575 first endpoints were observed and 145 participants were right-censored. Through 24 months, estimated mean event-free time was 2.08 months longer with intervention than with standard management (95% CI, 0.92 to 3.23). The effect changed substantially over time, however: the adjusted HR was 0.5 during the first 12 months and 1.83 afterward; the KM curves crossed later, and the PH test also provided evidence against proportionality. A single follow-up-wide HR or log-rank p-value should therefore not be treated as the complete conclusion. These results come from randomized simulated data with a known mechanism and cannot be generalized as evidence of a real intervention effect.

9.1.1 Auditable reporting template

In [target population], follow-up began at [time zero] and continued until [event or observation endpoint]. The primary estimand was [RMST difference/fixed-time risk/other] through $\tau$=[value]. We included [n] people, observed [event count] events, and right-censored [censoring count] people. The estimate for [group] relative to [comparison group] was [effect estimate and 95% CI]; survival probability at [time point] was [value]. Checks of [PH, censoring, functional form, competing risks] showed [result]. Because of [specific design or data limitation], the result should be interpreted as [description, association, conditional prediction, or causal effect under additional assumptions].

10 Common Errors at a Glance

Error	Why it is wrong	Better approach
Describe censoring as “no event”	The outcome after censoring is unknown	Report the observed time known to have been exceeded
Reverse the 0/1 event coding	Curves and model direction become completely reversed	Check frequencies against the source definition before modeling
Show a KM curve without numbers at risk	Reliability near the tail cannot be judged	Report a risk table and confidence intervals
Claim curves are the same because log-rank p>0.05	Failure to reject is not equivalence, and crossing effects may cancel	Report effect estimates and assess heterogeneity over time
Treat the HR as a fixed-time risk ratio	A conditional instantaneous rate differs from a cumulative probability	Also report fixed-time survival or risk
Claim cox.zph p>0.05 proves PH	The test may have low power	Combine the test with plots, numbers at risk, and subject knowledge
Switch automatically to AFT when PH fails	AFT has its own time-scaling assumption	Diagnose each model and choose based on the estimand
Ignore entry	Future enrollees are put into early risk sets	Use Surv(entry, exit, event)
Treat “ever treated later” as a baseline variable	Introduces future information and immortal time bias	Use start–stop data with time-updated treatment
Use 1-KM after censoring competing death	Overestimates the real-world probability of the event of interest	Use the Aalen–Johansen cumulative incidence function
Extrapolate beyond risk-set support	Results are driven by the modeled tail	Mark the observed range and perform sensitivity analyses
Give a causal conclusion from a significant adjusted HR	Design, confounding, and censoring may still cause bias	State the identification conditions and target population

11 Exercises and Answers

1. A participant is lost to follow-up at 8 months without a prior event. What is known about the event time?
2. The two KM curves cross at 15 months. Is it sufficient to report one Cox HR?
3. Why must an RMST analysis prespecify a common cutoff $\tau$?
4. A registry participant enrolls 4 months after diagnosis. How should Surv() be constructed?
5. When death prevents recurrence, why can recurrence probability not be represented by 1-KM after treating death as censored?
6. Medication begins during follow-up. Should “ever medicated” be coded as a baseline variable?

Show exercise answers

1. Only that the true event time exceeds 8 months; it cannot be recorded as never occurring.
2. No. Assess PH and report an estimand suited to the question, such as time-specific effects, fixed-time results, or RMST.
3. RMST is the area under the curve through $\tau$; changing $\tau$ changes the estimand, and both groups must have follow-up support through that cutoff.
4. If time is measured from diagnosis, use Surv(entry, exit, event) with entry=4 months.
5. A person who experiences a competing death can no longer recur. Treating death as ordinary censoring assumes that recurrence remains possible and usually overestimates cumulative incidence.
6. No. That would use future information. Treat medication status as a time-varying covariate and address possible time-varying confounding.

12 Quick Reference

12.1 Core formulas and R entry points

Target	Formula or code	Interpretation
Survival function	$S(t)=P(T>t)$	Probability of remaining event-free through time $t$
Cumulative hazard	$H(t)=-\log S(t)$	Hazard accumulated over time
KM	$\prod_{t_j\le t}(1-d_j/n_j)$	Nonparametric survival estimate under right censoring
Nelson–Aalen	$\sum_{t_j\le t}d_j/n_j$	Nonparametric cumulative hazard estimate
RMST	$\int_0^\tau S(t)dt$	Mean event-free time through $\tau$
KM curve	survfit(Surv(time, event) ~ group, data=d)	Report numbers at risk as well
Log-rank	survdiff(Surv(time, event) ~ group, data=d)	A test, not an effect estimate
Delayed entry	Surv(entry, exit, event)	Enter the risk set only after entry
Time updating	Surv(start, stop, event)	Use the contemporaneous covariate value in each interval
Competing risks	Surv(time, factor_status)	survfit() returns pstate

12.2 Final checklist

• Are the population, time zero, event, time scale, competing events, and observation endpoint explicit?
• Were the primary estimand and time point or RMST cutoff prespecified?
• Were event coding, duplicate records, entry, temporal ordering, and missingness audited?
• Are KM estimates, confidence intervals, numbers at risk, and censoring information shown together?
• Are tests accompanied by effect estimates rather than only p-values?
• Were PH, continuous-variable functional forms, influential observations, and censoring assumptions assessed?
• Are time-varying effects distinguished from time-varying covariates?
• Under competing risks, is the hazard or cumulative-incidence scale aligned with the question?
• Is extrapolation restricted to data-supported ranges and accompanied by sensitivity analyses?
• Does causal or predictive language match the design, validation, and identification assumptions?

Next directions for study

Next topics include flexible parametric survival models, spline-based time effects, adjusted RMST, inverse-probability-of-censoring weighting, semiparametric AFT models, competing-risk regression, joint longitudinal–survival models, recurrent events, frailty, multistate regression, dynamic prediction, and internal and external validation.

sessionInfo()

## R version 4.6.1 (2026-06-24)
## Platform: aarch64-apple-darwin23
## Running under: macOS Tahoe 26.5.1
##
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/4.6/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.6/Resources/lib/libRlapack.dylib;  LAPACK version 3.12.1
##
## locale:
## [1] C.UTF-8/C.UTF-8/C.UTF-8/C/C.UTF-8/C.UTF-8
##
## time zone: America/Edmonton
## tzcode source: internal
##
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base
##
## loaded via a namespace (and not attached):
##  [1] digest_0.6.39   R6_2.6.1        fastmap_1.2.0   Matrix_1.7-5
##  [5] xfun_0.60       lattice_0.22-9  splines_4.6.1   cachem_1.1.0
##  [9] knitr_1.51      htmltools_0.5.9 rmarkdown_2.31  lifecycle_1.0.5
## [13] cli_3.6.6       grid_4.6.1      sass_0.4.10     jquerylib_0.1.4
## [17] compiler_4.6.1  tools_4.6.1     evaluate_1.0.5  bslib_0.12.0
## [21] survival_3.8-6  yaml_2.3.12     rlang_1.3.0     jsonlite_2.0.0

For teaching only: This module introduces the description, estimation, diagnosis, and reporting of survival data. It does not replace a study-specific statistical analysis plan, data audit, or consultation with a biostatistician. All data and numerical results are simulated.