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. The third column returned by the helper is Greenwood’s variance sum, used in Exercise 1 to show uncertainty.

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)\). This is the likelihood used by the parametric AFT example in Exercise 3.

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.

Exercises

1. Survival at a horizon, with uncertainty. Compare latent sample survival, Kaplan–Meier, and the two naive estimates. Show how much follow-up remains at each horizon.

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)
T = g.weibull(1.5, 20000) * 100
C = g.uniform(0, 260, 20000)
obs, event = np.minimum(T, C), (T <= C).astype(int)
s = km(obs, event)
for h in (50, 100, 150, 200):
    if h > obs.max():
        raise ValueError("Horizon is beyond observed follow-up")
    rows = s[s[:, 0] <= h]
    shat, v = rows[-1, 1:] if len(rows) else (1.0, 0.0)
    if 0 < shat < 1:
        q = np.log(-np.log(shat))
        se = np.sqrt(v) / abs(np.log(shat))
        lo, hi = np.exp(-np.exp(q + 1.96*se)), np.exp(-np.exp(q - 1.96*se))
    else:
        lo, hi = np.nan, np.nan
    print(f"t={h:3d} latent {(T>h).mean():.4f} KM {shat:.4f} "
          f"95% pointwise [{lo:.4f},{hi:.4f}] at risk {(obs>=h).sum()}")
    print(f"      naive all {(obs>h).mean():.4f} event-only {(obs[event==1]>h).mean():.4f}")
# t= 50 latent 0.7003 KM 0.7015 95% pointwise [0.6946,0.7082] at risk 11318
#       naive all 0.5659 event-only 0.5971
# t=100 latent 0.3673 KM 0.3678 95% pointwise [0.3600,0.3756] at risk 4596
#       naive all 0.2298 event-only 0.2308
# t=150 latent 0.1567 KM 0.1595 95% pointwise [0.1528,0.1664] at risk 1364
#       naive all 0.0682 event-only 0.0605
# t=200 latent 0.0590 KM 0.0626 95% pointwise [0.0572,0.0683] at risk 297
#       naive all 0.0149 event-only 0.0097

2. Unknown horizon labels and probability error. Compare the fraction censored before a horizon with the actual label error, then inspect both ranking and average predicted probability in a simulated audit set.

import numpy as np
from sklearn.linear_model import LogisticRegression
from sklearn.metrics import roc_auc_score

h = 100
for mechanism in ("independent", "shared score", "x0 strongly"):
    g = np.random.default_rng(0)
    X = g.normal(size=(30000, 4))
    b = np.array([0.8, -0.5, 0.3, 0.0])
    T = g.weibull(1.5, 30000) * 100 * np.exp(-X @ b)
    base_C = g.uniform(0, 260, 30000)
    C = base_C if mechanism == "independent" else base_C * np.exp(
        -0.8 * (X @ b) if mechanism == "shared score" else -2.5 * X[:, 0])
    obs, event = np.minimum(T, C), (T <= C).astype(int)
    truth = (T <= h).astype(int)
    unknown = (event == 0) & (obs < h)
    label = ((obs <= h) & (event == 1)).astype(int)
    tr, te = np.arange(20000), np.arange(20000, 30000)
    print(f"{mechanism}: unknown {unknown[te].mean():.3f} "
          f"false zeros {(truth[te] != label[te]).mean():.3f} true event rate {truth[te].mean():.3f}")
    for name, rows, target in (("label zero", tr, label),
                               ("drop unknown", tr[~unknown[tr]], label),
                               ("latent-label fit", tr, truth)):
        model = LogisticRegression(max_iter=2000).fit(X[rows], target[rows])
        prob = model.predict_proba(X[te])[:, 1]
        print(f"  {name:16s} AUC {roc_auc_score(truth[te], prob):.4f} mean prediction {prob.mean():.3f}")
    print(f"  generating score AUC {roc_auc_score(truth[te], X[te] @ b):.4f}")
# independent: unknown 0.249 false zeros 0.098 true event rate 0.602
#   label zero       AUC 0.8885 mean prediction 0.509
#   drop unknown     AUC 0.8885 mean prediction 0.637
#   latent-label fit AUC 0.8884 mean prediction 0.605
#   generating score AUC 0.8885
# shared score: unknown 0.233 false zeros 0.140 true event rate 0.602
#   label zero       AUC 0.8882 mean prediction 0.467
#   drop unknown     AUC 0.8884 mean prediction 0.633
#   latent-label fit AUC 0.8884 mean prediction 0.605
#   generating score AUC 0.8885
# x0 strongly: unknown 0.385 false zeros 0.287 true event rate 0.602
#   label zero       AUC 0.6274 mean prediction 0.316
#   drop unknown     AUC 0.8866 mean prediction 0.642
#   latent-label fit AUC 0.8884 mean prediction 0.605
#   generating score AUC 0.8885

3. Compare a censored AFT fit with naive regression. Use the same Weibull model family as the generator, fixing its shape at 1.5 to isolate the treatment of censoring. Compare time-scale slopes, not Cox log hazard ratios.

import numpy as np
from scipy.optimize import minimize
from sklearn.linear_model import LinearRegression

def fit_aft(X, time, event, shape=1.5):
    D = np.column_stack((np.ones(len(X)), X))
    logt = np.log(time)
    def loss_grad(theta):
        eta = D @ theta
        z = np.exp(shape * (logt - eta))
        loss = np.mean(z - event * (np.log(shape) - logt + shape * (logt - eta)))
        gradient = shape * D.T @ (event - z) / len(time)
        return loss, gradient
    initial = np.r_[np.log(np.median(time)), np.zeros(X.shape[1])]
    fit = minimize(loss_grad, initial, jac=True, method="BFGS", options={"gtol":1e-7})
    if not fit.success:
        raise RuntimeError(fit.message)
    return fit.x

beta = np.array([0.8, -0.5, 0.3, 0.0])
g = np.random.default_rng(0)
X = g.normal(size=(30000, 4))
T = g.weibull(1.5, 30000) * 100 * np.exp(-X @ beta)
u = g.uniform(size=30000)
print("censored method     x0       x1       x2")
for upper in (1e9, 260, 60):
    C = upper * u
    time, event = np.minimum(T, C), (T <= C).astype(int)
    ols = LinearRegression().fit(X, np.log(time)).coef_
    aft = fit_aft(X, time, event)[1:]
    for name, coef in (("OLS", ols), ("AFT", aft)):
        print(f"{1-event.mean():.4f}   {name:3s} " + " ".join(f"{v:8.4f}" for v in coef[:3]))
print("true AFT slopes", -beta[:3])
# censored method     x0       x1       x2
# 0.0000   OLS  -0.7955   0.5020  -0.3004
# 0.0000   AFT  -0.7973   0.5029  -0.3009
# 0.3977   OLS  -0.4791   0.3114  -0.1777
# 0.3977   AFT  -0.8070   0.5054  -0.3020
# 0.7640   OLS  -0.1893   0.1296  -0.0667
# 0.7640   AFT  -0.8086   0.5168  -0.3085
# true AFT slopes [-0.8  0.5 -0.3]

Further reading

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