Sweet Foods Liking and Diet Quality

Load packages and call data

Decoding surveys

## Liking score
liking.df <- qualtrics %>% 
  select("ID", "FL-01_3":"FL-60_3") %>%
  mutate(veg.liking = rowMeans(select(., "FL-01_3":"FL-06_3"), na.rm = TRUE), 
         veg.liking = (veg.liking+100)/2, 
         fruit.liking = rowMeans(select(., "FL-07_3":"FL-11_3"), na.rm = TRUE), 
         fruit.liking = (fruit.liking+100)/2, 
         saltyfat.liking = rowMeans(select(., "FL-12_3":"FL-15_3"), na.rm = TRUE), 
         saltyfat.liking = (saltyfat.liking+100)/2,
         hfprotein.liking = rowMeans(select(., "FL-16_3":"FL-20_3"), na.rm = TRUE),   
         hfprotein.liking = (hfprotein.liking+100)/2, 
         alcohol.liking = rowMeans(select(., "FL-22_3":"FL-25_3"), na.rm = TRUE), 
         alcohol.liking = (alcohol.liking+100)/2,
         wg.liking = rowMeans(select(., "FL-26_3":"FL-29_3"), na.rm = TRUE), 
         wg.liking = (wg.liking+100)/2,
         carbs.liking = rowMeans(select(., "FL-30_3":"FL-34_3"), na.rm = TRUE), 
         carbs.liking = (carbs.liking+100)/2,
         healthyfat.liking = rowMeans(select(., "FL-35_3":"FL-38_3"), na.rm = TRUE), 
         healthyfat.liking = (healthyfat.liking+100)/2,
        sweets.liking = rowMeans(select(., "FL-39_3":"FL-43_3"), na.rm = TRUE),
        sweets.liking = (sweets.liking+100)/2,
        ssb.liking = rowMeans(select(., "FL-44_3":"FL-47_3"), na.rm = TRUE),
        ssb.liking = (ssb.liking+100)/2,
        sfbl.liking = rowMeans(select(., "FL-39_3":"FL-47_3"), na.rm = TRUE),
        sfbl.liking = (sfbl.liking+100)/2,
         unhealthyfat.liking = rowMeans(select(., "FL-48_3":"FL-52_3"), na.rm = TRUE),
        unhealthyfat.liking = (unhealthyfat.liking+100)/2,
         protein.liking = rowMeans(select(., "FL-53_3":"FL-57_3"), na.rm = TRUE), 
        protein.liking = (protein.liking+100)/2,
         spicy.liking = rowMeans(select(., "FL-58_3":"FL-60_3"), na.rm = TRUE),
        spicy.liking = (spicy.liking+100)/2,) %>%
  select(c("ID", ends_with(".liking"))) 

## Three factor questionnaire
TFQ.df <- qualtrics %>% 
  select("ID", "TFQ-01":"TFQ-21") %>%
  mutate(uncontrolled = rowSums(select(., "TFQ-03", "TFQ-06", "TFQ-08", "TFQ-09", "TFQ-12", "TFQ-13", "TFQ-15", "TFQ-19", "TFQ-20"), na.rm = TRUE)) %>%  
  mutate(uncontrolled = ((uncontrolled-1)/36)*100, na.rm = TRUE) %>%
  mutate(restraint = rowSums(select(., "TFQ-01", "TFQ-05", "TFQ-11"), na.rm = TRUE)) %>%
  mutate(restraint = ((restraint-1)/12)*100, na.rm = TRUE) %>%
  mutate(emotional = rowSums(select(., "TFQ-02", "TFQ-04", "TFQ-07", "TFQ-10", "TFQ-14", "TFQ-16"), na.rm = TRUE)) %>%
  mutate(emotional = ((emotional-1)/24)*100, na.rm = TRUE) %>%
    select(c("ID", "uncontrolled", "restraint", "emotional")) 

## Body appreciation scale 
BAS.df <- qualtrics %>% 
  select("ID", starts_with("BAS")) %>%
  mutate(BAS = rowSums(select(., 2:11), na.rm = TRUE))%>%
  select(c("ID", "BAS")) 

## Intuitive eating scale 
IES.df <- qualtrics %>% 
  select("ID", starts_with("IES")) %>%
  mutate(IES = (rowSums(select(., 2:24), na.rm = TRUE))/23) %>%
  select(c("ID", "IES"))

