Introduction

Aim

This analysis was conducted as part of the broader GLYSIMI study, which investigates the role of the human small intestinal microbiota in shaping postprandial glucose responses (PPGR) to standardized meals.

Within this project, the focus was on characterising individual glucose response patterns. Specifically, the analysis aimed to assess (1) the consistency of glycaemic responses within individuals to repeated consumption of the same meal, and (2) difference in responses to different food products across participants.

Standardized Test Meals

Participants consumed a set of standardized breakfast and dinner meals designed to elicit postprandial glucose responses.

Breakfast meals (each containing ~50 g carbohydrates): - Whole grain wheat bread with cream cheese (~390 kcal; WG_Bread) - Refined wheat bread with cream cheese (~350 kcal; R_Bread) - Mixed dried fruit (banana and apple, 50:50 w/w) with yogurt drink (~275 kcal; Fruit_Yoghurt) - Gingerbread with yogurt drink (~285 kcal; Gingerbread_Yoghurt)

Dinner meals (~68 g carbohydrates): - Refined wheat pasta with tomato sauce, soy-based minced meat replacer, and vegetables (~680 kcal; Pasta) - Whole grain rice with soy-based meat replacer, ketjap sauce, and vegetables (~735 kcal; Nasi)

In addition, participants consumed a glucose control drink containing 50 g of glucose.

Methods

Continuous glucose measurements were obtained using the FreeStyle Libre sensor system. These data were integrated with dietary intake records collected via the Traqq application, where participants logged meal timing. An additional dataset provided information on the allocation of specific food products to participants across study days.

The analysis focused on quantifying postprandial responses using incremental area under the curve (iAUC) above baseline. In addition, glucose response curves were explored and compared across individuals and meal types to assess patterns in glycaemic variability.

Preparation

Attaching necessary libraries

knitr::opts_chunk$set(echo = TRUE)
library(tidyverse)

## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.2.0     ✔ readr     2.2.0
## ✔ forcats   1.0.1     ✔ stringr   1.6.0
## ✔ ggplot2   4.0.2     ✔ tibble    3.3.1
## ✔ lubridate 1.9.5     ✔ tidyr     1.3.2
## ✔ purrr     1.2.1     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

library(ggpubr)
library(pals)
library(kableExtra)

## 
## Attaching package: 'kableExtra'
## 
## The following object is masked from 'package:dplyr':
## 
##     group_rows

library(rstatix)

## 
## Attaching package: 'rstatix'
## 
## The following object is masked from 'package:stats':
## 
##     filter

library(here)

## here() starts at C:/Users/malin009/OneDrive - Wageningen University & Research/Git/CGM_data

library(dplyr)
library(data.table)

## 
## Attaching package: 'data.table'
## 
## The following objects are masked from 'package:lubridate':
## 
##     hour, isoweek, isoyear, mday, minute, month, quarter, second, wday,
##     week, yday, year
## 
## The following objects are masked from 'package:dplyr':
## 
##     between, first, last
## 
## The following object is masked from 'package:purrr':
## 
##     transpose

library(lubridate)
library(janitor)

## 
## Attaching package: 'janitor'
## 
## The following object is masked from 'package:rstatix':
## 
##     make_clean_names
## 
## The following objects are masked from 'package:stats':
## 
##     chisq.test, fisher.test

library(readxl)
library(ggsci)
library(hms)

## 
## Attaching package: 'hms'
## 
## The following object is masked from 'package:lubridate':
## 
##     hms

library(fuzzyjoin)

## Warning: package 'fuzzyjoin' was built under R version 4.5.3

library(ggrepel)
library(DescTools)

## Warning: package 'DescTools' was built under R version 4.5.3

## 
## Attaching package: 'DescTools'
## 
## The following object is masked from 'package:data.table':
## 
##     %like%

library(lme4)

## Loading required package: Matrix
## 
## Attaching package: 'Matrix'
## 
## The following objects are masked from 'package:tidyr':
## 
##     expand, pack, unpack

library(purrr)
library(broom.mixed)

## Warning: package 'broom.mixed' was built under R version 4.5.3

library(lmerTest)

## Warning: package 'lmerTest' was built under R version 4.5.3

## 
## Attaching package: 'lmerTest'
## 
## The following object is masked from 'package:lme4':
## 
##     lmer
## 
## The following object is masked from 'package:stats':
## 
##     step

library(patchwork)

Setting the work directory

here::i_am("Postprandial_GM.Rmd")

## here() starts at C:/Users/malin009/OneDrive - Wageningen University & Research/Git/CGM_data

script_path <- here::here()
cat("The project root is located at:", script_path, "\n")

## The project root is located at: C:/Users/malin009/OneDrive - Wageningen University & Research/Git/CGM_data

setwd(script_path)

Loading CGM data

Continuous glucose monitoring (CGM) data were imported from multiple raw .csv files. An initial file was inspected to determine relevant columns and file structure, which informed the selection of variables and removal of non-data rows.

Participant identifiers were extracted from file names and appended to each dataset prior to merging into a single dataframe.

Before converting variables to numeric and datetime formats, the structure of the Timestamp column was explored by quantifying the occurrence of specific characters (e.g. dashes, colons, slashes). This step ensured consistent parsing of date-time values.

Glucose values were converted to numeric format (handling comma decimal separators), and timestamps were parsed into standard datetime objects. As historical measurements were included in the raw exports, data prior to 2022 were excluded from further analysis.

Missing CGM observations were excluded to ensure that baseline glucose values reflected the closest valid preprandial measurement.

# Read all files from the raw_data folder (in work directory) and select only columns of interest, remove first rows in each file

files <- list.files(path = "./raw_data", recursive=TRUE,  full.names = TRUE)
my_cols<- c("Timestamp", "glucose mmol/l")
tbl_fread <- map(files, ~fread(.x, skip = "Apparat", select = my_cols))

# the following code adds a column in each df in a list with a name of participant

part_id <- tools::file_path_sans_ext(basename(files))
tbl_fread <- map2(tbl_fread, part_id, ~mutate(.x, part_id = .y))

# join data frames from list
tbl_fread_df <- bind_rows(tbl_fread)


# explore Timestamp column
characters_in_string<- tbl_fread_df %>% mutate(dashes =str_count(Timestamp, pattern="-"),
                                                   colon =str_count(Timestamp, pattern=":"),
                                                   slash =str_count(Timestamp, pattern="/")) %>% select(dashes, colon, slash) %>%
  pivot_longer(cols = c("dashes", "colon", "slash"), names_to = "character", values_to = "number")


ggplot(characters_in_string, aes(number))+
  geom_histogram()+
  facet_wrap(~character)+
  theme_bw()

## `stat_bin()` using `bins = 30`. Pick better value `binwidth`.

# looks like date string has always 1 colon, 2 dashes and 0 slashes, therefore the date looks like: DD/MM/YYYY HH:MM

# changing Timestamp string into date format and glucose measurement into numeric

tbl_fread_df<- tbl_fread_df %>%
  mutate(Timestamp = parse_date_time(Timestamp, "dmy HM"), CGM_glucose = as.numeric(gsub(",", ".", `glucose mmol/l`)))%>%
  filter(year(Timestamp) >= 2022) %>% select(-`glucose mmol/l`)

# deleting timepoint with no glucose measurement

tbl_fread_df <- tbl_fread_df %>%
  drop_na(CGM_glucose)

Loading dietary data

Dietary intake data consisted of two sources:
(1) data exported from the Traqq application, where participants recorded meal timing, and
(2) a study protocol defining the planned food products assigned to each participant.

The Traqq dataset was originally in Dutch. After initial inspection, relevant columns were selected and renamed for clarity. Meal types were derived from session labels (e.g. breakfast vs dinner), based on consultation with colleagues involved in the study setup.

Meal timing information was stored separately as date and time (in seconds), requiring conversion into a unified datetime format to enable integration with CGM data.

Additional variables describing meal completion and quantity consumed were consolidated from separate columns (breakfast vs dinner) into unified fields for downstream analysis.

The code below shows the selection and transformation of variables from traqq application used in further analysis.

