Forecasting glucose from CGM data

import numpy as np

def synthetic_cgm(days=14, seed=0):
    """A crude 5-minute CGM trace: baseline, meal rises, slow decay, overnight dips and sensor noise."""
    rng = np.random.default_rng(seed)
    steps = days * 288                                   # 288 readings per day
    glucose = np.empty(steps)
    g, meal_effect = 130.0, 0.0
    for t in range(steps):
        minute = (t % 288) * 5
        if minute in (8 * 60, 13 * 60, 20 * 60) and rng.random() < 0.95:
            meal_effect += rng.uniform(60, 160)          # a meal of varying size
        if minute == 3 * 60 and rng.random() < 0.4:
            meal_effect -= rng.uniform(120, 170)         # an occasional overnight low
        meal_effect *= 0.965                             # absorbed and handled over a few hours
        g += 0.03 * (125 - g) + 0.04 * meal_effect + rng.normal(0, 1.5)
        glucose[t] = g
    return np.clip(glucose + rng.normal(0, 4, steps), 40, 400)   # sensor noise and range limits

cgm = synthetic_cgm()
in_range = np.mean((cgm >= 70) & (cgm <= 180))
print(f"{len(cgm)} readings over 14 days; mean {cgm.mean():.0f} mg/dL; "
      f"time in range {in_range:.0%}; below 70: {np.mean(cgm < 70):.1%}; above 180: {np.mean(cgm > 180):.0%}")

# 4032 readings over 14 days; mean 152 mg/dL; time in range 79%; below 70: 2.9%; above 180: 18%

from sklearn.linear_model import Ridge

HISTORY, HORIZON = 6, 6                       # use the last 30 minutes to predict 30 minutes ahead
X = np.array([cgm[t - HISTORY:t] for t in range(HISTORY, len(cgm) - HORIZON)])
y = np.array([cgm[t + HORIZON - 1] for t in range(HISTORY, len(cgm) - HORIZON)])
split = 10 * 288 - HISTORY                    # train on days 1-10, test on days 11-14 (temporal split)
X_tr, X_te, y_tr, y_te = X[:split], X[split:], y[:split], y[split:]

forecasts = {
    "Persistence (last value)": X_te[:, -1],
    "Linear trend (last 15 min)": X_te[:, -1] + (X_te[:, -1] - X_te[:, -4]) / 3 * HORIZON,
    "Ridge regression on last 30 min": Ridge(alpha=1.0).fit(X_tr, y_tr).predict(X_te),
}
for name, pred in forecasts.items():
    rmse = np.sqrt(np.mean((pred - y_te) ** 2))
    mard = np.mean(np.abs(pred - y_te) / y_te)
    print(f"{name:<32} RMSE {rmse:5.1f} mg/dL   MARD {mard:5.1%}")

# Persistence (last value)         RMSE  10.2 mg/dL   MARD  5.6%
# Linear trend (last 15 min)       RMSE  16.7 mg/dL   MARD  9.3%
# Ridge regression on last 30 min  RMSE   9.3 mg/dL   MARD  5.1%

low = y_te < 70                                          # a low 30 minutes later
print(f"{low.sum()} test time points are followed by a low (< 70 mg/dL) 30 minutes later")
for name, pred in forecasts.items():
    alert = pred < 80                                    # alert if the forecast is below 80 mg/dL
    sens = np.mean(alert[low])
    false_alerts_per_day = np.sum(alert & ~low) / 4      # four test days
    print(f"{name:<32} lows caught {sens:4.0%}   false alerts per day {false_alerts_per_day:4.1f}")

# 56 test time points are followed by a low (< 70 mg/dL) 30 minutes later
# Persistence (last value)         lows caught  93%   false alerts per day  9.0
# Linear trend (last 15 min)       lows caught  82%   false alerts per day  9.8
# Ridge regression on last 30 min  lows caught  89%   false alerts per day  6.5

new_low = low & (X_te[:, -1] >= 70)                      # not low yet, but will be in 30 minutes
print(f"{new_low.sum()} of these lows start after a normal current reading")
for name, pred in forecasts.items():
    print(f"{name:<32} new lows predicted in advance {np.mean((pred < 80)[new_low]):4.0%}")

# 16 of these lows start after a normal current reading
# Persistence (last value)         new lows predicted in advance  75%
# Linear trend (last 15 min)       new lows predicted in advance  69%
# Ridge regression on last 30 min  new lows predicted in advance  62%