4.2 - Naive Bayes classifier for emergency triage

The second algorithm the Minermont team studies for triage is the Naive Bayes classifier, a probabilistic method based on Bayes' theorem. Martha's notes (section 4.5) develop its mathematics; in this notebook we implement it by hand, compare it with scikit-learn and see where its "naive" assumption helps and where it hurts:

1. Import the libraries.
2. Bayes' theorem and the naive assumption.
3. The triage data.
4. A prediction by hand, step by step.
5. Gaussian Naive Bayes, and why it fails here.
6. Counting frequencies: Naive Bayes with categories.
7. When the naive assumption breaks.
8. Final evaluation.

Step 1: Import the libraries

---

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn.model_selection import train_test_split
from sklearn.naive_bayes import GaussianNB, CategoricalNB
from sklearn.metrics import accuracy_score, recall_score, classification_report, ConfusionMatrixDisplay

Step 2: Bayes' theorem and the naive assumption

---

For a patient with features $X = (x_1, \dots, x_d)$ we want the probability of each urgency level $C_k$. Bayes' theorem gives it:

$$P(C_k \mid X) = \frac{P(X \mid C_k)\, P(C_k)}{P(X)} \qquad\Longrightarrow\qquad P(C_k \mid X) \propto P(X \mid C_k)\, P(C_k)$$

• $P(C_k)$ is the prior: how frequent the level is before looking at the patient.
• $P(X \mid C_k)$ is the likelihood: how probable these features are for a patient of that level.
• $P(X)$ is the same for every level, so to choose the level we can ignore it.

Estimating $P(x_1, \dots, x_d \mid C_k)$ for every combination of values would need an enormous amount of data. The naive assumption is that, once the level is known, the features are independent of each other, so the likelihood becomes a product:

$$P(X \mid C_k) = \prod_{j=1}^{d} P(x_j \mid C_k) \qquad\qquad \text{Predicted level} = \arg\max_k \; P(C_k) \prod_{j=1}^{d} P(x_j \mid C_k)$$

Step 3: The triage data

---

The emergency records of a real hospital are confidential, so we generate synthetic patients that behave like the ones in the book. The next cell contains the generator. You do not need to follow it line by line: what matters is that the urgency level comes from a hidden protocol (an early-warning score based on the vital signs, plus some clinical warning signs) and from the nurse's own judgement, which sometimes moves the level one step. The models never see that rule: they have to discover it from the examples.

Column	Meaning
age	Age in years (1–97; about one patient in eight is a child)
male	1 if the patient is a man, 0 if a woman
systolic_bp	Systolic blood pressure (mmHg)
heart_rate	Heart rate (beats per minute)
resp_rate	Respiratory rate (breaths per minute)
spo2	Oxygen saturation (%)
temperature	Temperature (°C)
consciousness	Level of consciousness: 0 alert, 1 responds to voice, 2 responds only to pain, 3 unresponsive
pain	Pain reported by the patient (0–10)
chest_pain, dyspnea, trauma, vomiting	Presenting symptoms: chest pain, difficulty breathing, injury, vomiting (1 = yes)
level	Target: urgency level, from 1 (red, immediate) to 5 (blue, non-urgent)

