Survival Analysis in Epidemiology & Causal Inference Frameworks

Let $T \ge 0$ denote a non-negative continuous random variable representing the time from a well-defined baseline origin (such as date of cancer diagnosis, cardiac surgery, or study enrollment) until the occurrence of a pre-specified clinical event of interest (such as death, relapse, or recovery). ### The Four Fundamental Mathematical Representations The probability distribution of $T$ is uniquely characterized by four mathematically equivalent functions: 1. **Cumulative Distribution Function** ($F(t)$): The probability that the event occurs on or before time $t$: $$F(t) = P(T \le t) = \int_0^t f(u) \, du$$ 2. **Survival Function** ($S(t)$): The probability that an individual survives longer than time $t$ (remains event-free): $$S(t) = P(T > t) = 1 - F(t) = \int_t^\infty f(u) \, du$$ Properties: $S(0) = 1$, $\lim_{t \to \infty} S(t) = 0$, and $S(t)$ is non-increasing. 3. **Probability Density Function** ($f(t)$): The unconditional rate of event occurrence per unit time: $$f(t) = \frac{dF(t)}{dt} = -\frac{dS(t)}{dt}$$ 4. **Hazard Function** ($h(t)$ or $\lambda(t)$): The instantaneous rate of event occurrence at time $t$, conditional upon having survived up to time $t$: $$h(t) = \lim_{\Delta t \to 0} \frac{P(t \le T < t + \Delta t \mid T \ge t)}{\Delta t} = \frac{f(t)}{S(t)} = -\frac{1}{S(t)} \frac{dS(t)}{dt} = -\frac{d}{dt} \ln S(t)$$ Integrating the hazard function yields the **Cumulative Hazard Function** ($H(t)$): $$H(t) = \int_0^t h(u) \, du = -\ln S(t)$$ Exponentiating yields the universal mathematical link connecting survival and hazard: $$S(t) = \exp(-H(t)) = \exp\left( -\int_0^t h(u) \, du \right)$$ ### Censoring and Truncation Mechanisms Survival data are almost never fully observed; they are subject to incomplete observation mechanisms: - **Right Censoring**: The most prevalent mechanism. The true event time $T_i$ is known only to exceed a recorded censoring time $C_i$ ($T_i > C_i$). Right censoring occurs when a participant reaches the administrative study end date without experiencing the event, migrates away, or dies from an unrelated competing cause. The observed data pair is $(Y_i, \delta_i)$, where: $$Y_i = \min(T_i, C_i), \qquad \delta_i = I(T_i \le C_i) \in \{0, 1\}$$ - **Left Truncation (Delayed Entry)**: Occurs when individuals enter observation at an age or time $L_i > 0$ subsequent to the natural origin. Individuals who experience the event prior to $L_i$ are never observed, introducing severe survival selection unless conditioned upon via staggered risk sets ($T_i > L_i$). - **Informative vs Independent Censoring**: Standard survival methods assume **independent (non-informative) censoring**: conditional on observed covariates $\mathbf{X}$, the event time $T$ and censoring time $C$ are statistically independent ($T \perp C \mid \mathbf{X}$).

When the mathematical functional form of the hazard rate is unknown, imposing parametric distributions (such as the exponential or Weibull) risks severe model misspecification. In 1958, Edward L. Kaplan and Paul Meier formulated the non-parametric Product-Limit Estimator.

Mathematical Derivation of the Kaplan-Meier Estimator

Suppose that in a sample of $n$ individuals, distinct failure events occur at $k$ ordered calendar times:

At each distinct event time $t_j$:

• $d_j$: The number of individuals who experience the event at exactly time $t_j$.
• $Y_j$: The number of individuals who are alive, uncensored, and under active observation immediately prior to time $t_j$ (the risk set at $t_j$).
• $c_j$: The number of individuals censored in the half-open interval $[t_j, t_{j+1})$.

The conditional probability of surviving past $t_j$, given survival up to $t_j$, is non-parametrically estimated as:

