📊 Tutorial: CGM Data Analysis with JupyterHealth

from enum import Enum

import cgmquantify
import pandas as pd
from jupyterhealth_client import Code, JupyterHealthClient

CGM = Code.BLOOD_GLUCOSE

pd.options.mode.chained_assignment = None

jh_client = JupyterHealthClient()

# construct a dict of all organizations
org_dict = {org['id']: org for org in jh_client.list_organizations()}
# populate a tree of children, from each organization's `partOf` reference to its parent
# the root org has id=0

for org in org_dict.values():
    parent_id = org['partOf']
    if parent_id is not None:
        parent = org_dict[parent_id]
        parent.setdefault("child_ids", []).append(org['id'])
    

def print_org(org, indent=''):
    print(f"{indent}[{org['id']}] {org['name']}")
    for org_id in org.get('child_ids', []):
        print_org(org_dict[org_id], indent=" " * len(indent) + " ⮑")


for org_id in org_dict[0]['child_ids']:
    print_org(org_dict[org_id])

for study in jh_client.list_studies():
    print(f"  - [{study['id']}] {study['name']} org:{study['organization']['name']}")

study_id = 30011
study = jh_client.get_study(study_id)
study

from dateutil.parser import parse as parse_date

# show all the patients with study data I have access to:
print(f"Patients in study [{study['id']}] {study['name']}")

for patient in jh_client.list_patients():
    consents = jh_client.get_patient_consents(patient['id'])
    consented = {
        _study['id']: _study
        for _study in consents['studies']
    }
    pending_consent = {
        _study['id']: _study
        for _study in consents['studiesPendingConsent']
    }
    if study_id not in consented and study_id not in pending_consent:
        # ignore patients not in the study
        continue
    print(f"[{patient['id']}] {patient['nameFamily']}, {patient['nameGiven']}: {patient['telecomEmail']}")
    if study_id in consented:
        study_consent = consented[study_id]
        for scope_consent in study_consent['scopeConsents']:
            date = parse_date(scope_consent['consentedTime']).date()
            print(f" - {scope_consent['code']['text']}, consented {date}")
    if study_id in pending_consent:
        study_pending_consent = pending_consent[study_id]
        for scope_consent in study_pending_consent['pendingScopeConsents']:
            print(f" - {scope_consent['code']['text']} (not consented)")

jh_client.get_patient_consents(40001)

jh_client.get_patient_consents(40054)

# pick patient id from above
patient_id = 40072

df = jh_client.list_observations_df(patient_id=patient_id, study_id=study_id, limit=10_000, code=CGM)

df

df.iloc[0]

# Reduce data to relevant subset for cgm
# make sure we're getting the units we expect
assert (df.blood_glucose_unit == 'MGDL').all()
# reduce data to interesting columns
cgm = df.loc[
    :,
    [
        "blood_glucose_value",
        "effective_time_frame_date_time_local",
    ],
]
# ensure sorted by date
cgm = cgm.sort_values("effective_time_frame_date_time_local")

# Mean of glucose
mean_glucose = cgm["blood_glucose_value"].mean()

# Standard Deviation (SD) of glucose
std_glucose = cgm["blood_glucose_value"].std()

# Coefficient of Variation (CV) = (SD / Mean) * 100
cv_glucose = (std_glucose / mean_glucose) * 100

# Print results
print(f"Mean Glucose: {mean_glucose:.2f} mg/dL")
print(f"Standard Deviation (SD): {std_glucose:.2f} mg/dL")
print(f"Coefficient of Variation (CV): {cv_glucose:.2f}%")

# Compute Glucose Management Indicator (GMI)
gmi = 3.31 + (0.02392 * mean_glucose)

# Print result
print(f"Glucose Management Indicator (GMI): {gmi:.2f}%")

# --- GRI COMPONENTS ---

# Total CGM readings
total_readings = len(cgm)
#print(f"Total CGM Readings: {total_readings}")

