V The V Lab

Accelerated Failure Time and Cox Proportional Hazards Models

A Reproducible R Guide from Censoring and Survival Functions to Model Choice, Diagnosis, and Interpretation

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, and statistics learners who have studied regression fundamentals

Study timeApproximately 120–180 minutes

PrerequisitesConfidence intervals, regression coefficients, and basic R syntax

About the tutorial data This tutorial uses the survival package and a fixed random seed to generate a simulated follow-up cohort. Every record, effect, and result is for teaching only, contains no identifiable personal health information, and must not be treated as clinical or causal evidence for a real population.

How to use this tutorial

This module develops one question throughout: Is an intervention associated with a later occurrence of the outcome? We first study the structure of time-to-event outcomes, then answer the question on two distinct but complementary scales:

• the Cox proportional hazards model (Cox PH) addresses relative instantaneous event rates, with the hazard ratio (HR) as its principal effect measure;
• the accelerated failure time model (AFT) addresses how much event time is stretched or compressed, with the time ratio (TR) as its principal effect measure.

Read the concepts in sequence, run the code, interpret the output, and then open the knowledge checks. Code is shown by default and can be folded with the controls at the top of the page.

Learning objectives

After completing this module, you should be able to:

• define the time origin, event, follow-up time, censoring, and risk set correctly;
• distinguish the survival, hazard, and cumulative hazard functions;
• fit and interpret HRs from Cox PH models and TRs from AFT models;
• explain the difference between the Cox partial likelihood and the full likelihood of a parametric AFT model;
• assess proportional hazards, functional form, influential observations, and AFT distributional assumptions;
• understand the special conversion between HRs and TRs under a Weibull model;
• choose a method according to the research question, the data-supported time range, and model assumptions;
• report absolute survival probabilities, relative effects, uncertainty, diagnostic findings, and limitations.

1 Survival Data: Events, Time, and Censoring

1.1 Why ordinary regression is not enough

An outcome in survival analysis includes not only whether an event occurs, but also when it occurs. When neither of two participants has an observed event, 2 months and 24 months of follow-up provide different amounts of information. An analysis restricted to participants with events would systematically discard the event-free time contributed by censored participants.

Time

From a prespecified origin to the event or last observation.

Event

Defined by reproducible rules with clinical or public health meaning.

Censoring

The event time is known only to exceed an observed time; this does not mean “no event.”

Common examples include time from diagnosis to death, treatment initiation to recurrence, discharge to readmission, or enrollment to discontinuation of a health behavior. “Survival” is a historical term; the event does not have to be death.

1.2 Define time zero, the time scale, and the event first

Every analysis should state at least the following:

Component	Question that must be answered	Definition in this module
Time zero	When does risk begin?	Date intervention or standard management begins
Time scale	Days, months, age, or calendar time?	Months since management began
Event	What counts as an event, and can it occur only once?	First occurrence of the simulated study endpoint
Competing event	Can another event prevent the outcome of interest?	Not simulated, for teaching simplicity
Observation endpoint	When does follow-up stop?	Event, loss to follow-up, or administrative cutoff

Inconsistent definitions of time zero can introduce immortal time bias. For example, assigning only people who survive long enough to receive treatment to the treatment group while starting their follow-up at an earlier diagnosis date creates an artificial period during which death could not have occurred for that group.

1.3 Right censoring and its key assumption

If a participant has not experienced the event by the last observation, we know only that the true event time $T$ exceeds the censoring time $C$. We record:

$\tilde T=\min(T,C)$, together with $\delta=I(T\le C)$

Here, $\delta=1$ indicates an observed event and $\delta=0$ indicates right censoring. Other common observation structures include:

• left censoring: the event occurred before a known time, but its exact time is unknown;
• interval censoring: the event is known only to have occurred between two examinations;
• left truncation/delayed entry: an individual enters the risk set only after surviving and meeting entry conditions; this is not left censoring.

The standard Cox and survreg() examples here focus on right censoring. They generally require the censoring mechanism to carry no additional information about the potential event time after conditioning on model covariates. If participants whose condition is deteriorating rapidly are more likely to be lost to follow-up, and that deterioration is not adequately recorded in the model, treating loss to follow-up as independent censoring can introduce bias.

1.4 Survival, hazard, and cumulative hazard functions

Let $T$ be a continuous event time:

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=-\log S(t)$

The hazard function is the instantaneous event rate immediately after a time among individuals who remain at risk. It is not a probability over a fixed interval and need not lie between 0 and 1. The survival probability $S(t)$ is the probability of remaining event-free from time zero through time $t$.

Hazard is not risk An HR compares conditional instantaneous event rates; a risk ratio compares cumulative event probabilities by a specified time. Their numerical values usually differ even when the HR is constant over time.

1.5 Constructing a survival outcome in R

Surv(time, event) stores right-censoring information using a follow-up-time column and an event-indicator column. This module explicitly uses 1 = event and 0 = censored rather than relying on automatic interpretation of character or factor states.

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

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

data_preview <- transform(
  head(survival_data[, c(
    "participant_id", "time_months", "event",
    "treatment", "age", "severity", "biomarker"
  )]),
  event = ifelse(event == 1, "Event", "Censored"),
  treatment = ifelse(treatment == "Intervention", "Intervention", "Standard management"),
  severity = c("Mild" = "Mild", "Moderate" = "Moderate", "Severe" = "Severe")[
    as.character(severity)
  ]
)

knitr::kable(
  data_preview,
  col.names = c(
    "Participant", "Observed time (months)", "Observed endpoint",
    "Management strategy", "Age", "Baseline severity", "Standardized biomarker"
  ),
  caption = "First six rows of the simulated follow-up data"
)

First six rows of the simulated follow-up data

Participant	Observed time (months)	Observed endpoint	Management strategy	Age	Baseline severity	Standardized biomarker
S001	2.17	Event	Standard management	72	Mild	-1.54
S002	3.05	Event	Intervention	47	Severe	0.61
S003	11.72	Censored	Standard management	52	Mild	0.26
S004	28.41	Censored	Intervention	37	Mild	0.51
S005	13.97	Event	Standard management	43	Mild	-0.55
S006	3.18	Event	Standard management	64	Mild	-0.70

followup_summary <- data.frame(
  Sample_size = nrow(survival_data),
  Events = sum(survival_data$event),
  Censored = sum(survival_data$event == 0),
  Censoring_percentage = pct(mean(survival_data$event == 0)),
  Median_observed_months = median(survival_data$time_months)
)
knitr::kable(
  followup_summary,
  digits = 2,
  col.names = c(
    "Sample size", "Events", "Censored",
    "Censoring percentage", "Median observed time (months)"
  ),
  caption = "Overview of follow-up completeness"
)

Overview of follow-up completeness

Sample size	Events	Censored	Censoring percentage	Median observed time (months)
650	416	234	36.0%	11.1

Longer observed time does not necessarily imply longer true event time because censored event times are unknown. A data description should report the sample size, event count, censoring count, time range, and number at risk in key groups.

1.6 Kaplan–Meier curves: inspect the data before modeling

The Kaplan–Meier (KM) estimator updates survival probability at each event time using the risk set:

$\widehat S(t)=\prod_{t_j\le t}\left(1-\dfrac{d_j}{n_j}\right)$

Here, $d_j$ is the number of events at time $t_j$, and $n_j$ is the number still at risk immediately before that time. Censoring is not counted as an event, but a participant leaves subsequent risk sets after being censored.

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

plot(
  km_fit,
  col = c(palette_surv["orange"], palette_surv["teal"]),
  lwd = 2.3,
  mark.time = TRUE,
  conf.int = FALSE,
  xlab = "Months since management began",
  ylab = "Estimated probability of remaining event-free",
  xlim = c(0, 32),
  ylim = c(0, 1),
  las = 1
)
legend(
  "topright",
  legend = c("Standard management", "Intervention"),
  col = c(palette_surv["orange"], palette_surv["teal"]),
  lwd = 2.3,
  bty = "n"
)

