Sensory Nutrition Lab R Group Starter Code

Load packages and call data

• Load packages and call Qualtrics data file
• Removed the participants who 1) took less than 15 minutes to complete the survey, and 2) responded that they like “seeing a mouse in their house”/like “being late for an important date”/like “the smell of garbage”/dislike “going on vacation”
• Abbreviations: BAS = body appreciation score; BAScat = body appreciation score category, by median; CR = cognitive restraint; CRcat = cognitive restraint, by tertiles; EBs = eating behaviors; EM = emotional eating; EMcat = emotional eating category, by tertiles; SFBL = sweet foods and beverages liking; SFBLcat = sweet foods and beverages liking, by tertiles; sHEI = health eating index; UE = uncontrolled eating; UEcat = uncontrolled eating category, by tertiles; UHF = unhealthy fat liking; UHFcat = unhealthy fat liking, by tertiles

pacman::p_load(corrplot, MASS, readxl, writexl, ggsci, ggpubr, ggrepel, msm, ltm, ppcor, nnet, rstatix, FactoMineR, factoextra, purrr, multcomp, tidyr, GGally, ggcorrplot, ppcor, tidyverse)

## load surveys 
setwd("/Users/maycheung/Documents/R/SNL_R/Data") ## Please note that you will need to set the working directory based on where your data is stored on your computer
EB.survey.df <- read_excel("EB.qualtrics_29MAY25.xlsx")
SL.survey.df <- read_excel("SL.qualtrics_29MAY25.xlsx")

## join surveys
common_cols <- intersect(names(EB.survey.df), names(SL.survey.df))

EB.survey.df <- EB.survey.df %>% select(all_of(common_cols))
SL.survey.df <- SL.survey.df %>% select(all_of(common_cols))

survey.df <- bind_rows(EB.survey.df, SL.survey.df)

survey.df <- survey.df %>% 
  filter(!row_number() %in% 1) %>%
  mutate(across(-c("DM-05", "DM-05_7_TEXT", "DM-06"), as.numeric)) %>%
    filter(Duration > 720) %>%
  filter("FL-61_3" >= 0 | "FL-62_3" >= 0 | "FL-65_3" <= 0 | "FL-66_3" >= 0 ) %>%
  rename(
    age = "DM-01", height = "DM-02", weight = "DM-03", sex = "DM-04", race = "DM-05", orace = "DM-05_7_TEXT", ethn = "DM-06", usb = "DM-07", education = "DM-09", pinc = "DM-10", hinc = "DM-11", exec = "DM-12", smoke = "DM-13", diag = "DM-14") %>%   
  mutate(usb = case_when(
    usb == 1 ~ "U.S. Born",
    usb == 2 ~ "Immigrant")) %>%
  mutate(race = case_when(
    str_detect(race, ",") ~ "Mixed/Other",
    race == 1 & ethn == 2 ~ "African",
    race == 2 & ethn == 2 ~ "Mixed/Other",
    race == 3 & ethn == 2 ~ "Asian",
    race == 4 & ethn == 2 ~ "European",
    race == 5 & ethn == 2 ~ "NorAfrican",
    race == 6 & ethn == 2 ~ "Islander",
    race == 7 & ethn == 2 ~ "Mixed/Other",
    race >= 1 & ethn == 1 ~ "Hispanic")) 

survey.df <- survey.df %>% 
  rownames_to_column(., var = "ID") %>%
  filter(sex == 2)

Decoding survey

All scores and dietary intake estimates are calculated in this section

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

## Three factor questionnaire
TFQ.df <- survey.df %>% 
  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")) 

