Healthcare capacity & clinical risk — survival analysis meets bed planning

Machine Learning in Operations Research · from time-to-event to how many beds

The staffing notebook modelled the arrival side of hospital capacity (how many patients show up). This one models the service side and the clinical risk that drives it: how long patients stay, how they leave (discharged alive vs died), who bounces back (30-day readmission), and how the resulting length-of-stay distribution sets how many beds a ward needs. It runs the arc survival analysis → competing risks → Bayesian survival → risk scoring → capacity, and it deliberately extends the mortality/longevity-forecasting toolkit into a private-sector operations problem.

Where the queueing notebook assumed a fixed service time, here we derive the length-of-stay distribution from real patient data — the exact input a capacity model needs — and then propagate it to bed occupancy. Staffing modelled demand; this models service; together they are the full picture.

Data: the UCI Diabetes 130-US hospitals set — 101,766 real inpatient encounters (1999–2008) with length of stay, discharge disposition (including in-hospital death), and 30-day readmission.

import os, re
os.environ.setdefault("KMP_DUPLICATE_LIB_OK", "TRUE")
import numpy as np, pandas as pd, matplotlib.pyplot as plt, warnings
from pathlib import Path
warnings.filterwarnings("ignore")
BLUE, RED, GREEN, ORANGE, GREY, PURP = "#2b6cb0","#c53030","#2f855a","#dd6b20","#a0aec0","#6b46c1"
here = Path.cwd(); DATA = next(p for p in [here/"data", here.parent/"data"] if (p/"diabetes.parquet").exists())
df = pd.read_parquet(DATA/"diabetes.parquet")

# --- outcomes ---
DEATH = {11, 19, 20, 21}                                   # 'Expired' discharge dispositions
df["LOS"] = df["time_in_hospital"].astype(float)           # length of stay (days), capped at 14
df["died"] = df["discharge_disposition_id"].isin(DEATH).astype(int)
df["discharged"] = 1 - df["died"]                          # 1 = alive discharge (event for LOS survival)
df["cr_event"] = np.where(df["died"]==1, 2, 1)             # competing risks: 1=discharge, 2=death
df["age_num"] = df["age"].map(lambda a: np.mean([int(x) for x in re.findall(r"\d+", str(a))]) if re.findall(r"\d+", str(a)) else np.nan)
df["readmit30"] = (df["readmitted"] == "<30").astype(int)
print(f"{len(df):,} encounters | LOS mean {df.LOS.mean():.1f}d (max {df.LOS.max():.0f}) | "
      f"in-hospital death {df.died.mean()*100:.1f}% | 30-day readmission {df.readmit30.mean()*100:.1f}%")

101,766 encounters | LOS mean 4.4d (max 14) | in-hospital death 1.6% | 30-day readmission 11.2%

1 · The data — length of stay, outcomes, readmission

Each row is an inpatient stay. Three clinically and operationally central quantities: the length of stay (how long a bed is occupied), the exit type (discharged alive vs died in hospital — competing ways to leave), and whether the patient was readmitted within 30 days.

Data dictionary — the fields we use

The record has 50 columns; below are the ones this notebook uses, grouped by what they actually measure — a distinction that matters for interpretation.

Outcome & timing

field	meaning
time_in_hospital (LOS)	length of stay, in days (1–14; the dataset caps at 14)
discharge_disposition_id	coded outcome/destination at discharge; codes 11, 19, 20, 21 = "Expired" (in-hospital death) — these define the competing death event
readmitted	<30 (readmitted within 30 days), >30, or NO; we model the <30 event
age	10-year bracket (e.g. [70-80)); we use the bracket midpoint

This-encounter activity — counts during the current stay

field	meaning
num_lab_procedures	lab tests run during the encounter
num_procedures	non-lab procedures performed
num_medications	distinct medications administered
number_diagnoses	diagnoses recorded for the encounter

Prior-year utilization — visits in the 12 months before this admission (history/comorbidity signals, not current activity)

field	meaning
number_inpatient	prior inpatient (hospital) admissions
number_emergency	prior emergency-department visits
number_outpatient	prior outpatient visits

The prior-year counts are the strongest risk signals: a patient with several recent inpatient/ER visits is chronically sicker, which surfaces both as longer stays (the Cox and Bayesian models) and higher readmission risk — the same underlying frailty showing up on the service side and the risk side. (admission_type_id — emergency/urgent/elective — and the ICD-9 diagnosis codes diag_1..3 are also in the record; we keep the models to the numeric utilization/activity fields for clarity.)