app_data<- read.csv("part_portfolio_app.csv")
app_data <- app_data %>% select(-1)

app_data<- app_data %>% select(Subject: Ontbijt.volledig) %>%
  mutate(meal = case_when(Session.Name == "Nuttigen maaltijd?" ~ "Dinner",
                          Session.Name == "Activeiten vandaag?" ~ "Breakfast",
                          TRUE ~ NA_character_)) %>%
  unite(col = "time", c(Hoe.laat.maaltijd, Hoe.laat.ontbijt), na.rm = TRUE) %>%
    mutate(time = as_hms(as.numeric(time)),  # seconds → hms
           Datetime = as.POSIXct(ymd(Date)) + as.numeric(time),
           Date = ymd(Date),
           breakfast_duration = hms::as_hms(Duur.ontbijt),
           completed = coalesce(Maaltijd.volledig, Ontbijt.volledig),
           amount = coalesce(Hoeveel.maaltijd..gram., Hoeveel.ontbijt))%>%
  select(-Session.Name) %>%
  select(Subject, Date, time, Datetime, breakfast_duration, meal, completed, amount) %>%
  drop_na(time)

The code below demonstrates the import of the study protocol containing information on allocated test products, followed by merging with Traqq application data to generate a complete dietary dataset.

products<- read.csv("products_portfolio.csv") %>% select(-1) %>%
  mutate(Date = dmy(Date))
eaten_products<- app_data %>% left_join(products, by = c("Subject", "Date", "meal"))

Data integration

Merging CGM measurements with dietary data

Extracting baseline glucose measurements

To accurately quantify postprandial responses, baseline glucose was defined as the closest CGM measurement preceding each meal. This avoids bias introduced by postprandial glucose increases when no pre-meal measurement is available.

baseline <- tbl_fread_df %>%
  fuzzy_left_join(eaten_products, by = c("part_id" = "Subject", "Timestamp" = "Datetime"),
    match_fun = list(`==`, `<=`)) %>%  # ONLY values before meal
  drop_na(meal) %>% # to ensure only meal moments are considered
  group_by(part_id, Datetime) %>%
  slice_max(Timestamp, n = 1, with_ties = FALSE) %>%  # closest before
  ungroup() %>%
  rename(baseline_glucose = CGM_glucose)

Before baseline measurement will be joined with 2h frame CGM measurement the window time before meal consumption will be calculated to evaluate acceptability of this gap. 30 minute gap will be still accepted regarding that the sensor saves measurements every 30 minutes and sometime participants might have forgotten to scan it with the application..

baseline <- baseline %>%
  mutate(baseline_gap = as.numeric(difftime(Datetime, Timestamp, units = "mins")),
         baseline_quality = case_when(baseline_gap <= 30 ~ "good",
                                 baseline_gap <= 45 ~ "acceptable",
                                 TRUE ~ "poor"))

baseline %>% filter(baseline_quality != "good") %>% group_by(product, baseline_quality) %>%count()

## # A tibble: 4 × 3
## # Groups:   product, baseline_quality [4]
##   product  baseline_quality     n
##   <chr>    <chr>            <int>
## 1 Nasi     acceptable           1
## 2 Pasta    acceptable           1
## 3 Pasta    poor                 1
## 4 WG_Bread poor                 1

baseline %>% filter(baseline_quality != "good") %>% select(part_id, product, Date, baseline_quality, baseline_gap)

## # A tibble: 4 × 5
##   part_id product  Date       baseline_quality baseline_gap
##   <chr>   <chr>    <date>     <chr>                   <dbl>
## 1 P03     Pasta    2022-02-10 acceptable                 31
## 2 P03     Pasta    2022-02-13 poor                       61
## 3 P09     Nasi     2022-02-23 acceptable                 32
## 4 P10     WG_Bread 2022-02-07 poor                       70

In two cases, the temporal gap between meal onset and baseline glucose measurement exceeded one hour, which may compromise baseline estimation. This occurred for the Pasta meal in participant P03 and the WG_Bread meal in participant P10. These time series should therefore be interpreted with caution and may be excluded from downstream analyses if considered biologically unreliable.

Extracting 2 hour window glycemic response

Glucose response curves were extracted over a 2-hour window following each meal. Due to the discrete nature of CGM measurements, exact alignment with meal timing was not always possible; All available measurements within the defined window and after start of meal eating were included.

response <- tbl_fread_df %>%
  fuzzy_inner_join(eaten_products, by = c("part_id" = "Subject",
                                          "Timestamp" = "Datetime"),
    match_fun = list(`==`, function(x, y) x >= y & x <= y + lubridate::hours(2)))

Preparing dataframe with baseline, response measurement and meal names

full_response<- baseline %>% select(-baseline_gap) %>% rename(CGM_glucose = baseline_glucose) %>%
       bind_rows(response) %>% arrange(Timestamp) %>%
  group_by(part_id, Date, product) %>% fill(baseline_quality) %>%
  mutate(time_from_meal = as.numeric(difftime(Timestamp, Datetime, units = "mins"))) %>% ungroup

Quality control

To assess the quality of postprandial glucose response curves, a quality control procedure was conducted at the meal-event level. For each participant and meal occasion, the number of available CGM measurements was calculated, together with the timing of the first and last measurements within the 2-hour postprandial window as well as the coverage of the 2-hour window.

Time of first and last measurement, coverage ratio

qc_response <- full_response %>%
  group_by(part_id, Date, product) %>%
  summarise(first_measurement = min(time_from_meal),
            n_measurements = n(),
            last_measurement = max(time_from_meal),
            coverage_ratio = (last_measurement - first_measurement) / 120,
            .groups = "drop") %>%
  mutate(coverage_description = case_when(coverage_ratio<=0.2 ~ "no",
                                          coverage_ratio<=0.6 ~ "poor",
                                          coverage_ratio<0.8 ~ "to examine",
                                          coverage_ratio>=0.8 ~ "good")) %>%
  mutate(coverage_description = fct_relevel(as.factor(coverage_description), c("no", "poor", "to examine", "good")))

Number of measurements per meal

Exploratory visualisations were used to assess variability in measurement density and temporal coverage across participants and meal types.

qc_response %>% ggplot(aes(x = first_measurement, y = last_measurement, shape = coverage_description, fill = n_measurements)) +
  geom_point(size = 3, alpha = 0.6, aes(x = first_measurement, y = last_measurement, shape = coverage_description, fill = n_measurements)) +
  scale_shape_manual(values = c(25, 24, 23, 21)) +
  geom_text_repel(data=(qc_response %>%filter(n_measurements <3 | coverage_ratio<0.8)), size = 3, aes(label = part_id))+
  scale_fill_gsea() +
  facet_wrap(~product)+
  theme_bw()

Visual inspection of postprandial glucose response curves indicated good data quality for the following meals: Fruit_Yoghurt, Gingerbread_Yoghurt, WG_Bread, and the glucose drink control.

In contrast, several response curves (Pasta, Nasi, R_Bread, and one WG_Bread observation) exhibited incomplete temporal coverage, with the final available CGM measurement frequently occurring approximately 70–80 minutes after meal consumption.

qc_response %>% filter(coverage_ratio<0.8 | n_measurements <6 | first_measurement<=-30 | last_measurement < 60) %>% arrange(coverage_ratio) %>% kbl() %>%
  kableExtra::kable_paper(font_size = 12) %>%
  scroll_box(width = "800px", height = "500px")