def generate_patients(n, seed):
    """Synthetic emergency patients at Minermont: one row per patient, with the urgency level assigned at triage."""
    rng = np.random.default_rng(seed)
    age = np.clip(np.round(np.where(rng.random(n) < 0.12, rng.uniform(1, 17, n), rng.normal(52, 20, n))), 1, 97)
    male = rng.binomial(1, 0.5, n)
    s = rng.gamma(0.9, 0.55, n)                        # hidden severity of the episode (never shown to the model)
    infection = rng.binomial(1, 0.12 + 0.1 * (s > 1), n)
    trauma = rng.binomial(1, 0.2, n)
    chest_pain = rng.binomial(1, np.clip(0.04 + 0.0025 * age, 0, 0.35), n)
    dyspnea = rng.binomial(1, np.clip(0.06 + 0.25 * s, 0, 0.9), n)
    vomiting = rng.binomial(1, 0.12, n)
    shock = (rng.random(n) < 0.35) & (s > 1.2)
    hr = 76 + 16 * s + 12 * infection + rng.normal(0, 10, n)
    brady = (rng.random(n) < 0.08) & (s > 1.3)         # a few severe patients have a very slow heart rate
    hr = np.where(brady, rng.normal(42, 6, n), hr)
    rr = 15 + 4.5 * s + 3 * dyspnea + rng.normal(0, 2.2, n)
    spo2 = np.clip(98.2 - 2.2 * s * (0.5 + dyspnea) + rng.normal(0, 1.1, n), 70, 100)
    sbp = np.clip(128 + 0.25 * (age - 50) + rng.normal(0, 16, n) - 28 * s * shock, 60, 230)
    temp = 36.8 + rng.normal(0, 0.35, n) + infection * rng.normal(1.6, 0.6, n)
    p_consc = np.clip(0.02 * s ** 2.2, 0, 0.9)
    consc = np.where(rng.random(n) < p_consc, rng.choice([1, 2, 3], n, p=[0.55, 0.3, 0.15]), 0)
    pain = np.clip(np.round(rng.normal(2.5 + 3.5 * trauma + 3 * chest_pain + 0.8 * s, 1.8, n)), 0, 10).astype(int)
    hr, rr, sbp = np.round(hr).astype(int), np.round(rr).astype(int), np.round(sbp).astype(int)
    temp, spo2 = np.round(temp, 1), np.round(spo2).astype(int)

    # The hidden protocol behind the labels, used only to generate the data:
    # an early-warning score in which a vital sign scores points when it is too low OR too high...
    score = (np.select([rr <= 8, rr <= 11, rr <= 20, rr <= 24], [3, 1, 0, 2], 3)
             + np.select([spo2 <= 91, spo2 <= 93, spo2 <= 95], [3, 2, 1], 0)
             + np.select([sbp <= 90, sbp <= 100, sbp <= 110, sbp <= 219], [3, 2, 1, 0], 3)
             + np.select([hr <= 40, hr <= 50, hr <= 90, hr <= 110, hr <= 130], [3, 1, 0, 1, 2], 3)
             + np.where(consc > 0, 3, 0)
             + np.select([temp <= 35, temp <= 36, temp <= 38, temp <= 39], [3, 1, 0, 1], 2))
    # ...plus clinical discriminators, from the least to the most urgent level
    level = np.full(n, 5)
    level = np.where((score >= 1) | (pain >= 4) | (trauma == 1), 4, level)
    level = np.where((score >= 4) | (pain >= 7) | (temp >= 39) | ((chest_pain == 1) & (age < 45)), 3, level)
    level = np.where((score >= 6) | ((chest_pain == 1) & (age >= 45)) | (spo2 <= 91), 2, level)
    level = np.where((score >= 9) | (consc >= 2) | (sbp <= 80), 1, level)
    # and the nurse's own judgement, which sometimes moves the level one step
    level = np.clip(level + rng.choice([-1, 0, 1], n, p=[0.05, 0.9, 0.05]), 1, 5)

    return pd.DataFrame({
        'age': age.astype(int), 'male': male, 'systolic_bp': sbp, 'heart_rate': hr, 'resp_rate': rr,
        'spo2': spo2, 'temperature': temp, 'consciousness': consc, 'pain': pain, 'chest_pain': chest_pain,
        'dyspnea': dyspnea, 'trauma': trauma, 'vomiting': vomiting, 'level': level,
    })

df = generate_patients(4000, seed=42)
feature_names = [c for c in df.columns if c != 'level']
level_names = ['1 Red', '2 Orange', '3 Yellow', '4 Green', '5 Blue']
X = df[feature_names]
y = df['level']

X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.25, random_state=42, stratify=y)
print(f"Training set: {len(y_train)} patients, test set: {len(y_test)} patients\n")
print("Patients per urgency level:")
print(pd.DataFrame({'patients': y.value_counts().sort_index(),
                    'share': y.value_counts(normalize=True).sort_index().round(3)},).set_axis(level_names).to_string())

Training set: 3000 patients, test set: 1000 patients

Patients per urgency level:
          patients  share