# Define time-in-range categories
VLow = (cgm["blood_glucose_value"] < 54).sum() / total_readings * 100
Low = ((cgm["blood_glucose_value"] >= 54) & (cgm["blood_glucose_value"] < 70)).sum() / total_readings * 100
High = ((cgm["blood_glucose_value"] > 180) & (cgm["blood_glucose_value"] <= 250)).sum() / total_readings * 100
VHigh = (cgm["blood_glucose_value"] > 250).sum() / total_readings * 100

# Calculate GRI components
hypo_component = VLow + (0.8 * Low)
hyper_component = VHigh + (0.5 * High)
gri = (3.0 * hypo_component) + (1.6 * hyper_component)
gri = min(gri, 100)  # Cap GRI at 100

# Print results
print("\n--- GRI COMPONENTS ---")
print(f"Level 2 Hypoglycemia (<54 mg/dL): {VLow:.2f}%")
print(f"Level 1 Hypoglycemia (54-69 mg/dL): {Low:.2f}%")
print(f"Level 1 Hyperglycemia (181-250 mg/dL): {High:.2f}%")
print(f"Level 2 Hyperglycemia (>250 mg/dL): {VHigh:.2f}%")
print(f"Hypoglycemia Component: {hypo_component:.2f}%")
print(f"Hyperglycemia Component: {hyper_component:.2f}%")
print(f"Glycemia Risk Index (GRI): {gri:.2f}")

# GRI GRID

import matplotlib.pyplot as plt
import numpy as np
from matplotlib.colors import ListedColormap
from matplotlib.patches import Patch

# Define the plotting area
fig, ax = plt.subplots(figsize=(8, 6))

# Create a meshgrid for the background
x = np.linspace(0, 40, 400)
y = np.linspace(0, 60, 400)
X, Y = np.meshgrid(x, y)

# Define risk zones: diagonal bands (Zone A to E)
zone = np.zeros_like(X)
zone_bounds = [(0, 20), (20, 40), (40, 60), (60, 80), (80, 100)]

for i, (low, high) in enumerate(zone_bounds):
    mask = (3.0 * X + 1.6 * Y >= low) & (3.0 * X + 1.6 * Y < high)
    zone[mask] = i + 1

# Define zone colors (Zone A to E: green to red)
zone_colors = ['#d2fbd4', '#f0fcb2', '#fff4a3', '#ffd0a1', '#ff9d9d']
zone_labels = ['Zone A', 'Zone B', 'Zone C', 'Zone D', 'Zone E']
zone_cmap = ListedColormap(zone_colors)

# Show background zones
ax.contourf(X, Y, zone, levels=[0.5,1.5,2.5,3.5,4.5,5.5], colors=zone_colors, alpha=0.5)

# Plot GRI point (replace with your values)
hypo_component
hyper_component
gri

ax.scatter(hypo_component, hyper_component, color='black', s=120, edgecolors='white', zorder=10)
ax.text(hypo_component + 0.5, hyper_component, f"GRI = {gri:.1f}", fontsize=10)

# Labels and formatting
ax.set_xlabel("Hypoglycemia Component (%)")
ax.set_ylabel("Hyperglycemia Component (%)")
ax.set_title("GRI Grid with Color-Zoned Risk Levels")
ax.set_xlim(0, 40)
ax.set_ylim(0, 60)
ax.grid(True)

# Create a legend for the zones
legend_patches = [Patch(color=zone_colors[i], label=zone_labels[i]) for i in range(len(zone_labels))]
ax.legend(handles=legend_patches, title="GRI Zones", loc="upper right", frameon=True)

plt.show()

# CGM and AGP Visualization

# 1. Imports
import datetime

import matplotlib.pyplot as plt
import pandas as pd

# 2. Enums and Utility Mappings


class Category(Enum):
    very_low = "Very Low"
    low = "Low"
    target_range = "Target Range"
    high = "High"
    very_high = "Very High"


category_colors = {
    Category.very_low.value: "#a00",
    Category.low.value: "#f44",
    Category.target_range.value: "#CDECCD",
    Category.high.value: "#FDBE85",
    Category.very_high.value: "#FD8D3C",
}