part_id	Date	product	first_measurement	n_measurements	last_measurement	coverage_ratio	coverage_description
P03	2022-02-10	Pasta	-31	1	-31	0.0000000	no
P03	2022-02-13	Pasta	-61	1	-61	0.0000000	no
P03	2022-02-14	Pasta	0	2	0	0.0000000	no
P03	2022-02-15	Nasi	-28	1	-28	0.0000000	no
P03	2022-02-16	Nasi	-22	1	-22	0.0000000	no
P06	2022-02-10	Pasta	-5	1	-5	0.0000000	no
P09	2022-02-23	Nasi	-32	1	-32	0.0000000	no
P06	2022-02-14	Nasi	-2	2	13	0.1250000	no
P09	2022-03-01	Nasi	-14	2	1	0.1250000	no
P01	2022-02-09	R_Bread	-11	4	34	0.3750000	poor
P09	2022-02-27	Pasta	-6	4	39	0.3750000	poor
P09	2022-03-02	Nasi	-8	5	52	0.5000000	poor
P06	2022-02-09	R_Bread	-2	6	73	0.6250000	to examine
P09	2022-02-24	Pasta	-4	6	71	0.6250000	to examine
P06	2022-02-16	R_Bread	-3	7	87	0.7500000	to examine
P10	2022-02-17	R_Bread	-6	7	84	0.7500000	to examine
P09	2022-02-25	glucose_drink	-2	5	108	0.9166667	good
P01	2022-02-08	Pasta	-26	5	107	1.1083333	good
P10	2022-02-07	WG_Bread	-70	4	111	1.5083333	good

Based on this quality assessment, response curves with insufficient postprandial coverage were excluded from downstream analyses. Additionally, postprandial response curves with baseline glucose measurements obtained more than 1 hour before meal consumption were excluded from further analysis. This decision was made to ensure comparability of iAUC estimates across meals and participants.

Generating a list of measurements for inclusion in data analysis

The code below generates a vector of participant–date–product combinations that were removed from the full response dataset due to poor temporal coverage.

records_to_exclude<- qc_response %>% filter(coverage_ratio<0.8 | first_measurement < -60) %>%
 arrange(coverage_ratio) %>%
  mutate(exclude = paste0(part_id, "_", Date, "_", product)) %>% pull(exclude)

Preparing final dataframe with good quality data

Below, the final analytical dataset is prepared following quality control and exclusion procedures. Data quality was reassessed through additional visual inspection of postprandial glucose response curves.

Subsequently, repeated consumption occasions of the same test products were identified by assigning consecutive meal numbers within each participant, as most products (except the glucose drink control) were consumed multiple times.

full_response_qc<- full_response %>% mutate(exclude = paste0(part_id, "_", Date, "_", product)) %>%
  filter(!exclude %in% records_to_exclude) %>% select(-exclude)

full_response_qc %>% group_by(part_id, Date, product) %>%
  summarise(first_measurement = min(time_from_meal),
            n_measurements = n(),
            last_measurement = max(time_from_meal),
            coverage_ratio = (last_measurement - first_measurement) / 120,
            .groups = "drop") %>%
  mutate(coverage_description = case_when(coverage_ratio<=0.2 ~ "no",
                                          coverage_ratio<=0.6 ~ "poor",
                                          coverage_ratio<0.8 ~ "to examine",
                                          coverage_ratio>=0.8 ~ "good")) %>%
  mutate(coverage_description = fct_relevel(as.factor(coverage_description), c("no", "poor", "to examine", "good"))) %>%
  ggplot(aes(x = first_measurement, y = last_measurement, shape = coverage_description, fill = n_measurements)) +
  geom_point(size = 3, alpha = 0.6, aes(x = first_measurement, y = last_measurement, shape = coverage_description, fill = n_measurements)) +
  scale_shape_manual(values = c(25, 24, 23, 21)) +
  scale_fill_gsea() +
  facet_wrap(~product)+
  theme_bw()

## Warning: There was 1 warning in `mutate()`.
## ℹ In argument: `coverage_description = fct_relevel(...)`.
## Caused by warning:
## ! 3 unknown levels in `f`: no, poor, and to examine

full_response_qc<- full_response_qc %>%
       arrange(part_id, product, Date) %>%
       group_by(part_id, product) %>% mutate(consumption = dense_rank(Date)) %>%
  ungroup %>%
  unique %>%
  mutate(completed = if_else(completed == "1", "yes", "no")) %>%
  mutate(completed = fct_relevel(as.factor(completed), c("yes", "no"))) %>%
  mutate(consumption = case_when(consumption == "1" ~ "first",
                                 consumption == "2" ~ "second",
                                 consumption == "3" ~ "third",
                                 consumption == "4" ~ "fourth",
                                 consumption == "5" ~ "fifth"))

Downstream analysis

Exploration of glycemic response

Below, several postprandial glycemic response curves are presented alongside the reference response to the glucose drink. The first panel shows breakfasts with a caloric content of approximately 280 kcal, the second panel includes breakfasts ranging from 350–390 kcal, and the third panel displays all dinner meals.

# plot caloric content of approximately 280 kcal
full_response_qc %>% filter(product %in% c("Fruit_Yoghurt", "Gingerbread_Yoghurt", "glucose_drink")) %>%
  mutate(product = fct_relevel(as.factor(product), c("Fruit_Yoghurt", "Gingerbread_Yoghurt", "glucose_drink"))) %>%
  ggplot(aes(x=time_from_meal, y = CGM_glucose, group = interaction(part_id, product, consumption)))+
  geom_point(size = 3, color = "black", aes(shape = completed, fill = product))+
  scale_shape_manual(values = c(21, 24), name = "completed")+
  geom_smooth(se = FALSE, aes(color = product, linetype = consumption))+
  scale_color_manual(values = pal_bmj("default")(9)[c(1,2,9)])+
  scale_fill_manual(values = pal_bmj("default")(9)[c(1,2,9)])+
  facet_wrap(~part_id, scales = "free")+
  labs(x = "time from consumption [min]", y = "glucose concentration [mmol/l]", title = "Individual glycemic response to breakfast meal", subtitle = "fruit with yoghurt, gingerbread with yoghurt and standard glucose drink")+
  guides(color = guide_legend(ncol = 1, override.aes = list(shape = 21)),
         linetype = guide_legend(ncol = 1, override.aes = list(color = "black")),
         shape = guide_legend(ncol = 1)) +
  theme_bw()+
  theme(legend.box = "horizontal", legend.position = c(0.75, 0.2))

## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'

# plot caloric content of 350-390 kcal
full_response_qc %>% filter(product %in% c("WG_Bread", "R_Bread", "glucose_drink")) %>%
  mutate(product = fct_relevel(as.factor(product), c("WG_Bread", "R_Bread", "glucose_drink")))%>%
  ggplot(aes(x=time_from_meal, y = CGM_glucose, group = interaction(part_id, product, consumption)))+
  geom_point(size = 3, color = "black", aes(shape = completed, fill = product))+
  scale_shape_manual(values = c(21, 24), name = "completed")+
  geom_smooth(se = FALSE, aes(color = product, linetype = consumption))+
  scale_color_manual(values = pal_bmj("default")(9)[c(4,5,9)])+
  scale_fill_manual(values = pal_bmj("default")(9)[c(4,5,9)])+
  facet_wrap(~part_id, scales = "free")+
  labs(x = "time from consumption [min]", y = "glucose concentration [mmol/l]", title = "Individual glycemic response to breakfast meal", subtitle = "whole grain with cheese spread, refined bread with cheeses spread and standard glucose drink")+
  guides(color = guide_legend(ncol = 1, override.aes = list(shape = 21)),
         linetype = guide_legend(ncol = 1, override.aes = list(color = "black")),
         shape = guide_legend(ncol = 1)) +
  theme_bw()+
  theme(legend.box = "horizontal", legend.position = c(0.75, 0.2))

## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'

# plot dinner
full_response_qc %>% filter(product %in% c("Pasta", "Nasi", "glucose_drink")) %>%
  mutate(product = fct_relevel(as.factor(product), c("Pasta", "Nasi", "glucose_drink")))%>%
  mutate(consumption = fct_relevel(as.factor(consumption), c("first", "second", "third", "fourth", "fifth")))%>%
  ggplot(aes(x=time_from_meal, y = CGM_glucose, group = interaction(part_id, product, consumption)))+
  geom_point(size = 2, color = "black", aes(shape = completed, fill = product))+
  scale_shape_manual(values = c(21, 24), name = "completed")+
  geom_smooth(se = FALSE, aes(color = product, linetype = consumption))+
  scale_color_manual(values = pal_bmj("default")(9)[c(7,8,9)])+
  scale_fill_manual(values = pal_bmj("default")(9)[c(7,8,9)])+
  scale_linetype_manual(values = c("solid", "longdash", "dashed", "dotdash", "dotted"))+
  facet_wrap(~part_id, scales = "free")+
  labs(x = "time from consumption [min]", y = "glucose concentration [mmol/l]", title = "Individual glycemic response to dinner", subtitle = "nasi, pasta and standard glucose drink")+
   guides(color = guide_legend(ncol = 1, override.aes = list(shape = 21)),
         linetype = guide_legend(ncol = 1, keywidth = unit(2, "cm"), override.aes = list(color = "black")),
         shape = guide_legend(ncol = 1)) +
  theme_bw()+
  theme(legend.box = "horizontal", legend.position = c(0.75, 0.17))

## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'

Evaluation of iAUC

Calculation of iAUC

To quantify the postprandial response to each meal and consumption occasion, the incremental area under the curve (iAUC) was calculated. This approach was chosen because it adjusts for baseline glucose levels and focuses on glucose excursions above baseline, without penalizing decreases below baseline. Postprandial response curves were aligned to meal start time, and only measurements collected after meal consumption contributed to iAUC calculation.

Since the DescTools package does not provide a direct function for calculating incremental area under the curve (iAUC), a custom implementation was used. Glucose values were first adjusted relative to baseline glucose concentration. During iAUC calculation, only glucose values above baseline contributed to the final area, whereas values below baseline were set to zero before applying trapezoidal integration. This approach ensures that decreases in glucose concentration below baseline do not artificially reduce the estimated postprandial glycemic response.

The code below shows the iAUC calculation. To account for variability in observation windows, all curves were standardized to a fixed 2-hour (120 minutes) postprandial period. Although the target window was 2 hours, actual coverage varied slightly depending on exact meal timing and sensor scanning frequency. Standardization ensures comparability of iAUC values across participants and consumption occasions.

full_response_qc<- full_response_qc %>%
  group_by(part_id, product, consumption) %>%
  arrange(Timestamp, .by_group = TRUE) %>%
  mutate(time_from_baseline = as.numeric(difftime(Timestamp, first(Timestamp), units = "mins")),
         glucose_adj = CGM_glucose - first(CGM_glucose),
         glucose_increment = pmax(glucose_adj, 0)) %>%
  ungroup


parameters_2h<- full_response_qc %>%  group_by(part_id, product, consumption) %>%
       arrange(Timestamp, .by_group = TRUE) %>%
       mutate(base_glc = first(CGM_glucose),
            peak_glc = max(CGM_glucose),
            min_glc = min(CGM_glucose)) %>% ungroup %>%
       arrange(Timestamp) %>%
       filter(time_from_meal >=0) %>%
       group_by(part_id, product, consumption, base_glc, peak_glc, min_glc) %>%
       arrange(Timestamp, .by_group = TRUE) %>%
       summarize(duration=last(time_from_meal) - first(time_from_meal),
     iAUC = AUC(x = time_from_meal,
        y = glucose_increment,
        method = "trapezoid"), .groups = "drop") %>%
  mutate(iAUC_norm = iAUC*120/duration)

Visualisations of calculated iAUC

The visualization below illustrates the approach used for incremental area under the curve (iAUC) calculation using the WG_Bread response as an example. Each panel represents one participant and shows the postprandial glucose response measured using continuous glucose monitoring (CGM).

The green vertical line indicates the time of meal consumption, whereas the red dashed horizontal line represents baseline glucose concentration determined from the closest CGM measurement prior to meal intake. Blue dashed vertical lines indicate the beginning and end of the response window used for iAUC calculation.

The shaded blue area represents glucose excursions above baseline that contributed to the final iAUC value. Glucose concentrations below baseline were excluded from the calculation by setting negative baseline-adjusted values to zero prior to trapezoidal integration.

response_WG_Bread_iAUC_calc<- full_response_qc %>%
  group_by(part_id, product, consumption) %>%
       arrange(Timestamp, .by_group = TRUE) %>%
       mutate(base_glc = first(CGM_glucose),
            peak_glc = max(CGM_glucose),
            min_glc = min(CGM_glucose)) %>% ungroup %>%
  filter(product == "WG_Bread" & consumption == "first")

i_response_WG_Bread_iAUC_calc <- response_WG_Bread_iAUC_calc%>%
  group_by(part_id, product, consumption) %>%
       arrange(Timestamp, .by_group = TRUE) %>%
                filter(time_from_meal>=0) %>%
  mutate(first_meas = first(time_from_meal),
         last_meas = last(time_from_meal),
         glucose_iAUC = pmax(CGM_glucose, base_glc)) %>% ungroup

response_WG_Bread_iAUC_calc %>%
  ggplot(aes(x=time_from_meal, y = CGM_glucose, group = part_id, product))+
  geom_point(size = 3, color = "black", fill = "#7D5CC6FF")+
  geom_line(color = "#7D5CC6FF")+
  geom_vline(aes(xintercept=0), linetype="solid", size = 1, color = "green4")+
  geom_vline(data = i_response_WG_Bread_iAUC_calc, aes(xintercept=first_meas), linetype="dashed", linewidth = 1, color = "dodgerblue1")+
  geom_vline(data = i_response_WG_Bread_iAUC_calc, aes(xintercept=last_meas), linetype="dashed", linewidth = 1, color = "dodgerblue1")+
  geom_hline(aes(yintercept=base_glc), linetype="dashed", linewidth = 1, color = "red3")+
  geom_ribbon(data = i_response_WG_Bread_iAUC_calc,
              aes(ymin = base_glc, ymax = glucose_iAUC),
              fill = "dodgerblue1",alpha = 0.3)+
  facet_wrap(~part_id, scales = "free")+
  labs(x = "time from consumption [min]", y = "glucose concentration [mmol/l]", title = "iAUC calculation", subtitle = "WG_Bread shown as an example.\nGreen line = meal consumption time; blue dashed lines = iAUC calculation window;\nred dashed line = baseline glucose concentration.")+
  theme_bw()+
  theme(legend.box = "horizontal", legend.position = c(0.75, 0.2))

## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once per session.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

Comparison of calculated parameters between consumption occassions

Normality of data distribution

First, it was assessed whether the same meal elicits different glycemic responses across repeated consumption occasions, while accounting for participant-specific baseline differences using linear mixed-effects models with participant-level random effects.

Prior to model fitting, distributional assumptions were evaluated. Normality was assessed using visual inspection of Q–Q plots, supported by Shapiro–Wilk tests. Given the sensitivity of formal tests to sample size, graphical diagnostics were prioritized for inference.

Due to fundamentally different statistical properties, integrated postprandial response (iAUC) was analysed separately from instantaneous glucose metrics (baseline, peak, and minimum glucose levels). This ensured appropriate scaling and interpretability of distributional diagnostics.

plot_data_consumption <- parameters_2h %>% select(-c(duration, iAUC)) %>%
   filter(product != "glucose_drink") %>%
  pivot_longer(cols = c("base_glc", "peak_glc", "min_glc", "iAUC_norm"), values_to = "value", names_to = "parameter") 

stat_consumption_df <- plot_data_consumption %>%
  group_by(parameter, product) %>%
  summarise(p = shapiro.test(value)$p.value,
    label = paste0("Shapiro p = ", signif(p, 2)),
    .groups = "drop")

plot_data_consumption %>%
  filter(parameter == "iAUC_norm") %>%
  ggplot(aes(value)) +
  geom_histogram(aes(y = after_stat(density)), bins = 25, fill = "dodgerblue4", color = "black", alpha = 0.8) +
  geom_density(color = "red", linewidth = 1) +
  geom_text(data = (stat_consumption_df %>% filter(parameter == "iAUC_norm")),
    aes(x = -Inf, y = Inf, label = label), hjust = -0.1, vjust = 1.2, inherit.aes = FALSE, size = 4) +
  facet_wrap(~product) +
  theme_bw() +
  labs(title = "Distribution of iAUC across meals", x = "iAUC (normalized)", y = "Density")