By the chain rule of conditional probabilities, the overall probability of surviving past any arbitrary time $t$ is the product of all conditional survival probabilities up to time $t$:

The Kaplan-Meier curve is a step function that starts at $\widehat{S}(0) = 1.0$ and steps downward exclusively at each observed failure time $t_j$, with vertical step magnitude proportional to $d_j / Y_j$. Censored observations do not trigger steps; instead, they reduce the size of the risk set $Y_{j+1}$ for subsequent intervals.

Greenwood's Formula for Sampling Variance

Major Greenwood (1926) derived the asymptotic sampling variance of the product-limit survival estimator using the delta method on the log-transformed survival function $\ln \widehat{S}(t) = \sum_{j : t_j \le t} \ln(1 - d_j / Y_j)$. Treating each conditional failure count as an independent binomial realization:

Applying the delta method again to $\exp(\ln \widehat{S}(t))$ yields Greenwood’s Formula:

The standard Greenwood 95% Wald confidence interval is $\widehat{S}(t) \pm 1.96 \sqrt{\widehat{\text{Var}}(\widehat{S}(t))}$. To prevent confidence intervals from overshooting the unit boundaries $[0, 1]$, modern software applies the log-log transformation $g(S) = \ln(-\ln S)$.

In randomized oncology trials and comparative cohort studies, a central task is testing whether two or more survival curves differ significantly across the entire follow-up history ($H_0: S_1(t) = S_2(t) \quad \forall t$). ### The Mantel-Cox Log-Rank Test Nathan Mantel (1966) and David Cox (1972) formulated the **Log-Rank Test** by constructing a $2 \times 2$ contingency table at each distinct failure time $t_j$: Under the null hypothesis of equal hazard rates ($H_0: h_1(t) = h_2(t)$), the number of events in Group 1 at time $t_j$ follows a hypergeometric distribution conditional upon marginal totals. The expected number of events in Group 1 is: $$E_{1j} = \mathbb{E}[d_{1j} \mid H_0] = d_j \frac{Y_{1j}}{Y_j}$$ The conditional hypergeometric variance is: $$V_{1j} = \text{Var}(d_{1j} \mid H_0) = \frac{Y_{1j} Y_{2j} d_j (Y_j - d_j)}{Y_j^2 (Y_j - 1)}$$ Summing over all distinct failure times $t_1, \dots, t_k$: $$O_1 = \sum_{j=1}^k d_{1j}, \qquad E_1 = \sum_{j=1}^k E_{1j}, \qquad V_1 = \sum_{j=1}^k V_{1j}$$ The **Log-Rank Test Statistic** is: $$\chi^2_{\text{log-rank}} = \frac{(O_1 - E_1)^2}{V_1}$$ Under $H_0$, $\chi^2_{\text{log-rank}}$ asymptotically follows a chi-square distribution with $1$ degree of freedom. An approximate estimate of the Hazard Ratio is obtained as: $$\widehat{HR} \approx \frac{O_1 / E_1}{O_2 / E_2}$$ The log-rank test places equal weight on events across all time points, making it the most statistically powerful test under proportional hazards alternatives.

While the log-rank test evaluates binary or categorical group differences, it cannot adjust for continuous confounders or evaluate multiple prognostic variables simultaneously. In 1972, Sir David Cox revolutionized survival analysis by formulating the Proportional Hazards Regression Model.

Model Specification

The Cox model postulates that the hazard of disease for an individual with covariate vector $\mathbf{x}_i = (x_{i1}, \dots, x_{ip})^T$ is the product of an unspecified non-parametric baseline hazard and an exponential parametric function:

where:

• $h_0(t)$: The baseline hazard function, representing the hazard for an individual whose covariates are all zero ($\mathbf{x} = \mathbf{0}$). It is left completely unspecified, requiring no parametric distributional assumptions!
• $\exp(\mathbf{x}_i^T \boldsymbol{\beta})$: The relative hazard (risk multiplier), which is strictly parametric and time-invariant.