## Adult eating behavior questionnaire
AEBQ.df <- survey.df %>% 
  select("ID", "AEBQ-01":"AEBQ-35") %>%
  mutate(EF = rowSums(select(., "AEBQ-01", "AEBQ-03", "AEBQ-04"), na.rm = TRUE)) %>%
  mutate(EOE = rowSums(select(., "AEBQ-05", "AEBQ-08", "AEBQ-10", "AEBQ-16"), na.rm = TRUE)) %>%
  mutate(EUE = rowSums(select(., "AEBQ-15", "AEBQ-20", "AEBQ-27", "AEBQ-35"), na.rm = TRUE)) %>%
  mutate(FF = rowSums(select(., "AEBQ-02", "AEBQ-07", "AEBQ-12", "AEBQ-19", "AEBQ-24"), na.rm = TRUE)) %>%
  mutate(FR = rowSums(select(., "AEBQ-13", "AEBQ-22", "AEBQ-33"), na.rm = TRUE)) %>%
  mutate(SE = rowSums(select(., "AEBQ-14", "AEBQ-25", "AEBQ-26", "AEBQ-29"), na.rm = TRUE)) %>%
  mutate(H = rowSums(select(., "AEBQ-06", "AEBQ-09", "AEBQ-28", "AEBQ-32", "AEBQ-34"), na.rm = TRUE)) %>%
  mutate(SR = rowSums(select(., "AEBQ-11", "AEBQ-23", "AEBQ-30", "AEBQ-31"), na.rm = TRUE)) %>%
  select(c("ID", "EF":"SR")) 

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

## sHEI score
HEI.df <- survey.df %>%
  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

survey.df$ID <- as.character(survey.df$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)
AEBQ.df$ID <- as.character(AEBQ.df$ID)

full.df <- survey.df  %>%
  filter(age >= 18) %>%
  mutate(BMI = (weight/(height^2)*703)) %>%
  filter(BMI < 60) %>%
  select(c("ID", "sex", "age", "BMI", "race", "usb", "education", "hinc", "exec")) %>%
  mutate(agecat = case_when(
    age >= 18 & age < 26 ~ "1",
    age >= 26 & age < 36 ~ "2",
    age >= 36 & age < 66 ~ "3",
    age >= 66 ~ "4")) %>%
  mutate(BMIcat = case_when(
      BMI < 18.5 ~ "under",
      BMI >= 18.5 & BMI  < 25 & race != "Asian" ~ "healthy",
      BMI >= 25 & BMI < 30 & race != "Asian" ~ "overweight",
      BMI > 30 & race != "Asian" ~ "obese",
      BMI >= 18.5 & BMI  < 23 & race == "Asian" ~ "healthy",
      BMI >= 23 & BMI < 25 & race == "Asian" ~ "overweight",
      BMI > 25 & race == "Asian" ~ "obese"))  %>%
  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(AEBQ.df, by = "ID") %>% 
  mutate(BAScat = ntile(BAS, 2)) %>%
  mutate(UEcat = ntile(uncontrolled, 3)) %>%
  mutate(CRcat = ntile(restraint, 3)) %>%
  mutate(EMcat = ntile(emotional, 3)) %>%
  mutate(SFBLcat = ntile(sfbl.liking, 3)) %>%
  mutate(SuFatcat = ntile(sufat.liking, 3)) %>%
  mutate(UHFcat = ntile(unhealthyfat.liking, 3)) %>%
  mutate(HFcat = ntile(healthyfat.liking, 3)) %>%
  mutate(sugarcat = ntile(sugar.intake, 3)) %>%
  mutate(DQ = ntile(sHEI, 3))

full.df$group <- paste(full.df$race, full.df$usb)

Demographics

## Counts - by migration
full.df %>% 
  count(group) 

# A tibble: 12 × 2
   group                     n
   <chr>                 <int>
 1 African Immigrant        47
 2 African U.S. Born        95
 3 Asian Immigrant          62
 4 Asian U.S. Born          66
 5 European Immigrant        7
 6 European U.S. Born       29
 7 Hispanic Immigrant       13
 8 Hispanic U.S. Born       19
 9 Mixed/Other Immigrant     3
10 Mixed/Other U.S. Born    11
11 NorAfrican Immigrant      4
12 NorAfrican U.S. Born     11