plot_data_consumption %>%
  filter(parameter != "iAUC_norm") %>%
  ggplot(aes(value)) +
  geom_histogram(aes(y = after_stat(density)), bins = 25, fill = "dodgerblue2", color = "black", alpha = 0.8) +
  geom_density(color = "red") +
  geom_text(data = (stat_consumption_df %>% filter(parameter != "iAUC_norm")),
    aes(x = -Inf, y = Inf, label = label), hjust = -0.1, vjust = 1.2, inherit.aes = FALSE, size = 4) +
  facet_grid(cols = vars(product), rows = vars(parameter), scales = "free") +
  theme_bw()

ggplot(plot_data_consumption, aes(sample = value)) +
  stat_qq(size = 2, alpha = 0.7) +
  stat_qq_line(color = "red") +
  facet_grid(cols = vars(product), rows = vars(parameter), scales = "free") +
  theme_bw() +
  labs(title = "Q-Q plots for glycemic response parameters")

Histogram and Q–Q plots indicated that all parameters, apart from iAUC, were approximately normally distributed. Therefore, iAUC will be log and square root-transformed to reduce skewness, and distributional assumptions will be reassessed after transformation.

stat_consumption_df_iAUC_trans<- parameters_2h %>%
  mutate(iAUC_log = log1p(iAUC_norm),
         iAUC_sqrt = sqrt(iAUC_norm)) %>%
  pivot_longer(cols = c("iAUC_log", "iAUC_sqrt"), values_to = "value", names_to = "iAUC_trans") %>%
  group_by(product, iAUC_trans) %>%
  summarise(p = shapiro.test(value)$p.value,
    label = paste0("Shapiro p = ", signif(p, 2)),
    .groups = "drop")

For future analyses, the transformation strategy will be tailored to the distributional characteristics of each meal type. A square-root transformation of iAUC will be applied for dinner meals and bread-based breakfasts, whereas a log transformation will be used for glucose drink and yoghurt-based breakfasts, due to their more pronounced right-skewed distributions.

Comparison between consumption moments

Statistics

parameters_2h_transformed<- parameters_2h %>% select(-c(duration, iAUC)) %>% mutate(iAUC_log = log1p(iAUC_norm),
         iAUC_sqrt = sqrt(iAUC_norm)) %>%
  mutate(iAUC_transformed = case_when(
      product %in% c("Pasta", "Nasi", "WG_Bread", "R_Bread") ~ iAUC_sqrt,
      product %in% c("Fruit_Yoghurt", "Gingerbread_Yoghurt", "glucose_drink") ~ iAUC_log,
      TRUE ~ NA_real_))

consumption_models<- parameters_2h_transformed %>% 
  filter(product != "glucose_drink") %>%
  mutate(consumption = fct_relevel(as.factor(consumption), c("first", "second", "third", "fourth", "fifth"))) %>%
  pivot_longer(cols = c("base_glc", "peak_glc", "min_glc", "iAUC_transformed"), values_to = "value", names_to = "parameter") %>% group_by(parameter, product) %>%
   nest() %>%
  mutate(model = map(data, ~ lmer(value ~ consumption + (1|part_id), data = .x)),
     tidy_results = map(model,~ broom.mixed::tidy(.x, effects = "fixed")))%>%
  select(parameter, product, tidy_results) %>%
  unnest(tidy_results)

consumption_models<- consumption_models %>% filter(term != "(Intercept)")  %>% mutate(
    label = paste0("p = ", round(p.value, 3)),
    text_color = ifelse(p.value < 0.05, "red", "black"))
  
consumption_models %>% arrange(p.value) %>% kbl() %>%
  kableExtra::kable_paper(font_size = 12) %>%
  scroll_box(width = "800px", height = "500px")

parameter	product	effect	term	estimate	std.error	statistic	df	p.value	label	text_color
iAUC_transformed	Nasi	fixed	consumptionthird	-6.1513251	2.6046365	-2.3616828	21.392847	0.0277403	p = 0.028	red
base_glc	Pasta	fixed	consumptionfourth	1.2586786	0.5844212	2.1537182	28.202484	0.0399541	p = 0.04	red
base_glc	Nasi	fixed	consumptionthird	1.2522227	0.5777195	2.1675271	20.901264	0.0418929	p = 0.042	red
base_glc	WG_Bread	fixed	consumptionsecond	0.2006307	0.0977390	2.0527185	7.130190	0.0784781	p = 0.078	black
iAUC_transformed	R_Bread	fixed	consumptionsecond	-1.6340990	0.8795486	-1.8578837	5.753646	0.1146482	p = 0.115	black
peak_glc	Pasta	fixed	consumptionfourth	0.6873727	0.4248813	1.6177995	27.814802	0.1169918	p = 0.117	black
peak_glc	Pasta	fixed	consumptionfifth	0.6540393	0.4248813	1.5393462	27.814802	0.1350200	p = 0.135	black
base_glc	Nasi	fixed	consumptionsecond	0.8118857	0.5289063	1.5350276	20.252687	0.1402523	p = 0.14	black
iAUC_transformed	WG_Bread	fixed	consumptionsecond	-2.1820941	1.4827926	-1.4716111	8.179538	0.1785236	p = 0.179	black
base_glc	Fruit_Yoghurt	fixed	consumptionsecond	-0.2300000	0.1605892	-1.4322259	9.000000	0.1858809	p = 0.186	black
min_glc	Nasi	fixed	consumptionthird	0.4306475	0.3156705	1.3642312	20.110048	0.1875690	p = 0.188	black
min_glc	Fruit_Yoghurt	fixed	consumptionsecond	-0.2300000	0.1626516	-1.4140651	9.000000	0.1909894	p = 0.191	black
peak_glc	Pasta	fixed	consumptionthird	0.4790661	0.3680378	1.3016764	27.455091	0.2038437	p = 0.204	black
iAUC_transformed	Fruit_Yoghurt	fixed	consumptionsecond	0.2578203	0.1882431	1.3696136	9.000000	0.2040061	p = 0.204	black
iAUC_transformed	Pasta	fixed	consumptionthird	1.9914792	1.9061546	1.0447627	28.405791	0.3049425	p = 0.305	black
peak_glc	Nasi	fixed	consumptionsecond	0.4297878	0.4123316	1.0423354	19.670490	0.3098983	p = 0.31	black
base_glc	Nasi	fixed	consumptionfourth	0.5947847	0.6098738	0.9752586	21.060445	0.3405018	p = 0.341	black
peak_glc	Pasta	fixed	consumptionsecond	0.3400000	0.3543977	0.9593741	27.234324	0.3458107	p = 0.346	black
iAUC_transformed	Pasta	fixed	consumptionfourth	-2.0801632	2.1719616	-0.9577348	29.929874	0.3458695	p = 0.346	black
peak_glc	WG_Bread	fixed	consumptionsecond	-0.3551710	0.3584243	-0.9909233	7.382764	0.3530828	p = 0.353	black
min_glc	Pasta	fixed	consumptionthird	0.2723426	0.3063728	0.8889256	27.514270	0.3817491	p = 0.382	black
iAUC_transformed	Gingerbread_Yoghurt	fixed	consumptionsecond	-0.1124059	0.1446809	-0.7769230	9.000000	0.4571374	p = 0.457	black
iAUC_transformed	Nasi	fixed	consumptionfourth	-1.8476534	2.7395099	-0.6744467	21.941505	0.5070719	p = 0.507	black
peak_glc	Nasi	fixed	consumptionthird	-0.3055662	0.4528623	-0.6747442	19.927186	0.5075919	p = 0.508	black
peak_glc	R_Bread	fixed	consumptionsecond	-0.3506697	0.5107990	-0.6865122	5.872951	0.5185609	p = 0.519	black
peak_glc	Fruit_Yoghurt	fixed	consumptionsecond	0.1300000	0.1961009	0.6629241	9.000000	0.5239924	p = 0.524	black
base_glc	Pasta	fixed	consumptionfifth	0.3253453	0.5844212	0.5566965	28.202484	0.5821259	p = 0.582	black
iAUC_transformed	Nasi	fixed	consumptionsecond	-1.3021006	2.4107747	-0.5401171	20.123794	0.5950438	p = 0.595	black
peak_glc	Nasi	fixed	consumptionfourth	-0.2277007	0.4785869	-0.4757770	19.965160	0.6393978	p = 0.639	black
peak_glc	Gingerbread_Yoghurt	fixed	consumptionsecond	-0.1900000	0.3973384	-0.4781819	9.000000	0.6439221	p = 0.644	black
min_glc	Nasi	fixed	consumptionfourth	0.1543236	0.3335580	0.4626591	20.162940	0.6485629	p = 0.649	black
iAUC_transformed	Pasta	fixed	consumptionsecond	0.7548602	1.8471750	0.4086566	27.807755	0.6859218	p = 0.686	black
iAUC_transformed	Pasta	fixed	consumptionfifth	0.8388478	2.1719616	0.3862167	29.929874	0.7020682	p = 0.702	black
min_glc	WG_Bread	fixed	consumptionsecond	0.0285707	0.0772040	0.3700684	7.050758	0.7221977	p = 0.722	black
min_glc	Nasi	fixed	consumptionsecond	0.1021737	0.2876614	0.3551875	19.788989	0.7262068	p = 0.726	black
base_glc	Pasta	fixed	consumptionsecond	0.1600000	0.4889668	0.3272206	27.332949	0.7459940	p = 0.746	black
min_glc	Pasta	fixed	consumptionfourth	-0.0878739	0.3534534	-0.2486153	27.943403	0.8054765	p = 0.805	black
base_glc	R_Bread	fixed	consumptionsecond	-0.0393056	0.1602149	-0.2453305	6.177904	0.8141477	p = 0.814	black
base_glc	Pasta	fixed	consumptionthird	0.1028173	0.5072099	0.2027115	27.647978	0.8408449	p = 0.841	black
min_glc	Pasta	fixed	consumptionfifth	0.0454594	0.3534534	0.1286151	27.943403	0.8985835	p = 0.899	black
min_glc	Pasta	fixed	consumptionsecond	0.0300000	0.2951385	0.1016472	27.258084	0.9197808	p = 0.92	black
base_glc	Gingerbread_Yoghurt	fixed	consumptionsecond	-0.0100000	0.1779201	-0.0562050	9.000000	0.9564065	p = 0.956	black
min_glc	R_Bread	fixed	consumptionsecond	-0.0018960	0.1211812	-0.0156463	6.062898	0.9880187	p = 0.988	black
min_glc	Gingerbread_Yoghurt	fixed	consumptionsecond	0.0000000	0.1282359	0.0000000	9.000000	1.0000000	p = 1	black

