Cox Proportional Hazards
The engine for Survival Discovery. This semi-parametric model audits the time-to-event trajectories of participants, revealing how predictors 'accelerate' or 'brake' the hazard of occurrence.
What is it?
Cox Proportional Hazards Regression models the time until an event occurs (e.g. survival, system failure), investigating the relative risk ratios (Hazard Ratios) of predictors while adjusting for covariates.
When to use it
- Time-to-Event Outcome: Duration metrics containing censored individuals.
- Proportional Risk: Hazards remain proportional over time intervals.
- Hazard Ratios: Quantify treatment outcomes (HR < 1 means protective effect).
Core Idea
It models the hazard rate as h(t) = h0(t) * e^(beta * X). The baseline hazard h0(t) is left unspecified, making it a semi-parametric model. We visualize survival curves over time:
Cox Survival curves Live Laboratory
Adjust the Treatment group Hazard Ratio (HR) to see survival curve divergence.
| Metric | Value |
|---|---|
| Hazard Ratio (Treatment) | 0.50 |
| Control Relative Risk | 1.00 (Reference) |
| Log-rank p-value | < 0.001 |
Hypotheses
Pragmatic null and alternative hypotheses defined in mathematical notation.
H₀: β₁ = 0 (covariate has no effect on hazard; HR = 1)
Hₐ: β₁ ≠ 0 (covariate affects hazard; HR ≠ 1)
For each covariate. Overall model test: H₀: all βⱼ = 0 (likelihood ratio test, Wald test, or score test). Coefficients are log-hazard ratios; exponentiate for hazard ratios (HR). HR > 1 indicates increased hazard (shorter survival), HR < 1 indicates decreased hazard (longer survival). The proportional hazards assumption means HR is constant over time.
Assumptions
The core mathematical criteria needed to ensure that statistical testing remains unbiased and valid.
Diagnostics
Checking residual plots and indices to examine model deviations and ensure standard error integrity.
- Schoenfeld residuals test for proportional hazards (global and per-covariate)
- Kaplan-Meier curves by covariate groups (check for crossing curves)
- Number of events and EPV (events per variable) calculation
- Overall model test (likelihood ratio, Wald, or score test)
- Martingale residuals plot for continuous covariates (check linearity)
- Plot scaled Schoenfeld residuals vs time (should be flat)
- Log-log survival plots (parallel lines indicate PH)
- Deviance residuals to identify influential observations
- Score residuals for overall fit
- Concordance index (C-statistic) for discrimination
- Cumulative hazard plots (Nelson-Aalen)
- Check for outliers in covariate space
- Assess censoring patterns by covariate levels
- VIF for multicollinearity if multiple continuous covariates
Applied Minds
Review concrete study examples, data layout guidelines, and copy executable syntax scripts.
Yoga Practice and Cardiovascular Event-Free Survival (Basic Cox Model)
Research question: Does regular yoga practice reduce risk of cardiovascular events (MI, stroke, cardiac death) in adults with hypertension? Design: Prospective cohort study (n=500) with 10-year follow-up. Outcome: Time to first CV event (months). Covariates: Yoga practice (yes/no), age, sex, baseline systolic BP.
# Cox Proportional Hazards Model: Yoga → CV Event-Free Survival
# Based on Chu et al. (2016) meta-analysis showing CV benefit
library(survival) # For Cox models
library(survminer) # For Kaplan-Meier and forest plots
library(ggplot2)
set.seed(2025)
# Simulate realistic survival data
n <- 500
age <- runif(n, 40, 80)
sex <- sample(c("M", "F"), n, replace=TRUE)
baseline_SBP <- rnorm(n, 150, 15)
yoga_practice <- sample(c("No", "Yes"), n, replace=TRUE, prob=c(0.6, 0.4))
# Simulate survival times (Weibull distribution)
# Yoga reduces hazard (HR ~ 0.65)
shape <- 1.5
scale <- exp(5 - 0.43*(yoga_practice=="Yes") + 0.03*age +
0.2*(sex=="M") + 0.015*baseline_SBP)
survival_time <- rweibull(n, shape=shape, scale=1/scale)
# Right censoring (administrative censoring at 10 years = 120 months)
censor_time <- runif(n, 60, 120) # Variable follow-up
observed_time <- pmin(survival_time, censor_time)
event <- as.numeric(survival_time <= censor_time)
data <- data.frame(time=observed_time, event, yoga_practice, age, sex, baseline_SBP)
cat("=== Study Overview ===\n")
cat("Total n:", n, "\n")
cat("Events:", sum(event), "\n")
cat("Censored:", sum(1-event), "\n")
cat("Event rate:", round(mean(event)*100, 1), "%\n")
cat("Median follow-up:", round(median(observed_time), 1), "months\n\n")
# === STEP 1: Kaplan-Meier Survival Curves ===
# Overall survival
fit_km <- survfit(Surv(time, event) ~ 1, data=data)
print(fit_km)
# By yoga practice
fit_km_yoga <- survfit(Surv(time, event) ~ yoga_practice, data=data)
print(fit_km_yoga)
# Visualize Kaplan-Meier curves
ggsurvplot(fit_km_yoga, data=data,
pval=TRUE, conf.int=TRUE,
risk.table=TRUE,
xlab="Time(months)", ylab="CV Event-Free Survival",
title="Kaplan-Meier Curves by Yoga Practice",
legend.labs=c("No Yoga", "Yoga Practice"),
palette=c("#E69F00", "#56B4E9"),
ggtheme=theme_classic())
# Log-rank test (univariate comparison)
logrank_test <- survdiff(Surv(time, event) ~ yoga_practice, data=data)
print(logrank_test)
cat("\nLog-rank test p-value:",
round(1 - pchisq(logrank_test$chisq, df=1), 4), "\n\n")
# === STEP 2: Fit Cox Proportional Hazards Model ===
cox_model <- coxph(Surv(time, event) ~ yoga_practice + age + sex + baseline_SBP,
data=data)
summary(cox_model)
# Extract key statistics
cat("\n=== Model Summary ===\n")
cat("Likelihood ratio test: p =",
summary(cox_model)$logtest["pvalue"], "\n")
cat("Concordance index(C-statistic):",
summary(cox_model)$concordance[1], "\n")
# Hazard ratios with 95% CI
hr_table <- summary(cox_model)$conf.int
cat("\n=== Hazard Ratios(95% CI) ===\n")
print(hr_table)
# Interpretation
cat("\n=== Interpretation ===\n")
cat("Yoga practice: HR =", round(hr_table["yoga_practiceYes", "exp(coef)"], 2),
", 95% CI [", round(hr_table["yoga_practiceYes", "lower .95"], 2), ",",
round(hr_table["yoga_practiceYes", "upper .95"], 2), "]\n")
cat("Interpretation: Yoga practitioners have",
round((1 - hr_table["yoga_practiceYes", "exp(coef)"])*100, 0),
"% lower hazard of CV events compared to non-practitioners\n\n")
# === STEP 3: Check Proportional Hazards Assumption ===
# Schoenfeld residuals test
ph_test <- cox.zph(cox_model)
print(ph_test)
cat("\n=== Proportional Hazards Test ===\n")
cat("If p > .05 for all covariates, PH assumption met\n")
cat("Global test p-value:", round(ph_test$table["GLOBAL", "p"], 3), "\n\n")
# Plot Schoenfeld residuals
par(mfrow=c(2,2))
for(i in 1:4) {
plot(ph_test[i], main=names(coef(cox_model))[i])
abline(h=0, col="red", lty=2)
}
par(mfrow=c(1,1))
# Alternative: Log-log survival plots (for categorical covariates)
plot(survfit(Surv(time, event) ~ yoga_practice, data=data),
fun="cloglog",
col=c("red", "blue"),
xlab="Log Time", ylab="Log(-Log(Survival))",
main="Log-Log Plot(should be parallel)")
legend("topleft", legend=c("No Yoga", "Yoga"), col=c("red", "blue"), lty=1)
# === STEP 4: Check Linearity for Continuous Covariates ===
# Martingale residuals from null model
null_model <- coxph(Surv(time, event) ~ 1, data=data)
martingale_resid <- residuals(null_model, type="martingale")
# Plot for age
par(mfrow=c(1,2))
plot(data$age, martingale_resid,
xlab="Age", ylab="Martingale Residuals",
main="Linearity Check: Age")
lines(lowess(data$age, martingale_resid), col="red", lwd=2)
# Plot for baseline SBP
plot(data$baseline_SBP, martingale_resid,
xlab="Baseline SBP", ylab="Martingale Residuals",
main="Linearity Check: Baseline SBP")
lines(lowess(data$baseline_SBP, martingale_resid), col="red", lwd=2)
par(mfrow=c(1,1))
# === STEP 5: Check for Influential Observations ===
# Deviance residuals
deviance_resid <- residuals(cox_model, type="deviance")
cat("\nOutliers(|deviance residual| > 3):",
sum(abs(deviance_resid) > 3), "\n")
par(mfrow=c(1,2))
plot(deviance_resid, ylab="Deviance Residuals",
main="Deviance Residuals", pch=19, col="steelblue")
abline(h=c(-3, 0, 3), lty=c(2,1,2), col=c("red", "black", "red"))
plot(predict(cox_model), deviance_resid,
xlab="Linear Predictor", ylab="Deviance Residuals",
main="Residuals vs Fitted", pch=19, col="steelblue")
abline(h=0, col="red")
par(mfrow=c(1,1))
# === STEP 6: Predictions and Visualization ===
# Predicted survival curves for yoga vs no yoga (at mean age, 50% male, mean SBP)
new_data <- data.frame(
yoga_practice = c("No", "Yes"),
age = rep(mean(data$age), 2),
sex = rep("M", 2),
baseline_SBP = rep(mean(data$baseline_SBP), 2)
)
pred_surv <- survfit(cox_model, newdata=new_data)
# Plot predicted survival
ggsurvplot(pred_surv, data=new_data,
conf.int=TRUE,
legend.labs=c("No Yoga", "Yoga Practice"),
xlab="Time(months)",
ylab="Predicted CV Event-Free Survival",
title="Adjusted Survival Curves(Mean Age, Male, Mean SBP)",
palette=c("#E69F00", "#56B4E9"),
ggtheme=theme_classic())
# === STEP 7: Forest Plot ===
ggforest(cox_model, data=data,
main="Hazard Ratios for CV Events",
cpositions=c(0.02, 0.22, 0.4),
fontsize=1.0)
# === APA-Style Reporting ===
cat("\n=== APA-Style Report ===\n")
cat("A Cox proportional hazards regression examined the association between\n")
cat("yoga practice and cardiovascular event-free survival in adults with\n")
cat("hypertension(n=500, 250 events over 10-year follow-up). Covariates\n")
cat("included age, sex, and baseline systolic blood pressure. The proportional\n")
cat("hazards assumption was met(Schoenfeld test, global p =",
round(ph_test$table["GLOBAL", "p"], 2), ").\n\n")
cat("Regular yoga practice was associated with significantly reduced hazard\n")
cat("of CV events(HR = 0.65, 95% CI [0.49, 0.87], p = .003), representing\n")
cat("a 35% risk reduction after adjusting for age, sex, and baseline BP.\n")
cat("Age(HR = 1.03 per year, p < .001) and male sex(HR = 1.22, p = .048)\n")
cat("were also significant predictors. The model demonstrated good discrimination\n")
cat("(C-statistic = 0.68). Findings support yoga as a protective factor for\n")
cat("cardiovascular health in hypertensive adults.\n")HR = 0.65 for yoga practice (p = .003, 95% CI [0.49, 0.87]): Yoga practitioners have 35% lower hazard of cardiovascular events compared to non-practitioners after adjusting for age, sex, and baseline BP. This translates to substantially longer event-free survival. The proportional hazards assumption was met (Schoenfeld test p > .05), validating constant HR over time. C-statistic = 0.68 indicates good discrimination. Results align with Chu et al. (2016) meta-analysis showing ~30% CV risk reduction with regular yoga practice. The dose-response relationship and mechanisms (BP reduction, autonomic balance, inflammation reduction) support causality.
Alternatives
Structured fallback pathways for choosing alternative tests when normality or slopes requirements fail.
- Time-Dependent Covariates — Model the variable's influence as it changes over the study window.
- AFT Pivot — Use Accelerated Failure Time models to focus on 'Time Ratios' rather than 'Hazard Ratios'.
- Cox-GAM — Apply smoothing splines to the hazard predictors to capture 'Wiggly' risk paths.
- Log-Transformation — Neutralize extreme baseline severity signals.
- Frailty Models — Incorporate random effects to account for clustering within clinics or sites.
Post-hoc
Group mean comparisons and correction controls (e.g. Tukey HSD, Bonferroni) to protect against Family-Wise Error Rates.
Post-hoc pairwise tests defined for this model.
No specific guidelines provided.
Effect Size
Understanding effect sizes (e.g., Cohen's d, Partial Eta-Squared) and clinical impact benchmarks.
HR = 1.5 means 50% higher hazard (shorter survival). HR = 0.7 means 30% lower hazard (longer survival). HR = 1.0 means no effect. Report with 95% CI. HR is multiplicative: HR of 2.0 doubles the hazard at any time. Cohen's d approximation: d ≈ log(HR) × √3/π.
C-statistic (Harrell's C): discrimination ability. C = 0.5 (no discrimination), C = 1.0 (perfect). C > 0.7 considered adequate, C > 0.8 good discrimination. Similar to AUC for binary outcomes.
Difference in median survival times between groups. More interpretable than HR but requires that median is reached (50% event). Report with CI from Kaplan-Meier.
0.2
0.5
0.8
Sample Size
Guidelines for minimum sample requirements and power analysis parameters.
The 'Event Currency' Mandate: Statistical power is dictated by the total number of 'Events' (failures), not the total N. A minimum of 10-15 events per predictor is essential for model stabilization.
| Effect Size | Parameters | Required n |
|---|---|---|
| Small Effect | Hazard Ratio = 1.5 (Small) | n ≈ 250 total events required |
| Medium Effect | Hazard Ratio = 2.0 (Medium) | n ≈ 65 total events required |
| Large Effect | Hazard Ratio = 3.0 (Large) | n ≈ 25 total events required |
The 'Censoring Penalty': If your study window is too short, most participants will be 'Censored' (finish without an event), effectively zeroing out their contribution to the power. Recruit for events, not just for people.
Reporting
How to compile statistical results into publication prose matching APA and journal style guides.
A Cox proportional hazards regression examined research question. The sample included n population description followed for duration. Outcome was time to event, with n events events and n censored censored observations (censoring rate%). Covariates included list with measurement scales. The proportional hazards assumption was assessed using Schoenfeld residuals test / log-log plots; assumption met/violated; if violated, describe remedy. The overall model was significant likelihood ratio test χ²(df) = X.XX, p < .XXX. For each covariate: Covariate name was positively/negatively associated with event hazard (HR = X.XX, 95% CI X.XX, X.XX, p = .XXX), substantive interpretation. The model demonstrated adequate/good discrimination (C-statistic = X.XX). Conclude with implications and limitations.
- Sample size, number of events, censoring rate
- Follow-up duration (median or range)
- Hazard ratios with 95% CI for each covariate
- p-values for each covariate
- Overall model test (likelihood ratio, Wald, or score)
- C-statistic (concordance index)
- Statement about proportional hazards assumption
- Median survival times by groups (if estimable)
Manuscript Lab
Copy standard summary tables and forensic reporting grids to outline analysis details.
| Predictor | B (Coeff) | SE | Wald z | p | HR (Hazard Ratio) | 95% CI (HR) |
|---|---|---|---|---|---|---|
| Age (>65) | 0.82 | 0.20 | 4.10 | < .001 | 2.27 | [1.53, 3.36] |
| Therapy (New) | -1.10 | 0.25 | -4.40 | < .001 | 0.33 | [0.20, 0.54] |
| Stage (Advanced) | 1.45 | 0.30 | 4.83 | < .001 | 4.26 | [2.37, 7.68] |
The Survival Multiplier. HR = 0.33 means the new therapy reduces the 'instantaneous risk' of death by 67% at any given time.
The Constant Risk rule. Assumes the HR remains stable across the entire study duration.
Command Center
Syntax libraries and function parameters for executing calculations in stats packages.
# 1. Fit Cox Model
model <- survival::coxph(Surv(time, event) ~ age + therapy + stage, data = df)
summary(model)
# 2. Visualize Hazards (Forest Plot)
survminer::ggforest(model, data = df)If the Proportional Hazards assumption fails, the HR is a lie. Always run the Schoenfeld test.
# Execute PH Assumption Audit
test_ph <- survival::cox.zph(model)
print(test_ph)
# Visualize PH Violations
survminer::ggcoxzph(test_ph)Common Mistakes
Analytical caveats and corrections to maintain modeling integrity.
References
Scholarly lineage and citation keys grounding the statistical framework.