# From GitHub, for the latest version (CRAN release can lag)
install.packages("pak")
pak::pak("irinagain/iglu")
install.packages(c("dplyr", "quarto"))Continuous Glucose Monitoring (CGM) Tutorial
Digital Health Technologies: Data, Methods and Applications with R — CGM Module
CGMs are small wearable devices that measure interstitial glucose levels continuously throughout the day, with some monitors taking measurements as often as every 5 minutes. Data from these monitors provide a detailed quantification of the variation in glucose levels during the course of the day, and thus CGMs play an increasing role in clinical practice. For more on CGMs, see for example Rodbard (2016) “Continuous Glucose Monitoring: A Review of Successes, Challenges, and Opportunities.”.
This tutorial will focus on visualizing raw CGM data, calculating consensus metrics of glycemic control and some additional CGM metrics, and visualizing the data using the Ambulatory Glucose Profile (AGP) based on R package iglu Broll, Urbanek, Buchanan, Chun, Muschelli, Punjabi and Gaynanova (2021) “Interpreting blood glucose data with R package iglu.” PLoS One, 16(4), e0248560, with expanded functionality described in Chun, Fernandes and Gaynanova (2024) “An Update on the iglu Software Package for Interpreting Continuous Glucose Monitoring Data.” Diabetes Technology & Therapeutics, 26(12), 939-950.
We will then move beyond standard metrics into functional and distributional approaches to CGM data, sometimes referred to as “CGM Data Analysis 2.0” Klonoff, Bergenstal, Cengiz et al. (2025) “CGM Data Analysis 2.0: Functional Data Pattern Recognition and Artificial Intelligence Applications.” Journal of Diabetes Science and Technology.
1. Prerequisites and data format
Prerequisites
This tutorial relies on R package iglu with all its dependencies and also utilizes R package dplyr
Alternatively, you can follow via Shiny App here
Parts 5 and 6 rely on additional, more specialized packages, all available from CRAN. You will not need these until we get there. Those parts also load pre-processed Hall et al. (2018) data from hall_data.Rda, built by the companion script prepare_hall_data.qmd - run that once (it’s the slow part) before rendering Parts 5-6 of this tutorial. That script has its own, slightly different package needs (notably DBI/RSQLite, for reading the clinical SQLite database - not needed here, since this tutorial only ever reads the already-processed .Rda); see the top of prepare_hall_data.qmd for its install code.
# For Part 5 (Functional Data Analysis)
install.packages("refund")
install.packages("tidyr")
install.packages("stringr")
install.packages("patchwork")
# For Part 6 (Distributional analysis)
install.packages("biosensors.usc")
pak::pak("https://github.com/alexandercoulter/fastfrechet")
pak::pak("IrinaStatsLab/OptiThresholds/optithresholdr")Data format
Continuous Glucose Monitoring (CGM) data typically has the following three columns, often provided as .csv file.
Glucose level measurement [in mg/dL] (
"gl")Timestamp for glucose measurement (
"time")Subject identification (
"id")
iglu is designed to work with mg/dL unit of glucose measurements (US standard). In Europe and other countries, it’s more common to measure glucose in mmol/L, with conversion as follows: \[
mg/dL = 18.018 \times mmol/L
\]
2. CGM data visualization and data quality checks
Data visualization
We start with example_data_5_subject, 5 min frequency Dexcom G4 CGM data from 5 subjects with type 2 diabetes not on insulin therapy, part of a larger study analyzed in Gaynanova et al. (2020). example_data_1_subject is Subject 1 from this same study, provided separately for single-subject demos.
Dexcom G4 CGM measurements with 5 min frequency of a subject with type 2 diabetes.
plot_glu(example_data_1_subject)
We can also look at all subjects at once (all with type 2 diabetes)
plot_glu(example_data_5_subject)
Horizontal lines by default are at 70 and 180 mg/dL, thresholds representing target glucose levels (in range) for most patients with diabetes.
Data availability check
Based on current consensus recommendations Battelino et al. (2023), 14 days of monitoring with 70% non-missing data is a minimal data availability standard for outpatient studies to represent glycemic exposure over 3 month period. We can check data availability using active_percent
data_check = active_percent(example_data_5_subject)
data_check# A tibble: 5 × 5
id active_percent ndays start_date end_date
<fct> <dbl> <drtn> <dttm> <dttm>
1 Subject 1 79.8 12.7 days 2015-06-06 17:50:27 2015-06-19 09:59:36
2 Subject 2 58.9 16.7 days 2015-02-24 17:31:29 2015-03-13 10:38:01
3 Subject 3 92.1 5.8 days 2015-03-10 16:36:26 2015-03-16 11:11:05
4 Subject 4 98.7 12.9 days 2015-03-13 13:44:09 2015-03-26 11:01:58
5 Subject 5 95.8 10.6 days 2015-02-28 17:40:06 2015-03-11 09:04:28To check for consensus recommendation, we can look at number of days (ndays) times percent non-missing (active_percent), and compare with \(14 * 0.7 = 9.8\) threshold.
data_check %>%
dplyr::mutate(threshold = 14 * 0.7,
actual = ndays * active_percent/100,
meets_consensus = actual >= threshold) %>%
dplyr::select(id, ndays, active_percent, actual, threshold, meets_consensus)# A tibble: 5 × 6
id ndays active_percent actual threshold meets_consensus
<fct> <drtn> <dbl> <drtn> <dbl> <lgl>
1 Subject 1 12.7 days 79.8 10.139825 days 9.8 TRUE
2 Subject 2 16.7 days 58.9 9.838463 days 9.8 TRUE
3 Subject 3 5.8 days 92.1 5.343389 days 9.8 FALSE
4 Subject 4 12.9 days 98.7 12.729760 days 9.8 TRUE
5 Subject 5 10.6 days 95.8 10.152259 days 9.8 TRUE All subjects meet the threshold except Subject 3. For illustrative purposes, we would still keep that subject in later analyses, but in practice, one might want to exclude that subject or do sensitivity checks.
Tip
Whether the 14-day/70% rule is adequate depends on what you are trying to represent and which metric you need. Recent work by Kok, Williamson, Lee and Gaynanova (2026) “Impact of missing data and monitoring duration on downstream analyses in continuous glucose monitoring.” Diabetes Care, shows that the required monitoring duration varies by metric, and that closer to a month of data may be needed when the goal is to approximate a longer-term (e.g. 90-day, HbA1c-comparable) glycemic exposure rather than a shorter snapshot.
Caution
The examples above are curated study data, used more or less “as-is.” CGM data pulled from research data warehouses or biobanks is typically messier - duplicate readings, device/format heterogeneity, and other artifacts that can silently bias downstream metrics if not addressed before analysis. See Williamson, Lee and Gaynanova (2026) “A Processing Algorithm to Address Real-World Data Quality Issues With Continuous Glucose Monitoring Data.” Journal of Diabetes Science and Technology, 20(4), 1376-1380, for a processing algorithm addressing this.
3. CGM metrics: consensus and beyond
Consensus CGM metrics described in Battelino et al. (2023)
- Time-in-Range (TIR, 70-180 mg/dL), Time in Hypoglycemia (Level 1 and Level 2), Time in Hyperglycemia (Level 1 and Level 2)
- CV (Coefficient of Variation)
- Mean glucose and GMI (Glucose Management Indicator)
- GRI (Glucose Risk Index)
- Glycemic Episodes
Time in range metrics
Time in range (TIR)
Most common and accepted metric as treatment target. Calculated as % of time spend within fixed thresholds of \([70, 180]\) mg/dL.
in_range_percent(example_data_5_subject,
target_ranges = list(c(70, 180)))# A tibble: 5 × 2
id in_range_70_180
<fct> <dbl>
1 Subject 1 91.7
2 Subject 2 26.4
3 Subject 3 81.3
4 Subject 4 95.1
5 Subject 5 62.1A typical goal is have TIR over 70%. Again, we see that Subjects 2 and 5 are outside of this goal. Subjects without diabetes typically have over 95% TIR.
Time in range can also be judged from the plots.
plot_glu(example_data_5_subject, LLTR = 80, ULTR = 180)
Sometimes, a tight time in range (TTIR) is used (e.g. pregnancy) which corresponds to [70, 140] mg/dL thresholds.
Glycemic thresholds
It is typical to divide the whole range of measurements into time spent within prespecified thresholds
- Level 2 Hypoglycemia [< 54 mg/dL]
- Level 1 Hypoglycemia [54 - 70 mg/dL]
- In-range [70 - 180 mg/dL]
- Level 1 Hyperglycemia [180 - 250 mg/dL]
- Level 2 Hyperglycemia [> 250 mg/dL]
The sum across ranges is 100%, giving rise to barplot. Sometimes, Levels 1 and 2 are combined, giving rise to just 3 areas.
plot_glu(example_data_5_subject %>%
dplyr::filter(id == "Subject 5"))
plot_ranges(example_data_5_subject %>%
dplyr::filter(id == "Subject 5"))
For Subject 5, there is significant time in Hyperglycemia, and no time in Hypoglycemia. In general, Hypoglycemia is more prominent in subjects with Type 1 diabetes.
The ranges can be evaluated separately with any thresholds
below_percent(example_data_5_subject %>%
dplyr::filter(id == "Subject 5"))# A tibble: 1 × 3
id below_54 below_70
<fct> <dbl> <dbl>
1 Subject 5 0 0.103above_percent(example_data_5_subject %>%
dplyr::filter(id == "Subject 5"))# A tibble: 1 × 4
id above_140 above_180 above_250
<fct> <dbl> <dbl> <dbl>
1 Subject 5 69.8 37.8 11.3Glycemic variability metrics
Coefficient of variation
CV is a global measure of variability (mean/sd). A typical treatment target is below 36%.
cv_glu(example_data_5_subject)# A tibble: 5 × 2
id CV
<fct> <dbl>
1 Subject 1 26.9
2 Subject 2 24.0
3 Subject 3 29.1
4 Subject 4 22.4
5 Subject 5 33.5All 5 subjects meet this treatment goal.
Other variability measures
Different types of standard deviation (SD), between days, within time points, all are highly correlated but can be computed at once.
sd_measures(example_data_5_subject)# A tibble: 5 × 7
id SDw SDhhmm SDwsh SDdm SDb SDbdm
<fct> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 Subject 1 26.4 19.6 6.54 16.7 27.9 24.0
2 Subject 2 36.7 22.8 7.62 52.0 48.0 35.9
3 Subject 3 42.9 14.4 9.51 12.4 42.8 42.5
4 Subject 4 24.5 12.9 6.72 16.9 25.5 22.0
5 Subject 5 50.0 29.6 12.8 23.3 50.3 45.9MAGE (Mean Amplitude of Glycemic Excursions)
Uses automatic peak and nadir identification algorithm to compute amplitude of excursions, and then takes an average amplitude for all excursions above 1 standard deviation.
mage(example_data_5_subject %>%
dplyr::filter(id == "Subject 3"), plot = TRUE)
# numeric
mage(example_data_5_subject)# A tibble: 5 × 2
# Rowwise:
id MAGE
<fct> <dbl>
1 Subject 1 72.4
2 Subject 2 118.
3 Subject 3 116.
4 Subject 4 70.9
5 Subject 5 142. Rate of change
Rate of change is a discrete approximation of 1st derivative associated with glucose curve (discretized using 15 min). Unlike glucose measurements themselves, rate of change is symmetric, and standard deviation (SD) of rate of change can be used a measure of local glucose variability.
# visual
hist_roc(example_data_1_subject)
# numeric
sd_roc(example_data_5_subject)# A tibble: 5 × 2
id sd_roc
<fct> <dbl>
1 Subject 1 0.620
2 Subject 2 0.642
3 Subject 3 0.831
4 Subject 4 0.617
5 Subject 5 1.05 Other metrics and all together
Mean and GMI
Mean values match our intuition from data visualization, all are in mg/dL.
mean_glu(example_data_5_subject)# A tibble: 5 × 2
id mean
<fct> <dbl>
1 Subject 1 124.
2 Subject 2 218.
3 Subject 3 154.
4 Subject 4 130.
5 Subject 5 175.GMI is a deterministic transformation of mean on HbA1c scale \[ GMI = 3.31 + 0.02392 \times \text{mean glucose} \]
gmi(example_data_5_subject)# A tibble: 5 × 2
id GMI
<fct> <dbl>
1 Subject 1 6.27
2 Subject 2 8.54
3 Subject 3 6.99
4 Subject 4 6.41
5 Subject 5 7.49HbA1c is a measure of average glucose over the past 3 months
Pre-diabetes - A1c of 5.7%-6.4%; Diabetes - A1c> 6.5%
Typical treatment goal: A1c < 7%
Based on the above, Subjects 2 and 5 have the highest glucose levels and outside of treatment goal.
Note
The linear GMI formula above can over- or under-estimate HbA1c. An updated, non-linear GMI (uGMI) has recently been proposed to better align with HbA1c: Bergenstal, Xu, Dunn et al. (2026) “Updated glucose management indicator (GMI) better aligns with HbA1c than current GMI: implications for clinical practice and reporting.” Diabetologia, 69, 2182-2188.
GRI (Glucose Risk Index)
More recent measure, based on Principal Component Analysis of clinicians’ ratings of CGM profiles. The final formula is based on percentages within each level - attempt to arrive at one-number summary of glycemic control that is reflective of both hypo and hyperglycemia. \[ GRI = 3 \times \text{ Lv2 Hypo} + 2.4 \times \text{ Lv1 Hypo} + 0.8 \times \text{ Lv1 Hyper} + 1.6 \times \text{ Lv2 Hyper} \]
gri(example_data_5_subject)# A tibble: 5 × 2
id GRI
<fct> <dbl>
1 Subject 1 7.19
2 Subject 2 79.7
3 Subject 3 20.0
4 Subject 4 4.38
5 Subject 5 39.5 GRI \(=0\) indicates time-in-range of 100%. Maximum allowable GRI is 100%.
Glycemic Episodes
Glycemic episodes are defined based on consecutive measurements within a certain range.
- High Glucose (Level 1) - \(>180\), at least 15 consecutive min, episode ends when \(\geq 15\) consecutive min of values \(< 180\)
- Very High Glucose (Level 2) - \(> 250\), at least 15 consecutive min, episode ends when \(\geq 15\) consecutive min of values \(< 250\)
Similar for low glucose and very low glucose.
Dexcom CGMs have 5 measurement frequency, hence require 3 consecutive readings.
We illustrate with example_data_hall, Dexcom G4 5 min frequency CGM data from 19 subjects with pre-diabetes and type 2 diabetes from Hall et al. (2018).
epicalc_profile(example_data_hall %>%
dplyr::filter(id == "2133-039"))
Alternative numeric output directly
episode_calculation(example_data_hall %>%
dplyr::filter(id == "2133-039"))# A tibble: 7 × 7
id type level avg_ep_per_day avg_ep_duration avg_ep_gl total_episodes
<chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl>
1 2133-039 hypo lv1 1.33 49 63.8 10
2 2133-039 hypo lv2 0.133 15 51.4 1
3 2133-039 hypo extend… 0 0 NA 0
4 2133-039 hyper lv1 0.266 42.5 188. 2
5 2133-039 hyper lv2 0 0 NA 0
6 2133-039 hypo lv1_ex… 1.20 45 64.2 9
7 2133-039 hyper lv1_ex… 0.266 42.5 188. 2
Caution
Episode counts are more sensitive to algorithmic choices than they first appear - e.g. how missing-data gaps within a candidate episode are handled, or whether short dips below threshold are allowed to interrupt it. Different software can disagree substantially on the resulting counts for the same raw data. See Gaynanova and Lee (2025) “When Algorithms Diverge: Quantification of Glycemic Episodes from Continuous Glucose Monitor Data.” Diabetes Technology & Therapeutics, 27(6), 500-502.
All consensus metrics at once
All consensus metrics can be computed at once with all_metrics call.
all_metrics(example_data_5_subject,
metrics_to_include = "consensus_only")# A tibble: 5 × 18
id below_54 below_70 in_range_70_180 above_180 above_250 SD mean CV
<fct> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 Subje… 0 0.137 91.7 8.20 0.377 33.3 124. 26.9
2 Subje… 0 0 26.4 73.6 26.1 52.4 218. 24.0
3 Subje… 0 0.326 81.3 18.3 5.68 44.8 154. 29.1
4 Subje… 0.0546 0.273 95.1 4.61 0 29.1 130. 22.4
5 Subje… 0 0.103 62.1 37.8 11.3 58.6 175. 33.5
# ℹ 9 more variables: active_percent <dbl>, ndays <drtn>, start_date <dttm>,
# end_date <dttm>, in_range_70_140 <dbl>, GMI <dbl>, GRI <dbl>,
# total_extended_hypo_episodes <dbl>, total_extended_hyper_episodes <dbl>Out of all these metrics, TIR is considered the 1st default and most commonly used metric.
Metrics heatmap
Multiple additional metrics of glycemic control based on CGM data exist. These can be visualized as a heatmap for a particular study.
cluster_out = metrics_heatmap(data = example_data_hall)
Tip
Heatmap can be useful to narrow down the metrics list to a few that offer the most distinct information. Based on experience across multiple datasets, a complementary shortlist is: TIR for overall control, TBR for hypoglycemia, GMI (tends to track well with TAR) for average exposure, CV for overall variability, and MAGE and SD of rate of change for more local variability.
4. Advanced visualization of CGM data
Now that Part 3 has introduced both the glycemic threshold levels (Level 2 Hypoglycemia, Level 1 Hypoglycemia, In-range, Level 1 Hyperglycemia, Level 2 Hyperglycemia) and glycemic variability metrics, we can make use of visualizations that build on both.
AGP (Ambulatory Glucose Profile)
Most consensus metrics are typically summarized together in a single Ambulatory Glucose Profile. Internally, agp() composes a stats table, a range/TIR plot, and the percentile-band plot into one figure.
agp(example_data_1_subject, daily = FALSE)
Lasagna plots
An alternative way to visualize CGM data across many days (or many subjects) at once is a lasagna plot Swihart et al. (2010), where each row is a day (or subject) and color encodes glycemic range.
plot_lasagna_1subject(example_data_1_subject,
color_scheme = "red-orange")
The plot can also be time-sorted to highlight average 24-hour patterns.
plot_lasagna_1subject(example_data_1_subject,
lasagnatype = 'timesorted',
color_scheme = "red-orange")
Lasagna plots also make missing-data handling directly visible, since unfilled gaps show up as blank tiles. Like most iglu functions, plot_lasagna_1subject() linearly interpolates across gaps up to inter_gap minutes (default 45); larger gaps are left as missing rather than filled in. Widening that window fills in more of the plot, but at the cost of interpolating over longer stretches of genuinely unobserved data - worth doing deliberately, not by leaving the default unexamined (see also the data-quality callouts in Part 2).
With the default inter_gap = 45, some gaps remain blank:
plot_lasagna_1subject(example_data_1_subject,
color_scheme = "red-orange",
inter_gap = 45)
With inter_gap = 300 (5 hours), those same gaps are now filled in by interpolation:
plot_lasagna_1subject(example_data_1_subject,
color_scheme = "red-orange",
inter_gap = 300)
5. Functional data analysis
Standard CGM metrics summarize each subject’s glucose profile into one or a handful of numbers. Functional Data Analysis (FDA) approaches instead treat the CGM trajectory itself (or a meaningfully aligned segment of it, e.g. by sleep onset or meal time) as the unit of analysis, preserving its temporal shape.
For illustration, we use the full glycemic-response-to-meals dataset from Hall et al. (2018) (“S6 Data” of that paper - CC BY 4.0 licensed), rather than the 3-subject demo bundled with iglu (example_meals_hall). This file (HallFullData/pbio.2005143.s015.tsv) covers 30 subjects and 176 standardized-meal windows (3 meal types - Cereal Flakes “CF”, Peanut Butter Sandwich “PB”, Protein Bar “Bar” - each replicated twice), with glucose sampled from 30 minutes before to 2.5 hours after each meal.
Constructing meal-aligned glucose curves
Building meal_curves from the raw S6 file involves several cleaning steps - parsing timestamps, handling the “Low” sensor-floor convention, splitting meal type/replicate out of the Meal column, deduplicating repeated readings, and computing time relative to meal onset (each window starts 30 minutes before the meal, so onset = first sampled time + 30 minutes; the small number of meals with a substantially incomplete window are dropped). That wrangling is slow to redo on every render, so it now lives in the companion script prepare_hall_data.qmd, which saves the result (along with everything Part 6 needs below) to hall_data.Rda. Here we just load it:
load("hall_data.Rda")
dplyr::n_distinct(meal_curves$id)[1] 30dplyr::n_distinct(paste(meal_curves$id, meal_curves$Meal))[1] 160Now the illustration analogous to plot_glu, but with the x-axis re-centered on meal onset instead of clock time, one (semi-transparent) line per meal, faceted by meal type:
ggplot(meal_curves, aes(x = rel_time_hr, y = gl, group = interaction(id, Meal))) +
geom_vline(xintercept = 0, linetype = "dashed", color = "grey30") +
geom_line(alpha = 0.3, color = "steelblue") +
facet_wrap(~ meal_type) +
labs(x = "Time relative to meal start (hours)", y = "Glucose (mg/dL)",
title = "Meal-aligned glucose curves, Hall et al. (2018) full cohort") +
theme_minimal()
Each curve is now a comparable functional observation on a common domain (rel_time_hr), which is the input PCA on aligned glucose curves and Function-on-scalar regression below will build on.
PCA on aligned glucose curves
To keep this a simple one-curve-per-subject illustration (rather than the full multilevel/repeated-measures structure - see the note at the end of this section), prepare_hall_data.qmd restricts to the Cereal Flakes (“CF”) meal type and keeps each subject’s earliest CF occurrence only, so every subject contributes exactly one curve - already loaded above as the fpca_mat/subject_ids matrix (subjects x relative-time grid).
library(refund)
library(tidyr)
library(stringr)
library(patchwork)
dim(fpca_mat)[1] 29 37We fit the FPCA (Functional Principal Component Analysis) model to these meal-aligned glucose curves using refund::fpca.face(). Here we fix npc = 3 to keep a fixed, simple 3-component decomposition for interpretation below. Because our functional domain here is short (37 time points, -30 to +150 minutes), we also reduce knots well below the fpca.face default of 35, which was tuned for much longer/denser curves.
cf_fpca = fpca.face(fpca_mat, npc = 3, knots = 10)
cf_fpca$npc[1] 3cf_fpca$pve[1] 0.9581952Next, we plot the estimated eigenfunctions and a scree plot of the proportion of variance explained by each component.
eigs_df = as.data.frame(cf_fpca$efunctions)
colnames(eigs_df) = paste0("Phi", seq_len(cf_fpca$npc))
eigs_df$rel_min = as.numeric(stringr::str_remove(colnames(fpca_mat), "t_"))
eigen_plot = eigs_df %>%
tidyr::pivot_longer(starts_with("Phi"),
names_to = "component",
values_to = "value") %>%
dplyr::filter(component %in% paste0("Phi", 1:3)) %>%
ggplot(aes(x = rel_min / 60, y = value, colour = component)) +
geom_vline(xintercept = 0, linetype = "dashed", color = "grey60") +
geom_line() +
labs(x = "Time relative to meal start (hours)", y = "Eigenfunction value",
title = "Leading FPCs (CF meal, first occurrence per subject)") +
theme_minimal() +
theme(legend.position = c(0.75, 0.85))
evalues = cf_fpca$evalues
var_explained = evalues / sum(evalues)
scree_plot = tibble::tibble(k = seq_along(evalues),
FVE = var_explained,
cumFVE = cumsum(var_explained)) %>%
ggplot(aes(x = k, y = FVE)) +
geom_point() +
geom_line() +
labs(x = "Component", y = "Percent variance explained",
title = "Scree plot") +
theme_minimal()
eigen_plot + scree_plot
The eigenfunctions above are on an abstract, hard-to-interpret scale. A more directly interpretable view is to show the overall mean glucose curve, and how adding/subtracting each PC shifts it - i.e. what a “high score” vs. “low score” on that PC actually looks like in mg/dL.
pc_effect_df = dplyr::bind_rows(lapply(seq_len(cf_fpca$npc), function(k) {
effect = sqrt(cf_fpca$evalues[k]) * cf_fpca$efunctions[, k]
tibble::tibble(rel_min = eigs_df$rel_min,
component = paste0("PC", k),
`Mean` = cf_fpca$mu,
`Mean + 1 SD` = cf_fpca$mu + effect,
`Mean - 1 SD` = cf_fpca$mu - effect)
}))
pc_effect_df %>%
tidyr::pivot_longer(cols = c(`Mean`, `Mean + 1 SD`, `Mean - 1 SD`),
names_to = "curve", values_to = "gl") %>%
dplyr::mutate(curve = factor(curve, levels = c("Mean - 1 SD", "Mean", "Mean + 1 SD"))) %>%
ggplot(aes(x = rel_min / 60, y = gl, color = curve, linetype = curve)) +
geom_vline(xintercept = 0, linetype = "dotted", color = "grey60") +
geom_line(linewidth = 0.9) +
facet_wrap(~ component) +
scale_color_manual(values = c("Mean - 1 SD" = "#D55E00", "Mean" = "black", "Mean + 1 SD" = "#0072B2")) +
scale_linetype_manual(values = c("Mean - 1 SD" = "dashed", "Mean" = "solid", "Mean + 1 SD" = "dashed")) +
labs(x = "Time relative to meal start (hours)", y = "Glucose (mg/dL)",
color = NULL, linetype = NULL,
title = "Mean CF-meal response, and the effect of each PC") +
theme_minimal()
Each panel now reads directly in glucose units: the black line is the average CF-meal response across subjects, and the colored lines show how a subject scoring high (blue) vs. low (orange) on that particular PC deviates from the average - e.g. a PC largely shifting the post-meal peak height/timing, vs. one shifting the baseline/recovery level.
Note
This uses one curve per subject for simplicity. The data actually has up to 6 curves per subject (3 meal types x 2 replicates each - see previous section), a repeated-measures/multilevel structure. A natural extension is multilevel FPCA across all meal curves, see for example Gaynanova, Punjabi and Crainiceanu (2022) “Modeling continuous glucose monitoring (CGM) data during sleep.” Biostatistics, 23(1), 223-239 (developed for repeated within-subject sleep periods, directly applicable to repeated within-subject meal challenges here).
Function-on-scalar regression
We use the same curves as the PCA above (fpca_mat/subject_ids, CF meal, first occurrence per subject), with A1C as a single covariate pulled from the study’s clinical data (“S5 Data” of Hall et al. (2018), a SQLite database) and already joined into subject_covs (loaded above from hall_data.Rda). Note that this file’s userID column mixes two different ID formats across its 57 rows; conveniently, the CF-meal subjects we already have curves for all use the same "2133-###" format as our data, so no ID remapping was needed for them.
We use A1C (HbA1c) as a single, clinically standard continuous covariate. Two subjects are missing A1C, so we restrict to the 27 subjects with a value.
refund::pffr() fits a penalized function-on-scalar regression, estimating a smooth coefficient function \(\beta(t)\) - i.e. how much a one-unit (1 percentage point) increase in A1C shifts the expected glucose curve at each relative-time point \(t\), rather than a single number.
rel_min_grid = as.numeric(stringr::str_remove(colnames(fpca_mat), "t_"))
complete = !is.na(subject_covs$A1C)
sum(complete)[1] 27fosr_data = list(Y = fpca_mat[complete, , drop = FALSE],
A1C = subject_covs$A1C[complete])
fosr_fit = pffr(Y ~ A1C, yind = rel_min_grid, data = fosr_data)
summary(fosr_fit)
Family: gaussian
Link function: identity
Formula:
Y ~ A1C
Constant coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -150.16 80.05 -1.876 0.061 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Smooth terms & functional coefficients:
edf Ref.df F p-value
Intercept(yindex) 11.081 19.000 3.198 0.0332 *
A1C(yindex) 4.002 4.157 2.828 0.0226 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
R-sq.(adj) = 0.492 Deviance explained = 50%
-REML score = 4554.8 Scale est. = 796.13 n = 953 (27 x 37)
Sandwich correction (cluster) applied to covariance matrices.plot(fosr_fit, select = 2, rug = FALSE,
ylab = "A1C", ylim = c(-50, 150),
main = "A1C effect on CF-meal glucose response", xlab = "Meal time")
This is mgcv’s native plotting output (via pffr’s underlying GAM machinery), not ggplot2 like the rest of this tutorial - select = 2 shows only the estimated A1C coefficient function (dropping the functional-intercept panel), with a pointwise confidence band. A band that excludes zero at a given relative time indicates A1C has a detectable effect on glucose at that point in the meal response - e.g. only during the post-meal peak, or persisting through recovery.
Note
This is a minimal illustration with 1 covariate and 27 subjects. Sergazinov, Leroux, Cui, Crainiceanu, Aurora, Punjabi and Gaynanova (2023) Biometrics, 79(4), 3873-3882, develop fast function-on-scalar regression inference for the repeated-measures case (multiple curves per subject), directly applicable if using all meal replicates rather than one curve per subject.
6. Distributional analysis
Rather than fixing thresholds (70/180 mg/dL, etc.) or aligning curves in time, distributional approaches represent each subject’s CGM data as a full probability distribution of glucose values (a “glucodensity”), providing a threshold-independent summary that simultaneously captures the mean, TIR, and other consensus metrics Matabuena, Petersen, Vidal and Gude (2021) “Glucodensities: A new representation of glucose profiles using distributional data analysis.” Statistical Methods in Medical Research, 30(6), 1445-1464.
For this section we use all subjects from the full Hall et al. (2018) cohort, ignoring meal timing entirely: each subject’s entire monitoring-period glucose stream becomes their empirical distribution. This is “S1 Data” of the Hall paper (raw CGM recordings for all 57 subjects, not the 19-subject example_data_hall subset bundled with iglu, and not restricted to the meal windows used in Part 5). Reading and cleaning this ~105K-row file, and pulling clinical covariates for all 57 subjects, is handled by prepare_hall_data.qmd as well (same script as Part 5). Loading again here defensively, in case you’re jumping straight to Part 6:
load("hall_data.Rda")
dplyr::n_distinct(hall_full$id)[1] 57Clustering using biosensors package
Unlike fpca.face/pffr, biosensors.usc::load_data() expects its inputs as CSV files - raw long-format data (id, time, value) and a covariates file (id plus any covariates). Building the non-parametric density/quantile estimates for all 57 subjects is the slowest step in this whole section, so prepare_hall_data.qmd writes those CSVs, calls load_data() once, and saves the result as hall_biosensor (already loaded above).
library(biosensors.usc)
names(hall_biosensor)[1] "data" "densities" "quantiles" "variables"load_data() converts the raw readings into both a non-parametric density estimate ($density) and an empirical quantile estimate ($quantiles) per subject. clustering() then performs energy-distance clustering with the 2-Wasserstein metric on these distributional representations, and plots the resulting clusters. We ask for 3 clusters, matching the number of categories in the Hall study’s own glucotype classification (low/moderate/severe). clustering() has no seed argument of its own but uses randomized restarts internally, so we set R’s seed beforehand for reproducible cluster labels (cluster 1/2/3 assignments would otherwise not be guaranteed stable across re-runs):
set.seed(1)
clus = clustering(hall_biosensor, clusters = 3)
cluster_assignments = clustering_prediction(clus, hall_biosensor$quantiles$data)
tibble::tibble(id = hall_biosensor$variables$id,
cluster = cluster_assignments) %>%
dplyr::left_join(clinical_all, by = "id") %>%
dplyr::count(cluster, glucotype) %>%
tidyr::pivot_wider(names_from = glucotype, values_from = n, values_fill = 0)# A tibble: 3 × 4
cluster low moderate severe
<int> <int> <int> <int>
1 1 6 7 0
2 2 0 3 22
3 3 0 19 0Regression using fastFrechet
fastfrechet::frechetreg_univar2wass() performs Fréchet regression Petersen and Müller (2019) “Fréchet regression for random objects with Euclidean predictors.” Annals of Statistics, 47(2), 691-719 - when the response is a full distribution per subject (represented as an empirical quantile function, EQF, on a shared probability grid) rather than a single number - a distributional generalization of linear regression. The fastfrechet package itself is described in Coulter, Lee and Gaynanova (2025) “fastfrechet: An R package for fast implementation of Fréchet regression with distributional responses.” Journal of Open Source Software, 10(109), 7925. As with Part 5, we use A1C as the covariate.
library(fastfrechet)
m = 100
mseq = seq(0.5 / m, 1 - 0.5 / m, length.out = m)
subject_eqf = hall_full %>%
dplyr::group_by(id) %>%
dplyr::summarize(qf = list(stats::quantile(gl, probs = mseq, na.rm = TRUE)), .groups = "drop")
Y = do.call(rbind, subject_eqf$qf)
rownames(Y) = subject_eqf$id
dim(Y)[1] 57 100subject_a1c = tibble::tibble(id = subject_eqf$id) %>%
dplyr::left_join(clinical_all %>% dplyr::select(id, A1C), by = "id")
complete = !is.na(subject_a1c$A1C)
sum(complete)[1] 55X = matrix(subject_a1c$A1C[complete], ncol = 1)
colnames(X) = "A1C"
Y_complete = Y[complete, , drop = FALSE]Because most CGMs have 40 as lower measurement range, we constrain the fitted quantile functions with lower = 40. Rather than only recovering each observed subject’s own fitted distribution, we also predict fitted distributions (Z) at three representative A1C values (10th/50th/90th percentile among our subjects), to see how the whole distribution shape shifts with A1C:
a1c_grid = stats::quantile(X[, "A1C"], probs = c(0.1, 0.5, 0.9))
Z = matrix(a1c_grid, ncol = 1)
colnames(Z) = "A1C"
frechet_fit = frechetreg_univar2wass(X = X, Y = Y_complete, Z = Z, lower = 40)
dim(frechet_fit$Qhat)[1] 3 100tibble::tibble(p = rep(mseq, times = nrow(Z)),
A1C = factor(rep(round(a1c_grid, 1), each = m)),
glucose = as.vector(t(frechet_fit$Qhat))) %>%
ggplot(aes(x = p, y = glucose, color = A1C)) +
geom_line(linewidth = 0.9) +
labs(x = "Probability", y = "Glucose (mg/dL)",
color = "A1C (%)",
title = "Fitted glucose quantile functions by A1C level") +
theme_minimal()
The x-axis is probability rather than time, and each line is the entire fitted glucose distribution (not just its mean) for a subject at that A1C level. Compare the three curves to see whether higher A1C shifts the whole distribution upward, or mainly the upper tail (i.e. only the already-high glucose values) - both patterns are physiologically plausible, and this plot answers it empirically for this cohort.
Note
This is Fréchet regression without variable selection - a direct analogue of simple linear regression, but for distributional responses. Coulter, Aurora, Punjabi and Gaynanova (2025) develop the fast variable selection extension (FRiSO_univar2wass() and related functions in fastfrechet) for choosing among many candidate covariates simultaneously, which would be the natural next step using the fuller set of covariates.
Optimal / data-driven thresholds
Park, Kok and Gaynanova (2026) “Beyond fixed thresholds: optimizing summaries of wearable device data via piecewise linearization of quantile functions.” Statistics in Medicine, 45(15-17), e70646, asks whether the consensus 70/180 mg/dL cutoffs are actually the best choice for summarizing a given population into time-in-range-style bins, or whether other cutoffs would better preserve the information in the original distributions. The method (and its R port, OptiThresholdR) works directly from each subject’s raw measurements - no meal alignment or manual quantile construction needed - so we apply it straight to hall_full, the same 57-subject pooled data used above.
library(OptiThresholdR)
dist = as_distribution(hall_full, range = c(40, 400))We ask for K = 2 fully data-driven thresholds (matching the “one cutoff for hypoglycemia, one for hyperglycemia” structure of the consensus 70/180 pair), using loss1 - the distribution-preservation criterion appropriate here since we’re treating this as one cohort, not comparing subgroups.
opt_fit = optimal_thresholds(dist, K = 2, loss = "loss1", seed = 1)
opt_fit$cutoffs[1] 83.18805 133.34972opt_fit$objective[1] 110.4828For comparison, we evaluate that same loss1 criterion at the consensus 70/180 cutoffs, by asking for zero additional free thresholds beyond those two fixed ones:
fixed_fit = optimal_thresholds(dist, K = 0, fixed = c(70, 180), loss = "loss1", seed = 1)
fixed_fit$objective[1] 468.259tibble::tibble(
thresholds = c("Data-driven (K = 2)", "Consensus (70, 180)"),
cutoffs = c(paste(round(opt_fit$cutoffs, 1), collapse = ", "), "70, 180"),
loss1_objective = c(opt_fit$objective, fixed_fit$objective)
)# A tibble: 2 × 3
thresholds cutoffs loss1_objective
<chr> <chr> <dbl>
1 Data-driven (K = 2) 83.2, 133.3 110.
2 Consensus (70, 180) 70, 180 468.ggplot(hall_full, aes(x = gl, group = id)) +
geom_density(color = "steelblue", alpha = 0.15, linewidth = 0.4) +
geom_vline(xintercept = opt_fit$cutoffs, color = "#D55E00", linewidth = 1) +
geom_vline(xintercept = c(70, 180), color = "black", linetype = "dashed", linewidth = 1) +
labs(x = "Glucose (mg/dL)", y = "Density",
title = "Individual glucose distributions, all 57 subjects overlaid",
subtitle = "Orange solid = data-driven (K=2, loss1); black dashed = consensus (70, 180)") +
theme_minimal()
Overlaying all 57 individual curves shows both where subjects agree (peaks lining up) and where they differ (spread/multiple modes across subjects) - and lets you see whether the data-driven cutoffs (orange) look like they’re doing a better job of catching where individual distributions actually separate than 70/180 (black) does.
Tip
We used loss1 (single-population, distribution-preservation) since all 57 subjects are being pooled together. If the goal were instead to find thresholds that best separate known subgroups (e.g. by diagnosis or glucotype), loss2 (distance preservation) is the more appropriate criterion - see the package documentation for details. The semi-supervised fixed argument used above for evaluation can also be used the other way around, to keep one clinically-mandated cutoff while optimizing the rest.
7. Conclusion
This tutorial walked through a recommended pipeline for analyzing CGM data, from raw trace to modern extensions:
- Visualize and check data quality (Part 2) - plot the raw trace, and check data availability/monitoring duration before trusting any downstream metric.
- Compute consensus metrics (Part 3) - TIR, CV, mean/GMI, GRI, and glycemic episodes, all in one call via
all_metrics(), plus the wider metric libraryigluoffers when consensus metrics alone are insufficient for your question. - Use richer visual summaries (Part 4) - AGP and lasagna plots, once the underlying threshold levels and variability metrics are understood.
- Go beyond single-number summaries (Parts 5-6) - when the research question calls for it, treat the glucose trajectory as a curve (FDA) or as a full distribution (glucodensities), rather than compressing it into fixed predefined thresholds.
A few themes worth carrying forward:
- Data quality first. How much data you need depends on what you’re trying to represent and which metric you need - 14 days at 70% non-missing is a minimum for many outpatient settings, but longer monitoring may be required to stabilize some metrics or to approximate longer-term exposure (Part 2). Similarly, shorter monitoring may be sufficient for some use cases, e.g. inpatient settings.
- No single metric tells the whole story. Consensus metrics are interpretable and clinically validated, but they were developed with type 1 diabetes in mind and compress away a lot of structure; the metrics heatmap (Part 3) and the FDA/distributional approaches (Parts 5-6) are ways of recovering some of that structure. See also Gaynanova (2022) “Digital biomarkers of glucose control - reproducibility challenges and opportunities.” ASA Biopharmaceutical Report, 29(1), 21-26, on the reproducibility challenges this metric heterogeneity creates.
- Match the tool to the question. Fixed thresholds and TIR are the right lens for clinical treatment targets; functional and distributional methods are the right lens for research questions about when and how glucose varies, not just how much time is spent where.
Where to go next:
- iglu website - full metric/function reference and vignettes (MAGE, AGP and episodes, lasagna plots)
- Awesome-CGM - curated public CGM datasets for practicing on your own
- Gaynanova lab CGM page - papers and software referenced throughout this tutorial
- GlucoBench - if your interest is glucose forecasting (machine learning/AI), which this tutorial did not cover