The results indicate that only baseline glucose prior to the third occasion of Nasi consumption, the fourth occasion of Pasta consumption, and iAUC during the third Nasi consumption differed significantly from the first consumption occasion.

Visualisation

Below also a plot with such visualization is presented.

# Since participant P06 had no glycemic response data for R_Bread consumption, participant-specific colors were fixed across plots to maintain consistent color mapping and enable a common legend.

all_participants <- c("P01","P02","P03","P04","P05",
                      "P06","P07","P08","P09","P10")

plot_data_consumption <- parameters_2h_transformed %>%
  mutate(consumption = fct_relevel(as.factor(consumption),
                                   c("first","second","third","fourth","fifth"))) %>%
  mutate(part_id = factor(part_id, levels = all_participants)) %>%
  mutate(product = fct_relevel(as.factor(product), c("Fruit_Yoghurt","Gingerbread_Yoghurt","WG_Bread","R_Bread","Nasi", "Pasta", "glucose_drink"))) %>%
  mutate(meal_set = case_when(product %in% c("Fruit_Yoghurt","Gingerbread_Yoghurt") ~ "Breakfast_280kcal",
                              product %in% c("WG_Bread","R_Bread") ~ "Breakfast_380kcal",
                              product %in% c("Nasi", "Pasta") ~ "Dinner",
                              TRUE ~ NA_character_)) %>%
  pivot_longer(cols = c(base_glc, peak_glc, min_glc, iAUC_transformed),
    names_to = "parameter",
    values_to = "value")

participant_ids <- sort(unique(parameters_2h_transformed$part_id))

participant_colors <- setNames(
  glasbey(length(participant_ids)),
  participant_ids)


breakfast_set1_consumption_plot<- plot_data_consumption %>%
    filter(meal_set == "Breakfast_280kcal") %>%
    ggplot(aes(x = consumption, y = value)) +
    geom_boxplot(fill = "grey90", outlier.shape = NA, width = 0.5) +
    geom_jitter(aes(fill = part_id),
                color = "black",
                shape = 21,
                width = 0.2,
                size = 2.5,
                alpha = 0.8) +
    geom_line(aes(group = part_id, color = part_id), alpha = 0.4) +
    stat_summary(fun = mean, geom = "point", shape = 23, fill = "dodgerblue1", size = 2) +
    stat_summary(fun = mean, geom = "text",
                 color = "dodgerblue3", vjust = -2,
                 aes(label = round(after_stat(y), 1))) +
  geom_text(data = (consumption_models %>% filter(term != "(Intercept)") %>%
                      filter(product %in% c("Fruit_Yoghurt","Gingerbread_Yoghurt"))),
    aes(x = -Inf, y = Inf, label = paste0("p = ", round(p.value, 3))), hjust = -0.1, vjust = 1.2, inherit.aes = FALSE, size = 4) +  
  scale_fill_manual(values = participant_colors, drop = FALSE) +
    scale_color_manual(values = participant_colors, drop = FALSE) +
    facet_grid(rows = vars(parameter), cols = vars(product), scales = "free")+
    theme_bw() +
    theme(axis.text.x = element_text(angle = 45, hjust = 1), axis.title.x = element_blank(), axis.title.y = element_blank())

breakfast_set2_consumption_plot<- plot_data_consumption %>%
    filter(meal_set == "Breakfast_380kcal") %>%
    ggplot(aes(x = consumption, y = value)) +
    geom_boxplot(fill = "grey90", outlier.shape = NA, width = 0.5) +
    geom_jitter(aes(fill = part_id),
                color = "black",
                shape = 21,
                width = 0.2,
                size = 2.5,
                alpha = 0.8) +
    geom_line(aes(group = part_id, color = part_id), alpha = 0.4) +
    stat_summary(fun = mean, geom = "point", shape = 23, fill = "dodgerblue1", size = 2) +
    stat_summary(fun = mean, geom = "text",
                 color = "dodgerblue3", vjust = -2,
                 aes(label = round(after_stat(y), 1))) +
  geom_text(data = (consumption_models %>% filter(term != "(Intercept)") %>%
                      filter(product %in% c("WG_Bread","R_Bread"))),
    aes(x = -Inf, y = Inf, label = paste0("p = ", round(p.value, 3))), hjust = -0.1, vjust = 1.2, inherit.aes = FALSE, size = 4) +  
  scale_fill_manual(values = participant_colors, drop = FALSE) +
    scale_color_manual(values = participant_colors, drop = FALSE) +
    facet_grid(rows = vars(parameter), cols = vars(product), scales = "free")+
    theme_bw() +
    theme(axis.text.x = element_text(angle = 45, hjust = 1), axis.title.x = element_blank(), axis.title.y = element_blank())