Kaplan–Meier survival curves by management strategy. Short vertical marks indicate censoring; fewer participants remain at risk late in follow-up, so uncertainty is usually greater.

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

km_table <- data.frame(
  Management = ifelse(
    km_selected$strata == "treatment=Intervention",
    "Intervention",
    "Standard management"
  ),
  Time_months = km_selected$time,
  Survival_probability = km_selected$surv,
  CI_lower = km_selected$lower,
  CI_upper = km_selected$upper
)

knitr::kable(
  km_table,
  digits = 3,
  col.names = c(
    "Management strategy", "Time (months)", "Survival probability",
    "95% CI lower", "95% CI upper"
  ),
  caption = "KM survival probabilities and 95% confidence intervals at prespecified times"
)

KM survival probabilities and 95% confidence intervals at prespecified times

Management strategy	Time (months)	Survival probability	95% CI lower	95% CI upper
Standard management	6	0.804	0.758	0.842
Standard management	12	0.529	0.472	0.582
Standard management	18	0.312	0.257	0.370
Intervention	6	0.847	0.802	0.883
Intervention	12	0.619	0.561	0.672
Intervention	18	0.440	0.378	0.500

km_risk_summary <- summary(
  km_fit,
  times = c(0, 6, 12, 18, 24)
)
km_risk_long <- data.frame(
  Management = ifelse(
    km_risk_summary$strata == "treatment=Intervention",
    "Intervention",
    "Standard management"
  ),
  Time_months = km_risk_summary$time,
  Number_at_risk = km_risk_summary$n.risk
)
km_risk_table <- xtabs(
  Number_at_risk ~ Management + Time_months,
  data = km_risk_long
)
names(dimnames(km_risk_table)) <- c("Management strategy", "Time (months)")

knitr::kable(
  km_risk_table,
  caption = "Numbers remaining at risk immediately before selected KM time points"
)

Numbers remaining at risk immediately before selected KM time points

	0	6	12	18	24
Intervention	308	261	159	87	32
Standard management	342	275	139	61	15

KM curves are unadjusted descriptions. Differences between curves may reflect management strategy, age, severity, or other baseline differences. Crossing curves, systematic changes in separation over time, or sparse late follow-up can all affect subsequent model choice and interpretation.

KM estimation also relies on censoring carrying no additional prognostic information within groups. If a group’s curve never falls to 0.50, its KM median event time is not reached; the last observed time must not be reported as the median.

Check your understanding: does a censored participant have “no event”?

A participant is lost to follow-up after 10 months without a prior outcome. Can the event time be recorded as infinity or the outcome as never occurring?

Answer: No. We know only that the true event time exceeds 10 months. Right-censoring methods retain the information from those 10 months without assuming that the event can never occur later.

2 The Cox Proportional Hazards Model

2.1 What question does the model answer?

The Cox PH model links covariates to the conditional hazard function:

$h(t\mid X)=h_0(t)\exp(\beta_1X_1+\cdots+\beta_pX_p)$

• $h_0(t)$ is the unknown baseline hazard function and may vary freely over time;
• $\exp(\beta_j)$ is the HR associated with a one-unit increase in $X_j$, holding other covariates fixed;
• “proportional hazards” means that the ratio of two hazard functions does not change over time.

For a binary intervention variable:

$HR=\dfrac{h(t\mid \text{intervention})}{h(t\mid \text{standard management})}=\exp(\beta_{\text{intervention}})$

An HR below 1 indicates a lower instantaneous event rate in the intervention group among comparable participants who remain at risk at each time. It does not directly quantify the reduction in event probability and is not a multiplier for survival time.

2.2 How partial likelihood avoids specifying the baseline hazard

At every observed event time, the Cox model compares the covariates of the participant who has the event with those of everyone in the risk set at that time. Partial likelihood estimates $\beta$ primarily from relative information about who experiences the event next, without first imposing a parametric distribution on $h_0(t)$.

Ignoring tied events for the moment, the partial likelihood is:

$$L_p(\beta)= \prod_{i:\delta_i=1} \frac{\exp(X_i^\mathsf{T}\beta)} {\sum_{j\in R(t_i)}\exp(X_j^\mathsf{T}\beta)},$$

Here, $R(t_i)$ is the set of participants still at risk immediately before event time $t_i$. The numerator corresponds to the participant who actually has the event, and the denominator sums the relative hazards of everyone who could have had it at that time.

Two consequences follow:

1. the Cox model is flexible about the shape of the baseline hazard;
2. its coefficients are estimated by partial likelihood, whereas a parametric AFT model’s AIC is based on a full event-time likelihood, so the two AIC values should not be compared directly to decide which model is better.

When multiple events are recorded at the same time, tied events occur. Times in this module are rounded to two decimal places, so we explicitly use the commonly applied Efron approximation. If time is inherently recorded in coarse discrete intervals and ties are very frequent, reconsider both the time representation and a discrete-time model.

2.3 Fitting an adjusted Cox model

Age is divided by 10 and centered at 55 years so that its HR compares a 10-year increase and the reference age used for prediction is easy to interpret.

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

cox_summary <- summary(cox_fit)
cox_terms <- rownames(cox_summary$coefficients)
cox_labels <- c(
  treatmentIntervention = "Intervention vs standard management",
  age10 = "Age (per 10-year increase)",
  severityModerate = "Moderate vs mild severity",
  severitySevere = "Severe vs mild severity",
  biomarker = "Biomarker (per 1-SD increase)"
)

cox_results <- data.frame(
  Term = unname(cox_labels[cox_terms]),
  HR = cox_summary$coefficients[, "exp(coef)"],
  `95% CI lower` = cox_summary$conf.int[, "lower .95"],
  `95% CI upper` = cox_summary$conf.int[, "upper .95"],
  `p-value` = format_p(cox_summary$coefficients[, "Pr(>|z|)"]),
  check.names = FALSE
)

knitr::kable(
  cox_results,
  digits = 3,
  align = c("l", "r", "r", "r", "r"),
  caption = "Adjusted Cox proportional hazards model"
)

Adjusted Cox proportional hazards model

	Term	HR	95% CI lower	95% CI upper	p-value
treatmentIntervention	Intervention vs standard management	0.669	0.549	0.815	<0.001
age10	Age (per 10-year increase)	1.194	1.096	1.301	<0.001
severityModerate	Moderate vs mild severity	1.467	1.182	1.821	<0.001
severitySevere	Severe vs mild severity	2.334	1.791	3.043	<0.001
biomarker	Biomarker (per 1-SD increase)	0.998	0.901	1.104	0.962

cox_treatment_hr <- unname(
  cox_summary$coefficients["treatmentIntervention", "exp(coef)"]
)
cox_treatment_ci <- unname(
  cox_summary$conf.int[
    "treatmentIntervention",
    c("lower .95", "upper .95")
  ]
)

Holding age, baseline severity, and biomarker level fixed, the estimated HR for intervention versus standard management is 0.67 (95% confidence interval: 0.55 to 0.82). If proportional hazards and the other model assumptions hold, participants in the intervention group who remain at risk have an instantaneous event rate approximately 66.9% of that in the standard-management group. This is not direct evidence that event risk is reduced by 33.1%, nor is it a causal effect.

2.3.1 Adjusted survival curves

The HR is a relative measure. To support decisions, report absolute survival probabilities at clinically or publicly meaningful times as well. The following example fixes a reference individual at age 55, mild baseline severity, and a biomarker value of 0.

reference_profiles <- data.frame(
  treatment = factor(
    c("Standard", "Intervention"),
    levels = levels(survival_data$treatment)
  ),
  age10 = c(0, 0),
  severity = factor(
    c("Mild", "Mild"),
    levels = levels(survival_data$severity)
  ),
  biomarker = c(0, 0)
)

