Cox Proportional Hazards: A Method Guide

A Worked Example With Made-Up Data

Question: In a made-up study, can a Cox model find a set effect on time to an event?

Inputs: Simulated data, seed 17. The script makes 300 cases in two groups. It builds in a hazard ratio of 0.70. It draws event and censor times on their own, with base rates of 0.08 and 0.04.

Method: The script fits a Cox model by its partial likelihood. It gets a Wald 95% range from the model's information. It repeats the fit across seeds 0 to 999 to show the spread. The code uses NumPy and SciPy.

Result: The fit gave a hazard ratio of 0.694, with a 95% range from 0.525 to 0.918, from 201 events. Across seeds, the median was 0.695; the 5th to 95th percentile ran from 0.543 to 0.884. The event count had a median of 188, with a 5th to 95th percentile from 174 to 201. The range covers the built-in 0.70 and the median sits at it, which is expected from how the data were built. This checks the code, not real people or firms.

Limit: The code sets the effect to stay fixed over time, so the proportional hazards assumption holds by construction. The check for that assumption was not run. It does not show if that rule holds in real data.

import numpy as np
from scipy.optimize import minimize_scalar

TRUE_BETA = np.log(0.70)
BASE_HAZARD = 0.08
CENSOR_HAZARD = 0.04
N = 300


def fit(seed):
    rng = np.random.default_rng(seed)
    x = rng.binomial(1, 0.5, N)
    event_time = rng.exponential(1 / (BASE_HAZARD * np.exp(TRUE_BETA * x)))
    censor_time = rng.exponential(1 / CENSOR_HAZARD, N)
    time = np.minimum(event_time, censor_time)
    event = event_time <= censor_time

    def neg_loglik(beta):
        total = 0.0
        for t in time[event]:
            risk = time >= t
            z = beta * x[risk]
            m = np.max(z)
            total += beta * x[(time == t) & event].sum() - (m + np.log(np.exp(z - m).sum()))
        return -total

    result = minimize_scalar(neg_loglik, bounds=(-4, 4), method="bounded")
    beta = result.x
    info = 0.0
    for t in time[event]:
        risk_x = x[time >= t]
        weights = np.exp(beta * risk_x)
        mean = np.sum(weights * risk_x) / np.sum(weights)
        variance = np.sum(weights * (risk_x - mean) ** 2) / np.sum(weights)
        info += variance
    se = 1 / np.sqrt(info)
    hr = np.exp(beta)
    lo, hi = np.exp(beta - 1.96 * se), np.exp(beta + 1.96 * se)
    return hr, lo, hi, int(event.sum())


seed = 17
hr, lo, hi, events = fit(seed)
runs = np.array([fit(s) for s in range(1000)])
hrs = runs[:, 0]
event_counts = runs[:, 3]
print(f"Simulated data, seed {seed}; subjects={N}; one binary group covariate")
print("Built-in effect: hazard ratio=0.70; event hazard base=0.08; censor hazard=0.04")
print(f"One draw: events={events}; fitted hazard ratio={hr:.3f}; 95% CI=({lo:.3f}, {hi:.3f})")
print(f"Across 1,000 seeds: fitted hazard ratio median={np.median(hrs):.3f}; 5th-95th percentile=({np.percentile(hrs, 5):.3f}, {np.percentile(hrs, 95):.3f})")
print(f"Across 1,000 seeds: event count median={np.median(event_counts):.0f}; 5th-95th percentile=({np.percentile(event_counts, 5):.0f}, {np.percentile(event_counts, 95):.0f})")
print("Diagnostic: proportional hazards was built into the simulation; this run does not test that assumption.")

Printed output from running the code above:

Simulated data, seed 17; subjects=300; one binary group covariate
Built-in effect: hazard ratio=0.70; event hazard base=0.08; censor hazard=0.04
One draw: events=201; fitted hazard ratio=0.694; 95% CI=(0.525, 0.918)
Across 1,000 seeds: fitted hazard ratio median=0.695; 5th-95th percentile=(0.543, 0.884)
Across 1,000 seeds: event count median=188; 5th-95th percentile=(174, 201)
Diagnostic: proportional hazards was built into the simulation; this run does not test that assumption.

1. Introduction

1.1 Problem Statement

