Fat Liking Study

Questions

Plant-based milk is a growing marker, especially in young women. Although generally considered a healthier alternative compared to cow’s milk by consumers, some plant-based milk may contain high levels of fat and sugar, which may off-set some of the benefits of switching from cow’s milk to plant-based milk. The goal of this study is to test whether liking for fat and sugar may interact with plant-based milk consumption to affect diet quality and body composition.

H1. Does H2.

Load packages and call data

Reliability

# Define your function
reliability <- function(data, xvar, yvar, title) {
  ggplot(data, aes_string(x = xvar, y = yvar, label = "ID")) +
    geom_point(alpha = 0.3) +
    geom_abline(intercept = 0, slope = 1, color = "red", linetype = "dashed") +
    stat_cor() +
    ylim(0, 100) +
    xlim(0, 100) +
    xlab("Replicate 1") +
    ylab("Replicate 2") +
    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")
    )
}

# geom_label_repel() 


## liking 

liking.plots <- list(
  reliability(oatmilk.liking.wide, "liking_S0F0.1", "liking_S0F0.2", "Control (water)"),
  reliability(oatmilk.liking.wide, "liking_S1F0.1", "liking_S1F0.2", "14% sucrose 0% fat"),
  reliability(oatmilk.liking.wide, "liking_S1F1.1", "liking_S1F1.2", "14% sucrose 2% fat"),
  reliability(oatmilk.liking.wide, "liking_S1F2.1", "liking_S1F2.2", "14% sucrose 6% fat"),
  reliability(oatmilk.liking.wide, "liking_S1F3.1", "liking_S1F3.2", "14% sucrose 12% fat"),
  reliability(oatmilk.liking.wide, "liking_S2F0.1", "liking_S2F0.2", "28% sucrose 0% fat "),
  reliability(oatmilk.liking.wide, "liking_S2F1.1", "liking_S2F1.2", "28% sucrose 2% fat"),
  reliability(oatmilk.liking.wide, "liking_S2F2.1", "liking_S2F2.2", "28% sucrose 6% fat"),
  reliability(oatmilk.liking.wide, "liking_S2F3.1", "liking_S2F3.2", "28% sucrose 12% fat")
)

# liking plot grid
liking.reliability <- ((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]]))

print(liking.reliability)

ggsave("liking.reliability.png", height = 8, width = 8)

# 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]])

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,
        sufat.liking = rowMeans(select(., c("FL-39_3","FL-40_3", "FL-43_3")), na.rm = TRUE),
        sufat.liking = (sufat.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)

### Missing stress survey data

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") %>%
  left_join(bia, by  = "ID") %>%
  mutate(BMIcat = case_when(
      BMI < 18.5 ~ "under",
      BMI >= 18.5 & BMI  < 25 ~ "normal",
      BMI >= 25 & BMI < 30 ~ "over",
      BMI > 30  ~ "obese"))  %>%
  mutate(sugarcat = ntile(sugar.intake, 3)) %>%
  mutate(DQ = ntile(sHEI, 3)) %>%
  mutate(sweetPref = case_when(
      sweetpref.atF2 == 1 ~ "Low",
      sweetpref.atF2 == 2 ~ "High")) %>%
  mutate(fatPref = case_when(
    fatpref.atS1 == 1 ~ "Low",
    fatpref.atS1 == 2 ~ "Medium",
    fatpref.atS1 == 3 ~ "High")) %>%
    mutate(sfblCat = ntile(sfbl.liking, 3),
         sfblCat = factor(sfblCat, levels = c(1, 2, 3), labels = c("Low", "Medium", "High"))) %>%
    mutate(healthyfatCat = ntile(healthyfat.liking, 3),
         healthyfatCat = factor(healthyfatCat, levels = c(1, 2, 3), labels = c("Low", "Medium", "High"))) %>%
      mutate(unhealthyfatCat = ntile(unhealthyfat.liking, 3),
         unhealthyfatCat = factor(unhealthyfatCat, levels = c(1, 2, 3), labels = c("Low", "Medium", "High"))) %>%
      mutate(visceralCat = ntile(visceral_fat, 3),
         visceralCat = factor(visceralCat, levels = c(1, 2, 3), labels = c("Low", "Medium", "High"))) %>%
      mutate(pbfCat = ntile(PBF, 3),
      pbfCat = factor(pbfCat, levels = c(1, 2, 3), labels = c("Low", "Medium", "High"))) %>%
      mutate(leanmassCat = ntile(lean_mass, 3),
      leanmassCat = factor(leanmassCat, levels = c(1, 2, 3), labels = c("Low", "Medium", "High"))) %>%    
      mutate(ffmCat = ntile(FFM, 3),
      ffmCat = factor(ffmCat, levels = c(1, 2, 3), labels = c("Low", "Medium", "High")))

## Joining sweetPref with oat milk dfs
sweetPref <- full.df %>% select(ID, sweetPref)

oatmilk.liking.long <- oatmilk.liking.long %>%  
  left_join(sweetPref, by = "ID") 
oatmilk.sweet.long <- oatmilk.sweet.long %>%  
  left_join(sweetPref, by = "ID") 
oatmilk.creamy.long <- oatmilk.creamy.long %>%  
  left_join(sweetPref, by = "ID") 
oatmilk.bitter.long <- oatmilk.bitter.long %>%  
  left_join(sweetPref, by = "ID") 

oatmilk.liking.wide <- oatmilk.liking.wide %>%  
  left_join(sweetPref, by = "ID") 
oatmilk.sweet.wide <- oatmilk.sweet.wide %>%  
  left_join(sweetPref, by = "ID") 
oatmilk.creamy.wide <- oatmilk.creamy.wide %>%  
  left_join(sweetPref, by = "ID") 
oatmilk.bitter.wide <- oatmilk.bitter.wide %>%  
  left_join(sweetPref, by = "ID") 

## Joining fatPref with oat milk dfs
fatPref <- full.df %>% select(ID, fatPref)

oatmilk.liking.long <- oatmilk.liking.long %>%  
  left_join(fatPref, by = "ID") 
oatmilk.sweet.long <- oatmilk.sweet.long %>%  
  left_join(fatPref, by = "ID") 
oatmilk.creamy.long <- oatmilk.creamy.long %>%  
  left_join(fatPref, by = "ID") 
oatmilk.bitter.long <- oatmilk.bitter.long %>%  
  left_join(fatPref, by = "ID") 

oatmilk.liking.wide <- oatmilk.liking.wide %>%  
  left_join(fatPref, by = "ID") 
oatmilk.sweet.wide <- oatmilk.sweet.wide %>%  
  left_join(fatPref, by = "ID") 
oatmilk.creamy.wide <- oatmilk.creamy.wide %>%  
  left_join(fatPref, by = "ID") 
oatmilk.bitter.wide <- oatmilk.bitter.wide %>%  
  left_join(fatPref, by = "ID") 

full.df$sweetPref <- factor(full.df$sweetPref , levels = c("Low", "High"))
full.df$fatPref <- factor(full.df$fatPref , levels = c("Low", "Medium", "High"))
full.df$prefcal.cat <- factor(full.df$prefcal.cat , levels = c("1", "2", "3"))

Preference data frames

pref.df <- oatmilk.liking.wide %>%
  select(ID, 
         # Rep 1
         liking_S1F0.1, liking_S2F0.1, liking_S1F1.1, liking_S1F2.1, liking_S1F3.1,
         liking_S1F2.1, liking_S2F2.1,
         # Rep 2
         liking_S1F0.2, liking_S2F0.2, liking_S1F1.2, liking_S1F2.2, liking_S1F3.2,
         liking_S1F2.2, liking_S2F2.2,
         # Overall (averaged across reps)
         liking_S1F0, liking_S2F0, liking_S1F1, liking_S1F2, liking_S1F3,
         liking_S2F2,
         # Other variables
         pref.cal, prefcal.cat) %>%
  rowwise() %>%
  mutate(
    # Fat preference rep 1 (at S1)
    fatPref.1 = case_when(
      `liking_S1F1.1` == max(c(`liking_S1F1.1`, `liking_S1F2.1`, `liking_S1F3.1`), na.rm = TRUE) ~ "low",
      `liking_S1F2.1` == max(c(`liking_S1F1.1`, `liking_S1F2.1`, `liking_S1F3.1`), na.rm = TRUE) ~ "medium",
      `liking_S1F3.1` == max(c(`liking_S1F1.1`, `liking_S1F2.1`, `liking_S1F3.1`), na.rm = TRUE) ~ "high"
    ),
    
    # Fat preference rep 2 (at S1)
    fatPref.2 = case_when(
      `liking_S1F1.2` == max(c(`liking_S1F1.2`, `liking_S1F2.2`, `liking_S1F3.2`), na.rm = TRUE) ~ "low",
      `liking_S1F2.2` == max(c(`liking_S1F1.2`, `liking_S1F2.2`, `liking_S1F3.2`), na.rm = TRUE) ~ "medium",
      `liking_S1F3.2` == max(c(`liking_S1F1.2`, `liking_S1F2.2`, `liking_S1F3.2`), na.rm = TRUE) ~ "high"
    ),
    
    # Fat preference overall (at S1) - using averaged liking scores
    fatPref = case_when(
      `liking_S1F1` == max(c(`liking_S1F1`, `liking_S1F2`, `liking_S1F3`), na.rm = TRUE) ~ "low",
      `liking_S1F2` == max(c(`liking_S1F1`, `liking_S1F2`, `liking_S1F3`), na.rm = TRUE) ~ "medium",
      `liking_S1F3` == max(c(`liking_S1F1`, `liking_S1F2`, `liking_S1F3`), na.rm = TRUE) ~ "high"
    ),
    
    # Sweet preference rep 1 (at F2)
    sweetPref.1 = case_when(
      `liking_S1F2.1` == max(c(`liking_S1F2.1`, `liking_S2F2.1`), na.rm = TRUE) ~ "low",
      `liking_S2F2.1` == max(c(`liking_S1F2.1`, `liking_S2F2.1`), na.rm = TRUE) ~ "high"
    ),
    
    # Sweet preference rep 2 (at F2)
    sweetPref.2 = case_when(
      `liking_S1F2.2` == max(c(`liking_S1F2.2`, `liking_S2F2.2`), na.rm = TRUE) ~ "low",
      `liking_S2F2.2` == max(c(`liking_S1F2.2`, `liking_S2F2.2`), na.rm = TRUE) ~ "high"
    ),
    
    # Sweet preference overall (at F2) - using averaged liking scores
    sweetPref = case_when(
      `liking_S1F2` == max(c(`liking_S1F2`, `liking_S2F2`), na.rm = TRUE) ~ "low",
      `liking_S2F2` == max(c(`liking_S1F2`, `liking_S2F2`), na.rm = TRUE) ~ "high"
    )
  ) %>%
  ungroup()

pref.df <- pref.df %>%
  mutate(
    # Relevel fat preferences (low, medium, high)
    fatPref = factor(fatPref, levels = c("low", "medium", "high")),
    fatPref.1 = factor(fatPref.1, levels = c("low", "medium", "high")),
    fatPref.2 = factor(fatPref.2, levels = c("low", "medium", "high")),
    
    # Relevel sweet preferences (low, high)
    sweetPref = factor(sweetPref, levels = c("low", "high")),
    sweetPref.1 = factor(sweetPref.1, levels = c("low", "high")),
    sweetPref.2 = factor(sweetPref.2, levels = c("low", "high"))
  )

Liking clusters

# Step 1: Prepare PCA input data
pca.df <- full.df %>%
  select(ID, liking_S0F0, liking_S1F0, liking_S2F0, liking_S1F1, liking_S1F2, liking_S1F3, liking_S2F1, liking_S2F2, liking_S2F3) %>%
  drop_na()

# Step 2: Run PCA
res.pca <- pca.df %>%
  select(-ID) %>%
  prcomp(., scale. = TRUE)

# Step 3: Get PCA coordinates
pca_coords <- res.pca$x

# Step 4: K-means clustering
km.res <- kmeans(pca_coords, centers = 3, nstart = 25)

    # Elbow and silhouette for K-means
    fviz_nbclust(pca_coords, kmeans, method = "wss")

    fviz_nbclust(pca_coords, kmeans, method = "silhouette")

# Step 5: Hierarchical clustering

dist_matrix <- dist(pca_coords)
hc.res <- hclust(dist_matrix, method = "ward.D2")
hc.groups <- cutree(hc.res, k = 3)

    # Plot dendrogram
    plot(hc.res, labels = FALSE, hang = -1, main = "Dendrogram")

# Step 6: Plot PCA biplot with K-means clusters
fviz_pca_biplot(res.pca,
                habillage = as.factor(km.res$cluster),
                addEllipses = TRUE,
                ellipse.level = 0.95,
                palette = "jco",
                repel = TRUE,
                label = "var") +
  ggtitle("PCA Biplot with K-means Clusters")

Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
ℹ Please use `linewidth` instead.
ℹ The deprecated feature was likely used in the ggpubr package.
  Please report the issue at <https://github.com/kassambara/ggpubr/issues>.

# Step 7: Plot PCA biplot with Hierarchical clusters
fviz_pca_biplot(res.pca,
                habillage = as.factor(hc.groups),
                addEllipses = TRUE,
                ellipse.level = 0.95,
                palette = "npg",
                repel = TRUE,
                label = "var") +
  ggtitle("PCA Biplot with Hierarchical Clusters")

## See who is in the cluster (hierarchical)
pca_coords$ID <- as.character(pca.df$ID)

Warning in pca_coords$ID <- as.character(pca.df$ID): Coercing LHS to a list

pca_coords$cluster <- as.factor(hc.groups)
pca_coords <- as.data.frame(pca_coords)

pca_coords %>%  
  group_by(cluster) %>%
  summarize(participants = paste(ID, collapse = " "))

# A tibble: 3 × 2
  cluster participants                                                          
  <fct>   <chr>                                                                 
1 1       3 8 7 11 18 19 21 31 36 38 39 50 54 61 62 68 76 81 101 114 122 130 13…
2 2       4 5 10 12 13 14 15 23 24 25 26 27 28 30 32 33 35 37 40 41 45 47 48 51…
3 3       6 9 16 17 20 22 29 34 42 43 44 46 49 56 57 69 70 80 88 90 96 99 102 1…

cluster <- pca_coords %>% select(ID, cluster)

cluster %>%
  group_by(cluster) %>%
  count()

# A tibble: 3 × 2
# Groups:   cluster [3]
  cluster     n
  <fct>   <int>
1 1          30
2 2         113
3 3          37

## Joining cluster with oat milk liking df 
full.df <- full.df  %>%
    left_join(cluster, by = "ID") 

oatmilk.liking.long$ID <- as.factor(oatmilk.liking.long$ID)

oatmilk.liking.long <- oatmilk.liking.long %>%  
  left_join(cluster, by = "ID") 

oatmilk.liking.wide <- oatmilk.liking.wide %>%  
  left_join(cluster, by = "ID") 

Fat liking slopes

## Calculate the slopes for the low sweet samples 
slopes.s1.df <- oatmilk.liking.long %>%
  filter(sample %in% c("S0F0", "S1F1", "S1F2", "S1F3")) %>%
  mutate(fat.level = case_when(
    str_detect(sample, "S0F0") ~ 0, 
    str_detect(sample, "S1F1") ~ 0.33,
    str_detect(sample, "S1F2") ~ 0.67,
    str_detect(sample, "S1F3") ~ 1.0
  )) 

slopes.s1.df <- slopes.s1.df %>%
  filter(fat.level == c(0, 0.33, 0.67, 1.0)) %>%
  filter(!is.na(liking)) %>%
  group_by(ID) %>%
  nest() %>%
  mutate(model = map(data, ~lm(liking ~ fat.level, data = .x)),
         slope = map_dbl(model, ~coef(.)[2])) %>%
  rename(fatliking.s1 = slope) %>%
  select(c(ID, fatliking.s1)) 

slopes.s1.df <- slopes.s1.df  %>%
  mutate(fatPhen.s1 = case_when(
    fatliking.s1 <= -15 ~ 1,
    fatliking.s1 > -15 & fatliking.s1 < 15 ~ 2,
    fatliking.s1 >= 15 ~ 3 )) 

## Calculate the slopes for the high sweet samples 
slopes.s2.df <- oatmilk.liking.long %>%
  filter(sample %in% c("S0F0", "S2F1", "S2F2", "S2F3")) %>%
  mutate(fat.level = case_when(
    str_detect(sample, "S0F0") ~ 0, 
    str_detect(sample, "S2F1") ~ 0.33,
    str_detect(sample, "S2F2") ~ 0.67,
    str_detect(sample, "S2F3") ~ 1.0
  )) 

slopes.s2.df <- slopes.s2.df %>%
  filter(fat.level == c(0, 0.33, 0.67, 1.0)) %>%
  filter(!is.na(liking)) %>%
  group_by(ID) %>%
  nest() %>%
  mutate(model = map(data, ~lm(liking ~ fat.level, data = .x)),
         slope = map_dbl(model, ~coef(.)[2])) %>%
  rename(fatliking.s2 = slope) %>%
  select(c(ID, fatliking.s2)) 

slopes.s2.df <- slopes.s2.df  %>%
  mutate(fatPhen.s2 = case_when(
    fatliking.s2 <= -15 ~ 1,
    fatliking.s2 > -15 & fatliking.s2 < 15 ~ 2,
    fatliking.s2 >= 15 ~ 3 )) 

### Joining the dfs to create a slopes dataframe
fatlikingslopes.df <- slopes.s1.df  %>%
  select(c("ID", "fatliking.s1", "fatPhen.s1")) %>%
  left_join(slopes.s2.df, by = "ID") 

fatlikingslopes.df$ID <- as.factor(fatlikingslopes.df$ID)


## Joining fating liking slopes df with oat milk liking df 
full.df <- full.df  %>%
    left_join(fatlikingslopes.df, by = "ID") 

oatmilk.liking.long$ID <- as.factor(oatmilk.liking.long$ID)

oatmilk.liking.long <- oatmilk.liking.long %>%  
  left_join(fatlikingslopes.df, by = "ID")


oatmilk.liking.wide <- oatmilk.liking.wide %>%  
  left_join(fatlikingslopes.df, by = "ID")

Sweet liking slopes

## Calculate the slopes for the fat-free samples 
slopes.f0.df <- oatmilk.liking.long %>%
  filter(sample %in% c("S0F0", "S1F0", "S2F0")) %>%
  mutate(sugar.level = case_when(
    str_detect(sample, "S0F0") ~ 0,
    str_detect(sample, "S1F0") ~ 0.5,
    str_detect(sample, "S2F0") ~ 1.0,
  )) 

slopes.f0.df <- slopes.f0.df %>%
  filter(sugar.level == c(0, 0.5, 1.0)) %>%
  filter(!is.na(liking)) %>%
  group_by(ID) %>%
  nest() %>%
  mutate(model = map(data, ~lm(liking ~ sugar.level, data = .x)),
         slope = map_dbl(model, ~coef(.)[2])) %>%
  rename(sweetliking.f0 = slope) %>%
  select(c(ID, sweetliking.f0)) 

slopes.f0.df <- slopes.f0.df  %>%
  mutate(sweetPhen.f0 = case_when(
    sweetliking.f0 <= -15 ~ 1,
    sweetliking.f0 > -15 & sweetliking.f0 < 15 ~ 2,
    sweetliking.f0 >= 15 ~ 3 )) 

## Calculate the slopes for the low fat samples 
slopes.f1.df <- oatmilk.liking.long %>%
  filter(sample %in% c("S0F0", "S1F1", "S2F1")) %>%
  mutate(sugar.level = case_when(
    str_detect(sample, "S0F0") ~ 0,
    str_detect(sample, "S1F1") ~ 0.5,
    str_detect(sample, "S2F1") ~ 1.0,
  )) 

slopes.f1.df <- slopes.f1.df %>%
  filter(sugar.level == c(0, 0.5, 1.0)) %>%
  filter(!is.na(liking)) %>%
  group_by(ID) %>%
  nest() %>%
  mutate(model = map(data, ~lm(liking ~ sugar.level, data = .x)),
         slope = map_dbl(model, ~coef(.)[2])) %>%
  rename(sweetliking.f1 = slope) %>%
  select(c(ID, sweetliking.f1)) 

slopes.f1.df <- slopes.f1.df  %>%
  mutate(sweetPhen.f1 = case_when(
    sweetliking.f1 <= -15 ~ 1,
    sweetliking.f1 > -15 & sweetliking.f1 < 15 ~ 2,
    sweetliking.f1 >= 15 ~ 3 )) 

## Calculate the slopes for the middle fat samples 
slopes.f2.df <- oatmilk.liking.long %>%
  filter(sample %in% c("S0F0", "S1F2", "S2F2")) %>%
  mutate(sugar.level = case_when(
    str_detect(sample, "S0F0") ~ 0, 
    str_detect(sample, "S1F2") ~ 0.5,
    str_detect(sample, "S2F2") ~ 1.,
  )) 

slopes.f2.df <- slopes.f2.df %>%
  filter(sugar.level == c(0, 0.5, 1.0)) %>%
  filter(!is.na(liking)) %>%
  group_by(ID) %>%
  nest() %>%
  mutate(model = map(data, ~lm(liking ~ sugar.level, data = .x)),
         slope = map_dbl(model, ~coef(.)[2])) %>%
  rename(sweetliking.f2 = slope) %>%
  select(c(ID, sweetliking.f2)) 

slopes.f2.df <- slopes.f2.df  %>%
  mutate(sweetPhen.f2 = case_when(
    sweetliking.f2 <= -15 ~ 1,
    sweetliking.f2 > -15 & sweetliking.f2 < 15 ~ 2,
    sweetliking.f2 >= 15 ~ 3 )) 

## Calculate the slopes for the high fat samples 
slopes.f3.df <- oatmilk.liking.long %>%
  filter(sample %in% c("S0F0", "S1F3", "S2F3")) %>%
  mutate(sugar.level = case_when(
    str_detect(sample, "S0F0") ~ 0,
    str_detect(sample, "S1F3") ~ 0.5,
    str_detect(sample, "S2F3") ~ 1.0,
  )) 

slopes.f3.df <- slopes.f3.df %>%
  filter(sugar.level == c(0, 0.5, 1.0)) %>%
  filter(!is.na(liking)) %>%
  group_by(ID) %>%
  nest() %>%
  mutate(model = map(data, ~lm(liking ~ sugar.level, data = .x)),
         slope = map_dbl(model, ~coef(.)[2])) %>%
  rename(sweetliking.f3 = slope) %>%
  select(c(ID, sweetliking.f3)) 

slopes.f3.df <- slopes.f3.df  %>%
  mutate(sweetPhen.f3 = case_when(
    sweetliking.f3 <= -15 ~ 1,
    sweetliking.f3 > -15 & sweetliking.f3 < 15 ~ 2,
    sweetliking.f3 >= 15 ~ 3 )) 


### Joining the dfs to create a slopes dataframe
sweetlikingslopes.df <- slopes.f0.df  %>%
  select(c("ID", "sweetliking.f0", "sweetPhen.f0")) %>%
    left_join(slopes.f1.df, by = "ID") %>%
    left_join(slopes.f2.df, by = "ID") %>%
    left_join(slopes.f3.df, by = "ID")

sweetlikingslopes.df$ID <- as.factor(sweetlikingslopes.df$ID)

## Joining sweet liking slopes with oat milk liking df 
full.df <- full.df  %>%
    left_join(sweetlikingslopes.df, by = "ID") 

oatmilk.liking.long$ID <- as.factor(oatmilk.liking.long$ID)

oatmilk.liking.long <- oatmilk.liking.long %>%  
  left_join(sweetlikingslopes.df, by ="ID") 

oatmilk.liking.wide <- oatmilk.liking.wide %>%  
  left_join(sweetlikingslopes.df, by ="ID") 

Participants flow, between measures - Sankey diagram

## Sweet preference flow: sweetPref → sweetPref.1 → sweetPref.2
sweet_flow <- pref.df %>%
  select(ID, sweetPref, sweetPref.1, sweetPref.2) %>%
  count(sweetPref, sweetPref.1, sweetPref.2) %>%
  filter(!is.na(sweetPref) & !is.na(sweetPref.1) & !is.na(sweetPref.2))

sweet_flow_graph <- ggplot(sweet_flow, 
       aes(axis1 = sweetPref, axis2 = sweetPref.1, axis3 = sweetPref.2, y = n)) +
  geom_alluvium(aes(fill = sweetPref), width = 1/12, alpha = 0.7) +
  geom_stratum(width = 1/5, fill = "white", color = "grey") +
  geom_text(stat = "stratum", aes(label = after_stat(stratum)), angle = 90, size = 4, fontface = "bold") +
  scale_x_discrete(limits = c("Overall", "Rep 1", "Rep 2"), expand = c(.05, .05)) +
  scale_fill_manual(values = c("low" = "#FFCA28", "high" = "#F44336")) +
  labs(title = " ",
       y = "Number of Participants",
       fill = "Sweet Preference") +
  theme_minimal() +
  theme(legend.position = "none")

print(sweet_flow_graph)

Warning in to_lodes_form(data = data, axes = axis_ind, discern =
params$discern): Some strata appear at multiple axes.
Warning in to_lodes_form(data = data, axes = axis_ind, discern =
params$discern): Some strata appear at multiple axes.
Warning in to_lodes_form(data = data, axes = axis_ind, discern =
params$discern): Some strata appear at multiple axes.

ggsave("sweet_flow_graph.png", height = 4, width = 5)