cox_reference_curves <- survival::survfit(
  cox_fit,
  newdata = reference_profiles
)

plot(
  cox_reference_curves,
  col = c(palette_surv["orange"], palette_surv["teal"]),
  lwd = 2.3,
  conf.int = FALSE,
  xlab = "Months since management began",
  ylab = "Adjusted probability of remaining event-free",
  xlim = c(0, 32),
  ylim = c(0, 1),
  las = 1
)
legend(
  "topright",
  legend = c("Standard management", "Intervention"),
  col = c(palette_surv["orange"], palette_surv["teal"]),
  lwd = 2.3,
  bty = "n"
)

Adjusted survival curves from the Cox model for the same reference covariate profile. The absolute level of each curve depends on the estimated baseline survival function.

These are conditional predictions and represent only the covariate profiles supplied in newdata. If the target is population-average survival, prespecify a standard population and average predictions across its members. Do not automatically describe a “typical individual” curve as a population curve.

2.4 What assumptions does the Cox model require?

Assumption	Meaning	Common assessment
Proportional hazards	The HR for each covariate remains constant over time	Schoenfeld residual plots, cox.zph(), stratified curves, and subject-matter knowledge
Functional form of continuous variables	The form in the linear predictor is correct; for example, age is approximately linear on the log-hazard scale	Martingale residuals, splines, and prespecified nonlinear terms
Conditional independent censoring	After conditioning on covariates, censoring no longer predicts event time	Compare loss-to-follow-up patterns, perform sensitivity analyses, and improve data collection
Dependence structure is handled	Clustering, recurrent events, or multicenter correlation cannot be treated as independent	Robust variance, frailty, multilevel, or recurrent-event methods
Covariate measurement and temporal order are appropriate	Post-baseline information must not be substituted incorrectly for baseline values	Protocol, data-dictionary, and timestamp review

2.4.1 Assessing proportional hazards with Schoenfeld residuals

cox.zph() assesses whether scaled Schoenfeld residuals show a systematic trend over time. A small p-value suggests that an effect may vary with time; a large p-value means only that the current data provide no strong evidence against PH, not that exact proportionality has been proven.

cox_ph_test <- survival::cox.zph(cox_fit, transform = "km")

cox_ph_table <- data.frame(
  Test = c(
    "Management strategy", "Age", "Baseline severity",
    "Biomarker", "Global test"
  ),
  Chi_square = cox_ph_test$table[, "chisq"],
  df = cox_ph_test$table[, "df"],
  `p-value` = format_p(cox_ph_test$table[, "p"]),
  check.names = FALSE
)