Survival analysis is the statistical framework for data where the outcome is the time until an event occurs and some subjects have not yet had the event when observation ends. The Cox proportional hazards model, introduced by David Cox in 1972, is the most widely taught semi-parametric approach for relating covariates to the hazard of that event. It handles censored data, estimates covariate effects without specifying a baseline hazard, and reports effects as hazard ratios.

The same properties that make the model easy to fit make it easy to fit carelessly. Software returns hazard ratios and p-values whether or not the assumptions behind them hold. This guide explains those assumptions and the standard checks for them.

Typical settings include time to patient relapse in medicine, time to component failure in reliability engineering, and time to cancellation in a subscription business. The mathematics is identical in each; only the meaning of the event and of censoring changes.

1.2 Scope and Objectives

This guide presents no original research or survey of practice. It includes a simulated teaching example, not evidence from real people or firms. Its aims are to:

• State the model and the partial likelihood used to fit it
• Explain how to interpret hazard ratios and confidence intervals
• Describe the proportional hazards assumption and how to diagnose violations
• Describe time-varying covariates, stratification and tie handling
• List well-known pitfalls and give generic guidance for avoiding them

Numbers appear only in the formulas and in one worked example that is labelled as made up. Published figures are attributed to the cited sources.

2. Background

2.1 The Cox Proportional Hazards Model

The Cox model specifies the hazard function for individual i at time t as:

The model is semi-parametric: it estimates covariate effects without requiring a parametric form for the baseline hazard. Covariate effects are reported as hazard ratios, HR = exp(β).

The proportional hazards assumption, that hazard ratios are constant over time, is the model's central requirement. When it holds and covariate values stay fixed, the hazard for any individual stays in a fixed ratio to the hazard for any other individual throughout follow-up; with time-dependent covariates the ratio can change as the covariates change. When it fails, the single hazard ratio the model reports is a weighted average of effects that differ across time.

2.2 The Partial Likelihood

Cox's insight was to estimate β without ever estimating h₀(t). At each observed event time t_i, consider the risk set R(t_i), the subjects still under observation just before t_i. The probability that the subject who actually had the event is the one who did, given that exactly one event occurred, depends only on β:

2.3 Interpreting Hazard Ratios

A hazard ratio of 1 means no association. A hazard ratio above 1 means a higher instantaneous event rate for a one-unit increase in the covariate (or for the exposed group relative to the reference group), and below 1 means a lower rate. A hazard ratio is a ratio of rates, not a ratio of probabilities or of survival times. A confidence interval that excludes 1 corresponds to a Wald test rejecting the null at the matching level. Interpretation for a continuous covariate is per unit of that covariate as coded, so rescaling the covariate changes the reported HR.

2.4 Alternatives and Extensions

The Cox model is one tool among several. Parametric models (Weibull, exponential, log-normal) specify the baseline hazard and allow direct prediction of survival times. Accelerated failure time models describe covariate effects on the time scale. Stratified Cox models allow separate baseline hazards per stratum. Models with time-varying coefficients relax proportionality; time-dependent covariates handle covariates that change over time, which is a separate issue. Competing-risks models apply when other events preclude the event of interest. For the nonparametric starting point, see the Kaplan-Meier whitepaper.

3. Method Overview

3.1 How This Guide Is Organized

This section describes the standard workflow for fitting and checking a Cox model. It is a description of established practice from the textbooks and papers in the references.

1. Define the event, the time origin and the time scale, and code the censoring indicator.
2. Plot Kaplan-Meier curves by the main groups to see the data before modeling.
3. Choose covariates and check that the number of events supports the number of parameters.
4. Fit the model, stating the tie-handling method.
5. Check proportional hazards, functional form and influential observations.
6. If the assumption fails, stratify, add time-varying terms, or change the model.
7. Report hazard ratios with confidence intervals, the event count, and the diagnostics performed.

3.2 Data Considerations

Censoring mechanisms: Right censoring (follow-up ends before the event) is the standard case. Left and interval censoring need different methods. The Cox model assumes censoring is non-informative: the reason a subject is censored is unrelated to their event risk once the covariates are accounted for.

Tied event times: Ties arise when events are recorded at coarse granularity, such as days. The exact methods handle tied events either by averaging over their possible orderings or by summing over the possible sets of subjects who had the event, and both are costly with many ties, so approximations such as the Breslow and Efron methods are used.

Time scale: Time since entry, calendar time and age are all possible time scales. The choice changes the risk sets and can change whether proportionality looks plausible.

