Survival Analysis Fundamentals

This page covers the statistical theory behind the Survival Analysis tabs. See that page for usage instructions.

Time-to-Event Data and Censoring

Survival analysis is a set of methods for analyzing time until an event occurs. Despite the name "survival," the event need not be death — it can be machine failure, customer churn, time to recidivism, or any event of interest.

The defining feature of survival data is censoring. Subjects who did not experience the event during the observation period (e.g., patients still alive at the end of a clinical trial, or patients lost to follow-up) carry only incomplete information: "the event had not occurred by at least this time."

Simply excluding censored observations biases the analysis toward subjects who experienced the event sooner, underestimating survival times. Treating censored observations as "no event" overstates survival times since the true event time is unknown. Survival analysis methods are designed to handle censoring properly, provided that censoring is independent of event occurrence (non-informative censoring). When censoring is related to the likelihood of the event — for example, when patients drop out due to worsening side effects — the KM estimator and Cox model estimates become biased. MIDAS handles right censoring only (censoring due to end of observation or loss to follow-up). Left censoring and interval censoring are not supported.

Why Ordinary Regression Fails

Without censoring, survival times could be analyzed as a response variable in ordinary regression. But censored data provides inequality information — "the true value is at least as large as the observed value" — and the usual residual ( $y_i - \hat{y}_i$ ) cannot be defined. Survival analysis incorporates this inequality into the likelihood function, correctly accounting for censoring.

Survival Function and Hazard Function

The distribution of survival time $T$ is characterized by two functions.

The survival function $S(t) = P(T > t)$ is the probability of not having experienced the event by time $t$ . It starts at $S(0) = 1$ and decreases monotonically over time.

The hazard function $h(t)$ is the instantaneous rate of event occurrence at time $t$ , given survival up to that point:

h(t) = \lim_{\Delta t \to 0} \frac{P(t \le T < t + \Delta t \mid T \ge t)}{\Delta t}

The hazard is a rate (per unit time), not a probability, so it can exceed 1. The survival and hazard functions are related by $S(t) = \exp\left(-\int_0^t h(u)\,du\right)$ ; knowing one determines the other.

Kaplan-Meier Estimator

The Kaplan-Meier estimator is a nonparametric estimator of the survival function. It makes no distributional assumptions, estimating $S(t)$ directly from observed event times.

Let the distinct event times be $t_1 < t_2 < \cdots < t_k$ , with $n_i$ subjects at risk and $d_i$ events at each time $t_i$ :

\hat{S}(t) = \prod_{t_i \le t} \left(1 - \frac{d_i}{n_i}\right)

This cumulatively multiplies the "survival fraction" at each event time. Censoring is reflected through changes in the risk set: when a subject is censored, they leave the risk set but are not counted as an event.

Under non-informative censoring, the KM estimator is a consistent estimator of $S(t)$ . The variance is estimated using Greenwood's formula, derived via the delta method:

\widehat{\operatorname{Var}}[\hat{S}(t)] = \hat{S}(t)^2 \sum_{t_i \le t} \frac{d_i}{n_i(n_i - d_i)}

Confidence intervals are constructed from this variance. MIDAS uses the log transformation method by default, computing $\exp(\log \hat{S}(t) \pm z \cdot \text{SE} / \hat{S}(t))$ . The log transformation prevents the interval from falling outside $[0, 1]$ .

Log-rank Test

The log-rank test is a nonparametric test for whether two or more groups have equal hazard functions.

At each event time $t_i$ , the observed event count $d_{ij}$ and expected event count $e_{ij} = n_{ij} \cdot d_i / n_i$ are computed for group $j$ , where $n_{ij}$ is the size of group $j$ 's risk set at time $t_i$ . Under the null hypothesis, events at each time point are expected to be distributed according to each group's share of the risk set. This allocation follows a hypergeometric distribution, from which the variance is also derived. For two groups, since $O_1 + O_2 = \sum d_i$ and $E_1 + E_2 = \sum d_i$ , it follows that $O_1 - E_1 = -(O_2 - E_2)$ , so only one group's information is needed. The test statistic is:

\chi^2 = \frac{(O_1 - E_1)^2}{V_1}

where $O_1 = \sum_i d_{1i}$ (total observed events in group 1), $E_1 = \sum_i e_{1i}$ (total expected events), and $V_1 = \sum_i \frac{n_{1i} \cdot n_{2i} \cdot d_i (n_i - d_i)}{n_i^2(n_i - 1)}$ (variance). The denominator is the variance derived from the hypergeometric distribution, not the expected value $E_1$ . For three or more groups, this extends to the quadratic form $\chi^2 = \mathbf{U}' \mathbf{V}^{-1} \mathbf{U}$ using the variance-covariance matrix. Under the null hypothesis, this statistic approximately follows a chi-squared distribution with $g - 1$ degrees of freedom ( $g$ = number of groups).