dinner_consumption_plot<- plot_data_consumption %>%
    filter(meal_set == "Dinner") %>%
    ggplot(aes(x = consumption, y = value)) +
    geom_boxplot(fill = "grey90", outlier.shape = NA, width = 0.5) +
    geom_jitter(aes(fill = part_id),
                color = "black",
                shape = 21,
                width = 0.2,
                size = 2.5,
                alpha = 0.8) +
    geom_line(aes(group = part_id, color = part_id), alpha = 0.4) +
    stat_summary(fun = mean, geom = "point", shape = 23, fill = "dodgerblue1", size = 2) +
    stat_summary(fun = mean, geom = "text",
                 color = "dodgerblue3", vjust = -2,
                 aes(label = round(after_stat(y), 1))) +
  geom_text(data = (consumption_models %>%
                      filter(product %in% c("Nasi","Pasta")) %>%
                      filter(term == "consumptionsecond")),
    aes(x = -Inf, y = Inf, label = paste0("p_second = ", round(p.value, 3))), hjust = -0.1, vjust = 1.2, inherit.aes = FALSE, size = 4, color = consumption_models %>%
    filter(product %in% c("Nasi","Pasta")) %>%
    filter(term == "consumptionsecond") %>%
    pull(text_color)) +  
  geom_text(data = (consumption_models %>%
                      filter(product %in% c("Nasi","Pasta")) %>%
                      filter(term == "consumptionthird")),
    aes(x = -Inf, y = Inf, label = paste0("p_third = ", round(p.value, 3))), hjust = -0.1, vjust = 2.5, inherit.aes = FALSE, size = 4, color = consumption_models %>%
    filter(product %in% c("Nasi","Pasta")) %>%
    filter(term == "consumptionthird") %>%
    pull(text_color)) + 
  geom_text(data = (consumption_models %>% 
                      filter(product %in% c("Nasi","Pasta")) %>%
                      filter(term == "consumptionfourth")),
    aes(x = -Inf, y = Inf, label = paste0("p_fourth = ", round(p.value, 3))), hjust = -0.1, vjust = 3.7, inherit.aes = FALSE, size = 4, color = consumption_models %>%
    filter(product %in% c("Nasi","Pasta")) %>%
    filter(term == "consumptionfourth") %>%
    pull(text_color)) +  
  scale_fill_manual(values = participant_colors, drop = FALSE) +
    scale_color_manual(values = participant_colors, drop = FALSE) +
    facet_grid(rows = vars(parameter), cols = vars(product), scales = "free")+
    theme_bw() +
    theme(axis.text.x = element_text(angle = 45, hjust = 1), axis.title.x = element_blank(), axis.title.y = element_blank())


print(breakfast_set1_consumption_plot)

print(breakfast_set2_consumption_plot)

print(dinner_consumption_plot)

Comparison between different products

This section compares glycemic response parameters between products with similar caloric value and meal type. Products were grouped into three meal sets: lower-calorie breakfasts (~280 kcal), bread-based breakfasts (~350–390 kcal), and dinner meals. For each parameter (baseline glucose, peak glucose, minimum glucose, and transformed iAUC), linear mixed-effects models were fitted with product as a fixed effect and participant ID as a random effect to account for repeated measurements within individuals.

Model results were extracted and p-values corresponding to product differences were visualized directly on the plots. Boxplots together with individual observations were used to present the distribution of glycemic response parameters across meal sets. Mean values for each product are additionally displayed on the plots. Significant product differences (p < 0.05) are highlighted in red.

product_models<- parameters_2h_transformed %>%
  mutate(meal_set = case_when(product %in% c("Fruit_Yoghurt","Gingerbread_Yoghurt") ~ "Breakfast_280kcal",
                              product %in% c("WG_Bread","R_Bread") ~ "Breakfast_380kcal",
                              product %in% c("Nasi", "Pasta") ~ "Dinner",
                              TRUE ~ NA_character_)) %>%
  filter(product != "glucose_drink") %>%
  mutate(consumption = fct_relevel(as.factor(consumption), c("first", "second", "third", "fourth", "fifth"))) %>%
  pivot_longer(cols = c("base_glc", "peak_glc", "min_glc", "iAUC_transformed"), values_to = "value", names_to = "parameter") %>%
  group_by(parameter, meal_set) %>%
   nest() %>%
  mutate(model = map(data, ~ lmer(value ~ product + (1|part_id), data = .x)),
     tidy_results = map(model,~ broom.mixed::tidy(.x, effects = "fixed")))%>%
  select(parameter, meal_set, tidy_results) %>%
  unnest(tidy_results)

product_models<- product_models %>% filter(term != "(Intercept)")  %>% mutate(
    label = paste0("p = ", round(p.value, 3)),
    text_color = ifelse(p.value < 0.05, "red", "black"))
  
product_models %>% arrange(p.value) %>% kbl() %>%
  kableExtra::kable_paper(font_size = 12) %>%
  scroll_box(width = "800px", height = "500px")

parameter	meal_set	effect	term	estimate	std.error	statistic	df	p.value	label	text_color
peak_glc	Breakfast_280kcal	fixed	productGingerbread_Yoghurt	1.2700000	0.2127894	5.9683436	29.00000	0.0000017	p = 0	red
iAUC_transformed	Breakfast_280kcal	fixed	productGingerbread_Yoghurt	0.5415323	0.1150014	4.7089185	29.00000	0.0000569	p = 0	red
peak_glc	Dinner	fixed	productPasta	-0.5081198	0.1888682	-2.6903401	62.19093	0.0091542	p = 0.009	red
base_glc	Breakfast_380kcal	fixed	productWG_Bread	0.1848718	0.1046232	1.7670258	23.24711	0.0903515	p = 0.09	black
iAUC_transformed	Dinner	fixed	productPasta	-1.9002372	1.1669401	-1.6283931	63.44545	0.1083977	p = 0.108	black
iAUC_transformed	Breakfast_380kcal	fixed	productWG_Bread	-1.2356188	0.7909588	-1.5621784	23.49331	0.1316171	p = 0.132	black
base_glc	Breakfast_280kcal	fixed	productGingerbread_Yoghurt	0.1300000	0.1105004	1.1764660	29.00000	0.2489757	p = 0.249	black
min_glc	Breakfast_380kcal	fixed	productWG_Bread	0.0754855	0.0735237	1.0266827	23.09903	0.3152018	p = 0.315	black
peak_glc	Breakfast_380kcal	fixed	productWG_Bread	-0.2370783	0.2595725	-0.9133412	23.26611	0.3704313	p = 0.37	black
min_glc	Dinner	fixed	productPasta	-0.0589175	0.1435923	-0.4103108	62.21280	0.6829876	p = 0.683	black
min_glc	Breakfast_280kcal	fixed	productGingerbread_Yoghurt	0.0150000	0.1024737	0.1463790	29.00000	0.8846350	p = 0.885	black
base_glc	Dinner	fixed	productPasta	0.0357259	0.2612722	0.1367383	62.25436	0.8916789	p = 0.892	black

plot_data_consumption %>%
    filter(product != "glucose_drink") %>%
    ggplot(aes(x = meal_set, y = value)) +
    geom_boxplot(aes(fill = product), outlier.shape = NA, width = 0.5, alpha = 0.5) +
    geom_jitter(aes(fill = product),
                color = "black", shape = 21, size = 2.5, alpha = 0.8, position = position_jitterdodge(jitter.width = 0.2,
    dodge.width = 0.5)) +
    stat_summary(aes(fill = product), fun = mean, geom = "point", shape = 23, size = 4, position = position_dodge(width = 0.5)) +
    stat_summary(fun = mean, geom = "label", vjust = -2, position = position_dodge(width = 0.5), aes(color = product, label = round(after_stat(y), 1))) +
  geom_text(data = product_models,
    aes(x = meal_set, y = Inf, label = paste0("p = ", round(p.value, 4))), hjust = 0.5, vjust = 3, inherit.aes = FALSE, size = 4, color = product_models$text_color) +  
    scale_color_manual(values = pal_bmj("default")(9)[c(1,2,4,5,7,8)])+
  scale_fill_manual(values = pal_bmj("default")(9)[c(1,2,4,5,7,8)])+
    facet_wrap(~parameter, scales = "free")+
  guides(color = "none")+
    theme_bw() +
    theme(axis.text.x = element_text(angle = 45, hjust = 1), axis.title.x = element_blank(), axis.title.y = element_blank())

The plot indicate a significantly stronger glycemic response after consumption of gingerbread with yoghurt compared to fruit with yoghurt, as reflected by both higher iAUC values and higher peak glucose concentrations. No statistically significant differences were observed for the remaining compared product pairs.

Comparison between each breakfast and standard glucose drink

Linear mixed-effects models were used to assess differences in glycemic response parameters between food products. The analysis included baseline glucose, peak glucose, minimum glucose, and incremental area under the curve (iAUC), which were reshaped into a long format for modelling.