1 Red          158  0.040
2 Orange       530  0.132
3 Yellow       780  0.195
4 Green       1671  0.418
5 Blue         861  0.215

Step 4: A prediction by hand

---

Before using all the features, let's do the calculation with only two symptoms. A patient arrives with chest pain and difficulty breathing. From the training data we estimate, for each level, the prior and the two likelihoods, exactly as the book describes: the number of patients of the level with the symptom, divided by the number of patients of the level.

prior = y_train.value_counts(normalize=True).sort_index()
p_chest = X_train.groupby(y_train)['chest_pain'].mean()      # P(chest pain | level)
p_dysp = X_train.groupby(y_train)['dyspnea'].mean()        # P(dyspnea | level)
product = prior * p_chest * p_dysp                            # numerator of Bayes' theorem

example = pd.DataFrame({'P(level)': prior, 'P(chest pain | level)': p_chest, 'P(dyspnea | level)': p_dysp,
                        'product': product, 'P(level | symptoms)': product / product.sum()})
print(example.set_axis(level_names).round(3).to_string())

          P(level)  P(chest pain | level)  P(dyspnea | level)  product  P(level | symptoms)
1 Red        0.039                  0.280               0.432    0.005                0.116
2 Orange     0.133                  0.761               0.271    0.027                0.671
3 Yellow     0.195                  0.202               0.212    0.008                0.204
4 Green      0.418                  0.006               0.154    0.000                0.009
5 Blue       0.215                  0.000               0.088    0.000                0.000

Level 2 (orange) is by far the most probable (0.671): chest pain is very common among Level 2 patients. Level 1 gets only 0.116, even though difficulty breathing is more frequent there than in any other level, because it is a rare level (its prior is 0.039). And look at Level 5: no Level 5 patient in the training set had chest pain, so its product is exactly zero and the level is ruled out completely. We will come back to these zeros.

Step 5: Gaussian Naive Bayes

---

For numerical features, such as the heart rate, the book explains that we usually assume a normal (Gaussian) distribution within each level: we estimate its mean $\mu$ and variance $\sigma^2$, and the likelihood of a value is the height of the bell curve,

$$P(x_j \mid C_k) = \frac{1}{\sqrt{2\pi\sigma_{kj}^2}} \exp\left(-\frac{(x_j - \mu_{kj})^2}{2\sigma_{kj}^2}\right)$$

Multiplying 13 small numbers gives numbers so tiny that the computer rounds them to zero, so in practice we add logarithms instead of multiplying probabilities. The maximum does not change. To turn the scores back into probabilities we divide by their sum, which is the denominator $P(X)$ of Bayes' theorem.

levels = np.sort(y_train.unique())
Xtr = X_train.to_numpy(dtype=float)
ytr = y_train.to_numpy()

log_prior = np.log([np.mean(ytr == c) for c in levels])
means = np.array([Xtr[ytr == c].mean(axis=0) for c in levels])                 # levels x features
epsilon = 1e-9 * Xtr.var(axis=0).max()                                         # the tiny safety margin scikit-learn adds
variances = np.array([Xtr[ytr == c].var(axis=0) for c in levels]) + epsilon

def gaussian_log_likelihood(X_new, columns):
    j = [feature_names.index(c) for c in columns]
    x = np.asarray(X_new[columns], dtype=float)[:, None, :]                   # patients x 1 x features
    ll = -0.5 * np.log(2 * np.pi * variances[:, j]) - (x - means[:, j]) ** 2 / (2 * variances[:, j])
    return ll.sum(axis=2)                                                      # patients x levels

def to_probabilities(log_scores):
    p = np.exp(log_scores - log_scores.max(axis=1, keepdims=True))
    return p / p.sum(axis=1, keepdims=True)                                    # divide by P(X)

gauss_scores = log_prior + gaussian_log_likelihood(X_test, feature_names)
pred_gauss = levels[gauss_scores.argmax(axis=1)]

sk_gauss = GaussianNB().fit(X_train, y_train)
print(f"Same prediction as scikit-learn's GaussianNB: {np.sum(pred_gauss == sk_gauss.predict(X_test))} of {len(y_test)}")
print(f"Test accuracy: {accuracy_score(y_test, pred_gauss):.3f}")
print("Recall per level:", dict(zip(level_names, recall_score(y_test, pred_gauss, average=None).round(2).tolist())))