## Counts - by age
full.df %>% 
  group_by(agecat) %>%
  count() 

# A tibble: 4 × 2
# Groups:   agecat [4]
  agecat     n
  <chr>  <int>
1 1        279
2 2         52
3 3         34
4 4          2

Data distribution

## Histogram function
histfunc <- function(data, x, title) {
  ggplot(data = data, aes(x = {{x}})) +
    geom_histogram(aes(y=..density..), colour="black", fill="white")+
    geom_density(alpha=.2, fill="#FF6666") +
    theme_light()
}

## QQ plot function
qqfunc <- function(data, sample, title) {
  ggplot(data = data, aes(sample = {{sample}})) +
  stat_qq() +
  stat_qq_line() +
    ggtitle(title) +
    theme_light()
}

## Age
histfunc(full.df, age, "Age histogram") 

qqfunc(full.df, age, "Age QQ plot")

shapiro.test(log10(full.df$age)) 

    Shapiro-Wilk normality test

data:  log10(full.df$age)
W = 0.80782, p-value < 2.2e-16

## BMI
histfunc(full.df, BMI, "BMI histogram")

qqfunc(full.df, BMI, "BMI QQ plot") 

shapiro.test(full.df$BMI) 

    Shapiro-Wilk normality test

data:  full.df$BMI
W = 0.91554, p-value = 1.743e-13

## Body appreciation score
histfunc(full.df, BAS, "BAS histogram")

qqfunc(full.df, BAS, "BAS QQ plot") 

shapiro.test(full.df$BAS) 

    Shapiro-Wilk normality test

data:  full.df$BAS
W = 0.97869, p-value = 3.034e-05

## Uncontrolled eating
histfunc(full.df, uncontrolled, "Uncontrolled eating histogram")

qqfunc(full.df, uncontrolled, "Uncontrolled eating QQ plot") 

shapiro.test(full.df$uncontrolled) 

    Shapiro-Wilk normality test

data:  full.df$uncontrolled
W = 0.98701, p-value = 0.002269

## Restrictive eating
histfunc(full.df, restraint, "Restrictive eating histogram")

qqfunc(full.df, restraint, "Restrictive eating QQ plot") 

shapiro.test(full.df$restraint) 

    Shapiro-Wilk normality test

data:  full.df$restraint
W = 0.95009, p-value = 8.336e-10

## Emotional eating
histfunc(full.df, emotional, "Emotional eating histogram")

qqfunc(full.df, emotional, "Emotional eating QQ plot") 

shapiro.test(full.df$emotional) 

    Shapiro-Wilk normality test

data:  full.df$emotional
W = 0.9469, p-value = 3.32e-10

## Sweet foods and beverages liking 
histfunc(full.df, sfbl.liking, "Sweet liking histogram")

qqfunc(full.df, sfbl.liking, "Sweet liking index QQ plot") 

shapiro.test(full.df$sfbl.liking)

    Shapiro-Wilk normality test

data:  full.df$sfbl.liking
W = 0.97939, p-value = 4.215e-05

## Unhealthy fat liking
histfunc(full.df, unhealthyfat.liking, "Unhealthy fat liking histogram")

qqfunc(full.df, unhealthyfat.liking, "Unhealthy fat liking QQ plot") 

shapiro.test(full.df$unhealthyfat.liking) 

    Shapiro-Wilk normality test

data:  full.df$unhealthyfat.liking
W = 0.98085, p-value = 8.595e-05

## sHEI 
histfunc(full.df, sHEI, "Healthy eating index histogram")

qqfunc(full.df, sHEI, "Healthy eating index QQ plot") 

shapiro.test(full.df$sHEI) 

    Shapiro-Wilk normality test

data:  full.df$sHEI
W = 0.99297, p-value = 0.08353

Correlations - overall corrplot