## sHEI score
HEI.df <- qualtrics %>%
  select("ID", "sex", "Fruit":"Water") %>%
  mutate(fruit.1 = case_when(
    Fruit == 1 ~ 0,
    Fruit == 2 ~ 2, 
    Fruit == 3 ~ 3.5,
    Fruit >= 4 ~ 5)) %>%
  mutate(fruit.2 = case_when (
    Juice == 1 ~ 0,
    Juice == 2 ~ 2,
    Juice == 3 ~ 3.5,
    Juice >= 4 ~ 5)) %>% 
  mutate(tfruitHEI = case_when(
    fruit.1 + fruit.2 <5 ~ 0,
    fruit.1 + fruit.2 >= 5 ~ 5)) %>%
  mutate(wfruitHEI = case_when(
    Fruit == 1 ~ 0,
    Fruit == 2 ~ 2.5, 
    Fruit >= 3 ~ 5)) %>%
  mutate(vegHEI = case_when(
    `Green veg` == 1 ~ 1.6,
    `Green veg` == 2 & `Starchy veg` >= 2 ~ 2.46,
    `Green veg` >= 2 & `Starchy veg` >= 2 ~ 3.24,
    `Green veg` >= 2 & `Starchy veg` == 1 ~ 3.56)) %>%
  mutate(bean.1 = case_when(
    `Green veg` == 1 ~ 0,
    `Green veg` >= 2 ~5)) %>%
  mutate(bean.2 = case_when(
    Beans == 1 ~ 0,
    Beans >= 2 ~ 5)) %>%
  mutate(beanHEI = case_when(
    bean.1 + bean.2 < 5 ~ 0,
    bean.1 + bean.2 >= 5 ~ 5)) %>%
  mutate(wgHEI = case_when(
    `Whole grains.day` == 1 ~ 0.51,
    sex == 1 & `Whole grains.day` >= 2 ~ 2.97,
    sex == 2 & `Whole grains.day` >= 2 & `Whole grains.day` <= 3 ~ 5.20,
    sex == 2 & `Whole grains.day` >= 4 ~ 6.94)) %>%
  mutate(dairyHEI = case_when(
    sex == 1 & `Milk.day` <= 3 ~ 3.22,
    sex == 2 & `Milk.day` <= 3 & `Low-fat milk.day` == 1 ~ 3.32,
    sex == 2 & `Milk.day` <= 3 & `Low-fat milk.day` >= 2 ~ 4.81,
    `Milk.day` >= 4 ~ 6.51)) %>%
  mutate(tproteinHEI = case_when(
    sex == 1 & `Seafood.freq` <= 4 ~ 4.11,
    sex == 1 & `Seafood.freq` >= 5 ~ 4.98,
    sex == 1 & is.na(`Seafood.freq`) ~ 4.11,
    sex == 2 ~ 4.97)) %>%
  mutate(sfpproteinHEI = case_when(
    sex == 1 & Nuts <= 2 ~ 0.49,
    sex == 2 & Nuts <= 2 ~ 1.50,
    Nuts >= 3 ~ 4.20)) %>%
  mutate(fatHEI = case_when(
    `Milk.day` >= 4 ~ 2.56,
    `Saturated fats` >= 2 & `Saturated fats` <= 3 & `Milk.freq` >= 1 & `Low-fat milk.freq` <= 2 ~ 2.63,
    `Saturated fats` >= 2 & `Saturated fats` <= 3 & is.na(`Milk.freq`) | is.na(`Low-fat milk.freq`) ~ 2.63,
    `Saturated fats` >= 2 & `Saturated fats` <= 3 & `Milk.freq` >= 1 & `Low-fat milk.freq` >= 3 ~ 4.54,
    `Saturated fats` == 1 & `Milk.freq` >= 1 ~ 5.93,
    `Saturated fats` == 1 & is.na(`Milk.freq`) ~ 5.93)) %>%
  mutate(grainHEI = case_when(
    `Green veg` == 1 ~ 2.13,
    `Grains.day` >= 3 & `Seafood.day` >= 2 & `Green veg` >= 2 ~ 2.27,
    `Grains.day` >= 3 & Nuts >= 1 & Nuts <= 2 & `Seafood.day` == 1 & `Green veg` >= 2 ~ 4.73,
    `Grains.day` >= 3 & Nuts >= 3 & `Seafood.day` == 1 & `Green veg` >= 2 ~ 8.45,
    `Grains.day` >= 1 & `Grains.day` <= 2 & `Green veg` >= 2 ~ 9.25)) %>%
  mutate(sodiumHEI = case_when(
    Fruit >= 1 & Fruit <= 2 & `Grains.day` >= 3 & Water == 3 ~ 0.70,
    Fruit >= 3 & `Grains.day` >= 3 & Water == 3 ~ 2.30,
    `Grains.day` >= 3 & Water >= 1 & Water <= 2 ~ 4.94,
    `Grains.day` >= 1 & `Grains.day` <= 2 ~ 6.07)) %>%
  mutate(ssb.cal = case_when(
    `SSB.day` == 1  ~ 0,
    `SSB.day` == 2 ~ 156,
    `SSB.day` == 3 ~ 312,
    `SSB.day` == 4 ~ 468,
    `SSB.day` == 5 ~ 624,
    `SSB.day` == 6 ~ 780,
    `SSB.day` == 7 ~ 936),
    sugar.cal = case_when(
    `Added sugars` == 1 ~ 130,
    `Added sugars` == 2 ~ 260,
    `Added sugars` == 3 ~ 520)) %>%
  mutate(sugar.intake = (ssb.cal + sugar.cal)) %>%
  mutate(sugarHEI = case_when(
    sugar.intake <= 130 ~ 10,
    sugar.intake > 130 & sugar.cal < 520 ~ 5,
    sugar.intake >= 520 ~ 0)) %>%
  mutate(sugar.tsp = case_when(
    `SSB.freq` <= 4 ~ 13.26,
    `SSB.freq` >= 5 & `SSB.freq` <=6 ~ 16.00,
    `SSB.day` ==2 ~ 16.00,
    `SSB.day` >=3 ~ 26.87)) %>%
  mutate(satfatHEI = case_when(
    `SSB.day` >= 3 ~ 1.82,
    `SSB.day` <= 2 & `Grains.day` <= 2 ~ 3.20,
    `SSB.day` <= 2 & `Grains.day` >= 3 & Nuts <= 2 ~ 4.64,
    `SSB.day` <= 2 & `Grains.day` >= 3 & Nuts >= 3 ~ 6.56)) %>%
  select("ID", "sugar.intake", "sugar.tsp", ends_with("HEI")) %>%
  mutate(sHEI = rowSums(select(., "tfruitHEI":"satfatHEI"), na.rm = TRUE)) %>%
  mutate(SuFatHEI = rowSums(select(., c("sugarHEI", "satfatHEI")))) 