The model is semi-parametric because it pairs an infinite-dimensional non-parametric baseline hazard $h_0(t)$ with finite-dimensional regression coefficients $\boldsymbol{\beta}$.

The Proportional Hazards Property

Consider two individuals with covariate vectors $\mathbf{x}_1$ and $\mathbf{x}_0$. The ratio of their hazards is:

Notice that the baseline hazard $h_0(t)$ cancels out completely! The Hazard Ratio ($HR$) is strictly constant over time: the hazard curves for any two individuals remain strictly proportional across the entire observational horizon.

Maximum Partial Likelihood Estimation

Cox’s mathematical genius lay in demonstrating that $\boldsymbol{\beta}$ can be estimated without knowing or estimating $h_0(t)$ through the method of Partial Likelihood.

Suppose failure events occur at $k$ distinct ordered times $t_1 < \dots < t_k$. At failure time $t_j$, let $\mathcal{R}(t_j)$ denote the risk set of all individuals alive and under observation immediately prior to $t_j$. Conditional on an event occurring at $t_j$, the probability that the event occurs to the specific individual $j$ who actually failed is:

Canceling $h_0(t_j)$:

The Partial Likelihood is the product across all $k$ failure events:

Maximizing the partial log-likelihood $\ell_{\text{partial}}(\boldsymbol{\beta})$ via Newton-Raphson yields estimators that are asymptotically unbiased, normally distributed, and achieve full semiparametric efficiency.

The validity of the Cox model depends on the strict satisfaction of the Proportional Hazards (PH) assumption. If the true hazard ratio varies over time (for example, if a surgical intervention carries high immediate perioperative mortality but strong long-term survival benefit), the baseline model fails.

Schoenfeld Residuals

David Schoenfeld (1982) formulated residuals defined exclusively at observed failure times. For covariate $k$ and failure event $j$:

The Schoenfeld residual represents the difference between the actual covariate value of the failing individual and the risk-weighted average of that covariate across the active risk set at that moment.

Under the proportional hazards assumption:

• Schoenfeld residuals have an expected value of zero ($\mathbb{E}[r_{jk}] = 0$) across all failure times.
• A plot of scaled Schoenfeld residuals against time (or log-time) should trace a flat, horizontal line with slope zero.
• Grambsch and Therneau (1994) formulated a formal score test of non-zero slope: if $p < 0.05$, the proportional hazards assumption is rejected.

When the proportional hazards assumption is violated, several analytical remedies exist in contemporary survival epidemiology:

1. Stratified Cox Models: If a categorical covariate violates proportional hazards, one can stratify on that variable. In a stratified Cox model, the baseline hazard function is allowed to vary arbitrarily across strata $s \in \{1, \dots, S\}$, while the regression coefficients $\boldsymbol{\beta}$ for the other covariates remain common across all strata:

The total partial log-likelihood is simply the sum of the stratum-specific partial log-likelihoods $\ell(\boldsymbol{\beta}) = \sum_{s=1}^S \ell_s(\boldsymbol{\beta})$. This strategy completely eliminates confounding or non-proportionality by the stratifying factor without requiring any parametric modeling of its time-dependent dynamics.

2. Time-Varying Coefficients and Step Functions: If an exposure of primary interest exhibits non-proportional hazards, one can model the regression coefficient as an explicit function of time $\beta(t) = \beta + \theta g(t)$, where common choices for $g(t)$ include $\ln(t)$, $t$, or piecewise step functions over clinically meaningful follow-up intervals ($0-1$ year, $1-5$ years, and $>5$ years). Testing whether $\theta = 0$ provides a direct likelihood-ratio or Wald test of proportionality.

3. Restricted Mean Survival Time (RMST): When proportional hazards fails entirely (for instance, when survival curves cross due to early surgical mortality followed by late oncologic curative cure), the hazard ratio becomes an unstable, follow-up-dependent summary. In such settings, the difference in Restricted Mean Survival Time up to a fixed milestone $\tau$, defined as $\Delta \text{RMST}(\tau) = \int_0^{\tau} S_1(t)\,dt - \int_0^{\tau} S_0(t)\,dt$, provides an unconfounded, model-free metric of absolute life expectancy gained.