Warning in to_lodes_form(data = data, axes = axis_ind, discern =
params$discern): Some strata appear at multiple axes.
Warning in to_lodes_form(data = data, axes = axis_ind, discern =
params$discern): Some strata appear at multiple axes.
Warning in to_lodes_form(data = data, axes = axis_ind, discern =
params$discern): Some strata appear at multiple axes.

# Sweet preference stability
pref.df %>%
  filter(!is.na(sweetPref) & !is.na(sweetPref.1) & !is.na(sweetPref.2)) %>%
  filter(sweetPref == sweetPref.1 & sweetPref.1 == sweetPref.2) %>%
  count(sweetPref)

# A tibble: 2 × 2
  sweetPref     n
  <fct>     <int>
1 low          33
2 high        111

pref.df %>%
  filter(!is.na(sweetPref) & !is.na(sweetPref.1) & !is.na(sweetPref.2)) %>%
  mutate(stable = sweetPref == sweetPref.1 & sweetPref.1 == sweetPref.2) %>%
  summarise(n_stable = sum(stable),
            total = n(),
            percent_stable = round(n_stable / total * 100, 1))

# A tibble: 1 × 3
  n_stable total percent_stable
     <int> <int>          <dbl>
1      144   180             80

## Fat preference flow: fatPref → fatPref.1 → fatPref.2
fat_flow <- pref.df %>%
  select(ID, fatPref, fatPref.1, fatPref.2) %>%
  count(fatPref, fatPref.1, fatPref.2) %>%
  filter(!is.na(fatPref) & !is.na(fatPref.1) & !is.na(fatPref.2))

fat_flow_graph <- ggplot(fat_flow, 
       aes(axis1 = fatPref, axis2 = fatPref.1, axis3 = fatPref.2, y = n)) +
  geom_alluvium(aes(fill = fatPref), width = 1/12, alpha = 0.7) +
  geom_stratum(width = 1/5, fill = "white", color = "grey") +
  geom_text(stat = "stratum", aes(label = after_stat(stratum)), angle = 90, size = 4, fontface = "bold") +
  scale_x_discrete(limits = c("Overall", "Rep 1", "Rep 2"), expand = c(.05, .05)) +
  scale_fill_manual(values = c("low" = "#ffb300", "medium" = "#00897b", "high" = "#e53935")) +
  labs(title = " ",
       y = "Number of Participants",
       fill = "Fat Preference") +
  theme_minimal() +
  theme(legend.position = "none")

print(fat_flow_graph)

Warning in to_lodes_form(data = data, axes = axis_ind, discern =
params$discern): Some strata appear at multiple axes.
Warning in to_lodes_form(data = data, axes = axis_ind, discern =
params$discern): Some strata appear at multiple axes.
Warning in to_lodes_form(data = data, axes = axis_ind, discern =
params$discern): Some strata appear at multiple axes.

ggsave("fat_flow_graph.png", height = 4, width = 5)

Warning in to_lodes_form(data = data, axes = axis_ind, discern =
params$discern): Some strata appear at multiple axes.
Warning in to_lodes_form(data = data, axes = axis_ind, discern =
params$discern): Some strata appear at multiple axes.
Warning in to_lodes_form(data = data, axes = axis_ind, discern =
params$discern): Some strata appear at multiple axes.

# Combine side by side
combined_flow_graph <- sweet_flow_graph + fat_flow_graph +
  plot_annotation(title = " ")

print(combined_flow_graph)

Warning in to_lodes_form(data = data, axes = axis_ind, discern =
params$discern): Some strata appear at multiple axes.
Warning in to_lodes_form(data = data, axes = axis_ind, discern =
params$discern): Some strata appear at multiple axes.
Warning in to_lodes_form(data = data, axes = axis_ind, discern =
params$discern): Some strata appear at multiple axes.
Warning in to_lodes_form(data = data, axes = axis_ind, discern =
params$discern): Some strata appear at multiple axes.
Warning in to_lodes_form(data = data, axes = axis_ind, discern =
params$discern): Some strata appear at multiple axes.
Warning in to_lodes_form(data = data, axes = axis_ind, discern =
params$discern): Some strata appear at multiple axes.

ggsave("combined_flow_graph.png", plot = combined_flow_graph, height = 4, width = 8)

Warning in to_lodes_form(data = data, axes = axis_ind, discern =
params$discern): Some strata appear at multiple axes.
Warning in to_lodes_form(data = data, axes = axis_ind, discern =
params$discern): Some strata appear at multiple axes.
Warning in to_lodes_form(data = data, axes = axis_ind, discern =
params$discern): Some strata appear at multiple axes.
Warning in to_lodes_form(data = data, axes = axis_ind, discern =
params$discern): Some strata appear at multiple axes.
Warning in to_lodes_form(data = data, axes = axis_ind, discern =
params$discern): Some strata appear at multiple axes.
Warning in to_lodes_form(data = data, axes = axis_ind, discern =
params$discern): Some strata appear at multiple axes.

# Fat preference stability
pref.df %>%
  filter(!is.na(fatPref) & !is.na(fatPref.1) & !is.na(fatPref.2)) %>%
  filter(fatPref == fatPref.1 & fatPref.1 == fatPref.2) %>%
  count(fatPref)

# A tibble: 3 × 2
  fatPref     n
  <fct>   <int>
1 low        14
2 medium     20
3 high       27

pref.df %>%
  filter(!is.na(fatPref) & !is.na(fatPref.1) & !is.na(fatPref.2)) %>%
  mutate(stable = fatPref == fatPref.1 & fatPref.1 == fatPref.2) %>%
  summarise(n_stable = sum(stable),
            total = n(),
            percent_stable = round(n_stable / total * 100, 1))

# A tibble: 1 × 3
  n_stable total percent_stable
     <int> <int>          <dbl>
1       61   180           33.9

Participants flow, across sweet/fat levels - Sankey diagram

## Fat preference flow: fatPref → fatPref.1 → fatPref.2
fat_flow <- full.df %>%
  select(ID, fatpref.atS1, fatpref.atS2) %>%
  filter(!is.na(fatpref.atS1) & !is.na(fatpref.atS2)) %>%
  mutate(
    fatpref.atS1 = case_when(
      fatpref.atS1 == 1 ~ "Low",
      fatpref.atS1 == 2 ~ "Medium",
      fatpref.atS1 == 3 ~ "High"
    ),
    fatpref.atS2 = case_when(
      fatpref.atS2 == 1 ~ "Low",
      fatpref.atS2 == 2 ~ "Medium",
      fatpref.atS2 == 3 ~ "High"
    )
  ) %>%
  count(fatpref.atS1, fatpref.atS2) %>%
    mutate(
    fatpref.atS1 = factor(fatpref.atS1, levels = c("Low", "Medium", "High")),
    fatpref.atS2 = factor(fatpref.atS2, levels = c("Low", "Medium", "High"))
  )

ggplot(fat_flow, 
       aes(axis1 = fatpref.atS1, axis2 = fatpref.atS2, y = n)) +
  geom_alluvium(aes(fill = fatpref.atS1), width = 1/12, alpha = 0.7) +
  geom_stratum(width = 1/5, fill = "white", color = "grey") +
  geom_text(stat = "stratum", aes(label = after_stat(stratum)), size = 4, fontface = "bold") +
  scale_x_discrete(limits = c("Low Sweet", "High Sweet"), expand = c(.05, .05)) +
  scale_fill_manual(values = c("Low" = "#ffb300", "Medium" = "#00897b", "High" = "#e53935")) +
  labs(title = "Fat Preference Flow: Low Sweet → High Sweet",
       y = "Number of Participants",
       fill = "Fat Preference") +
  theme_minimal() +
  theme(legend.position = "right")

Warning in to_lodes_form(data = data, axes = axis_ind, discern =
params$discern): Some strata appear at multiple axes.
Warning in to_lodes_form(data = data, axes = axis_ind, discern =
params$discern): Some strata appear at multiple axes.
Warning in to_lodes_form(data = data, axes = axis_ind, discern =
params$discern): Some strata appear at multiple axes.

## Sweet preference flow: sweetPref → sweetPref.1 → sweetPref.2
sweet_flow <- full.df %>%
  select(ID, sweetpref.atF1, sweetpref.atF2, sweetpref.atF3) %>%
  filter(!is.na(sweetpref.atF1) & !is.na(sweetpref.atF2) & !is.na(sweetpref.atF3)) %>%
  mutate(
    sweetpref.atF1 = case_when(
      sweetpref.atF1 == 1 ~ "Low",
      sweetpref.atF1 == 2 ~ "High",
    ),
    sweetpref.atF2 = case_when(
      sweetpref.atF2 == 1 ~ "Low",
      sweetpref.atF2 == 2 ~ "High",
    ),
      sweetpref.atF3 = case_when(
      sweetpref.atF3 == 1 ~ "Low",
      sweetpref.atF3 == 2 ~ "High",
    )
  ) %>%
  count(sweetpref.atF1, sweetpref.atF2, sweetpref.atF3) %>%
    mutate(
    sweetpref.atF1 = factor(sweetpref.atF1, levels = c("Low", "High")),
    sweetpref.atF2 = factor(sweetpref.atF2, levels = c("Low", "High")),
    sweetpref.atF3 = factor(sweetpref.atF3, levels = c("Low", "High"))
  )

ggplot(sweet_flow, 
       aes(axis1 = sweetpref.atF1, axis2 = sweetpref.atF2, axis3 = sweetpref.atF3, y = n)) +
  geom_alluvium(aes(fill = sweetpref.atF1), width = 1/12, alpha = 0.7) +
  geom_stratum(width = 1/5, fill = "white", color = "grey") +
  geom_text(stat = "stratum", aes(label = after_stat(stratum)), size = 4, fontface = "bold") +
  scale_x_discrete(limits = c("Low Fat", "Medium Fat", "High Fat"), expand = c(.05, .05)) +
  scale_fill_manual(values = c("Low" = "#FFCA28", "High" = "#F44336")) +
  labs(title = "Sweet Preference Flow: Low Fat → Medium Fat → High Fat",
       y = "Number of Participants",
       fill = "Sweet Preference") +
  theme_minimal() +
  theme(legend.position = "right")

Warning in to_lodes_form(data = data, axes = axis_ind, discern =
params$discern): Some strata appear at multiple axes.
Warning in to_lodes_form(data = data, axes = axis_ind, discern =
params$discern): Some strata appear at multiple axes.
Warning in to_lodes_form(data = data, axes = axis_ind, discern =
params$discern): Some strata appear at multiple axes.

Dose response - liking

## Make sweet liking dose dataframe
sweet.liking.dose <- oatmilk.liking.long %>%
  group_by(sample, sweetPref) %>%
  summarize(mean = mean(liking), 
            sd = sd(liking),
            n = n(),           
            sem = sd/sqrt(n)) 

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

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

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

sweet.liking.dose$sweetPref <- factor(sweet.liking.dose$sweetPref , levels = c("Low", "High"))

## sweet liking by sweet preference
sweet.liking.dose %>%
  filter(!sugar.level %in% "S0") %>%
  filter(!fat.level %in% "F0") %>%
  ggplot(aes(x = sugar.level, y = mean, group = sweetPref, color = sweetPref)) +
  geom_line()+
  geom_hline(yintercept = 50, linetype = "dashed", color = "gray50") + 
  geom_point() +
  scale_color_manual(values = c("#00897b",  "#F44336"), name = "Sweet Preference") +
  geom_errorbar(aes(ymin = mean - sem, ymax = mean + sem), width = 0.2) +
  labs(y = "Liking") +
  facet_wrap(~fat.level) +
  ylim(0, 100) +
  theme_classic() 

## Make fat liking dose dataframe
fat.liking.dose <- oatmilk.liking.long %>%
  group_by(sample, fatPref) %>%
  summarize(mean = mean(liking), 
            sd = sd(liking),
            n = n(),           
            sem = sd/sqrt(n)) 

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

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

fat.liking.dose$fatPref <- factor(fat.liking.dose$fatPref , levels = c("Low", "Medium", "High"))

## fat liking by fat preference
fat.liking.dose %>%
  filter(!sugar.level %in% "S0") %>%
  filter(!fat.level %in% "F0") %>%
  ggplot(aes(x = fat.level, y = mean, group = fatPref, color = fatPref)) +
  geom_line()+
  geom_hline(yintercept = 50, linetype = "dashed", color = "gray50") + 
  geom_point() +
  scale_color_manual(values = c("#00897b", "#ffb300",  "#F44336"), name = "Fat Preference") +
  geom_errorbar(aes(ymin = mean - sem, ymax = mean + sem), width = 0.2) +
  labs(y = "Liking") +
  facet_wrap(~sugar.level) +
  ylim(0, 100) +
  theme_classic() 

## By individual
oatmilk.liking.long %>%
  group_by(ID, sample) %>%
    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() +
  labs(y = "Liking") +
  scale_color_manual(values = c("#FFCA28", "#F44336")) +
  facet_wrap(~ID)+
  theme_classic()

Dose response, by sweet preference

## Make sweet liking dose dataframe for sweet intensity
sweetPref.sweet.dose <- oatmilk.sweet.long %>%
  separate(sample, into = c("sugar.level", "fat.level"), sep = 2) %>%
  select(-fatPref) %>%
  group_by(sugar.level, sweetPref) %>%
  summarize(mean = mean(sweet), 
            sd = sd(sweet) ,
            n = n(),           
            sem = sd/sqrt(n)) 

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

sweetPref.sweet.dose$sugar.level <- factor(sweetPref.sweet.dose$sugar.level, levels = c("S0", "S1", "S2"))
sweetPref.sweet.dose$sweetPref <- factor(sweetPref.sweet.dose$sweetPref, levels = c("Low","High"))

## Testing the difference in sweet intensity perception between sweet liking phenotypes
sweetint.sweetPref.aov <- oatmilk.sweet.long %>%
    separate(sample, into = c("sugar.level", "fat.level"), sep = 2) %>%
  aov(sweet ~ sweetPref*sugar.level, .)
summary(sweetint.sweetPref.aov)

                        Df Sum Sq Mean Sq F value  Pr(>F)    
sweetPref                1   6270    6270  16.872 4.2e-05 ***
sugar.level              2 583372  291686 784.938 < 2e-16 ***
sweetPref:sugar.level    2   1630     815   2.193   0.112    
Residuals             1614 599768     372                    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

sweetint.sweetPref.emm <- emmeans(sweetint.sweetPref.aov, ~ sweetPref | sugar.level)
sweetint.sweetPref.pairwise <- pairs(sweetint.sweetPref.emm, adjust = "tukey") %>%
  as.data.frame()

# Prepare significance annotations
sig_data <- sweetint.sweetPref.pairwise %>%
  mutate(p.signif = case_when(p.value < 0.001 ~ "***", p.value < 0.01 ~ "**", p.value < 0.05 ~ "*", TRUE ~ "ns")) %>%
  filter(p.value < 0.05)  # Keep only significant comparisons

# Get max y value for positioning asterisks
max_y <- max(sweetPref.sweet.dose$mean + sweetPref.sweet.dose$sem, na.rm = TRUE)

# Add y position for each sugar level
sig_data <- sig_data %>%
  mutate(y.position = max_y + 5)  # Adjust height as needed

## Sweetness intensity dose response by sweetness preference
sweetPref.sweet <- sweetPref.sweet.dose %>%
  ggplot(aes(x = sugar.level, y = mean, group = sweetPref, color = sweetPref, shape = sweetPref)) +
  geom_line()+
  geom_point(size = 0.75) +
  scale_color_manual(values = c("#00897b", "#e53935"), name = "Sweet Preference", labels = c("Low", "High")) +
  scale_shape_manual(values = c(16, 15), name = "Sweet Preference", labels = c("Low", "High")) +  
  geom_errorbar(aes(ymin = mean - sem, ymax = mean + sem), width = 0.1) +
    geom_text(data = sig_data, aes(x = sugar.level, y = y.position, label = p.signif), color = "black", size = 5, inherit.aes = FALSE) +
    scale_x_discrete(labels = c("0%", "14%", "28%")) +
  labs(y = "Sweetness intensity (0 - 100)", x = "Sugar Concentration") +
  ylim(0, 100) +
  theme_classic()

print(sweetPref.sweet)

ggsave("sweetPref.sweet.png", height = 3, width = 4)

## Make sweet liking dose dataframe for creaminess intensity
sweetPref.creamy.dose <- oatmilk.creamy.long %>%
  separate(sample, into = c("sugar.level", "fat.level"), sep = 2) %>%
  select(-fatPref) %>%
  group_by(fat.level, sweetPref) %>%
  summarize(mean = mean(creamy), 
            sd = sd(creamy) ,
            n = n(),           
            sem = sd/sqrt(n)) 

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

sweetPref.creamy.dose$fat.level <- factor(sweetPref.creamy.dose$fat.level, levels = c("F0", "F1", "F2", "F3"))
sweetPref.creamy.dose$sweetPref <- factor(sweetPref.creamy.dose$sweetPref, levels = c("Low","High"))

## Testing the difference in sweet intensity perception between sweet liking phenotypes
creamyint.sweetPref.aov <- oatmilk.creamy.long %>%
    separate(sample, into = c("sugar.level", "fat.level"), sep = 2) %>%
  aov(creamy ~ sweetPref*fat.level, .)
summary(creamyint.sweetPref.aov)

                      Df Sum Sq Mean Sq F value Pr(>F)    
sweetPref              1    232     232   1.008  0.316    
fat.level              3 208473   69491 301.893 <2e-16 ***
sweetPref:fat.level    3    123      41   0.178  0.912    
Residuals           1612 371057     230                   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

creamyint.sweetPref.emm <- emmeans(creamyint.sweetPref.aov, ~ sweetPref | fat.level)
creamyint.sweetPref.pairwise <- pairs(creamyint.sweetPref.emm, adjust = "tukey") %>%
  as.data.frame()

# Prepare significance annotations
sig_data <- creamyint.sweetPref.pairwise %>%
  mutate(p.signif = case_when(p.value < 0.001 ~ "***", p.value < 0.01 ~ "**", p.value < 0.05 ~ "*", TRUE ~ "ns")) %>%
  filter(p.value < 0.05)  # Keep only significant comparisons

# Get max y value for positioning asterisks
max_y <- max(sweetPref.creamy.dose$mean + sweetPref.creamy.dose$sem, na.rm = TRUE)

# Add y position for each sugar level
sig_data <- sig_data %>%
  mutate(y.position = max_y + 5)  # Adjust height as needed

## Sweetness intensity dose response by sweetness preference
sweetPref.creamy <- sweetPref.creamy.dose %>%
  ggplot(aes(x = fat.level, y = mean, group = sweetPref, color = sweetPref, shape = sweetPref)) +
  geom_line()+
  geom_point(size = 0.75) +
  scale_color_manual(values = c("#00897b", "#e53935"),name = "Sweet Preference", labels = c("Low", "High")) +
  scale_shape_manual(values = c(16, 15),name = "Sweet Preference", labels = c("Low", "High")) +   geom_errorbar(aes(ymin = mean - sem, ymax = mean + sem), width = 0.1) +
    geom_text(data = sig_data, aes(x = fat.level, y = y.position, label = p.signif), color = "black", size = 5, inherit.aes = FALSE) +
    scale_x_discrete(labels = c("0%", "2%", "6%", "12%")) +
  labs(y = "Creaminess intensity (0 - 100)", x = "Fat Concentration") +
  ylim(0, 100) +
  theme_classic()

print(sweetPref.creamy)

ggsave("sweetPref.creamy.png", height = 3, width = 4)

Dose response, by fat preference

## Make fat liking dose dataframe for sweet intensity
fatPref.sweet.dose <- oatmilk.sweet.long %>%
  separate(sample, into = c("sugar.level", "fat.level"), sep = 2) %>%
  select(-sweetPref) %>%
  group_by(sugar.level, fatPref) %>%
  summarize(mean = mean(sweet), 
            sd = sd(sweet),
            n = n(),           
            sem = sd/sqrt(n)) 

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

fatPref.sweet.dose$sugar.level <- factor(fatPref.sweet.dose$sugar.level, levels = c("S0", "S1", "S2"))
fatPref.sweet.dose$fatPref <- factor(fatPref.sweet.dose$fatPref, levels = c("Low", "Medium", "High"))

## Testing the difference in sweet intensity perception between fat liking phenotypes
sweetint.fatPref.aov <- oatmilk.sweet.long %>%
  separate(sample, into = c("sugar.level", "fat.level"), sep = 2) %>%
  aov(sweet ~ fatPref*sugar.level, .)

summary(sweetint.fatPref.aov)

                      Df Sum Sq Mean Sq F value Pr(>F)    
fatPref                2   1557     779   2.071  0.126    
sugar.level            2 583372  291686 775.676 <2e-16 ***
fatPref:sugar.level    4    308      77   0.205  0.936    
Residuals           1611 605801     376                   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

sweetint.fatPref.emm <- emmeans(sweetint.fatPref.aov, ~ fatPref | sugar.level)
sweetint.fatPref.pairwise <- pairs(sweetint.fatPref.emm, adjust = "tukey") %>%
  as.data.frame()

# Prepare significance annotations
sig_data <- sweetint.fatPref.pairwise %>%
  mutate(p.signif = case_when(p.value < 0.001 ~ "***", 
                               p.value < 0.01 ~ "**", 
                               p.value < 0.05 ~ "*", 
                               TRUE ~ "ns")) %>%
  filter(p.value < 0.05)  # Keep only significant comparisons

# Get max y value for positioning asterisks
max_y <- max(fatPref.sweet.dose$mean + fatPref.sweet.dose$sem, na.rm = TRUE)

# Add y position for each fat level
sig_data <- sig_data %>%
  mutate(y.position = max_y + 5)  # Adjust height as needed

## Fattiness intensity dose response by fattiness preference
fatPref.sweet <- fatPref.sweet.dose %>%
  ggplot(aes(x = sugar.level, y = mean, group = fatPref, color = fatPref, shape = fatPref)) +
  geom_line() +
  geom_point(size = 0.75) +
  scale_color_manual(values = c("#00897b", "#ffb300", "#e53935"), name = "Fat Preference", labels = c("Low", "Medium", "High")) +
  scale_shape_manual(values = c(16, 17, 15), name = "Fat Preference", labels = c("Low", "Medium", "High")) +  
  geom_errorbar(aes(ymin = mean - sem, ymax = mean + sem), width = 0.1) +
  geom_text(data = sig_data, aes(x = sugar.level, y = y.position, label = p.signif), 
            color = "black", size = 5, inherit.aes = FALSE) +
  labs(y = "Sweetness intensity (0 - 100)",
       x = "Sugar Concentration") +
    scale_x_discrete(labels = c("0%", "14%", "28")) +
  ylim(0, 100) +
  theme_classic()

print(fatPref.sweet)

ggsave("fatPref.sweet.png", height = 3, width = 4)

## Make fat liking dose dataframe for creamy intensity
fatPref.creamy.dose <- oatmilk.creamy.long %>%
  separate(sample, into = c("sugar.level", "fat.level"), sep = 2) %>%
  select(-sweetPref) %>%
  group_by(fat.level, fatPref) %>%
  summarize(mean = mean(creamy), 
            sd = sd(creamy),
            n = n(),           
            sem = sd/sqrt(n)) 

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

fatPref.creamy.dose$fat.level <- factor(fatPref.creamy.dose$fat.level, levels = c("F0", "F1", "F2", "F3"))
fatPref.creamy.dose$fatPref <- factor(fatPref.creamy.dose$fatPref, levels = c("Low", "Medium", "High"))

## Testing the difference in fat intensity perception between fat liking phenotypes
fatint.fatPref.aov <- oatmilk.creamy.long %>%
  separate(sample, into = c("sugar.level", "fat.level"), sep = 2) %>%
  aov(creamy ~ fatPref*fat.level, .)

summary(fatint.fatPref.aov)

                    Df Sum Sq Mean Sq F value   Pr(>F)    
fatPref              2   3916    1958   8.640 0.000185 ***
fat.level            3 208473   69491 306.659  < 2e-16 ***
fatPref:fat.level    6   3111     519   2.288 0.033311 *  
Residuals         1608 364385     227                     
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

fatint.fatPref.emm <- emmeans(fatint.fatPref.aov, ~ fatPref | fat.level)
fatint.fatPref.pairwise <- pairs(fatint.fatPref.emm, adjust = "tukey") %>%
  as.data.frame()

# Prepare significance annotations
sig_data <- fatint.fatPref.pairwise %>%
  mutate(p.signif = case_when(p.value < 0.001 ~ "***", 
                               p.value < 0.01 ~ "**", 
                               p.value < 0.05 ~ "*", 
                               TRUE ~ "ns")) %>%
  filter(p.value < 0.05)  # Keep only significant comparisons

# Get max y value for positioning asterisks
max_y <- max(fatPref.creamy.dose$mean + fatPref.creamy.dose$sem, na.rm = TRUE)

# Add y position for each fat level
sig_data <- sig_data %>%
  mutate(y.position = max_y + 5)  # Adjust height as needed

## Fattiness intensity dose response by fattiness preference
fatPref.creamy <- fatPref.creamy.dose %>%
  ggplot(aes(x = fat.level, y = mean, group = fatPref, color = fatPref, shape = fatPref)) +
  geom_line() +
  geom_point(size = 0.75) +
  scale_color_manual(values = c("#00897b", "#ffb300", "#e53935"), name = "Fat Preference", labels = c("Low", "Medium", "High")) +
  scale_shape_manual(values = c(16, 17, 15), name = "Fat Preference", labels = c("Low", "Medium", "High")) +    
  geom_errorbar(aes(ymin = mean - sem, ymax = mean + sem), width = 0.1) +
  geom_text(data = sig_data, aes(x = fat.level, y = y.position, label = p.signif), 
            color = "black", size = 5, inherit.aes = FALSE) +
  scale_x_discrete(labels = c("0%", "2%", "6%", "12%")) +
  labs(y = "Creaminess intensity (0 - 100)",
       x = "Fat Concentration") +
  ylim(0, 100) +
  theme_classic()

