Predictive Analytics for Transfusion Need¶
Setup¶
This notebook connects to a local MySQL database (BloodMatching) containing
patient, vitals, and lab tables. Credentials are loaded from a .env file
(not committed to version control — see .gitignore) rather than hardcoded,
so this notebook can be shared publicly without exposing database access.
To run this notebook yourself: create a .env file in this directory with
DB_HOST, DB_USER, DB_PASSWORD, and DB_NAME.
Project Overview¶
A small-data ML pipeline predicting blood transfusion need from routine patient vitals (heart rate, respiratory rate, oxygen saturation, blood pressure, sex) — built on a 206-row MIMIC-style demo dataset (14 positive cases)
Pipeline: SQL ingestion → data cleaning & aggregation → leakage-safe feature selection → SMOTE-balanced cross-validation (Random Forest vs. Logistic Regression) → threshold & calibration analysis → client-side model deployment.
Dataset¶
The Medical Information Mart for Intensive Care (MIMIC)-IV database is comprised of deidentified electronic health records for patients admitted to the Beth Israel Deaconess Medical Center https://physionet.org/content/mimic-iv-demo/2.2/
| Attribute | Value |
|---|---|
| Total rows (patient stays) | 206 |
Positive cases (needs_transfusion = 1) |
14 (~6.8%) |
| Negative cases | 192 (~93.2%) |
| Features used in final model | 6 |
| Source tables | patients, vital, lab_result |
| Database | MySQL (BloodMatching) |
Features(6, after leakage removal)¶
| Feature | Description | Type |
|---|---|---|
| Gender | Patient sex, label-encoded | Categorical |
| heart_rate_mean | Mean heart rate across the stay | Continuous |
| resprate_mean | Mean respiratory rate across the stay | Continuous |
| o2sat_mean | Mean oxygen saturation across the stay | Continuous |
| sbp_mean | Mean systolic blood pressure across the stay | Continuous |
| dbp_mean | Mean diastolic blood pressure across the stay | Continuous |
Target¶
needs_transfusion is a rule-derived synthetic label (hemoglobin < 7.5 OR (sbp_min < 90 AND heart_rate_max > 100)), not a real clinical transfusion-order event.
hemoglobin, sbp_min, and heart_rate_max are excluded from the feature set since they directly define the label.
Methodology¶
Ingestion — patients, vitals, and lab data loaded into a local MySQL database via SQLAlchemy, credentials handled through environment variables.
Cleaning & aggregation — vitals aggregated from repeated per-reading rows to one row per stay (fixing a many-to-many join bug in the process); numeric coercion, sentinel-value handling, and median imputation for missing values.
Leakage-safe feature selection — label-defining columns and everything derived from them excluded from X, after two separate rounds of leakage were traced and fixed.
Class imbalance handling — SMOTE applied inside each cross-validation fold via an imblearn.Pipeline, avoiding synthetic-sample leakage across folds.
Model comparison — Random Forest vs. Logistic Regression, both evaluated on identical stratified 5-fold splits, reporting mean ROC-AUC ± standard deviation rather than a single score.
Threshold & calibration — out-of-fold probabilities used to select an operating threshold via precision-recall curve (favoring recall, since missing a transfusion need is costlier than a false alarm), with a calibration curve checking whether predicted probabilities are trustworthy.
Deployment — final model refit on the full dataset, exported to JavaScript via m2cgen, running entirely client-side in the interactive demo.
Key Findings¶
| # | Finding | Impact |
|---|---|---|
| 1 | Label leakage found and fixed twice — once as a pandas formula, once baked into the SQL view itself | Both produced a meaningless perfect ROC-AUC of 1.0 until corrected |
| 2 | Database join-explosion bug: two one-to-many tables joined without a shared row key | Inflated 222 patients into 2M+ rows before being caught |
| 3 | Reproducibility bug: unseeded random hemoglobin generator | Dataset labels silently changed on every re-run until fixed |
| 4 | Logistic Regression matched/outperformed Random Forest in CV | Legitimate small-data finding, not a failed experiment |
| 5 | Real model shows non-monotonic behavior on heart rate | Concrete evidence of small-sample overfitting, disclosed openly in the demo |
| 6 | Final CV metric | [fill in mean ROC-AUC ± std once your reproducible run is locked] |
import sys
import pandas as pd
from sqlalchemy import create_engine
import numpy as np
import mysql.connector as db
import os
import matplotlib.pyplot as plt
import seaborn as sns
import plotly.graph_objects as go
from plotly.subplots import make_subplots
from dotenv import load_dotenv
from sklearn.preprocessing import LabelEncoder
from sklearn.model_selection import train_test_split
from sklearn.ensemble import RandomForestClassifier
from sklearn.linear_model import LogisticRegression
from sklearn.metrics import classification_report, confusion_matrix, roc_auc_score , precision_recall_curve
from sklearn.model_selection import StratifiedKFold, cross_val_score , cross_val_predict
from sklearn.calibration import calibration_curve
from imblearn.over_sampling import SMOTE
from imblearn.pipeline import Pipeline as ImbPipeline
import m2cgen as m2c
print(os.getcwd())
print(sys.executable)
found = load_dotenv()
print('Loaded:',found)
import os
from dotenv import load_dotenv
load_dotenv() # reads the .env file into environment variables
con = db.connect(
host=os.getenv("DB_HOST"),
user=os.getenv("DB_USER"),
password=os.getenv("DB_PASSWORD"),
database=os.getenv("DB_NAME")
)
engine = create_engine(
f'mysql+mysqlconnector://{os.getenv("DB_USER")}:{os.getenv("DB_PASSWORD")}@{os.getenv("DB_HOST")}/{os.getenv("DB_NAME")}'
)
Data Ingestion¶
Source data (patients.csv, vitals, lab values) is loaded into a local MySQL
database. Tables are written with if_exists='replace' rather than 'append' —
this makes re-running the notebook idempotent.
An earlier version of this pipeline
used 'append', which silently duplicated rows on every re-run and inflated the
dataset from ~222 patients to over 2 million rows before this was caught and fixed.
df_patients = pd.read_csv(r'/Users/samikshatripathy/Downloads/patients (1).csv')
df_patients = pd.DataFrame(df_patients)
df_patients.head()
df_patients = df_patients[['subject_id', 'stay_id','gender','race','intime']]
df_patients['intime'] = pd.to_datetime(df_patients['intime'], format='mixed')
df_patients.to_sql('patients',con = engine,if_exists = 'replace',index = False)
print(df_patients.shape)
print(df_patients['subject_id'].nunique())
print(df_patients['stay_id'].nunique())
dupe_ids = df_patients[df_patients.duplicated('subject_id', keep=False)]
print(dupe_ids.sort_values('subject_id').head(20))
df_patients.info()
df_vitals = pd.read_csv(r"/Users/samikshatripathy/Downloads/vital_sign.csv")
df_vitals = pd.DataFrame(df_vitals)
df_vitals = df_vitals.rename(columns={"heartrate": "heart_rate"})
df_vitals.head()
df_vitals.info
Why vitals are aggregated to one row per stay¶
stay_id is a primary key in patients (one row per stay) but a foreign key in
vital, where it legitimately repeats — one row per vitals reading, several
readings per stay. That's normal relational structure.
The bug: lab_result was originally built by copying vital's row structure
directly, giving it the same one-to-many repetition. When the view joined two
tables that were both one-to-many against the same stay_id, SQL paired every
matching row on both sides — a many-to-many join that multiplied ~222 patients
into millions of nonsensical rows (a vitals reading from 2am paired with a lab
value that was really copied from an unrelated 6pm reading).
Fix: both vital and lab_result are aggregated to exactly one row per
stay_id before the join, so each row in the final dataset represents one
genuine patient encounter — which is also the right unit for this prediction
task in the first place.
df_vitals_agg = df_vitals.groupby('stay_id').agg(
heart_rate_mean=('heart_rate', 'mean'),
heart_rate_max=('heart_rate', 'max'),
resprate_mean=('resprate', 'mean'),
o2sat_mean=('o2sat', 'mean'),
sbp_mean=('sbp', 'mean'),
sbp_min=('sbp', 'min'),
dbp_mean=('dbp', 'mean')
).reset_index()
print(df_vitals_agg.shape)
print(df_vitals_agg.columns.tolist())
df_vitals_agg.to_sql("vital",con=engine,if_exists ="replace",index =False)
def stimulate_hb(sbp, rng=np.random.default_rng(42)):
if sbp < 90:
return round(rng.uniform(6.0, 8.5), 1)
else:
return round(rng.uniform(11.0, 15.0), 1)
df_lab = df_vitals_agg[["stay_id", "sbp_min"]].copy()
df_lab["hemoglobin"] = df_lab['sbp_min'].apply(stimulate_hb)
df_lab[['stay_id', 'hemoglobin']].to_sql('lab_result', con=engine, if_exists='replace', index=False)
print(df_lab.shape)
Data Cleaning¶
Missing value handling follows these steps:
- Coerce to numeric —
pd.to_numeric(..., errors='coerce')turns any unparseable values intoNaN. - Sentinel values → NaN — zero/negative readings are physiologically impossible for vitals like heart rate or blood pressure, so they're treated as missing rather than real measurements. This is applied only to continuous vitals/lab columns, never to categorical or label columns (where 0 can be a valid value, not an error).
- Median imputation — missingness here is low (under 6% per column across ~638k source readings before aggregation), so simple median fill is a reasonable choice. Higher missingness would call for dropping the column or adding a missingness-indicator flag instead.
Gender is handled separately via label encoding, not the numeric-coercion
pipeline, since it's categorical, not continuous.
df = pd.read_sql("select * from vw_transfusion_input",engine)
print(df.head())
df.shape
print(df.isna().sum())
print(df.info())
df.dtypes
print(df['stay_id'].duplicated().sum())
print(df.columns.tolist())
print(df.shape)
print(df[["heart_rate_mean","resprate_mean","o2sat_mean","sbp_min","dbp_mean"]].isna().mean()*100)
numeric_col = df[["heart_rate_mean","resprate_mean","o2sat_mean","sbp_min","dbp_mean"]]
for col in numeric_col:
pd.to_numeric(df[col])
# Relacing the zeros/negatives with NAN values
for col in numeric_col:
df.loc[df[col]<=0,col] = np.nan
for col in numeric_col:
df[col] = df[col].fillna(df[col].median())
print(df.isna().sum())
# no NAN vales
print(df['Gender'].unique())
lb = LabelEncoder()
df['Gender'] = lb.fit_transform(df['Gender'].astype(str))
Label Construction and Data Leakage¶
needs_transfusion is defined by a rule in the vw_transfusion_input SQL view:
CASE WHEN (hemoglobin < 7.5) OR (sbp_min < 90 AND heart_rate_max > 100)
THEN 1 ELSE 0 ENDThis is a synthetic, rule-derived label for demo purposes — not a real clinical transfusion-order event. Two rounds of data leakage were found and fixed while building this pipeline:
- First leakage: an early version of this notebook additionally overwrote
needs_transfusionin pandas using a formula built directly fromhemoglobinandsbp— the same columns used as model features. This let the model trivially "predict" the label by recovering the formula, producing a meaningless perfect score. - Second leakage (same root cause, one layer deeper): after removing the
pandas formula, the model still scored a perfect 1.0 ROC-AUC. Tracing this
back showed the SQL view itself was built using the same kind of hard
threshold on
hemoglobin,sbp_min, andheart_rate_max— so the "real" label was still deterministic given those columns.
Fix: hemoglobin, sbp_min, and heart_rate_max — along with
shock_index_mean and pulse_pressure_mean, which are derived from those
same columns and therefore still leak the same information indirectly — are
excluded from the feature set X. The model is trained only on genuinely
independent signals: Gender, heart_rate_mean, resprate_mean, o2sat_mean,
sbp_mean, dbp_mean.
df['needs_transfusion'].value_counts()
df["shock_index_mean"] = df["heart_rate_mean"]/df["sbp_min"]
df["pulse_pressure_mean"] = df["sbp_min"] - df["dbp_mean"]
df.head()
print(df["shock_index_mean"].isna().sum())
print(df["pulse_pressure_mean"].isna().sum())
Exploratory Data Analysis (EDA)¶
type_ = ["no","yes"]
fig = make_subplots(rows =1,cols =1)
night_colour = ["#BCD9A2","#A82323"]
fig.add_trace(go.Pie(labels = type_ , values = df["needs_transfusion"].value_counts(), name = "Transfusion", marker_colors = night_colour ))
fig.update_traces(hole = 0.4,hoverinfo="label+percent+name",textfont_size = 16)
fig.update_layout(title_text="<b> Blood Transfusion Pie Chart <b>",
width = 500,
height = 500,
annotations = [dict(text='Tranfusion',x=0.5,y=0.5,showarrow = False)])
fig.show()
cols = ['heart_rate_mean', 'sbp_min', 'dbp_mean', 'o2sat_mean', 'resprate_mean', 'hemoglobin', 'shock_index_mean', 'pulse_pressure_mean']
fig = make_subplots(rows=2, cols=4, subplot_titles=cols)
for i, col in enumerate(cols):
row = i // 4 + 1
col_pos = i % 4 + 1
for outcome in [0, 1]:
fig.add_trace(
go.Box(y=df[df['needs_transfusion']==outcome][col], name=str(outcome), showlegend=(i==0)),
row=row, col=col_pos
)
fig.update_layout(height=700, title_text="Vitals by Transfusion Need")
fig.show()
import plotly.express as px
corr = df.corr(numeric_only=True).round(2)
fig = px.imshow(
corr,
text_auto=True, # shows numbers on each cell
color_continuous_scale='Viridis',
aspect='auto',
labels=dict(color='Correlation')
)
fig.update_layout(
width=800,
height=800,
font=dict(size=10), # shrink font so numbers fit in cells
title='Correlation Heatmap'
)
fig.update_traces(textfont_size=10) # controls the number size inside cells specifically
fig.show()
print(df[["heart_rate_mean","resprate_mean","o2sat_mean","sbp_min","dbp_mean"]].describe())
df.groupby('Gender')['needs_transfusion'].mean()
Handling Class Imbalance¶
needs_transfusion is rare (14 positive cases out of 206 rows, ~6.8%). SMOTE
generates synthetic minority-class examples by interpolating between real
positive cases, so the model isn't incentivized to just predict the majority
class for every row.
Important: SMOTE is applied inside the cross-validation pipeline
(imblearn.pipeline.Pipeline), not once on the whole dataset beforehand.
Applying it before splitting would let synthetic points generated from
training-fold data leak into test folds — the same category of mistake as the
label leakage above, just relocated into the evaluation step.
Smote fixes the training imbalance problem beacuse without it the model can get high accuracy by just predicting "no transfusion" for everyone, since that's right 95% of the time, and SMOTE stops that lazy shortcut by making the positive class numerically comparable during training.
The model can learn who needs tranfusion only set by those 14 real examples only, it dosn't matter how many synethetic points we get generated from them.
columns_to_drop = [
'needs_transfusion', 'stay_id', 'subject_id',
'hemoglobin', 'sbp_min', 'heart_rate_max', # directly in the label formula
'shock_index_mean', 'pulse_pressure_mean' # derived from sbp_min, so still leaky
]
X = df.drop(columns=columns_to_drop, errors='ignore')
Y = df['needs_transfusion']
print(X.columns.tolist())
Cross-Validation¶
A single 80/20 train/test split proved unreliable here — with only 14 positive cases total, one particular split can land 2-3 (or even 0-1) positive cases in the test set, making any single score highly sensitive to which rows happened to land where.
StratifiedKFold(n_splits=5) is used instead: the data is divided into 5 folds,
each preserving the same ~6.8% positive rate, and every row is used as test
data exactly once across the 5 rounds. Reporting the mean and standard
deviation across folds (rather than a single number) makes the uncertainty in
this estimate explicit, which matters given the small sample size.
cv = StratifiedKFold(n_splits = 5, shuffle = True, random_state = 42)
rf_pipeline = ImbPipeline([('smote',SMOTE(random_state = 42, k_neighbors = 3)),
('model',RandomForestClassifier(n_estimators=300, max_depth=5, random_state=42,n_jobs=-1))
])
rf_scores = cross_val_score(rf_pipeline,X,Y, cv = cv, scoring= 'roc_auc')
print("Per-fold ROC-AUC", np.round(rf_scores,3))
print(f"Mean:{rf_scores.mean():.3f} ± Std: {rf_scores.std():3f}")
oof_proba = cross_val_predict(rf_pipeline,X,Y,cv=cv,method='predict_proba')[:,1]
results_oof = pd.DataFrame({'actual': Y.values, 'predicted_proba':oof_proba})
print(results_oof.sort_values('predicted_proba', ascending=False).head(20))
''' This is a demo data the majority of classes where classified as not in a need of tranfusion so the accuracy of the model will be very high so we wont be using that as a feild for classification report'''
Model Comparison: Random Forest vs. Logistic Regression¶
Logistic Regression is included as an interpretable baseline. With this little
data, a simpler model with fewer parameters can generalize better than a more
flexible one — Random Forest has more room to fit noise in the training folds,
while Logistic Regression's simplicity acts as a form of regularization almost
by construction. Both models are evaluated on the same 5 folds (same cv
object passed to both), so any difference in score reflects the models
themselves, not random variation in how the data was split.
# cv already defined
logreg_pipeline = ImbPipeline([('smote',SMOTE(random_state=42, k_neighbors = 3)),
('model',LogisticRegression(max_iter=1000,random_state=42))
])
logreg_scores = cross_val_score(logreg_pipeline, X, Y, cv=cv, scoring='roc_auc')
print("Per-fold ROC-AUC:", np.round(logreg_scores, 3))
print(f"Mean: {logreg_scores.mean():.3f} ± Std: {logreg_scores.std():.3f}")
Model comparison (cross-validated, 5 folds):
| Model | Mean ROC-AUC | Std | Notes |
|---|---|---|---|
| Random Forest (+ SMOTE) | 0.791 |
0.063 |
Non-monotonic on heart rate (see findings) |
| Logistic Regression (+ SMOTE) | 0.812 |
0.178 |
Interpretable baseline; matched or outperformed RF |
Threshold Selection and Calibration¶
Out-of-fold predicted probabilities (cross_val_predict) are used here rather
than a single test set's predictions — this uses every row's prediction from
when it was held out across the 5 folds, giving a fuller and less noisy picture
than the 42-row single split used earlier.
The default 0.5 classification threshold isn't assumed to be correct. Instead, a precision-recall curve is used to choose a threshold deliberately: for a transfusion-need screening use case, missing a real positive case (false negative) is arguably costlier than an unnecessary follow-up check (false positive), so a threshold favoring recall is a defensible, explicit choice — not the default.
A calibration curve checks whether predicted probabilities are trustworthy in their own right (e.g. does "70% predicted risk" really correspond to about 70% of such cases being positive). With only 14 positive cases, expect this curve to look noisy rather than a smooth diagonal — that's an honest reflection of the sample size, not a modeling error.
precision , recall, thresholds = precision_recall_curve(Y,oof_proba)
plt.figure(figsize=(7,5))
plt.plot(thresholds, precision[:-1], label = 'Precision')
plt.plot(thresholds, recall[:-1], label = 'recall')
plt.xlabel('Threshold')
plt.ylabel('Score')
plt.title('Precision and Recall vs Tresholds')
plt.legend()
plt.show()
target_recall = 0.8
idx = np.where(recall >= target_recall)[0][-1]
print(f"Thresholds for >= {target_recall} recall: {thresholds[idx]:.3f}")
print(f"At that threshold — precision: {precision[idx]:.3f}, recall: {recall[idx]:.3f}")
rac_pos, mean_pred = calibration_curve(Y, oof_proba, n_bins=5, strategy='quantile')
plt.figure(figsize=(6,6))
plt.plot(mean_pred, rac_pos, marker='o', label='Random Forest')
plt.plot([0,1],[0,1], linestyle='--', color='gray', label='Perfectly calibrated')
plt.xlabel('Mean predicted probability')
plt.ylabel('Fraction of actual positives')
plt.title('Calibration Curve')
plt.legend()
plt.show()
Final Model and Deployment¶
Cross-validation is used for evaluating the model honestly — not for
producing the final deployed model. The final model is refit on the entire
available dataset (all 206 rows, SMOTE-balanced), since more training data is
better once evaluation is complete. This final model is exported to JavaScript
via m2cgen, allowing it to run entirely client-side in the accompanying
interactive demo (transfusion_dashboard.html) with no server or API required.
final_model = RandomForestClassifier(n_estimators=300, max_depth=5, random_state=42)
sm_final = SMOTE(random_state=42, k_neighbors=3)
X_res, Y_res = sm_final.fit_resample(X, Y)
final_model.fit(X_res, Y_res)
print("Features model expects:", final_model.feature_names_in_)
code = m2c.export_to_javascript(final_model)
with open("model.js", "w") as f:
f.write(code)
Limitations¶
- Sample size. 206 patient stays and 14 positive cases is small; all metrics here should be read as indicative, not definitive. Cross-validation standard deviations reflect this directly.
- Synthetic label.
needs_transfusionis rule-derived for this demo, not a real clinical outcome — this project demonstrates a leakage-aware, honestly evaluated ML methodology, not a validated clinical tool. - No additional real data was available for this iteration. A larger patient sample would likely improve reliability more than further modeling changes.
- Age is collected in the interactive demo for context but is not one of the model's trained features.
End.¶