Length of stay (LOS) — the pivotal quantity

LOS is the number of days a patient occupies an inpatient bed, from admission to discharge or death (the time_in_hospital field). It sits at the centre of this notebook because it carries two meanings at once:

• Clinically, it proxies illness severity and recovery — sicker, more complex patients stay longer (and the models below confirm it).
• Operationally, it is the service time of a bed: how long each admission ties up the resource. By Little's Law, average occupancy equals admission rate times mean length of stay, $\;\overline{\text{beds}}=\lambda\,\mathbb{E}[\text{LOS}]$, so LOS is the lever for capacity — shave the average and beds free up without turning anyone away (§6 shows a 10% LOS cut freeing ~9% of beds).

That is why we treat LOS as a time-to-event outcome, where the "event" is leaving the bed: survival methods give P(still in hospital after $t$ days), competing risks separate the two ways a stay ends (discharge vs death), and the fitted LOS distribution becomes the service-time input the capacity model needs — the bridge from survival analysis to operations.

One honest caveat: time_in_hospital is capped at 14 days, so longer stays are truncated — the LOS here is biased toward shorter stays and under-represents the long-stay tail that strains capacity most. Keep that in mind wherever we quote a mean LOS or a bed count.

fig, ax = plt.subplots(1, 3, figsize=(14, 4))
ax[0].hist(df.LOS, bins=np.arange(0.5,15.5,1), color=BLUE, rwidth=.9); ax[0].set_xlabel("length of stay (days)"); ax[0].set_ylabel("encounters")
ax[0].set_title(f"Length of stay (mean {df.LOS.mean():.1f}d)")
ex = pd.Series({"discharged alive": df.discharged.sum(), "died in hospital": df.died.sum()})
ax[1].bar(["discharged\nalive","died in\nhospital"], ex.values, color=[GREEN,RED]); ax[1].set_yscale("log"); ax[1].set_title("How stays end (log scale)")
for i,v in enumerate(ex.values): ax[1].text(i,v,f"{v:,}",ha="center",va="bottom",fontsize=8)
rr = df.groupby(pd.cut(df.age_num,[0,40,55,65,75,85,100]))["readmit30"].mean()*100
ax[2].bar([str(int(i.left))+"-"+str(int(i.right)) for i in rr.index], rr.values, color=PURP); ax[2].set_xlabel("age"); ax[2].set_ylabel("% readmitted <30d")
ax[2].set_title("30-day readmission by age"); ax[2].tick_params(axis="x", rotation=30)
fig.tight_layout(); plt.show()

2 · Survival analysis of length of stay

Length of stay is a time-to-event outcome: the "event" is leaving the bed. The Kaplan-Meier estimator gives the probability a patient is still in hospital after $t$ days (treating in-hospital death as censoring for the discharge event); a Cox proportional-hazards model then says which patient factors speed up or slow down discharge — the hazard ratio $>1$ means discharged sooner (shorter stay), $<1$ means stays longer. This is the same hazard-modelling toolkit used for mortality/longevity, pointed at bed-occupancy.

from lifelines import KaplanMeierFitter, CoxPHFitter
km = KaplanMeierFitter().fit(df.LOS, df.discharged, label="still in hospital")
covs = ["age_num","number_diagnoses","num_medications","number_inpatient","number_emergency","num_procedures"]
cox_df = df[covs + ["LOS","discharged"]].dropna()
cox = CoxPHFitter().fit(cox_df, duration_col="LOS", event_col="discharged")
fig, ax = plt.subplots(1, 2, figsize=(13, 4.5))
ax[0].step(km.survival_function_.index, km.survival_function_.values[:,0], where="post", color=BLUE, lw=2)
ax[0].set_xlabel("days since admission"); ax[0].set_ylabel("P(still in hospital)"); ax[0].set_title("Kaplan-Meier: length of stay")
ax[0].grid(alpha=.3)
hr = cox.hazard_ratios_.sort_values()
ax[1].errorbar(hr.values, range(len(hr)), xerr=None, fmt="o", color=BLUE)
lo = np.exp(cox.confidence_intervals_.iloc[:,0]).reindex(hr.index); hi = np.exp(cox.confidence_intervals_.iloc[:,1]).reindex(hr.index)
ax[1].hlines(range(len(hr)), lo.values, hi.values, color=BLUE, alpha=.5)
ax[1].axvline(1, color=RED, ls="--"); ax[1].set_yticks(range(len(hr))); ax[1].set_yticklabels(hr.index, fontsize=8)
ax[1].set_xlabel("hazard ratio for discharge  (>1 = shorter stay)"); ax[1].set_title("Cox model — drivers of length of stay")
fig.tight_layout(); plt.show()
print(f"median LOS (KM): {km.median_survival_time_:.0f} days | Cox concordance {cox.concordance_index_:.3f}")
print("longest-stay drivers (HR<1):", ", ".join(f"{k} {v:.2f}" for k,v in cox.hazard_ratios_.sort_values().head(3).items()))