print(fatPref.creamy)

ggsave("fatPref.creamy.png", height = 3, width = 4)

Dose response - bitter

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

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 = mean, group = sugar.level, color = sugar.level, shape = sugar.level)) +
  geom_line()+
  geom_point() +
  scale_color_manual(values = c("#FFCA28", "#F44336")) +
  geom_errorbar(aes(ymin = mean - sem, ymax = mean + sem), width = 0.05) +
  labs(y = "Bitter") +
  ylim(-10, 100) +
  theme_classic()

## By individual
oatmilk.bitter.long %>%
  group_by(ID, sample) %>%
    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()

Corrplot - all

## 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)
}

## survey corr
    ## other variables to potentially add to the corrplot (liking_S0F0, liking_S1F0, liking_S1F1, liking_S1F2, liking_S1F3, liking_S2F0, liking_S2F1, liking_S2F2, liking_S2F3, sweetliking.f0, sweetliking.f1, sweetliking.f2, sweetliking.f3)

all.corr.df <- full.df %>%
  select(age, BMI, WC, hinc, restraint, emotional, uncontrolled, sfbl.liking, plantmilk_freq, sugar.intake, sHEI, lean_mass, SMM, FFM, SMM_wt, trunk_lean_mass, BFM, PBF, trunk_fat_mass, visceral_fat, BMR, sweetpref.atF2, sweetpref.atF3, fatpref.atS1, pref.cal) %>%
  mutate(across(everything(), as.numeric)) %>%
  rename("household income" = hinc, "sugar intake" = sugar.intake, 
         "healthy eating index" = sHEI, "sweet foods liking" = sfbl.liking) 

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

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

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

No id variables; using all as measure variables

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

No id variables; using all as measure variables

Partial correlation

# Step 1: Prepare data
all.corr.df <- full.df %>%
  select(age, plantmilk_freq, BMI, WC, hinc, restraint, emotional, uncontrolled, 
         sfbl.liking, healthyfat.liking, unhealthyfat.liking, sugarHEI, fatHEI, sHEI, lean_mass, SMM, FFM, SMM_wt, 
         trunk_lean_mass, BFM, PBF, trunk_fat_mass, visceral_fat, BMR, sweetpref.atF2, fatpref.atS1, pref.cal) %>%
  mutate(across(everything(), as.numeric)) %>%
  rename("Household Income" = hinc, "Sugar Intake Score" = sugarHEI, "Fatty Acid Intake Score" = fatHEI, "Healthy Eating Index" = sHEI, "Sweet Foods Liking" = sfbl.liking, "Healthy Fat Liking" = healthyfat.liking, "Unhealthy Fat Liking" = unhealthyfat.liking, "Sweet Preference" = sweetpref.atF2, "Fat Preference" = fatpref.atS1, 
         "Restrictive Eating" = restraint, "Emotional Eating" = emotional, "Uncontrolled Eating" = uncontrolled,
         "Waist Circumference" = WC ,"Lean Mass" = lean_mass, "Skeletal Muscle Mass" = SMM, "Seletal Muscle Mass, Weight Adjusted" = SMM_wt, "Fat-Free Mass" = FFM,  "Trunk Lean Mass" = trunk_lean_mass, "Body Fat Mass" = BFM, "Body Fat Percentage" = PBF, "Trunk Fat Mass" = trunk_fat_mass, "Visceral Fat" = visceral_fat, "Basal Metabolic Rate" = BMR)

# Variables of interest (excluding both covariates)
vars_of_interest <- c("Sweet Preference", "Fat Preference", 
                      "Sweet Foods Liking", "Healthy Fat Liking", "Unhealthy Fat Liking",
                      "Restrictive Eating", "Emotional Eating", "Uncontrolled Eating", 
                      "Sugar Intake Score", "Fatty Acid Intake Score", "Healthy Eating Index", 
                      "BMI", "Waist Circumference", "Lean Mass", "Skeletal Muscle Mass", "Seletal Muscle Mass, Weight Adjusted", "Fat-Free Mass", "Trunk Lean Mass", "Body Fat Mass", "Body Fat Percentage", "Trunk Fat Mass", "Visceral Fat", "Basal Metabolic Rate")

# Combine with BOTH control variables
all_vars <- c(vars_of_interest, "age", "plantmilk_freq")
data_subset <- all.corr.df[, all_vars]

# Compute correlation matrix with pairwise deletion
cor_matrix <- cor(data_subset, use = "pairwise.complete.obs")

# Compute partial correlations from the correlation matrix
partial_cor_full <- cor2pcor(cor_matrix)

# Extract just the correlations among your variables of interest
# (exclude the last 2 rows/columns which are age and plantmilk_freq)
n_vars <- length(vars_of_interest)
partial_cor_matrix <- partial_cor_full[1:n_vars, 1:n_vars]

# Add variable names
rownames(partial_cor_matrix) <- vars_of_interest
colnames(partial_cor_matrix) <- vars_of_interest

# Updated function to compute p-values with TWO control variables
compute_partial_cor_pvalues_2covs <- function(data, vars_interest, control_vars) {
  n_vars <- length(vars_interest)
  p_matrix <- matrix(1, nrow = n_vars, ncol = n_vars)
  
  for (i in 1:(n_vars-1)) {
    for (j in (i+1):n_vars) {
      # Get complete cases for this specific set of variables (2 vars + 2 controls)
      subset_data <- data[, c(vars_interest[i], vars_interest[j], control_vars)]
      subset_complete <- na.omit(subset_data)
      
      if (nrow(subset_complete) > 4) {  # Need at least 5 observations for 2 covariates
        # Compute partial correlation for this pair, controlling for BOTH covariates
        pcor_result <- pcor.test(subset_complete[,1], subset_complete[,2], 
                                 subset_complete[, 3:4])  # Both control variables
        p_matrix[i, j] <- pcor_result$p.value
        p_matrix[j, i] <- pcor_result$p.value
      }
    }
  }
  
  return(p_matrix)
}

# Compute p-values (this may take a moment)
p_values <- compute_partial_cor_pvalues_2covs(data_subset, vars_of_interest, 
                                               c("age", "plantmilk_freq"))
rownames(p_values) <- vars_of_interest
colnames(p_values) <- vars_of_interest

# Apply FDR correction
lower_tri_indices <- lower.tri(p_values)
p_values_vector <- p_values[lower_tri_indices]
p_adjusted <- p.adjust(p_values_vector, method = "fdr")

# Create adjusted p-value matrix
p_adjusted_matrix <- matrix(1, nrow = n_vars, ncol = n_vars)
p_adjusted_matrix[lower_tri_indices] <- p_adjusted
p_adjusted_matrix <- t(p_adjusted_matrix)
p_adjusted_matrix[lower_tri_indices] <- p_adjusted
diag(p_adjusted_matrix) <- 0

rownames(p_adjusted_matrix) <- vars_of_interest
colnames(p_adjusted_matrix) <- vars_of_interest

# Plot
corr <- ggcorrplot(partial_cor_matrix, 
           type = "lower", 
           p.mat = p_adjusted_matrix, 
           insig = "blank", 
           lab = TRUE, 
           lab_size = 2.5, 
           colors = c("#E46726", "white", "#6D9EC1"),
           title = "")

print(corr)

ggsave("corr.png", height = 10, width = 10)

Partial correlation - Julia version

# Step 1: Prepare data
all.corr.df <- full.df %>%
  select(age, plantmilk_freq, BMI, WC, hinc, sweetpref.atF2, fatpref.atS1, sugarHEI, fatHEI, sHEI) %>%
  mutate(across(everything(), as.numeric)) %>%
  rename("Sugar Intake Score" = sugarHEI, "Fatty Acid Intake Score" = fatHEI, "Healthy Eating Index" = sHEI, "Sweet Preference" = sweetpref.atF2, "Fat Preference" = fatpref.atS1, "Waist Circumference" = WC )

# Variables of interest (excluding both covariates)
vars_of_interest <- c("Sweet Preference", "Fat Preference", 
                      "Sugar Intake Score", "Fatty Acid Intake Score", "Healthy Eating Index", 
                      "BMI", "Waist Circumference")

# Combine with BOTH control variables
all_vars <- c(vars_of_interest, "age", "plantmilk_freq")
data_subset <- all.corr.df[, all_vars]

# Compute correlation matrix with pairwise deletion
cor_matrix <- cor(data_subset, use = "pairwise.complete.obs")

# Compute partial correlations from the correlation matrix
partial_cor_full <- cor2pcor(cor_matrix)

# Extract just the correlations among your variables of interest
# (exclude the last 2 rows/columns which are age and plantmilk_freq)
n_vars <- length(vars_of_interest)
partial_cor_matrix <- partial_cor_full[1:n_vars, 1:n_vars]

# Add variable names
rownames(partial_cor_matrix) <- vars_of_interest
colnames(partial_cor_matrix) <- vars_of_interest

# Updated function to compute p-values with TWO control variables
compute_partial_cor_pvalues_2covs <- function(data, vars_interest, control_vars) {
  n_vars <- length(vars_interest)
  p_matrix <- matrix(1, nrow = n_vars, ncol = n_vars)
  
  for (i in 1:(n_vars-1)) {
    for (j in (i+1):n_vars) {
      # Get complete cases for this specific set of variables (2 vars + 2 controls)
      subset_data <- data[, c(vars_interest[i], vars_interest[j], control_vars)]
      subset_complete <- na.omit(subset_data)
      
      if (nrow(subset_complete) > 4) {  # Need at least 5 observations for 2 covariates
        # Compute partial correlation for this pair, controlling for BOTH covariates
        pcor_result <- pcor.test(subset_complete[,1], subset_complete[,2], 
                                 subset_complete[, 3:4])  # Both control variables
        p_matrix[i, j] <- pcor_result$p.value
        p_matrix[j, i] <- pcor_result$p.value
      }
    }
  }
  
  return(p_matrix)
}

# Compute p-values (this may take a moment)
p_values <- compute_partial_cor_pvalues_2covs(data_subset, vars_of_interest, 
                                               c("age", "plantmilk_freq"))
rownames(p_values) <- vars_of_interest
colnames(p_values) <- vars_of_interest

# Apply FDR correction
lower_tri_indices <- lower.tri(p_values)
p_values_vector <- p_values[lower_tri_indices]
p_adjusted <- p.adjust(p_values_vector, method = "fdr")

# Create adjusted p-value matrix
p_adjusted_matrix <- matrix(1, nrow = n_vars, ncol = n_vars)
p_adjusted_matrix[lower_tri_indices] <- p_adjusted
p_adjusted_matrix <- t(p_adjusted_matrix)
p_adjusted_matrix[lower_tri_indices] <- p_adjusted
diag(p_adjusted_matrix) <- 0

rownames(p_adjusted_matrix) <- vars_of_interest
colnames(p_adjusted_matrix) <- vars_of_interest

# Plot
corr <- ggcorrplot(partial_cor_matrix, 
           type = "lower", 
           p.mat = p_adjusted_matrix, 
           insig = "blank", 
           lab = TRUE, 
           lab_size = 2.5, 
           colors = c("#E46726", "white", "#6D9EC1"),
           title = "")

print(corr)

ggsave("julia.corr.png", height = 5, width = 5)

Participant characteristics

full.df %>%
  summarise(agemean = mean(age), 
            agesd = sd(age),
            n = n(),           
            agesem = agesd/sqrt(n))

# A tibble: 1 × 4
  agemean agesd     n agesem
    <dbl> <dbl> <int>  <dbl>
1    21.5  3.22   180  0.240

full.df %>%
  count(race)

# A tibble: 12 × 2
   race      n
   <chr> <int>
 1 1        46
 2 1,2       1
 3 1,2,4     2
 4 1,3,4     1
 5 1,4       2
 6 1,7       1
 7 2         1
 8 3        61
 9 3,4       1
10 4        41
11 5        10
12 7        13

full.df %>%
  count(ethn)

# A tibble: 2 × 2
   ethn     n
  <dbl> <int>
1     1    33
2     2   147

full.df %>%
  summarise(BMImean = mean(BMI), 
            BMIsd = sd(BMI),
            n = n(),           
            BMIsem = BMIsd/sqrt(n))

# A tibble: 1 × 4
  BMImean BMIsd     n BMIsem
    <dbl> <dbl> <int>  <dbl>
1    25.1  6.66   180  0.497

full.df %>%
  count(BMIcat)

# A tibble: 4 × 2
  BMIcat     n
  <chr>  <int>
1 normal    99
2 obese     25
3 over      43
4 under     13

full.df %>%
  count(sweetPref)

# A tibble: 2 × 2
  sweetPref     n
  <fct>     <int>
1 Low          40
2 High        140

full.df %>%
  count(fatPref)

# A tibble: 3 × 2
  fatPref     n
  <fct>   <int>
1 Low        46
2 Medium     70
3 High       64

full.df %>%
  count(overall.pref)

# A tibble: 9 × 2
  overall.pref     n
  <chr>        <int>
1 S0F0            22
2 S1F0            15
3 S1F1            11
4 S1F2             9
5 S1F3            15
6 S2F0            10
7 S2F1            29
8 S2F2            39
9 S2F3            30

full.df %>%
  group_by(sfblCat) %>%
  summarize(meanLiking = mean(sfbl.liking), sdLiking = sd(sfbl.liking))

# A tibble: 3 × 3
  sfblCat meanLiking sdLiking
  <fct>        <dbl>    <dbl>
1 Low           56.0     8.80
2 Medium        69.8     3.14
3 High          82.2     4.70

full.df %>%
  group_by(unhealthyfatCat) %>%
  summarize(meanLiking = mean(unhealthyfat.liking), sdLiking = sd(unhealthyfat.liking))

# A tibble: 3 × 3
  unhealthyfatCat meanLiking sdLiking
  <fct>                <dbl>    <dbl>
1 Low                   49.8     8.13
2 Medium                63.3     2.81
3 High                  75.8     6.22

full.df %>%
  group_by(healthyfatCat) %>%
  summarize(meanLiking = mean(healthyfat.liking), sdLiking = sd(healthyfat.liking))

# A tibble: 3 × 3
  healthyfatCat meanLiking sdLiking
  <fct>              <dbl>    <dbl>
1 Low                 53.4     7.97
2 Medium              68.5     3.47
3 High                82.8     6.11

Group difference, 1-way ANOVA demographics

## Participant counts
table(full.df$sweetPref)

 Low High 
  40  140 

chisq.test(table(full.df$sweetPref))

    Chi-squared test for given probabilities

data:  table(full.df$sweetPref)
X-squared = 55.556, df = 1, p-value = 9.085e-14

table(full.df$fatPref)

   Low Medium   High 
    46     70     64 

chisq.test(table(full.df$fatPref))

    Chi-squared test for given probabilities

data:  table(full.df$fatPref)
X-squared = 5.2, df = 2, p-value = 0.07427

## Age
sweet.age.aov <- full.df %>%
  aov(age ~ sweetPref , .)
summary(sweet.age.aov)

             Df Sum Sq Mean Sq F value Pr(>F)
sweetPref     1    2.5    2.54   0.243  0.622
Residuals   178 1858.4   10.44               

full.df %>% 
  group_by(sweetPref) %>%
  summarize(mean = mean(age), sd = sd(age) , n = n(), sem = sd/sqrt(n))

# A tibble: 2 × 5
  sweetPref  mean    sd     n   sem
  <fct>     <dbl> <dbl> <int> <dbl>
1 Low        21.3  2.44    40 0.386
2 High       21.6  3.42   140 0.289

fat.age.aov <- full.df %>%
  aov(age ~ fatPref , .)
summary(fat.age.aov)

             Df Sum Sq Mean Sq F value Pr(>F)
fatPref       2      5   2.499   0.238  0.788
Residuals   177   1856  10.485               

full.df %>% 
  group_by(fatPref) %>%
  summarize(mean = mean(age), sd = sd(age) , n = n(), sem = sd/sqrt(n))

# A tibble: 3 × 5
  fatPref  mean    sd     n   sem
  <fct>   <dbl> <dbl> <int> <dbl>
1 Low      21.6  2.53    46 0.373
2 Medium   21.3  3.14    70 0.375
3 High     21.7  3.75    64 0.469

## BMI
sweet.BMI.aov <- full.df %>%
  aov(BMI ~ sweetPref , .)
summary(sweet.BMI.aov)

             Df Sum Sq Mean Sq F value Pr(>F)  
sweetPref     1    133  133.31   3.037 0.0831 .
Residuals   178   7813   43.89                 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

full.df %>% 
  group_by(sweetPref) %>%
  summarize(mean = mean(BMI), sd = sd(BMI) , n = n(), sem = sd/sqrt(n))

# A tibble: 2 × 5
  sweetPref  mean    sd     n   sem
  <fct>     <dbl> <dbl> <int> <dbl>
1 Low        23.5  4.80    40 0.759
2 High       25.6  7.05   140 0.596

fat.BMI.aov <- full.df %>%
  aov(BMI ~ fatPref , .)
summary(fat.BMI.aov)

             Df Sum Sq Mean Sq F value Pr(>F)
fatPref       2    174   87.11   1.984  0.141
Residuals   177   7772   43.91               

full.df %>% 
  group_by(fatPref) %>%
  summarize(mean = mean(BMI), sd = sd(BMI) , n = n(), sem = sd/sqrt(n))

# A tibble: 3 × 5
  fatPref  mean    sd     n   sem
  <fct>   <dbl> <dbl> <int> <dbl>
1 Low      26.7  8.05    46 1.19 
2 Medium   24.2  5.24    70 0.627
3 High     25.0  6.86    64 0.857

## WC
sweet.WC.aov <- full.df %>%
  aov(WC ~ sweetPref , .)
summary(sweet.WC.aov)

             Df Sum Sq Mean Sq F value Pr(>F)
sweetPref     1     49   48.84   1.649  0.201
Residuals   178   5273   29.62               

full.df %>% 
  group_by(sweetPref) %>%
  summarize(mean = mean(WC), sd = sd(WC) , n = n(), sem = sd/sqrt(n))

# A tibble: 2 × 5
  sweetPref  mean    sd     n   sem
  <fct>     <dbl> <dbl> <int> <dbl>
1 Low        30.0  4.57    40 0.723
2 High       31.3  5.66   140 0.479

fat.WC.aov <- full.df %>%
  aov(WC ~ fatPref , .)
summary(fat.WC.aov)

             Df Sum Sq Mean Sq F value Pr(>F)
fatPref       2     84   41.93   1.417  0.245
Residuals   177   5238   29.59               

full.df %>% 
  group_by(fatPref) %>%
  summarize(mean = mean(WC), sd = sd(WC) , n = n(), sem = sd/sqrt(n))

# A tibble: 3 × 5
  fatPref  mean    sd     n   sem
  <fct>   <dbl> <dbl> <int> <dbl>
1 Low      32.1  6.55    46 0.965
2 Medium   30.4  4.66    70 0.557
3 High     30.8  5.37    64 0.671

## Education
table(full.df$sweetPref, full.df$education)

      
        4  5  6  8  9 10
  Low   0  7 22  1  8  2
  High  1 24 57 18 36  4

chisq.test(full.df$sweetPref, full.df$education)

Warning in stats::chisq.test(x, y, ...): Chi-squared approximation may be
incorrect

    Pearson's Chi-squared test

data:  full.df$sweetPref and full.df$education
X-squared = 5.7405, df = 5, p-value = 0.3323

table(full.df$fatPref, full.df$education)

        
          4  5  6  8  9 10
  Low     0  5 23  8  9  1
  Medium  1 13 33  3 16  4
  High    0 13 23  8 19  1

chisq.test(full.df$fatPref, full.df$education)

Warning in stats::chisq.test(x, y, ...): Chi-squared approximation may be
incorrect

    Pearson's Chi-squared test

data:  full.df$fatPref and full.df$education
X-squared = 12.651, df = 10, p-value = 0.2438

## Household income
table(full.df$sweetPref, full.df$hinc)

      
        1  2  3  4  5  6  7  8
  Low   2  4  6  9  8  3  4  4
  High  3  9 20 35 25 15 19 14

chisq.test(full.df$sweetPref, full.df$hinc)

Warning in stats::chisq.test(x, y, ...): Chi-squared approximation may be
incorrect

    Pearson's Chi-squared test

data:  full.df$sweetPref and full.df$hinc
X-squared = 2.2642, df = 7, p-value = 0.9438

table(full.df$fatPref, full.df$hinc)

        
          1  2  3  4  5  6  7  8
  Low     2  2  8  8  8  5  9  4
  Medium  2  5 12 20 10  7  6  8
  High    1  6  6 16 15  6  8  6

chisq.test(full.df$fatPref, full.df$hinc)

Warning in stats::chisq.test(x, y, ...): Chi-squared approximation may be
incorrect

    Pearson's Chi-squared test

data:  full.df$fatPref and full.df$hinc
X-squared = 9.3825, df = 14, p-value = 0.8058

## PB-milk frequency
sweet.pbmilk.aov <- full.df %>%
  aov(plantmilk_freq ~ sweetPref, .)
summary(sweet.pbmilk.aov)

             Df Sum Sq Mean Sq F value Pr(>F)  
sweetPref     1   16.7  16.677   5.469 0.0205 *
Residuals   178  542.8   3.049                 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

full.df %>% 
  group_by(sweetPref) %>%
  summarize(mean = mean(plantmilk_freq), sd = sd(plantmilk_freq) , n = n(), sem = sd/sqrt(n))

# A tibble: 2 × 5
  sweetPref  mean    sd     n   sem
  <fct>     <dbl> <dbl> <int> <dbl>
1 Low        3.62  1.85    40 0.292
2 High       2.89  1.72   140 0.145

fat.pbmilk.aov <- full.df %>%
  aov(plantmilk_freq ~ fatPref , .)
summary(fat.pbmilk.aov)

             Df Sum Sq Mean Sq F value  Pr(>F)   
fatPref       2   35.9  17.926    6.06 0.00285 **
Residuals   177  523.6   2.958                   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

emmeans(fat.pbmilk.aov, pairwise ~ fatPref, adjust = "tukey")

$emmeans
 fatPref emmean    SE  df lower.CL upper.CL
 Low       2.74 0.254 177     2.24     3.24
 Medium    2.71 0.206 177     2.31     3.12
 High      3.66 0.215 177     3.23     4.08

Confidence level used: 0.95 

$contrasts
 contrast      estimate    SE  df t.ratio p.value
 Low - Medium    0.0248 0.326 177   0.076  0.9968
 Low - High     -0.9171 0.332 177  -2.759  0.0175
 Medium - High  -0.9420 0.297 177  -3.167  0.0051

P value adjustment: tukey method for comparing a family of 3 estimates 

emmeans(fat.pbmilk.aov, ~fatPref)

 fatPref emmean    SE  df lower.CL upper.CL
 Low       2.74 0.254 177     2.24     3.24
 Medium    2.71 0.206 177     2.31     3.12
 High      3.66 0.215 177     3.23     4.08

Confidence level used: 0.95 

full.df %>% 
  group_by(fatPref) %>%
  summarize(mean = mean(plantmilk_freq), sd = sd(plantmilk_freq) , n = n(), sem = sd/sqrt(n))

# A tibble: 3 × 5
  fatPref  mean    sd     n   sem
  <fct>   <dbl> <dbl> <int> <dbl>
1 Low      2.74  1.73    46 0.255
2 Medium   2.71  1.64    70 0.196
3 High     3.66  1.79    64 0.224

## Emotional Eating
sweet.emo.aov <- full.df %>%
  aov(emotional ~ sweetPref , .)
summary(sweet.emo.aov)

             Df Sum Sq Mean Sq F value Pr(>F)
sweetPref     1     88    88.0   0.215  0.644
Residuals   178  72894   409.5               

full.df %>% 
  group_by(sweetPref) %>%
  summarize(mean = mean(emotional), sd = sd(emotional) , n = n(), sem = sd/sqrt(n))

# A tibble: 2 × 5
  sweetPref  mean    sd     n   sem
  <fct>     <dbl> <dbl> <int> <dbl>
1 Low        48.2  21.6    40  3.42
2 High       49.9  19.8   140  1.68

fat.emo.aov <- full.df %>%
  aov(emotional ~ fatPref , .)
summary(fat.emo.aov)

             Df Sum Sq Mean Sq F value Pr(>F)
fatPref       2    194    97.2   0.236   0.79
Residuals   177  72788   411.2               

full.df %>% 
  group_by(fatPref) %>%
  summarize(mean = mean(emotional), sd = sd(emotional) , n = n(), sem = sd/sqrt(n))

# A tibble: 3 × 5
  fatPref  mean    sd     n   sem
  <fct>   <dbl> <dbl> <int> <dbl>
1 Low      51.3  19.2    46  2.83
2 Medium   48.7  21.5    70  2.57
3 High     49.2  19.6    64  2.46

## Uncontrolled Eating
sweet.unc.aov <- full.df %>%
  aov(uncontrolled ~ sweetPref , .)
summary(sweet.unc.aov)

             Df Sum Sq Mean Sq F value Pr(>F)
sweetPref     1    191   191.4   1.118  0.292
Residuals   178  30459   171.1               