The log-rank test is a member of the weighted log-rank test family. It has maximum power compared to other weighted variants (such as Wilcoxon-type tests) when the hazard ratio is constant over time (under the proportional hazards assumption). Power decreases when survival curves cross.

Cox Proportional Hazards Model

Model Formulation

The Cox (1972) proportional hazards model is a semiparametric model that estimates the effect of covariates on hazard:

h(t \mid X) = h_0(t) \exp(\beta_1 X_1 + \beta_2 X_2 + \cdots + \beta_p X_p)

$h_0(t)$ is the baseline hazard (the hazard when all covariates are zero), and $\exp(\beta_j)$ is the hazard ratio for a one-unit increase in covariate $X_j$ .

It is called "semiparametric" because $\beta$ is estimated parametrically, but no functional form is specified for $h_0(t)$ . This removes the need to assume a distribution for the baseline hazard. After estimating $\beta$ , the baseline hazard $h_0(t)$ can be estimated nonparametrically using the Breslow estimator, which in turn allows computing the survival function $S(t|X)$ for specific covariate values. MIDAS outputs the baseline cumulative hazard $H_0(t)$ and the adjusted survival curve $S(t|X)$ at user-specified covariate values.

The covariates $X$ in this model are fixed values for each subject throughout the observation period. Handling covariates that change over time (time-varying covariates) requires extensions that MIDAS does not currently support.

Proportional Hazards Assumption

The core assumption is that covariate effects are constant over time. That is, the hazard ratio $h(t \mid X_1) / h(t \mid X_2) = \exp(\beta'(X_1 - X_2))$ does not depend on $t$ .

When this assumption is violated (e.g., a treatment effect that fades over time), the estimated $\beta$ represents a weighted average of the time-varying effect, with weights that depend on the risk set composition and baseline hazard, making interpretation difficult (Struthers & Kalbfleisch, 1986).

Schoenfeld Residuals

Schoenfeld residuals are used to assess the proportional hazards assumption. They are defined at each event time $t_{(i)}$ for each covariate $j$ :

r_{ij} = X_{ij} - \bar{X}_j(t_{(i)})

$X_{ij}$ is the value of covariate $j$ for the subject who experienced the event. $\bar{X}_j(t_{(i)})$ is the weighted mean of covariate $j$ over the risk set, defined as:

\bar{X}_j(t_{(i)}) = \frac{\sum_{k \in \mathcal{R}(t_{(i)})} X_{kj} \exp(X_k'\hat\beta)}{\sum_{k \in \mathcal{R}(t_{(i)})} \exp(X_k'\hat\beta)}

Here $k$ ranges over subjects in the risk set $\mathcal{R}(t_{(i)})$ , $X_k$ is the covariate vector for subject $k$ , and $X_k'\hat\beta$ is the linear predictor. The weight $\exp(X_k'\hat\beta)$ is subject-specific — subjects with higher hazard contribute more — and is the same for all covariates $j$ .

The sum of Schoenfeld residuals equals the score function (the gradient of the log partial likelihood with respect to $\beta$ ) (Schoenfeld, 1982). At the MLE $\hat\beta$ the score function is theoretically zero, so the sum of residuals is also zero. In practice, Newton-Raphson terminates after finitely many iterations, so the sum is zero only up to convergence tolerance.

Scaled Schoenfeld residuals adjust the raw residuals by the variance-covariance matrix:

r^*_i = d \cdot \hat{V} \, r_i + \hat\beta

where $d$ is the total number of events, $\hat{V} = \hat{I}(\hat\beta)^{-1}$ is the estimated variance-covariance matrix, and $r_i$ and $r^*_i$ are the raw and scaled residual vectors at event time $i$ . The $j$ -th component of $r^*_i$ can be interpreted as an estimate of $\beta_j(t_{(i)})$ . Under proportional hazards, $r^*_i$ shows no systematic trend over time (Grambsch & Therneau, 1994).

MIDAS displays the following diagnostics (usage):

Grambsch-Therneau test: Tests the association between scaled Schoenfeld residuals and time. Reports per-covariate tests and a global test
Scaled Schoenfeld residual plots: Plots $r^*_{ij}$ against time with a LOESS smooth
log(-log(S(t))) plot: Plots Kaplan-Meier estimates as $\log(-\log(S(t)))$ versus $\log(t)$ by group. Under proportional hazards, the curves should be approximately parallel

Partial Likelihood

Cox model parameters are estimated using partial likelihood. For subject $i$ who experienced an event at time $t_{(i)}$ , consider the conditional probability that subject $i$ — among all subjects still at risk at that time $\mathcal{R}(t_{(i)})$ — is the one who experiences the event:

L(\beta) = \prod_{i:\text{event}} \frac{\exp(X_i'\beta)}{\sum_{j \in \mathcal{R}(t_{(i)})} \exp(X_j'\beta)}

Each factor corresponds to the conditional probability $h(t_{(i)}|X_i) / \sum_j h(t_{(i)}|X_j)$ within the risk set at time $t_{(i)}$ . Substituting $h(t|X) = h_0(t)\exp(X'\beta)$ , the $h_0(t_{(i)})$ terms cancel between numerator and denominator, so estimating $\beta$ does not require knowing $h_0(t)$ . Although the partial likelihood is not a full likelihood, it has been shown to yield estimators with the same asymptotic properties as maximum likelihood — consistency and asymptotic normality (Cox, 1975).

Interpreting Hazard Ratios

$\exp(\beta_j)$ is interpreted as the hazard ratio (HR):

HR > 1: A one-unit increase in $X_j$ increases the hazard by $(\text{HR} - 1) \times 100\%$
HR < 1: The hazard decreases by $(1 - \text{HR}) \times 100\%$
HR = 1: $X_j$ has no effect on the hazard

When the confidence interval for the hazard ratio does not include 1, the data are inconsistent with the null hypothesis of no effect. The width of the confidence interval reflects estimation precision: a narrow interval indicates a more precise estimate, while a wide interval indicates limited information from the data. Hazard ratios directly convey the direction and magnitude of the effect, making them more informative than the p-value alone.

Test Statistics

Three test statistics are reported for the null hypothesis that all $\beta = 0$ . Each is asymptotically $\chi^2$ with $p$ (number of covariates) degrees of freedom, and they converge to the same value in large samples.

Likelihood Ratio Test

\Lambda = 2\bigl[\ell(\hat\beta) - \ell(0)\bigr]

Based on the difference in log partial likelihood between the null model ( $\beta = 0$ ) and the fitted model. It reflects the global shape of the likelihood surface rather than relying on the local approximations used by the Wald and Score tests.

Wald Test

W = \hat\beta' \, \hat{I}(\hat\beta) \, \hat\beta

Based on the MLE $\hat\beta$ and the estimated information matrix $\hat{I}(\hat\beta)$ . This is a local approximation using the curvature of the log-likelihood at the MLE. The per-covariate p-values in the coefficients table are the univariate version of this test: $z_j = \hat\beta_j / SE(\hat\beta_j)$ . Comparing the overall Wald statistic with individual $z$ values can reveal inconsistencies — for example, individually significant covariates that are not jointly significant — which may indicate multicollinearity or estimation instability.

The weakness of the Wald test is that it fails to capture asymmetry of the likelihood surface when $\hat\beta$ is large.

Score Test

S = U(0)' \, I(0)^{-1} \, U(0)

$U(0)$ is the score function (gradient of the log partial likelihood) evaluated at $\beta = 0$ , and $I(0)$ is the information matrix at that point. Since it does not require $\hat\beta$ , it can be computed even when Newton-Raphson fails to converge.

The typical non-convergence scenario is when a covariate nearly perfectly predicts events. In a small clinical trial where events concentrate in the treatment group, $\hat\beta \to \infty$ and the iterations diverge. Because the model parameterizes the hazard ratio as $\exp(\beta)$ , an extremely large effect causes the coefficient to diverge. The likelihood ratio and Wald tests cannot be computed in this case, but the Score test can still provide evidence that covariates are associated with survival. However, the effect size cannot be estimated. Firth correction and exact methods can address this, but MIDAS does not currently support them.

References

Cox, D. R. (1972). Regression models and life-tables. Journal of the Royal Statistical Society: Series B, 34(2), 187-220. https://www.jstor.org/stable/2985181
Kaplan, E. L., & Meier, P. (1958). Nonparametric estimation from incomplete observations. Journal of the American Statistical Association, 53(282), 457-481. https://www.jstor.org/stable/2281868
Cox, D. R. (1975). Partial likelihood. Biometrika, 62(2), 269-276. https://www.jstor.org/stable/2335362
Struthers, C. A., & Kalbfleisch, J. D. (1986). Misspecified proportional hazard models. Biometrika, 73(2), 363-369. https://www.jstor.org/stable/2336212
Schoenfeld, D. (1982). Partial residuals for the proportional hazards regression model. Biometrika, 69(1), 239-241. https://www.jstor.org/stable/2335876
Grambsch, P. M., & Therneau, T. M. (1994). Proportional hazards tests and diagnostics based on weighted residuals. Biometrika, 81(3), 515-526. https://www.jstor.org/stable/2337123