When the PH assumption is violated, epidemiologists resolve the non-proportionality by:

1. Including time-varying covariates ($x_j \times \ln t$).
2. Fitting Stratified Cox Models, allowing the baseline hazard $h_{0g}(t)$ to differ across strata.
3. Fitting Accelerated Failure Time (AFT) models.

Observational epidemiology has increasingly embraced the formal mathematical language of Causal Inference, rooted in Donald Rubin’s Potential Outcomes model and Judea Pearl’s structural causal frameworks.

The Three Identifiability Conditions

To identify the Average Causal Effect ($ACE = \mathbb{E}[Y(1)] - \mathbb{E}[Y(0)]$) from observational data, three foundational conditions must hold:

1. Consistency: An individual's observed outcome under observed exposure $A = a$ equals their potential counterfactual outcome under that exposure level: $Y = Y(a)$. This requires well-defined interventions with no hidden variations of treatment.

2. Conditional Exchangeability (No Unmeasured Confounding): Conditional on a measured set of baseline covariates $\mathbf{L}$, the potential outcomes are statistically independent of the treatment assignment:

3. Positivity: Every individual in the population has a non-zero probability of receiving every treatment level:

Marginal Structural Models and Inverse Probability of Treatment Weighting (IPTW)

When an observational study involves time-varying exposures and time-varying confounders that are themselves affected by prior exposure (such as CD4 counts in longitudinal HIV studies), standard multivariable regression models fail catastrophically: conditioning on the time-varying confounder introduces collider bias, while omitting it leaves confounding unadjusted.

James Robins (2000) resolved this dilemma through Inverse Probability of Treatment Weighting (IPTW). The propensity score is the conditional probability of receiving treatment given baseline covariates:

Each participant is assigned a statistical weight equal to the inverse of the probability of receiving the treatment they actually received:

Weighting creates a pseudo-population in which the association between confounders $\mathbf{L}$ and exposure $A$ is completely severed. In this re-weighted pseudo-population, a simple unadjusted regression of outcome $Y$ on exposure $A$ (a Marginal Structural Model) yields a mathematically unconfounded, asymptotically unbiased estimate of the true average causal effect.

While the Cox proportional hazards model is the most widely applied survival model in medical literature, two critical alternative frameworks are essential for comprehensive epidemiological mastery:

Accelerated Failure Time (AFT) Models

Under the Accelerated Failure Time (AFT) formulation, the effect of an exposure is modeled not as a multiplicative hazard shift, but as a direct acceleration or deceleration of the biological survival time clock:

where $\varepsilon_i$ is a random error disturbance with a pre-specified distribution (e.g., standard extreme value yielding the Weibull AFT, logistic yielding log-logistic AFT, or Gaussian yielding log-normal AFT). The exponentiated coefficient $\exp(\beta_1)$ represents the Acceleration Factor ($\gamma$):

• If $\gamma > 1.0$: The exposure prolongs time-to-event (decelerates aging/progression); the median survival time is multiplied by $\gamma$.
• If $\gamma < 1.0$: The exposure accelerates time-to-failure.

Unlike hazard ratios, acceleration factors are easily communicated to non-statistical clinical audiences as percentage changes in survival time.

Competing Risks and Fine-Gray Subdistribution Hazards

In aging populations, individuals are simultaneously at risk of multiple mutually exclusive terminating events (e.g., death from cardiovascular disease, death from prostate cancer, or fatal motor vehicle trauma). Treating competing deaths as simple independent censoring violates non-informative censoring, producing severe upward bias in the cumulative incidence of the primary event.

Robert Gray (1988) and Jason Fine and Robert Gray (1999) resolved this by modeling the Cumulative Incidence Function ($CIF_k(t) = P(T \le t, \text{Cause} = k)$) directly through the Subdistribution Hazard:

In the Fine-Gray subdistribution risk set, an individual who dies from a competing cause (e.g., heart disease at month 12) remains in the denominator risk set indefinitely for the cancer outcome, preventing the false inflation of cancer cumulative incidence. Fitting proportional subdistribution hazards regression models ensures mathematically sound risk prediction in the presence of competing terminal events.

In Section 8.4, we evaluated the Basic Reproduction Number $\mathcal{R}_0 = \beta / \gamma$ for the simple classical three-compartment SIR model. However, real infectious disease pathogens involve complex multi-compartment transmission pathways—such as latent incubation states ($E$ in SEIR models), asymptomatic transmission compartments ($A$), hospitalized isolation states ($H$), and cross-species vector compartments (mosquitoes in dengue and malaria).

For high-dimensional, heterogeneous compartmental transmission models, calculating $\mathcal{R}_0$ via scalar algebraic manipulation is prone to severe errors. In 1990, Odo Diekmann, Hans Heesterbeek, and Johan Metz introduced the rigorous mathematical standard: The Next Generation Matrix (NGM).

1. Classification of Compartments and the Disease-Free Equilibrium

Let the continuous population be divided into $n$ distinct epidemiological compartments $\mathbf{x} = [x_1, x_2, \dots, x_n]^T$. Sort the compartments such that the first $m$ compartments ($x_1, \dots, x_m$) are infected compartments (individuals currently harboring pathogen), while the remaining $n - m$ compartments are uninfected (e.g., susceptible $S$, recovered $R$).

The disease system evolves according to the coupled ordinary differential equations:

where:

• $\mathcal{F}_i(\mathbf{x})$: The rate of appearance of new infections into compartment $i$.
• $\mathcal{V}_i(\mathbf{x}) = \mathcal{V}_i^-(\mathbf{x}) - \mathcal{V}_i^+(\mathbf{x})$: The net rate of transfer of individuals out of compartment $i$ via all other biological transitions (recovery, progression to subsequent disease stages, natural mortality, and disease-induced mortality).

Let $\mathbf{x}_0 = [0, \dots, 0, x_{m+1}^*, \dots, x_n^*]^T$ denote the Disease-Free Equilibrium (DFE), where all infected states are identically zero.

2. Jacobian Linearization and the Next Generation Operator

Linearizing the infected subsystem about the disease-free equilibrium $\mathbf{x}_0$ yields two $m \times m$ Jacobian matrices:

The matrices admit transparent biological interpretations:

• Element $F_{ij}$ represents the rate at which an infected individual in compartment $j$ generates new infections entering compartment $i$.
• Element $(V^{-1})_{jk}$ represents the average expected duration of time that an individual who initially enters infected compartment $k$ will spend residing in infected compartment $j$ throughout their entire infective lifespan.

The product matrix is the Next Generation Matrix:

Element $K_{ik} = (\mathbf{F} \mathbf{V}^{-1})_{ik}$ represents the expected total number of secondary infections in compartment $i$ generated throughout the entire infectious lifetime of a single primary infected individual originally introduced into compartment $k$.

3. Theorem on $\mathcal{R}_0$ and Epidemic Stability

Diekmann and Heesterbeek proved the foundational stability theorem: Theorem: The Basic Reproduction Number $\mathcal{R}_0$ is exactly the Spectral Radius (dominant eigenvalue in absolute value) of the Next Generation Matrix:

Furthermore:

• If $\mathcal{R}_0 < 1$, the disease-free equilibrium $\mathbf{x}_0$ is locally asymptotically stable: any small viral introduction decays exponentially to extinction.
• If $\mathcal{R}_0 > 1$, the disease-free equilibrium is unstable, guaranteeing that a pathogen introduction triggers an exponential epidemic wave.

By converting complex multi-stage differential systems into eigenvalue problems, the Next Generation Matrix provides mathematical epidemiologists with an analytical tool for evaluating vaccine coverage thresholds, quarantine protocols, and vector control programs.