Same prediction as scikit-learn's GaussianNB: 1000 of 1000
Test accuracy: 0.477
Recall per level: {'1 Red': 0.65, '2 Orange': 0.67, '3 Yellow': 0.25, '4 Green': 0.25, '5 Blue': 0.98}

Both implementations agree on every patient, but the result is poor: an accuracy of 0.477, much lower than KNN. The recall per level shows the problem: almost every Level 5 patient is found (0.98), but only a quarter of the Level 3 and Level 4 patients. Why? Look at what the model has learned about chest pain:

j = feature_names.index('chest_pain')
print(pd.DataFrame({'mean = P(chest pain)': means[:, j], 'variance': variances[:, j]}, index=level_names).to_string(float_format='{:.3g}'.format))

          mean = P(chest pain)  variance
1 Red                     0.28     0.201
2 Orange                 0.761     0.182
3 Yellow                 0.202     0.161
4 Green                0.00559   0.00556
5 Blue                       0  5.36e-07

Chest pain is a yes/no column, and no Level 5 patient has it. The Gaussian of Level 5 therefore has a mean of 0 and a variance of practically zero: a spike, enormously tall at 0. Every patient without chest pain, which is most of them, receives a huge likelihood for Level 5 from this single column, and it outweighs everything else. The problem is not Bayes' theorem but how the feature is described: a bell curve is not a sensible model for a symptom that is either present or absent.

Step 6: Counting frequencies

---

For categorical features, the book estimates the likelihood by counting: the number of patients of the level with that value, divided by the number of patients of the level. We can use the same idea for the vital signs if we first cut each one into intervals. We use 6 intervals per measurement, with the same number of training patients in each (the cut points are computed on the training set only).

One detail: if a value never appears in a level, its frequency is 0, and a single zero in the product would rule the level out completely. Laplace smoothing adds one imaginary patient to every count ($\alpha = 1$), so no probability is exactly zero:

$$P(x_j = v \mid C_k) = \frac{\text{patients of level } k \text{ with } x_j = v \;+\; \alpha}{\text{patients of level } k \;+\; \alpha \cdot (\text{number of possible values of } x_j)}$$

continuous = ['age', 'systolic_bp', 'heart_rate', 'resp_rate', 'spo2', 'temperature']
categorical = [c for c in feature_names if c not in continuous]
cut_points = {c: np.unique(np.quantile(X_train[c], np.linspace(0, 1, 7)[1:-1])) for c in continuous}

def to_intervals(X_):
    B = X_.copy()
    for c in continuous:
        B[c] = np.digitize(X_[c], cut_points[c])      # 0, 1, 2, ... = interval number
    return B.astype(int)

B_train, B_test = to_intervals(X_train), to_intervals(X_test)
n_values = {c: len(cut_points[c]) + 1 for c in continuous}
n_values.update({'male': 2, 'consciousness': 4, 'pain': 11, 'chest_pain': 2, 'dyspnea': 2, 'trauma': 2, 'vomiting': 2})

def count_tables(B, y, alpha=1):
    """P(value | level) for every feature, by counting, with Laplace smoothing."""
    tables = {}
    for col in B.columns:
        t = np.zeros((len(levels), n_values[col]))
        for k, c in enumerate(levels):
            counts = np.bincount(B[col][y == c], minlength=n_values[col])
            t[k] = (counts + alpha) / (counts.sum() + alpha * n_values[col])
        tables[col] = t
    return tables

tables = count_tables(B_train, ytr)

def count_log_likelihood(B, columns):
    return sum(np.log(tables[col][:, B[col].to_numpy()]).T for col in columns)   # patients x levels

print("Heart rate intervals (bpm):", cut_points['heart_rate'])
print(pd.DataFrame(tables['heart_rate'], index=level_names,
                   columns=[f'interval {i}' for i in range(n_values['heart_rate'])]).round(3).to_string())

Heart rate intervals (bpm): [72. 78. 84. 90. 98.]
          interval 0  interval 1  interval 2  interval 3  interval 4  interval 5