median LOS (KM): 4 days | Cox concordance 0.689
longest-stay drivers (HR<1): number_diagnoses 0.95, number_inpatient 0.95, num_medications 0.95

3 · Competing risks — discharged alive vs died

A patient leaves a bed in one of two competing ways: discharged alive or died in hospital. Treating death as mere "censoring" (as plain KM does) overstates the chance of discharge, because some of those "censored" patients were never going to be discharged — they died. The correct object is the cumulative incidence function (CIF): the probability of each exit type by day $t$, accounting for the other. We estimate both with the Aalen-Johansen estimator.

from lifelines import AalenJohansenFitter
aj_d = AalenJohansenFitter(calculate_variance=False).fit(df.LOS, df.cr_event, event_of_interest=1)
aj_x = AalenJohansenFitter(calculate_variance=False).fit(df.LOS, df.cr_event, event_of_interest=2)
cif_d = aj_d.cumulative_density_; cif_x = aj_x.cumulative_density_
fig, ax = plt.subplots(1, 2, figsize=(13, 4.3))
ax[0].step(cif_d.index, cif_d.values[:,0], where="post", color=GREEN, lw=2, label="discharged alive")
ax[0].step(cif_x.index, cif_x.values[:,0], where="post", color=RED, lw=2, label="died in hospital")
ax[0].set_xlabel("days since admission"); ax[0].set_ylabel("cumulative incidence"); ax[0].set_title("Competing-risks CIF (whole cohort)"); ax[0].legend(fontsize=8); ax[0].grid(alpha=.3)
# death CIF by age band (higher-risk groups)
for band, c in zip([(0,55),(55,75),(75,100)], [GREEN,ORANGE,RED]):
    m = df.age_num.between(*band)
    aj = AalenJohansenFitter(calculate_variance=False).fit(df.LOS[m], df.cr_event[m], event_of_interest=2)
    ax[1].step(aj.cumulative_density_.index, aj.cumulative_density_.values[:,0], where="post", color=c, lw=2, label=f"age {band[0]}-{band[1]}")
ax[1].set_xlabel("days since admission"); ax[1].set_ylabel("cumulative incidence of death"); ax[1].set_title("In-hospital death CIF by age"); ax[1].legend(fontsize=8); ax[1].grid(alpha=.3)
fig.tight_layout(); plt.show()
print(f"by day 14: P(discharged alive)={cif_d.values[-1,0]:.3f}, P(died)={cif_x.values[-1,0]:.3f}")

by day 14: P(discharged alive)=0.984, P(died)=0.016

4 · Bayesian survival — length of stay with quantified uncertainty

A parametric Weibull model of length of stay, fit by MCMC (NumPyro), gives a full posterior over covariate effects — credible intervals, not just point estimates — and a generative LOS distribution we can sample for capacity planning. Deaths are treated as right-censored for the discharge time (they exit the bed without being "discharged"). This is the Bayesian survival modelling from the broader portfolio, applied to bed-occupancy. (Fit on a subsample for speed.)