## Putting scores together
qualtrics$ID <- as.character(qualtrics$ID)
anthro$ID <- as.character(anthro$ID)
oatmilk.liking.wide$ID <- as.character(oatmilk.liking.wide$ID) 
oatmilk.sweet.wide$ID <- as.character(oatmilk.sweet.wide$ID) 
oatmilk.creamy.wide$ID <- as.character(oatmilk.creamy.wide$ID) 
oatmilk.bitter.wide$ID <- as.character(oatmilk.bitter.wide$ID) 
liking.df$ID <- as.character(liking.df$ID)
HEI.df$ID <- as.character(HEI.df$ID)
TFQ.df$ID <- as.character(TFQ.df$ID)
BAS.df$ID <- as.character(BAS.df$ID)
IES.df$ID <- as.character(IES.df$ID)

full.df <- qualtrics  %>%
  select(c("ID", "age", "race", "ethn", "usb", "education", "hinc", "exec", "plantmilk", "plantmilk_freq")) %>%
  left_join(anthro, by = "ID") %>%
  left_join(oatmilk.liking.wide, by = "ID") %>%
  left_join(oatmilk.sweet.wide, by = "ID") %>%
  left_join(oatmilk.creamy.wide, by = "ID") %>%
  left_join(oatmilk.bitter.wide, by = "ID") %>%
  left_join(liking.df, by = "ID") %>%
  left_join(HEI.df, by = "ID") %>%
  left_join(TFQ.df, by = "ID") %>%
  left_join(BAS.df, by = "ID") %>%
  left_join(IES.df, by = "ID") %>%
  mutate(BMIcat = case_when(
      BMI < 18.5 ~ "0",
      BMI >= 18.5 & BMI  < 25 & race != "Asian" ~ "1",
      BMI >= 25 & BMI < 30 & race != "Asian" ~ "2",
      BMI > 30 & race != "Asian" ~ "3",
      BMI >= 18.5 & BMI  < 23 & race == "Asian" ~ "1",
      BMI >= 23 & BMI < 25 & race == "Asian" ~ "2",
      BMI > 25 & race == "Asian" ~ "3"))  %>%
  mutate(sugarcat = ntile(sugar.intake, 3)) %>%
  mutate(DQ = ntile(sHEI, 3)) 