1 Red          0.089       0.032       0.065       0.081       0.153       0.581
2 Orange       0.141       0.121       0.178       0.151       0.158       0.250
3 Yellow       0.112       0.110       0.161       0.159       0.210       0.249
4 Green        0.147       0.151       0.177       0.172       0.215       0.137
5 Blue         0.209       0.253       0.255       0.213       0.049       0.021

The table shows what a bell curve cannot: the probability of each heart rate interval within each level, without assuming any shape. Most Level 1 patients are in the fastest interval (0.581), but for them the slowest interval (0.089) is more frequent than any of the intermediate ones: some critical patients have a heart rate that is too slow. Level 5 patients, instead, are concentrated below 90 bpm.

We compare three versions: all Gaussian; the book's recipe (Gaussian for the six measurements and counting for the categorical features); and counting for everything, with the measurements cut into intervals. We also check the last one against scikit-learn's CategoricalNB.

versions = {
    'all Gaussian': gauss_scores,
    'Gaussian + counting (book)': log_prior + gaussian_log_likelihood(X_test, continuous) + count_log_likelihood(X_test, categorical),
    'counting with intervals': log_prior + count_log_likelihood(B_test, feature_names),
}
rows = []
for name, scores in versions.items():
    pred = levels[scores.argmax(axis=1)]
    rows.append({'version': name, 'accuracy': accuracy_score(y_test, pred),
                 **dict(zip([f'recall {n}' for n in level_names], recall_score(y_test, pred, average=None)))})
print(pd.DataFrame(rows).set_index('version').round(3).to_string())

pred_counts = levels[versions['counting with intervals'].argmax(axis=1)]
sk_cat = CategoricalNB(alpha=1, min_categories=[n_values[c] for c in feature_names]).fit(B_train, y_train)
print(f"\nSame prediction as scikit-learn's CategoricalNB: {np.sum(pred_counts == sk_cat.predict(B_test))} of {len(y_test)}")

                            accuracy  recall 1 Red  recall 2 Orange  recall 3 Yellow  recall 4 Green  recall 5 Blue
version                                                                                                            
all Gaussian                   0.477          0.65            0.667            0.246           0.249          0.981
Gaussian + counting (book)     0.742          0.65            0.697            0.605           0.801          0.795
counting with intervals        0.761          0.65            0.667            0.513           0.904          0.786

Same prediction as scikit-learn's CategoricalNB: 1000 of 1000

The book's recipe, a Gaussian only for the six measurements and counting for the categorical features, jumps from 0.477 to 0.742. Counting everything with intervals goes a little further, to 0.761, because the intervals capture shapes that a bell curve cannot, like the slow heart rates of some critical patients. Our implementation and scikit-learn's CategoricalNB give the same prediction for all 1000 test patients.

Step 7: When the naive assumption breaks

---

Grace objects in the book that symptoms are not independent. In our patients, the respiratory rate, the oxygen saturation and the difficulty breathing all come from the same underlying problem. Within a single level they are clearly correlated:

print("Correlations among Level 2 patients (training set):")
print(X_train[y_train == 2][['resp_rate', 'spo2', 'dyspnea', 'heart_rate']].corr().round(2).to_string())

Correlations among Level 2 patients (training set):
            resp_rate  spo2  dyspnea  heart_rate
resp_rate        1.00 -0.68     0.59        0.37
spo2            -0.68  1.00    -0.61       -0.28
dyspnea          0.59 -0.61     1.00        0.22
heart_rate       0.37 -0.28     0.22        1.00

Naive Bayes multiplies their likelihoods as if they were three separate pieces of evidence, when to a large extent they tell the same story. The effect is easiest to see with an extreme case: we give the model the same oxygen saturation column several times. A copy adds no information at all, but Naive Bayes counts it again. Counting the same evidence $n$ times is the same as adding its log-likelihood $n$ times:

base_scores = versions['counting with intervals']
spo2_term = count_log_likelihood(B_test, ['spo2'])

for copies in [1, 3, 5]:
    scores = base_scores + (copies - 1) * spo2_term
    pred = levels[scores.argmax(axis=1)]
    confidence = to_probabilities(scores).max(axis=1).mean()
    print(f"SpO2 counted x{copies}: accuracy {accuracy_score(y_test, pred):.3f}, "
          f"mean probability of the chosen level {confidence:.3f}")

