Survival Analysis Censoring Kaplan-Meier Estimator Cox Model

# Survival Analysis: Censoring, Kaplan-Meier Estimator & Cox Model

## Introduction & Motivation

Survival analysis models time-to-event data (duration until failure/death). Censoring: incomplete observation. Kaplan-Meier (product-limit) nonparametrically estimates survival curve; handles censoring via conditional probability. Cox proportional hazards regresses on covariates while accounting for censoring. Critical for medical trials, reliability engineering, customer churn.

Motivation: Standard regression assumes complete observation. Censoring (dropout, end-of-study) common in practice. Ignoring censoring biases estimates.

Applications: Medical trials (time-to-recovery), reliability (time-to-failure), customer lifetime value, unemployment duration.

---

## Core Concepts & Theory

### Survival Function

S(t) = P(T > t). Probability of surviving past time t.

### Hazard Function

h(t) = lim_{Δt→0} P(t < T ≤ t+Δt | T > t) / Δt. Instantaneous failure rate.

### Censoring

Right censoring: observation time c < T (unknown T). Kaplan-Meier handles via product-limit.

---

## Mathematical Formulation

Kaplan-Meier survival estimate:
$$S(t) = \prod_{t_i \leq t} \left(1 - \frac{d_i}{n_i} ight)$$

where d_i = events at t_i, n_i = at-risk at t_i.

Cox partial likelihood:
$$\ell(\beta) = \sum_{i: \delta_i = 1} \left[\beta^T x_i - \log \sum_{j \in R(t_i)} \exp(\beta^T x_j) ight]$$

where R(t_i) = risk set at time t_i.

Hazard ratio:
$$ ext{HR} = \exp(\beta)$$

---

## Advanced Theory & Extensions

### Competing Risks

Multiple event types; cumulative incidence accounts for competing events.

### Time-Dependent Covariates

Covariates change over time; stratified or time-interaction models.

### Frailty Models

Random intercept for hierarchical data; shared frailty for clusters.

---

## Computational Considerations

Kaplan-Meier: O(n log n) sorting + O(n) computation.

Cox partial likelihood: O(iterations × n × p) via Newton-Raphson; scales poorly with p.

Stratified Cox: Condition on strata; handles non-proportional hazards locally.

---

## Practical Implementation Strategies

### Checking Proportional Hazards

Schoenfeld residuals: plot vs. time; should be flat.

### Handling Ties

Breslow, Efron approximations in partial likelihood.

### Confidence Intervals

95% CI: exp(β ± 1.96 * SE(β)). Large-sample normal approximation.

---

## Benchmark Datasets & Evaluation

WHAS500: Worcester Heart Attack Study; 500 patients, 1-year mortality.

Veteran's Lung Cancer: 137 patients, treatment comparison.

Metrics: C-index (concordance), log-rank test (curve comparison), AUC at t.

---

## Key Challenges & Limitations

### Non-Proportional Hazards

Cox assumes constant hazard ratio; violated if lines cross. Use stratified models.

### Sparse Events

Few events limit precision; use Firth's bias-corrected partial likelihood.

### Competing Risks

Ignoring competitors overestimates event probability.

---

## Hyperparameter Tuning

Smoothing (smoothed baseline hazard): Bandwidth via cross-validation.

Frailty variance prior: Weakly informative to avoid overfitting.

Stratification: Balance bias-variance.

---

## Real-World Applications & Case Studies

Clinical Trials: FDA requires survival curves; intent-to-treat with censoring.

Organ Transplant: Waiting list survival; right-censoring from transplant/censoring.

Insurance: Mortality tables; age-stratified Cox model.

---

## Integration with Other Methods

Survival + Machine Learning → Random Survival Forests, gradient-boosted survival.

Survival + Neural Networks → DeepHit (multi-task learning for competing risks).

---

## Summary & Key Takeaways

Survival analysis quantifies time-to-event distributions, accounting for censoring via Kaplan-Meier (nonparametric) or Cox regression (semiparametric).

Principles:
1. Censoring is missing data; must be handled correctly.
2. Kaplan-Meier uses product-limit formula.
3. Cox assumes proportional hazards; verify via residuals.
4. Hazard ratio exp(β) is interpretable effect size.
5. Log-rank test compares survival curves.

---

---

## Appendix: Practical Labs

### Lab 1: Kaplan-Meier Estimator

import numpy as np
import pandas as pd

def kaplan_meier(T, event):
 """
 T: time-to-event
 event: 1 = event, 0 = censored
 """
 data = pd.DataFrame({'T': T, 'event': event})
 data = data.sort_values('T')
 
 S = 1.0
 survival_times = []
 survival_probs = []
 
 for t in data['T'].unique():
 subset = data[data['T'] <= t]
 at_risk = len(data[data['T'] >= t])
 events = len(subset[subset['event'] == 1])
 
 if at_risk > 0:
 S *= (1 - events / at_risk)
 survival_times.append(t)
 survival_probs.append(S)
 
 return np.array(survival_times), np.array(survival_probs)