knitr::kable(
  cox_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	1.170	1	0.279
age10	Age	0.230	1	0.631
severity	Baseline severity	0.072	2	0.965
biomarker	Biomarker	1.305	1	0.253
GLOBAL	Global test	2.922	5	0.712

old_par <- par(mfrow = c(1, 2), mar = c(4.4, 4.3, 2.6, 1))

plot(
  cox_ph_test,
  var = 1,
  resid = TRUE,
  se = TRUE,
  col = palette_surv["teal"]
)
abline(h = 0, lty = 2, col = palette_surv["gray"])
title("Management strategy")

plot(
  cox_ph_test,
  var = 2,
  resid = TRUE,
  se = TRUE,
  col = palette_surv["blue"]
)
abline(h = 0, lty = 2, col = palette_surv["gray"])
title("Age (per 10 years)")

Scaled Schoenfeld residual diagnostics for management strategy and age. A smooth curve that clearly departs from a horizontal line suggests that the corresponding log HR may change over time.

par(old_par)

The global test p-value in these simulated data is 0.712. Together with the plots, there is no clear evidence against PH, but this conclusion remains limited by the event count, follow-up range, and statistical power. In particular, a model should not be selected mechanically from multiple covariate-specific p-values.

2.4.2 Assessing the functional form of a continuous variable

If the true age effect is curved, forcing a linear age term can distort both the HR and covariate adjustment. One exploratory approach is to fit a model without age and then plot Martingale residuals against age; a smooth curve suggests a potentially useful functional form. A formal analysis should prespecify flexible forms such as restricted cubic splines using subject-matter knowledge and avoid repeated data-driven searching.

cox_without_age <- survival::coxph(
  survival::Surv(time_months, event) ~
    treatment + severity + biomarker,
  data = survival_data,
  ties = "efron"
)

martingale_age <- residuals(cox_without_age, type = "martingale")

plot(
  survival_data$age,
  martingale_age,
  pch = 16,
  cex = 0.55,
  col = grDevices::adjustcolor(palette_surv["navy"], alpha.f = 0.28),
  xlab = "Age (years)",
  ylab = "Martingale residual",
  las = 1
)
abline(h = 0, lty = 2, col = palette_surv["gray"])
lines(
  lowess(survival_data$age, martingale_age, f = 0.65),
  col = palette_surv["vermillion"],
  lwd = 2.4
)

Age plotted against Martingale residuals from a Cox model that omits age. The smooth line explores a possible functional form for age and is not a formal significance test.

2.4.3 Assessing unusual and influential observations

Deviance residuals can help identify observations whose event patterns are fitted poorly, while dfbeta residuals approximate how deleting one observation would affect each coefficient. They are starting points for investigation, not rules for automatic deletion.

cox_dfbeta <- residuals(cox_fit, type = "dfbeta")
colnames(cox_dfbeta) <- names(coef(cox_fit))
max_dfbeta <- apply(abs(cox_dfbeta), 2, max)

influence_table <- data.frame(
  Coefficient = unname(cox_labels[colnames(cox_dfbeta)]),
  Maximum_absolute_dfbeta = max_dfbeta,
  check.names = FALSE
)

knitr::kable(
  influence_table,
  digits = 3,
  col.names = c("Coefficient", "Maximum absolute dfbeta"),
  caption = "Largest absolute dfbeta residual observed for each coefficient"
)

Largest absolute dfbeta residual observed for each coefficient

	Coefficient	Maximum absolute dfbeta
treatmentIntervention	Intervention vs standard management	0.015
age10	Age (per 10-year increase)	0.016
severityModerate	Moderate vs mild severity	0.018
severitySevere	Severe vs mild severity	0.039
biomarker	Biomarker (per 1-SD increase)	0.009

After finding an influential observation, verify its time, event code, covariates, and eligibility, and report sensitivity analyses with and without it when appropriate. Deleting an observation merely because it does not conform to the model creates a new source of bias.

2.5 What to do when proportional hazards does not hold

The response should follow the scientific question, not just a diagnostic p-value:

• if a covariate is included only for adjustment and its HR is not needed, use strata() to allow different baseline hazards across its levels;
• if the effect of the primary exposure varies over time, explicitly include an exposure-by-time interaction and report $HR(t)$;
• if risk at a fixed time or average event-free time is more important, report adjusted survival probabilities or restricted mean survival time (RMST);
• if nonproportional hazards reflect a different time-scale mechanism, consider a suitable AFT model, flexible parametric model, or piecewise model;
• if apparent change occurs only at a few late times with very small risk sets, first assess whether the estimates are sufficiently stable.

The code below demonstrates syntax only. The function inside tt() must be designed in advance for the scientific problem; $\log(t+1)$ is not a universal solution.

cox_time_varying <- survival::coxph(
  survival::Surv(time_months, event) ~
    treatment_num + age10 + severity + biomarker +
    tt(treatment_num),
  data = survival_data,
  ties = "efron",
  tt = function(x, t, ...) x * log(t + 1)
)

A time-varying effect is different from a time-varying covariate “The intervention HR changes over time” describes a time-varying effect. “Blood pressure is repeatedly updated during follow-up” describes a time-varying covariate. The latter generally requires multiple rows per participant in start–stop format with Surv(start, stop, event), and every measurement must be ordered correctly in time.

Check your understanding: how should HR = 0.70 be described?

Can we say, “The intervention reduces 12-month event risk by 30%”?

Answer: Not directly. An HR of 0.70 means that, among participants who remain at risk and have the same model covariates, the intervention group has approximately 0.70 times the instantaneous event rate if PH holds. A 12-month risk difference or risk ratio must be calculated separately from the corresponding survival probabilities.

3 The Accelerated Failure Time Model

3.1 Shifting from event rates to event times

An AFT model directly describes log event time:

$\log(T)=\beta_0+\beta_1X_1+\cdots+\beta_pX_p+\sigma\varepsilon$

An equivalent time-scaling expression is:

$S(t\mid X)=S_0\{t\exp(-X^\mathsf{T}\beta)\}$

For a one-unit increase in $X_j$:

$TR=\exp(\beta_j)$

In the standard AFT structure, the TR is a common multiplier for every conditional event-time quantile:

• $TR>1$: the time scale is stretched, so events tend to occur later;
• $TR<1$: the time scale is compressed, so events tend to occur earlier;
• $TR=1$: the modeled event-time scale is unchanged.

For example, TR = 1.30 means that, holding other covariates fixed, the time required to reach the same event-time quantile is multiplied by 1.30. If the conditional median event time under standard management is 10 months, the corresponding conditional median under intervention is 13 months. This does not mean that risk is reduced by 30%.

A parametric AFT model constructs a full likelihood from the event density for observed events and the information that censored participants remained event-free through their observation times:

$$L(\theta)= \prod_i f(\tilde T_i\mid X_i;\theta)^{\delta_i} S(\tilde T_i\mid X_i;\theta)^{1-\delta_i}.$$

Thus, an event record contributes the density $f$, while a right-censored record contributes the survival probability $S$. The chosen distribution determines both contributions as well as tail behavior.

3.2 A parametric AFT model requires an event-time distribution

The parametric AFT models in this tutorial use a complete event-time distribution. survreg() specifies an error distribution for $\log(T)$ in location–scale form. Semiparametric AFT methods based on rank estimation also exist and do not require the same full parametric distribution. Common survreg() choices include:

survreg() distribution	Distribution of $T$	Typical hazard shape	Also satisfies PH?
exponential	Exponential	Constant hazard	Yes
weibull	Weibull	Monotonically increasing or decreasing	Yes
lognormal	Log-normal	Often increases and then decreases	Usually no
loglogistic	Log-logistic	Can increase and then decrease; heavier tail	Usually no

Distribution choice should integrate the underlying mechanism, KM curves, residuals, calibration, AIC, and the time range over which predictions are needed. A lower AIC indicates only a better relative tradeoff between fit and complexity among the candidate models on the same data. It neither proves that a distribution is true nor protects long-term extrapolation.

3.3 Fitting a Weibull AFT model

aft_weibull <- survival::survreg(
  survival::Surv(time_months, event) ~
    treatment + age10 + severity + biomarker,
  data = survival_data,
  dist = "weibull"
)

aft_summary <- summary(aft_weibull)
aft_terms <- names(coef(aft_weibull))[-1]
aft_labels <- cox_labels

aft_results <- data.frame(
  Term = unname(aft_labels[aft_terms]),
  `Log-time coefficient` = coef(aft_weibull)[aft_terms],
  TR = exp(coef(aft_weibull)[aft_terms]),
  `95% CI lower` = exp(
    coef(aft_weibull)[aft_terms] -
      1.96 * aft_summary$table[aft_terms, "Std. Error"]
  ),
  `95% CI upper` = exp(
    coef(aft_weibull)[aft_terms] +
      1.96 * aft_summary$table[aft_terms, "Std. Error"]
  ),
  `p-value` = format_p(aft_summary$table[aft_terms, "p"]),
  check.names = FALSE
)

knitr::kable(
  aft_results,
  digits = 3,
  align = c("l", "r", "r", "r", "r", "r"),
  caption = "Adjusted Weibull AFT model"
)

Adjusted Weibull AFT model

	Term	Log-time coefficient	TR	95% CI lower	95% CI upper	p-value
treatmentIntervention	Intervention vs standard management	0.267	1.306	1.146	1.487	<0.001
age10	Age (per 10-year increase)	-0.118	0.889	0.840	0.941	<0.001
severityModerate	Moderate vs mild severity	-0.253	0.777	0.673	0.897	<0.001
severitySevere	Severe vs mild severity	-0.556	0.573	0.481	0.683	<0.001
biomarker	Biomarker (per 1-SD increase)	0.004	1.004	0.938	1.074	0.912

aft_treatment_beta <- unname(
  coef(aft_weibull)["treatmentIntervention"]
)
aft_treatment_se <- unname(
  aft_summary$table["treatmentIntervention", "Std. Error"]
)
aft_treatment_tr <- exp(aft_treatment_beta)
aft_treatment_ci <- exp(
  aft_treatment_beta + c(-1, 1) * 1.96 * aft_treatment_se
)
aft_weibull_shape <- 1 / aft_weibull$scale

Holding age, baseline severity, and biomarker level fixed, the estimated TR for intervention versus standard management is 1.31 (95% confidence interval: 1.15 to 1.49). The model therefore estimates every conditional event-time quantile to be 1.31 times that under standard management, an extension of approximately 30.6%. This statement depends on the Weibull distribution, a constant acceleration factor, correct covariate forms, independent censoring, and the other model assumptions.

3.3.1 The survreg() scale is not the Weibull shape

This is one of the most common parameterization pitfalls. Under R’s survreg(dist="weibull") parameterization:

$\text{Weibull shape}=\dfrac{1}{\texttt{survreg scale}}$

The survreg scale in this example is 0.663, corresponding to a Weibull shape of approximately 1.508. A shape above 1 indicates that the baseline hazard increases over time, a shape of 1 gives the exponential distribution, and a shape below 1 indicates a decreasing baseline hazard.

Let $\eta=X^\mathsf{T}\beta$. R’s Weibull AFT parameterization corresponds to:

$S(t\mid X)=\exp\left[-\left\{t/\exp(\eta)\right\}^{\kappa}\right]$, where $\kappa=1/\sigma$

Do not copy parameters named scale or shape directly between software packages. Before using another function, determine whether it uses a Weibull proportional-hazards parameterization, an AFT parameterization, or a different parameterization.

3.4 Predicting event-time quantiles with an AFT model

The next example predicts the 25th, 50th, and 75th percentiles of event time for the same two reference covariate profiles used for the Cox curves. Here, the 25th percentile is the time by which 25% of conditional event times are expected to have occurred.

aft_quantile_prediction <- predict(
  aft_weibull,
  newdata = reference_profiles,
  type = "quantile",
  p = c(0.25, 0.50, 0.75),
  se.fit = TRUE
)

aft_quantiles <- aft_quantile_prediction$fit
aft_quantile_se <- aft_quantile_prediction$se.fit
quantile_labels <- c("25th percentile", "Median", "75th percentile")

aft_quantile_table <- data.frame(
  Management = rep(c("Standard management", "Intervention"), times = 3),
  Event_time_quantile = rep(quantile_labels, each = 2),
  `Estimated time (months)` = as.vector(aft_quantiles),
  `95% CI lower` = pmax(
    0,
    as.vector(aft_quantiles - 1.96 * aft_quantile_se)
  ),
  `95% CI upper` = as.vector(
    aft_quantiles + 1.96 * aft_quantile_se
  ),
  check.names = FALSE
)

knitr::kable(
  aft_quantile_table,
  digits = 2,
  col.names = c(
    "Management strategy", "Event-time quantile", "Estimated time (months)",
    "95% CI lower", "95% CI upper"
  ),
  caption = "Conditional event-time quantiles and approximate 95% confidence intervals from the Weibull AFT model"
)

Conditional event-time quantiles and approximate 95% confidence intervals from the Weibull AFT model

Management strategy	Event-time quantile	Estimated time (months)	95% CI lower	95% CI upper
Standard management	25th percentile	8.29	7.27	9.32
Intervention	25th percentile	10.83	9.42	12.24
Standard management	Median	14.86	13.20	16.52
Intervention	Median	19.40	17.02	21.78
Standard management	75th percentile	23.53	20.80	26.26
Intervention	75th percentile	30.72	26.75	34.68

These intervals are Wald approximations based on the fitted model’s asymptotic covariance and do not include uncertainty from distribution selection. Because censoring increases late in follow-up, higher quantiles may extend beyond the best-supported range of the data. Numerical precision in the prediction table should not exceed what the data and model can support.

3.5 Comparing candidate AFT distributions

AIC values are comparable only when every candidate model uses the same participants, outcome definition, and covariates. The same right-censoring likelihood is retained here as well.

aft_exponential <- update(aft_weibull, dist = "exponential")
aft_lognormal <- update(aft_weibull, dist = "lognormal")
aft_loglogistic <- update(aft_weibull, dist = "loglogistic")

aft_aic <- AIC(
  aft_exponential,
  aft_weibull,
  aft_lognormal,
  aft_loglogistic
)

aft_aic_table <- data.frame(
  Distribution = c("Exponential", "Weibull", "Log-normal", "Log-logistic"),
  Parameters = aft_aic[, "df"],
  AIC = aft_aic[, "AIC"],
  `Difference from minimum AIC` = aft_aic[, "AIC"] - min(aft_aic[, "AIC"]),
  check.names = FALSE
)

knitr::kable(
  aft_aic_table[order(aft_aic_table$AIC), ],
  digits = 1,
  caption = "Candidate parametric AFT models using the same data and covariates"
)

Candidate parametric AFT models using the same data and covariates

	Distribution	Parameters	AIC	Difference from minimum AIC
2	Weibull	7	3182	0.0
4	Log-logistic	7	3206	24.2
3	Log-normal	7	3240	57.7
1	Exponential	6	3264	81.9

The Weibull model has the lowest AIC, consistent with the data-generating mechanism, but a real analysis does not reveal the “correct answer.” Distributional diagnostics should be evaluated alongside the clinical process and predictive calibration. If several models agree during observed follow-up but diverge substantially when extrapolated, treat distribution choice as a sensitivity analysis rather than reporting only the most favorable prediction.

3.6 Assessing overall fit with Cox–Snell residuals

For a Weibull AFT model, the Cox–Snell residual for an individual at observed time $\tilde T_i$ is the fitted cumulative hazard:

$r_i=\widehat H_i(\tilde T_i)=\exp\left\{\dfrac{\log(\tilde T_i)-\widehat\mu_i}{\widehat\sigma}\right\}$

If the model is well calibrated overall, the cumulative hazard of these residuals should approximately follow $y=x$. Because the residuals are still censored, their cumulative hazard must be estimated with survival methods.

aft_linear_predictor <- predict(aft_weibull, type = "lp")
cox_snell_residual <- exp(
  (log(survival_data$time_months) - aft_linear_predictor) /
    aft_weibull$scale
)

cox_snell_fit <- survival::survfit(
  survival::Surv(cox_snell_residual, survival_data$event) ~ 1
)
cox_snell_cumhaz <- -log(cox_snell_fit$surv)
diagnostic_limit <- unname(
  quantile(cox_snell_fit$time, probs = 0.95)
)

plot(
  cox_snell_fit$time,
  cox_snell_cumhaz,
  type = "s",
  lwd = 2.3,
  col = palette_surv["teal"],
  xlim = c(0, diagnostic_limit),
  ylim = c(0, diagnostic_limit),
  xlab = "Cox–Snell residual",
  ylab = "Estimated cumulative hazard of residuals",
  las = 1
)
abline(
  a = 0,
  b = 1,
  lty = 2,
  lwd = 2,
  col = palette_surv["vermillion"]
)
legend(
  "topleft",
  legend = c("Estimated curve", "Ideal reference line"),
  col = c(palette_surv["teal"], palette_surv["vermillion"]),
  lty = c(1, 2),
  lwd = 2.2,
  bty = "n"
)

Cox–Snell residual diagnostic for the Weibull AFT model. An estimated cumulative hazard close to the 45-degree line is consistent with adequate overall distributional fit; late deviations are often affected by dwindling risk sets.

A Cox–Snell plot is a global check and may conceal a misspecified functional form for one covariate. Also compare observed and predicted survival within groups, inspect deviance or response residuals, assess continuous-variable forms, and internally validate key predictions. A deviation caused by a few observations in the tail should not be interpreted in the same way as systematic deviation across the main follow-up range.

3.7 What assumptions does an AFT model require?

• Correct distributional form: the selected Weibull, log-normal, or other distribution must reasonably describe conditional event times;
• Constant acceleration factor: a covariate stretches or compresses the entire conditional time distribution by the same factor;
• Correct covariate functional form: the form of a continuous variable on the log-time scale must be reasonable;
• Conditional independent censoring: after conditioning on model information, censoring carries no further information about the potential event time;
• Correct observation structure: clustering, recurrent events, competing risks, and delayed entry require suitable methods;
• Defensible extrapolation: predictions beyond observed follow-up are driven primarily by assumptions about the distributional tail.

Statistical significance cannot choose the effect scale Do not report an HR merely because the Cox p-value is smaller, and do not ignore distributional fit merely because an AFT TR is more intuitive. The research question, target estimand, and assumptions should determine the primary model before the results are known.

Check your understanding: what does TR = 1.25 mean?

How should a TR of 1.25 for a binary exposure be interpreted in a standard AFT model?

Answer: Holding other covariates fixed and assuming the AFT model holds, every quantile of the exposed group’s conditional event-time distribution is 1.25 times the corresponding reference-group quantile, so the time scale is extended by approximately 25%. This is not an HR of 0.75 and does not directly mean a 25% reduction in risk at a fixed time.

4 Relationships, Differences, and Choice Between AFT and Cox PH Models

4.1 The two models do not negate each other

Feature	Cox PH	Parametric AFT
Principal question	How much does the instantaneous event rate change relatively?	How much is event time stretched or compressed?
Principal effect measure	HR $=\exp(\beta_{\text{PH}})$	TR $=\exp(\beta_{\text{AFT}})$
Baseline structure	No parametric shape specified for $h_0(t)$	A complete time distribution must be selected
Likelihood information	Partial likelihood estimates relative-hazard coefficients	Full likelihood estimates location and scale jointly
Core effect assumption	The HR remains constant over time	The time distribution is scaled by a constant factor
Absolute prediction	Estimated baseline survival supports prediction within the observed range	Quantiles and survival probabilities can be predicted directly
Extrapolation	Generally unsuitable beyond the last event time	Computable, but highly dependent on the tail distribution
Interpretive advantage	Common in clinical literature; flexible baseline hazard	“How much earlier or later?” is often easier to communicate

Cox PH and AFT models do not reduce to a simple contrast in which semiparametric models are always robust and parametric models are always risky. A Cox model still requires proportional hazards, functional forms, and the censoring structure to be appropriate. When its distribution is reasonable, an AFT model can be efficient and provide a direct time-scale interpretation.

4.2 Weibull is the intersection of the two model families

The Weibull distribution belongs to both the PH and AFT families. Let the AFT survreg scale be $\sigma$ and the Weibull shape be $\kappa=1/\sigma$. For the same covariate, the coefficients then satisfy:

$\beta_{\text{PH}}=-\kappa\beta_{\text{AFT}}$

$HR=\exp(-\kappa\beta_{\text{AFT}})=TR^{-\kappa}$

aft_implied_hr <- exp(
  -aft_weibull_shape * aft_treatment_beta
)

bridge_table <- data.frame(
  Source = c("Direct Cox PH estimate", "Converted from Weibull AFT"),
  Intervention_vs_standard_HR = c(
    cox_treatment_hr,
    aft_implied_hr
  )
)

knitr::kable(
  bridge_table,
  digits = 3,
  col.names = c("Source", "HR: intervention vs standard management"),
  caption = "Comparison of HRs when a Weibull model satisfies both PH and AFT"
)

Comparison of HRs when a Weibull model satisfies both PH and AFT

Source	HR: intervention vs standard management
Direct Cox PH estimate	0.669
Converted from Weibull AFT	0.669

The two estimates are close in this example because the data were simulated from a Weibull mechanism; this should not be expected for every dataset. Log-normal and log-logistic AFT models generally do not produce a constant HR and cannot use this conversion. Even when the same data approximately satisfy both structures, HRs and TRs still answer different questions.

4.3 A practical selection process

1. Write down the target estimand first. Does the decision-maker need risk at a fixed time, an HR, a TR, median event time, or an RMST difference?
2. Plot KM curves and numbers at risk. Examine curve shapes, crossings, censoring, and the time range supported by the data.
3. Specify the primary model from the question. Do not choose an interpretation scale after inspecting p-values.
4. Assess core assumptions. Assess PH for Cox; assess time scaling and distributional calibration for AFT; assess functional forms and censoring for both.
5. Report absolute quantities. Even when the primary effect is an HR or TR, give a survival probability at a prespecified time or an event-time quantile.
6. Perform sensitivity analyses. Compare reasonable distributions, time-varying effects, nonlinear forms, and potentially informative censoring.
7. Limit extrapolation and causal language. Model fit does not automatically make the study design support causal interpretation.

4.3.1 When Cox PH may be preferable

• the research question and field conventions explicitly focus on conditional instantaneous event rates;
• no parametric distribution is desired for the baseline hazard;
• the primary focus is a relative effect during observed follow-up;
• proportional hazards is approximately reasonable over the important time range.

4.3.2 When AFT may be preferable

• “How much is the event delayed or advanced?” is the primary scientific question;
• a reasonable distribution agrees with the mechanism, plots, and diagnostics;
• conditional event-time quantiles or well-validated parametric predictions are needed;
• PH is inappropriate and constant time scaling is better supported.

4.3.3 When neither may be the best primary choice

• curves cross clearly and the effect direction changes over time;
• competing events are common and the target is cumulative incidence;
• recurrent events, multistate transitions, or complex clustering are present;
• the primary target is average event-free time within a fixed window, for which RMST may be more direct;
• censoring depends strongly on unobserved health status, making a simple independent-censoring model indefensible.

Check your understanding: can Cox AIC and Weibull AFT AIC be compared?

Both models were fitted to the same data. Can the model with the smaller AIC be selected?

Answer: Not directly. Standard Cox coefficients come from a partial likelihood, whereas a parametric AFT model’s AIC comes from a full event-time likelihood; their likelihood bases differ. AIC can be compared across full-likelihood parametric models fitted to the same data and outcome, but diagnostics and the research question still matter.

5 Common Complexities and Pitfalls

5.1 Delayed entry, clustering, and recurrent events

• Delayed entry: in a Cox model, a participant contributes to the risk set only after entering the study; use coxph(Surv(entry, exit, event) ~ ...). Ignoring entry conditions can cause selection bias. Standard survreg() does not support this start–stop input, so a parametric AFT analysis with left truncation requires another implementation with the appropriate likelihood;
• center or household clustering: depending on the target, use cluster-robust standard errors, a shared frailty, or a multilevel model;
• recurrent events: when the same person can experience the outcome multiple times, event gaps, the overall event process, and terminal events require explicit definitions;
• multistate processes: transitions such as “healthy → hospitalized → dead” cannot be represented fully by a single endpoint.

Changing only the standard errors cannot repair an incorrect risk set, time origin, or target estimand.

5.2 Competing risks

If death prevents recurrence, death is a competing event for recurrence. When a competing death is treated as ordinary censoring:

• a cause-specific Cox model can estimate the cause-specific hazard;
• however, $1-\widehat S(t)$ is generally not the cumulative incidence of the event of interest;
• if the target is real-world cumulative incidence, use a cumulative incidence function and specify whether the effect scale is cause-specific or subdistribution-based.

Treating a competing event as censoring is not a computational error, but it changes the estimand and must align with the research question.

5.3 Missing data and measurement time

A complete-case analysis can be unbiased only under missingness conditions appropriate to the target analysis. A multiple-imputation model should include outcome information, follow-up information, important auxiliary variables, and variables related to missingness. Treating a covariate measured after the event as a baseline adjustment variable can control a mediator, introduce collider bias, or violate temporal order.

5.4 Causal interpretation requires additional conditions

An adjusted HR or TR may still be affected by unmeasured confounding, selection bias, measurement error, and model specification. If the target is a causal effect, specify:

• the intervention, comparator, and target population;
• the start of follow-up, grace period, and treatment strategy;
• baseline confounders and their causal justification;
• handling of loss to follow-up and treatment switching;
• whether the target effect is conditional or marginal;
• identification conditions such as consistency, exchangeability, and positivity.

Adjustment variables should be selected from the research question, temporal order, and causal structure, not screened by univariable p-values. Continuous variables should not be dichotomized arbitrarily merely to create “significant groups.”

Model accuracy is not decision fairness High-risk predictions may reflect access to care, opportunities for diagnosis, or structural inequities rather than biological risk alone. When reporting a model, examine variable meanings, censoring and calibration across groups, potential harms, and how results will be used.

6 Guided Mini Case Study

6.1 Research question and analysis plan

Research question: In the simulated cohort, how is the intervention associated with time to the first occurrence of the study endpoint?

Before examining the results, prespecify:

• primary Cox estimand: the intervention HR adjusted for age, severity, and biomarker level;
• primary AFT estimand: the Weibull TR with the same adjustment variables;
• absolute outcome: the 12-month probability of remaining event-free for a 55-year-old reference individual with mild severity and a biomarker value of 0;
• diagnostics: PH assessment, AFT distribution comparison, and a Cox–Snell plot;
• interpretive limits: intervention is not randomized and all data are simulated.

6.2 Step 1: Bring descriptive, relative, and absolute results together

km_median <- summary(km_fit)$table

case_relative <- data.frame(
  Model = c("Cox PH", "Weibull AFT"),
  Effect_measure = c("HR", "TR"),
  Estimate = c(cox_treatment_hr, aft_treatment_tr),
  `95% CI lower` = c(cox_treatment_ci[1], aft_treatment_ci[1]),
  `95% CI upper` = c(cox_treatment_ci[2], aft_treatment_ci[2]),
  check.names = FALSE
)

case_medians <- data.frame(
  Management = c("Standard management", "Intervention"),
  `Unadjusted KM median time (months)` = km_median[, "median"],
  `95% CI lower` = km_median[, "0.95LCL"],
  `95% CI upper` = km_median[, "0.95UCL"],
  check.names = FALSE
)

knitr::kable(
  case_relative,
  digits = 3,
  col.names = c("Model", "Effect measure", "Estimate", "95% CI lower", "95% CI upper"),
  caption = "Two adjusted relative-effect scales"
)

Two adjusted relative-effect scales

Model	Effect measure	Estimate	95% CI lower	95% CI upper
Cox PH	HR	0.669	0.549	0.815
Weibull AFT	TR	1.306	1.146	1.487

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

Unadjusted KM median event time by management strategy

	Management strategy	Unadjusted KM median time (months)	95% CI lower	95% CI upper
treatment=Standard	Standard management	13.0	11.2	14.1
treatment=Intervention	Intervention	15.7	14.3	18.2

An unadjusted KM median and an adjusted AFT conditional median are not the same estimand. The former describes the observed covariate mixture within each group; the latter fixes or conditions on covariates in the model.

6.3 Step 2: Calculate absolute 12-month survival probabilities

cox_12_summary <- summary(
  cox_reference_curves,
  times = 12
)
cox_survival_12 <- as.numeric(cox_12_summary$surv)
cox_survival_12_lower <- as.numeric(cox_12_summary$lower)
cox_survival_12_upper <- as.numeric(cox_12_summary$upper)

aft_reference_lp <- predict(
  aft_weibull,
  newdata = reference_profiles,
  type = "lp"
)
aft_survival_12 <- exp(
  -exp(
    (log(12) - aft_reference_lp) /
      aft_weibull$scale
  )
)

# Use simulation to propagate the large-sample joint covariance of survreg parameters.
set.seed(20260814)
n_parameter_draws <- 4000
aft_parameter_mean <- c(
  coef(aft_weibull),
  log(aft_weibull$scale)
)
standard_normal_draws <- matrix(
  rnorm(n_parameter_draws * length(aft_parameter_mean)),
  nrow = n_parameter_draws
)
aft_parameter_draws <- sweep(
  standard_normal_draws %*% chol(aft_weibull$var),
  MARGIN = 2,
  STATS = aft_parameter_mean,
  FUN = "+"
)

n_aft_coefficients <- length(coef(aft_weibull))
aft_beta_draws <- aft_parameter_draws[
  , seq_len(n_aft_coefficients), drop = FALSE
]
aft_sigma_draws <- exp(
  aft_parameter_draws[, n_aft_coefficients + 1]
)
aft_reference_matrix <- model.matrix(
  delete.response(terms(aft_weibull)),
  data = reference_profiles
)
aft_lp_draws <- aft_beta_draws %*% t(aft_reference_matrix)
aft_standardized_time <- sweep(
  log(12) - aft_lp_draws,
  MARGIN = 1,
  STATS = aft_sigma_draws,
  FUN = "/"
)
aft_survival_draws <- exp(-exp(aft_standardized_time))
aft_survival_12_ci <- apply(
  aft_survival_draws,
  MARGIN = 2,
  FUN = quantile,
  probs = c(0.025, 0.975)
)

all_survival_12 <- c(cox_survival_12, aft_survival_12)
all_survival_12_lower <- c(
  cox_survival_12_lower,
  aft_survival_12_ci[1, ]
)
all_survival_12_upper <- c(
  cox_survival_12_upper,
  aft_survival_12_ci[2, ]
)

absolute_results <- data.frame(
  Model = rep(c("Cox PH", "Weibull AFT"), each = 2),
  Management = rep(c("Standard management", "Intervention"), times = 2),
  `12-month survival probability` = all_survival_12,
  `95% CI lower` = all_survival_12_lower,
  `95% CI upper` = all_survival_12_upper,
  `12-month cumulative event probability` = 1 - all_survival_12,
  check.names = FALSE
)

knitr::kable(
  absolute_results,
  digits = 3,
  col.names = c(
    "Model", "Management strategy", "12-month survival probability",
    "95% CI lower", "95% CI upper", "12-month cumulative event probability"
  ),
  caption = "Conditional 12-month model predictions and 95% confidence intervals for the reference profile"
)

Conditional 12-month model predictions and 95% confidence intervals for the reference profile

Model	Management strategy	12-month survival probability	95% CI lower	95% CI upper	12-month cumulative event probability
Cox PH	Standard management	0.611	0.557	0.670	0.389
Cox PH	Intervention	0.719	0.671	0.771	0.281
Weibull AFT	Standard management	0.605	0.552	0.657	0.395
Weibull AFT	Intervention	0.715	0.666	0.759	0.285

For this reference individual, the Cox model estimates that 12-month survival probability increases from 61.1% under standard management to 71.9% under intervention. This is a conditional model prediction, not a marginal risk from a randomized trial, and it does not imply that every participant would receive the same benefit.

The Cox intervals reflect asymptotic model uncertainty. The AFT intervals propagate the joint covariance of coefficients and log scale through 4,000 draws from a normal parameter approximation. Both condition on the selected model form and omit uncertainty from distribution selection, unmeasured confounding, and prediction error in new data.

6.4 Step 3: Compare model predictions within the observed range

prediction_times <- seq(0.1, 30, by = 0.1)
aft_curve_matrix <- sapply(
  aft_reference_lp,
  function(mu) {
    exp(
      -exp(
        (log(prediction_times) - mu) /
          aft_weibull$scale
      )
    )
  }
)

plot(
  cox_reference_curves,
  col = c(palette_surv["orange"], palette_surv["teal"]),
  lwd = 2.4,
  lty = 1,
  conf.int = FALSE,
  xlab = "Months since management began",
  ylab = "Predicted probability of remaining event-free",
  xlim = c(0, 30),
  ylim = c(0, 1),
  las = 1
)
lines(
  prediction_times,
  aft_curve_matrix[, 1],
  col = palette_surv["orange"],
  lwd = 2.4,
  lty = 2
)
lines(
  prediction_times,
  aft_curve_matrix[, 2],
  col = palette_surv["teal"],
  lwd = 2.4,
  lty = 2
)
legend(
  "topright",
  legend = c(
    "Standard management: Cox", "Intervention: Cox",
    "Standard management: AFT", "Intervention: AFT"
  ),
  col = c(
    palette_surv["orange"], palette_surv["teal"],
    palette_surv["orange"], palette_surv["teal"]
  ),
  lty = c(1, 1, 2, 2),
  lwd = 2.3,
  bty = "n"
)

Survival predictions from Cox PH and Weibull AFT models for the same reference individual. Color identifies management strategy and line type identifies model; agreement over the observed range does not guarantee agreement in long-term extrapolation.

The two predictions are close in this example, as expected under the Weibull data-generating mechanism. In a real study, compare calibration only where enough participants remain at risk. Do not allow curves to extend automatically into poorly supported long-term follow-up and then interpret them with certainty.

6.5 Step 4: Write an auditable result

Among 650 simulated participants, 416 first endpoints were observed and 234 participants were right-censored. After adjustment for age, baseline severity, and biomarker level, intervention was associated with a lower conditional instantaneous event rate (Cox HR = 0.67, 95% CI 0.55–0.82). The Weibull AFT model estimated that conditional event times were multiplied by 1.31 (95% CI 1.15–1.49). PH residuals showed no clear global departure, Weibull had the lowest AIC among candidate parametric distributions, and the Cox–Snell curve remained close to the reference line over the main range. Because management strategy was not randomized and neither independent censoring nor model form can be verified completely, these associations should not be interpreted as real-world causal treatment effects.

6.5.1 Reporting template

In [target population], participants were followed from [time zero] to [explicit event] for [duration]. After adjustment for [prespecified covariates], the [HR/TR] for [exposure] versus [comparator] was [estimate] (95% CI: [lower, upper]). At [prespecified time or covariate profile], the model-estimated survival probabilities were [values]. Diagnostics for [PH/AFT distribution/censoring/functional form] showed [findings]. Because of [specific design or data limitation], the result should be interpreted as [an association/a conditional prediction/a causal effect under additional conditions].

Mini-case reflection: why report absolute survival probabilities as well?

If both the HR and TR have confidence intervals, why also provide 12-month survival probabilities?

Answer: A relative effect does not reveal the baseline event level. The same HR can correspond to very different absolute benefits in low- and high-risk populations, and a TR does not directly provide event probability at a fixed time. Absolute outcomes align more closely with many decisions, but the corresponding population, covariate profile, and time point must be stated.

7 Common Errors at a Glance

Common statement or practice	Problem	Better approach
Treating event=0 as the event	Reverses the directions of event and censoring	Verify explicitly that 1 = event and 0 = censored, and inspect raw counts
“HR 0.70 means a 30% lower one-year risk”	Treats an instantaneous rate ratio as a cumulative risk ratio	Name the HR correctly and report one-year survival probability or risk separately
Treating exp(survreg coefficient) as an HR	An exponentiated AFT coefficient is a TR	Interpret it on the time scale; convert only under a compatible Weibull model
Treating survreg$scale as the Weibull shape	R uses the reciprocal parameterization	Use shape = 1 / fit$scale
Declaring PH valid because cox.zph p > 0.05	Failure to reject is not proof	Combine residual plots, power, numbers at risk, and subject-matter knowledge
Automatically switching to AFT when PH fails	Failure of PH does not guarantee constant time scaling	Assess the AFT distribution and acceleration-factor assumption independently
Comparing Cox and AFT AIC values	Partial- and full-likelihood bases differ	Compare only candidate parametric models that share the same full-likelihood basis
Assuming robust standard errors repair the model	A variance correction cannot fix non-PH, nonlinearity, or an incorrect risk set	Correct the model structure and perform targeted sensitivity analyses
Predicting beyond the last observed time without disclosure	Results are driven mainly by tail-distribution assumptions	Mark the observed range clearly and compare extrapolations across reasonable distributions
Reporting only an HR or TR	Omits the baseline level and decision-relevant absolute quantities	Also report survival probability, risk, or a time quantile at a prespecified time

8 Quick Reference

8.1 Core formulas

Concept	Formula	Principal interpretation
Survival function	$S(t)=P(T>t)$	Probability of remaining event-free through $t$
Cumulative hazard	$H(t)=-\log S(t)$	Accumulation of instantaneous hazard over time
Cox PH	$h(t\mid X)=h_0(t)e^{X\beta}$	$e^{\beta}$ is a conditional HR
AFT	$\log T=X\beta+\sigma\varepsilon$	$e^{\beta}$ is a conditional TR
Weibull parameter bridge	$\kappa=1/\sigma$	The reciprocal of survreg scale is shape
Weibull conversion	$HR=TR^{-\kappa}$	Valid only for compatible Weibull PH/AFT models
Cumulative event probability at a fixed time	$F(t)=1-S(t)$	Probability of an event by $t$ in the absence of competing risks

8.2 Common R code

Goal	Code pattern
Construct a right-censored outcome	survival::Surv(time, event)
Kaplan–Meier curve	survival::survfit(survival::Surv(time, event) ~ group, data = d)
Cox PH	survival::coxph(survival::Surv(time, event) ~ x1 + x2, data = d)
PH diagnostic	survival::cox.zph(cox_fit)
Adjusted Cox curve	survival::survfit(cox_fit, newdata = profiles)
Weibull AFT	survival::survreg(survival::Surv(time, event) ~ x1 + x2, data = d, dist = "weibull")
AFT time ratios	exp(coef(aft_fit)[names(coef(aft_fit)) != "(Intercept)"])
AFT conditional quantile	predict(aft_fit, newdata = profiles, type = "quantile", p = 0.5)
Candidate parametric-model AIC	AIC(aft_weibull, aft_lognormal)

8.3 Pre-modeling checklist

• Are the time origin, event, time scale, and observation endpoint explicit?
• Is the event truly coded 1 and censoring coded 0?
• Are delayed entry, competing events, recurrent events, or clustering present?
• What are the event count, censoring proportion, and number at risk in each group?
• Is the primary target an HR, TR, fixed-time risk, quantile, or RMST?
• Were the units, centering, and functional forms of continuous variables prespecified?
• Were unsupported arbitrary cutpoints for continuous variables avoided?
• Were adjustment variables measured before time zero, and do they have a defensible causal basis?
• Are there enough events to support the fitted parameters and interactions?

8.4 Post-fit diagnostic checklist

8.4.1 Cox PH

• examine KM curves, log-minus-log plots, and numbers at risk;
• examine global and covariate-specific cox.zph() results and residual plots;
• assess continuous-variable nonlinearity, unusual observations, and dfbeta influence;
• evaluate clustering, tied events, and the censoring mechanism;
• assess absolute predictions and calibration at prespecified times.

8.4.2 AFT

• compare scientifically reasonable candidate distributions rather than searching every available distribution;
• assess Cox–Snell residuals and agreement between observed and predicted survival within groups;
• assess continuous-variable functional forms on the log-time scale;
• evaluate sensitivity of key conclusions to the distributional tail and censoring assumptions;
• state the survreg distribution, scale, and software parameterization explicitly;
• separate extrapolated results clearly from the actual observed range.

8.5 Plain-language glossary

Term	Meaning
Risk set	People still observed and event-free immediately before an event time
Right censoring	Knowing only that event time is later than the last observation
Delayed entry	A participant begins contributing to the risk set after time zero
Survival probability	Probability of remaining free of the target event through a specified time
Hazard function	Instantaneous event rate immediately afterward among people still event-free
HR	Ratio of two conditional instantaneous event rates
TR	Multiplier for conditional event-time quantiles in an AFT model
Proportional hazards	An HR that remains constant over the analysis time range
Acceleration factor	Multiplier that stretches or compresses the event-time scale in an AFT model
Baseline hazard	Hazard function in a Cox model when all covariates are at their reference values
Partial likelihood	Cox method that uses event ordering and risk sets to estimate relative-hazard coefficients
Extrapolation	Prediction beyond the time range actually supported by the data

9 Final Knowledge Check

1. Does a right-censored participant still contribute information to the risk set before censoring?
2. Does HR = 0.60 necessarily imply a 12-month risk ratio of 0.60?
3. Why does Cox PH not require choosing a Weibull or log-normal distribution for $h_0(t)$?
4. Does a global cox.zph() p-value of 0.40 prove that PH holds?
5. What does $\exp(\beta)=1.40$ mean in an AFT model?
6. If survreg(dist="weibull") reports scale = 0.80, what is the Weibull shape?
7. Can a TR from a log-normal AFT model usually be converted to a constant HR?
8. Why should AIC not be used to compare a standard Cox model directly with a parametric AFT model?
9. If a competing death is treated as ordinary censoring, is $1-\widehat S(t)$ still the cumulative incidence of recurrence?
10. Does a statistically significant multivariable-adjusted HR automatically have a causal interpretation?

Show final answers

1. Yes. Information that the participant remained event-free before the censoring time contributes to the corresponding risk sets.
2. No. An HR is a conditional instantaneous rate ratio, whereas a fixed-time risk ratio compares cumulative probabilities.
3. Cox partial likelihood estimates $\beta$ through risk-set comparisons at event times without parameterizing the shape of the baseline hazard.
4. No. It means only that the current data provide no strong evidence against PH; power, plots, and subject-matter knowledge still matter.
5. Holding other covariates fixed and assuming the AFT model holds, every conditional event-time quantile is multiplied by 1.40, extending the time scale by approximately 40%.
6. Shape = $1/0.80=1.25$.
7. Usually not. A log-normal AFT model generally does not satisfy constant proportional hazards.
8. Cox uses a partial likelihood, while a parametric AFT model uses a full event-time likelihood, so their AIC values have different likelihood bases.
9. Usually not. A cumulative incidence function is needed to express real-world cumulative probability in the presence of competing risks.
10. No. Appropriate study design, temporal ordering, control of confounding, measurement, censoring, and identification assumptions are still required.

Next steps

Further topics include restricted mean survival time, flexible parametric survival models, splines and time-varying effects, time-varying covariates, competing risks, multistate models, recurrent events, frailty, causal survival analysis, weighting methods for informative censoring, multiple imputation, 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  stats4_4.6.1   
## [13] lifecycle_1.0.5 cli_3.6.6       grid_4.6.1      sass_0.4.10    
## [17] jquerylib_0.1.4 compiler_4.6.1  tools_4.6.1     evaluate_1.0.5 
## [21] bslib_0.12.0    survival_3.8-6  yaml_2.3.12     rlang_1.3.0    
## [25] jsonlite_2.0.0

Teaching use: This module supports learning about estimation, diagnosis, and interpretation of survival models. It does not replace a statistical analysis plan, data audit, or professional biostatistical consultation for a specific study. All data and numerical results are simulated.