# 3. Data Preparation
cgm_plot = cgm.copy()
cgm_plot["effective_time_frame_date_time_local"] = pd.to_datetime(
    cgm_plot["effective_time_frame_date_time_local"]
)
cgm_plot["date"] = cgm_plot["effective_time_frame_date_time_local"].dt.date.astype(str)
cgm_plot["hour"] = (
    cgm_plot["effective_time_frame_date_time_local"].dt.hour +
    cgm_plot["effective_time_frame_date_time_local"].dt.minute / 60
)

# 4. Classification Function


def classify_glucose(row):
    if row.blood_glucose_value < 54:
        return Category.very_low.value
    elif row.blood_glucose_value < 70:
        return Category.low.value
    elif row.blood_glucose_value < 180:
        return Category.target_range.value
    elif row.blood_glucose_value < 250:
        return Category.high.value
    else:
        return Category.very_high.value


cgm_plot["category"] = cgm_plot.apply(classify_glucose, axis=1)

# 5. Helper to Add Background Zones


def add_glucose_zones(ax):
    ax.axhspan(0, 54, facecolor="#a00", alpha=0.1)
    ax.axhspan(54, 70, facecolor="#f44", alpha=0.1)
    ax.axhspan(70, 180, facecolor="#CDECCD", alpha=0.3)
    ax.axhspan(180, 250, facecolor="#FDBE85", alpha=0.2)
    ax.axhspan(250, 400, facecolor="#FD8D3C", alpha=0.2)

# 6. Daily Subplots
cols = 5
dates = sorted(cgm_plot["date"].unique())[:10]
rows = (len(dates) + cols - 1) // cols
fig, axes = plt.subplots(rows, cols, figsize=(cols * 5, rows * 3), sharey=True)
axes = axes.flatten()
grouped = cgm_plot.groupby("date")

for i, date in enumerate(dates):
    ax = axes[i]
    day_data = grouped.get_group(date).sort_values("hour")
    add_glucose_zones(ax)

    smoothed = day_data["blood_glucose_value"].rolling(window=5, center=True, min_periods=1).median()
    ax.plot(day_data["hour"], smoothed, color="green", label="Median Glucose")

    # Fill high
    ax.fill_between(
        day_data["hour"], 180, day_data["blood_glucose_value"],
        where=day_data["blood_glucose_value"] > 180,
        color="#FDAE61", alpha=0.5
    )
    # Plot low
    low_mask = day_data["blood_glucose_value"] < 70
    ax.plot(
        day_data["hour"][low_mask],
        day_data["blood_glucose_value"][low_mask],
        color="#E41A1C", linewidth=1.5
    )

    weekday = datetime.datetime.strptime(date, "%Y-%m-%d").strftime("%A")
    ax.set_xlim(0, 24)
    ax.set_xticks([12])
    ax.set_xticklabels(["12 PM"])
    ax.set_ylim(0, 300)
    ax.set_yticks([54,70, 180, 250, 300])
    ax.set_yticklabels(["", "70", "180", "250", "300"])
    ax.set_title(f"{weekday}\n{date}")
    ax.axvline(x=12, color="gray", linestyle="--", linewidth=1, alpha=0.6)
    ax.grid(True)
    if i == 0:
        ax.legend(fontsize=8, loc="upper right", frameon=False)

for j in range(i + 1, len(axes)):
    axes[j].axis("off")

fig.suptitle("Daily CGM Profiles with Glucose Zones", fontsize=16, y=0.94)
plt.tight_layout()
plt.subplots_adjust(
    top=0.75,      # space for the title
    hspace=0.75,    # vertical space between rows (lower = tighter)
    wspace=0    # horizontal space between columns (try 0 for near-touching)
)
plt.show()