Covariate types: Continuous or categorical, constant over time or time-varying, linear or nonlinear in the log hazard. Each affects which diagnostics matter.

3.3 Diagnostic Tools

Proportional hazards: The Grambsch and Therneau test correlates scaled Schoenfeld residuals with a function of time, globally and per covariate. Plots of scaled Schoenfeld residuals against time show the direction of any drift. Log-log survival plots by group are a visual check for categorical covariates, and a covariate-by-time interaction can be tested with a likelihood ratio test.

Functional form: Martingale residuals plotted against a continuous covariate suggest whether a transformation such as a spline is needed.

Outliers and influence: Deviance residuals flag poorly fit observations. Score residuals are each observation's contribution to the fitting equations, and dfbeta approximates how much each coefficient would change if that observation were removed.

Fit and discrimination: The concordance index summarizes how well the risk score orders subjects by event time. Calibration should be checked separately when predictions are used.

In R, the survival package provides coxph, cox.zph and residual methods; the lifelines package offers comparable tools in Python.

4. Key Points

5. Worked Example (Illustrative)

5.1 Reading a Hazard Ratio

Suppose a subscription service fits a Cox model for time to cancellation with one binary covariate, onboarding call (1 = received a call, 0 = did not), and obtains a coefficient β = −0.357.

• The hazard ratio is HR = exp(−0.357) ≈ 0.70.
• Interpretation: at any given time, customers who received a call cancel at about 0.70 times the rate of otherwise similar customers who did not, which is a 30% lower instantaneous rate.
• Suppose the standard error of β is 0.12. The 95% confidence interval on the log scale is −0.357 ± 1.96 × 0.12 = (−0.592, −0.122), which gives an HR interval of (exp(−0.592), exp(−0.122)) ≈ (0.55, 0.89).

This says nothing about whether the call caused the difference. A Cox model estimates association; causal reading needs a design that supports it, such as randomization or a credible adjustment set.

5.2 What a Proportionality Violation Looks Like

Suppose that, in the same invented example, the scaled Schoenfeld residuals for onboarding call trend upward with time, and the Schoenfeld test returns a small p-value. That pattern would mean the effect of the call is strongest early and fades, so HR ≈ 0.70 is an average over a changing effect. A reasonable response is to fit a model in which the coefficient depends on time, for example β(t) = β₀ + β₁·log(t), and report the hazard ratio at several time points instead of one number.

5.3 Events per Variable

Suppose the dataset has 40 cancellations and the analyst wants to include 8 covariates. That is 40 / 8 = 5 events per variable, below the common heuristic of 10. The options are to drop covariates on subject-matter grounds, fit a penalized model, or wait for more follow-up, as described in Point 3.

6. Practical Considerations

Computation: Exact tie handling can be slow when there are many ties in large datasets. Efron is the usual compromise, and discrete-time models are an alternative framework.

Multiple testing: A thorough set of diagnostics runs many tests. Read them together, and prefer plots alongside p-values.

Interpretability: Stratified and time-varying models are harder to summarize than one hazard ratio. Decide in advance how the results will be communicated to non-statistical readers, and show survival curves where possible, since a curve is easier to read than a ratio of rates.

Prediction versus explanation: A Cox model used for prediction needs calibration and validation on data not used for fitting. A Cox model used for explanation needs a defensible set of adjustment variables.

Documentation: Record the time origin, event definition, censoring rules, tie method, covariate coding, diagnostics run and any model changes made after diagnostics, so the analysis can be reproduced.

7. Recommendations

The recommendations below are general guidance drawn from the standard references. They are not conclusions from any data analyzed for this paper.

8. Conclusion

The Cox proportional hazards model is a flexible and widely used tool for time-to-event data. Its output is only as meaningful as the assumptions behind it: constant hazard ratios over time, non-informative censoring, enough events for the number of parameters, a stated treatment of ties, and a functional form that matches the data.

This paper has explained the model, its partial likelihood, the interpretation of hazard ratios, and the standard diagnostics and remedies, and it has illustrated them with a worked example using made-up numbers. It reports no measured results. Readers who want empirical evidence on any point should consult the primary sources in the references.

The practical takeaway is a short checklist: define the event and censoring carefully, count events against parameters, state the tie method, test proportional hazards, look at residuals, and report the number of events with the hazard ratios.