full.df %>% 
  group_by(sweetPref) %>%
  summarize(mean = mean(uncontrolled), sd = sd(uncontrolled) , n = n(), sem = sd/sqrt(n))

# A tibble: 2 × 5
  sweetPref  mean    sd     n   sem
  <fct>     <dbl> <dbl> <int> <dbl>
1 Low        50.1  11.1    40  1.76
2 High       52.6  13.6   140  1.15

fat.unc.aov <- full.df %>%
  aov(uncontrolled ~ fatPref , .)
summary(fat.unc.aov)

             Df Sum Sq Mean Sq F value Pr(>F)  
fatPref       2    850   425.1   2.525 0.0829 .
Residuals   177  29800   168.4                 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

full.df %>% 
  group_by(fatPref) %>%
  summarize(mean = mean(uncontrolled), sd = sd(uncontrolled) , n = n(), sem = sd/sqrt(n))

# A tibble: 3 × 5
  fatPref  mean    sd     n   sem
  <fct>   <dbl> <dbl> <int> <dbl>
1 Low      55.4  13.6    46  2.01
2 Medium   49.8  12.5    70  1.50
3 High     52.1  13.0    64  1.62

## Restrictive Eating
sweet.res.aov <- full.df %>%
  aov(restraint ~ sweetPref , .)
summary(sweet.res.aov)

             Df Sum Sq Mean Sq F value Pr(>F)
sweetPref     1    863   863.3   2.037  0.155
Residuals   178  75442   423.8               

full.df %>% 
  group_by(sweetPref) %>%
  summarize(mean = mean(restraint), sd = sd(restraint) , n = n(), sem = sd/sqrt(n))

# A tibble: 2 × 5
  sweetPref  mean    sd     n   sem
  <fct>     <dbl> <dbl> <int> <dbl>
1 Low        50.2  22.6    40  3.57
2 High       44.9  20.0   140  1.69

fat.res.aov <- full.df %>%
  aov(restraint ~ fatPref , .)
summary(fat.res.aov)

             Df Sum Sq Mean Sq F value Pr(>F)
fatPref       2    883   441.5   1.036  0.357
Residuals   177  75423   426.1               

full.df %>% 
  group_by(fatPref) %>%
  summarize(mean = mean(restraint), sd = sd(restraint) , n = n(), sem = sd/sqrt(n))

# A tibble: 3 × 5
  fatPref  mean    sd     n   sem
  <fct>   <dbl> <dbl> <int> <dbl>
1 Low      46.0  18.7    46  2.75
2 Medium   43.7  21.2    70  2.53
3 High     48.8  21.4    64  2.68

full.df %>% 
  group_by(fatPref) %>%
  summarize(mean = mean(PBF, na.rm = TRUE), 
            sd = sd(PBF, na.rm = TRUE), 
            n = n())

# A tibble: 3 × 4
  fatPref  mean    sd     n
  <fct>   <dbl> <dbl> <int>
1 Low      36.3  8.66    46
2 Medium   33.3  7.21    70
3 High     32.8  7.25    64

full.df %>% 
  group_by(fatPref) %>%
  summarize(mean = mean(visceral_fat, na.rm = TRUE), 
            sd = sd(visceral_fat, na.rm = TRUE), 
            n = n())

# A tibble: 3 × 4
  fatPref  mean    sd     n
  <fct>   <dbl> <dbl> <int>
1 Low     124.   68.3    46
2 Medium   97.1  50.5    70
3 High     97.2  49.0    64

Group difference, 2-way ANCOVA, taste on diet

## sHEI
sHEI.aov <- full.df %>%
  aov(sHEI ~ sweetPref*fatPref + age + plantmilk_freq, data = .)
summary(sHEI.aov)

                   Df Sum Sq Mean Sq F value  Pr(>F)   
sweetPref           1    708   707.9   9.770 0.00208 **
fatPref             2    686   342.9   4.733 0.00998 **
age                 1    125   125.1   1.727 0.19055   
plantmilk_freq      1     21    21.4   0.296 0.58736   
sweetPref:fatPref   2     52    26.0   0.359 0.69909   
Residuals         172  12462    72.5                   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

eta_squared(sHEI.aov)

        sweetPref           fatPref               age    plantmilk_freq 
      0.050368944       0.048797627       0.008903117       0.001523876 
sweetPref:fatPref 
      0.003698606 

emmeans(sHEI.aov, pairwise ~ sweetPref, adjust = "tukey")

NOTE: Results may be misleading due to involvement in interactions

$emmeans
 sweetPref emmean    SE  df lower.CL upper.CL
 Low         49.8 1.430 172     46.9     52.6
 High        44.2 0.737 172     42.7     45.7

Results are averaged over the levels of: fatPref 
Confidence level used: 0.95 

$contrasts
 contrast   estimate   SE  df t.ratio p.value
 Low - High     5.55 1.62 172   3.423  0.0008

Results are averaged over the levels of: fatPref 

sweet.emm <- emmeans(sHEI.aov, ~ sweetPref, adjust = "tukey") %>% as.data.frame()

NOTE: Results may be misleading due to involvement in interactions

Note: adjust = "tukey" was changed to "sidak"
because "tukey" is only appropriate for one set of pairwise comparisons

emmeans(sHEI.aov, ~sweetPref)

NOTE: Results may be misleading due to involvement in interactions

 sweetPref emmean    SE  df lower.CL upper.CL
 Low         49.8 1.430 172     46.9     52.6
 High        44.2 0.737 172     42.7     45.7

Results are averaged over the levels of: fatPref 
Confidence level used: 0.95 

emmeans(sHEI.aov, pairwise ~ fatPref, adjust = "tukey")

NOTE: Results may be misleading due to involvement in interactions

$emmeans
 fatPref emmean   SE  df lower.CL upper.CL
 Low       46.9 1.47 172     44.0     49.8
 Medium    45.0 1.15 172     42.8     47.3
 High      49.0 1.52 172     46.0     52.0

Results are averaged over the levels of: sweetPref 
Confidence level used: 0.95 

$contrasts
 contrast      estimate   SE  df t.ratio p.value
 Low - Medium      1.84 1.86 172   0.986  0.5865
 Low - High       -2.11 2.13 172  -0.989  0.5848
 Medium - High    -3.95 1.91 172  -2.062  0.1010

Results are averaged over the levels of: sweetPref 
P value adjustment: tukey method for comparing a family of 3 estimates 

fat.emm <- emmeans(sHEI.aov, ~ fatPref) %>% as.data.frame()

NOTE: Results may be misleading due to involvement in interactions

emmeans(sHEI.aov, ~fatPref)

NOTE: Results may be misleading due to involvement in interactions

 fatPref emmean   SE  df lower.CL upper.CL
 Low       46.9 1.47 172     44.0     49.8
 Medium    45.0 1.15 172     42.8     47.3
 High      49.0 1.52 172     46.0     52.0

Results are averaged over the levels of: sweetPref 
Confidence level used: 0.95 

## Graph

# Sweet preference significance
sweet.sig <- full.df %>%
  tukey_hsd(sHEI ~ sweetPref) %>%
  add_significance()

# Fat preference significance
fat.sig <- full.df %>%
  tukey_hsd(sHEI ~ fatPref) %>%
  add_significance()

# Sweet preference plot
sweet.plot <- ggplot(sweet.emm, aes(x = sweetPref, y = emmean, fill = sweetPref)) +
  geom_bar(stat = "identity", width = 0.5, alpha = 0.8) +
  geom_errorbar(aes(ymin = lower.CL, ymax = upper.CL), width = 0.2) +
  stat_pvalue_manual(sweet.sig, hide.ns = TRUE, label = "p.adj.signif",
                     label.size = 10, bracket.size = 1, tip.length = 0.05,
                     y.position = c(58)) +
  scale_fill_manual(values = c("Low" = "#00897b", "High" = "#e53935")) +  # low = teal, high = red
  labs(x = "Sweet Preference", y = "Healthy Eating Index (sHEI)", title = " ") +
  ylim(0, 70) +
  theme_classic() +
  theme(legend.position = "none",
        axis.text = element_text(size = 14),
        axis.title = element_text(size = 16))

# Fat preference plot
fat.plot <- ggplot(fat.emm, aes(x = fatPref, y = emmean, fill = fatPref)) +
  geom_bar(stat = "identity", width = 0.5, alpha = 0.8) +
  geom_errorbar(aes(ymin = lower.CL, ymax = upper.CL), width = 0.2) +
  stat_pvalue_manual(fat.sig, hide.ns = TRUE, label = "p.adj.signif",
                     label.size = 10, bracket.size = 1, tip.length = 0.05,
                     y.position = c(60)) +
  scale_fill_manual(values = c("Low" = "#00897b", "Medium" = "#ffb300", "High" = "#e53935")) +  # low = teal, medium = yellow, high = red
  labs(x = "Fat Preference", y = "Healthy Eating Index (sHEI)", title = " ") +
  ylim(0, 70) +
  theme_classic() +
  theme(legend.position = "none",
        axis.text = element_text(size = 14),
        axis.title = element_text(size = 16))

# Stitch together
sweet.plot + fat.plot

ggsave("taste.sHEI.png", height = 5, width = 6)


## Fat HEI
fatHEI.aov <- full.df %>%
  aov(fatHEI ~ sweetPref*fatPref + age + plantmilk_freq, .)
summary(fatHEI.aov)

                   Df Sum Sq Mean Sq F value Pr(>F)
sweetPref           1   3.59   3.586   2.052  0.154
fatPref             2   3.03   1.513   0.865  0.423
age                 1   2.22   2.221   1.271  0.261
plantmilk_freq      1   3.72   3.716   2.126  0.147
sweetPref:fatPref   2   1.41   0.703   0.402  0.670
Residuals         172 300.61   1.748               

summary.aov(fatHEI.aov)

                   Df Sum Sq Mean Sq F value Pr(>F)
sweetPref           1   3.59   3.586   2.052  0.154
fatPref             2   3.03   1.513   0.865  0.423
age                 1   2.22   2.221   1.271  0.261
plantmilk_freq      1   3.72   3.716   2.126  0.147
sweetPref:fatPref   2   1.41   0.703   0.402  0.670
Residuals         172 300.61   1.748               

eta_squared(fatHEI.aov)

        sweetPref           fatPref               age    plantmilk_freq 
      0.011399460       0.009617428       0.007061332       0.011812661 
sweetPref:fatPref 
      0.004467402 

emmeans(fatHEI.aov, pairwise ~ sweetPref, adjust = "tukey")

NOTE: Results may be misleading due to involvement in interactions

$emmeans
 sweetPref emmean    SE  df lower.CL upper.CL
 Low         3.56 0.222 172     3.12     4.00
 High        3.34 0.115 172     3.12     3.57

Results are averaged over the levels of: fatPref 
Confidence level used: 0.95 

$contrasts
 contrast   estimate    SE  df t.ratio p.value
 Low - High    0.216 0.252 172   0.859  0.3917

Results are averaged over the levels of: fatPref 

emmeans(fatHEI.aov, ~sweetPref)

NOTE: Results may be misleading due to involvement in interactions

 sweetPref emmean    SE  df lower.CL upper.CL
 Low         3.56 0.222 172     3.12     4.00
 High        3.34 0.115 172     3.12     3.57

Results are averaged over the levels of: fatPref 
Confidence level used: 0.95 

emmeans(fatHEI.aov, pairwise ~ fatPref, adjust = "tukey")

NOTE: Results may be misleading due to involvement in interactions

$emmeans
 fatPref emmean    SE  df lower.CL upper.CL
 Low       3.22 0.229 172     2.77     3.67
 Medium    3.63 0.178 172     3.28     3.99
 High      3.50 0.236 172     3.03     3.97

Results are averaged over the levels of: sweetPref 
Confidence level used: 0.95 

$contrasts
 contrast      estimate    SE  df t.ratio p.value
 Low - Medium    -0.416 0.289 172  -1.439  0.3235
 Low - High      -0.282 0.331 172  -0.850  0.6726
 Medium - High    0.135 0.297 172   0.454  0.8929

Results are averaged over the levels of: sweetPref 
P value adjustment: tukey method for comparing a family of 3 estimates 

emmeans(fatHEI.aov, ~fatPref)

NOTE: Results may be misleading due to involvement in interactions

 fatPref emmean    SE  df lower.CL upper.CL
 Low       3.22 0.229 172     2.77     3.67
 Medium    3.63 0.178 172     3.28     3.99
 High      3.50 0.236 172     3.03     3.97

Results are averaged over the levels of: sweetPref 
Confidence level used: 0.95 

## Graph

# Fat HEI significance
sweet.fatHEI.sig <- full.df %>%
  tukey_hsd(fatHEI ~ sweetPref) %>%
  add_significance()

fat.fatHEI.sig <- full.df %>%
  tukey_hsd(fatHEI ~ fatPref) %>%
  add_significance()

# Get emmeans
sweet.fatHEI.emm <- emmeans(fatHEI.aov, ~ sweetPref) %>% as.data.frame()

NOTE: Results may be misleading due to involvement in interactions

fat.fatHEI.emm <- emmeans(fatHEI.aov, ~ fatPref) %>% as.data.frame()

NOTE: Results may be misleading due to involvement in interactions

# Sweet preference plot
sweet.fatHEI.plot <- ggplot(sweet.fatHEI.emm, aes(x = sweetPref, y = emmean, fill = sweetPref)) +
  geom_bar(stat = "identity", width = 0.5, alpha = 0.8) +
  geom_errorbar(aes(ymin = lower.CL, ymax = upper.CL), width = 0.2) +
  scale_fill_manual(values = c("Low" = "#00897b", "High" = "#e53935")) +
  labs(x = "Sweet Preference", y = "Fat Intake Score", title = " ") +
  ylim(0, 5) +
  theme_classic() +
  theme(legend.position = "none",
        axis.text = element_text(size = 14),
        axis.title = element_text(size = 16))

# Fat preference plot
fat.fatHEI.plot <- ggplot(fat.fatHEI.emm, aes(x = fatPref, y = emmean, fill = fatPref)) +
  geom_bar(stat = "identity", width = 0.5, alpha = 0.8) +
  geom_errorbar(aes(ymin = lower.CL, ymax = upper.CL), width = 0.2) +
  scale_fill_manual(values = c("Low" = "#00897b", "Medium" = "#ffb300", "High" = "#e53935")) +
  labs(x = "Fat Preference", y = "Fat Intake Score", title = " ") +
  ylim(0, 5) +
  theme_classic() +
  theme(legend.position = "none",
        axis.text = element_text(size = 14),
        axis.title = element_text(size = 16))

# Stitch together
sweet.fatHEI.plot + fat.fatHEI.plot + plot_layout(axis_titles = "collect")

ggsave("taste.fatHEI.png", height = 5, width = 6)


## Sugar HEI
sugarHEI.aov <- full.df %>%
  aov(sugarHEI ~ sweetPref*fatPref + age + plantmilk_freq, .)
summary(sugarHEI.aov)

                   Df Sum Sq Mean Sq F value  Pr(>F)   
sweetPref           1   61.9   61.91   8.891 0.00328 **
fatPref             2    2.3    1.17   0.167 0.84606   
age                 1   11.6   11.63   1.671 0.19789   
plantmilk_freq      1   24.0   23.97   3.442 0.06526 . 
sweetPref:fatPref   2    7.2    3.62   0.520 0.59547   
Residuals         172 1197.8    6.96                   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

summary.aov(sugarHEI.aov)

                   Df Sum Sq Mean Sq F value  Pr(>F)   
sweetPref           1   61.9   61.91   8.891 0.00328 **
fatPref             2    2.3    1.17   0.167 0.84606   
age                 1   11.6   11.63   1.671 0.19789   
plantmilk_freq      1   24.0   23.97   3.442 0.06526 . 
sweetPref:fatPref   2    7.2    3.62   0.520 0.59547   
Residuals         172 1197.8    6.96                   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

eta_squared(sugarHEI.aov)

        sweetPref           fatPref               age    plantmilk_freq 
      0.047449251       0.001785972       0.008916406       0.018370496 
sweetPref:fatPref 
      0.005550011 

emmeans(sugarHEI.aov, pairwise ~ sweetPref, adjust = "tukey")

NOTE: Results may be misleading due to involvement in interactions

$emmeans
 sweetPref emmean    SE  df lower.CL upper.CL
 Low         5.39 0.442 172     4.52     6.27
 High        3.66 0.229 172     3.21     4.11

Results are averaged over the levels of: fatPref 
Confidence level used: 0.95 

$contrasts
 contrast   estimate    SE  df t.ratio p.value
 Low - High     1.73 0.502 172   3.446  0.0007

Results are averaged over the levels of: fatPref 

emmeans(sugarHEI.aov, ~sweetPref)

NOTE: Results may be misleading due to involvement in interactions

 sweetPref emmean    SE  df lower.CL upper.CL
 Low         5.39 0.442 172     4.52     6.27
 High        3.66 0.229 172     3.21     4.11

Results are averaged over the levels of: fatPref 
Confidence level used: 0.95 

emmeans(sugarHEI.aov, pairwise ~ fatPref, adjust = "tukey")

NOTE: Results may be misleading due to involvement in interactions

$emmeans
 fatPref emmean    SE  df lower.CL upper.CL
 Low       4.77 0.457 172     3.87     5.67
 Medium    4.21 0.355 172     3.51     4.91
 High      4.61 0.472 172     3.68     5.54

Results are averaged over the levels of: sweetPref 
Confidence level used: 0.95 

$contrasts
 contrast      estimate    SE  df t.ratio p.value
 Low - Medium     0.557 0.578 172   0.963  0.6011
 Low - High       0.159 0.661 172   0.240  0.9687
 Medium - High   -0.398 0.594 172  -0.670  0.7812

Results are averaged over the levels of: sweetPref 
P value adjustment: tukey method for comparing a family of 3 estimates 

emmeans(sugarHEI.aov, ~fatPref)

NOTE: Results may be misleading due to involvement in interactions

 fatPref emmean    SE  df lower.CL upper.CL
 Low       4.77 0.457 172     3.87     5.67
 Medium    4.21 0.355 172     3.51     4.91
 High      4.61 0.472 172     3.68     5.54

Results are averaged over the levels of: sweetPref 
Confidence level used: 0.95 

## Graph

# Sugar HEI significance
sweet.sugarHEI.sig <- full.df %>%
  tukey_hsd(sugarHEI ~ sweetPref) %>%
  add_significance()

fat.sugarHEI.sig <- full.df %>%
  tukey_hsd(sugarHEI ~ fatPref) %>%
  add_significance()

# Get emmeans
sweet.sugarHEI.emm <- emmeans(sugarHEI.aov, ~ sweetPref) %>% as.data.frame()

NOTE: Results may be misleading due to involvement in interactions

fat.sugarHEI.emm <- emmeans(sugarHEI.aov, ~ fatPref) %>% as.data.frame()

NOTE: Results may be misleading due to involvement in interactions

# Sweet preference plot
sweet.sugarHEI.plot <- ggplot(sweet.sugarHEI.emm, aes(x = sweetPref, y = emmean, fill = sweetPref)) +
  geom_bar(stat = "identity", width = 0.5, alpha = 0.8) +
  geom_errorbar(aes(ymin = lower.CL, ymax = upper.CL), width = 0.2) +
  stat_pvalue_manual(sweet.sugarHEI.sig, hide.ns = TRUE, label = "p.adj.signif",
                     label.size = 10, bracket.size = 1, tip.length = 0.05,
                     y.position = c(7.5)) +
  scale_fill_manual(values = c("Low" = "#00897b", "High" = "#e53935")) +
  labs(x = "Sweet Preference", y = "Sugar Intake Score", title = " ") +
  ylim(0, 8) +
  theme_classic() +
  theme(legend.position = "none",
        axis.text = element_text(size = 14),
        axis.title = element_text(size = 16))

# Fat preference plot
fat.sugarHEI.plot <- ggplot(fat.sugarHEI.emm, aes(x = fatPref, y = emmean, fill = fatPref)) +
  geom_bar(stat = "identity", width = 0.5, alpha = 0.8) +
  geom_errorbar(aes(ymin = lower.CL, ymax = upper.CL), width = 0.2) +
  scale_fill_manual(values = c("Low" = "#00897b", "Medium" = "#ffb300", "High" = "#e53935")) +
  labs(x = "Fat Preference", y = "Sugar Intake Score", title = " ") +
  ylim(0, 8) +
  theme_classic() +
  theme(legend.position = "none",
        axis.text = element_text(size = 14),
        axis.title = element_text(size = 16))

# Stitch together
sweet.sugarHEI.plot + fat.sugarHEI.plot + plot_layout(axis_titles = "collect")

ggsave("taste.sugarHEI.png", height = 5, width = 6)

## Sugar intake
sugarintake.aov <- full.df %>%
  aov(sugar.intake ~ sweetPref*fatPref + age + plantmilk_freq, .)
summary(sugarintake.aov)

                   Df   Sum Sq Mean Sq F value Pr(>F)  
sweetPref           1   194270  194270   2.797 0.0962 .
fatPref             2    42961   21480   0.309 0.7344  
age                 1        3       3   0.000 0.9945  
plantmilk_freq      1    96892   96892   1.395 0.2392  
sweetPref:fatPref   2    42487   21243   0.306 0.7369  
Residuals         172 11944385   69444                 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

summary.aov(sugarintake.aov)

                   Df   Sum Sq Mean Sq F value Pr(>F)  
sweetPref           1   194270  194270   2.797 0.0962 .
fatPref             2    42961   21480   0.309 0.7344  
age                 1        3       3   0.000 0.9945  
plantmilk_freq      1    96892   96892   1.395 0.2392  
sweetPref:fatPref   2    42487   21243   0.306 0.7369  
Residuals         172 11944385   69444                 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

eta_squared(sugarintake.aov)

        sweetPref           fatPref               age    plantmilk_freq 
     1.576738e-02      3.486779e-03      2.720832e-07      7.863963e-03 
sweetPref:fatPref 
     3.448333e-03 

emmeans(sugarintake.aov, pairwise ~ sweetPref, adjust = "tukey")

NOTE: Results may be misleading due to involvement in interactions

$emmeans
 sweetPref emmean   SE  df lower.CL upper.CL
 Low          400 44.2 172      313      487
 High         499 22.8 172      454      544

Results are averaged over the levels of: fatPref 
Confidence level used: 0.95 

$contrasts
 contrast   estimate   SE  df t.ratio p.value
 Low - High    -99.3 50.2 172  -1.980  0.0493

Results are averaged over the levels of: fatPref 

emmeans(sugarintake.aov, ~sweetPref)

NOTE: Results may be misleading due to involvement in interactions

 sweetPref emmean   SE  df lower.CL upper.CL
 Low          400 44.2 172      313      487
 High         499 22.8 172      454      544

Results are averaged over the levels of: fatPref 
Confidence level used: 0.95 

emmeans(sugarintake.aov, pairwise ~ fatPref, adjust = "tukey")

NOTE: Results may be misleading due to involvement in interactions

$emmeans
 fatPref emmean   SE  df lower.CL upper.CL
 Low        421 45.6 172      331      511
 Medium     482 35.5 172      412      552
 High       446 47.1 172      353      539

Results are averaged over the levels of: sweetPref 
Confidence level used: 0.95 

$contrasts
 contrast      estimate   SE  df t.ratio p.value
 Low - Medium     -61.5 57.7 172  -1.066  0.5364
 Low - High       -25.1 66.0 172  -0.381  0.9232
 Medium - High     36.4 59.3 172   0.614  0.8129

Results are averaged over the levels of: sweetPref 
P value adjustment: tukey method for comparing a family of 3 estimates 

emmeans(sugarintake.aov, ~fatPref)

NOTE: Results may be misleading due to involvement in interactions

 fatPref emmean   SE  df lower.CL upper.CL
 Low        421 45.6 172      331      511
 Medium     482 35.5 172      412      552
 High       446 47.1 172      353      539

Results are averaged over the levels of: sweetPref 
Confidence level used: 0.95 

Group difference, 3-way ANCOVA, liking on diet

## sHEI
sHEI.aov <- full.df %>%
  aov(sHEI ~ sfblCat*healthyfatCat*unhealthyfatCat + age + plantmilk_freq, data = .)
summary(sHEI.aov)

                                       Df Sum Sq Mean Sq F value   Pr(>F)    
sfblCat                                 2   1285   642.7   9.721 0.000107 ***
healthyfatCat                           2    940   470.0   7.108 0.001122 ** 
unhealthyfatCat                         2    274   136.9   2.070 0.129750    
age                                     1     27    27.0   0.408 0.524045    
plantmilk_freq                          1      0     0.3   0.004 0.948604    
sfblCat:healthyfatCat                   4    318    79.4   1.201 0.312881    
sfblCat:unhealthyfatCat                 4    412   103.1   1.559 0.188162    
healthyfatCat:unhealthyfatCat           4    188    47.1   0.712 0.584624    
sfblCat:healthyfatCat:unhealthyfatCat   8    626    78.2   1.183 0.312950    
Residuals                             151   9984    66.1                     
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

eta_squared(sHEI.aov)

                              sfblCat                         healthyfatCat 
                         9.146692e-02                          6.687799e-02 
                      unhealthyfatCat                                   age 
                         1.947678e-02                          1.918593e-03 
                       plantmilk_freq                 sfblCat:healthyfatCat 
                         1.961296e-05                          2.259638e-02 
              sfblCat:unhealthyfatCat         healthyfatCat:unhealthyfatCat 
                         2.933729e-02                          1.340711e-02 
