Cox Proportional Hazard Model

Motivating Example: Founder vs. Professional CEO Tenure#

A management researcher wants to know whether founder-CEOs stay in their positions longer than professional (externally hired) CEOs. She collects data on 1,200 CEO spells at publicly traded firms between 1995 and 2015, recording the date each CEO took office and, if applicable, the date they left.

Here is the problem: 40% of the CEOs in her sample are still in their position at the end of the study period in 2015. These cases are right-censored observations — she knows the CEO was still active at the end of 2015, but she does not know when they will eventually leave.

If she simply compares average observed tenure between the two groups, she will underestimate average tenure for both groups (because she is treating the end of the observation window as though it were the event date). Worse, if founder-CEOs are disproportionately censored (because they tend to stay longer), the bias will be asymmetric — she will underestimate founder-CEO tenure more than professional-CEO tenure, potentially masking the very difference she wants to detect.

She cannot simply drop the censored observations either. Doing so would select on the outcome: she would be left with only CEOs who departed, which is a non-random subset. If the censored CEOs are systematically different (e.g., better performing), dropping them introduces selection bias.

The Cox proportional hazard model (Cox, 1972) solves this problem. It models the instantaneous rate of departure (the hazard rate) as a function of covariates, while properly accounting for censored observations. Censored CEOs contribute information up to the point they are last observed — they were "at risk" of departing during all the time they were observed, and the model uses this information without requiring knowledge of when they eventually leave.

What the Cox Model Does#

The Cox proportional hazard model estimates hazard ratios — the multiplicative effect of covariates on the instantaneous rate of experiencing an event. The model is:

where:

• $h(t | X_i)$ is the hazard rate for individual $i$ at time $t$
• $h_0(t)$ is the baseline hazard — the hazard when all covariates equal zero
• $\exp(X_i'\beta)$ is the multiplicative shift due to covariates

The critical feature of the Cox model is that it is semiparametric: the baseline hazard $h_0(t)$ is left completely unspecified. You do not need to assume that it follows any particular distribution (exponential, Weibull, log-normal, etc.). The model only estimates how covariates shift the hazard, not the shape of the hazard itself.

Key Concepts#

• Survival function $S(t) = P(T > t)$: the probability of surviving (not experiencing the event) beyond time $t$. The Kaplan-Meier estimator provides a nonparametric estimate of this function (Kaplan & Meier, 1958).

• Hazard rate $h(t)$: the instantaneous risk of the event at time $t$, conditional on having survived to $t$. Unlike a probability, the hazard rate can exceed 1 because it is a rate (events per unit time), not a probability.

• Right censoring: when the event has not occurred by the end of the observation period. The Cox model handles this by including censored observations in the risk set up to the censoring time but not requiring them to contribute an event.

• Proportional hazards: the assumption that the ratio of hazards for any two individuals is constant over time. Proportional hazards is the defining assumption of the Cox model.

How It Differs from OLS#

The key difference from OLS is that the Cox model correctly handles censored observations. In a regression of observed duration on covariates, OLS treats censored observations as though their observed time is the true duration — biasing coefficients downward. The Cox model instead uses the partial likelihood, which conditions on the observed ordering of events and does not require knowing the full distribution of event times.

When to Use the Cox Model#

• Your outcome is time-to-event: time to firm exit, CEO departure, technology adoption, employee turnover, patent citation, project completion, loan default
• You have right-censored observations: subjects who have not experienced the event by the end of your study
• You want to estimate how covariates affect the rate of the event rather than a binary yes/no outcome
• You do not want to assume a particular parametric form for the baseline hazard

When NOT to Use the Cox Model#

• Your outcome is binary with no meaningful time dimension — use logit/probit
• Your outcome is a continuous variable that is not a duration — use OLS
• The proportional hazards assumption is violated and cannot be fixed by stratification — consider accelerated failure time (AFT) models
• You have multiple competing event types (e.g., CEO can leave voluntarily, be fired, or retire) and you want event-specific effects — use competing risks models (Fine & Gray, 1999)

When to Use (Detailed)#

1. Your outcome is a duration with censoring. Time to firm exit, CEO tenure, time to technology adoption, employee turnover, time to loan default, patent lifetime — any setting where you observe a time-to-event and some subjects have not yet experienced the event.

2. You want to estimate covariate effects without specifying the baseline hazard. If you do not have a strong prior about the shape of the hazard over time (whether it increases, decreases, or is constant), the semiparametric Cox model is the safe choice.

3. You want to compare groups controlling for covariates. The Cox model provides hazard ratios that quantify how much faster or slower the event occurs for one group relative to another, holding other factors constant.

When NOT to Use (Detailed)#

1. The proportional hazards assumption is badly violated. If a treatment effect wears off sharply over time, the Cox model produces a single "average" hazard ratio that may not describe any actual time period well. Consider:

   • Stratified Cox model: allows different baseline hazards for different groups
   • Time-varying coefficients: interact the covariate with time
   • Accelerated failure time (AFT) models: model log(time) directly

2. You have interval-censored data. If you only know that the event occurred between two time points (e.g., between annual surveys), use discrete-time hazard models rather than the standard Cox model.

3. Multiple event types compete. If a CEO can leave voluntarily, be fired, or retire, you must choose an appropriate competing risks approach. To estimate cause-specific hazard ratios (how covariates affect the instantaneous rate among those still event-free), fit separate Cox models for each event type, censoring the others — standard Cox software handles this directly. To estimate effects on the cumulative incidence of a specific event type, use the Fine-Gray subdistribution hazard model (Fine & Gray, 1999).

4. You need to estimate the baseline hazard. The Cox model does not directly estimate $h_0(t)$. If you need the baseline hazard for prediction or simulation, consider parametric alternatives (Weibull, Gompertz).

Connection to Other Methods#

The Cox model sits within a broader ecosystem of methods for different outcome types:

• Logit/Probit: models whether the event occurred (binary 0/1) but ignores when. Use logit when the time dimension is not meaningful. Use Cox when timing matters. Shumway (2001) shows that switching from static logit to a hazard model significantly improves bankruptcy prediction accuracy.

• OLS on duration: regressing observed duration on covariates ignores censoring. Censored observations are treated as though their observed time is the true event time, biasing all coefficients toward zero (attenuation bias). This approach is never appropriate when censoring is present.

• Competing risks: with multiple event types, two approaches target different estimands (Fine & Gray, 1999). Cause-specific hazard models estimate how covariates affect the instantaneous event rate among those still event-free (fit with standard Cox, censoring competing events). Fine-Gray models estimate effects on the cumulative incidence function (subdistribution hazard). The right choice depends on whether you want to understand cause-specific rates or absolute risks.

• Parametric survival models: Weibull, exponential, Gompertz, and log-normal models fully specify the baseline hazard. They are more efficient than Cox when the distributional assumption is correct, but biased when it is wrong.

• Fixed effects: in panel settings with repeated spells, stratified Cox models can be used to absorb time-invariant unobserved heterogeneity (analogous to fixed effects).

Common Confusions#

For the Cox model to provide valid inference, three key assumptions must hold.

Assumption 1: Proportional Hazards#

Plain language: The effect of a covariate on the hazard rate is constant over time. If founder-CEOs have a 30% lower departure hazard than professional CEOs at year 1, they also have a 30% lower hazard at year 5, year 10, and so on.

Formally: $h(t | X_i) = h_0(t) \cdot \exp(X_i'\beta)$, where $\beta$ does not depend on $t$.

This property means the hazard ratio for any two individuals is:

The baseline hazard $h_0(t)$ cancels, and the ratio is constant over time. If the proportional hazards assumption fails — for example, if a treatment effect wears off over time — the Cox model estimates a weighted average of the time-varying effect, which may be misleading (Grambsch & Therneau, 1994).

Assumption 2: Non-Informative Censoring#

Plain language: The reason an observation is censored is unrelated to the likelihood of the event occurring. In the CEO example, this assumption means that the end of the study period, firm going private, or data unavailability is not systematically related to whether the CEO would have left soon.

Formally: $T \perp C | X$, where $T$ is the true event time and $C$ is the censoring time.

This assumption is violated, for example, if healthier patients selectively drop out of a clinical trial, or if firms with troubled CEOs are more likely to be acquired (and thus disappear from the sample).

Assumption 3: Correct Specification (for Causal Interpretation)#

Plain language: For the Cox model coefficients to have a causal interpretation, the covariates must be exogenous — the same zero conditional mean assumption required for OLS. If there are unobserved confounders that affect both the covariates and the hazard, the coefficient estimates are biased.

Formally: Conditional on $X_i$, the event times are independent of any unobserved factors that affect the hazard. This requirement is the analog of the exogeneity condition in the proportional hazards framework.

This requirement is the same exogeneity concern as in any regression model. The Cox model does not solve endogeneity — it handles censoring and the functional form for duration data.

Adjust the hazard ratio and censoring rate to see how the Kaplan-Meier survival curves diverge between treatment and control groups. Higher censoring makes the curves noisier and the log-rank test less powerful.

Explore how the proportional hazards assumption breaks down when the founder CEO advantage erodes over time:

Partial Likelihood Derivation#

Cox Regression with Diagnostics#

1# Requires: survival
2# survival: R's core package for survival analysis (Therneau & Grambsch)
3library(survival)
4
5# --- Step 1: Kaplan-Meier survival curves ---
6# Visualize non-parametric survival estimates by group before modeling
7# Surv() creates a survival object: (time, event indicator)
8km_fit <- survfit(Surv(tenure, departed) ~ founder_ceo, data = df)
9plot(km_fit, col = c("blue", "red"), lwd = 2,
10   xlab = "Years", ylab = "Survival probability",
11   main = "CEO Tenure by Type")
12legend("topright", c("Professional CEO", "Founder CEO"),
13     col = c("blue", "red"), lwd = 2)
14
15# Log-rank test: non-parametric test for equality of survival curves
16# H0: survival functions are the same across groups
17survdiff(Surv(tenure, departed) ~ founder_ceo, data = df)
18
19# --- Step 2: Cox proportional hazards regression ---
20# coxph() estimates the semiparametric Cox model: h(t|X) = h0(t) * exp(X*beta)
21# ties = "efron": Efron approximation for tied event times (more accurate than Breslow)
22cox_fit <- coxph(Surv(tenure, departed) ~ founder_ceo + firm_size +
23                 roa + industry,
24               data = df,
25               ties = "efron")
26summary(cox_fit)
27
28# Hazard ratios with 95% confidence intervals
29# exp(coef) = hazard ratio: HR > 1 means higher hazard (shorter survival)
30exp(cbind(HR = coef(cox_fit), confint(cox_fit)))
31
32# --- Step 3: Test the proportional hazards assumption ---
33# cox.zph() tests whether Schoenfeld residuals trend with time
34# H0: hazard ratios are constant over time (PH holds)
35# A significant p-value indicates PH violation for that covariate
36ph_test <- cox.zph(cox_fit)
37print(ph_test)
38plot(ph_test)  # Flat line = PH holds; trend = PH violated
39
40# --- Step 4: Predicted survival curves ---
41# Generate survival curves at specific covariate values for comparison
42# Setting covariates to median creates a "representative" profile
43newdata <- data.frame(founder_ceo = c(0, 1),
44                    firm_size = median(df$firm_size),
45                    roa = median(df$roa),
46                    industry = "Manufacturing")
47surv_pred <- survfit(cox_fit, newdata = newdata)
48plot(surv_pred, col = c("blue", "red"), lwd = 2,
49   xlab = "Years", ylab = "Survival probability")
50legend("topright", c("Professional CEO", "Founder CEO"),
51     col = c("blue", "red"), lwd = 2)

F.1 Schoenfeld Residuals (Proportional Hazards Test)#

The most important diagnostic for the Cox model. Grambsch and Therneau (1994) proposed testing the PH assumption by examining whether Schoenfeld residuals trend with time. Under the null hypothesis of proportional hazards, the scaled Schoenfeld residuals for each covariate should show no systematic pattern over time.

• Global test: tests whether any covariate violates PH (reported by cox.zph() in R, estat phtest in Stata, check_assumptions() in lifelines)
• Covariate-specific test: tests each covariate individually
• Visual inspection: plot scaled Schoenfeld residuals against time. A flat line (zero slope) supports PH; a trend suggests violation

If the PH assumption is violated for a specific covariate:

1. Stratify: coxph(Surv(time, event) ~ x1 + strata(x2), data = df) — allows different baseline hazards by strata
2. Time interaction: include $x \times \ln(t)$ or $x \times t$ to allow the effect to vary over time
3. Split the time axis: estimate separate models for early and late periods

F.2 Cox-Snell Residuals (Overall Fit)#

Cox-Snell residuals assess overall model fit. If the model is correct, these residuals should follow a unit exponential distribution. Plot the Nelson-Aalen cumulative hazard of the Cox-Snell residuals against the residuals themselves — the plot should follow a 45-degree line. Systematic departures indicate poor overall fit.

F.3 Log-Log Survival Plot#

Plot $\ln(-\ln(\hat{S}(t)))$ versus $\ln(t)$ for different groups. Under proportional hazards, these curves should be approximately parallel. Crossing or converging curves indicate PH violation. This graphical diagnostic is especially useful for categorical covariates.

F.4 Martingale Residuals (Functional Form)#

Martingale residuals assess whether continuous covariates enter the model with the correct functional form (Lin et al., 1993). Plot martingale residuals from a null model (no covariates) against each covariate. A nonlinear pattern suggests the covariate should be transformed (log, polynomial) or categorized.

F.5 Deviance Residuals (Outliers)#

Deviance residuals are a normalized transformation of martingale residuals that are more symmetrically distributed around zero. Observations with large positive deviance residuals experienced the event "too soon" relative to the model's prediction; large negative values experienced it "too late" or were censored unexpectedly early.

1library(survival)
2
3cox_fit <- coxph(Surv(tenure, departed) ~ founder_ceo + firm_size + roa,
4               data = df, ties = "efron")
5
6# F.1 Schoenfeld residuals — PH test
7ph_test <- cox.zph(cox_fit)
8print(ph_test)          # Global and per-covariate tests
9par(mfrow = c(1, 3))
10plot(ph_test)            # Scaled Schoenfeld residuals vs time
11
12# F.2 Cox-Snell residuals — overall fit
13cs_resid <- df$departed - resid(cox_fit, type = "martingale")
14surv_cs <- survfit(Surv(cs_resid, df$departed) ~ 1)
15plot(surv_cs, fun = "cumhaz",
16   xlab = "Cox-Snell residuals",
17   ylab = "Cumulative hazard",
18   main = "Cox-Snell Residual Plot")
19abline(0, 1, col = "red", lty = 2)
20
21# F.3 Log-log survival plot
22km <- survfit(Surv(tenure, departed) ~ founder_ceo, data = df)
23plot(km, fun = "cloglog",
24   xlab = "ln(time)", ylab = "ln(-ln(S(t)))",
25   col = c("blue", "red"), main = "Log-Log Plot")
26legend("topleft", c("Professional", "Founder"),
27     col = c("blue", "red"), lwd = 2)
28
29# F.4 Martingale residuals — functional form
30null_fit <- coxph(Surv(tenure, departed) ~ 1, data = df)
31mart_resid <- resid(null_fit, type = "martingale")
32plot(df$firm_size, mart_resid,
33   xlab = "Firm size", ylab = "Martingale residuals",
34   main = "Functional Form Check")
35lines(lowess(df$firm_size, mart_resid), col = "red", lwd = 2)

Hazard Ratios vs. Coefficients#

The Cox model estimates $\beta$ coefficients, but results are typically reported as hazard ratios (HR) $\exp(\beta)$:

For a continuous covariate: $HR = 1.05$ means a one-unit increase in $X$ is associated with a 5% increase in the instantaneous event rate, holding other covariates constant.

For a binary covariate: $HR = 0.65$ means the group with $X = 1$ has a hazard that is 65% of the reference group's hazard — a 35% lower event rate at every point in time.

What to Report in a Table#

A well-reported Cox regression table should include:

1. Hazard ratios (not just coefficients) with 95% confidence intervals
2. Number of subjects and number of events (not just total N)
3. Median follow-up time or total person-time at risk
4. Proportion censored
5. PH test results (Schoenfeld global test p-value)
6. Tie-handling method (Efron vs. Breslow)
7. Standard error type (robust/sandwich if used)

H.1 Concept Checks#

H.2 Guided Exercise#

H.3 Error Detective#

H.4 You Are the Referee#

• Parametric models (Weibull, exponential, Gompertz): when you need to estimate the baseline hazard — for prediction, simulation, or forecasting. More efficient than Cox when the distributional assumption is correct. Weibull is the simplest two-parameter family (can model increasing or decreasing hazards). See Therneau and Grambsch (2000) for extensions.

• Accelerated failure time (AFT) models: when the proportional hazards assumption is violated. AFT models parameterize the effect of covariates as a multiplicative shift in the time scale, not the hazard scale: $\ln(T) = X'\gamma + \sigma \varepsilon$. This formulation is analogous to OLS on log(duration) but with proper censoring handling. Common distributions: log-normal, log-logistic, Weibull.

• Discrete-time hazard models: when time is measured in discrete intervals (e.g., annual employment spells, quarterly observations) or when you have interval-censored data. These models are logit or complementary log-log specifications applied to person-period data (Singer & Willett, 2003).

• Competing risks (Fine-Gray): when multiple event types exist and you want to estimate effects on the cumulative incidence of a specific event, rather than the cause-specific hazard (Fine & Gray, 1999). The Fine-Gray model accounts for the fact that experiencing a competing event removes the subject from risk for the event of interest.

• Frailty models (random effects): when you have clustered data (e.g., employees within firms, patients within hospitals) and want to account for unobserved heterogeneity. A frailty term adds a random effect to the Cox model (Therneau & Grambsch, 2000).