# 7. AGP (Ambulatory Glucose Profile)
cgm_plot["minute_of_day"] = (
    cgm_plot["effective_time_frame_date_time_local"].dt.hour * 60 +
    cgm_plot["effective_time_frame_date_time_local"].dt.minute
)
cgm_plot["time_bin"] = (cgm_plot["minute_of_day"] // 5) * 5

agp_summary = cgm_plot.groupby("time_bin")["blood_glucose_value"].quantile([0.05, 0.25, 0.5, 0.75, 0.95]).unstack()
agp_summary.columns = ["p5", "p25", "p50", "p75", "p95"]
agp_summary = agp_summary.reset_index()
agp_summary["hour"] = agp_summary["time_bin"] / 60

plt.figure(figsize=(12, 6))
add_glucose_zones(plt)
plt.fill_between(agp_summary["hour"], agp_summary["p5"], agp_summary["p95"], color="#D8E3E7", alpha=0.5, label="5–95% Range")
plt.fill_between(agp_summary["hour"], agp_summary["p25"], agp_summary["p75"], color="#B0C4DE", alpha=0.6, label="25–75% Range")
plt.plot(agp_summary["hour"], agp_summary["p50"], color="green", linewidth=2, label="Median")

plt.title("Ambulatory Glucose Profile (AGP)")
plt.xlabel("Time of Day")
plt.ylabel("Glucose (mg/dL)")
plt.xticks([0, 6, 12, 18, 24], ["12 AM", "6 AM", "12 PM", "6 PM", "12 AM"])
plt.xlim(0, 24)
plt.ylim(0, 300)
plt.yticks([54, 70, 180, 250, 300])
plt.grid(True)
plt.legend(loc="upper right", fontsize=10, frameon=True, framealpha=0.9)
plt.tight_layout()
plt.show()

# Count total readings
total = len(cgm_plot)

# Compute % time in each range
tir_data = {
    "Very Low (<54)": (cgm_plot["blood_glucose_value"] < 54).sum() / total * 100,
    "Low (54–69)": ((cgm_plot["blood_glucose_value"] >= 54) & (cgm_plot["blood_glucose_value"] < 70)).sum() / total * 100,
    "Target (70–180)": ((cgm_plot["blood_glucose_value"] >= 70) & (cgm_plot["blood_glucose_value"] <= 180)).sum() / total * 100,
    "High (181–250)": ((cgm_plot["blood_glucose_value"] > 180) & (cgm_plot["blood_glucose_value"] <= 250)).sum() / total * 100,
    "Very High (>250)": (cgm_plot["blood_glucose_value"] > 250).sum() / total * 100,
}

# Plot
plt.figure(figsize=(8, 5))
tir_colors = [
    "#a00", "#f44", "#CDECCD", "#FDBE85", "#FD8D3C"
]
bars = plt.barh(list(tir_data.keys()), list(tir_data.values()), color=tir_colors)

plt.xlabel("Time in Range (%)")
plt.title("Time in Glucose Ranges")
plt.xlim(0, 100)
plt.grid(axis="x", linestyle="--", alpha=0.5)

# Add percentage labels
for bar in bars:
    width = bar.get_width()
    plt.text(width + 1, bar.get_y() + bar.get_height()/2,
             f"{width:.1f}%", va='center', fontsize=9)

plt.tight_layout()
plt.show()

# Define the same time-in-range calculation
tir_data = {
    "Very Low": (cgm_plot["blood_glucose_value"] < 54).sum(),
    "Low": ((cgm_plot["blood_glucose_value"] >= 54) & (cgm_plot["blood_glucose_value"] < 70)).sum(),
    "Target": ((cgm_plot["blood_glucose_value"] >= 70) & (cgm_plot["blood_glucose_value"] <= 180)).sum(),
    "High": ((cgm_plot["blood_glucose_value"] > 180) & (cgm_plot["blood_glucose_value"] <= 250)).sum(),
    "Very High": (cgm_plot["blood_glucose_value"] > 250).sum(),
}

total = sum(tir_data.values())
tir_pct = {k: v / total * 100 for k, v in tir_data.items()}

# Colors that match your other plots
tir_colors = {
    "Very Low": "#a00",
    "Low": "#f44",
    "Target": "#CDECCD",
    "High": "#FDBE85",
    "Very High": "#FD8D3C",
}

# Plot the vertical stacked bar
plt.figure(figsize=(2, 6))

bottom = 0
for label in ["Very Low", "Low", "Target", "High", "Very High"]:
    height = tir_pct[label]
    plt.bar(0, height, bottom=bottom, color=tir_colors[label], width=0.5, edgecolor='white')
    # Add percent labels
    if height > 3:  # skip labeling tiny slivers
        plt.text(0.6, bottom + height / 2, f"{height:.0f}%", va='center', fontsize=9)
    bottom += height

# Formatting
plt.xlim(-0.5, 2)
plt.ylim(0, 100)
plt.xticks([])
plt.ylabel("Time in Range (%)")
plt.title("Glucose Zones")

# Optional: add reference lines
for y in [70, 180, 250]:  # these are just visual references, can be removed
    plt.axhline(y=y, color='gray', linestyle='--', linewidth=0.5, alpha=0.3)

plt.box(False)
plt.tight_layout()
plt.show()

df['Time'] = df.effective_time_frame_date_time_local
df['Glucose'] = df.blood_glucose_value
df['Day'] = df['Time'].dt.date
cgmquantify.plotglucosebounds(df)

cgm.head()

# Calculate differences between consecutive glucose values
cgm["glucose_diff"] = cgm["blood_glucose_value"].diff().abs()

# Define a threshold for significant glucose excursions (e.g., 1 SD)
threshold = std_glucose

# Compute MAGE: Mean of all large excursions above the threshold
mage = cgm[cgm["glucose_diff"] > threshold]["glucose_diff"].mean()
print(f"Mean Amplitude of Glycemic Excursions (MAGE): {mage:.2f} mg/dL")

# Extract hour and minute from the timestamp
cgm["time_of_day"] = cgm["effective_time_frame_date_time_local"].dt.strftime("%H:%M")

# Compute the mean glucose value for each time of day across all days
mean_per_time = cgm.groupby("time_of_day")["blood_glucose_value"].mean()

# Compute MODD as the mean of absolute day-to-day differences
modd = mean_per_time.diff().abs().mean()
print(f"Mean of Daily Differences (MODD): {modd:.2f} mg/dL")

all_study_data = jh_client.list_observations_df(study_id=study_id, limit=50_000, code=CGM)
all_study_data

all_study_data["patient_id"] = all_study_data["subject_reference"].str.replace("Patient/", "")

all_cgm = all_study_data.loc[
    :,
    [
        "patient_id",
        "blood_glucose_value",
        "effective_time_frame_date_time_local",
    ],
]
# ensure sorted by date
all_cgm = all_cgm.sort_values("effective_time_frame_date_time_local")
all_cgm.head()

all_cgm.groupby("patient_id").blood_glucose_value.mean().sort_values().plot(kind="bar");
plt.title("Mean blood glucose")
plt.ylabel("mg/dL");

def compute_gmi(mean_glucose):
    return 3.31 + (0.02392 * mean_glucose)

all_cgm.groupby("patient_id").blood_glucose_value.mean().apply(compute_gmi).sort_values().plot(kind="bar")
plt.title("Glucose Management Indicator (GMI)");

def compute_time_in_range(df):
    tir_data = {
        "Very Low": (df["blood_glucose_value"] < 54).sum(),
        "Low": ((df["blood_glucose_value"] >= 54) & (df["blood_glucose_value"] < 70)).sum(),
        "Target": ((df["blood_glucose_value"] >= 70) & (df["blood_glucose_value"] <= 180)).sum(),
        "High": ((df["blood_glucose_value"] > 180) & (df["blood_glucose_value"] <= 250)).sum(),
        "Very High": (df["blood_glucose_value"] > 250).sum(),
    }
    return pd.Series(tir_data) * 100 / len(df)

ranges = all_cgm.groupby("patient_id").apply(compute_time_in_range, include_groups=False)
ranges

ranges.plot(kind='bar', stacked=True, color=tir_colors)
plt.legend(bbox_to_anchor=(1, 1), reverse=True)
plt.ylim(0, 100)
plt.title("% time in range");