sfblCat:healthyfatCat:unhealthyfatCat 
                         4.451822e-02 

emmeans(sHEI.aov, pairwise ~ sfblCat, adjust = "tukey")

NOTE: Results may be misleading due to involvement in interactions

$emmeans
 sfblCat emmean   SE  df lower.CL upper.CL
 Low       50.0 1.21 151     47.6     52.4
 Medium    43.6 1.10 151     41.4     45.8
 High      43.7 1.19 151     41.3     46.0

Results are averaged over the levels of: healthyfatCat, unhealthyfatCat 
Confidence level used: 0.95 

$contrasts
 contrast      estimate   SE  df t.ratio p.value
 Low - Medium    6.3866 1.65 151   3.868  0.0005
 Low - High      6.3088 1.70 151   3.708  0.0009
 Medium - High  -0.0778 1.61 151  -0.048  0.9987

Results are averaged over the levels of: healthyfatCat, unhealthyfatCat 
P value adjustment: tukey method for comparing a family of 3 estimates 

emmeans(sHEI.aov, ~sfblCat)

NOTE: Results may be misleading due to involvement in interactions

 sfblCat emmean   SE  df lower.CL upper.CL
 Low       50.0 1.21 151     47.6     52.4
 Medium    43.6 1.10 151     41.4     45.8
 High      43.7 1.19 151     41.3     46.0

Results are averaged over the levels of: healthyfatCat, unhealthyfatCat 
Confidence level used: 0.95 

emmeans(sHEI.aov, pairwise ~ healthyfatCat, adjust = "tukey")

NOTE: Results may be misleading due to involvement in interactions

$emmeans
 healthyfatCat emmean   SE  df lower.CL upper.CL
 Low             43.8 1.23 151     41.3     46.2
 Medium          44.4 1.17 151     42.1     46.7
 High            49.2 1.16 151     46.9     51.4

Results are averaged over the levels of: sfblCat, unhealthyfatCat 
Confidence level used: 0.95 

$contrasts
 contrast      estimate   SE  df t.ratio p.value
 Low - Medium    -0.614 1.69 151  -0.364  0.9295
 Low - High      -5.409 1.73 151  -3.128  0.0059
 Medium - High   -4.795 1.65 151  -2.906  0.0117

Results are averaged over the levels of: sfblCat, unhealthyfatCat 
P value adjustment: tukey method for comparing a family of 3 estimates 

emmeans(sHEI.aov, ~healthyfatCat)

NOTE: Results may be misleading due to involvement in interactions

 healthyfatCat emmean   SE  df lower.CL upper.CL
 Low             43.8 1.23 151     41.3     46.2
 Medium          44.4 1.17 151     42.1     46.7
 High            49.2 1.16 151     46.9     51.4

Results are averaged over the levels of: sfblCat, unhealthyfatCat 
Confidence level used: 0.95 

emmeans(sHEI.aov, pairwise ~ unhealthyfatCat, adjust = "tukey")

NOTE: Results may be misleading due to involvement in interactions

$emmeans
 unhealthyfatCat emmean   SE  df lower.CL upper.CL
 Low               47.3 1.18 151     44.9     49.6
 Medium            44.4 1.08 151     42.3     46.6
 High              45.6 1.25 151     43.1     48.0

Results are averaged over the levels of: sfblCat, healthyfatCat 
Confidence level used: 0.95 

$contrasts
 contrast      estimate   SE  df t.ratio p.value
 Low - Medium      2.84 1.60 151   1.773  0.1823
 Low - High        1.70 1.74 151   0.973  0.5952
 Medium - High    -1.14 1.65 151  -0.691  0.7688

Results are averaged over the levels of: sfblCat, healthyfatCat 
P value adjustment: tukey method for comparing a family of 3 estimates 

emmeans(sHEI.aov, ~healthyfatCat)

NOTE: Results may be misleading due to involvement in interactions

 healthyfatCat emmean   SE  df lower.CL upper.CL
 Low             43.8 1.23 151     41.3     46.2
 Medium          44.4 1.17 151     42.1     46.7
 High            49.2 1.16 151     46.9     51.4

Results are averaged over the levels of: sfblCat, unhealthyfatCat 
Confidence level used: 0.95 

# Significance tests
sfbl.sHEI.sig <- full.df %>%
  tukey_hsd(sHEI ~ sfblCat) %>%
  add_significance()

healthyfat.sHEI.sig <- full.df %>%
  tukey_hsd(sHEI ~ healthyfatCat) %>%
  add_significance()

unhealthyfat.sHEI.sig <- full.df %>%
  tukey_hsd(sHEI ~ unhealthyfatCat) %>%
  add_significance()

# Get emmeans
sfbl.sHEI.emm <- emmeans(sHEI.aov, ~ sfblCat) %>% as.data.frame()

NOTE: Results may be misleading due to involvement in interactions

healthyfat.sHEI.emm <- emmeans(sHEI.aov, ~ healthyfatCat) %>% as.data.frame()

NOTE: Results may be misleading due to involvement in interactions

unhealthyfat.sHEI.emm <- emmeans(sHEI.aov, ~ unhealthyfatCat) %>% as.data.frame()

NOTE: Results may be misleading due to involvement in interactions

# sfblCat plot
sfbl.sHEI.plot <- ggplot(sfbl.sHEI.emm, aes(x = sfblCat, y = emmean, fill = sfblCat)) +
  geom_bar(stat = "identity", width = 0.5, alpha = 0.8) +
  geom_errorbar(aes(ymin = lower.CL, ymax = upper.CL), width = 0.2) +
  stat_pvalue_manual(sfbl.sHEI.sig, hide.ns = TRUE, label = "p.adj.signif",
                     label.size = 10, bracket.size = 1, tip.length = 0.05,
                     y.position = c(58, 68)) +
  scale_fill_manual(values = c("Low" = "#00897b", "Medium" = "#ffb300", "High" = "#e53935")) +
  labs(x = "Sweet Fat Beverage Liking", y = "Healthy Eating Index (sHEI)", title = " ") +
  ylim(0, 75) +
  theme_classic() +
  theme(legend.position = "none",
        axis.text = element_text(size = 14),
        axis.title = element_text(size = 16))

# healthyfatCat plot
healthyfat.sHEI.plot <- ggplot(healthyfat.sHEI.emm, aes(x = healthyfatCat, y = emmean, fill = healthyfatCat)) +
  geom_bar(stat = "identity", width = 0.5, alpha = 0.8) +
  geom_errorbar(aes(ymin = lower.CL, ymax = upper.CL), width = 0.2) +
  stat_pvalue_manual(healthyfat.sHEI.sig, hide.ns = TRUE, label = "p.adj.signif",
                     label.size = 10, bracket.size = 1, tip.length = 0.05,
                     y.position = c(60)) +
  scale_fill_manual(values = c("Low" = "#00897b", "Medium" = "#ffb300", "High" = "#e53935")) +
  labs(x = "Healthy Fat Liking", y = " ", title = " ") +
  ylim(0, 75) +
  theme_classic() +
  theme(legend.position = "none",
        axis.text = element_text(size = 14),
        axis.title = element_text(size = 16))

# unhealthyfatCat plot
unhealthyfat.sHEI.plot <- ggplot(unhealthyfat.sHEI.emm, aes(x = unhealthyfatCat, y = emmean, fill = unhealthyfatCat)) +
  geom_bar(stat = "identity", width = 0.5, alpha = 0.8) +
  geom_errorbar(aes(ymin = lower.CL, ymax = upper.CL), width = 0.2) +
  scale_fill_manual(values = c("Low" = "#00897b", "Medium" = "#ffb300", "High" = "#e53935")) +
  labs(x = "Unhealthy Fat Liking", y = " ", title = " ") +
  ylim(0, 75) +
  theme_classic() +
  theme(legend.position = "none",
        axis.text = element_text(size = 14),
        axis.title = element_text(size = 16))

# Stitch together
sfbl.sHEI.plot + healthyfat.sHEI.plot + unhealthyfat.sHEI.plot +
  plot_layout(axis_titles = "collect")

ggsave("liking.sHEI.png", height = 5, width = 9)

## Fat HEI
fatHEI.aov <- full.df %>%
  aov(fatHEI ~ sfblCat*healthyfatCat*unhealthyfatCat + age + plantmilk_freq, .)
summary(fatHEI.aov)

                                       Df Sum Sq Mean Sq F value   Pr(>F)    
sfblCat                                 2   9.20   4.600   2.857 0.060557 .  
healthyfatCat                           2   6.02   3.009   1.869 0.157874    
unhealthyfatCat                         2  27.06  13.532   8.404 0.000346 ***
age                                     1   2.70   2.695   1.674 0.197694    
plantmilk_freq                          1   1.77   1.769   1.099 0.296187    
sfblCat:healthyfatCat                   4   2.14   0.534   0.332 0.856432    
sfblCat:unhealthyfatCat                 4   7.80   1.950   1.211 0.308379    
healthyfatCat:unhealthyfatCat           4   3.65   0.913   0.567 0.687045    
sfblCat:healthyfatCat:unhealthyfatCat   8  11.10   1.387   0.861 0.550610    
Residuals                             151 243.13   1.610                     
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

eta_squared(fatHEI.aov)

                              sfblCat                         healthyfatCat 
                          0.029247148                           0.019130726 
                      unhealthyfatCat                                   age 
                          0.086035818                           0.008568904 
                       plantmilk_freq                 sfblCat:healthyfatCat 
                          0.005624827                           0.006787946 
              sfblCat:unhealthyfatCat         healthyfatCat:unhealthyfatCat 
                          0.024801561                           0.011607088 
sfblCat:healthyfatCat:unhealthyfatCat 
                          0.035272870 

emmeans(fatHEI.aov, pairwise ~ sfblCat, adjust = "tukey")

NOTE: Results may be misleading due to involvement in interactions

$emmeans
 sfblCat emmean    SE  df lower.CL upper.CL
 Low       3.45 0.190 151     3.07     3.82
 Medium    3.62 0.172 151     3.28     3.96
 High      3.09 0.186 151     2.73     3.46

Results are averaged over the levels of: healthyfatCat, unhealthyfatCat 
Confidence level used: 0.95 

$contrasts
 contrast      estimate    SE  df t.ratio p.value
 Low - Medium    -0.174 0.258 151  -0.677  0.7775
 Low - High       0.354 0.266 151   1.335  0.3783
 Medium - High    0.529 0.252 151   2.101  0.0931

Results are averaged over the levels of: healthyfatCat, unhealthyfatCat 
P value adjustment: tukey method for comparing a family of 3 estimates 

emmeans(fatHEI.aov, ~sfblCat)

NOTE: Results may be misleading due to involvement in interactions

 sfblCat emmean    SE  df lower.CL upper.CL
 Low       3.45 0.190 151     3.07     3.82
 Medium    3.62 0.172 151     3.28     3.96
 High      3.09 0.186 151     2.73     3.46

Results are averaged over the levels of: healthyfatCat, unhealthyfatCat 
Confidence level used: 0.95 

emmeans(fatHEI.aov, pairwise ~ healthyfatCat, adjust = "tukey")

NOTE: Results may be misleading due to involvement in interactions

$emmeans
 healthyfatCat emmean    SE  df lower.CL upper.CL
 Low             3.26 0.191 151     2.89     3.64
 Medium          3.66 0.182 151     3.30     4.02
 High            3.24 0.180 151     2.88     3.60

Results are averaged over the levels of: sfblCat, unhealthyfatCat 
Confidence level used: 0.95 

$contrasts
 contrast      estimate    SE  df t.ratio p.value
 Low - Medium   -0.3952 0.263 151  -1.502  0.2928
 Low - High      0.0231 0.270 151   0.086  0.9960
 Medium - High   0.4183 0.257 151   1.625  0.2384

Results are averaged over the levels of: sfblCat, unhealthyfatCat 
P value adjustment: tukey method for comparing a family of 3 estimates 

emmeans(fatHEI.aov, ~healthyfatCat)

NOTE: Results may be misleading due to involvement in interactions

 healthyfatCat emmean    SE  df lower.CL upper.CL
 Low             3.26 0.191 151     2.89     3.64
 Medium          3.66 0.182 151     3.30     4.02
 High            3.24 0.180 151     2.88     3.60

Results are averaged over the levels of: sfblCat, unhealthyfatCat 
Confidence level used: 0.95 

emmeans(fatHEI.aov, pairwise ~ unhealthyfatCat, adjust = "tukey")

NOTE: Results may be misleading due to involvement in interactions

$emmeans
 unhealthyfatCat emmean    SE  df lower.CL upper.CL
 Low               3.88 0.184 151     3.52     4.25
 Medium            3.21 0.169 151     2.87     3.54
 High              3.07 0.195 151     2.68     3.46

Results are averaged over the levels of: sfblCat, healthyfatCat 
Confidence level used: 0.95 

$contrasts
 contrast      estimate    SE  df t.ratio p.value
 Low - Medium     0.677 0.250 151   2.708  0.0205
 Low - High       0.813 0.272 151   2.989  0.0091
 Medium - High    0.136 0.258 151   0.528  0.8579

Results are averaged over the levels of: sfblCat, healthyfatCat 
P value adjustment: tukey method for comparing a family of 3 estimates 

emmeans(fatHEI.aov, ~unhealthyfatCat)

NOTE: Results may be misleading due to involvement in interactions

 unhealthyfatCat emmean    SE  df lower.CL upper.CL
 Low               3.88 0.184 151     3.52     4.25
 Medium            3.21 0.169 151     2.87     3.54
 High              3.07 0.195 151     2.68     3.46

Results are averaged over the levels of: sfblCat, healthyfatCat 
Confidence level used: 0.95 

## Graph

# Significance tests
sfbl.fatHEI.sig <- full.df %>%
  tukey_hsd(fatHEI ~ sfblCat) %>%
  add_significance()

healthyfat.fatHEI.sig <- full.df %>%
  tukey_hsd(fatHEI ~ healthyfatCat) %>%
  add_significance()

unhealthyfat.fatHEI.sig <- full.df %>%
  tukey_hsd(fatHEI ~ unhealthyfatCat) %>%
  add_significance()

# Get emmeans
sfbl.fatHEI.emm <- emmeans(fatHEI.aov, ~ sfblCat) %>% as.data.frame()

NOTE: Results may be misleading due to involvement in interactions

healthyfat.fatHEI.emm <- emmeans(fatHEI.aov, ~ healthyfatCat) %>% as.data.frame()

NOTE: Results may be misleading due to involvement in interactions

unhealthyfat.fatHEI.emm <- emmeans(fatHEI.aov, ~ unhealthyfatCat) %>% as.data.frame()

NOTE: Results may be misleading due to involvement in interactions

# sfblCat plot
sfbl.fatHEI.plot <- ggplot(sfbl.fatHEI.emm, aes(x = sfblCat, y = emmean, fill = sfblCat)) +
  geom_bar(stat = "identity", width = 0.5, alpha = 0.8) +
  geom_errorbar(aes(ymin = lower.CL, ymax = upper.CL), width = 0.2) +
  scale_fill_manual(values = c("Low" = "#00897b", "Medium" = "#ffb300", "High" = "#e53935")) +
  labs(x = "Sweet Fat Beverage Liking", y = "Fat HEI", title = " ") +
  ylim(0, 6) +
  theme_classic() +
  theme(legend.position = "none",
        axis.text = element_text(size = 14),
        axis.title = element_text(size = 16))

# healthyfatCat plot
healthyfat.fatHEI.plot <- ggplot(healthyfat.fatHEI.emm, aes(x = healthyfatCat, y = emmean, fill = healthyfatCat)) +
  geom_bar(stat = "identity", width = 0.5, alpha = 0.8) +
  geom_errorbar(aes(ymin = lower.CL, ymax = upper.CL), width = 0.2) +
  scale_fill_manual(values = c("Low" = "#00897b", "Medium" = "#ffb300", "High" = "#e53935")) +
  labs(x = "Healthy Fat Liking", y = " ", title = " ") +
  ylim(0, 6) +
  theme_classic() +
  theme(legend.position = "none",
        axis.text = element_text(size = 14),
        axis.title = element_text(size = 16))

# unhealthyfatCat plot
unhealthyfat.fatHEI.plot <- ggplot(unhealthyfat.fatHEI.emm, aes(x = unhealthyfatCat, y = emmean, fill = unhealthyfatCat)) +
  geom_bar(stat = "identity", width = 0.5, alpha = 0.8) +
  geom_errorbar(aes(ymin = lower.CL, ymax = upper.CL), width = 0.2) +
  stat_pvalue_manual(unhealthyfat.fatHEI.sig, hide.ns = TRUE, label = "p.adj.signif",
                     label.size = 10, bracket.size = 1, tip.length = 0.05,
                     y.position = c(4.7, 5.5)) +
  scale_fill_manual(values = c("Low" = "#00897b", "Medium" = "#ffb300", "High" = "#e53935")) +
  labs(x = "Unhealthy Fat Liking", y = " ", title = " ") +
  ylim(0, 6) +
  theme_classic() +
  theme(legend.position = "none",
        axis.text = element_text(size = 14),
        axis.title = element_text(size = 16))

# Stitch together
sfbl.fatHEI.plot + healthyfat.fatHEI.plot + unhealthyfat.fatHEI.plot +
  plot_layout(axis_titles = "collect")

ggsave("liking.fatHEI.png", height = 5, width = 9)

## Sugar HEI
sugarHEI.aov <- full.df %>%
  aov(sugarHEI ~ sfblCat*healthyfatCat*unhealthyfatCat + age + plantmilk_freq, .)
summary(sugarHEI.aov)

                                       Df Sum Sq Mean Sq F value   Pr(>F)    
sfblCat                                 2  243.6  121.81  22.654 2.49e-09 ***
healthyfatCat                           2    3.2    1.59   0.296  0.74402    
unhealthyfatCat                         2   46.6   23.29   4.332  0.01481 *  
age                                     1    1.6    1.60   0.297  0.58660    
plantmilk_freq                          1   37.3   37.33   6.943  0.00929 ** 
sfblCat:healthyfatCat                   4   41.9   10.47   1.947  0.10555    
sfblCat:unhealthyfatCat                 4   40.7   10.18   1.894  0.11441    
healthyfatCat:unhealthyfatCat           4    4.5    1.11   0.207  0.93402    
sfblCat:healthyfatCat:unhealthyfatCat   8   73.6    9.20   1.711  0.10010    
Residuals                             151  811.9    5.38                     
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

summary.aov(sugarHEI.aov)

                                       Df Sum Sq Mean Sq F value   Pr(>F)    
sfblCat                                 2  243.6  121.81  22.654 2.49e-09 ***
healthyfatCat                           2    3.2    1.59   0.296  0.74402    
unhealthyfatCat                         2   46.6   23.29   4.332  0.01481 *  
age                                     1    1.6    1.60   0.297  0.58660    
plantmilk_freq                          1   37.3   37.33   6.943  0.00929 ** 
sfblCat:healthyfatCat                   4   41.9   10.47   1.947  0.10555    
sfblCat:unhealthyfatCat                 4   40.7   10.18   1.894  0.11441    
healthyfatCat:unhealthyfatCat           4    4.5    1.11   0.207  0.93402    
sfblCat:healthyfatCat:unhealthyfatCat   8   73.6    9.20   1.711  0.10010    
Residuals                             151  811.9    5.38                     
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

eta_squared(sugarHEI.aov)

                              sfblCat                         healthyfatCat 
                          0.186695051                           0.002441545 
                      unhealthyfatCat                                   age 
                          0.035704680                           0.001223653 
                       plantmilk_freq                 sfblCat:healthyfatCat 
                          0.028609153                           0.032089391 
              sfblCat:unhealthyfatCat         healthyfatCat:unhealthyfatCat 
                          0.031212442                           0.003417822 
sfblCat:healthyfatCat:unhealthyfatCat 
                          0.056397138 

emmeans(sugarHEI.aov, pairwise ~ sfblCat, adjust = "tukey")

NOTE: Results may be misleading due to involvement in interactions

$emmeans
 sfblCat emmean    SE  df lower.CL upper.CL
 Low       5.80 0.346 151     5.11     6.48
 Medium    3.33 0.313 151     2.71     3.95
 High      3.36 0.339 151     2.69     4.03

Results are averaged over the levels of: healthyfatCat, unhealthyfatCat 
Confidence level used: 0.95 

$contrasts
 contrast      estimate    SE  df t.ratio p.value
 Low - Medium    2.4711 0.471 151   5.248  <.0001
 Low - High      2.4418 0.485 151   5.032  <.0001
 Medium - High  -0.0293 0.460 151  -0.064  0.9978

Results are averaged over the levels of: healthyfatCat, unhealthyfatCat 
P value adjustment: tukey method for comparing a family of 3 estimates 

emmeans(sugarHEI.aov, ~sfblCat)

NOTE: Results may be misleading due to involvement in interactions

 sfblCat emmean    SE  df lower.CL upper.CL
 Low       5.80 0.346 151     5.11     6.48
 Medium    3.33 0.313 151     2.71     3.95
 High      3.36 0.339 151     2.69     4.03

Results are averaged over the levels of: healthyfatCat, unhealthyfatCat 
Confidence level used: 0.95 

emmeans(sugarHEI.aov, pairwise ~ healthyfatCat, adjust = "tukey")

NOTE: Results may be misleading due to involvement in interactions

$emmeans
 healthyfatCat emmean    SE  df lower.CL upper.CL
 Low             3.71 0.350 151     3.02     4.40
 Medium          4.10 0.333 151     3.44     4.76
 High            4.67 0.330 151     4.02     5.32

Results are averaged over the levels of: sfblCat, unhealthyfatCat 
Confidence level used: 0.95 

$contrasts
 contrast      estimate    SE  df t.ratio p.value
 Low - Medium    -0.391 0.481 151  -0.812  0.6959
 Low - High      -0.958 0.493 151  -1.942  0.1305
 Medium - High   -0.567 0.470 151  -1.205  0.4519

Results are averaged over the levels of: sfblCat, unhealthyfatCat 
P value adjustment: tukey method for comparing a family of 3 estimates 

emmeans(sugarHEI.aov, ~healthyfatCat)

NOTE: Results may be misleading due to involvement in interactions

 healthyfatCat emmean    SE  df lower.CL upper.CL
 Low             3.71 0.350 151     3.02     4.40
 Medium          4.10 0.333 151     3.44     4.76
 High            4.67 0.330 151     4.02     5.32

Results are averaged over the levels of: sfblCat, unhealthyfatCat 
Confidence level used: 0.95 

emmeans(sugarHEI.aov, pairwise ~ unhealthyfatCat, adjust = "tukey")

NOTE: Results may be misleading due to involvement in interactions

$emmeans
 unhealthyfatCat emmean    SE  df lower.CL upper.CL
 Low               5.06 0.337 151     4.39     5.72
 Medium            3.97 0.309 151     3.36     4.58
 High              3.46 0.357 151     2.76     4.17

Results are averaged over the levels of: sfblCat, healthyfatCat 
Confidence level used: 0.95 

$contrasts
 contrast      estimate    SE  df t.ratio p.value
 Low - Medium     1.089 0.457 151   2.384  0.0480
 Low - High       1.593 0.497 151   3.202  0.0047
 Medium - High    0.503 0.472 151   1.067  0.5360

Results are averaged over the levels of: sfblCat, healthyfatCat 
P value adjustment: tukey method for comparing a family of 3 estimates 

emmeans(sugarHEI.aov, ~unhealthyfatCat)

NOTE: Results may be misleading due to involvement in interactions

 unhealthyfatCat emmean    SE  df lower.CL upper.CL
 Low               5.06 0.337 151     4.39     5.72
 Medium            3.97 0.309 151     3.36     4.58
 High              3.46 0.357 151     2.76     4.17

Results are averaged over the levels of: sfblCat, healthyfatCat 
Confidence level used: 0.95 

## Graph
# Significance tests
sfbl.sugarHEI.sig <- full.df %>%
  tukey_hsd(sugarHEI ~ sfblCat) %>%
  add_significance()

healthyfat.sugarHEI.sig <- full.df %>%
  tukey_hsd(sugarHEI ~ healthyfatCat) %>%
  add_significance()

unhealthyfat.sugarHEI.sig <- full.df %>%
  tukey_hsd(sugarHEI ~ unhealthyfatCat) %>%
  add_significance()

# Get emmeans
sfbl.sugarHEI.emm <- emmeans(sugarHEI.aov, ~ sfblCat) %>% as.data.frame()

NOTE: Results may be misleading due to involvement in interactions

healthyfat.sugarHEI.emm <- emmeans(sugarHEI.aov, ~ healthyfatCat) %>% as.data.frame()

NOTE: Results may be misleading due to involvement in interactions

unhealthyfat.sugarHEI.emm <- emmeans(sugarHEI.aov, ~ unhealthyfatCat) %>% as.data.frame()

NOTE: Results may be misleading due to involvement in interactions

# sfblCat plot
sfbl.sugarHEI.plot <- ggplot(sfbl.sugarHEI.emm, aes(x = sfblCat, y = emmean, fill = sfblCat)) +
  geom_bar(stat = "identity", width = 0.5, alpha = 0.8) +
  geom_errorbar(aes(ymin = lower.CL, ymax = upper.CL), width = 0.2) +
  stat_pvalue_manual(sfbl.sugarHEI.sig, hide.ns = TRUE, label = "p.adj.signif",
                     label.size = 10, bracket.size = 1, tip.length = 0.05,
                     y.position = c(7, 8)) +
  scale_fill_manual(values = c("Low" = "#00897b", "Medium" = "#ffb300", "High" = "#e53935")) +
  labs(x = "Sweet Fat Beverage Liking", y = "Sugar HEI", title = " ") +
  ylim(0, 9) +
  theme_classic() +
  theme(legend.position = "none",
        axis.text = element_text(size = 14),
        axis.title = element_text(size = 16))

