Survival Analysis

A subscriber who is still active after 40 days has supplied useful information: their cancellation time exceeds 40 days. Recording 40 as the cancellation time discards that distinction; removing the subscriber discards the information entirely. Survival analysis uses observed event times together with these incomplete observations to estimate when an event occurs.

Record both follow-up time and event status

Choose a time origin and an event before building the table: for example, days from subscription to first cancellation. Let \(T\) be the event time and \(C\) the time observation stops. Right-censored data contain \(Y=\min(T,C)\) and \(\delta=1\{T\leq C\}\). If \(\delta=0\), we know \(T>C\), not when the event eventually happens. “Survival” means remaining event-free and applies to equipment failures or customer cancellation as well as lifetimes.

The survival function \(S(t)=P(T>t)\) answers a horizon question: what fraction remain event-free after \(t\)? For continuous event times, the hazard \(h(t)=f(t)/S(t)\) is an instantaneous event rate among those still event-free. Over a short interval of length \(\Delta t\), their event probability is approximately \(h(t)\Delta t\). A hazard is measured per unit time and is not itself a probability. Its cumulative integral \(H(t)\) gives \(S(t)=\exp[-H(t)]\).

The examples assume independent subjects, entry at time zero, and right censoring. Delayed entry requires different risk sets; interval censoring records a range containing the event time. A competing event that prevents the event of interest requires distinguishing cause-specific hazards from the probability of that event occurring. Treating competing events as ordinary censoring in 1−Kaplan–Meier does not generally estimate that probability.

Build a survival curve from risk sets

At an observed event time \(t_j\), let \(n_j\) be the number still under observation and event-free just before that time, and \(d_j\) the number of events. Kaplan–Meier multiplies the estimated conditional survival fractions: \[\widehat S(t)=\prod_{t_j\leq t}\left(1-\frac{d_j}{n_j}\right).\] Censoring reduces future risk sets without creating an event. If events and censoring share a recorded time, the code includes those censored at that time in the event risk set.

For four subjects, suppose events occur at days 2 and 5, while the other two subjects are censored at days 4 and 6. At day 2, survival becomes 3/4. At day 5, only two subjects remain at risk, so survival becomes (3/4)×(1/2)=3/8. The day-4 censoring changes the denominator of the second step; it does not create a drop of its own.

import numpy as np

def km(time, event):
    order = np.argsort(time)
    t, e = np.asarray(time)[order], np.asarray(event)[order]
    times, first, counts = np.unique(t, return_index=True, return_counts=True)
    deaths = np.add.reduceat(e, first)
    risk = len(t) - np.r_[0, np.cumsum(counts[:-1])]
    keep = deaths > 0
    times, deaths, risk = times[keep], deaths[keep], risk[keep]
    survival = np.cumprod(1 - deaths / risk)
    increments = np.full(len(times), np.inf)
    np.divide(deaths, risk * (risk - deaths), out=increments, where=risk > deaths)
    greenwood = np.cumsum(increments)
    return np.column_stack((times, survival, greenwood))

g = np.random.default_rng(0)
n = 20000
T = g.weibull(1.5, n) * 100
C = g.uniform(0, 260, n)
obs, event = np.minimum(T, C), (T <= C).astype(int)
s = km(obs, event)
reached = s[s[:, 1] <= 0.5]
median = reached[0, 0] if len(reached) else np.nan
print(f"latent sample mean {T.mean():.2f} median {np.median(T):.2f} censored {1-event.mean():.4f}")
print(f"observed-time mean {obs.mean():.2f} median {np.median(obs):.2f}")
print(f"event-only mean {obs[event == 1].mean():.2f} median {np.median(obs[event == 1]):.2f}")
print(f"Kaplan-Meier median {median:.2f}")
print("four-subject example", km(np.array([2,4,5,6]), np.array([1,0,1,0]))[:, :2])
# latent sample mean 89.88 median 77.89 censored 0.3425
# observed-time mean 67.64 median 57.77
# event-only mean 69.13 median 60.69
# Kaplan-Meier median 78.07
# four-subject example [[2.    0.75 ]
#  [5.    0.375]]

The observed-time mean is 67.64 and the event-only mean is 69.13, compared with the latent sample mean 89.88. Their medians are 57.77 and 60.69, while Kaplan–Meier gives 78.07 against the latent sample median 77.89. The latent values are available only because this is a simulation; they are sample summaries, not exact population quantities.

Because \(Y\leq T\) for every row, treating follow-up as the event time cannot increase the empirical mean. Under independent censoring, longer event times are less likely to be observed, so retaining only observed events also favors shorter times. Censored subjects are not necessarily the subjects with the longest actual lifetimes: someone censored early may have an event soon afterward.

For the pooled Kaplan–Meier curve here, censoring is independent of event time. If independence holds only conditional on features \(X\), an unadjusted pooled curve need not be valid; stratification or an appropriate conditional model may be needed. Dropout related to an impending event through unobserved information can violate the assumption. Tail estimates also need adequate follow-up: beyond the observed support, drawing a flat line does not establish continued survival.