import jax, jax.numpy as jnp, numpyro, numpyro.distributions as dist
from numpyro.infer import MCMC, NUTS
sub = df[covs + ["LOS","discharged"]].dropna().sample(6000, random_state=0)
Xs = sub[covs].to_numpy(float); mu_, sd_ = Xs.mean(0), Xs.std(0)+1e-9; Xz = (Xs-mu_)/sd_
y = sub.LOS.to_numpy(float); delta = sub.discharged.to_numpy(float)      # 1=discharge observed, 0=censored(death)
def model(X, y, delta):
    k = numpyro.sample("k", dist.LogNormal(0., 0.5))                     # Weibull shape
    b0 = numpyro.sample("b0", dist.Normal(1., 1.))
    beta = numpyro.sample("beta", dist.Normal(0., 1.).expand([X.shape[1]]))
    lam = jnp.exp(b0 + X @ beta)                                         # per-patient scale
    d = dist.Weibull(scale=lam, concentration=k)
    ll = jnp.where(delta==1, d.log_prob(y), -(y/lam)**k)                 # logpdf if discharged, log-survival if censored
    numpyro.factor("lik", ll.sum())
mcmc = MCMC(NUTS(model), num_warmup=500, num_samples=600, num_chains=1, progress_bar=False)
mcmc.run(jax.random.PRNGKey(0), jnp.array(Xz), jnp.array(y), jnp.array(delta)); post = mcmc.get_samples()
beta = np.asarray(post["beta"])
fig, ax = plt.subplots(figsize=(9, 4.2))
med = np.median(beta,0); lo,hi = np.percentile(beta,[5,95],0); order=np.argsort(med)
ax.errorbar(med[order], range(len(covs)), xerr=[med[order]-lo[order], hi[order]-med[order]], fmt="o", color=PURP)
ax.axvline(0, color=RED, ls="--"); ax.set_yticks(range(len(covs))); ax.set_yticklabels([covs[i] for i in order], fontsize=8)
ax.set_xlabel("posterior effect on log-LOS scale  (>0 = longer stay)"); ax.set_title("Bayesian Weibull AFT — length-of-stay drivers (90% credible)")
fig.tight_layout(); plt.show()
print(f"Weibull shape k posterior median {np.median(post['k']):.2f} (>1 -> discharge hazard rises with days already stayed)")

Weibull shape k posterior median 1.71 (>1 -> discharge hazard rises with days already stayed)

5 · Clinical risk — who bounces back (30-day readmission)

Readmission within 30 days is the operational cost of premature or poorly-supported discharge. We train a gradient-boosted classifier to score readmission risk from the encounter record, and — because the class is imbalanced (~11%) — judge it by AUC and calibration rather than accuracy.

import lightgbm as lgb
from sklearn.model_selection import train_test_split
from sklearn.metrics import roc_auc_score, roc_curve
from sklearn.calibration import calibration_curve
feats = ["age_num","LOS","num_medications","number_diagnoses","number_inpatient","number_emergency",
         "number_outpatient","num_lab_procedures","num_procedures"]
M = df[feats+["readmit30"]].dropna()
Xtr,Xte,ytr,yte = train_test_split(M[feats], M["readmit30"], test_size=0.25, random_state=0, stratify=M["readmit30"])
clf = lgb.LGBMClassifier(n_estimators=300, learning_rate=0.05, num_leaves=31, verbosity=-1).fit(Xtr,ytr)
p = clf.predict_proba(Xte)[:,1]; auc = roc_auc_score(yte,p)
fig, ax = plt.subplots(1, 3, figsize=(15, 4))
fpr,tpr,_ = roc_curve(yte,p); ax[0].plot(fpr,tpr,color=BLUE,lw=2,label=f"AUC {auc:.3f}"); ax[0].plot([0,1],[0,1],"--",color=GREY)
ax[0].set_xlabel("false positive rate"); ax[0].set_ylabel("true positive rate"); ax[0].set_title("Readmission ROC"); ax[0].legend(fontsize=9)
frac,mean = calibration_curve(yte,p,n_bins=10); ax[1].plot(mean,frac,"o-",color=GREEN); ax[1].plot([0,1],[0,1],"--",color=GREY)
ax[1].set_xlabel("predicted risk"); ax[1].set_ylabel("observed rate"); ax[1].set_title("Calibration")
imp = pd.Series(clf.feature_importances_, index=feats).sort_values().tail(8)
ax[2].barh(imp.index, imp.values, color=ORANGE); ax[2].set_title("Top risk features (gain)"); ax[2].tick_params(axis="y", labelsize=8)
fig.tight_layout(); plt.show()
print(f"readmission model: AUC {auc:.3f} | base rate {yte.mean()*100:.1f}%")

readmission model: AUC 0.639 | base rate 11.2%