# healthyfatCat plot
healthyfat.sugarHEI.plot <- ggplot(healthyfat.sugarHEI.emm, aes(x = healthyfatCat, y = emmean, fill = healthyfatCat)) +
  geom_bar(stat = "identity", width = 0.5, alpha = 0.8) +
  geom_errorbar(aes(ymin = lower.CL, ymax = upper.CL), width = 0.2) +
  scale_fill_manual(values = c("Low" = "#00897b", "Medium" = "#ffb300", "High" = "#e53935")) +
  labs(x = "Healthy Fat Liking", y = " ", title = " ") +
  ylim(0, 9) +
  theme_classic() +
  theme(legend.position = "none",
        axis.text = element_text(size = 14),
        axis.title = element_text(size = 16))

# unhealthyfatCat plot
unhealthyfat.sugarHEI.plot <- ggplot(unhealthyfat.sugarHEI.emm, aes(x = unhealthyfatCat, y = emmean, fill = unhealthyfatCat)) +
  geom_bar(stat = "identity", width = 0.5, alpha = 0.8) +
  geom_errorbar(aes(ymin = lower.CL, ymax = upper.CL), width = 0.2) +
  stat_pvalue_manual(unhealthyfat.sugarHEI.sig, hide.ns = TRUE, label = "p.adj.signif",
                     label.size = 10, bracket.size = 1, tip.length = 0.05,
                     y.position = c(6.5)) +
  scale_fill_manual(values = c("Low" = "#00897b", "Medium" = "#ffb300", "High" = "#e53935")) +
  labs(x = "Unhealthy Fat Liking", y = " ", title = " ") +
  ylim(0, 9) +
  theme_classic() +
  theme(legend.position = "none",
        axis.text = element_text(size = 14),
        axis.title = element_text(size = 16))

# Stitch together
sfbl.sugarHEI.plot + healthyfat.sugarHEI.plot + unhealthyfat.sugarHEI.plot +
  plot_layout(axis_titles = "collect")

ggsave("liking.sugarHEI.png", height = 5, width = 9)


## Sugar intake
sugarintake.aov <- full.df %>%
  aov(sugar.intake ~ sfblCat*healthyfatCat*unhealthyfatCat + age + plantmilk_freq, .)
summary(sugarintake.aov)

                                       Df  Sum Sq Mean Sq F value   Pr(>F)    
sfblCat                                 2 1465981  732991  12.621 8.54e-06 ***
healthyfatCat                           2   69661   34830   0.600   0.5503    
unhealthyfatCat                         2  416338  208169   3.584   0.0301 *  
age                                     1  122798  122798   2.114   0.1480    
plantmilk_freq                          1  209658  209658   3.610   0.0593 .  
sfblCat:healthyfatCat                   4  288561   72140   1.242   0.2955    
sfblCat:unhealthyfatCat                 4  321415   80354   1.384   0.2423    
healthyfatCat:unhealthyfatCat           4  201351   50338   0.867   0.4854    
sfblCat:healthyfatCat:unhealthyfatCat   8  455937   56992   0.981   0.4528    
Residuals                             151 8769297   58075                     
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

summary.aov(sugarintake.aov)

                                       Df  Sum Sq Mean Sq F value   Pr(>F)    
sfblCat                                 2 1465981  732991  12.621 8.54e-06 ***
healthyfatCat                           2   69661   34830   0.600   0.5503    
unhealthyfatCat                         2  416338  208169   3.584   0.0301 *  
age                                     1  122798  122798   2.114   0.1480    
plantmilk_freq                          1  209658  209658   3.610   0.0593 .  
sfblCat:healthyfatCat                   4  288561   72140   1.242   0.2955    
sfblCat:unhealthyfatCat                 4  321415   80354   1.384   0.2423    
healthyfatCat:unhealthyfatCat           4  201351   50338   0.867   0.4854    
sfblCat:healthyfatCat:unhealthyfatCat   8  455937   56992   0.981   0.4528    
Residuals                             151 8769297   58075                     
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

eta_squared(sugarintake.aov)

                              sfblCat                         healthyfatCat 
                          0.118982339                           0.005653828 
                      unhealthyfatCat                                   age 
                          0.033790943                           0.009966599 
                       plantmilk_freq                 sfblCat:healthyfatCat 
                          0.017016319                           0.023420232 
              sfblCat:unhealthyfatCat         healthyfatCat:unhealthyfatCat 
                          0.026086753                           0.016342134 
sfblCat:healthyfatCat:unhealthyfatCat 
                          0.037004894 

emmeans(sugarintake.aov, pairwise ~ sfblCat, adjust = "tukey")

NOTE: Results may be misleading due to involvement in interactions

$emmeans
 sfblCat emmean   SE  df lower.CL upper.CL
 Low        341 36.0 151      270      412
 Medium     493 32.6 151      429      558
 High       586 35.3 151      516      655

Results are averaged over the levels of: healthyfatCat, unhealthyfatCat 
Confidence level used: 0.95 

$contrasts
 contrast      estimate   SE  df t.ratio p.value
 Low - Medium    -152.7 48.9 151  -3.121  0.0061
 Low - High      -245.1 50.4 151  -4.861  <.0001
 Medium - High    -92.4 47.8 151  -1.933  0.1331

Results are averaged over the levels of: healthyfatCat, unhealthyfatCat 
P value adjustment: tukey method for comparing a family of 3 estimates 

emmeans(sugarintake.aov, ~sfblCat)

NOTE: Results may be misleading due to involvement in interactions

 sfblCat emmean   SE  df lower.CL upper.CL
 Low        341 36.0 151      270      412
 Medium     493 32.6 151      429      558
 High       586 35.3 151      516      655

Results are averaged over the levels of: healthyfatCat, unhealthyfatCat 
Confidence level used: 0.95 

emmeans(sugarintake.aov, pairwise ~ healthyfatCat, adjust = "tukey")

NOTE: Results may be misleading due to involvement in interactions

$emmeans
 healthyfatCat emmean   SE  df lower.CL upper.CL
 Low              523 36.4 151      451      594
 Medium           475 34.6 151      407      544
 High             422 34.3 151      354      489

Results are averaged over the levels of: sfblCat, unhealthyfatCat 
Confidence level used: 0.95 

$contrasts
 contrast      estimate   SE  df t.ratio p.value
 Low - Medium      47.3 50.0 151   0.947  0.6114
 Low - High       100.9 51.2 151   1.968  0.1237
 Medium - High     53.5 48.9 151   1.095  0.5187

Results are averaged over the levels of: sfblCat, unhealthyfatCat 
P value adjustment: tukey method for comparing a family of 3 estimates 

emmeans(sugarintake.aov, ~healthyfatCat)

NOTE: Results may be misleading due to involvement in interactions

 healthyfatCat emmean   SE  df lower.CL upper.CL
 Low              523 36.4 151      451      594
 Medium           475 34.6 151      407      544
 High             422 34.3 151      354      489

Results are averaged over the levels of: sfblCat, unhealthyfatCat 
Confidence level used: 0.95 

emmeans(sugarintake.aov, pairwise ~ unhealthyfatCat, adjust = "tukey")

NOTE: Results may be misleading due to involvement in interactions

$emmeans
 unhealthyfatCat emmean   SE  df lower.CL upper.CL
 Low                396 35.0 151      327      465
 Medium             489 32.1 151      426      553
 High               534 37.1 151      461      607

Results are averaged over the levels of: sfblCat, healthyfatCat 
Confidence level used: 0.95 

$contrasts
 contrast      estimate   SE  df t.ratio p.value
 Low - Medium     -93.3 47.5 151  -1.964  0.1248
 Low - High      -137.8 51.7 151  -2.667  0.0230
 Medium - High    -44.6 49.0 151  -0.909  0.6353

Results are averaged over the levels of: sfblCat, healthyfatCat 
P value adjustment: tukey method for comparing a family of 3 estimates 

emmeans(sugarintake.aov, ~healthyfatCat)

NOTE: Results may be misleading due to involvement in interactions

 healthyfatCat emmean   SE  df lower.CL upper.CL
 Low              523 36.4 151      451      594
 Medium           475 34.6 151      407      544
 High             422 34.3 151      354      489

Results are averaged over the levels of: sfblCat, unhealthyfatCat 
Confidence level used: 0.95 

Group difference, 2-way ANCOVA, anthropometrics

## Lean mass
leanmass.aov <- full.df %>%
  aov(lean_mass ~ as.factor(sweetPref) * as.factor(fatPref) + age + plantmilk_freq, .)
summary(leanmass.aov)

                                         Df Sum Sq Mean Sq F value Pr(>F)  
as.factor(sweetPref)                      1   12.4   12.42   0.565 0.4533  
as.factor(fatPref)                        2   83.9   41.95   1.910 0.1519  
age                                       1   75.3   75.33   3.430 0.0661 .
plantmilk_freq                            1    1.5    1.46   0.067 0.7966  
as.factor(sweetPref):as.factor(fatPref)   2   86.8   43.39   1.975 0.1425  
Residuals                               142 3118.7   21.96                 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
30 observations deleted due to missingness

eta_squared(leanmass.aov)

                   as.factor(sweetPref)                      as.factor(fatPref) 
                           0.0036759507                            0.0248333697 
                                    age                          plantmilk_freq 
                           0.0222956756                            0.0004333148 
as.factor(sweetPref):as.factor(fatPref) 
                           0.0256833707 

## SMM
SMM.aov <- full.df %>%
  aov(SMM ~ as.factor(sweetPref) * as.factor(fatPref) + age + plantmilk_freq, .)
summary(SMM.aov)

                                         Df Sum Sq Mean Sq F value Pr(>F)  
as.factor(sweetPref)                      1     59    59.0   0.533 0.4667  
as.factor(fatPref)                        2    488   244.0   2.204 0.1141  
age                                       1    384   383.7   3.466 0.0647 .
plantmilk_freq                            1      8     8.1   0.073 0.7869  
as.factor(sweetPref):as.factor(fatPref)   2    461   230.5   2.082 0.1284  
Residuals                               142  15720   110.7                 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
30 observations deleted due to missingness

eta_squared(SMM.aov)

                   as.factor(sweetPref)                      as.factor(fatPref) 
                           0.0034444251                            0.0285081421 
                                    age                          plantmilk_freq 
                           0.0224141411                            0.0004741942 
as.factor(sweetPref):as.factor(fatPref) 
                           0.0269289267 

## SMM, weight adjusted
SMMwt.aov <- full.df %>%
    aov(SMM_wt ~ sweetPref * fatPref + age + plantmilk_freq, .)
summary(SMMwt.aov)

                   Df Sum Sq Mean Sq F value Pr(>F)  
sweetPref           1   55.5   55.54   3.465 0.0647 .
fatPref             2   63.0   31.51   1.966 0.1438  
age                 1   42.3   42.30   2.639 0.1065  
plantmilk_freq      1   17.3   17.29   1.079 0.3008  
sweetPref:fatPref   2    4.4    2.19   0.137 0.8725  
Residuals         142 2275.9   16.03                 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
30 observations deleted due to missingness

emmeans(SMMwt.aov, pairwise ~ sweetPref, adjust = "tukey")

NOTE: Results may be misleading due to involvement in interactions

$emmeans
 sweetPref emmean    SE  df lower.CL upper.CL
 Low         36.8 0.788 142     35.3     38.4
 High        35.5 0.379 142     34.7     36.2

Results are averaged over the levels of: fatPref 
Confidence level used: 0.95 

$contrasts
 contrast   estimate    SE  df t.ratio p.value
 Low - High     1.33 0.881 142   1.516  0.1318

Results are averaged over the levels of: fatPref 

emmeans(SMMwt.aov, ~sweetPref)

NOTE: Results may be misleading due to involvement in interactions

 sweetPref emmean    SE  df lower.CL upper.CL
 Low         36.8 0.788 142     35.3     38.4
 High        35.5 0.379 142     34.7     36.2

Results are averaged over the levels of: fatPref 
Confidence level used: 0.95 

eta_squared(SMMwt.aov)

        sweetPref           fatPref               age    plantmilk_freq 
      0.022592417       0.025634297       0.017206650       0.007031179 
sweetPref:fatPref 
      0.001780440 

## FFM
FFM.aov <- full.df %>%
    aov(FFM ~ as.factor(sweetPref)*as.factor(fatPref) + age + plantmilk_freq, .)
summary(FFM.aov)

                                         Df Sum Sq Mean Sq F value Pr(>F)  
as.factor(sweetPref)                      1    189   189.2   0.619 0.4327  
as.factor(fatPref)                        2   1291   645.7   2.112 0.1247  
age                                       1   1045  1044.6   3.417 0.0666 .
plantmilk_freq                            1     13    13.3   0.043 0.8352  
as.factor(sweetPref):as.factor(fatPref)   2   1226   613.0   2.005 0.1384  
Residuals                               142  43409   305.7                 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
30 observations deleted due to missingness

eta_squared(FFM.aov)

                   as.factor(sweetPref)                      as.factor(fatPref) 
                           0.0040112449                            0.0273773704 
                                    age                          plantmilk_freq 
                           0.0221439057                            0.0002814142 
as.factor(sweetPref):as.factor(fatPref) 
                           0.0259892585 

## Trunk lean mass
trunklean.aov <- full.df %>%
    aov(trunk_lean_mass ~ as.factor(sweetPref)*as.factor(fatPref) + age + plantmilk_freq, .)
summary(trunklean.aov)

                                         Df Sum Sq Mean Sq F value Pr(>F)  
as.factor(sweetPref)                      1     56   55.69   0.874 0.3515  
as.factor(fatPref)                        2    324  161.89   2.540 0.0825 .
age                                       1    287  287.02   4.503 0.0356 *
plantmilk_freq                            1      0    0.19   0.003 0.9562  
as.factor(sweetPref):as.factor(fatPref)   2    318  158.85   2.492 0.0864 .
Residuals                               142   9052   63.74                 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
30 observations deleted due to missingness

eta_squared(trunklean.aov)

                   as.factor(sweetPref)                      as.factor(fatPref) 
                           5.548988e-03                            3.226186e-02 
                                    age                          plantmilk_freq 
                           2.859941e-02                            1.922404e-05 
as.factor(sweetPref):as.factor(fatPref) 
                           3.165650e-02 

## BFM
BFMfat.aov <- full.df %>%
    aov(BFM ~ as.factor(sweetPref)*as.factor(fatPref) + age + plantmilk_freq, .)
summary(BFMfat.aov)

                                         Df Sum Sq Mean Sq F value Pr(>F)  
as.factor(sweetPref)                      1   1695  1694.7   2.719 0.1014  
as.factor(fatPref)                        2   5052  2526.0   4.053 0.0194 *
age                                       1   2934  2934.2   4.708 0.0317 *
plantmilk_freq                            1    511   510.9   0.820 0.3668  
as.factor(sweetPref):as.factor(fatPref)   2   1027   513.7   0.824 0.4407  
Residuals                               142  88508   623.3                 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
30 observations deleted due to missingness

emmeans(BFMfat.aov, pairwise ~ fatPref, adjust = "tukey")

NOTE: Results may be misleading due to involvement in interactions

$emmeans
 fatPref emmean   SE  df lower.CL upper.CL
 Low       55.7 5.63 142     44.6     66.9
 Medium    46.7 3.65 142     39.5     53.9
 High      43.6 4.75 142     34.2     53.0

Results are averaged over the levels of: sweetPref 
Confidence level used: 0.95 

$contrasts
 contrast      estimate   SE  df t.ratio p.value
 Low - Medium      9.04 6.70 142   1.349  0.3708
 Low - High       12.16 7.45 142   1.632  0.2356
 Medium - High     3.12 6.04 142   0.517  0.8633

Results are averaged over the levels of: sweetPref 
P value adjustment: tukey method for comparing a family of 3 estimates 

emmeans(BFMfat.aov, ~ fatPref)

NOTE: Results may be misleading due to involvement in interactions

 fatPref emmean   SE  df lower.CL upper.CL
 Low       55.7 5.63 142     44.6     66.9
 Medium    46.7 3.65 142     39.5     53.9
 High      43.6 4.75 142     34.2     53.0

Results are averaged over the levels of: sweetPref 
Confidence level used: 0.95 

eta_squared(BFMfat.aov)

                   as.factor(sweetPref)                      as.factor(fatPref) 
                            0.016993469                             0.050658979 
                                    age                          plantmilk_freq 
                            0.029422240                             0.005123125 
as.factor(sweetPref):as.factor(fatPref) 
                            0.010302228 

## PBF
PBFfat.aov <- full.df %>%
      aov(PBF ~ as.factor(sweetPref)*as.factor(fatPref) + age + plantmilk_freq, .)
summary(PBFfat.aov)

                                         Df Sum Sq Mean Sq F value Pr(>F)  
as.factor(sweetPref)                      1    192  191.97   3.418 0.0666 .
as.factor(fatPref)                        2    275  137.64   2.451 0.0899 .
age                                       1    197  197.04   3.509 0.0631 .
plantmilk_freq                            1     37   37.10   0.661 0.4177  
as.factor(sweetPref):as.factor(fatPref)   2     36   17.87   0.318 0.7280  
Residuals                               142   7975   56.16                 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
30 observations deleted due to missingness

emmeans(PBFfat.aov, pairwise ~ fatPref, adjust = "tukey")

NOTE: Results may be misleading due to involvement in interactions

$emmeans
 fatPref emmean   SE  df lower.CL upper.CL
 Low       35.3 1.69 142     32.0     38.6
 Medium    32.9 1.10 142     30.8     35.1
 High      31.6 1.43 142     28.8     34.4

Results are averaged over the levels of: sweetPref 
Confidence level used: 0.95 

$contrasts
 contrast      estimate   SE  df t.ratio p.value
 Low - Medium      2.37 2.01 142   1.180  0.4668
 Low - High        3.71 2.24 142   1.659  0.2246
 Medium - High     1.34 1.81 142   0.737  0.7419

Results are averaged over the levels of: sweetPref 
P value adjustment: tukey method for comparing a family of 3 estimates 

emmeans(PBFfat.aov, ~ fatPref)

NOTE: Results may be misleading due to involvement in interactions

 fatPref emmean   SE  df lower.CL upper.CL
 Low       35.3 1.69 142     32.0     38.6
 Medium    32.9 1.10 142     30.8     35.1
 High      31.6 1.43 142     28.8     34.4

Results are averaged over the levels of: sweetPref 
Confidence level used: 0.95 

eta_squared(PBFfat.aov)

                   as.factor(sweetPref)                      as.factor(fatPref) 
                            0.022035243                             0.031599473 
                                    age                          plantmilk_freq 
                            0.022618116                             0.004258686 
as.factor(sweetPref):as.factor(fatPref) 
                            0.004101613 

## Trunk fat mass
trunkfat.aov <- full.df %>%
    aov(trunk_fat_mass ~ as.factor(sweetPref)*as.factor(fatPref) + age + plantmilk_freq, .)
summary(trunkfat.aov)

                                         Df Sum Sq Mean Sq F value Pr(>F)  
as.factor(sweetPref)                      1    374   373.9   2.649 0.1058  
as.factor(fatPref)                        2    965   482.5   3.419 0.0355 *
age                                       1    868   868.1   6.151 0.0143 *
plantmilk_freq                            1     85    84.9   0.601 0.4393  
as.factor(sweetPref):as.factor(fatPref)   2    287   143.5   1.017 0.3644  
Residuals                               142  20041   141.1                 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
30 observations deleted due to missingness

emmeans(trunkfat.aov, pairwise ~ fatPref, adjust = "tukey")

NOTE: Results may be misleading due to involvement in interactions

$emmeans
 fatPref emmean   SE  df lower.CL upper.CL
 Low       26.6 2.68 142     21.3     31.9
 Medium    22.7 1.74 142     19.3     26.1
 High      21.0 2.26 142     16.5     25.4

Results are averaged over the levels of: sweetPref 
Confidence level used: 0.95 

$contrasts
 contrast      estimate   SE  df t.ratio p.value
 Low - Medium      3.87 3.19 142   1.215  0.4465
 Low - High        5.60 3.55 142   1.581  0.2573
 Medium - High     1.73 2.87 142   0.602  0.8192

Results are averaged over the levels of: sweetPref 
P value adjustment: tukey method for comparing a family of 3 estimates 

emmeans(trunkfat.aov, ~ fatPref)

NOTE: Results may be misleading due to involvement in interactions

 fatPref emmean   SE  df lower.CL upper.CL
 Low       26.6 2.68 142     21.3     31.9
 Medium    22.7 1.74 142     19.3     26.1
 High      21.0 2.26 142     16.5     25.4

Results are averaged over the levels of: sweetPref 
Confidence level used: 0.95 

eta_squared(trunkfat.aov)

                   as.factor(sweetPref)                      as.factor(fatPref) 
                            0.016529873                             0.042663900 
                                    age                          plantmilk_freq 
                            0.038380162                             0.003752572 
as.factor(sweetPref):as.factor(fatPref) 
                            0.012686597 

## Visceral fat mass
visceralfat.aov <- full.df %>%
    aov(visceral_fat ~ as.factor(sweetPref)*as.factor(fatPref) + age + plantmilk_freq, .)
summary(visceralfat.aov)

                                         Df Sum Sq Mean Sq F value Pr(>F)  
as.factor(sweetPref)                      1   8347    8347   2.937 0.0887 .
as.factor(fatPref)                        2  17440    8720   3.069 0.0496 *
age                                       1  17505   17505   6.160 0.0142 *
plantmilk_freq                            1   3233    3233   1.138 0.2879  
as.factor(sweetPref):as.factor(fatPref)   2   4475    2237   0.787 0.4570  
Residuals                               142 403510    2842                 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
30 observations deleted due to missingness

emmeans(visceralfat.aov, pairwise ~ fatPref, adjust = "tukey")

NOTE: Results may be misleading due to involvement in interactions

$emmeans
 fatPref emmean    SE  df lower.CL upper.CL
 Low      112.9 12.00 142     89.1      137
 Medium    96.2  7.79 142     80.8      112
 High      89.3 10.10 142     69.3      109

Results are averaged over the levels of: sweetPref 
Confidence level used: 0.95 

$contrasts
 contrast      estimate   SE  df t.ratio p.value
 Low - Medium      16.7 14.3 142   1.165  0.4760
 Low - High        23.6 15.9 142   1.481  0.3029
 Medium - High      6.9 12.9 142   0.535  0.8543

Results are averaged over the levels of: sweetPref 
P value adjustment: tukey method for comparing a family of 3 estimates 

emmeans(visceralfat.aov, ~ fatPref)

NOTE: Results may be misleading due to involvement in interactions

 fatPref emmean    SE  df lower.CL upper.CL
 Low      112.9 12.00 142     89.1      137
 Medium    96.2  7.79 142     80.8      112
 High      89.3 10.10 142     69.3      109

Results are averaged over the levels of: sweetPref 
Confidence level used: 0.95 

eta_squared(visceralfat.aov)

                   as.factor(sweetPref)                      as.factor(fatPref) 
                            0.018365332                             0.038369994 
                                    age                          plantmilk_freq 
                            0.038513969                             0.007113967 
as.factor(sweetPref):as.factor(fatPref) 
                            0.009845245 

## BMR
BMRfat.aov <- full.df %>%
      aov(BMR ~ as.factor(sweetPref) + as.factor(fatPref) + age + plantmilk_freq, .)
summary(BMRfat.aov)

                      Df  Sum Sq Mean Sq F value Pr(>F)  
as.factor(sweetPref)   1   18157   18157   0.610 0.4362  
as.factor(fatPref)     2  124126   62063   2.084 0.1281  
age                    1  100880  100880   3.388 0.0677 .
plantmilk_freq         1    1314    1314   0.044 0.8339  
Residuals            144 4287735   29776                 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
30 observations deleted due to missingness

Group difference, 3-way ANCOVA, anthropometrics

## Lean mass
leanmass.aov <- full.df %>%
  aov(lean_mass ~ sfblCat*healthyfatCat*unhealthyfatCat + age + plantmilk_freq, .)
summary(leanmass.aov)

                                       Df Sum Sq Mean Sq F value Pr(>F)  
sfblCat                                 2  116.6   58.31   2.774 0.0664 .
healthyfatCat                           2   15.8    7.92   0.377 0.6870  
unhealthyfatCat                         2    4.7    2.35   0.112 0.8943  
age                                     1  128.3  128.26   6.101 0.0149 *
plantmilk_freq                          1    0.1    0.10   0.005 0.9438  
sfblCat:healthyfatCat                   4  118.7   29.68   1.412 0.2342  
sfblCat:unhealthyfatCat                 4   63.6   15.89   0.756 0.5560  
healthyfatCat:unhealthyfatCat           4  170.1   42.53   2.023 0.0955 .
sfblCat:healthyfatCat:unhealthyfatCat   8  216.7   27.09   1.289 0.2556  
Residuals                             121 2543.9   21.02                 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
30 observations deleted due to missingness

eta_squared(leanmass.aov)

                              sfblCat                         healthyfatCat 
                         3.452001e-02                          4.686505e-03 
                      unhealthyfatCat                                   age 
                         1.392044e-03                          3.796242e-02 
                       plantmilk_freq                 sfblCat:healthyfatCat 
                         3.100631e-05                          3.513683e-02 
              sfblCat:unhealthyfatCat         healthyfatCat:unhealthyfatCat 
                         1.881507e-02                          5.034900e-02 
