📊 JupyterHealth Demo: CGM Data Analysis (Researcher View)

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()

org_dict = {org['id']: org for org in jh_client.list_organizations()}
# children is not populated
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])

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

# show all the patients with study data I have access to:
print("Patients with data I have access to:")

for patient in jh_client.list_patients():
    consents = jh_client.get_patient_consents(patient['id'])
    print(f"[{patient['id']}] {patient['nameFamily']}, {patient['nameGiven']}: {patient['telecomEmail']}")
    for study in consents['studies']:
        print(f"  - [{study['id']}] {study['name']}")
    for study in consents['studiesPendingConsent']:
        print(f"  - (not consented) [{study['id']}] {study['name']}")
    if not consents['studies'] and not consents['studiesPendingConsent']:
        print("  (no studies)") 

# pick patient id, study id from above
study_id = 30006
patient_id = 40012

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

df.head()

df.iloc[0]

# Reduce data to relevant subset for cgm
assert (df.blood_glucose_unit == 'mg/dL').all()
# reduce data
cgm = df.loc[
    df.resource_type == CGM.value,
    [
        "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))
bars = plt.barh(list(tir_data.keys()), list(tir_data.values()), color=[
    "#a00", "#f44", "#CDECCD", "#FDBE85", "#FD8D3C"
])

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")

from cgmsandbox import (
    CGMViewer,
    MageJumpOverlay,
    MeanGlucoseOverlay,
    SleepCompositionExtension,
    SleepWindowOverlay,
    TimeInRangeOverlay,
    WakeupGlucoseOverlay,
    load_sleep_nights,
    load_step_count,
)

viewer = CGMViewer(
    source="client",    # read the JupyterHealth Exchange data frame format
    client_df=df,       # that frame: this participant's CGM, from Step 2
    gl_range=(0, 250),  # y-axis bounds in mg/dL, fitted to this participant
)

# Each entry adds one layer to the same trace. They are constructed but not yet
# applied, so the cell first shows the raw signal.
layers = [
    TimeInRangeOverlay(),  # colours the trace by in/below/above range
    MeanGlucoseOverlay(source="client", scope="daily"),  # mean of the visible window
]

# Uncomment these two lines and re-run to apply them:
# for layer in layers:
#     viewer.add_overlay(layer)

# One cell, two scales: switch "daily" to "full" and re-run to lay the whole record
# out as consecutive weeks instead of a single day.
viewer.show(view_mode="daily")

# Sleep is stored under IEEE 1752 rather than Open mHealth.
IEEE_W3ID = "https://w3id.org/ieee1752"

sleep_raw = jh_client.list_observations_df(
    patient_id=patient_id,
    study_id=study_id,
    code=f"{IEEE_W3ID}|{Code.SLEEP_STAGE_SUMMARY.value}",
    limit=10_000,
)

nights = load_sleep_nights(sleep_raw)
nights.head()

viewer = CGMViewer(source="client", client_df=df, gl_range=(0, 250))

viewer.add_overlay(SleepWindowOverlay(nights))
viewer.add_overlay(TimeInRangeOverlay())
viewer.add_extensions(SleepCompositionExtension(nights))

# The sleep windows also shade the trace in the full view. The per-night panel does
# not, because extensions draw beneath a single day, not beneath the multi-week layout.
viewer.show(view_mode="daily")

viewer = CGMViewer(source="client", client_df=df, gl_range=(0, 250))

# The same 1 SD threshold Step 3 used, so the jumps marked here are the ones MAGE averaged.
viewer.add_overlay(MageJumpOverlay(threshold=std_glucose))

# One marker per wake time, taken from the sleep nights loaded above. 7 hours is a lower bar
# than the overlay's 8 hour default, which would filter out this participant's shorter nights.
viewer.add_overlay(WakeupGlucoseOverlay(source="client", nights=nights, min_sleep_hours=7.0))

viewer.show(view_mode="daily")

# Your turn: fill in the two methods below, then run the cell.
# The query has the same shape as the sleep one above. omh:step-count:3.0 is an
# Open mHealth code, so it needs no system|code prefix.
step_raw = jh_client.list_observations_df(
    patient_id=patient_id,
    study_id=study_id,
    code="omh:step-count:3.0",
    limit=10_000,
)
steps = load_step_count(step_raw)


class StepCountOverlay:
    """Shade each visible day by how many steps were recorded that day."""

    def __init__(self, steps):
        # TODO 1: keep what draw() will need. The busiest day is a convenient
        #   denominator for turning raw counts into a shading scale.
        ...

    def draw(self):
        # TODO 2: walk viewer.iter_axes_by_time() and shade each day that overlaps
        #   the visible window. viewer.scale(daily_value, full_value) picks the
        #   right value at either scale, and viewer.add_overlay() applies it.
        ...

# Solution. A daily step total is a per-day scalar with no home on a mg/dL axis, so
# shading intensity carries it: a darker band means a busier day.


# The filled-in overlay the scaffold above asks for.
class StepCountOverlay:  # noqa: F811
    """Shade each visible day by how many steps were recorded that day."""

    def __init__(self, steps):
        self.steps = steps
        self.busiest = max(steps["steps"].max(), 1)

    def draw(self):
        for ax, start, end in self.viewer.iter_axes_by_time():
            for _, row in self.steps.iterrows():
                left, right = max(row["start"], start), min(row["end"], end)
                if left >= right:
                    continue
                ax.axvspan(
                    left,
                    right,
                    color="#f0a500",
                    alpha=0.05 + 0.25 * (row["steps"] / self.busiest),
                    zorder=0,
                )


viewer = CGMViewer(source="client", client_df=df, gl_range=(0, 250))
viewer.add_overlay(StepCountOverlay(steps))
# Full view, so every day's shading is visible at once rather than a single day's.
viewer.show(view_mode="full")