i = np.where((y_test.to_numpy() == 2) & (B_test['spo2'].to_numpy() == 0))[0][0]
print("\nA Level 2 patient with low saturation:")
print(X_test.iloc[[i]].to_string(index=False))
for copies in [1, 2, 3, 5]:
    p = to_probabilities(base_scores[[i]] + (copies - 1) * spo2_term[[i]])[0]
    print(f"  SpO2 counted x{copies}: " + ", ".join(f"{n} {v:.3f}" for n, v in zip(level_names, p)))

SpO2 counted x1: accuracy 0.761, mean probability of the chosen level 0.710
SpO2 counted x3: accuracy 0.725, mean probability of the chosen level 0.737
SpO2 counted x5: accuracy 0.701, mean probability of the chosen level 0.759

A Level 2 patient with low saturation:
 age  male  systolic_bp  heart_rate  resp_rate  spo2  temperature  consciousness  pain  chest_pain  dyspnea  trauma  vomiting
  85     1          144          39         23    94         36.1              0     2           0        1       0         1
  SpO2 counted x1: 1 Red 0.320, 2 Orange 0.179, 3 Yellow 0.204, 4 Green 0.286, 5 Blue 0.011
  SpO2 counted x2: 1 Red 0.644, 2 Orange 0.150, 3 Yellow 0.136, 4 Green 0.070, 5 Blue 0.000
  SpO2 counted x3: 1 Red 0.847, 2 Orange 0.082, 3 Yellow 0.059, 4 Green 0.011, 5 Blue 0.000
  SpO2 counted x5: 1 Red 0.976, 2 Orange 0.016, 3 Yellow 0.008, 4 Green 0.000, 5 Blue 0.000

Each extra copy makes the model more confident and less accurate: the accuracy falls from 0.761 to 0.701 while the mean probability of the chosen level rises from 0.710 to 0.759. For the Level 2 patient, the probability of Level 1 goes from 0.320 to 0.976 just by repeating the same measurement. The correlated vital signs of real patients produce a milder version of the same effect: Naive Bayes can be very sure and still be wrong.

Step 8: Final evaluation

---

The best version is the one that counts frequencies with intervals. This is its detailed evaluation on the test set:

print(classification_report(y_test, pred_counts, target_names=level_names, digits=3))
ConfusionMatrixDisplay.from_predictions(y_test, pred_counts, display_labels=level_names, cmap='Blues', xticks_rotation=30)
plt.xlabel('Predicted level')
plt.ylabel('Real level')
plt.title('Naive Bayes (counting with intervals) on the test set')
plt.show()

              precision    recall  f1-score   support

       1 Red      0.406     0.650     0.500        40
    2 Orange      0.693     0.667     0.680       132
    3 Yellow      0.714     0.513     0.597       195
     4 Green      0.789     0.904     0.843       418
      5 Blue      0.889     0.786     0.835       215

    accuracy                          0.761      1000
   macro avg      0.698     0.704     0.691      1000
weighted avg      0.768     0.761     0.758      1000

Naive Bayes reaches 0.761 on the test set, better than KNN, with a model that trains in a fraction of a second. But look at the red level: it finds 26 of the 40 Level 1 patients (recall 0.650), and when it says Level 1 it is right less than half of the time (precision 0.406).

Conclusion

---

• Naive Bayes combines the prior of each level with the likelihood of each feature, using Bayes' theorem and the assumption that the features are independent within each level.
• Training it is just computing frequencies, means and variances: it is very fast and needs little data. We add logarithms instead of multiplying, and Laplace smoothing avoids zero probabilities.
• How we model each feature matters a lot: a Gaussian for a yes/no symptom, or for a vital sign that is dangerous both when it is too low and too high, is a bad description, and counting frequencies works much better.
• When several features tell the same story, Naive Bayes counts it several times and becomes overconfident. That is exactly Grace's objection, and the reason why, as Victor concludes in the book, it is probably not the most suitable model for the complexity of triage. The next notebook studies decision trees.