# Data
T = np.array([5, 6, 6, 7, 8, 9, 10])
event = np.array([1, 1, 0, 1, 0, 1, 1])

times, probs = kaplan_meier(T, event)
print(f"Survival times: {times}")
print(f"Survival probs: {probs}")
assert len(times) == len(probs), "Should match"
assert (probs >= 0).all() and (probs <= 1).all(), "Probs should be in [0, 1]"
print("✓ Kaplan-Meier working")

if __name__ == "__main__":
 print("Lab 1: Kaplan-Meier - PASSED")

### Lab 2: Log-Rank Test

import numpy as np
from scipy.stats import chi2

def log_rank_test(T1, event1, T2, event2):
 """Compare two survival curves"""
 # Merge data
 all_times = np.concatenate([T1, T2])
 all_events = np.concatenate([event1, event2])
 group = np.concatenate([np.zeros_like(T1), np.ones_like(T2)])
 
 # Compute test statistic
 O1 = 0 # Observed events in group 1
 E1 = 0 # Expected events
 V = 0 # Variance
 
 for t in np.unique(all_times):
 mask = all_times == t
 n1 = (group[mask] == 0).sum()
 n2 = (group[mask] == 1).sum()
 d = all_events[mask].sum()
 
 if n1 + n2 > 0:
 O1 += (all_events[mask] * (group[mask] == 0)).sum()
 E1 += d * n1 / (n1 + n2)
 V += d * (n1 + n2 - d) * n1 * n2 / ((n1 + n2)**2 * (n1 + n2 - 1) + 1e-8)
 
 # Chi-squared statistic
 Z = (O1 - E1)**2 / (V + 1e-8)
 pvalue = 1 - chi2.cdf(Z, df=1)
 
 return Z, pvalue

T1 = np.array([1, 3, 5, 6, 8])
event1 = np.array([1, 1, 0, 1, 1])
T2 = np.array([2, 4, 6, 7, 9])
event2 = np.array([1, 1, 1, 0, 1])

Z, pvalue = log_rank_test(T1, event1, T2, event2)
print(f"Log-rank test: Z={Z:.4f}, p-value={pvalue:.4f}")
assert 0 <= pvalue <= 1, "P-value should be in [0, 1]"
print("✓ Log-rank test working")

if __name__ == "__main__":
 print("Lab 2: Log-Rank - PASSED")

### Lab 3: Cox Proportional Hazards

import numpy as np
from scipy.optimize import minimize

def cox_partial_likelihood(beta, X, T, event):
 """Negative partial likelihood for Cox regression"""
 n = len(T)
 risk_score = np.exp(X @ beta)
 
 # Sort by time
 sort_idx = np.argsort(T)
 T_sorted = T[sort_idx]
 event_sorted = event[sort_idx]
 risk_sorted = risk_score[sort_idx]
 
 nll = 0
 for i in range(n):
 if event_sorted[i] == 1:
 risk_set = risk_sorted[i:]
 nll -= np.log(risk_sorted[i] / np.sum(risk_set) + 1e-8)
 
 return nll / n

# Data
np.random.seed(42)
n = 50
X = np.column_stack([np.random.randn(n), np.random.randn(n)])
beta_true = np.array([0.5, -0.3])
T = np.random.exponential(scale=10, size=n)
event = np.random.binomial(1, 0.7, n)

# Fit
beta_init = np.zeros(2)
result = minimize(cox_partial_likelihood, beta_init, args=(X, T, event), method='BFGS')
beta_hat = result.x

print(f"Estimated beta: {beta_hat}")
assert len(beta_hat) == 2, "Should estimate 2 coefficients"
assert np.isfinite(beta_hat).all(), "Should be finite"
print("✓ Cox model working")

if __name__ == "__main__":
 print("Lab 3: Cox - PASSED")

### Lab 4: Hazard Ratio Interpretation

import numpy as np

def hazard_ratio(beta):
 """Hazard ratio from log-scale coefficient"""
 return np.exp(beta)

def hazard_ratio_ci(beta, se_beta, ci=0.95):
 """95% CI for hazard ratio"""
 z_crit = 1.96
 ci_log = [beta - z_crit * se_beta, beta + z_crit * se_beta]
 return np.exp(ci_log[0]), np.exp(ci_log[1])

# Example: treatment effect
beta_treatment = 0.3
se_treatment = 0.1

hr = hazard_ratio(beta_treatment)
ci_low, ci_high = hazard_ratio_ci(beta_treatment, se_treatment)

print(f"Hazard ratio: {hr:.4f}, 95% CI: [{ci_low:.4f}, {ci_high:.4f}]")
assert hr > 0, "HR should be positive"
assert ci_low < hr < ci_high, "HR should be in CI"
assert ci_low > 0, "CI should be positive"
print("✓ Hazard ratio working")

if __name__ == "__main__":
 print("Lab 4: HR - PASSED")