Survival Analysis — Time-to-Event in Security#
Security events have duration. Dwell time, time-to-patch, time-to-detect — these are all time-to-event problems, and they share a statistical structure that ordinary summary statistics handle poorly. Censoring is natural: some vulnerabilities are not patched by the end of an observation window, and dropping them biases every estimate. Means mislead here mainly because of censoring, not skew — the median from a Kaplan-Meier curve is almost always more honest.
For a full treatment of time-to-patch survival analysis applied to real vulnerability data, see the Elastic Security Labs article Time-to-Patch Metrics: A Survival Analysis Approach Using Qualys and Elastic (Laura Voicu and Clement Fouque, 2025).
1 — Setup#
import numpy as np
import matplotlib.pyplot as plt
from decision_security.synth import make_rng, survival_times
from decision_security.survival import km_estimator, nelson_aalen
from decision_security.viz import plot_km
rng = make_rng(42)
plt.rcParams.update({
"font.family": "serif",
"font.size": 10,
"axes.labelsize": 11,
"axes.titlesize": 12,
"xtick.labelsize": 9,
"ytick.labelsize": 9,
"legend.fontsize": 9,
"figure.dpi": 150,
"axes.spines.top": False,
"axes.spines.right": False,
})
PRIMARY = "#1A1A1A"
ACCENT = "#E74C3C"
DARK_BG = "#34495E"
LIGHT_GRAY = "#95A5A6"
MED_GRAY = "#7F8C8D"
VERY_LIGHT = "#BDC3C7"
2 — Generating Survival Data#
We simulate 200 time-to-patch observations from a Weibull distribution (shape k=1.5, scale λ=30 days) with right-censoring at 60 days. Censoring means the vulnerability was still unpatched when we stopped watching.
times, events = survival_times("weibull", 200, rng=rng, k=1.5, lam=30, censor_at=60)
n_observed = int(events.sum())
n_censored = len(events) - n_observed
print("Censor at 60 days (light censoring):")
print(f" Observed (patched): {n_observed}")
print(f" Censored (unpatched): {n_censored}")
print(f" Censoring rate: {n_censored / len(events):.1%}")
times_hc, events_hc = survival_times("weibull", 200, rng=rng, k=1.5, lam=30, censor_at=30)
n_obs_hc = int(events_hc.sum())
n_cens_hc = len(events_hc) - n_obs_hc
print("\nCensor at 30 days (heavier censoring):")
print(f" Observed (patched): {n_obs_hc}")
print(f" Censored (unpatched): {n_cens_hc}")
print(f" Censoring rate: {n_cens_hc / len(events_hc):.1%}")
naive_mean_hc = times_hc[events_hc == 1].mean()
print(f"\n Naive mean of observed-only times: {naive_mean_hc:.1f} days <- biased low")
Censor at 60 days (light censoring):
Observed (patched): 193
Censored (unpatched): 7
Censoring rate: 3.5%
Censor at 30 days (heavier censoring):
Observed (patched): 122
Censored (unpatched): 78
Censoring rate: 39.0%
Naive mean of observed-only times: 15.6 days <- biased low
3 — Kaplan-Meier Curve#
The Kaplan-Meier estimator gives a non-parametric estimate of the survival function S(t) — the probability that a vulnerability remains unpatched at time t. The horizontal line at S(t) = 0.5 marks the median survival time.
from scipy.stats import norm
def km_greenwood(times, events, conf=0.95):
"""Kaplan-Meier with Greenwood CI, as step arrays (t, S, lower, upper).
Hand-rolled because km_estimator does not return the at-risk counts the CI needs."""
times = np.asarray(times, float)
events = np.asarray(events, int)
order = np.argsort(times)
times, events = times[order], events[order]
event_times = np.unique(times[events == 1])
z = norm.ppf(0.5 + conf / 2)
S, cum = 1.0, 0.0 # survival, running Greenwood sum
t_out, s_out, lo, hi = [0.0], [1.0], [1.0], [1.0]
for t in event_times:
at_risk = int((times >= t).sum())
d = int(((times == t) & (events == 1)).sum())
if at_risk == 0:
continue
S *= 1 - d / at_risk
if at_risk > d:
cum += d / (at_risk * (at_risk - d))
se = S * np.sqrt(cum)
t_out.append(t); s_out.append(S)
lo.append(max(0.0, S - z * se)); hi.append(min(1.0, S + z * se))
return np.array(t_out), np.array(s_out), np.array(lo), np.array(hi)
def km_median(t_arr, s_arr):
"""Smallest t with S(t) <= 0.5, or None if S(t) never crosses 0.5."""
below = np.where(s_arr <= 0.5)[0]
return float(t_arr[below[0]]) if len(below) else None
g_t, g_s, g_lo, g_hi = km_greenwood(times, events)
med = km_median(g_t, g_s)
fig, ax = plt.subplots(figsize=(8, 4))
ax.step(g_t, g_s, where="post", color=DARK_BG, label="KM estimate")
ax.fill_between(g_t, g_lo, g_hi, step="post", alpha=0.2, color=DARK_BG,
label="Greenwood 95% CI")
ax.axhline(0.5, color=ACCENT, linestyle="--", linewidth=0.8, label="S(t) = 0.5")
if med is not None:
ax.axvline(med, color=LIGHT_GRAY, linestyle=":", linewidth=0.8)
ax.set_xlabel("Days since disclosure")
ax.set_ylabel("S(t) — proportion unpatched")
ax.set_title("Kaplan-Meier: Time to Patch (Greenwood 95% CI)")
ax.legend()
plt.tight_layout()
plt.show()
# Under heavy censoring S(t) may never reach 0.5.
print(f"KM median: {'not reached' if med is None else f'{med:.0f} days'}"
" (consistent under non-informative censoring, unlike either mean)")
print(f"Naive mean of all recorded times: {times.mean():.1f} days")
print(f"Mean of observed (patched) times only: {times[events == 1].mean():.1f} days")
KM median: 23 days (consistent under non-informative censoring, unlike either mean)
Naive mean of all recorded times: 26.2 days
Mean of observed (patched) times only: 24.9 days
4 — Comparing Two Groups#
Risk-based patching prioritises critical vulnerabilities over medium ones. If the policy works, the two groups should have visibly different survival curves. Here we simulate “critical” (faster patching, k=2.0, λ=14) and “medium” (slower, k=1.2, λ=45) populations.
from scipy.stats import logrank, CensoredData
# One advancing RNG for both groups, not two fixed seeds.
rng_groups = make_rng(7)
t_crit, e_crit = survival_times("weibull", 200, rng=rng_groups, k=2.0, lam=14, censor_at=60)
t_med, e_med = survival_times("weibull", 200, rng=rng_groups, k=1.2, lam=45, censor_at=60)
gc_t, gc_s, gc_lo, gc_hi = km_greenwood(t_crit, e_crit)
gm_t, gm_s, gm_lo, gm_hi = km_greenwood(t_med, e_med)
fig, ax = plt.subplots(figsize=(8, 4))
ax.step(gc_t, gc_s, where="post", color=ACCENT, label="Critical (k=2.0, λ=14)")
ax.fill_between(gc_t, gc_lo, gc_hi, step="post", alpha=0.18, color=ACCENT)
ax.step(gm_t, gm_s, where="post", color=DARK_BG, label="Medium (k=1.2, λ=45)")
ax.fill_between(gm_t, gm_lo, gm_hi, step="post", alpha=0.18, color=DARK_BG)
ax.axhline(0.5, color=LIGHT_GRAY, linestyle="--", linewidth=0.8)
ax.set_xlabel("Days since disclosure")
ax.set_ylabel("S(t) — proportion unpatched")
ax.set_title("KM Curves: Critical vs. Medium (Greenwood 95% CI)")
ax.legend()
plt.tight_layout()
plt.show()
med_c = km_median(gc_t, gc_s)
med_m = km_median(gm_t, gm_s)
print("Per-group summaries:")
print(f" Critical: KM median {'not reached' if med_c is None else f'{med_c:.0f} days'}, "
f"censoring rate {1 - e_crit.mean():.1%}")
print(f" Medium: KM median {'not reached' if med_m is None else f'{med_m:.0f} days'}, "
f"censoring rate {1 - e_med.mean():.1%}")
crit_cd = CensoredData.right_censored(t_crit, e_crit == 0)
med_cd = CensoredData.right_censored(t_med, e_med == 0)
lr = logrank(crit_cd, med_cd)
print(f"\nLog-rank test (Critical vs Medium):")
print(f" z = {lr.statistic:.3f} chi2 = z^2 = {lr.statistic**2:.3f} p-value = {lr.pvalue:.2e}")
if lr.pvalue < 0.05:
print(" The curves differ significantly (p < 0.05).")
else:
print(" No significant difference (p >= 0.05).")
Per-group summaries:
Critical: KM median 12 days, censoring rate 0.0%
Medium: KM median 31 days, censoring rate 29.0%
Log-rank test (Critical vs Medium):
z = 13.861 chi2 = z^2 = 192.140 p-value = 1.08e-43
The curves differ significantly (p < 0.05).
5 — Cumulative Hazard#
The Nelson-Aalen estimator gives the cumulative hazard H(t). The hazard rate is the instantaneous risk of the event at time t, given survival to t. An increasing slope means the patch rate speeds up the longer a vulnerability stays open — which is what a Weibull shape k > 1 (here k=2.0) produces.
na_t, na_h = nelson_aalen(t_crit, e_crit)
fig, ax = plt.subplots(figsize=(8, 4))
ax.step(na_t, na_h, where="post", color=DARK_BG)
ax.set_xlabel("Days since disclosure")
ax.set_ylabel("H(t) — cumulative hazard")
ax.set_title("Nelson-Aalen: Cumulative Hazard (Critical Vulns)")
plt.tight_layout()
plt.show()
# S(t) = exp(-H(t)), so H = ln(2) marks the point where 50% are patched.
h_end, t_end = na_h[-1], na_t[-1]
idx_half = int(np.searchsorted(na_h, np.log(2)))
print(f"Day {t_end:.0f}: H(t) = {h_end:.2f}, S(t) = exp(-H) ≈ {np.exp(-h_end):.2f} still unpatched.")
if idx_half < len(na_t):
print(f"H(t) crosses ln(2) ≈ 0.69 (50% patched) around day {na_t[idx_half]:.0f}.")
Day 34: H(t) = 5.88, S(t) = exp(-H) ≈ 0.00 still unpatched.
H(t) crosses ln(2) ≈ 0.69 (50% patched) around day 12.
Pitfalls#
Censoring is information, not missing data. Dropping censored observations biases survival estimates downward — you are throwing away evidence that the event hadn’t happened yet.
Median survival time is not the mean. Means are dragged by the right tail; the median is the time at which S(t) crosses 0.5.
“Days since disclosure” is not “days of exposure.” Clock-start matters for time-to-patch metrics. A vulnerability disclosed publicly on day 0 may have been silently exploited for months. Survival analysis answers the question you ask it, just make sure the clock measures what you actually care about.
# Pitfall #1 on the heavier-censored data from Section 2.
kg_t, kg_s, *_ = km_greenwood(times_hc, events_hc)
km_med_kept = km_median(kg_t, kg_s)
dropped = times_hc[events_hc == 1]
naive_median = float(np.median(dropped))
print("Heavier-censored dataset (censor at 30 days):")
print(f" Drop censored rows -> naive median = {naive_median:.1f} days (too low)")
print(f" Keep them (KM) -> KM median = "
f"{'not reached' if km_med_kept is None else f'{km_med_kept:.0f} days'}")
Heavier-censored dataset (censor at 30 days):
Drop censored rows -> naive median = 15.1 days (too low)
Keep them (KM) -> KM median = 23 days