## All corr DF
corr.df <- full.df %>%
  select(age, BMI, education, hinc, exec, BAS, uncontrolled, restraint, emotional, sugar.intake, sHEI) 

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

p.mat <- cor_pmat(corr.df) 

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

## Liking corr DF
liking.corr.df <- full.df %>%
  select(veg.liking, fruit.liking, saltyfat.liking, hfprotein.liking, alcohol.liking, wg.liking, carbs.liking, fat.liking, healthyfat.liking, unhealthyfat.liking, sfbl.liking, protein.liking, spicy.liking, sHEI) 

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

p.mat <- cor_pmat(liking.corr.df) 

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

Differences - by DQ

## Age by DQ
age.DQ.aov <- full.df %>%
  aov(age ~ as.factor(DQ), .)
summary(age.DQ.aov)

               Df Sum Sq Mean Sq F value Pr(>F)  
as.factor(DQ)   2    553  276.48   3.355  0.036 *
Residuals     364  30000   82.42                 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

TukeyHSD(age.DQ.aov)

  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = age ~ as.factor(DQ), data = .)

$`as.factor(DQ)`
        diff       lwr      upr     p adj
2-1 1.576969 -1.152934 4.306872 0.3634200
3-1 3.003199  0.273296 5.733101 0.0269961
3-2 1.426230 -1.309239 4.161698 0.4379872

## BMI by DQ
BMI.sHEI.aov <- full.df %>%
  aov(BMI ~ as.factor(DQ), .)
summary(BMI.sHEI.aov)

               Df Sum Sq Mean Sq F value Pr(>F)
as.factor(DQ)   2    171   85.29   1.449  0.236
Residuals     364  21426   58.86               

TukeyHSD(BMI.sHEI.aov)

  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = BMI ~ as.factor(DQ), data = .)

$`as.factor(DQ)`
         diff       lwr       upr     p adj
2-1 -1.631331 -3.938369 0.6757077 0.2205103
3-1 -0.507317 -2.814356 1.7997214 0.8628579
3-2  1.124014 -1.187728 3.4357556 0.4875028

## BAS by DQ
BAS.sHEI.aov <- full.df %>%
  aov(BAS ~ as.factor(DQ), .)
summary(BAS.sHEI.aov)

               Df Sum Sq Mean Sq F value Pr(>F)
as.factor(DQ)   2     87   43.48   0.622  0.537
Residuals     364  25443   69.90               

TukeyHSD(BAS.sHEI.aov)

  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = BAS ~ as.factor(DQ), data = .)

$`as.factor(DQ)`
         diff       lwr      upr     p adj
2-1 0.6564041 -1.857610 3.170418 0.8123122
3-1 1.1891910 -1.324823 3.703205 0.5065341
3-2 0.5327869 -1.986353 3.051927 0.8724448

## Uncontrolled eating by DQ
UE.sHEI.aov <- full.df %>%
  aov(uncontrolled ~ as.factor(DQ), .)
summary(UE.sHEI.aov)

               Df Sum Sq Mean Sq F value Pr(>F)
as.factor(DQ)   2    819   409.6    2.24  0.108
Residuals     364  66563   182.9               

TukeyHSD(UE.sHEI.aov)

  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = uncontrolled ~ as.factor(DQ), data = .)

$`as.factor(DQ)`
         diff       lwr       upr     p adj
2-1 -2.195418 -6.261731 1.8708945 0.4127486
3-1 -3.629844 -7.696157 0.4364683 0.0910836
3-2 -1.434426 -5.509029 2.6401766 0.6855900

## Restrictive by DQ
CR.sHEI.aov <- full.df %>%
  aov(restraint ~ as.factor(DQ), .)
summary(CR.sHEI.aov)

               Df Sum Sq Mean Sq F value  Pr(>F)   
as.factor(DQ)   2   5044  2522.2   7.008 0.00103 **
Residuals     364 131002   359.9                   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

TukeyHSD(CR.sHEI.aov)

  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = restraint ~ as.factor(DQ), data = .)