The estimated median is the first event time where the curve reaches 0.5. If it never does, the median is not reached; the code returns NaN instead of inventing a value. Likewise, the unrestricted mean requires the full tail. A restricted mean survival time integrates the curve only to a stated, supported horizon.

What naive regression estimates

import numpy as np
from sklearn.linear_model import LinearRegression

beta = np.array([0.8, -0.5, 0.3, 0.0])
print(f"{'censoring rate':>15} {'x0':>9} {'x1':>9} {'x2':>9} {'attenuation':>12}")
for upper in (1e9, 400, 260, 160, 100, 60):
    g = np.random.default_rng(0); n = 30000
    X = g.normal(size=(n, 4))
    T = g.weibull(1.5, n) * (100 * np.exp(-X @ beta))     # accelerated failure time
    C = g.uniform(0, upper, n)
    obs, event = np.minimum(T, C), (T <= C).astype(int)
    m = LinearRegression().fit(X, np.log(obs))            # regress on what you observe
    print(f"{1 - event.mean():15.4f} " + " ".join(f"{c:9.4f}" for c in m.coef_[:3])
          + f" {np.mean(m.coef_[:3] / (-beta[:3])):12.4f}")
print(f"{'true -beta':>15} " + " ".join(f"{-b:9.4f}" for b in beta[:3]) + f" {1.0:12.4f}")
#  censoring rate        x0        x1        x2  attenuation
#          0.0000   -0.7955    0.5020   -0.3004       1.0000
#          0.2956   -0.5601    0.3601   -0.2089       0.7055
#          0.3977   -0.4791    0.3114   -0.1777       0.6046
#          0.5253   -0.3778    0.2490   -0.1391       0.4780
#          0.6476   -0.2814    0.1885   -0.1023       0.3566
#          0.7640   -0.1893    0.1296   -0.0667       0.2394
#      true -beta   -0.8000    0.5000   -0.3000       1.0000

The generator is an accelerated failure time (AFT) model: \(T=100\exp(-X\beta)W\), where \(W\) is Weibull with shape 1.5 and unit scale. Taking logs gives \(\log T=\log100-X\beta+\log W\), so the population slopes of log event time are \(-\beta\). Regressing \(\log\min(T,C)\) instead compresses long durations. The displayed “attenuation” is the average estimated-to-true ratio for the three nonzero slopes.

At 29.6% censoring that ratio is about 0.71; at 76.4% it is about 0.24. These approximately shared factors occur in this Gaussian-feature, independent-censoring setup. They are not a correction rule for arbitrary covariate distributions or censoring mechanisms. The zero coefficient is excluded from the ratio calculation. Similar coefficient ratios also do not guarantee identical rankings of every pair in a finite sample.

Fit the likelihood for the observation you have

With conditionally independent censoring and a censoring mechanism that can be ignored for inference about the event model, a row contributes \[f(Y\mid X)^{\delta}S(Y\mid X)^{1-\delta}.\] An observed event contributes a density at its time; a censored row contributes the probability of surviving beyond its follow-up. Thus a subject censored at day 40 contributes \(S(40\mid X)\), not the event density \(f(40\mid X)\).

A Cox model instead writes \(h(t\mid X)=h_0(t)\exp(X\gamma)\). The baseline hazard \(h_0(t)\) is left unspecified when estimating \(\gamma\). With no tied event times, each observed event contributes \(\exp(X_i\gamma)/\sum_{j\in R(t_i)}\exp(X_j\gamma)\) to the partial likelihood. The numerator is the relative hazard of the subject who had the event; the denominator sums relative hazards over everyone at risk. Censored follow-up times are needed to form those risk sets.

A one-unit increase in feature \(x_j\) multiplies the hazard by \(e^{\gamma_j}\), holding other features fixed. In an AFT model, a slope \(\theta_j\) multiplies time quantiles by \(e^{\theta_j}\). These are different quantities. In this particular Weibull generator, \(S(t\mid X)=\exp[-(t/100)^{1.5}\exp(1.5X\beta)]\), so it also satisfies proportional hazards with \(\gamma=1.5\beta\). The next block fits a Cox model and compares log hazard ratios with their generating values.

import numpy as np
from statsmodels.duration.hazard_regression import PHReg

g = np.random.default_rng(42)
X_cox = g.normal(size=(3000, 4))
beta_cox = np.array([0.8, -0.5, 0.3, 0.0])
T_cox = g.weibull(1.5, 3000) * 100 * np.exp(-X_cox @ beta_cox)
C_cox = g.uniform(0, 260, 3000)
Y_cox = np.minimum(T_cox, C_cox)
d_cox = (T_cox <= C_cox).astype(int)
fit_cox = PHReg(Y_cox, X_cox, status=d_cox).fit()
print("feature true log-HR fitted log-HR SE")
for j in range(4):
    print(f"x{j} {1.5*beta_cox[j]:10.3f} {fit_cox.params[j]:13.3f} {fit_cox.bse[j]:.3f}")
