Files
glucose-dashboard/reporting.py
T
willemjan 374ead70e5 Replace monthly GMI bar chart with daily scatter + rolling avg + linear trend/CI
Adds daily_gmi() and gmi_trend(): per-day GMI as scatter points, a 7-day
rolling average line, and a linear regression trend line with a 95%
confidence band (normal approximation, no scipy dependency, consistent
with the existing preference for a lean container image).

Recomputed from the already-fetched period data on each page load; the
existing 1-hour Streamlit cache on load_range_data() already avoids
repeated Nightscout calls for the same period, so no separate storage
layer is needed at the current max. 90-day report window.

Closes gitea #19.
2026-08-01 22:57:54 +02:00

517 lines
20 KiB
Python

"""
Rapportage-logica voor de klinische rapportagepagina's.
Zelfde opzet als data_processing.py: pure functies, geen Streamlit-imports, dus
los te unit-testen. Zie docs/requirements-reporting.md voor de klinische
definities waarop dit gebaseerd is (Battelino et al. 2019, Bergenstal et al.
2018).
"""
from __future__ import annotations
from itertools import pairwise
import numpy as np
import pandas as pd
from nightscout_client import INSULIN_EVENT_TYPES
# Klinische drempels in mmol/L (Battelino et al. 2019, Diabetes Care 42(8):1593-1603).
VERY_LOW = 3.0
LOW = 3.9
HIGH = 10.0
VERY_HIGH = 13.9
MINUTES_PER_DAY = 24 * 60
def gmi_percent(mean_mmol: float) -> float:
"""Geschat HbA1c (Glucose Management Indicator) uit gemiddelde glucose.
GMI% = 3.31 + 0.4311 * gemiddelde(mmol/L) — Bergenstal et al. 2018,
Diabetes Care 41(11):2275-2280. Dit is GEEN vervanging voor een labwaarde
(kan tot ~0.5 procentpunt afwijken door o.a. rode-bloedcel-factoren) en moet
als zodanig gelabeld worden in de UI.
"""
return round(3.31 + 0.4311 * mean_mmol, 1)
GMI_TREND_COLUMNS = ["date", "gmi", "rolling_avg", "trend", "trend_lower", "trend_upper"]
def daily_gmi(glucose: pd.DataFrame) -> pd.DataFrame:
"""Geschat HbA1c (GMI) per kalenderdag — scatterpunten voor de trendgrafiek (gitea #19)."""
if glucose.empty:
return pd.DataFrame(columns=["date", "gmi"])
df = glucose.copy()
df["date"] = df["timestamp"].dt.date
daily_mean = df.groupby("date")["mmol"].mean()
result = daily_mean.apply(gmi_percent).reset_index()
result.columns = ["date", "gmi"]
return result.sort_values("date").reset_index(drop=True)
def gmi_trend(glucose: pd.DataFrame, rolling_window_days: int = 7) -> pd.DataFrame:
"""Dagelijkse GMI-punten plus voortschrijdend gemiddelde en een lineaire
trendlijn met 95%-betrouwbaarheidsband (gitea #19).
Gebruikt een normaalbenadering (z=1.96) voor de betrouwbaarheidsband i.p.v.
een exacte t-verdeling — dat zou een scipy-dependency toevoegen aan een
verder bewust lichte container (dezelfde afweging als het niet gebruiken
van kaleido in pdf_export.py). Bij de hier gebruikelijke steekproefgroottes
(14-90 dagpunten) is het verschil met de exacte t-kritieke waarde
verwaarloosbaar voor een visuele trendindicatie.
Elke keer opnieuw berekend uit de al-opgehaalde periode-data (geen aparte
database) — de bestaande 1-uurs Streamlit-cache op load_range_data()
voorkomt al herhaalde Nightscout-calls voor dezelfde periode. Bij een
veel langere trend (jaren i.p.v. de huidige max. 90 dagen) zou een aparte
opslag wél de moeite waard worden.
"""
daily = daily_gmi(glucose)
if len(daily) < 2:
return daily.reindex(columns=GMI_TREND_COLUMNS)
daily = daily.copy()
daily["rolling_avg"] = daily["gmi"].rolling(window=rolling_window_days, min_periods=1).mean()
x = np.arange(len(daily), dtype=float)
y = daily["gmi"].to_numpy()
slope, intercept = np.polyfit(x, y, 1)
fitted = intercept + slope * x
residuals = y - fitted
degrees_of_freedom = max(len(x) - 2, 1)
residual_std = np.sqrt(np.sum(residuals**2) / degrees_of_freedom)
x_mean = x.mean()
sum_sq_x = np.sum((x - x_mean) ** 2)
standard_error = (
np.zeros_like(x)
if sum_sq_x == 0
else residual_std * np.sqrt(1 / len(x) + (x - x_mean) ** 2 / sum_sq_x)
)
z_95 = 1.96
daily["trend"] = fitted
daily["trend_lower"] = fitted - z_95 * standard_error
daily["trend_upper"] = fitted + z_95 * standard_error
return daily[GMI_TREND_COLUMNS]
def percent_cv(glucose: pd.DataFrame) -> float:
"""Glucosevariabiliteit als %CV (standaarddeviatie / gemiddelde * 100)."""
if glucose.empty:
return 0.0
mean = glucose["mmol"].mean()
if mean == 0:
return 0.0
return round(glucose["mmol"].std() / mean * 100, 1)
def tir_breakdown_5band(glucose: pd.DataFrame) -> dict[str, float]:
"""Percentage metingen in elk van de 5 klinische glucosebanden.
Banden: very_low (<3.0), low (3.0-3.9), in_range (3.9-10.0),
high (10.0-13.9), very_high (>13.9). Sluit exact aan op de 5-bands AGP-
indeling uit de consensus.
"""
if glucose.empty:
return {"very_low": 0.0, "low": 0.0, "in_range": 0.0, "high": 0.0, "very_high": 0.0}
def _band(mmol: float) -> str:
if mmol < VERY_LOW:
return "very_low"
if mmol < LOW:
return "low"
if mmol <= HIGH:
return "in_range"
if mmol <= VERY_HIGH:
return "high"
return "very_high"
counts = glucose["mmol"].apply(_band).value_counts(normalize=True) * 100
return {
band: round(counts.get(band, 0.0), 1)
for band in ("very_low", "low", "in_range", "high", "very_high")
}
def minutes_per_day(percentage: float) -> float:
"""Reken een percentage van de tijd om naar minuten per dag (ruwe vuistregel)."""
return round(percentage / 100 * MINUTES_PER_DAY, 0)
def sensor_active_percent(glucose: pd.DataFrame, start, end, expected_interval_minutes: int = 5) -> float:
"""Percentage van de periode waarin daadwerkelijk metingen binnenkwamen.
Vuistregel op basis van het verwachte meetinterval (CGM's meten doorgaans
elke 5 minuten) — geen exacte sensor-uptime-registratie, maar voldoende om
een periode met te veel hiaten te signaleren (zie FR9 / NFR "minimaal
geldig rapport").
"""
total_minutes = (pd.Timestamp(end) - pd.Timestamp(start)).total_seconds() / 60
if total_minutes <= 0:
return 0.0
expected_readings = total_minutes / expected_interval_minutes
if expected_readings <= 0:
return 0.0
return round(min(len(glucose) / expected_readings, 1.0) * 100, 1)
def completeness_warning(glucose: pd.DataFrame, start, end) -> str | None:
"""Bouw een waarschuwingstekst als er te weinig data is voor de gekozen periode.
Pure functie zodat zowel de interactieve pagina's (ui_common.py, als
st.warning) als het PDF-rapport (pdf_export.py, als tekstregel) exact
dezelfde beoordeling en bewoording gebruiken — een rapport dat naar het
diabetesteam gaat moet dezelfde disclaimer bevatten als wat de gebruiker
zelf op het scherm ziet.
Onderscheidt twee oorzaken:
- een structureel gat aan het BEGIN van de periode (de vroegste
binnengekomen meting ligt later dan de gekozen startdatum) — data is
simpelweg niet opgehaald/beschikbaar voor dat stuk, geen sensorhiaat.
- een lagere dekking dan de klinische AGP-ondergrens (Danne et al. 2017)
zonder zo'n gat aan het begin — verspreide hiaten.
Returns:
None als er geen waarschuwing nodig is; anders de te tonen tekst.
"""
if glucose.empty:
return "Geen glucosedata gevonden in deze periode."
start_date = pd.Timestamp(start).date()
end_date = pd.Timestamp(end).date()
active_pct = sensor_active_percent(glucose, start, end)
earliest = glucose["timestamp"].min().date()
missing_days_at_start = max((earliest - start_date).days, 0)
if missing_days_at_start > 0:
return (
f"Nightscout heeft pas data vanaf {earliest:%d-%m-%Y} teruggegeven — "
f"de eerste {missing_days_at_start} dag(en) van de gekozen periode "
f"({start_date:%d-%m-%Y} t/m {end_date:%d-%m-%Y}) ontbreken volledig, niet alleen "
f"losse sensorhiaten. Sensordekking over de wél opgehaalde periode: "
f"{active_pct:.0f}%. Onderstaande cijfers zijn dus berekend over een "
"kortere periode dan geselecteerd."
)
if active_pct < 70:
# ">=" i.p.v. "≥": deze tekst komt ook in het PDF-rapport terecht
# (pdf_export.py), en fpdf2's ingebouwde Helvetica-font ondersteunt
# alleen Latin-1 — zie gitea issue #11 voor het opruimen hiervan
# (Unicode-lettertype embedden) zodra dat de moeite waard is.
return (
f"Slechts {active_pct:.0f}% sensordekking in deze periode (streefwaarde "
">=70% voor een klinisch betrouwbaar rapport, Danne et al. 2017). De "
"onderstaande cijfers zijn mogelijk minder betrouwbaar door de vele hiaten."
)
return None
def _state_blocks(df: pd.DataFrame, out: pd.Series) -> list[dict]:
"""Groepeer aaneengesloten metingen met dezelfde in/out-status.
Eindtijd van elk blok = starttijd van het volgende blok (de status geldt
tot de eerstvolgende meting anders aangeeft); het laatste blok behoudt
zijn eigen laatste meting als eind.
"""
block_id = (out != out.shift()).cumsum()
blocks = [
{
"is_out": bool(out.loc[group.index[0]]),
"start": group["timestamp"].iloc[0],
"end": group["timestamp"].iloc[-1],
}
for _, group in df.groupby(block_id)
]
for i in range(len(blocks) - 1):
blocks[i]["end"] = blocks[i + 1]["start"]
return blocks
def _maybe_add_episode(events: list[dict], start, end, min_duration_minutes: int) -> None:
duration = (end - start).total_seconds() / 60
if duration >= min_duration_minutes:
events.append({"start": start, "end": end, "duration_minutes": round(duration, 1)})
def _scan_episodes(blocks: list[dict], min_duration_minutes: int) -> list[dict]:
"""Statelijke doorloop van de blokken volgens de persistentieregel.
Een lopende episode wordt alleen verlengd voorbij een kort in-bereik-blok
als er daarna weer een out-blok komt (het korte herstel wordt dan met
terugwerkende kracht "geabsorbeerd"). Een kort slot-blok waarna de data
simpelweg ophoudt telt dus NIET automatisch mee als voortzetting — anders
zou een niet-bevestigd herstel de duur van het event kunstmatig oprekken.
"""
events: list[dict] = []
episode_start = None
episode_end = None
for block in blocks:
if block["is_out"]:
episode_start = episode_start or block["start"]
episode_end = block["end"]
continue
if episode_start is None:
continue
recovery_minutes = (block["end"] - block["start"]).total_seconds() / 60
if recovery_minutes >= min_duration_minutes:
_maybe_add_episode(events, episode_start, episode_end, min_duration_minutes)
episode_start = None
episode_end = None
# anders: kortstondig herstel — episode blijft open, wordt pas
# verlengd als er weer een out-blok volgt (hierboven).
if episode_start is not None:
_maybe_add_episode(events, episode_start, episode_end, min_duration_minutes)
return events
def detect_events(
glucose: pd.DataFrame,
threshold: float,
direction: str,
min_duration_minutes: int = 15,
) -> pd.DataFrame:
"""Detecteer discrete hypo-/hyperglykemie-events volgens de persistentieregel.
Een event begint zodra glucose ononderbroken voorbij de drempel is
geweest voor >= min_duration_minutes, en eindigt pas zodra glucose
ononderbroken hersteld is voor >= min_duration_minutes — een korte terugval
binnen die hersteltijd beëindigt het event dus niet (Battelino et al.
2019). `direction` is "below" (hypo) of "above" (hyper).
Geeft een DataFrame terug met kolommen 'start', 'end', 'duration_minutes'.
"""
if glucose.empty:
return pd.DataFrame(columns=["start", "end", "duration_minutes"])
df = glucose.sort_values("timestamp").reset_index(drop=True)
if direction == "below":
out = df["mmol"] < threshold
elif direction == "above":
out = df["mmol"] > threshold
else:
raise ValueError("direction moet 'below' of 'above' zijn")
blocks = _state_blocks(df, out)
events = _scan_episodes(blocks, min_duration_minutes)
return pd.DataFrame(events, columns=["start", "end", "duration_minutes"])
def summarize_events(events: pd.DataFrame, period_days: float) -> dict[str, float]:
"""Samenvattende statistieken over een set gedetecteerde events."""
if events.empty or period_days <= 0:
return {
"count": 0,
"avg_duration_minutes": 0.0,
"longest_duration_minutes": 0.0,
"total_minutes": 0.0,
"events_per_week": 0.0,
}
total_minutes = events["duration_minutes"].sum()
return {
"count": len(events),
"avg_duration_minutes": round(events["duration_minutes"].mean(), 1),
"longest_duration_minutes": round(events["duration_minutes"].max(), 1),
"total_minutes": round(total_minutes, 1),
"events_per_week": round(len(events) / period_days * 7, 1),
}
AGP_PERCENTILE_COLUMNS = ["minute_of_day", "p5", "p25", "p50", "p75", "p95"]
def agp_percentile_bands(glucose: pd.DataFrame, bin_minutes: int = 30) -> pd.DataFrame:
"""Aggregeer glucosewaarden per tijdstip-van-de-dag naar percentielbanden.
Groepeert alle metingen in de periode naar hun tijd-van-de-dag (in bins
van `bin_minutes`), en berekent per bin de 5e/25e/50e/75e/95e percentiel
— de standaard AGP-percentielgrafiek over een samengesteld 24-uursbeeld
(Danne et al. 2017, Diabetes Care 40(12):1631-1640).
"""
if glucose.empty:
return pd.DataFrame(columns=AGP_PERCENTILE_COLUMNS)
df = glucose.copy()
minute_of_day = df["timestamp"].dt.hour * 60 + df["timestamp"].dt.minute
df["bin"] = (minute_of_day // bin_minutes) * bin_minutes
grouped = df.groupby("bin")["mmol"].quantile([0.05, 0.25, 0.5, 0.75, 0.95]).unstack()
grouped.columns = ["p5", "p25", "p50", "p75", "p95"]
grouped = grouped.reset_index().rename(columns={"bin": "minute_of_day"})
return grouped[AGP_PERCENTILE_COLUMNS].sort_values("minute_of_day").reset_index(drop=True)
def overnight_daytime_split(
glucose: pd.DataFrame, overnight_start: int = 0, overnight_end: int = 6
) -> dict[str, dict]:
"""Tijd-in-bereik apart voor de nacht (overnight_start-overnight_end uur) en overdag."""
empty = tir_breakdown_5band(glucose)
if glucose.empty:
return {"overnight": empty, "daytime": empty}
hour = glucose["timestamp"].dt.hour
is_overnight = (hour >= overnight_start) & (hour < overnight_end)
return {
"overnight": tir_breakdown_5band(glucose[is_overnight]),
"daytime": tir_breakdown_5band(glucose[~is_overnight]),
}
WEEKDAY_NAMES = ["Maandag", "Dinsdag", "Woensdag", "Donderdag", "Vrijdag", "Zaterdag", "Zondag"]
def day_of_week_pattern(glucose: pd.DataFrame) -> pd.DataFrame:
"""Gemiddelde glucose en tijd-in-bereik per weekdag."""
columns = ["weekday", "mean_mmol", "in_range_pct", "count"]
if glucose.empty:
return pd.DataFrame(columns=columns)
df = glucose.copy()
df["weekday"] = df["timestamp"].dt.dayofweek
rows = []
for weekday in range(7):
day_df = df[df["weekday"] == weekday]
if day_df.empty:
rows.append(
{"weekday": WEEKDAY_NAMES[weekday], "mean_mmol": None, "in_range_pct": None, "count": 0}
)
continue
tir = tir_breakdown_5band(day_df)
rows.append(
{
"weekday": WEEKDAY_NAMES[weekday],
"mean_mmol": round(day_df["mmol"].mean(), 1),
"in_range_pct": tir["in_range"],
"count": len(day_df),
}
)
return pd.DataFrame(rows, columns=columns)
SMB_TYPE = "SMB"
BASAL_EVENT_TYPE = "Temp Basal"
def bolus_stats(treatments: pd.DataFrame) -> dict:
"""Samenvatting van bolus-insuline, opgesplitst naar type.
AndroidAPS in SMB-looping registreert een automatische micro-bolus niet
als eigen `eventType`, maar als `eventType == "Correction Bolus"` met
`type == "SMB"` — zonder dit onderscheid werden alle SMB's stilzwijgend
als "Correction Bolus" (handmatige correctie) meegeteld, wat gitea #20
meldde als "bolus stats empty/onduidelijk". Geverifieerd rechtstreeks
tegen de productiedata: bij deze gebruiker is vrijwel elke "Correction
Bolus" in werkelijkheid een automatische SMB.
Let op: dit betreft alleen bolusinsuline (voeding/correctie/combo/SMB).
Basaalinsuline zit hier niet in — zie basal_stats().
"""
if treatments.empty:
return {"total_units": 0.0, "count": 0, "by_type": {}}
bolus = treatments[treatments["eventType"].isin(INSULIN_EVENT_TYPES) & treatments["insulin"].notna()].copy()
if bolus.empty:
return {"total_units": 0.0, "count": 0, "by_type": {}}
is_smb = (bolus["eventType"] == "Correction Bolus") & (bolus.get("type", "") == SMB_TYPE)
bolus["category"] = bolus["eventType"].where(~is_smb, SMB_TYPE)
by_type = bolus.groupby("category")["insulin"].sum().round(1).to_dict()
return {
"total_units": round(bolus["insulin"].sum(), 1),
"count": len(bolus),
"by_type": by_type,
}
def basal_stats(treatments: pd.DataFrame, period_days: float) -> dict:
"""Totale basaalinsuline uit Temp Basal-treatments (rate x duur), zonder pompprofiel.
AndroidAPS zet in gesloten-lus-modus vrijwel continu Temp Basal-treatments
met een `absolute` snelheid (E/uur) en `duration` (minuten); de som van
afgeleverde eenheden daaruit benadert de daadwerkelijke basaaldosis
nauwkeurig genoeg zonder het pompprofiel te hoeven ophalen (gitea #20).
Periodes zonder actieve temp basal (dan geldt het profiel-basaaltempo)
worden hierin gemist — dit is dus een ondergrens, geen exacte totale-
dagdosis.
"""
if treatments.empty or period_days <= 0:
return {"total_units": 0.0, "avg_units_per_day": 0.0}
temp_basal = treatments[
(treatments["eventType"] == BASAL_EVENT_TYPE)
& treatments["absolute"].notna()
& treatments["duration"].notna()
]
if temp_basal.empty:
return {"total_units": 0.0, "avg_units_per_day": 0.0}
total_units = (temp_basal["absolute"] * temp_basal["duration"] / 60).sum()
return {
"total_units": round(total_units, 1),
"avg_units_per_day": round(total_units / period_days, 1),
}
SITE_INTERVAL_COLUMNS = [
"start",
"end",
"duration_days",
"mean_mmol",
"in_range_pct",
"total_bolus_units",
"avg_bolus_units_per_day",
"total_carbs_g",
]
def _site_interval_row(glucose: pd.DataFrame, treatments: pd.DataFrame, start, end) -> dict:
duration_days = max((end - start).total_seconds() / 86400, 1e-9)
window_glucose = glucose[(glucose["timestamp"] >= start) & (glucose["timestamp"] < end)]
window_treatments = (
treatments[(treatments["timestamp"] >= start) & (treatments["timestamp"] < end)]
if not treatments.empty
else treatments
)
bolus = bolus_stats(window_treatments)
total_carbs = window_treatments["carbs"].fillna(0).sum() if not treatments.empty else 0.0
tir = tir_breakdown_5band(window_glucose)
return {
"start": start,
"end": end,
"duration_days": round(duration_days, 1),
"mean_mmol": round(window_glucose["mmol"].mean(), 1) if not window_glucose.empty else None,
"in_range_pct": tir["in_range"],
"total_bolus_units": bolus["total_units"],
"avg_bolus_units_per_day": round(bolus["total_units"] / duration_days, 1),
"total_carbs_g": round(total_carbs, 0),
}
def site_interval_stats(
glucose: pd.DataFrame, treatments: pd.DataFrame, site_changes: pd.DataFrame
) -> pd.DataFrame:
"""Statistieken per interval tussen opeenvolgende pompwissels (FR8).
Voor elk interval sinds een pompwissel: duur in dagen ("leeftijd" van de
pompsite bij het einde van het interval), gemiddelde glucose,
tijd-in-bereik, en bolus-insuline/koolhydraten in dat interval — als
benadering van of insulinebehoefte samenhangt met de leeftijd van de
pompsite. Zie bolus_stats() voor de kanttekening dat basaalinsuline hier
niet in zit.
"""
if site_changes.empty or glucose.empty:
return pd.DataFrame(columns=SITE_INTERVAL_COLUMNS)
change_times = sorted(site_changes["timestamp"])
boundaries = [*change_times, glucose["timestamp"].max()]
rows = [_site_interval_row(glucose, treatments, start, end) for start, end in pairwise(boundaries)]
return pd.DataFrame(rows, columns=SITE_INTERVAL_COLUMNS)