Evaluating Survival Predictions › Implementation in Python
Machine Learning from ScratchAll chapters Evaluating Survival Predictions › Implementation in PythonSection view · Sections in this chapter7 / 10 ↓
Reading tools Definitions, links & notes
Choose your expertise level to adjust how many terms are explained. Beginners see more tooltips, experts see fewer to maintain reading flow. Hover over underlined terms for instant definitions.
Article links
Make inline references clickable
The six-unit cohort validates the arithmetic, but it cannot tell us how the evaluation behaves as a statistical estimator. We therefore simulate a realistic equipment cohort with a known generating mechanism, fit three models on a training partition, estimate the censoring weights on that same training partition, and score the fitted models on a test partition they never saw. The generating Weibull AFT model also has proportional hazards, so both fitted model families have appropriate population forms. Their finite-sample estimates and held-out metrics can still differ.
Simulating a censored equipment cohort
PythonIn [22]Expand
rng = np.random.default_rng(20240501)
n_units = 900
temp_z = rng.normal(size=n_units)
cycles_z = rng.normal(size=n_units)
preventive = rng.binomial(1, 0.40, size=n_units)
log_scale = np.log(1500.0) - 0.35 * temp_z - 0.45 * cycles_z + 0.30 * preventive
rho_true = 1.6
t_failure = np.exp(log_scale) * (-np.log(rng.uniform(size=n_units))) ** (
1.0 / rho_true
)
t_dropout = rng.exponential(scale=2500.0, size=n_units)
t_admin = 3600.0
t_censor = np.minimum(t_dropout, t_admin)
equipment = pd.DataFrame(
{
"temp_z": temp_z,
"cycles_z": cycles_z,
"preventive": preventive,
"time": np.minimum(t_failure, t_censor),
"event": (t_failure <= t_censor).astype(int),
}
)The failure times follow with and a standard Gumbel random variable for minima, so hotter and busier units fail sooner and preventive maintenance extends life. Equivalently, the shape parameter is , matching in the simulation. Censoring comes from two sources: an exponential dropout process with a 2,500-hour mean, representing units that are retired or lose their sensor feed, and an administrative cut at 3,600 hours when the study closes. Both censoring sources are generated independently of failure times and prediction covariates, which supports the pooled censoring weights used below. The controlled simulation lets us check fitting and evaluation separately under stated assumptions. Agreement here does not establish validity under a different censoring mechanism or population.
ConsoleOut [23]Expand
temp_z cycles_z preventive time event -0.326967 1.040655 0 145.142787 1 -0.974315 0.307145 0 1002.431850 0 0.494588 0.895817 0 432.954549 1 units : 900 observed failures : 490 (54.4%) censored observations : 410 (45.6%) median follow-up : 784 h units censored at the 3600 h cut: 25
Roughly 45 percent of this cohort is censored. That overall fraction alone does not determine the weight stability or the bias of a naive rule; inspect the censoring survival at the scored horizons.
Splitting, standardising, and fitting the three models
The split is part of the evaluation design, not an afterthought. We keep a training partition for fitting, a validation partition for any model-choice step, and a test partition reserved for a final, prespecified evaluation phase.
PythonIn [24]Expand
features = ["temp_z", "cycles_z", "preventive"]
continuous = ["temp_z", "cycles_z"]
order_idx = rng.permutation(len(equipment))
train = equipment.iloc[order_idx[:540]].reset_index(drop=True)
val = equipment.iloc[order_idx[540:720]].reset_index(drop=True)
test = equipment.iloc[order_idx[720:]].reset_index(drop=True)
mu_train = train[continuous].mean()
sd_train = train[continuous].std(ddof=0)
def design(frame):
"""Standardise continuous columns with training statistics; keep the indicator raw."""
scaled = (frame[continuous] - mu_train) / sd_train
return np.column_stack(
[
scaled.to_numpy(dtype=float),
frame["preventive"].to_numpy(dtype=float),
]
)
X_train, X_val, X_test = design(train), design(val), design(test)
t_train = train["time"].to_numpy(dtype=float)
e_train = train["event"].to_numpy(dtype=float)
t_val = val["time"].to_numpy(dtype=float)
e_val = val["event"].to_numpy(dtype=float)
t_test = test["time"].to_numpy(dtype=float)
e_test = test["event"].to_numpy(dtype=float)ConsoleOut [25]Expand
partition units events censoring rate
train 540 290 0.463
validation 180 94 0.478
test 180 106 0.411Standardisation uses training means and standard deviations only, so no test-set information leaks into the design matrix. The preventive-maintenance indicator stays at 0/1 because rescaling a binary column would make its coefficient awkward to interpret.
The Cox partial likelihood is evaluated with Breslow's handling of tied event times, in a single function that returns the log-likelihood, its gradient, and its Hessian simultaneously.
PythonIn [26]Expand
def cox_pl(beta, X, time, event):
"""Breslow partial log-likelihood, gradient, and Hessian for Cox regression."""
order = np.argsort(time, kind="mergesort")
X_s, t_s, d_s = X[order], time[order], event[order]
eta = X_s @ beta
w = np.exp(eta)
risk_sum = np.cumsum(w[::-1])[::-1]
wx = w[:, None] * X_s
risk_x = np.cumsum(wx[::-1], axis=0)[::-1]
risk_xx = np.cumsum((wx[:, :, None] * X_s[:, None, :])[::-1], axis=0)[::-1]
ev = np.flatnonzero(d_s == 1)
starts = np.searchsorted(t_s, t_s[ev], side="left")
loglik = float(np.sum(eta[ev] - np.log(risk_sum[starts])))
mean_x = risk_x[starts] / risk_sum[starts][:, None]
grad = np.sum(X_s[ev] - mean_x, axis=0)
hess = -np.sum(
risk_xx[starts] / risk_sum[starts][:, None, None]
- mean_x[:, :, None] * mean_x[:, None, :],
axis=0,
)
return loglik, grad, hess
def fit_cox(X, time, event, ridge=1e-6, n_iter=40, tol=1e-10):
"""Newton-Raphson on the Breslow partial likelihood, with a tiny ridge term."""
beta = np.zeros(X.shape[1])
trace = []
for _ in range(n_iter):
loglik, grad, hess = cox_pl(beta, X, time, event)
trace.append(loglik)
step = np.linalg.solve(-hess + ridge * np.eye(X.shape[1]), grad)
beta = beta + step
if np.max(np.abs(step)) < tol:
break
return beta, np.array(trace)The gradient sums over event times. Each contribution compares the failing unit's covariates with the risk-weighted average among units still alive at that time. Here:
- : the covariate vector of the unit that failed at
- : the average of the risk-set covariates weighted by , namely
If failing units have higher temperatures than their risk-weighted risk-set averages, their contributions to the temperature score are positive. The joint Newton coefficient update also depends on the Hessian and the other covariates. The Hessian provides the curvature that turns this comparison into a Newton step.
PythonIn [27]Expand
def cox_baseline(X, time, event, beta):
"""Breslow baseline cumulative hazard, returned as a step function."""
order = np.argsort(time, kind="mergesort")
X_s, t_s, d_s = X[order], time[order], event[order]
risk_sum = np.cumsum(np.exp(X_s @ beta)[::-1])[::-1]
uniq, first, _ = np.unique(t_s, return_index=True, return_counts=True)
events_at = np.add.reduceat(d_s, first)
return uniq, np.cumsum(events_at / risk_sum[first])
def cox_survival(tau, X, beta, h0_times, h0_vals):
"""Cox survival probabilities at a single horizon tau for every row of X."""
h0 = step_value(tau, h0_times, h0_vals, default=0.0)
return np.exp(-h0 * np.exp(X @ beta))
def weibull_survival(tau, X, mu, gamma, rho):
"""Weibull AFT survival probabilities at a single horizon tau for every row of X."""
if tau == 0:
return np.ones(X.shape[0], dtype=float)
log_lam = mu + X @ gamma
return np.exp(-np.exp(rho * (np.log(tau) - log_lam)))
def fit_weibull_aft(X, time, event):
"""Maximum likelihood for a Weibull AFT model with independent censoring."""
time = np.asarray(time, dtype=float)
event = np.asarray(event, dtype=float)
log_t = np.log(time)
ev = event == 1
design_ev = np.column_stack([np.ones(ev.sum()), X[ev]])
init, *_ = np.linalg.lstsq(design_ev, log_t[ev], rcond=None)
theta0 = np.concatenate([init[:1], init[1:], [np.log(1.5)]])
def neg_loglik(theta):
mu_, gamma_ = theta[0], theta[1 : 1 + X.shape[1]]
rho_ = np.exp(theta[-1])
log_lam = mu_ + X @ gamma_
z = np.exp(rho_ * (log_t - log_lam))
ll = np.sum(
event * (theta[-1] + (rho_ - 1.0) * (log_t - log_lam) - log_lam - z)
)
ll -= np.sum((1.0 - event) * z)
return -ll
res = optimize.minimize(
neg_loglik,
theta0,
method="L-BFGS-B",
options={"maxiter": 5000, "ftol": 1e-12, "gtol": 1e-9},
)
return res.x[0], res.x[1 : 1 + X.shape[1]], float(np.exp(res.x[-1])), resThe AFT likelihood is built from the two contributions we derived earlier. A unit observed to fail at contributes its log density , and a censored unit contributes only , where and . The expression inside neg_loglik multiplies the first by the event indicator and adds the second contribution for the non-events, so one formula covers both cases.
PythonIn [28]Expand
beta_hat, cox_trace = fit_cox(X_train, t_train, e_train)
h0_times, h0_vals = cox_baseline(X_train, t_train, e_train, beta_hat)
mu_hat, gamma_hat, rho_hat, aft_res = fit_weibull_aft(X_train, t_train, e_train)
km_times, km_surv, _ = km_curve(t_train, e_train)ConsoleOut [29]Expand
Cox log hazard ratio AFT gamma implied by AFT: -rho * gamma temp_z 0.600 -0.366 0.581 cycles_z 0.690 -0.423 0.671 preventive -0.522 0.314 -0.498 AFT scale exp(mu) at the reference profile : 1586 h AFT shape rho : 1.587 Cox Newton iterations : 5 AFT optimizer converged : True
The fitted coefficient patterns are close in this seeded Weibull example. The Cox log hazard ratio for a covariate is close to times the AFT coefficient, since a Weibull hazard of the form can be written as a baseline hazard multiplied by , which is a proportional hazards model with . The fitted shape is close to the generating shape of 1.6. At zero standardized continuous covariates and preventive = 0, the fitted scale is about 1,586 hours, compared with a generating scale of about 1,552 hours at the training-mean covariate profile. The simulation's 1,500-hour baseline refers to zero raw covariates and preventive = 0. This agreement checks the fitted coefficient pattern in this seeded example; the derivative and reference-distribution checks below test additional calculations.
Computing the evaluation metrics
We estimate the censoring distribution once, on the training partition, and reuse it everywhere downstream. Fixing the horizons before touching the test set keeps the evaluation honest.
PythonIn [30]Expand
g_times, g_surv = censoring_curve(t_train, e_train)
horizons = np.array([800.0, 1500.0, 2500.0])
grid = np.linspace(300.0, 3000.0, 28)
tau_star = 3000.0 # preselected earlier-failure cutoff for IPCW concordance
g_at_horizons = step_value(horizons, g_times, g_surv)ConsoleOut [31]Expand
tau = 800 h G_hat(tau) = 0.709 weight 1/G_hat = 1.41 tau = 1500 h G_hat(tau) = 0.547 weight 1/G_hat = 1.83 tau = 2500 h G_hat(tau) = 0.378 weight 1/G_hat = 2.64 G_hat(3000 h) = 0.312 events at or before 2500 h in the test partition: 104 units still at risk at 2500 h in the test partition: 10
Both positivity requirements hold at all three horizons: , , and , each positive. At 2,500 hours, the retained test output reports 104 observed failures at or before the horizon and 10 units observed strictly beyond it. Positive censoring support makes the weights defined; it does not establish the precision of the score. Report these counts with uncertainty rather than treating a positive estimate as an adequacy threshold. At or beyond 3,600 hours, the administrative censoring jump makes , so this strict IPCW Brier implementation rejects the horizon. That support restriction does not make the outcome of a unit censored exactly at the cut unknown; it excludes an inclusive boundary estimator from this implementation.
PythonIn [32]Expand
def _metric_inputs(values, time, event):
time, event = _survival_inputs(time, event)
values = np.asarray(values, dtype=float)
if values.shape != time.shape or not np.all(np.isfinite(values)):
raise ValueError(
"predictions must be finite and aligned with observations"
)
return values, time, event
def _censoring_values(query, g_times, g_surv, side="right"):
g_times = np.asarray(g_times, dtype=float)
g_surv = np.asarray(g_surv, dtype=float)
if g_times.size == 0 or g_surv.shape != g_times.shape:
raise ValueError("a nonempty aligned censoring curve is required")
if (
not np.all(np.isfinite(g_surv))
or np.any((g_surv < 0) | (g_surv > 1))
or np.any(np.diff(g_surv) > 0)
):
raise ValueError(
"censoring survival must be finite, nonincreasing, and in [0,1]"
)
if np.any(np.asarray(query) > g_times[-1]):
raise ValueError(
"evaluation extends beyond observed censoring follow-up"
)
result = step_value(query, g_times, g_surv, side=side)
if np.any(np.asarray(result) <= 0):
raise ValueError(
"censoring survival must be positive at contributing times"
)
return result
def harrell_c(risk, time, event):
"""Comparable-pair concordance, giving half credit to tied risk scores."""
risk, time, event = _metric_inputs(risk, time, event)
usable = (event[:, None] == 1) & (time[:, None] < time[None, :])
if not usable.any():
raise ValueError("concordance requires at least one comparable pair")
higher = usable & (risk[:, None] > risk[None, :])
tied = usable & (risk[:, None] == risk[None, :])
return float((higher.sum() + 0.5 * tied.sum()) / usable.sum())
def uno_c(risk, time, event, g_times, g_surv, tau_star):
"""Truncated IPCW concordance; continuous censoring gives Uno's G^-2 weight."""
risk, time, event = _metric_inputs(risk, time, event)
if not np.isfinite(tau_star) or tau_star <= 0:
raise ValueError("tau_star must be a fixed positive finite horizon")
_censoring_values(tau_star, g_times, g_surv)
usable = (
(time[:, None] < time[None, :])
& (event[:, None] == 1)
& (time[:, None] < tau_star)
)
if not usable.any():
raise ValueError("truncated concordance requires a comparable pair")
earlier = usable.any(axis=1)
g_before = _censoring_values(time[earlier], g_times, g_surv, side="left")
g_after = _censoring_values(time[earlier], g_times, g_surv)
weights = np.zeros_like(time)
log_weights = -np.log(g_before) - np.log(g_after)
weights[earlier] = np.exp(log_weights - log_weights.max())
pair_weights = np.where(usable, weights[:, None], 0.0)
credit = (risk[:, None] > risk[None, :]) + 0.5 * (
risk[:, None] == risk[None, :]
)
return float((pair_weights * credit).sum() / pair_weights.sum())
def ipcw_brier(surv_pred, time, event, tau, g_times, g_surv):
"""IPCW Brier score: failure weights use G(t-), survivor weights G(tau)."""
surv_pred, time, event = _metric_inputs(surv_pred, time, event)
if (
np.any((surv_pred < 0) | (surv_pred > 1))
or not np.isfinite(tau)
or tau < 0
):
raise ValueError(
"survival predictions must be probabilities and tau nonnegative"
)
g_tau = float(_censoring_values(tau, g_times, g_surv))
failed = (time <= tau) & (event == 1)
alive = time > tau
g_before = _censoring_values(time[failed], g_times, g_surv, side="left")
terms = np.zeros_like(time)
terms[failed] = surv_pred[failed] ** 2 / g_before
terms[alive] = (1.0 - surv_pred[alive]) ** 2 / g_tau
return float(terms.mean())
def naive_brier(surv_pred, time, event, tau):
"""Brier score that treats every censored-before-tau unit as a survivor."""
time = np.asarray(time, dtype=float)
event = np.asarray(event, dtype=float)
outcome = 1.0 - ((time <= tau) & (event == 1)).astype(float)
return float(np.mean((outcome - np.asarray(surv_pred, dtype=float)) ** 2))
def trapezoid(y, x):
"""Trapezoid rule, written out so the integral is visible."""
return float(np.sum(0.5 * (y[1:] + y[:-1]) * np.diff(x)))
def predict_survival(tau, model):
"""Survival probabilities at one horizon for every test unit."""
if model == "Kaplan-Meier":
return np.full(len(test), float(step_value(tau, km_times, km_surv)))
if model == "Cox PH":
return cox_survival(tau, X_test, beta_hat, h0_times, h0_vals)
return weibull_survival(tau, X_test, mu_hat, gamma_hat, rho_hat)For Uno's ratio, the code evaluates relative row weights. It computes for failures with at least one eligible partner, then uses over those same contributing rows. Multiplying all contributing weights by one common positive factor leaves the numerator-to-denominator ratio unchanged. A row with no eligible partner does not set that scale. This avoids forming an underflowing probability product or an overflowing reciprocal. Extremely small relative weights can still round to zero in floating-point arithmetic.
The risk scores follow directly. The Cox model contributes its linear predictor, the AFT model contributes the negated log-scale term so that larger means higher risk, and the Kaplan–Meier fit contributes a constant, which is the honest summary of a model with no covariates.
PythonIn [33]Expand
model_names = ["Kaplan-Meier", "Cox PH", "Weibull AFT"]
risk_scores = {
"Kaplan-Meier": np.zeros(len(test)),
"Cox PH": X_test @ beta_hat,
"Weibull AFT": -(X_test @ gamma_hat),
}
surv_at = {}
for name in model_names:
for tau in horizons:
surv_at[(name, tau)] = predict_survival(tau, name)
val_surv_1500 = {
"Cox PH": cox_survival(1500.0, X_val, beta_hat, h0_times, h0_vals),
"Weibull AFT": weibull_survival(1500.0, X_val, mu_hat, gamma_hat, rho_hat),
}
val_brier = {
name: ipcw_brier(pred, t_val, e_val, 1500.0, g_times, g_surv)
for name, pred in val_surv_1500.items()
}The validation comparison is a small selection step: both covariate models are scored on the validation partition at the 1,500-hour horizon, and whichever is better would be the one we carry forward. The test partition is not consulted for this choice. The table below shows the result.
ConsoleOut [34]Expand
Validation IPCW Brier at tau = 1500 h Cox PH : 0.1611 Weibull AFT : 0.1616 Kaplan-Meier marginal baseline: 0.2498
PythonIn [35]Expand
rows = []
for name in model_names:
row = {
"model": name,
"Harrell C": harrell_c(risk_scores[name], t_test, e_test),
"Uno C (<3000 h)": uno_c(
risk_scores[name], t_test, e_test, g_times, g_surv, tau_star
),
}
for tau in horizons:
row[f"Brier {int(tau)}h"] = ipcw_brier(
surv_at[(name, tau)], t_test, e_test, tau, g_times, g_surv
)
rows.append(row)
metrics = pd.DataFrame(rows).set_index("model").round(4)ConsoleOut [36]Expand
Harrell C Uno C (<3000 h) Brier 800h Brier 1500h Brier 2500h model Kaplan-Meier 0.5000 0.5000 0.2614 0.2377 0.1366 Cox PH 0.7095 0.7110 0.1972 0.1968 0.0935 Weibull AFT 0.7100 0.7116 0.1991 0.1928 0.0889
The bars below compare the three chosen decision horizons within each model, using the rounded values in the table. The later horizon curves reuse those same horizons within a denser grid and add the naive censoring rule, so they show the path and the effect of outcome imputation rather than a separate evaluation sample.
VisualizationOut [38]Expand
PythonIn [39]Expand
brier_curves = {}
for name in model_names:
brier_curves[name] = np.array(
[
ipcw_brier(
predict_survival(t, name), t_test, e_test, t, g_times, g_surv
)
for t in grid
]
)
brier_curve_naive = np.array(
[
naive_brier(predict_survival(t, "Weibull AFT"), t_test, e_test, t)
for t in grid
]
)
ibs = {
name: trapezoid(brier_curves[name], grid) / (grid[-1] - grid[0])
for name in model_names
}ConsoleOut [40]Expand
Integrated Brier score over 300-3000 h Kaplan-Meier : 0.1986 Cox PH : 0.1513 Weibull AFT : 0.1482 Naive versus IPCW for the Weibull AFT model tau = 800 h IPCW 0.1991 naive 0.1767 tau = 1500 h IPCW 0.1928 naive 0.2144 tau = 2500 h IPCW 0.0889 naive 0.2330
In this test cohort, Kaplan–Meier lands at exactly 0.5 on both concordance measures because a constant score produces nothing but ties, and its IPCW Brier scores are 0.2614, 0.2377, and 0.1366 at 800, 1,500, and 2,500 hours. For comparison, Weibull AFT scores 0.1928 at 1,500 hours. A marginal curve can be roughly right on average while being useless for ranking. Both covariate models improve on it for discrimination. Report ranking and probability accuracy separately when the decision needs both. The naive Brier rule imputes survival to unknown horizon outcomes. Its difference from IPCW in this finite sample reflects both the effect of that imputation and sampling variation; it is not a direct measurement of the population bias.
VisualizationOut [42]Expand
The censoring curve supplies the weights for IPCW Brier error and truncated IPCW concordance. It declines gradually as units drop out, which inflates the weights modestly, and then it hits zero at the 3,600-hour administrative cut, where every remaining unit is censored at once. The jump records the end of observation. Units censored exactly at the cut are known to have survived past it, but the strict estimator here excludes them and has zero right-side censoring support at the cut. Follow-up alone does not reveal their subsequent failure times, so evaluating later outcomes would require additional assumptions or observation.
VisualizationOut [44]Expand
Three features of this figure carry the argument. First, the naive and IPCW curves differ because they use different rules for unknown horizon outcomes. Their difference in this finite sample combines censoring bias with sampling variation, and its direction can change with the horizon. Second, the two covariate models separate visibly from Kaplan–Meier in the middle of this plotted range. That visual separation does not establish where a comparison has greatest statistical power. Third, each covariate model has its largest plotted IPCW Brier value at 1,100 h. The IPCW curves have lower values toward the end of this window overall. Brier levels at different horizons reflect changing binary-outcome uncertainty as well as forecast performance, so the peak does not isolate the worst calibration. Near the end of follow-up, shrinking survivor counts and small censoring survival make the estimates sensitive to a few weighted contributions; report that support when interpreting the comparison.
PythonIn [45]Expand
def calibration_bins(surv_pred, time, event, tau, n_bins=5):
"""Group units by predicted survival and compare with Kaplan-Meier within each group."""
surv_pred = np.asarray(surv_pred, dtype=float)
order = np.argsort(surv_pred)
predicted, observed, lower, upper = [], [], [], []
for idx in np.array_split(order, n_bins):
u, s, v = km_curve(time[idx], event[idx])
obs = float(step_value(tau, u, s))
greenwood = float(step_value(tau, u, v, default=np.nan))
has_failure = np.any((event[idx] == 1) & (time[idx] <= tau))
se = (
float(np.sqrt(greenwood))
if has_failure and np.isfinite(greenwood) and greenwood > 0
else np.nan
)
predicted.append(surv_pred[idx].mean())
observed.append(obs)
lower.append(obs - 1.96 * se)
upper.append(obs + 1.96 * se)
return tuple(np.array(arr) for arr in (predicted, observed, lower, upper))
tau_cal = 1500.0
cal_cox = calibration_bins(
surv_at[("Cox PH", tau_cal)], t_test, e_test, tau_cal
)
cal_aft = calibration_bins(
surv_at[("Weibull AFT", tau_cal)], t_test, e_test, tau_cal
)
ici_cox = float(np.mean(np.abs(cal_cox[1] - cal_cox[0])))
ici_aft = float(np.mean(np.abs(cal_aft[1] - cal_aft[0])))ConsoleOut [46]Expand
Mean absolute calibration error at tau = 1500 h Cox PH : 0.1645 Weibull AFT : 0.1556
VisualizationOut [48]Expand
The mean absolute deviations are about 0.165 for Cox and 0.156 for Weibull AFT. Several quintiles sit appreciably below the diagonal; for example, the fourth Cox quintile predicts about 0.662 survival while its Kaplan–Meier estimate is about 0.302. The data were generated by a Weibull model and the Cox specification is correctly proportional. The error bars matter as much as the points. Each quintile contains 36 units, with 28, 23, 18, 15, and 9 observed failures by 1,500 hours in the Cox quintiles, so the observed value carries real uncertainty. For quintiles whose approximate pointwise 95 percent intervals cross the diagonal, the corresponding pointwise Wald comparison does not reject equality at the 5 percent level under the normal approximation, treating the fitted forecasts and binning as fixed. This does not establish calibration, demonstrate equivalence, or rule out smaller departures. These intervals are not a simultaneous test across quintiles. Show the intervals so that readers can assess the precision of the displayed estimates.
For this display, a bin's normal interval is unavailable when no failure is observed by the horizon or its Greenwood variance is nonpositive or nonfinite; the Kaplan–Meier point estimate is still shown. Passing these checks does not establish approximate coverage when failures are few.
Key Parameters
The list below covers evaluation settings and numerical controls used in this example.
horizons(): evaluation cut points. Fix them before seeing test results and tie them to a decision (inspection interval, warranty term, budget cycle). This strict IPCW Brier estimator requires . Report the number of observed failures by each horizon and the number of units observed strictly beyond it, and assess score uncertainty; positivity alone does not establish precision.- Censoring-weight estimation sample: this example fits the censoring curve on training data and keeps it fixed for evaluation. Evaluation-cohort estimation can also be valid under the censoring assumptions; include its estimation variation in uncertainty and check that the chosen censoring model describes the scored population.
n_bins: number of calibration bins (quintiles give 36 units per bin here). Inspect each bin's sample size, failures before the horizon, follow-up beyond it, censoring support, event-time risk sets, and Greenwood-based uncertainty interval before interpreting the point. Choose an acceptable interval width from the decision context. Fewer bins increase average group size but can hide within-bin miscalibration; an event-count target alone does not establish precision. The last Cox quintile has nine observed failures here.ridge: a diagonal damping term in each Newton solve (default: ). It stabilises the numerical step but does not penalise the objective or shrink the converged coefficients. An L2 penalty would also subtract from the score.- Convergence tolerance (
tol): maximum absolute Newton step below which fitting stops. The default value is ; this is a numerical stopping criterion and does not establish the statistical precision of the coefficient estimates.
Key Methods
km_curvereturns Kaplan–Meier survival and Greenwood variance at distinct observed times. Its undefined risk-exhausting variance is missing.censoring_curveestimates reverse Kaplan–Meier with original failures removed before a tied censoring jump.step_valueevaluates an aligned, sorted step table at query times;side="right"gives its right value andside="left"the value immediately before a jump.cox_plandfit_coxcompute the Breslow partial likelihood, gradient, and Hessian, then solve for the coefficients by Newton–Raphson.cox_baselineandcox_survivalturn the fitted coefficients into predicted survival probabilities through .fit_weibull_aftandweibull_survivaldo the same for the AFT parameterisation, optimising the censored log-likelihood with L-BFGS-B and converting into survival probabilities.harrell_ccomputes comparable-pair concordance from risk scores;uno_calso requires a censoring curve and an explicit supported truncation horizon.ipcw_brier,naive_brier, andtrapezoidproduce the IPCW Brier score, the naive survivor-imputed score, and the integrated summary.calibration_binsgroups units by predicted survival and compares each group's predictions with its Kaplan–Meier estimate.
Verifying the implementation
Three complementary checks compare selected implementation calculations with finite differences, explicit loops, and a reference distribution. The first compares the analytic Cox gradient with a finite-difference approximation. The second reproduces the six-unit hand arithmetic. The third compares our Weibull survival function with SciPy's reference implementation of the same distribution.
PythonIn [49]Expand
beta_check = np.array([0.10, -0.20, 0.05])
_, grad_analytic, _ = cox_pl(beta_check, X_train, t_train, e_train)
eps = 1e-6
grad_numeric = np.empty_like(beta_check)
for k in range(len(beta_check)):
bump = np.zeros_like(beta_check)
bump[k] = eps
plus = cox_pl(beta_check + bump, X_train, t_train, e_train)[0]
minus = cox_pl(beta_check - bump, X_train, t_train, e_train)[0]
grad_numeric[k] = (plus - minus) / (2 * eps)
grad_error = float(np.max(np.abs(grad_analytic - grad_numeric)))
tiny_c_check = harrell_c(tiny_risk, tiny_times, tiny_event)
tiny_brier_check = ipcw_brier(
tiny_pred, tiny_times, tiny_event, tau_tiny, tiny_g_times, tiny_g
)
scale_ref = np.exp(mu_hat + X_test @ gamma_hat)
scipy_surv = stats.weibull_min.sf(1500.0, c=rho_hat, scale=scale_ref)
ours_surv = weibull_survival(1500.0, X_test, mu_hat, gamma_hat, rho_hat)
surv_error = float(np.max(np.abs(scipy_surv - ours_surv)))ConsoleOut [50]Expand
max |analytic - numeric| Cox gradient : 3.27e-07 six-unit Harrell C : 0.8182 six-unit IPCW Brier at tau = 8 : 0.2534 max |Weibull survival - scipy| : 1.11e-16
At the tested coefficient vector, the printed maximum absolute gradient difference checks this derivative with one step size. Step size, function scale, and rounding can affect its error. The six-unit numbers reproduce and to the printed precision, so the vectorised functions agree with the loop-based definitions we wrote by hand. The Weibull survival probabilities agree with SciPy to machine precision at the evaluated inputs, supporting the stated and conversion. This finite comparison does not verify every possible input.
The Cox and transformed AFT coefficient estimates are close, as expected under this Weibull generating model. They need not be identical in a finite sample: Cox estimates its coefficients by partial likelihood without specifying the baseline, then reconstructs a Breslow baseline cumulative hazard from the risk sets and fitted coefficients. The AFT fit instead assumes a Weibull baseline and uses the full censored likelihood.
Uncertainty on the test metrics
A point estimate from a 180-unit test set needs uncertainty. Here we resample test units with replacement and take percentile intervals while holding the fitted models and training censoring curve fixed. These intervals describe test-cohort sampling variation conditional on those fitted objects; they do not include variation from refitting the models or censoring curve.
PythonIn [51]Expand
boot_rng = np.random.default_rng(1234)
n_boot = 200
pred_1500 = {name: surv_at[(name, 1500.0)] for name in model_names}
boot_records = []
for _ in range(n_boot):
pick = boot_rng.integers(0, len(test), len(test))
tb, eb = t_test[pick], e_test[pick]
for name in model_names:
boot_records.append(
{
"model": name,
"c_index": harrell_c(risk_scores[name][pick], tb, eb),
"brier": ipcw_brier(
pred_1500[name][pick], tb, eb, 1500.0, g_times, g_surv
),
}
)
boot = pd.DataFrame(boot_records)ConsoleOut [52]Expand
c_lower c_upper brier_lower brier_upper model Cox PH 0.6596 0.7599 0.1573 0.2397 Kaplan-Meier 0.5000 0.5000 0.2139 0.2659 Weibull AFT 0.6601 0.7603 0.1546 0.2347
The intervals are wide, as they should be. With 180 test units and roughly 41 percent of them censored, the two covariate-model concordance intervals each span about one tenth, and the two covariate models overlap substantially. Reporting a single number and declaring a winner would be overreach; reporting the interval makes the ambiguity part of the result.
Comments
Add to the conversation
A question, a correction, or an example from your own work.
No comments yet. Be the first to share your thoughts!