$`as.factor(DQ)`
        diff        lwr      upr     p adj
2-1 6.627349  0.9227595 12.33194 0.0179807
3-1 8.676529  2.9719399 14.38112 0.0011361
3-2 2.049180 -3.6670393  7.76540 0.6761188

p.sig <- full.df %>%
  tukey_hsd(restraint ~ as.factor(DQ)) %>% 
  add_significance() %>%
  rename(., DQ = term) 

full.df %>%
  group_by(DQ) %>%
  summarize(mean = mean(restraint), sd=sd(restraint), n = n()) %>%
  ggplot(aes(x = as.factor(DQ), y = mean, fill = as.factor(DQ))) +
  geom_bar(stat = "identity", position = position_dodge()) +
  geom_errorbar(aes(ymin = mean - sd/sqrt(n), ymax = mean + sd/sqrt(n), width = 0.25)) +
  theme_bw() +
  theme(legend.position = "none") +
  labs(x = "", y = "Restrictive eating score") +
      stat_pvalue_manual(p.sig, hide.ns = T, label = "p.adj.signif",label.size = 10, bracket.size = 1, tip.length = 0.05 , y.position = c(60, 70)) +
    scale_fill_npg() +
  ylim(0, 75) +
scale_x_discrete(labels = c("Low DQ", "Medium DQ", "High DQ")) +
    theme(axis.text = element_text(size = 24), axis.title = element_text(size = 24))

## Emotional eating by sHEI
EM.sHEI.aov <- full.df %>%
  aov(emotional ~ as.factor(DQ), .)
summary(EM.sHEI.aov)

               Df Sum Sq Mean Sq F value Pr(>F)
as.factor(DQ)   2   1627   813.3   2.128  0.121
Residuals     364 139136   382.2               

TukeyHSD(EM.sHEI.aov)

  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = emotional ~ as.factor(DQ), data = .)

$`as.factor(DQ)`
          diff        lwr      upr     p adj
2-1 -4.1000267  -9.979043 1.778989 0.2297055
3-1 -4.7489338 -10.627950 1.130082 0.1398438
3-2 -0.6489071  -6.539909 5.242095 0.9636384

## Sweet foods and beverages liking by sHEI
SFBL.sHEI.aov <- full.df %>%
  aov(sfbl.liking ~ as.factor(DQ), .)
summary(SFBL.sHEI.aov)

               Df Sum Sq Mean Sq F value   Pr(>F)    
as.factor(DQ)   2  12594    6297   7.963 0.000412 ***
Residuals     364 287864     791                     
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

TukeyHSD(SFBL.sHEI.aov)

  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = sfbl.liking ~ as.factor(DQ), data = .)

$`as.factor(DQ)`
          diff       lwr       upr     p adj
2-1  -5.022940 -13.47920  3.433324 0.3430707
3-1 -14.148688 -22.60495 -5.692424 0.0002902
3-2  -9.125748 -17.59925 -0.652244 0.0312795

p.sig <- full.df %>%
  tukey_hsd(sfbl.liking ~ as.factor(DQ)) %>% 
  add_significance() %>%
  rename(., DQ = term) 

full.df %>%
  group_by(DQ) %>%
  summarize(mean = mean(sfbl.liking), sd=sd(sfbl.liking), n = n()) %>%
  ggplot(aes(x = as.factor(DQ), y = mean, fill = as.factor(DQ))) +
  geom_bar(stat = "identity", position = position_dodge()) +
  geom_errorbar(aes(ymin = mean - sd/sqrt(n), ymax = mean + sd/sqrt(n), width = 0.25)) +
  theme_bw() +
  theme(legend.position = "none") +
  labs(x = "", y = "Sweet foods/beverages liking score") +
      stat_pvalue_manual(p.sig, hide.ns = T, label = "p.adj.signif",label.size = 10, bracket.size = 1, tip.length = 0.05 , y.position = c(60)) +
    scale_fill_npg() +
  ylim(0, 70) +