sfblCat:healthyfatCat:unhealthyfatCat 
                         6.415411e-02 

## SMM
SMM.aov <- full.df %>%
  aov(SMM ~ sfblCat*healthyfatCat*unhealthyfatCat + age + plantmilk_freq, .)
summary(SMM.aov)

                                       Df Sum Sq Mean Sq F value Pr(>F)  
sfblCat                                 2    534   267.2   2.501 0.0863 .
healthyfatCat                           2     94    46.9   0.439 0.6460  
unhealthyfatCat                         2     16     8.0   0.075 0.9279  
age                                     1    659   659.3   6.171 0.0144 *
plantmilk_freq                          1      1     1.1   0.010 0.9188  
sfblCat:healthyfatCat                   4    610   152.5   1.428 0.2289  
sfblCat:unhealthyfatCat                 4    294    73.6   0.689 0.6011  
healthyfatCat:unhealthyfatCat           4    880   219.9   2.058 0.0905 .
sfblCat:healthyfatCat:unhealthyfatCat   8   1102   137.8   1.289 0.2552  
Residuals                             121  12928   106.8                 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
30 observations deleted due to missingness

eta_squared(SMM.aov)

                              sfblCat                         healthyfatCat 
                         3.121615e-02                          5.473900e-03 
                      unhealthyfatCat                                   age 
                         9.352196e-04                          3.851186e-02 
                       plantmilk_freq                 sfblCat:healthyfatCat 
                         6.511145e-05                          3.564240e-02 
              sfblCat:unhealthyfatCat         healthyfatCat:unhealthyfatCat 
                         1.719565e-02                          5.138863e-02 
sfblCat:healthyfatCat:unhealthyfatCat 
                         6.438360e-02 

## SMM, weight adjusted
SMMwt.aov <- full.df %>%
    aov(SMM_wt ~ sfblCat*healthyfatCat*unhealthyfatCat + age + plantmilk_freq, .)
summary(SMMwt.aov)

                                       Df Sum Sq Mean Sq F value Pr(>F)  
sfblCat                                 2    0.3    0.15   0.009 0.9907  
healthyfatCat                           2   15.1    7.54   0.477 0.6219  
unhealthyfatCat                         2   51.9   25.96   1.641 0.1981  
age                                     1   58.9   58.94   3.725 0.0559 .
plantmilk_freq                          1   60.5   60.47   3.822 0.0529 .
sfblCat:healthyfatCat                   4   97.9   24.47   1.546 0.1931  
sfblCat:unhealthyfatCat                 4   20.3    5.07   0.320 0.8638  
healthyfatCat:unhealthyfatCat           4   17.9    4.48   0.283 0.8884  
sfblCat:healthyfatCat:unhealthyfatCat   8  221.1   27.64   1.747 0.0943 .
Residuals                             121 1914.5   15.82                 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
30 observations deleted due to missingness

eta_squared(SMMwt.aov)

                              sfblCat                         healthyfatCat 
                         0.0001203767                          0.0061379740 
                      unhealthyfatCat                                   age 
                         0.0211173084                          0.0239744699 
                       plantmilk_freq                 sfblCat:healthyfatCat 
                         0.0245992300                          0.0398105985 
              sfblCat:unhealthyfatCat         healthyfatCat:unhealthyfatCat 
                         0.0082486316                          0.0072891627 
sfblCat:healthyfatCat:unhealthyfatCat 
                         0.0899560964 

## FFM
FFM.aov <- full.df %>%
    aov(FFM ~ sfblCat*healthyfatCat*unhealthyfatCat + age + plantmilk_freq, .)
summary(FFM.aov)

                                       Df Sum Sq Mean Sq F value Pr(>F)  
sfblCat                                 2   1601   800.4   2.725 0.0695 .
healthyfatCat                           2    214   106.8   0.364 0.6959  
unhealthyfatCat                         2     42    20.8   0.071 0.9318  
age                                     1   1792  1792.0   6.101 0.0149 *
plantmilk_freq                          1      0     0.3   0.001 0.9735  
sfblCat:healthyfatCat                   4   1713   428.3   1.458 0.2191  
sfblCat:unhealthyfatCat                 4    845   211.3   0.720 0.5802  
healthyfatCat:unhealthyfatCat           4   2401   600.2   2.044 0.0925 .
sfblCat:healthyfatCat:unhealthyfatCat   8   3028   378.5   1.289 0.2555  
Residuals                             121  35538   293.7                 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
30 observations deleted due to missingness

eta_squared(FFM.aov)

                              sfblCat                         healthyfatCat 
                         3.393601e-02                          4.527565e-03 
                      unhealthyfatCat                                   age 
                         8.798484e-04                          3.798643e-02 
                       plantmilk_freq                 sfblCat:healthyfatCat 
                         6.881121e-06                          3.632064e-02 
              sfblCat:unhealthyfatCat         healthyfatCat:unhealthyfatCat 
                         1.791876e-02                          5.089422e-02 
sfblCat:healthyfatCat:unhealthyfatCat 
                         6.419422e-02 

## Trunk lean mass
trunklean.aov <- full.df %>%
    aov(trunk_lean_mass ~ sfblCat*healthyfatCat*unhealthyfatCat + age + plantmilk_freq, .)
summary(trunklean.aov)

                                       Df Sum Sq Mean Sq F value  Pr(>F)   
sfblCat                                 2    224   112.0   1.813 0.16755   
healthyfatCat                           2     47    23.3   0.377 0.68647   
unhealthyfatCat                         2      8     3.8   0.062 0.93967   
age                                     1    436   435.8   7.054 0.00897 **
plantmilk_freq                          1      2     2.1   0.034 0.85419   
sfblCat:healthyfatCat                   4    434   108.5   1.757 0.14200   
sfblCat:unhealthyfatCat                 4    132    33.0   0.535 0.71056   
healthyfatCat:unhealthyfatCat           4    550   137.5   2.226 0.07014 . 
sfblCat:healthyfatCat:unhealthyfatCat   8    727    90.9   1.470 0.17491   
Residuals                             121   7476    61.8                   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
30 observations deleted due to missingness

eta_squared(trunklean.aov)

                              sfblCat                         healthyfatCat 
                         0.0223255553                          0.0046465146 
                      unhealthyfatCat                                   age 
                         0.0007666182                          0.0434285564 
                       plantmilk_freq                 sfblCat:healthyfatCat 
                         0.0002088214                          0.0432638001 
              sfblCat:unhealthyfatCat         healthyfatCat:unhealthyfatCat 
                         0.0131652013                          0.0548162814 
sfblCat:healthyfatCat:unhealthyfatCat 
                         0.0724249743 

## BFM
BFMfat.aov <- full.df %>%
    aov(BFM ~ sfblCat*healthyfatCat*unhealthyfatCat + age + plantmilk_freq, .)
summary(BFMfat.aov)

                                       Df Sum Sq Mean Sq F value  Pr(>F)   
sfblCat                                 2    805     403   0.690 0.50339   
healthyfatCat                           2   1572     786   1.348 0.26370   
unhealthyfatCat                         2    933     467   0.800 0.45167   
age                                     1   4923    4923   8.442 0.00436 **
plantmilk_freq                          1   1750    1750   3.002 0.08573 . 
sfblCat:healthyfatCat                   4   3162     790   1.355 0.25349   
sfblCat:unhealthyfatCat                 4    235      59   0.101 0.98203   
healthyfatCat:unhealthyfatCat           4   3341     835   1.432 0.22736   
sfblCat:healthyfatCat:unhealthyfatCat   8  12441    1555   2.667 0.00987 **
Residuals                             121  70564     583                   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
30 observations deleted due to missingness

eta_squared(BFMfat.aov)

                              sfblCat                         healthyfatCat 
                          0.008073239                           0.015762453 
                      unhealthyfatCat                                   age 
                          0.009356850                           0.049364033 
                       plantmilk_freq                 sfblCat:healthyfatCat 
                          0.017552320                           0.031704164 
              sfblCat:unhealthyfatCat         healthyfatCat:unhealthyfatCat 
                          0.002355548                           0.033505184 
sfblCat:healthyfatCat:unhealthyfatCat 
                          0.124753855 

emmeans.df <- pairs(emmeans(BFMfat.aov, ~ sfblCat * healthyfatCat * unhealthyfatCat), adjust = "tukey") %>% as.data.frame()
emmeans(BFMfat.aov, ~ sfblCat * healthyfatCat * unhealthyfatCat)

 sfblCat healthyfatCat unhealthyfatCat emmean    SE  df lower.CL upper.CL
 Low     Low           Low               47.0 10.20 121    26.81     67.3
 Medium  Low           Low               83.3 14.20 121    55.23    111.5
 High    Low           Low               60.2  8.12 121    44.12     76.3
 Low     Medium        Low               53.0  9.14 121    34.89     71.1
 Medium  Medium        Low               44.3  8.77 121    26.97     61.7
 High    Medium        Low               91.5 17.30 121    57.22    125.7
 Low     High          Low               41.9 10.60 121    20.92     62.9
 Medium  High          Low               48.0 10.90 121    26.45     69.5
 High    High          Low               28.3 17.40 121    -6.16     62.8
 Low     Low           Medium            48.1  8.07 121    32.15     64.1
 Medium  Low           Medium            39.5  9.15 121    21.42     57.6
 High    Low           Medium            55.4 10.90 121    33.81     76.9
 Low     Medium        Medium            35.3  9.86 121    15.81     54.9
 Medium  Medium        Medium            77.3  9.87 121    57.73     96.8
 High    Medium        Medium            81.6 14.00 121    53.96    109.3
 Low     High          Medium            47.4 11.30 121    25.01     69.8
 Medium  High          Medium            43.1 10.90 121    21.63     64.6
 High    High          Medium            33.8 10.80 121    12.45     55.2
 Low     Low           High              42.1 17.30 121     7.80     76.4
 Medium  Low           High              33.0 10.90 121    11.42     54.5
 High    Low           High              65.7 11.10 121    43.86     87.6
 Low     Medium        High              35.9 14.00 121     8.20     63.5
 Medium  Medium        High              62.0  9.96 121    42.25     81.7
 High    Medium        High              46.4  7.66 121    31.27     61.6
 Low     High          High              49.8 12.30 121    25.38     74.3
 Medium  High          High              43.1  9.13 121    25.01     61.2
 High    High          High              49.0  8.12 121    32.96     65.1

Confidence level used: 0.95 

## PBF
PBFfat.aov <- full.df %>%
      aov(PBF ~ sfblCat*healthyfatCat*unhealthyfatCat + age + plantmilk_freq, .)
summary(PBFfat.aov)

                                       Df Sum Sq Mean Sq F value Pr(>F)  
sfblCat                                 2      3    1.28   0.024 0.9768  
healthyfatCat                           2     73   36.67   0.671 0.5129  
unhealthyfatCat                         2    136   68.09   1.247 0.2911  
age                                     1    281  281.50   5.154 0.0250 *
plantmilk_freq                          1    159  159.24   2.915 0.0903 .
sfblCat:healthyfatCat                   4    382   95.62   1.751 0.1433  
sfblCat:unhealthyfatCat                 4     45   11.31   0.207 0.9341  
healthyfatCat:unhealthyfatCat           4    110   27.54   0.504 0.7326  
sfblCat:healthyfatCat:unhealthyfatCat   8    912  114.04   2.088 0.0420 *
Residuals                             121   6609   54.62                 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
30 observations deleted due to missingness

eta_squared(PBFfat.aov)

                              sfblCat                         healthyfatCat 
                         0.0002948443                          0.0084191003 
                      unhealthyfatCat                                   age 
                         0.0156325845                          0.0323122865 
                       plantmilk_freq                 sfblCat:healthyfatCat 
                         0.0182783042                          0.0439049357 
              sfblCat:unhealthyfatCat         healthyfatCat:unhealthyfatCat 
                         0.0051913443                          0.0126471571 
sfblCat:healthyfatCat:unhealthyfatCat 
                         0.1047233866 

emmeans.df <- pairs(emmeans(PBFfat.aov, ~ sfblCat * healthyfatCat * unhealthyfatCat), adjust = "tukey") %>% as.data.frame()
emmeans(PBFfat.aov, ~ sfblCat * healthyfatCat * unhealthyfatCat)

 sfblCat healthyfatCat unhealthyfatCat emmean   SE  df lower.CL upper.CL
 Low     Low           Low               33.8 3.13 121     27.7     40.0
 Medium  Low           Low               41.0 4.35 121     32.4     49.6
 High    Low           Low               37.0 2.48 121     32.1     42.0
 Low     Medium        Low               36.3 2.80 121     30.8     41.9
 Medium  Medium        Low               31.7 2.68 121     26.4     37.0
 High    Medium        Low               45.6 5.29 121     35.1     56.0
 Low     High          Low               33.7 3.25 121     27.3     40.2
 Medium  High          Low               35.2 3.33 121     28.7     41.8
 High    High          Low               24.7 5.33 121     14.1     35.2
 Low     Low           Medium            33.7 2.47 121     28.8     38.5
 Medium  Low           Medium            31.0 2.80 121     25.5     36.6
 High    Low           Medium            34.8 3.33 121     28.2     41.4
 Low     Medium        Medium            29.4 3.02 121     23.5     35.4
 Medium  Medium        Medium            38.0 3.02 121     32.0     44.0
 High    Medium        Medium            40.5 4.28 121     32.0     48.9
 Low     High          Medium            34.8 3.46 121     27.9     41.6
 Medium  High          Medium            32.9 3.32 121     26.3     39.5
 High    High          Medium            28.9 3.31 121     22.4     35.5
 Low     Low           High              29.0 5.30 121     18.5     39.5
 Medium  Low           High              30.2 3.33 121     23.6     36.8
 High    Low           High              38.6 3.38 121     31.9     45.3
 Low     Medium        High              28.7 4.28 121     20.3     37.2
 Medium  Medium        High              37.8 3.05 121     31.7     43.8
 High    Medium        High              32.3 2.34 121     27.6     36.9
 Low     High          High              34.0 3.78 121     26.5     41.4
 Medium  High          High              32.8 2.80 121     27.2     38.3
 High    High          High              31.0 2.48 121     26.1     35.9

Confidence level used: 0.95 

## Trunk fat mass
trunkfat.aov <- full.df %>%
    aov(trunk_fat_mass ~ sfblCat*healthyfatCat*unhealthyfatCat + age + plantmilk_freq, .)
summary(trunkfat.aov)

                                       Df Sum Sq Mean Sq F value  Pr(>F)   
sfblCat                                 2    169    84.4   0.640 0.52916   
healthyfatCat                           2    206   103.0   0.781 0.46034   
unhealthyfatCat                         2    195    97.4   0.738 0.48011   
age                                     1   1281  1281.5   9.710 0.00229 **
plantmilk_freq                          1    308   308.3   2.336 0.12901   
sfblCat:healthyfatCat                   4    996   249.0   1.887 0.11711   
sfblCat:unhealthyfatCat                 4     69    17.1   0.130 0.97128   
healthyfatCat:unhealthyfatCat           4    748   187.0   1.417 0.23228   
sfblCat:healthyfatCat:unhealthyfatCat   8   2679   334.8   2.537 0.01373 * 
Residuals                             121  15969   132.0                   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
30 observations deleted due to missingness

eta_squared(trunkfat.aov)

                              sfblCat                         healthyfatCat 
                          0.007466046                           0.009110813 
                      unhealthyfatCat                                   age 
                          0.008614067                           0.056654525 
                       plantmilk_freq                 sfblCat:healthyfatCat 
                          0.013630414                           0.044027361 
              sfblCat:unhealthyfatCat         healthyfatCat:unhealthyfatCat 
                          0.003029634                           0.033077432 
sfblCat:healthyfatCat:unhealthyfatCat 
                          0.118422952 

emmeans.df <- pairs(emmeans(trunkfat.aov, ~ sfblCat * healthyfatCat * unhealthyfatCat), adjust = "tukey") %>% as.data.frame()

## Visceral fat mass
visceralfat.aov <- full.df %>%
    aov(visceral_fat ~ sfblCat*healthyfatCat*unhealthyfatCat + age + plantmilk_freq, .)
summary(visceralfat.aov)

                                       Df Sum Sq Mean Sq F value  Pr(>F)   
sfblCat                                 2   2737    1369   0.519 0.59635   
healthyfatCat                           2   2093    1046   0.397 0.67325   
unhealthyfatCat                         2   7962    3981   1.510 0.22503   
age                                     1  23311   23311   8.843 0.00355 **
plantmilk_freq                          1   9982    9982   3.787 0.05398 . 
sfblCat:healthyfatCat                   4  20414    5103   1.936 0.10877   
sfblCat:unhealthyfatCat                 4   1973     493   0.187 0.94474   
healthyfatCat:unhealthyfatCat           4  12595    3149   1.194 0.31679   
sfblCat:healthyfatCat:unhealthyfatCat   8  54468    6808   2.583 0.01223 * 
Residuals                             121 318977    2636                   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
30 observations deleted due to missingness

eta_squared(visceralfat.aov)

                              sfblCat                         healthyfatCat 
                          0.006022072                           0.004604400 
                      unhealthyfatCat                                   age 
                          0.017516801                           0.051287315 
                       plantmilk_freq                 sfblCat:healthyfatCat 
                          0.021962886                           0.044914110 
              sfblCat:unhealthyfatCat         healthyfatCat:unhealthyfatCat 
                          0.004340705                           0.027710329 
sfblCat:healthyfatCat:unhealthyfatCat 
                          0.119838581 

emmeans.df <- pairs(emmeans(visceralfat.aov, ~ sfblCat * healthyfatCat * unhealthyfatCat), adjust = "tukey") %>% as.data.frame()

## BMR
BMRfat.aov <- full.df %>%
      aov(BMR ~ sfblCat*healthyfatCat*unhealthyfatCat + age + plantmilk_freq, .)
summary(BMRfat.aov)

                                       Df  Sum Sq Mean Sq F value Pr(>F)  
sfblCat                                 2  154063   77031   2.730 0.0692 .
healthyfatCat                           2   20465   10233   0.363 0.6966  
unhealthyfatCat                         2    3996    1998   0.071 0.9317  
age                                     1  172733  172733   6.122 0.0147 *
plantmilk_freq                          1      36      36   0.001 0.9716  
sfblCat:healthyfatCat                   4  165521   41380   1.467 0.2165  
sfblCat:unhealthyfatCat                 4   80707   20177   0.715 0.5832  
healthyfatCat:unhealthyfatCat           4  229711   57428   2.035 0.0937 .
sfblCat:healthyfatCat:unhealthyfatCat   8  290937   36367   1.289 0.2555  
Residuals                             121 3414046   28215                 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
30 observations deleted due to missingness

Response Surface Methodology

#   filter(!sample %in% c('S0F0', 'S1F0', 'S2F0')) 

rsm.df <- oatmilk.liking.long %>%
  separate(sample, into = c("sugar.level", "fat.level"), sep = 2) %>%
  mutate(sugar.percent = case_when(
    sugar.level == "S0" ~ 0,
    sugar.level == "S1" ~ 14,
    sugar.level == "S2" ~ 28
  )) %>%
  mutate(fat.percent = case_when(
    fat.level == "F0" ~ 0,
    fat.level == "F1" ~ 2,
    fat.level == "F2" ~ 6,
    fat.level == "F3" ~ 12
    )) %>%
  select(c(sugar.percent, fat.percent, liking))

## Base RSM model
rsm.model <- rsm(liking ~ SO(sugar.percent, fat.percent), data = rsm.df)
summary(rsm.model)

Call:
rsm(formula = liking ~ SO(sugar.percent, fat.percent), data = rsm.df)

                            Estimate Std. Error t value  Pr(>|t|)    
(Intercept)               43.7416667  1.7826758 24.5371 < 2.2e-16 ***
sugar.percent             -1.0155909  0.2810951 -3.6130 0.0003120 ***
fat.percent                1.2492514  0.6738945  1.8538 0.0639531 .  
sugar.percent:fat.percent  0.0699972  0.0196481  3.5625 0.0003779 ***
sugar.percent^2            0.0385731  0.0089458  4.3119 1.717e-05 ***
fat.percent^2             -0.1625737  0.0418985 -3.8802 0.0001086 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Multiple R-squared:  0.09604,   Adjusted R-squared:  0.09324 
F-statistic:  34.3 on 5 and 1614 DF,  p-value: < 2.2e-16

Analysis of Variance Table

Response: liking
                                  Df Sum Sq Mean Sq F value    Pr(>F)
FO(sugar.percent, fat.percent)     2  52316 26158.2  45.729 < 2.2e-16
TWI(sugar.percent, fat.percent)    1  29501 29501.4  51.573 1.049e-12
PQ(sugar.percent, fat.percent)     2  16272  8136.0  14.223 7.530e-07
Residuals                       1614 923253   572.0                  
Lack of fit                        3  20435  6811.7  12.155 7.239e-08
Pure error                      1611 902818   560.4                  

Stationary point of response surface:
sugar.percent   fat.percent 
     8.096881      5.585192 

Eigenanalysis:
eigen() decomposition
$values
[1]  0.04448869 -0.16848927

$vectors
                    [,1]       [,2]
sugar.percent -0.9860143 -0.1666605
fat.percent   -0.1666605  0.9860143

persp(rsm.model, ~ sugar.percent + fat.percent, xlab = 
        c("Fat Concentration", "Sugar Concentration"), zlab = "Liking",  col = c("lightblue", "lightyellow", "pink"))

## Central Composite Design model 
oatmilk.coded <- coded.data(rsm.df, x1 ~(Sugar-110)/110, x2 ~(Fat-10)/10) 

Regression, fat as predictors

visceral.S1F1.lm <- lm(liking_S1F1 ~ visceral_fat + age + plantmilk_freq, full.df)
summary(visceral.S1F1.lm)

Call:
lm(formula = liking_S1F1 ~ visceral_fat + age + plantmilk_freq, 
    data = full.df)

Residuals:
    Min      1Q  Median      3Q     Max 
-41.134 -15.856  -1.394  11.344  61.773 

Coefficients:
               Estimate Std. Error t value Pr(>|t|)   
(Intercept)    32.35608   11.98948   2.699  0.00778 **
visceral_fat   -0.03759    0.03368  -1.116  0.26620   
age             0.08875    0.56534   0.157  0.87547   
plantmilk_freq  2.33937    1.04194   2.245  0.02626 * 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 22.05 on 146 degrees of freedom
  (30 observations deleted due to missingness)
Multiple R-squared:  0.04714,   Adjusted R-squared:  0.02756 
F-statistic: 2.408 on 3 and 146 DF,  p-value: 0.06959

visceral.S1F2.lm <- lm(liking_S1F2 ~ visceral_fat + age + plantmilk_freq, full.df)
summary(visceral.S1F2.lm)

Call:
lm(formula = liking_S1F2 ~ visceral_fat + age + plantmilk_freq, 
    data = full.df)

Residuals:
   Min     1Q Median     3Q    Max 
-50.35 -15.33   0.40  13.45  51.14 

Coefficients:
               Estimate Std. Error t value Pr(>|t|)    
(Intercept)    41.77491   11.67732   3.577 0.000471 ***
visceral_fat   -0.07318    0.03280  -2.231 0.027200 *  
age             0.05750    0.55062   0.104 0.916966    
plantmilk_freq  2.12163    1.01481   2.091 0.038291 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 21.48 on 146 degrees of freedom
  (30 observations deleted due to missingness)
Multiple R-squared:  0.06917,   Adjusted R-squared:  0.05005 
F-statistic: 3.617 on 3 and 146 DF,  p-value: 0.01475

PBF.S1F1.lm <- lm(liking_S1F1 ~ PBF + age + plantmilk_freq, full.df)
summary(PBF.S1F1.lm)

Call:
lm(formula = liking_S1F1 ~ PBF + age + plantmilk_freq, data = full.df)

Residuals:
    Min      1Q  Median      3Q     Max 
-44.217 -16.016  -1.066  10.880  62.096 

Coefficients:
               Estimate Std. Error t value Pr(>|t|)   
(Intercept)     41.1831    13.5138   3.047  0.00274 **
PBF             -0.3801     0.2398  -1.585  0.11513   
age              0.1017     0.5576   0.182  0.85549   
plantmilk_freq   2.2952     1.0356   2.216  0.02822 * 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 21.96 on 146 degrees of freedom
  (30 observations deleted due to missingness)
Multiple R-squared:  0.05526,   Adjusted R-squared:  0.03585 
F-statistic: 2.847 on 3 and 146 DF,  p-value: 0.03969

PBF.S1F2.lm <- lm(liking_S1F2 ~ PBF + age + plantmilk_freq, full.df)
summary(PBF.S1F2.lm)

Call:
lm(formula = liking_S1F2 ~ PBF + age + plantmilk_freq, data = full.df)

Residuals:
    Min      1Q  Median      3Q     Max 
-53.993 -15.320   0.042  13.153  50.142 

Coefficients:
                Estimate Std. Error t value Pr(>|t|)    
(Intercept)    53.793904  13.196588   4.076 7.49e-05 ***
PBF            -0.548132   0.234199  -2.340   0.0206 *  
age             0.007147   0.544511   0.013   0.9895    
plantmilk_freq  2.139319   1.011306   2.115   0.0361 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 21.44 on 146 degrees of freedom
  (30 observations deleted due to missingness)
Multiple R-squared:  0.07225,   Adjusted R-squared:  0.05318 
F-statistic:  3.79 on 3 and 146 DF,  p-value: 0.0118