6 · From length of stay to beds — the capacity link

Here survival meets operations. If a ward admits patients at rate $\lambda$/day and each occupies a bed for a length of stay drawn from the distribution we just modelled, the steady-state bed occupancy of an ample-bed ward (an $M/G/\infty$ queue) is Poisson with mean $\lambda\,\mathbb{E}[\text{LOS}]$ — a beautiful result: occupancy depends on the LOS mean alone (Little's law), not its shape. The beds needed to "almost never turn a patient away" is then a high quantile of that Poisson. We validate the analytic result by simulating with the real LOS draws, then show how a length-of-stay reduction (e.g. from better discharge planning / fewer readmissions) frees beds.

from scipy.stats import poisson
E_LOS = df.LOS.mean(); LAM = 12.0                          # 12 admissions/day (illustrative ward)
rho = LAM * E_LOS                                          # offered load = mean occupancy (Little's law)
def sim_occ(lam, los_pool, days=500, seed=0):
    rng = np.random.default_rng(seed); occ = np.zeros(days+20)
    for a in range(days):
        for l in rng.choice(los_pool, rng.poisson(lam)): occ[a:a+int(np.ceil(l))] += 1
    return occ[20:days]
occ = sim_occ(LAM, df.LOS.to_numpy())
beds95 = int(poisson.ppf(0.95, rho))
fig, ax = plt.subplots(1, 2, figsize=(13, 4.3))
ax[0].hist(occ, bins=30, density=True, color=BLUE, alpha=.6, label="simulated (real LOS)")
xs = np.arange(int(occ.min()), int(occ.max())+1); ax[0].plot(xs, poisson.pmf(xs, rho), "o-", color=RED, ms=4, label=f"Poisson(λ·E[LOS]={rho:.0f})")
ax[0].axvline(beds95, color=GREEN, ls="--", label=f"95% beds = {beds95}"); ax[0].set_xlabel("beds occupied"); ax[0].set_ylabel("density")
ax[0].set_title(f"Bed occupancy (λ={LAM:.0f}/day)"); ax[0].legend(fontsize=8)
cuts = np.array([0, .05, .10, .15]); beds = [int(poisson.ppf(0.95, LAM*E_LOS*(1-c))) for c in cuts]
ax[1].bar([f"-{int(c*100)}%" for c in cuts], beds, color=GREEN)
for i,b in enumerate(beds): ax[1].text(i,b,str(b),ha="center",va="bottom",fontsize=9)
ax[1].set_xlabel("length-of-stay reduction"); ax[1].set_ylabel("beds needed (95% service)"); ax[1].set_title("Shorter stays → fewer beds")
fig.tight_layout(); plt.show()
print(f"E[LOS] {E_LOS:.2f}d, λ {LAM:.0f}/day → mean occupancy {rho:.0f} beds; 95%-service needs {beds95} beds")
print(f"a 10% LOS reduction saves {beds[0]-beds[2]} beds ({(beds[0]-beds[2])/beds[0]*100:.0f}%)")

E[LOS] 4.40d, λ 12/day → mean occupancy 53 beds; 95%-service needs 65 beds
a 10% LOS reduction saves 6 beds (9%)

7 · Takeaways

• Length of stay is a survival problem. Kaplan-Meier and Cox turn the bed-occupancy question into time-to-event modelling — the same hazard toolkit used for mortality and longevity — and identify which patients occupy beds longest.
• Competing risks matter. Patients leave a bed by discharge or death; the cumulative-incidence view (Aalen-Johansen) gives the honest probability of each, which plain survival curves distort — and death risk rises sharply with age.
• Bayesian survival adds calibrated uncertainty on those drivers and a generative LOS distribution for planning.
• Clinical risk scoring (readmission) targets the interventions that shorten future stays.
• The capacity link closes the loop: the LOS distribution feeds an $M/G/\infty$ occupancy model, so bed requirements follow directly — and because occupancy depends on the LOS mean, a modest length-of-stay reduction (better discharge planning, fewer readmissions) translates straight into beds freed. This is the service side of the same capacity picture the staffing notebook approached from the demand side — survival analysis supplying the service-time distribution that queueing models need.

A fifth ML in Operations Research family, and a bridge from the portfolio's survival/Bayesian work to operations: predict time-to-event with uncertainty, then optimize the capacity decision it implies.