scale_x_discrete(labels = c("Low DQ", "Medium DQ", "High DQ")) +
      theme(axis.text = element_text(size = 24), axis.title = element_text(size = 24))

## Sweet and fatty foods liking by DQ
SuFat.DQ.aov <- full.df %>%
  aov(sufat.liking ~ as.factor(DQ), .)
summary(SuFat.DQ.aov)

               Df Sum Sq Mean Sq F value  Pr(>F)   
as.factor(DQ)   2   9551    4775   4.853 0.00832 **
Residuals     364 358175     984                   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

TukeyHSD(SuFat.DQ.aov)

  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = sufat.liking ~ as.factor(DQ), data = .)

$`as.factor(DQ)`
          diff       lwr       upr     p adj
2-1  -3.757497 -13.19012  5.675125 0.6169376
3-1 -12.198071 -21.63069 -2.765448 0.0070699
3-2  -8.440574 -17.89243  1.011279 0.0909176

full.df %>%
  group_by(DQ) %>%
  summarize(mean = mean(sufat.liking), sd=sd(sufat.liking), n = n()) %>%
  ggplot(aes(x = as.factor(DQ), y = mean, fill = as.factor(DQ))) +
  geom_bar(stat = "identity", position = position_dodge()) +
    geom_errorbar(aes(ymin = mean - sd/sqrt(n), ymax = mean + sd/sqrt(n), width = 0.25)) +
    stat_pvalue_manual(p.sig, hide.ns = F, label = "p.adj.signif", label.size = 8, bracket.size = 1, tip.length = 0.05, y.position = c(70, 80, 90)) +
    ylim(0, 95) +
  scale_x_discrete(labels = c("Low DQ", "Medium DQ", "High DQ")) +
  theme_bw()+
  labs(title = " ", x = " ", y = "Sweet and fatty liking") +
  theme(axis.text.x = element_text(size = 20), axis.text.y = element_text(size = 15), axis.title.y = element_text(size = 20), legend.position = "none") +
     scale_fill_npg() 

## Hight fat protein liking by DQ
HFP.sHEI.aov <- full.df %>%
  aov(hfprotein.liking ~ as.factor(DQ), .)
summary(HFP.sHEI.aov)

               Df Sum Sq Mean Sq F value  Pr(>F)    
as.factor(DQ)   2  22543   11272   9.581 8.8e-05 ***
Residuals     364 428214    1176                    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

TukeyHSD(HFP.sHEI.aov)

  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = hfprotein.liking ~ as.factor(DQ), data = .)

$`as.factor(DQ)`
          diff       lwr        upr     p adj
2-1 -10.143832 -20.45755  0.1698816 0.0550473
3-1 -19.172521 -29.48624 -8.8588069 0.0000472
3-2  -9.028689 -19.36343  1.3060525 0.1006310

p.sig <- full.df %>%
  tukey_hsd(hfprotein.liking ~ as.factor(DQ)) %>% 
  add_significance() %>%
  rename(., DQ = term) 

full.df %>%
  group_by(DQ) %>%
  summarize(mean = mean(hfprotein.liking), sd=sd(hfprotein.liking), n = n()) %>%
  ggplot(aes(x = as.factor(DQ), y = mean, fill = as.factor(DQ))) +
  geom_bar(stat = "identity", position = position_dodge()) +
  geom_errorbar(aes(ymin = mean - sd/sqrt(n), ymax = mean + sd/sqrt(n), width = 0.25)) +
  theme_bw() +
  theme(legend.position = "none") +
  labs(x = "", y = "High fat protein liking score") +
      stat_pvalue_manual(p.sig, hide.ns = T, label = "p.adj.signif",label.size = 10, bracket.size = 1, tip.length = 0.05 , y.position = c(65)) +
    scale_fill_npg() +
  ylim(0, 75) +