BFM.S1F1.lm <- lm(liking_S1F1 ~ PBF + age + plantmilk_freq, full.df)
summary(PBF.S1F1.lm)

Call:
lm(formula = liking_S1F1 ~ PBF + age + plantmilk_freq, data = full.df)

Residuals:
    Min      1Q  Median      3Q     Max 
-44.217 -16.016  -1.066  10.880  62.096 

Coefficients:
               Estimate Std. Error t value Pr(>|t|)   
(Intercept)     41.1831    13.5138   3.047  0.00274 **
PBF             -0.3801     0.2398  -1.585  0.11513   
age              0.1017     0.5576   0.182  0.85549   
plantmilk_freq   2.2952     1.0356   2.216  0.02822 * 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 21.96 on 146 degrees of freedom
  (30 observations deleted due to missingness)
Multiple R-squared:  0.05526,   Adjusted R-squared:  0.03585 
F-statistic: 2.847 on 3 and 146 DF,  p-value: 0.03969

BFM.S1F2.lm <- lm(liking_S1F2 ~ PBF + age + plantmilk_freq, full.df)
summary(PBF.S1F2.lm)

Call:
lm(formula = liking_S1F2 ~ PBF + age + plantmilk_freq, data = full.df)

Residuals:
    Min      1Q  Median      3Q     Max 
-53.993 -15.320   0.042  13.153  50.142 

Coefficients:
                Estimate Std. Error t value Pr(>|t|)    
(Intercept)    53.793904  13.196588   4.076 7.49e-05 ***
PBF            -0.548132   0.234199  -2.340   0.0206 *  
age             0.007147   0.544511   0.013   0.9895    
plantmilk_freq  2.139319   1.011306   2.115   0.0361 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 21.44 on 146 degrees of freedom
  (30 observations deleted due to missingness)
Multiple R-squared:  0.07225,   Adjusted R-squared:  0.05318 
F-statistic:  3.79 on 3 and 146 DF,  p-value: 0.0118

Group differences, demographics by body fat groups

## age
age.visceral.aov <- full.df %>%
  aov(age ~ visceralCat, data = .)
summary(age.visceral.aov)

             Df Sum Sq Mean Sq F value Pr(>F)  
visceralCat   2   56.4   28.19   2.612 0.0768 .
Residuals   147 1586.3   10.79                 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
30 observations deleted due to missingness

eta_squared(age.visceral.aov)

visceralCat 
 0.03431762 

emmeans(age.visceral.aov, pairwise ~ visceralCat, adjust = "tukey")

$emmeans
 visceralCat emmean    SE  df lower.CL upper.CL
 Low           21.1 0.465 147     20.2     22.0
 Medium        21.2 0.465 147     20.2     22.1
 High          22.4 0.465 147     21.5     23.4

Confidence level used: 0.95 

$contrasts
 contrast      estimate    SE  df t.ratio p.value
 Low - Medium     -0.04 0.657 147  -0.061  0.9980
 Low - High       -1.32 0.657 147  -2.009  0.1136
 Medium - High    -1.28 0.657 147  -1.948  0.1290

P value adjustment: tukey method for comparing a family of 3 estimates 

emmeans(age.visceral.aov, ~visceralCat)

 visceralCat emmean    SE  df lower.CL upper.CL
 Low           21.1 0.465 147     20.2     22.0
 Medium        21.2 0.465 147     20.2     22.1
 High          22.4 0.465 147     21.5     23.4

Confidence level used: 0.95 

## BMI
BMI.visceral.aov <- full.df %>%
  aov(BMI ~ visceralCat, data = .)
summary(BMI.visceral.aov)

             Df Sum Sq Mean Sq F value Pr(>F)    
visceralCat   2   3651  1825.4   79.74 <2e-16 ***
Residuals   147   3365    22.9                   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
30 observations deleted due to missingness

eta_squared(BMI.visceral.aov)

visceralCat 
  0.5203476 

emmeans(BMI.visceral.aov, pairwise ~ visceralCat, adjust = "tukey")

$emmeans
 visceralCat emmean    SE  df lower.CL upper.CL
 Low           20.1 0.677 147     18.8     21.4
 Medium        23.5 0.677 147     22.2     24.9
 High          31.8 0.677 147     30.5     33.2

Confidence level used: 0.95 

$contrasts
 contrast      estimate    SE  df t.ratio p.value
 Low - Medium     -3.44 0.957 147  -3.594  0.0013
 Low - High      -11.75 0.957 147 -12.281  <.0001
 Medium - High    -8.31 0.957 147  -8.687  <.0001

P value adjustment: tukey method for comparing a family of 3 estimates 

emmeans(BMI.visceral.aov, ~visceralCat)

 visceralCat emmean    SE  df lower.CL upper.CL
 Low           20.1 0.677 147     18.8     21.4
 Medium        23.5 0.677 147     22.2     24.9
 High          31.8 0.677 147     30.5     33.2

Confidence level used: 0.95 

## WC
WC.visceral.aov <- full.df %>%
  aov(WC ~ visceralCat, data = .)
summary(WC.visceral.aov)

             Df Sum Sq Mean Sq F value Pr(>F)    
visceralCat   2   2524  1261.9   87.83 <2e-16 ***
Residuals   147   2112    14.4                   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
30 observations deleted due to missingness

eta_squared(WC.visceral.aov)

visceralCat 
  0.5444272 

emmeans(WC.visceral.aov, pairwise ~ visceralCat, adjust = "tukey")

$emmeans
 visceralCat emmean    SE  df lower.CL upper.CL
 Low           26.6 0.536 147     25.5     27.7
 Medium        29.7 0.536 147     28.7     30.8
 High          36.4 0.536 147     35.4     37.5

Confidence level used: 0.95 

$contrasts
 contrast      estimate    SE  df t.ratio p.value
 Low - Medium     -3.14 0.758 147  -4.138  0.0002
 Low - High       -9.84 0.758 147 -12.974  <.0001
 Medium - High    -6.70 0.758 147  -8.836  <.0001

P value adjustment: tukey method for comparing a family of 3 estimates 

emmeans(WC.visceral.aov, ~visceralCat)

 visceralCat emmean    SE  df lower.CL upper.CL
 Low           26.6 0.536 147     25.5     27.7
 Medium        29.7 0.536 147     28.7     30.8
 High          36.4 0.536 147     35.4     37.5

Confidence level used: 0.95 

## Education
table(full.df$visceralCat, full.df$education)

        
          4  5  6  8  9 10
  Low     0 11 21  3 13  2
  Medium  0  8 20  6 16  0
  High    0  9 21  5 11  4

chisq.test(full.df$visceralCat, full.df$education)

Warning in stats::chisq.test(x, y, ...): Chi-squared approximation may be
incorrect

    Pearson's Chi-squared test

data:  full.df$visceralCat and full.df$education
X-squared = 6.4823, df = 8, p-value = 0.5934

## Household income
table(full.df$visceralCat, full.df$hinc)

        
          1  2  3  4  5  6  7  8
  Low     1  2  3 14 12  5 10  3
  Medium  2  3  6 13 10  7  5  4
  High    0  6 11 10  6  3  5  9

chisq.test(full.df$visceralCat, full.df$hinc)

Warning in stats::chisq.test(x, y, ...): Chi-squared approximation may be
incorrect

    Pearson's Chi-squared test

data:  full.df$visceralCat and full.df$hinc
X-squared = 19.941, df = 14, p-value = 0.132

## PB-milk frequency
pbmilk.visceral.aov <- full.df %>%
  aov(age ~ visceralCat, data = .)
summary(age.visceral.aov)

             Df Sum Sq Mean Sq F value Pr(>F)  
visceralCat   2   56.4   28.19   2.612 0.0768 .
Residuals   147 1586.3   10.79                 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
30 observations deleted due to missingness

eta_squared(pbmilk.visceral.aov)

visceralCat 
 0.03431762 

emmeans(pbmilk.visceral.aov, pairwise ~ visceralCat, adjust = "tukey")

$emmeans
 visceralCat emmean    SE  df lower.CL upper.CL
 Low           21.1 0.465 147     20.2     22.0
 Medium        21.2 0.465 147     20.2     22.1
 High          22.4 0.465 147     21.5     23.4

Confidence level used: 0.95 

$contrasts
 contrast      estimate    SE  df t.ratio p.value
 Low - Medium     -0.04 0.657 147  -0.061  0.9980
 Low - High       -1.32 0.657 147  -2.009  0.1136
 Medium - High    -1.28 0.657 147  -1.948  0.1290

P value adjustment: tukey method for comparing a family of 3 estimates 

emmeans(pbmilk.visceral.aov, ~visceralCat)

 visceralCat emmean    SE  df lower.CL upper.CL
 Low           21.1 0.465 147     20.2     22.0
 Medium        21.2 0.465 147     20.2     22.1
 High          22.4 0.465 147     21.5     23.4

Confidence level used: 0.95 

## Emotional Eating
emo.visceral.aov <- full.df %>%
  aov(emotional ~ visceralCat, .)
summary(emo.visceral.aov)

             Df Sum Sq Mean Sq F value Pr(>F)
visceralCat   2   1260   630.0   1.554  0.215
Residuals   147  59601   405.4               
30 observations deleted due to missingness

## Uncontrolled Eating
unc.visceral.aov <- full.df %>%
  aov(uncontrolled ~ visceralCat , .)
summary(unc.visceral.aov)

             Df Sum Sq Mean Sq F value Pr(>F)
visceralCat   2    564   281.8   1.658  0.194
Residuals   147  24994   170.0               
30 observations deleted due to missingness

## Restrictive Eating
res.visceral.aov <- full.df %>%
  aov(restraint ~ visceralCat , .)
summary(res.visceral.aov)

             Df Sum Sq Mean Sq F value Pr(>F)
visceralCat   2   1744   872.2    1.97  0.143
Residuals   147  65099   442.8               
30 observations deleted due to missingness

demo <- full.df %>%
  group_by(fatPref) %>%
  summarize(meanAge = mean(age), sdAge = sd(age), meanBMI = mean(BMI), sdBMI = sd(BMI), meanHEI = mean(sHEI), sdHEI = sd(sHEI), meanplantbasedmilk = mean(plantmilk_freq), sdplantbasedmilk = sd(plantmilk_freq))

Group differences, visceral fat on diet and taste liking

## visceral fat on fat preference
fatpref.visceral.multinom <- full.df %>%
  multinom(fatPref ~ visceralCat + age + plantmilk_freq, data = .)

# weights:  18 (10 variable)
initial  value 164.791843 
iter  10 value 153.178315
final  value 153.124425 
converged

tidy(fatpref.visceral.multinom, exponentiate = FALSE)  # coefficients (log odds)

# A tibble: 10 × 6
   y.level term              estimate std.error statistic p.value
   <chr>   <chr>                <dbl>     <dbl>     <dbl>   <dbl>
 1 Medium  (Intercept)        0.451      1.52      0.297  0.766  
 2 Medium  visceralCatMedium -0.375      0.557    -0.673  0.501  
 3 Medium  visceralCatHigh   -0.889      0.551    -1.61   0.107  
 4 Medium  age                0.0176     0.0699    0.253  0.801  
 5 Medium  plantmilk_freq     0.0723     0.137     0.526  0.599  
 6 High    (Intercept)       -0.156      1.55     -0.101  0.920  
 7 High    visceralCatMedium -0.439      0.582    -0.754  0.451  
 8 High    visceralCatHigh   -0.651      0.570    -1.14   0.253  
 9 High    age               -0.00474    0.0715   -0.0663 0.947  
10 High    plantmilk_freq     0.383      0.138     2.78   0.00539

tidy(fatpref.visceral.multinom, exponentiate = TRUE)   # odds ratios

# A tibble: 10 × 6
   y.level term              estimate std.error statistic p.value
   <chr>   <chr>                <dbl>     <dbl>     <dbl>   <dbl>
 1 Medium  (Intercept)          1.57     1.52      0.297  0.766  
 2 Medium  visceralCatMedium    0.687    0.557    -0.673  0.501  
 3 Medium  visceralCatHigh      0.411    0.551    -1.61   0.107  
 4 Medium  age                  1.02     0.0699    0.253  0.801  
 5 Medium  plantmilk_freq       1.07     0.137     0.526  0.599  
 6 High    (Intercept)          0.855    1.55     -0.101  0.920  
 7 High    visceralCatMedium    0.645    0.582    -0.754  0.451  
 8 High    visceralCatHigh      0.522    0.570    -1.14   0.253  
 9 High    age                  0.995    0.0715   -0.0663 0.947  
10 High    plantmilk_freq       1.47     0.138     2.78   0.00539

## visceral fat on sweet preference
sweetpref.visceral.multinom <- full.df %>%
  multinom(sweetPref ~ visceralCat + age + plantmilk_freq, data = .)

# weights:  6 (5 variable)
initial  value 103.972077 
iter  10 value 73.186727
iter  10 value 73.186727
iter  10 value 73.186727
final  value 73.186727 
converged

tidy(sweetpref.visceral.multinom, exponentiate = FALSE)  # coefficients (log odds)

# A tibble: 5 × 6
  y.level term              estimate std.error statistic p.value
  <chr>   <chr>                <dbl>     <dbl>     <dbl>   <dbl>
1 1       (Intercept)         1.17      1.42       0.822  0.411 
2 1       visceralCatMedium   0.438     0.480      0.914  0.361 
3 1       visceralCatHigh     0.808     0.531      1.52   0.128 
4 1       age                 0.0222    0.0667     0.333  0.739 
5 1       plantmilk_freq     -0.216     0.117     -1.85   0.0644

tidy(sweetpref.visceral.multinom, exponentiate = TRUE)   # odds ratios

# A tibble: 5 × 6
  y.level term              estimate std.error statistic p.value
  <chr>   <chr>                <dbl>     <dbl>     <dbl>   <dbl>
1 1       (Intercept)          3.22     1.42       0.822  0.411 
2 1       visceralCatMedium    1.55     0.480      0.914  0.361 
3 1       visceralCatHigh      2.24     0.531      1.52   0.128 
4 1       age                  1.02     0.0667     0.333  0.739 
5 1       plantmilk_freq       0.806    0.117     -1.85   0.0644

## visceral fat on sHEI
sHEI.visceral.aov <- full.df %>%
  aov(sHEI ~ visceralCat, data = .)
summary(sHEI.visceral.aov)

             Df Sum Sq Mean Sq F value Pr(>F)
visceralCat   2    114   57.22   0.741  0.479
Residuals   147  11355   77.25               
30 observations deleted due to missingness

eta_squared(sHEI.visceral.aov)

visceralCat 
0.009977517 

emmeans(sHEI.visceral.aov, pairwise ~ visceralCat, adjust = "tukey")

$emmeans
 visceralCat emmean   SE  df lower.CL upper.CL
 Low           46.7 1.24 147     44.2     49.1
 Medium        44.7 1.24 147     42.2     47.1
 High          45.0 1.24 147     42.6     47.5

Confidence level used: 0.95 

$contrasts
 contrast      estimate   SE  df t.ratio p.value
 Low - Medium     2.018 1.76 147   1.148  0.4864
 Low - High       1.625 1.76 147   0.925  0.6256
 Medium - High   -0.392 1.76 147  -0.223  0.9729

P value adjustment: tukey method for comparing a family of 3 estimates 

emmeans(sHEI.visceral.aov, ~visceralCat)

 visceralCat emmean   SE  df lower.CL upper.CL
 Low           46.7 1.24 147     44.2     49.1
 Medium        44.7 1.24 147     42.2     47.1
 High          45.0 1.24 147     42.6     47.5

Confidence level used: 0.95 

## Graph

# Fat preference
fatpref.or <- tidy(fatpref.visceral.multinom, exponentiate = TRUE, conf.int = TRUE) %>%
  filter(term %in% c("visceralCatMedium", "visceralCatHigh")) %>%
  mutate(outcome = "Fat Preference",
         comparison = y.level)

# Sweet preference - manually set y.level to "High"
sweetpref.or <- tidy(sweetpref.visceral.multinom, exponentiate = TRUE, conf.int = TRUE) %>%
  filter(term %in% c("visceralCatMedium", "visceralCatHigh")) %>%
  mutate(outcome = "Sweet Preference",
         comparison = "High",          # manually set to "High"
         y.level = "High")             # fix y.level too

# Combine
combined.or <- bind_rows(fatpref.or, sweetpref.or)

# Relabel term for y-axis
combined.or <- combined.or %>%
  mutate(term_label = recode(term,
                             "visceralCatMedium" = "Medium Visceral Fat",
                             "visceralCatHigh" = "High Visceral Fat"))

fatpref.or <- fatpref.or %>%
  mutate(term_label = recode(term,
                             "visceralCatMedium" = "Medium Visceral Fat",
                             "visceralCatHigh" = "High Visceral Fat"))

# Forest plot
ggplot(fatpref.or, aes(x = estimate, y = term_label, color = y.level)) +
  geom_point(size = 3, position = position_dodge(width = 0.5)) +          # dodge points
  geom_errorbarh(aes(xmin = conf.low, xmax = conf.high), height = 0.2,
                 position = position_dodge(width = 0.5)) +                 # dodge error bars
  geom_vline(xintercept = 1, linetype = "dashed", color = "red") +
  labs(x = "Odds Ratio (95% CI)",
       y = " ",
       color = "Preference Level",
       title = " ") +
  theme_classic() +
  theme(strip.text = element_text(size = 12, face = "bold"),
        legend.position = "right", 
        axis.text.x = element_text(size = 14),       # x-axis tick labels
        axis.text.y = element_text(size = 12),       # y-axis tick labels
        axis.title.x = element_text(size = 16),      # x-axis title
        axis.title.y = element_text(size = 16),      # y-axis title
        legend.text = element_text(size = 12),       # legend labels
        legend.title = element_text(size = 12))

Warning: `geom_errorbarh()` was deprecated in ggplot2 4.0.0.
ℹ Please use the `orientation` argument of `geom_errorbar()` instead.

`height` was translated to `width`.

ggsave("visceral.taste.png", height = 5, width = 6)

`height` was translated to `width`.

Group differences, body fat percent on diet and taste liking

## body fat percentage on fat preference
fatpref.pbf.multinom <- full.df %>%
  multinom(fatPref ~ pbfCat + age + plantmilk_freq, data = .)

# weights:  18 (10 variable)
initial  value 164.791843 
iter  10 value 153.557175
final  value 153.349142 
converged

tidy(fatpref.pbf.multinom, exponentiate = FALSE)  # coefficients (log odds)

# A tibble: 10 × 6
   y.level term            estimate std.error statistic p.value
   <chr>   <chr>              <dbl>     <dbl>     <dbl>   <dbl>
 1 Medium  (Intercept)     0.126       1.50     0.0836  0.933  
 2 Medium  pbfCatMedium    0.242       0.558    0.434   0.664  
 3 Medium  pbfCatHigh     -0.412       0.526   -0.784   0.433  
 4 Medium  age             0.0153      0.0689   0.222   0.824  
 5 Medium  plantmilk_freq  0.0752      0.136    0.552   0.581  
 6 High    (Intercept)    -0.385       1.55    -0.249   0.804  
 7 High    pbfCatMedium   -0.0659      0.575   -0.115   0.909  
 8 High    pbfCatHigh     -0.660       0.548   -1.20    0.229  
 9 High    age            -0.000304    0.0713  -0.00427 0.997  
10 High    plantmilk_freq  0.382       0.137    2.78    0.00540

tidy(fatpref.pbf.multinom, exponentiate = TRUE)   # odds ratios

# A tibble: 10 × 6
   y.level term           estimate std.error statistic p.value
   <chr>   <chr>             <dbl>     <dbl>     <dbl>   <dbl>
 1 Medium  (Intercept)       1.13     1.50     0.0836  0.933  
 2 Medium  pbfCatMedium      1.27     0.558    0.434   0.664  
 3 Medium  pbfCatHigh        0.662    0.526   -0.784   0.433  
 4 Medium  age               1.02     0.0689   0.222   0.824  
 5 Medium  plantmilk_freq    1.08     0.136    0.552   0.581  
 6 High    (Intercept)       0.681    1.55    -0.249   0.804  
 7 High    pbfCatMedium      0.936    0.575   -0.115   0.909  
 8 High    pbfCatHigh        0.517    0.548   -1.20    0.229  
 9 High    age               1.000    0.0713  -0.00427 0.997  
10 High    plantmilk_freq    1.47     0.137    2.78    0.00540

## body fat percentage on sweet preference
sweetpref.pbf.multinom <- full.df %>%
  multinom(sweetPref ~ pbfCat + age + plantmilk_freq, data = .)

# weights:  6 (5 variable)
initial  value 103.972077 
iter  10 value 73.759439
iter  10 value 73.759438
iter  10 value 73.759438
final  value 73.759438 
converged

tidy(sweetpref.pbf.multinom, exponentiate = FALSE)  # coefficients (log odds)

# A tibble: 5 × 6
  y.level term           estimate std.error statistic p.value
  <chr>   <chr>             <dbl>     <dbl>     <dbl>   <dbl>
1 1       (Intercept)      1.04      1.44       0.720  0.471 
2 1       pbfCatMedium     0.393     0.487      0.808  0.419 
3 1       pbfCatHigh       0.558     0.515      1.08   0.278 
4 1       age              0.0336    0.0668     0.503  0.615 
5 1       plantmilk_freq  -0.225     0.117     -1.93   0.0536

tidy(sweetpref.pbf.multinom, exponentiate = TRUE)   # odds ratios

# A tibble: 5 × 6
  y.level term           estimate std.error statistic p.value
  <chr>   <chr>             <dbl>     <dbl>     <dbl>   <dbl>
1 1       (Intercept)       2.82     1.44       0.720  0.471 
2 1       pbfCatMedium      1.48     0.487      0.808  0.419 
3 1       pbfCatHigh        1.75     0.515      1.08   0.278 
4 1       age               1.03     0.0668     0.503  0.615 
5 1       plantmilk_freq    0.798    0.117     -1.93   0.0536

## body fat on sHEI
sHEI.pbf.aov <- full.df %>%
  aov(sHEI ~ pbfCat, data = .)
summary(sHEI.pbf.aov)

             Df Sum Sq Mean Sq F value Pr(>F)
pbfCat        2     44   22.02   0.283  0.754
Residuals   147  11426   77.73               
30 observations deleted due to missingness

eta_squared(sHEI.pbf.aov)

     pbfCat 
0.003840229 

emmeans(sHEI.pbf.aov, pairwise ~ pbfCat, adjust = "tukey")

$emmeans
 pbfCat emmean   SE  df lower.CL upper.CL
 Low      46.0 1.25 147     43.6     48.5
 Medium   45.6 1.25 147     43.1     48.1
 High     44.7 1.25 147     42.3     47.2

Confidence level used: 0.95 

$contrasts
 contrast      estimate   SE  df t.ratio p.value
 Low - Medium     0.456 1.76 147   0.259  0.9638
 Low - High       1.308 1.76 147   0.742  0.7392
 Medium - High    0.851 1.76 147   0.483  0.8795

P value adjustment: tukey method for comparing a family of 3 estimates 

emmeans(sHEI.pbf.aov, ~pbfCat)

 pbfCat emmean   SE  df lower.CL upper.CL
 Low      46.0 1.25 147     43.6     48.5
 Medium   45.6 1.25 147     43.1     48.1
 High     44.7 1.25 147     42.3     47.2

Confidence level used: 0.95 

## Graph

# Fat preference
fatpref.or <- tidy(fatpref.pbf.multinom, exponentiate = TRUE, conf.int = TRUE) %>%
  filter(term %in% c("pbfCatMedium", "pbfCatHigh")) %>%
  mutate(outcome = "Fat Preference",
         comparison = y.level)

# Sweet preference - manually set y.level to "High"
sweetpref.or <- tidy(sweetpref.pbf.multinom, exponentiate = TRUE, conf.int = TRUE) %>%
  filter(term %in% c("pbfCatMedium", "pbfCatHigh")) %>%
  mutate(outcome = "Sweet Preference",
         comparison = "High",          # manually set to "High"
         y.level = "High")             # fix y.level too

# Combine
combined.or <- bind_rows(fatpref.or, sweetpref.or)

# Relabel term for y-axis
fatpref.or <- fatpref.or %>%
  mutate(term_label = recode(term,
                             "pbfCatMedium" = "Medium Percentage Body Fat",
                             "pbfCatHigh" = "High Percentage Body Fat"))

# Forest plot
ggplot(fatpref.or, aes(x = estimate, y = term_label, color = y.level)) +
  geom_point(size = 3, position = position_dodge(width = 0.5)) +          # dodge points
  geom_errorbarh(aes(xmin = conf.low, xmax = conf.high), height = 0.2,
                 position = position_dodge(width = 0.5)) +                 # dodge error bars
  geom_vline(xintercept = 1, linetype = "dashed", color = "red") +
  labs(x = "Odds Ratio (95% CI)",
       y = " ",
       color = "Preference Level",
       title = " ") +
  theme_classic() +
  theme(strip.text = element_text(size = 12, face = "bold"),
        legend.position = "right",         
        axis.text.x = element_text(size = 14),       # x-axis tick labels
        axis.text.y = element_text(size = 12),       # y-axis tick labels
        axis.title.x = element_text(size = 16),      # x-axis title
        axis.title.y = element_text(size = 16),      # y-axis title
        legend.text = element_text(size = 12),       # legend labels
        legend.title = element_text(size = 12))

`height` was translated to `width`.

ggsave("pbf.taste.png", height = 5, width = 6)

`height` was translated to `width`.