Reliability

# Define your function
reliability <- function(data, xvar, yvar, title) {
  ggplot(data, aes_string(x = xvar, y = yvar, label = "ID")) +
    geom_point() +
    geom_abline(intercept = 0, slope = 1, color = "red", linetype = "dashed") +
    geom_text_repel() +
    stat_cor() +
    ylim(0, 100) +
    xlim(0, 100) +
    ggtitle(title) +
    theme_light() +
    theme(
      axis.text = element_text(size = 10),
      axis.title = element_text(size = 12),
      plot.title = element_text(size = 14, face = "bold")
    )
}

## liking 

liking.plots <- list(
  reliability(oatmilk.liking.wide, "liking_S0F0.1", "liking_S0F0.2", "S0F0"),
  reliability(oatmilk.liking.wide, "liking_S1F0.1", "liking_S1F0.2", "S1F0"),
  reliability(oatmilk.liking.wide, "liking_S1F1.1", "liking_S1F1.2", "S1F1"),
  reliability(oatmilk.liking.wide, "liking_S1F2.1", "liking_S1F2.2", "S1F2"),
  reliability(oatmilk.liking.wide, "liking_S1F3.1", "liking_S1F3.2", "S1F3"),
  reliability(oatmilk.liking.wide, "liking_S2F0.1", "liking_S2F0.2", "S2F0"),
  reliability(oatmilk.liking.wide, "liking_S2F1.1", "liking_S2F1.2", "S2F1"),
  reliability(oatmilk.liking.wide, "liking_S2F2.1", "liking_S2F2.2", "S2F2"),
  reliability(oatmilk.liking.wide, "liking_S2F3.1", "liking_S2F3.2", "S2F3")
)

# liking plot grid
(liking.plots[[1]] + liking.plots[[2]] + liking.plots[[3]]) /
(liking.plots[[4]] + liking.plots[[5]] + liking.plots[[6]]) /
(liking.plots[[7]] + liking.plots[[8]] + liking.plots[[9]])

# sweet 
sweet.plots <- list(
  reliability(oatmilk.sweet.wide, "sweet_S0F0.1", "sweet_S0F0.2", "S0F0"),
  reliability(oatmilk.sweet.wide, "sweet_S1F0.1", "sweet_S1F0.2", "S1F0"),
  reliability(oatmilk.sweet.wide, "sweet_S1F1.1", "sweet_S1F1.2", "S1F1"),
  reliability(oatmilk.sweet.wide, "sweet_S1F2.1", "sweet_S1F2.2", "S1F2"),
  reliability(oatmilk.sweet.wide, "sweet_S1F3.1", "sweet_S1F3.2", "S1F3"),
  reliability(oatmilk.sweet.wide, "sweet_S2F0.1", "sweet_S2F0.2", "S2F0"),
  reliability(oatmilk.sweet.wide, "sweet_S2F1.1", "sweet_S2F1.2", "S2F1"),
  reliability(oatmilk.sweet.wide, "sweet_S2F2.1", "sweet_S2F2.2", "S2F2"),
  reliability(oatmilk.sweet.wide, "sweet_S2F3.1", "sweet_S2F3.2", "S2F3")
  )

# sweet plot grid
(sweet.plots[[1]] + sweet.plots[[2]] + sweet.plots[[3]]) /
(sweet.plots[[4]] + sweet.plots[[5]] + sweet.plots[[6]]) /
(sweet.plots[[7]] + sweet.plots[[8]] + sweet.plots[[9]])