# feature true log-HR fitted log-HR SE
# x0      1.200         1.177 0.030
# x1     -0.750        -0.776 0.027
# x2      0.450         0.473 0.025
# x3      0.000        -0.011 0.023

Cox coefficients alone do not give an absolute survival probability. Estimating the baseline cumulative hazard gives \(\widehat S(t\mid X)=\exp[-\widehat H_0(t)\exp(X\widehat\gamma)]\). Proportional hazards assumes that hazard ratios remain constant over time. Schoenfeld residuals compare the features of the subject who had an event with the model-weighted features of the risk set. Time trends in those residuals are one diagnostic; time-varying effects or another model may be needed when that assumption fails. Censoring-aware fitting still depends on model adequacy and follow-up support, and an adjusted association is not automatically causal.

Compare only event orderings that are known

import numpy as np

def build(seed=0, n=30000):
    g = np.random.default_rng(seed)
    X = g.normal(size=(n, 4))
    b = np.array([0.8, -0.5, 0.3, 0.0])
    T = g.weibull(1.5, n) * (100 * np.exp(-X @ b))
    C = g.uniform(0, 260, n)
    return X, T, np.minimum(T, C), (T <= C).astype(int), b

def cindex(time, event, score, cap=400000, seed=0):
    gg = np.random.default_rng(seed)
    i = gg.integers(0, len(time), cap)
    j = gg.integers(0, len(time), cap)
    distinct = i != j
    i, j = i[distinct], j[distinct]
    earlier_i = time[i] < time[j]
    earlier_j = time[j] < time[i]
    comparable = (earlier_i & (event[i] == 1)) | (earlier_j & (event[j] == 1))
    concordant = np.where(earlier_i, score[i] > score[j], score[j] > score[i])
    credit = concordant.astype(float) + 0.5 * (score[i] == score[j])
    value = credit[comparable].mean() if comparable.any() else np.nan
    return value, comparable.mean()

X, T, obs, event, b = build()
te = slice(20000, None)
risk = X[te] @ b
c, frac = cindex(obs[te], event[te], risk)
c_naive, frac_n = cindex(obs[te], np.ones(len(obs[te]), int), risk)
c_true, _ = cindex(T[te], np.ones(len(T[te]), int), risk)
print(f"comparable-pair C-index {c:.4f} usable pair fraction {frac:.4f}")
print(f"ignore-censoring index  {c_naive:.4f} usable pair fraction {frac_n:.4f}")
print(f"latent-time index      {c_true:.4f}")
# comparable-pair C-index 0.8017 usable pair fraction 0.6664
# ignore-censoring index  0.7021 usable pair fraction 1.0000
# latent-time index      0.7881

A concordance index measures ordering of event times. If one subject has an observed event at day 10 and another is observed through day 20, their order is known. If the first is censored at day 10, it is not known: their event might occur at day 12 or day 100. A higher risk score should predict the earlier event. Tied scores receive half credit; this illustrative implementation excludes tied observed times and samples pairs with replacement rather than enumerating them all. It returns NaN when no sampled pair is comparable.

The score uses the generating coefficients, so this is an evaluation demonstration, not performance of a trained model. Ignoring censoring yields 0.7021 against 0.7881 using latent event times. Censoring can make an apparent time ordering differ from the unobserved event ordering, causing either false discordance or false concordance.

The comparable-pair value, 0.8017, is higher here. Such indices can depend on the censoring distribution because it changes which event-time pairs are retained; those pairs need not be systematically easier in every setting. This single simulation also contains sampling variation. The usable fraction refers to pairs, not a fraction of subjects, and does not by itself quantify the bias.

For a target concordance over a stated time range, inverse-probability-of-censoring weighting (IPCW) reweights observable comparisons using an estimated censoring survival distribution. Valid weights require suitable independence assumptions and enough chance of follow-up across that range; feature-dependent censoring can require conditional weights. The horizon and weighting method are part of the reported metric. A time-dependent AUC instead compares cases and controls defined at a specific horizon; it is not identical to concordance.

Good ordering does not establish accurate survival probabilities. Evaluate those on held-out subjects using censoring-adjusted calibration checks and a time-dependent Brier score, which measures squared probability error with censoring weights. Keep preprocessing and model selection within training data, and use temporal or group splits when deployment requires them. No evaluation method recovers unsupported long-horizon outcomes without additional assumptions.

Further reading

• statsmodels survival and duration guide documents Kaplan–Meier estimates, survival uncertainty, and Cox fitting.
• scikit-survival evaluation guide compares concordance, IPCW, time-dependent AUC, and Brier scores.
• R survival project provides further material on competing risks and more general event histories.