For each glycemic parameter, a separate model was fitted with product as a fixed effect and participant ID as a random intercept to account for repeated measurements within individuals. The glucose drink served as the reference category.

Estimated fixed effects were extracted and used to obtain p-values for pairwise comparisons between each test product and the reference. Significance levels were encoded using standard notation (* p < 0.05, ** p < 0.01, *** p < 0.001).

These model-derived results were subsequently used for graphical annotation of statistically significant differences between products.

product_to_OGTT_models<- parameters_2h_transformed %>%
   filter(!product%in% c("Nasi", "Pasta")) %>%
  pivot_longer(cols = c("base_glc", "peak_glc", "min_glc", "iAUC_transformed"), values_to = "value", names_to = "parameter") %>%
  mutate(product = fct_relevel(as.factor(product), c("glucose_drink", "Fruit_Yoghurt","Gingerbread_Yoghurt", "WG_Bread","R_Bread"))) %>%
  group_by(parameter) %>%
   nest() %>%
  mutate(model = map(data, ~ lmer(value ~ product + (1|part_id), data = .x)),
     tidy_results = map(model,~ broom.mixed::tidy(.x, effects = "fixed")))%>%
  select(parameter, tidy_results) %>%
  unnest(tidy_results)

product_to_OGTT_models<- product_to_OGTT_models %>% filter(term != "(Intercept)")  %>% mutate(
    label = paste0("p = ", round(p.value, 3)), group1 = "glucose_drink") %>%
  rename("group2" = "term") %>%
  mutate(group2 = str_remove(group2, "product")) %>%
  mutate(p.label = case_when(p.value < 0.001 ~ "***",
                                  p.value < 0.01  ~ "**",
                                  p.value < 0.05  ~ "*",
                                  TRUE ~ "ns")) %>%
  mutate(xmin = 1, xmax = c(2,3,4,5))

y_positions <- plot_data_consumption %>%
  filter(!product %in% c("Nasi", "Pasta")) %>%
  group_by(parameter) %>%
  summarise(y.position = max(value, na.rm = TRUE) * 1.15)

product_to_OGTT_models <- left_join(product_to_OGTT_models, y_positions, by = "parameter") %>%
  group_by(parameter) %>%
  mutate(y.position = y.position + c(0,1,2,3)*0.05*max(y.position)) %>% ungroup

product_to_OGTT_models_sig <- product_to_OGTT_models %>%
  filter(p.label != "ns")
  
product_to_OGTT_models %>% arrange(p.value) %>% kbl() %>%
  kableExtra::kable_paper(font_size = 12) %>%
  scroll_box(width = "800px", height = "500px")

parameter	effect	group2	estimate	std.error	statistic	df	p.value	label	group1	p.label	xmin	xmax	y.position
peak_glc	fixed	Fruit_Yoghurt	-2.6950000	0.3048043	-8.8417395	70.01196	0.0000000	p = 0	glucose_drink	***	1	2	15.41000
iAUC_transformed	fixed	R_Bread	8.0215905	0.9294535	8.6304377	70.47392	0.0000000	p = 0	glucose_drink	***	1	5	31.50295
iAUC_transformed	fixed	WG_Bread	7.0758291	0.9053118	7.8159027	70.17883	0.0000000	p = 0	glucose_drink	***	1	4	30.13326
peak_glc	fixed	WG_Bread	-1.6150244	0.3111192	-5.1910155	70.02792	0.0000020	p = 0	glucose_drink	***	1	4	16.95100
peak_glc	fixed	Gingerbread_Yoghurt	-1.4250000	0.3048043	-4.6751313	70.01196	0.0000139	p = 0	glucose_drink	***	1	3	16.18050
peak_glc	fixed	R_Bread	-1.4106355	0.3197293	-4.4119686	70.06610	0.0000363	p = 0	glucose_drink	***	1	5	17.72150
min_glc	fixed	WG_Bread	0.2461201	0.1151450	2.1374806	70.01683	0.0360526	p = 0.036	glucose_drink		1	4	9.23450
min_glc	fixed	R_Bread	0.1998197	0.1183397	1.6885263	70.03592	0.0957585	p = 0.096	glucose_drink	ns	1	5	9.65425
min_glc	fixed	Gingerbread_Yoghurt	0.1700000	0.1128046	1.5070308	70.00885	0.1363021	p = 0.136	glucose_drink	ns	1	3	8.81475
iAUC_transformed	fixed	Fruit_Yoghurt	-1.2846222	0.8873017	-1.4477852	70.05358	0.1521379	p = 0.152	glucose_drink	ns	1	2	27.39387
min_glc	fixed	Fruit_Yoghurt	0.1550000	0.1128046	1.3740575	70.00885	0.1738072	p = 0.174	glucose_drink	ns	1	2	8.39500
base_glc	fixed	WG_Bread	0.1661695	0.1361955	1.2200813	70.03232	0.2265267	p = 0.227	glucose_drink	ns	1	4	9.36100
base_glc	fixed	Gingerbread_Yoghurt	0.1550000	0.1334295	1.1616623	70.01963	0.2493173	p = 0.249	glucose_drink	ns	1	3	8.93550
iAUC_transformed	fixed	Gingerbread_Yoghurt	-0.7430899	0.8873017	-0.8374716	70.05358	0.4051757	p = 0.405	glucose_drink	ns	1	3	28.76357
base_glc	fixed	Fruit_Yoghurt	0.0250000	0.1334295	0.1873649	70.01963	0.8519164	p = 0.852	glucose_drink	ns	1	2	8.51000
base_glc	fixed	R_Bread	0.0224388	0.1399686	0.1603133	70.06266	0.8730956	p = 0.873	glucose_drink	ns	1	5	9.78650

plot_data_consumption %>%
     filter(!product%in% c("Nasi", "Pasta")) %>%
    ggplot(aes(x = product, y = value)) +
    geom_boxplot(aes(fill = product), outlier.shape = NA, width = 0.5, alpha = 0.5) +
    geom_jitter(aes(fill = product),
                color = "black", shape = 21, size = 2.5, alpha = 0.8, width = 0.2) +
    stat_summary(aes(fill = product), fun = mean, geom = "point", shape = 23, size = 4) +
    stat_summary(fun = mean, geom = "label", vjust = -2, aes(color = product, label = round(after_stat(y), 1))) +
    stat_pvalue_manual(product_to_OGTT_models_sig, label = "p.label",xmin = "group1", xmax = "group2", y.position = "y.position", tip.length = 0.01, bracket.size = 0.5, size = 5)+
    scale_color_manual(values = pal_bmj("default")(9)[c(1,2,4,5, 9)])+
  scale_fill_manual(values = pal_bmj("default")(9)[c(1,2,4,5, 9)])+
    facet_wrap(~parameter, scales = "free")+
  guides(color = "none")+
    theme_bw() +
    theme(axis.text.x = element_text(angle = 45, hjust = 1), axis.title.x = element_blank(), axis.title.y = element_blank())

The highest peak in glycemic response was observed after consumption of the glucose drink. In contrast, the incremental area under the curve (iAUC) indicated a significantly greater overall glycemic response following bread-based breakfasts compared to the glucose drink.

Summary / Conclusions

Although postprandial glycemic curves over the 2-hour period showed minor visual differences, most derived glycemic parameters (peak glucose, minimum glucose, and iAUC) were not significantly different between breakfast meals; an exception was observed for rice (nasi)-based meals.
Gingerbread with yoghurt elicited a significantly higher glycemic response compared to fruit with yoghurt, as reflected by both iAUC and peak glucose.
Nasi consumption resulted in a significantly higher peak glucose compared to pasta, indicating a stronger acute glycemic excursion despite similar overall exposure.
The glucose drink produced the highest peak glucose concentration; however, bread-based meals showed a higher overall glycemic exposure (iAUC), indicating a more prolonged postprandial response despite lower peak values and faster return toward baseline.
Considerable inter-individual variability in glycemic responses was observed, which may be influenced by factors such as metabolic health, age, genetic background, and gut microbiome composition.

Postprandial glucose response

Anna Malinowska

2026-06-09