# creamy
creamy.plots <- list(
  reliability(oatmilk.creamy.wide, "creamy_S0F0.1", "creamy_S0F0.2", "S0F0"),
  reliability(oatmilk.creamy.wide, "creamy_S1F0.1", "creamy_S1F0.2", "S1F0"),
  reliability(oatmilk.creamy.wide, "creamy_S1F1.1", "creamy_S1F1.2", "S1F1"),
  reliability(oatmilk.creamy.wide, "creamy_S1F2.1", "creamy_S1F2.2", "S1F2"),
  reliability(oatmilk.creamy.wide, "creamy_S1F3.1", "creamy_S1F3.2", "S1F3"),
  reliability(oatmilk.creamy.wide, "creamy_S2F0.1", "creamy_S2F0.2", "S2F0"),
  reliability(oatmilk.creamy.wide, "creamy_S2F1.1", "creamy_S2F1.2", "S2F1"),
  reliability(oatmilk.creamy.wide, "creamy_S2F2.1", "creamy_S2F2.2", "S2F2"),
  reliability(oatmilk.creamy.wide, "creamy_S2F3.1", "creamy_S2F3.2", "S2F3")
)

# creamy plot grid
(creamy.plots[[1]] + creamy.plots[[2]] + creamy.plots[[3]]) /
(creamy.plots[[4]] + creamy.plots[[5]] + creamy.plots[[6]]) /
(creamy.plots[[7]] + creamy.plots[[8]] + creamy.plots[[9]])

# bitter
bitter.plots <- list(
  reliability(oatmilk.bitter.wide, "bitter_S0F0.1", "bitter_S0F0.2", "S0F0"),
  reliability(oatmilk.bitter.wide, "bitter_S1F0.1", "bitter_S1F0.2", "S1F0"),
  reliability(oatmilk.bitter.wide, "bitter_S1F1.1", "bitter_S1F1.2", "S1F1"),
  reliability(oatmilk.bitter.wide, "bitter_S1F2.1", "bitter_S1F2.2", "S1F2"),
  reliability(oatmilk.bitter.wide, "bitter_S1F3.1", "bitter_S1F3.2", "S1F3"),
  reliability(oatmilk.bitter.wide, "bitter_S2F0.1", "bitter_S2F0.2", "S2F0"),
  reliability(oatmilk.bitter.wide, "bitter_S2F1.1", "bitter_S2F1.2", "S2F1"),
  reliability(oatmilk.bitter.wide, "bitter_S2F2.1", "bitter_S2F2.2", "S2F2"),
  reliability(oatmilk.bitter.wide, "bitter_S2F3.1", "bitter_S2F3.2", "S2F3")
)

# bitter plot grid
(bitter.plots[[1]] + bitter.plots[[2]] + bitter.plots[[3]]) /
(bitter.plots[[4]] + bitter.plots[[5]] + bitter.plots[[6]]) /
(bitter.plots[[7]] + bitter.plots[[8]] + bitter.plots[[9]])

Counts

oatmilk.liking.wide %>%
  group_by(overall.pref) %>%
  count()

# A tibble: 9 × 2
# Groups:   overall.pref [9]
  overall.pref     n
  <chr>        <int>
1 S0F0             2
2 S1F0             4
3 S1F1             5
4 S1F2             2
5 S1F3             2
6 S2F0             1
7 S2F1             3
8 S2F2             4
9 S2F3             7

Participants flow - Sankey diagram

sankey <- oatmilk.liking.wide %>% 
  select(c(overall.pref, sweetpref.atF0, sweetpref.atF0, sweetpref.atF1, sweetpref.atF2, sweetpref.atF3, fatpref.atS1, fatpref.atS2)) %>%
  select(overall.pref, fatpref.atS1, fatpref.atS2) %>%
  mutate(id = row_number()) %>%
  pivot_longer(
    cols = c(overall.pref, fatpref.atS1, fatpref.atS2),
    names_to = "stage",
    values_to = "node"
  ) %>%
  mutate(stage = factor(stage, levels = c("overall.pref", "fatpref.atS1", "fatpref.atS2")))

sankey %>%
  ggplot(aes(x = stage, stratum = node, alluvium = id, fill = as.factor(node), label = node)) +
  geom_flow(stat = "alluvium", lode.guidance = "forward", alpha = 0.7) +
  geom_stratum(width = 1/6) +
  geom_text(stat = "stratum", size = 4) +
  scale_x_discrete(labels = c("Overall Preference", "Fat Preference\nat Sweet Level 1", "Fat Preference\nat Sweet Level 2"), expand = c(0.1, 0.1)) +
  labs(title = "Participant Flow: Overall Preference → Fat Preference at Sweet Level 1 → Fat Preference at Sweet Level 2",
       x = NULL, y = "Number of Participants") +
  theme_classic() +
  theme(legend.position = "none")