scale_x_discrete(labels = c("Low DQ", "Medium DQ", "High DQ")) +
      theme(axis.text = element_text(size = 24), axis.title = element_text(size = 24))

## Unhealthy fat by DQ
UHF.sHEI.aov <- full.df %>%
  aov(unhealthyfat.liking ~ as.factor(DQ), .)
summary(UHF.sHEI.aov)

               Df Sum Sq Mean Sq F value  Pr(>F)   
as.factor(DQ)   2   8593    4297   5.451 0.00465 **
Residuals     364 286927     788                   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

TukeyHSD(UHF.sHEI.aov)

  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = unhealthyfat.liking ~ as.factor(DQ), data = .)

$`as.factor(DQ)`
          diff       lwr       upr     p adj
2-1  -4.601539 -13.04403  3.840953 0.4058581
3-1 -11.756457 -20.19895 -3.313965 0.0032895
3-2  -7.154918 -15.61462  1.304787 0.1160097

full.df %>%
  group_by(DQ) %>%
  summarize(mean = mean(unhealthyfat.liking), sd=sd(unhealthyfat.liking), n = n()) %>%
  ggplot(aes(x = as.factor(DQ), y = mean, fill = as.factor(DQ))) +
  geom_bar(stat = "identity", position = position_dodge()) +
  labs(title = "Unhealthy fat liking and sHEI", x = "Diet Quality", y = "Unhealthy fat liking") +
  geom_errorbar(aes(ymin = mean - sd/sqrt(n), ymax = mean + sd/sqrt(n), width = 0.25)) +
  theme_bw() +
  scale_x_discrete(labels = c("Low DQ", "Medium DQ", "High DQ")) +
  scale_fill_npg() +
  ylim(0, 70) +
scale_x_discrete(labels = c("Low DQ", "Medium DQ", "High DQ")) +
    theme(axis.text = element_text(size = 24), axis.title = element_text(size = 24))

Scale for x is already present.
Adding another scale for x, which will replace the existing scale.

## Healthy fat by sHEI
HF.sHEI.aov <- full.df %>%
  aov(healthyfat.liking ~ as.factor(DQ), .)
summary(HF.sHEI.aov)

               Df Sum Sq Mean Sq F value   Pr(>F)    
as.factor(DQ)   2  16525    8262   9.246 0.000121 ***
Residuals     364 325267     894                     
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

TukeyHSD(HF.sHEI.aov)

  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = healthyfat.liking ~ as.factor(DQ), data = .)

$`as.factor(DQ)`
         diff        lwr      upr     p adj
2-1  7.982757 -1.0061086 16.97162 0.0933169
3-1 16.424014  7.4351482 25.41288 0.0000651
3-2  8.441257 -0.5659346 17.44845 0.0715820

p.sig <- full.df %>%
  tukey_hsd(healthyfat.liking ~ as.factor(DQ)) %>% 
  add_significance() %>%
  rename(., DQ = term) 

full.df %>%
  group_by(DQ) %>%
  summarize(mean = mean(healthyfat.liking), sd=sd(healthyfat.liking), n = n()) %>%
  ggplot(aes(x = as.factor(DQ), y = mean, fill = as.factor(DQ))) +
  geom_bar(stat = "identity", position = position_dodge()) +
  geom_errorbar(aes(ymin = mean - sd/sqrt(n), ymax = mean + sd/sqrt(n), width = 0.25)) +
  theme_bw() +
  theme(legend.position = "none") +
  labs(x = "", y = "Healthy fat liking score") +
      stat_pvalue_manual(p.sig, hide.ns = T, label = "p.adj.signif",label.size = 10, bracket.size = 1, tip.length = 0.05 , y.position = c(68)) +
    scale_fill_npg() +
  ylim(0, 70) +
scale_x_discrete(labels = c("Low DQ", "Medium DQ", "High DQ")) +
    theme(axis.text = element_text(size = 24), axis.title = element_text(size = 24))