Dose response - liking

## all participants
liking.dose <- oatmilk.liking.long %>%
  group_by(sample) %>%
  summarize(sem = sd(liking)/sqrt(30), liking = mean(liking))

liking.dose <- liking.dose %>%
  separate(sample, into = c("sugar.level", "fat.level"), sep = 2)  

liking.dose$fat.level <- factor(liking.dose$fat.level , levels = c("F0", "F1", "F2", "F3"))
liking.dose$sugar.level <- factor(liking.dose$sugar.level , levels = c("S0", "S1", "S2"))

## Overall
liking.dose %>%
  filter(sugar.level != "S0") %>%
  ggplot(aes(x = fat.level, y = liking, group = sugar.level, color = sugar.level, shape = sugar.level)) +
  geom_line()+
  geom_point() +
  scale_color_manual(values = c("#FFCA28", "#F44336")) +
  geom_errorbar(aes(ymin = liking - sem, ymax = liking + sem), width = 0.05) +
  ylim(0, 100) +
  theme_classic() 

## By individual
oatmilk.liking.long %>%
  group_by(ID, sample) %>%
  summarize(sem = sd(liking)/sqrt(30), liking = mean(liking)) %>%
    filter(sample != "S0F0") %>%
  separate(sample, into = c("sugar.level", "fat.level"), sep = 2)  %>%
  ggplot(aes(x = fat.level, y = liking, group = sugar.level, color = sugar.level)) +
  geom_point() +
  geom_line() +
  scale_color_manual(values = c("#FFCA28", "#F44336")) +
  facet_wrap(~ID)+
  theme_classic()

`summarise()` has grouped output by 'ID'. You can override using the `.groups`
argument.

Dose response - sweet

sweet.dose <- oatmilk.sweet.long %>%
  group_by(sample) %>%
  summarize(sem = sd(sweet)/sqrt(30), sweet = mean(sweet))

sweet.dose <- sweet.dose %>%
  separate(sample, into = c("sugar.level", "fat.level"), sep = 2)  

sweet.dose$fat.level <- factor(sweet.dose$fat.level , levels = c("F0", "F1", "F2", "F3"))
sweet.dose$sugar.level <- factor(sweet.dose$sugar.level , levels = c("S0", "S1", "S2"))

## Overall
sweet.dose %>%
  filter(sugar.level != "S0") %>%
  ggplot(aes(x = fat.level, y = sweet, group = sugar.level, color = sugar.level, shape = sugar.level)) +
  geom_line()+
  geom_point() +
  scale_color_manual(values = c("#FFCA28", "#F44336")) +
  geom_errorbar(aes(ymin = sweet - sem, ymax = sweet + sem), width = 0.05) +
  ylim(0, 100) +
  theme_classic()

## By individual
oatmilk.sweet.long %>%
  group_by(ID, sample) %>%
  summarize(sem = sd(sweet)/sqrt(25), sweet = mean(sweet)) %>%
    filter(sample != "S0F0") %>%
  separate(sample, into = c("sugar.level", "fat.level"), sep = 2)  %>%
  ggplot(aes(x = fat.level, y = sweet, group = sugar.level, color = sugar.level)) +
  geom_point() +
  geom_line() +
  scale_color_manual(values = c("#FFCA28", "#F44336")) +
  facet_wrap(~ID)+
  theme_classic()

`summarise()` has grouped output by 'ID'. You can override using the `.groups`
argument.

Dose response - creamy

creamy.dose <- oatmilk.creamy.long %>%
  group_by(sample) %>%
  summarize(sem = sd(creamy)/sqrt(30), creamy = mean(creamy))

creamy.dose <- creamy.dose %>%
  separate(sample, into = c("sugar.level", "fat.level"), sep = 2)  

creamy.dose$fat.level <- factor(creamy.dose$fat.level , levels = c("F0", "F1", "F2", "F3"))
creamy.dose$sugar.level <- factor(creamy.dose$sugar.level , levels = c("S0", "S1", "S2"))

## Overall
creamy.dose %>%
  filter(sugar.level != "S0") %>%
  ggplot(aes(x = fat.level, y = creamy, group = sugar.level, color = sugar.level, shape = sugar.level)) +
  geom_line()+
  geom_point() +
  scale_color_manual(values = c("#FFCA28", "#F44336")) +
  geom_errorbar(aes(ymin = creamy - sem, ymax = creamy + sem), width = 0.05) +
  ylim(0, 100) +
  theme_classic()

## By individual
oatmilk.creamy.long %>%
  group_by(ID, sample) %>%
  summarize(sem = sd(creamy)/sqrt(30), creamy = mean(creamy)) %>%
    filter(sample != "S0F0") %>%
  separate(sample, into = c("sugar.level", "fat.level"), sep = 2)  %>%
  ggplot(aes(x = fat.level, y = creamy, group = sugar.level, color = sugar.level)) +
  geom_point() +
  geom_line() +
  scale_color_manual(values = c("#FFCA28", "#F44336")) +
  facet_wrap(~ID)+
  theme_classic()

`summarise()` has grouped output by 'ID'. You can override using the `.groups`
argument.

Dose response - bitter

bitter.dose <- oatmilk.bitter.long %>%
  group_by(sample) %>%
  summarize(sem = sd(bitter)/sqrt(30), bitter = mean(bitter))

bitter.dose <- bitter.dose %>%
  separate(sample, into = c("sugar.level", "fat.level"), sep = 2)  

bitter.dose$fat.level <- factor(bitter.dose$fat.level , levels = c("F0", "F1", "F2", "F3"))
bitter.dose$sugar.level <- factor(bitter.dose$sugar.level , levels = c("S0", "S1", "S2"))

## Overall
bitter.dose %>%
  filter(sugar.level != "S0") %>%
  ggplot(aes(x = fat.level, y = bitter, group = sugar.level, color = sugar.level, shape = sugar.level)) +
  geom_line()+
  geom_point() +
  scale_color_manual(values = c("#FFCA28", "#F44336")) +
  geom_errorbar(aes(ymin = bitter - sem, ymax = bitter + sem), width = 0.05) +
  ylim(0, 100) +
  theme_classic()

## By individual
oatmilk.bitter.long %>%
  group_by(ID, sample) %>%
  summarize(sem = sd(bitter)/sqrt(30), bitter = mean(bitter)) %>%
    filter(sample != "S0F0") %>%
  separate(sample, into = c("sugar.level", "fat.level"), sep = 2)  %>%
  ggplot(aes(x = fat.level, y = bitter, group = sugar.level, color = sugar.level)) +
  geom_point() +
  geom_line() +
  scale_color_manual(values = c("#FFCA28", "#F44336")) +
  facet_wrap(~ID)+
  theme_classic()

`summarise()` has grouped output by 'ID'. You can override using the `.groups`
argument.

Corrplot - liking with demographics

## FDR adjustment
adjust.cor.pmat <- function(p.mat, method = "fdr") {
  p.mat[upper.tri(p.mat)] <- p.adjust(p.mat[upper.tri(p.mat)], method = method)
  p.mat[lower.tri(p.mat)] <- t(p.mat)[lower.tri(p.mat)]
  diag(p.mat) <- 0
  return(p.mat)
}

## All corr
corr.df <- full.df %>%
  select(age, BMI, WC, hinc, restraint, emotional, uncontrolled, sfbl.liking, plantmilk_freq, sugar.intake, sHEI, liking_S0F0, liking_S1F0, liking_S1F1, liking_S1F2, liking_S1F3, liking_S2F0, liking_S2F1, liking_S2F2, liking_S2F3, sweetpref.atF0, sweetpref.atF1, sweetpref.atF2, sweetpref.atF3, fatpref.atS1, fatpref.atS2) %>%
  rename("household income" = hinc, "sugar intake" = sugar.intake, 
         "healthy eating index" = sHEI, "sweet foods liking" = sfbl.liking) %>%
  mutate(across(everything(), as.numeric))

all.corr <- corr.df %>%
      cor(use = "pairwise.complete.obs")

p.mat <- cor_pmat(corr.df) %>%
  column_to_rownames()
p.mat.adj <- adjust.cor.pmat(p.mat, method = "fdr")

ggcorrplot(all.corr, type = "lower", p.mat = p.mat.adj, insig = "blank", lab = TRUE, lab_size = 2.5, colors = c("#E46726", "white", "#6D9EC1"))

No id variables; using all as measure variables