Sweet Taste Preference and Eating Behaviors on Diet

Study objectives

• Goal: To understand the influence of sweet taste liking and eating behaviors on dietary pattern

  • Objective 1: To determine whether sweet taste liking and eating behaviors are associated with snack choice
  • Objective 2: To determine whether sweet taste liking and eating behaviors are associated with added sugars intake
  • Objective 3: To determine whether sweet taste liking and eating behaviors are associated with diet quality (i.e., sHEI)

• Study design: Participants tasted a series of sucrose solutions to identify their preferred sucrose concentration (pref.mean). Participants also answered a series of surveys to assess their eating behaviors and diet quality (see survey variables below). Participant were asked to choose 1 snack (out of 12) to take home. The snacks can be divided into either 2 categories based on sweetness (not sweet vs sweet, SnackSwtCat), or caloric density (low vs high calories, SnackCalCat)

• Survey variable and their 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

Call data and decode surveys

setwd("/Users/maycheung/Documents/Analyses/SweetSnack")

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

## Demo df
raw.df <- read_excel("GSP.Master.xlsx") %>%
    filter(ID != "Test") %>%
  select(c("ID", "height.avg", "weight.avg", "BMI", "intensity", "liking", "pref.mean", "Snack", "SnackSwtCat", "SnackCalCat", "SnackCat")) %>% 
  mutate(PrefCat = case_when(
      pref.mean <= 0.35  ~ "disliker",
      pref.mean > 0.35 ~ "liker")) %>%
  mutate(across(-c("ID", "Snack", "SnackSwtCat", "SnackCalCat", "SnackCat", "PrefCat"), as.numeric))


## Survey df
survey.df <- read_excel("GSP.Qualtrics.xlsx") %>%
  filter(!row_number() %in% 1) %>%
  mutate(across(-c("DM-00", "DM-04", "DM-05", "DM-05_7_TEXT", "DM-06"), as.numeric)) %>%
  rename(ID = "DM-00", age = "DM-01", sex = "DM-04", race = "DM-05", orace = "DM-05_7_TEXT", ethn = "DM-06", education = "DM-09", hinc = "DM-11", pinc = "DM-12", smoke = "DM-13", exec = "DM-14") %>%  
  filter(!ID %in% c("Test-01", "Test-02", "Test-03")) %>%
  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")) 

## 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")))) %>%
  mutate(DQ = ntile(sHEI, 2)) 

## 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(restraint = rowSums(select(., "TFQ-01", "TFQ-05", "TFQ-11", "TFQ-17", "TFQ-18", "TFQ-21"), na.rm = TRUE)) %>%
  mutate(emotional = rowSums(select(., "TFQ-02", "TFQ-04", "TFQ-07", "TFQ-10", "TFQ-14", "TFQ-16"), 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")) 

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

## Putting scores together

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

full.df <- survey.df  %>%
  select(c("ID", "age", "race", "sex", "education", "hinc", "exec")) %>%
  left_join(raw.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") %>% 
  left_join(IES.df, by = "ID") %>% 
  mutate(SugarCat = ntile(sugar.intake, 3)) %>%
  mutate(BAScat = ntile(BAS, 2)) %>%
  mutate(UEcat = ntile(uncontrolled, 3)) %>%
  mutate(CRcat = ntile(restraint, 3)) %>%
  mutate(EMcat = ntile(emotional, 3)) %>%
  mutate(IEScat = ntile(IES, 2)) %>%
  mutate(SnackSwtCat = case_when(
      SnackSwtCat == "LS" ~ 1,
      SnackSwtCat == "HS" ~ 2)) %>%
    mutate(SnackCalCat = case_when(
      SnackCalCat == "LC" ~ 1,
      SnackCalCat == "HC" ~ 2)) %>%
  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"))  

Functions

• Useful functions (this part is hidden)

Data distribution

• I checked for data normality and I tried to transform the variables (e.g., log transform). It improved a little but it was not much. I recall there is a method called box-cox transformation but I have forgotten how to calculate the lambda (?).

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

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

shapiro.test(log10(full.df$age)) ## not normal

    Shapiro-Wilk normality test

data:  log10(full.df$age)
W = 0.81804, p-value = 1.683e-07

## pref.mean
full.df <- full.df %>% mutate(logPref = log10(pref.mean))
histfunc(full.df, logPref, "Sweet preference histogram") 

qqfunc(full.df, logPref, "Sweet preference QQ plot")

shapiro.test(full.df$logPref) ## not normal

    Shapiro-Wilk normality test

data:  full.df$logPref
W = 0.93718, p-value = 0.002565

## liking
full.df <- full.df %>% mutate(logLike = log10(liking+0.01))
histfunc(full.df, liking, "Liking histogram") 

qqfunc(full.df, liking, "Liking QQ plot")

shapiro.test(full.df$liking) ## not normal

    Shapiro-Wilk normality test

data:  full.df$liking
W = 0.95735, p-value = 0.02498

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

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

full.df <- full.df %>% mutate(logBMI = log10(BMI))
shapiro.test(full.df$logBMI) ## BMI not normal

    Shapiro-Wilk normality test

data:  full.df$logBMI
W = 0.93228, p-value = 0.001532

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

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

shapiro.test(full.df$BAS) ## normal

    Shapiro-Wilk normality test

data:  full.df$BAS
W = 0.97764, p-value = 0.2871

## IES eating
histfunc(full.df, IES, "Intuitive eating histogram")

qqfunc(full.df, IES, "Intuitive eating QQ plot") 

shapiro.test(full.df$IES) ## not normal

    Shapiro-Wilk normality test

data:  full.df$IES
W = 0.9619, p-value = 0.0431

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

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

shapiro.test(full.df$uncontrolled) ## normal

    Shapiro-Wilk normality test

data:  full.df$uncontrolled
W = 0.97901, p-value = 0.3353

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

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

shapiro.test(full.df$restraint) ## normal

    Shapiro-Wilk normality test

data:  full.df$restraint
W = 0.96796, p-value = 0.09

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

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

shapiro.test(full.df$emotional) ## not normal

    Shapiro-Wilk normality test

data:  full.df$emotional
W = 0.90727, p-value = 0.0001353

Counts

• To see how many people chose each snack category

## By sex
full.df %>%
  group_by(sex) %>%
  count()

# A tibble: 2 × 2
# Groups:   sex [2]
  sex       n
  <chr> <int>
1 1        28
2 2        37

## By race
full.df %>%
  group_by(race) %>%
  count()

# A tibble: 6 × 2
# Groups:   race [6]
  race            n
  <chr>       <int>
1 African        23
2 Asian          12
3 European       12
4 Hispanic        9
5 Mixed/Other     6
6 NorAfrican      3

## By hinc
full.df %>%
  group_by(hinc) %>%
  count()

# A tibble: 8 × 2
# Groups:   hinc [8]
   hinc     n
  <dbl> <int>
1     1     4
2     2     4
3     3     6
4     4    14
5     5    12
6     6    13
7     7     9
8     8     3

## By BMI
full.df %>%
  group_by(BMIcat) %>%
  count()

# A tibble: 4 × 2
# Groups:   BMIcat [4]
  BMIcat         n
  <chr>      <int>
1 healthy       26
2 obese         22
3 overweight    15
4 under          2

## By sweet liking phenotype
full.df %>%
  group_by(PrefCat) %>%
  count()

# A tibble: 2 × 2
# Groups:   PrefCat [2]
  PrefCat      n
  <chr>    <int>
1 disliker    32
2 liker       33

## By snack choice category
full.df %>%
  group_by(SnackCat) %>%
  count() %>%
  ungroup() %>%
  mutate(percent = (n/sum(n)) * 100)

# A tibble: 4 × 3
  SnackCat     n percent
  <chr>    <int>   <dbl>
1 HSHC        14   21.5 
2 HSLC        32   49.2 
3 LSHC        14   21.5 
4 LSLC         5    7.69

full.df %>%
  group_by(SnackSwtCat, SnackCalCat) %>%
  count() %>%
  ungroup() %>%
  mutate(percent = (n/sum(n)) * 100)

# A tibble: 4 × 4
  SnackSwtCat SnackCalCat     n percent
        <dbl>       <dbl> <int>   <dbl>
1           1           1     5    7.69
2           1           2    14   21.5 
3           2           1    30   46.2 
4           2           2    16   24.6 

chisq.test(full.df$SnackSwtCat, full.df$SnackCalCat)

    Pearson's Chi-squared test with Yates' continuity correction

data:  full.df$SnackSwtCat and full.df$SnackCalCat
X-squared = 6.6973, df = 1, p-value = 0.009656

full.df$SnackCalCat <- factor(full.df$SnackCalCat, levels = c("1", "2"), 
                  labels = c("Low calorie", "High calorie"))
  
full.df %>%
  group_by(SnackSwtCat, SnackCalCat) %>%
  count() %>%
  ggplot(aes(x = as.factor(SnackSwtCat), y = n, fill = as.factor(SnackSwtCat))) +
  geom_bar(stat = "identity", position = position_dodge()) +
  facet_wrap(. ~ as.factor(SnackCalCat)) +
  labs(x = " ", y = "Number of participants") +
  theme_bw() +
  theme(strip.text.x = element_text(size = 15, face = "bold")) + 
  theme(strip.text.x = element_text(size = 15, face = "bold")) + 
  theme(axis.text.x = element_text(size = 12)) +
  theme(axis.text.y = element_text(size = 15), axis.title=element_text(size=14)) +
  scale_x_discrete(labels=c("1" = "Low sweet", "2" = "High sweet")) +
  scale_fill_manual(values = c("#66BB6A", "#F44336"))

## By sugar cat
full.df %>%
  group_by(SugarCat) %>%
  select(c(SugarCat, sugar.intake)) %>%
  summarize(mean = mean(sugar.intake)/4, sd = sd(sugar.intake)/4)

# A tibble: 3 × 3
  SugarCat  mean    sd
     <int> <dbl> <dbl>
1        1  59.1  12.8
2        2 109.   18.9
3        3 206.   52.5

## By snack choice - individual
full.df %>%
  group_by(Snack) %>%
  count() %>%
  ungroup() %>%
  mutate(percent = (n/sum(n)) * 100) 

# A tibble: 10 × 3
   Snack                 n percent
   <chr>             <int>   <dbl>
 1 Almonds               1    1.54
 2 Apple                 1    1.54
 3 Beef/cheese stick     6    9.23
 4 Brownies              9   13.8 
 5 Cucumbers             5    7.69
 6 Fruit cup            27   41.5 
 7 Gummy bears           1    1.54
 8 Oreos                 6    9.23
 9 Potato Chips          4    6.15
10 Strawberries          5    7.69

full.df %>%
  group_by(Snack) %>%
  count() %>%
  ggplot(aes(Snack, n, fill = Snack)) +
  geom_bar(stat = "identity", position = position_dodge()) +
  labs(x = " ", y = "Number of participants ") +
        theme_bw() +
  ylim(0, 30)+
        theme(strip.text.x = element_text(size = 15, face = "bold")) + 
        theme(strip.text.x = element_text(size = 15, face = "bold")) + 
        theme(axis.text.x = element_text(angle = 45, size = 12, hjust = 1)) +
        theme(axis.text.y = element_text(size = 15), axis.title=element_text(size=14, vjust = 10)) +
  scale_fill_manual(values = c("#F57C00", "#B71C1C", "#4CC0E1", "#FDD835", "#00695C", "#66BB6A", "#283593", "#BDBDBD", "#8D6E63", "#E57373"))

Demographics

• Demographic table was generated and exported into excel (this part is hidden)

Correlations - overall corrplot

• Added sugars intake is associated with BMI, BAS, sweet intensity perception, and sweet preference
• sHEI is associated with household income and sweet preference

## Corr DF
corr.df <- full.df %>%
  select(age, BMI, hinc, uncontrolled, restraint, emotional, intensity, pref.mean, sugar.intake, sHEI) 

    ## All Participants
    all.corr <- corr.df %>%
      cor(use = "pairwise.complete.obs")
    
    all.res <- cor.mtest(corr.df, conf.level = 0.95)
    res.adj <- adjust_pvalue(all.res$p, method = "fdr")
    
    corrplot(all.corr, p.mat = res.adj, method = 'circle', tl.col = 'black', type = 'lower', insig='blank', addCoef.col ='black', number.cex = 0.7, order = 'original', diag=FALSE) 

Differences - by sweet preference

• This is not significant here, but I hypothesize (based on analyses from other independent data sets), that people who have a stronger sweet preference may exhibit more uncontrolled eating behaviors (binge eating)

## Sweet preference and age
t.test(age ~ as.factor(PrefCat), full.df)

    Welch Two Sample t-test

data:  age by as.factor(PrefCat)
t = -0.9142, df = 54.771, p-value = 0.3646
alternative hypothesis: true difference in means between group disliker and group liker is not equal to 0
95 percent confidence interval:
 -3.926927  1.466700
sample estimates:
mean in group disliker    mean in group liker 
              22.40625               23.63636 

full.df %>%
  group_by(PrefCat) %>%
  summarize(mean = mean(age), sd=sd(age), n = n()) %>%
  ggplot(aes(x = as.factor(PrefCat), y = mean, fill = as.factor(PrefCat))) +
  geom_bar(stat = "identity", position = position_dodge()) +
  labs(title = "Sweet preference and age", x = "Sweet preference", y = "age") +
  geom_errorbar(aes(ymin = mean - sd/sqrt(n), ymax = mean + sd/sqrt(n), width = 0.25)) +
  theme_bw() +
  scale_fill_manual(values = c("#66BB6A", "#F44336"))

## Sweet preference and BMI
t.test(BMI ~ as.factor(PrefCat), full.df)

    Welch Two Sample t-test

data:  BMI by as.factor(PrefCat)
t = -0.21944, df = 62.388, p-value = 0.827
alternative hypothesis: true difference in means between group disliker and group liker is not equal to 0
95 percent confidence interval:
 -4.590914  3.682584
sample estimates:
mean in group disliker    mean in group liker 
              27.69965               28.15382 

full.df %>%
  group_by(PrefCat) %>%
  summarize(mean = mean(BMI), sd=sd(BMI), n = n()) %>%
  ggplot(aes(x = as.factor(PrefCat), y = mean, fill = as.factor(PrefCat))) +
  geom_bar(stat = "identity", position = position_dodge()) +
  labs(title = "Sweet preference and BMI", x = "Sweet preference", y = "BMI") +
  geom_errorbar(aes(ymin = mean - sd/sqrt(n), ymax = mean + sd/sqrt(n), width = 0.25)) +
  theme_bw() +
  scale_fill_manual(values = c("#66BB6A", "#F44336"))

## Sweet preference and preferred sucrose concentration
t.test(pref.mean ~ as.factor(PrefCat), full.df)

    Welch Two Sample t-test

data:  pref.mean by as.factor(PrefCat)
t = -11.874, df = 45.036, p-value = 1.813e-15
alternative hypothesis: true difference in means between group disliker and group liker is not equal to 0
95 percent confidence interval:
 -0.5138395 -0.3648049
sample estimates:
mean in group disliker    mean in group liker 
             0.2024045              0.6417267 

## Sweet preference and intensity
t.test(intensity ~ as.factor(PrefCat), full.df)

    Welch Two Sample t-test

data:  intensity by as.factor(PrefCat)
t = 3.4558, df = 58.427, p-value = 0.001029
alternative hypothesis: true difference in means between group disliker and group liker is not equal to 0
95 percent confidence interval:
  4.955445 18.593797
sample estimates:
mean in group disliker    mean in group liker 
              79.56250               67.78788 

## Sweet preference and sugar intake
t.test(sugar.intake ~ as.factor(PrefCat), full.df)

    Welch Two Sample t-test

data:  sugar.intake by as.factor(PrefCat)
t = -2.7471, df = 55.348, p-value = 0.008098
alternative hypothesis: true difference in means between group disliker and group liker is not equal to 0
95 percent confidence interval:
 -308.61967  -48.28942
sample estimates:
mean in group disliker    mean in group liker 
              403.0000               581.4545 

full.df %>%
  group_by(PrefCat) %>%
  summarize(mean = mean(sugar.intake/4), sd=sd(sugar.intake/4), n = n()) %>%
  ggplot(aes(x = as.factor(PrefCat), y = mean, fill = as.factor(PrefCat))) +
  geom_bar(stat = "identity", position = position_dodge()) +
  geom_errorbar(aes(ymin = mean - sd/sqrt(n), ymax = mean + sd/sqrt(n), width = 0.25)) +
  labs(x = "", y = "Added sugars intake, grams") +
  theme(axis.text.x = element_text(size = 15)) +
  theme(axis.text.y = element_text(size = 15), axis.title=element_text(size=15)) +
  theme_bw() +
  scale_fill_manual(values = c("#66BB6A", "#F44336"))

## Sweet preference and sHEI
t.test(sHEI ~ as.factor(PrefCat), full.df)

    Welch Two Sample t-test

data:  sHEI by as.factor(PrefCat)
t = 1.3766, df = 58.098, p-value = 0.1739
alternative hypothesis: true difference in means between group disliker and group liker is not equal to 0
95 percent confidence interval:
 -1.340624  7.245378
sample estimates:
mean in group disliker    mean in group liker 
              46.10844               43.15606 

full.df %>%
  group_by(PrefCat) %>%
  summarize(mean = mean(sHEI), sd=sd(sHEI), n = n()) %>%
  ggplot(aes(x = as.factor(PrefCat), y = mean, fill = as.factor(PrefCat))) +
  geom_bar(stat = "identity", position = position_dodge()) +
  labs(title = "Sweet preference and sHEI", x = "Sweet preference", y = "sHEI") +
  geom_errorbar(aes(ymin = mean - sd/sqrt(n), ymax = mean + sd/sqrt(n), width = 0.25)) +
  theme_bw() +
  scale_fill_manual(values = c("#66BB6A", "#F44336"))

## Sweet preference and BAS
t.test(BAS ~ as.factor(PrefCat), full.df)

    Welch Two Sample t-test

data:  BAS by as.factor(PrefCat)
t = -0.12767, df = 62.43, p-value = 0.8988
alternative hypothesis: true difference in means between group disliker and group liker is not equal to 0
95 percent confidence interval:
 -4.337159  3.816326
sample estimates:
mean in group disliker    mean in group liker 
              34.40625               34.66667 

full.df %>%
  group_by(PrefCat) %>%
  summarize(mean = mean(BAS), sd=sd(BAS), n = n()) %>%
  ggplot(aes(x = as.factor(PrefCat), y = mean, fill = as.factor(PrefCat))) +
  geom_bar(stat = "identity", position = position_dodge()) +
  labs(title = "Sweet preference and BAS", x = "Sweet preference", y = "BAS") +
  geom_errorbar(aes(ymin = mean - sd/sqrt(n), ymax = mean + sd/sqrt(n), width = 0.25)) +
  theme_bw() +
  scale_fill_manual(values = c("#66BB6A", "#F44336"))

## Sweet preference and IES
t.test(IES ~ as.factor(PrefCat), full.df)

    Welch Two Sample t-test

data:  IES by as.factor(PrefCat)
t = 0.51995, df = 61.541, p-value = 0.605
alternative hypothesis: true difference in means between group disliker and group liker is not equal to 0
95 percent confidence interval:
 -0.2026562  0.3451134
sample estimates:
mean in group disliker    mean in group liker 
              3.480978               3.409750 

full.df %>%
  group_by(PrefCat) %>%
  summarize(mean = mean(IES), sd=sd(IES), n = n()) %>%
  ggplot(aes(x = as.factor(PrefCat), y = mean, fill = as.factor(PrefCat))) +
  geom_bar(stat = "identity", position = position_dodge()) +
  labs(title = "Sweet preference and IES", x = "Sweet preference", y = "IES") +
  geom_errorbar(aes(ymin = mean - sd/sqrt(n), ymax = mean + sd/sqrt(n), width = 0.25)) +
  theme_bw() +
  scale_fill_manual(values = c("#66BB6A", "#F44336"))

## Sweet preference and uncontrolled eating
t.test(uncontrolled ~ as.factor(PrefCat), full.df)

    Welch Two Sample t-test

data:  uncontrolled by as.factor(PrefCat)
t = -1.7291, df = 61.97, p-value = 0.08878
alternative hypothesis: true difference in means between group disliker and group liker is not equal to 0
95 percent confidence interval:
 -4.9063759  0.3552396
sample estimates:
mean in group disliker    mean in group liker 
              18.90625               21.18182 

full.df %>%
  group_by(PrefCat) %>%
  summarize(mean = mean(uncontrolled), sd=sd(uncontrolled), n = n()) %>%
  ggplot(aes(x = as.factor(PrefCat), y = mean, fill = as.factor(PrefCat))) +
  geom_bar(stat = "identity", position = position_dodge()) +
  labs(title = "Sweet preference and uncontrolled eating", x = "Sweet preference", y = "uncontrolled eating") +
  geom_errorbar(aes(ymin = mean - sd/sqrt(n), ymax = mean + sd/sqrt(n), width = 0.25)) +
  theme_bw() +
  scale_fill_manual(values = c("#66BB6A", "#F44336"))

## Sweet preference and restrictive eating
t.test(restraint ~ as.factor((PrefCat)), full.df)

    Welch Two Sample t-test

data:  restraint by as.factor((PrefCat))
t = -0.19961, df = 62.671, p-value = 0.8424
alternative hypothesis: true difference in means between group disliker and group liker is not equal to 0
95 percent confidence interval:
 -2.137822  1.749564
sample estimates:
mean in group disliker    mean in group liker 
              13.59375               13.78788 

full.df %>%
  group_by(PrefCat) %>%
  summarize(mean = mean(restraint), sd=sd(restraint), n = n()) %>%
  ggplot(aes(x = as.factor(PrefCat), y = mean, fill = as.factor(PrefCat))) +
  geom_bar(stat = "identity", position = position_dodge()) +
  labs(title = "Sweet preference and restrictive eating", x = "Sweet preference", y = "restrictive eating") +
  geom_errorbar(aes(ymin = mean - sd/sqrt(n), ymax = mean + sd/sqrt(n), width = 0.25)) +
  theme_bw() +
  scale_fill_manual(values = c("#66BB6A", "#F44336"))

## Sweet preference and emotional eating
t.test(emotional ~ as.factor(PrefCat), full.df)

    Welch Two Sample t-test

data:  emotional by as.factor(PrefCat)
t = 0.05463, df = 60.993, p-value = 0.9566
alternative hypothesis: true difference in means between group disliker and group liker is not equal to 0
95 percent confidence interval:
 -2.629777  2.777504
sample estimates:
mean in group disliker    mean in group liker 
              12.43750               12.36364 

full.df %>%
  group_by(PrefCat) %>%
  summarize(mean = mean(emotional), sd=sd(emotional), n = n()) %>%
  ggplot(aes(x = as.factor(PrefCat), y = mean, fill = as.factor(PrefCat))) +
  geom_bar(stat = "identity", position = position_dodge()) +
  labs(title = "Sweet preference and emotional eating", x = "Sweet preference", y = "emotional eating") +
  geom_errorbar(aes(ymin = mean - sd/sqrt(n), ymax = mean + sd/sqrt(n), width = 0.25)) +
  theme_bw() +
  scale_fill_manual(values = c("#66BB6A", "#F44336"))

Differences - by snack choice

## Snack choice and age
summary(aov(age ~ as.factor(SnackCat), full.df))

                    Df Sum Sq Mean Sq F value Pr(>F)
as.factor(SnackCat)  3     38   12.65   0.414  0.744
Residuals           61   1864   30.56               

## Snack choice and hinc
summary(aov(hinc ~ as.factor(SnackCat), full.df))

                    Df Sum Sq Mean Sq F value Pr(>F)
as.factor(SnackCat)  3    4.5   1.500   0.443  0.723
Residuals           61  206.5   3.385               

## Snack choice and sweet intensity
summary(aov(intensity ~ as.factor(SnackCat), full.df))

                    Df Sum Sq Mean Sq F value Pr(>F)
as.factor(SnackCat)  3    901   300.4   1.373   0.26
Residuals           61  13348   218.8               

## Snack choice and sweet preference
summary(aov(pref.mean ~ as.factor(SnackCat), full.df))

                    Df Sum Sq Mean Sq F value Pr(>F)
as.factor(SnackCat)  3  0.247 0.08239   1.164  0.331
Residuals           61  4.319 0.07080               

## Snack choice and BMI
summary(aov(BMI ~ as.factor(SnackCat), full.df))

                    Df Sum Sq Mean Sq F value Pr(>F)
as.factor(SnackCat)  3    214   71.49   1.041  0.381
Residuals           61   4191   68.71               

## Snack choice and sHEI
summary(aov(sHEI ~ as.factor(SnackCat), full.df))

                    Df Sum Sq Mean Sq F value Pr(>F)
as.factor(SnackCat)  3    460  153.37   2.109  0.108
Residuals           61   4437   72.74               

## Snack choice and BAS
summary(aov(BAS ~ as.factor(SnackCat), full.df))

                    Df Sum Sq Mean Sq F value Pr(>F)
as.factor(SnackCat)  3     99   33.10   0.483  0.695
Residuals           61   4177   68.47               

## Snack Choice and IES
summary(aov(IES ~ as.factor(SnackCat), full.df))

                    Df Sum Sq Mean Sq F value Pr(>F)
as.factor(SnackCat)  3  0.251 0.08373   0.267  0.849
Residuals           61 19.151 0.31394               

## Snack choice and uncontrolled eating
summary(aov(uncontrolled ~ as.factor(SnackCat), full.df))

                    Df Sum Sq Mean Sq F value Pr(>F)
as.factor(SnackCat)  3   16.1   5.375   0.177  0.911
Residuals           61 1849.6  30.322               

## Snack choice and restrictive eating
summary(aov(restraint ~ as.factor((SnackCat)), full.df))

                      Df Sum Sq Mean Sq F value Pr(>F)
as.factor((SnackCat))  3   47.6   15.87   1.047  0.378
Residuals             61  924.2   15.15               

## Snack choice and emotional eating
summary(aov(emotional ~ as.factor(SnackCat), full.df))

                    Df Sum Sq Mean Sq F value Pr(>F)
as.factor(SnackCat)  3   33.7   11.24   0.371  0.774
Residuals           61 1849.9   30.33               

Differences - by sugar intake

## Sugar intake and age
summary(aov(age ~ as.factor(SugarCat), full.df))

                    Df Sum Sq Mean Sq F value Pr(>F)
as.factor(SugarCat)  2     21   10.52   0.347  0.708
Residuals           62   1881   30.34               

full.df %>%
  group_by(SugarCat) %>%
  summarize(mean = mean(age), sd=sd(age), n = n()) %>%
  ggplot(aes(x = as.factor(SugarCat), y = mean, fill = as.factor(SugarCat))) +
  geom_bar(stat = "identity", position = position_dodge()) +
  labs(title = "Sugar intake and age", x = "Sugar intake", y = "age") +
  geom_errorbar(aes(ymin = mean - sd/sqrt(n), ymax = mean + sd/sqrt(n), width = 0.25)) +
  theme_bw() +
  scale_fill_manual(values = c("#66BB6A", "#FFCA28", "#F44336"))

## Sugar intake and BMI
summary(aov(BMI ~ as.factor(SugarCat), full.df))

                    Df Sum Sq Mean Sq F value Pr(>F)
as.factor(SugarCat)  2    131   65.49    0.95  0.392
Residuals           62   4275   68.94               

full.df %>%
  group_by(SugarCat) %>%
  summarize(mean = mean(BMI), sd=sd(BMI), n = n()) %>%
  ggplot(aes(x = as.factor(SugarCat), y = mean, fill = as.factor(SugarCat))) +
  geom_bar(stat = "identity", position = position_dodge()) +
  labs(title = "Sugar intake and BMI", x = "Sugar intake", y = "BMI") +
  geom_errorbar(aes(ymin = mean - sd/sqrt(n), ymax = mean + sd/sqrt(n), width = 0.25)) +
  theme_bw() +
  scale_fill_manual(values = c("#66BB6A", "#FFCA28", "#F44336"))

## Sugar intake and sweet preference
summary(aov(pref.mean ~ as.factor(SugarCat), full.df))

                    Df Sum Sq Mean Sq F value  Pr(>F)   
as.factor(SugarCat)  2  0.839  0.4196   6.981 0.00184 **
Residuals           62  3.726  0.0601                   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

full.df %>%
  group_by(SugarCat) %>%
  summarize(mean = mean(pref.mean), sd=sd(pref.mean), n = n()) %>%
  ggplot(aes(x = as.factor(SugarCat), y = mean, fill = as.factor(SugarCat))) +
  geom_bar(stat = "identity", position = position_dodge()) +
  labs(title = "Sugar intake and sweet preference", x = "Sugar intake", y = "Sweet preference") +
  geom_errorbar(aes(ymin = mean - sd/sqrt(n), ymax = mean + sd/sqrt(n), width = 0.25)) +
  theme_bw() +
  scale_fill_manual(values = c("#66BB6A", "#FFCA28", "#F44336"))

## Sugar intake and sHEI
summary(aov(sHEI ~ as.factor(SugarCat), full.df))

                    Df Sum Sq Mean Sq F value Pr(>F)
as.factor(SugarCat)  2     40   19.87   0.254  0.777
Residuals           62   4857   78.34               

full.df %>%
  group_by(SugarCat) %>%
  summarize(mean = mean(sHEI), sd=sd(sHEI), n = n()) %>%
  ggplot(aes(x = as.factor(SugarCat), y = mean, fill = as.factor(SugarCat))) +
  geom_bar(stat = "identity", position = position_dodge()) +
  labs(title = "Sugar intake and diet quality", x = "Sugar intake", y = "Diet quality") +
  geom_errorbar(aes(ymin = mean - sd/sqrt(n), ymax = mean + sd/sqrt(n), width = 0.25)) +
  theme_bw() +
  scale_fill_manual(values = c("#66BB6A", "#FFCA28", "#F44336"))

## Sugar intake and BAS
summary(aov(BAS ~ as.factor(SugarCat), full.df))

                    Df Sum Sq Mean Sq F value Pr(>F)  
as.factor(SugarCat)  2    475   237.6   3.876 0.0259 *
Residuals           62   3801    61.3                 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

full.df %>%
  group_by(SugarCat) %>%
  summarize(mean = mean(BAS), sd=sd(BAS), n = n()) %>%
  ggplot(aes(x = as.factor(SugarCat), y = mean, fill = as.factor(SugarCat))) +
  geom_bar(stat = "identity", position = position_dodge()) +
  labs(title = "Sugar intake and body appreciation", x = "Sugar intake", y = "Body appreciation") +
  geom_errorbar(aes(ymin = mean - sd/sqrt(n), ymax = mean + sd/sqrt(n), width = 0.25)) +
  theme_bw() +
  scale_fill_manual(values = c("#66BB6A", "#FFCA28", "#F44336"))

## Sugar intake and IES
summary(aov(IES ~ as.factor(SugarCat), full.df))

                    Df Sum Sq Mean Sq F value Pr(>F)
as.factor(SugarCat)  2  0.132 0.06608   0.213  0.809
Residuals           62 19.270 0.31080               

full.df %>%
  group_by(SugarCat) %>%
  summarize(mean = mean(IES), sd=sd(IES), n = n()) %>%
  ggplot(aes(x = as.factor(SugarCat), y = mean, fill = as.factor(SugarCat))) +
  geom_bar(stat = "identity", position = position_dodge()) +
  labs(title = "Sugar intake and intuitive eating", x = "Sugar intake", y = "Intuitive eating") +
  geom_errorbar(aes(ymin = mean - sd/sqrt(n), ymax = mean + sd/sqrt(n), width = 0.25)) +
  theme_bw() +
  scale_fill_manual(values = c("#66BB6A", "#FFCA28", "#F44336"))

## Sugar intake and uncontrolled eating
summary(aov(uncontrolled ~ as.factor(SugarCat), full.df))

                    Df Sum Sq Mean Sq F value Pr(>F)
as.factor(SugarCat)  2   36.9   18.44   0.625  0.539
Residuals           62 1828.9   29.50               

full.df %>%
  group_by(SugarCat) %>%
  summarize(mean = mean(uncontrolled), sd=sd(uncontrolled), n = n()) %>%
  ggplot(aes(x = as.factor(SugarCat), y = mean, fill = as.factor(SugarCat))) +
  geom_bar(stat = "identity", position = position_dodge()) +
  labs(title = "Sugar intake and uncontrolled eating", x = "Sugar intake", y = "Uncontrolled eating") +
  geom_errorbar(aes(ymin = mean - sd/sqrt(n), ymax = mean + sd/sqrt(n), width = 0.25)) +
  theme_bw() +
  scale_fill_manual(values = c("#66BB6A", "#FFCA28", "#F44336"))

## Sugar intake and restrictive eating
summary(aov(restraint ~ as.factor(SugarCat), full.df))

                    Df Sum Sq Mean Sq F value Pr(>F)
as.factor(SugarCat)  2    0.4   0.189   0.012  0.988
Residuals           62  971.5  15.669               

full.df %>%
  group_by(SugarCat) %>%
  summarize(mean = mean(restraint), sd=sd(restraint), n = n()) %>%
  ggplot(aes(x = as.factor(SugarCat), y = mean, fill = as.factor(SugarCat))) +
  geom_bar(stat = "identity", position = position_dodge()) +
  labs(title = "Sugar intake and restrictive eating", x = "Sugar intake", y = "Restrictive eating") +
  geom_errorbar(aes(ymin = mean - sd/sqrt(n), ymax = mean + sd/sqrt(n), width = 0.25)) +
  theme_bw() +
  scale_fill_manual(values = c("#66BB6A", "#FFCA28", "#F44336"))

## Sugar intake and emotional eating
summary(aov(emotional ~ as.factor(SugarCat), full.df))

                    Df Sum Sq Mean Sq F value Pr(>F)
as.factor(SugarCat)  2   49.1   24.54   0.829  0.441
Residuals           62 1834.5   29.59               

full.df %>%
  group_by(SugarCat) %>%
  summarize(mean = mean(emotional), sd=sd(emotional), n = n()) %>%
  ggplot(aes(x = as.factor(SugarCat), y = mean, fill = as.factor(SugarCat))) +
  geom_bar(stat = "identity", position = position_dodge()) +
  labs(title = "Sugar intake and emotional eating", x = "Sugar intake", y = "Emotional eating") +
  geom_errorbar(aes(ymin = mean - sd/sqrt(n), ymax = mean + sd/sqrt(n), width = 0.25)) +
  theme_bw() +
  scale_fill_manual(values = c("#66BB6A", "#FFCA28", "#F44336"))

Differences - by diet quality

## Diet quality and age
summary(aov(age ~ as.factor(DQ), full.df))

              Df Sum Sq Mean Sq F value Pr(>F)
as.factor(DQ)  1   10.4   10.38   0.346  0.559
Residuals     63 1891.6   30.02               

full.df %>%
  group_by(DQ) %>%
  summarize(mean = mean(age), sd=sd(age), n = n()) %>%
  ggplot(aes(x = as.factor(DQ), y = mean, fill = as.factor(DQ))) +
  geom_bar(stat = "identity", position = position_dodge()) +
  labs(title = "Diet quality and age", x = "Diet quality", y = "age") +
  geom_errorbar(aes(ymin = mean - sd/sqrt(n), ymax = mean + sd/sqrt(n), width = 0.25)) +
  theme_bw() +
  scale_fill_manual(values = c("#F44336", "#66BB6A"))

## Diet quality and BMI
summary(aov(BMI ~ as.factor(DQ), full.df))

              Df Sum Sq Mean Sq F value Pr(>F)
as.factor(DQ)  1     66   66.25   0.962   0.33
Residuals     63   4339   68.88               

full.df %>%
  group_by(DQ) %>%
  summarize(mean = mean(BMI), sd=sd(BMI), n = n()) %>%
  ggplot(aes(x = as.factor(DQ), y = mean, fill = as.factor(DQ))) +
  geom_bar(stat = "identity", position = position_dodge()) +
  labs(title = "Diet quality and BMI", x = "Diet quality", y = "BMI") +
  geom_errorbar(aes(ymin = mean - sd/sqrt(n), ymax = mean + sd/sqrt(n), width = 0.25)) +
  theme_bw() +
  scale_fill_manual(values = c("#F44336", "#66BB6A"))

## Diet quality and sweet preference
summary(aov(pref.mean ~ as.factor(DQ), full.df))

              Df Sum Sq Mean Sq F value Pr(>F)
as.factor(DQ)  1  0.027 0.02706   0.376  0.542
Residuals     63  4.539 0.07204               

full.df %>%
  group_by(DQ) %>%
  summarize(mean = mean(pref.mean), sd=sd(pref.mean), n = n()) %>%
  ggplot(aes(x = as.factor(DQ), y = mean, fill = as.factor(DQ))) +
  geom_bar(stat = "identity", position = position_dodge()) +
  labs(title = "Diet quality and sweet preference", x = "Diet quality", y = "Sweet preference") +
  geom_errorbar(aes(ymin = mean - sd/sqrt(n), ymax = mean + sd/sqrt(n), width = 0.25)) +
  theme_bw() +
  scale_fill_manual(values = c("#F44336", "#66BB6A"))

## Diet quality and sHEI
summary(aov(sHEI ~ as.factor(DQ), full.df))

              Df Sum Sq Mean Sq F value Pr(>F)    
as.factor(DQ)  1   3282    3282     128 <2e-16 ***
Residuals     63   1615      26                   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

full.df %>%
  group_by(DQ) %>%
  summarize(mean = mean(sHEI), sd=sd(sHEI), n = n()) %>%
  ggplot(aes(x = as.factor(DQ), y = mean, fill = as.factor(DQ))) +
  geom_bar(stat = "identity", position = position_dodge()) +
  labs(title = "Diet quality and diet quality", x = "Diet quality", y = "Diet quality") +
  geom_errorbar(aes(ymin = mean - sd/sqrt(n), ymax = mean + sd/sqrt(n), width = 0.25)) +
  theme_bw() +
  scale_fill_manual(values = c("#F44336", "#66BB6A"))

## Diet quality and BAS
summary(aov(BAS ~ as.factor(DQ), full.df))

              Df Sum Sq Mean Sq F value Pr(>F)
as.factor(DQ)  1    118   117.9   1.787  0.186
Residuals     63   4158    66.0               

full.df %>%
  group_by(DQ) %>%
  summarize(mean = mean(BAS), sd=sd(BAS), n = n()) %>%
  ggplot(aes(x = as.factor(DQ), y = mean, fill = as.factor(DQ))) +
  geom_bar(stat = "identity", position = position_dodge()) +
  labs(title = "Diet quality and body appreciation", x = "Diet quality", y = "Body appreciation") +
  geom_errorbar(aes(ymin = mean - sd/sqrt(n), ymax = mean + sd/sqrt(n), width = 0.25)) +
  theme_bw() +
  scale_fill_manual(values = c("#F44336", "#66BB6A"))

## Diet quality and IES
summary(aov(IES ~ as.factor(DQ), full.df))

              Df Sum Sq Mean Sq F value Pr(>F)
as.factor(DQ)  1  0.165  0.1646   0.539  0.466
Residuals     63 19.237  0.3054               

full.df %>%
  group_by(DQ) %>%
  summarize(mean = mean(IES), sd=sd(IES), n = n()) %>%
  ggplot(aes(x = as.factor(DQ), y = mean, fill = as.factor(DQ))) +
  geom_bar(stat = "identity", position = position_dodge()) +
  labs(title = "Diet quality and intuitive eating", x = "Diet quality", y = "Intuitive eating") +
  geom_errorbar(aes(ymin = mean - sd/sqrt(n), ymax = mean + sd/sqrt(n), width = 0.25)) +
  theme_bw() +
  scale_fill_manual(values = c("#F44336", "#66BB6A"))

## Diet quality and uncontrolled eating
summary(aov(uncontrolled ~ as.factor(DQ), full.df))

              Df Sum Sq Mean Sq F value Pr(>F)
as.factor(DQ)  1    6.1   6.117   0.207   0.65
Residuals     63 1859.6  29.518               

full.df %>%
  group_by(DQ) %>%
  summarize(mean = mean(uncontrolled), sd=sd(uncontrolled), n = n()) %>%
  ggplot(aes(x = as.factor(DQ), y = mean, fill = as.factor(DQ))) +
  geom_bar(stat = "identity", position = position_dodge()) +
  labs(title = "Diet quality and uncontrolled eating", x = "Diet quality", y = "Uncontrolled eating") +
  geom_errorbar(aes(ymin = mean - sd/sqrt(n), ymax = mean + sd/sqrt(n), width = 0.25)) +
  theme_bw() +
  scale_fill_manual(values = c("#F44336", "#66BB6A"))

## Diet quality and restrictive eating
summary(aov(restraint ~ as.factor(DQ), full.df))

              Df Sum Sq Mean Sq F value Pr(>F)
as.factor(DQ)  1    5.2   5.158   0.336  0.564
Residuals     63  966.7  15.344               

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()) +
  labs(title = "Diet quality and restrictive eating", x = "Diet quality", y = "Restrictive eating") +
  geom_errorbar(aes(ymin = mean - sd/sqrt(n), ymax = mean + sd/sqrt(n), width = 0.25)) +
  theme_bw() +
  scale_fill_manual(values = c("#F44336", "#66BB6A"))

## Diet quality and emotional eating
summary(aov(emotional ~ as.factor(DQ), full.df))

              Df Sum Sq Mean Sq F value Pr(>F)
as.factor(DQ)  1    0.9   0.889    0.03  0.864
Residuals     63 1882.7  29.884               

full.df %>%
  group_by(DQ) %>%
  summarize(mean = mean(emotional), sd=sd(emotional), n = n()) %>%
  ggplot(aes(x = as.factor(DQ), y = mean, fill = as.factor(DQ))) +
  geom_bar(stat = "identity", position = position_dodge()) +
  labs(title = "Diet quality and emotional eating", x = "Diet quality", y = "Emotional eating") +
  geom_errorbar(aes(ymin = mean - sd/sqrt(n), ymax = mean + sd/sqrt(n), width = 0.25)) +
  theme_bw() +
  scale_fill_manual(values = c("#F44336", "#66BB6A"))

Logistic regression - snack choice

## The effect of sweet preference on snack choice
    ## Sweet snack
    tidy(glm(as.factor(SnackSwtCat) ~ log10(pref.mean), family = binomial, full.df), exponentiate = TRUE, conf.int = TRUE)

# A tibble: 2 × 7
  term             estimate std.error statistic p.value conf.low conf.high
  <chr>               <dbl>     <dbl>     <dbl>   <dbl>    <dbl>     <dbl>
1 (Intercept)          2.58     0.501     1.89   0.0581    1.00       7.29
2 log10(pref.mean)     1.15     0.885     0.155  0.877     0.196      6.56

    ## Snack calories
    tidy(glm(as.factor(SnackCalCat) ~ log10(pref.mean), family = binomial, full.df), exponentiate = TRUE, conf.int = TRUE)

# A tibble: 2 × 7
  term             estimate std.error statistic p.value conf.low conf.high
  <chr>               <dbl>     <dbl>     <dbl>   <dbl>    <dbl>     <dbl>
1 (Intercept)          1.09     0.454     0.189   0.850    0.445      2.68
2 log10(pref.mean)     1.67     0.815     0.630   0.529    0.341      8.57

## The effect of BMI on snack choice
    ## Sweet snack
    tidy(glm(as.factor(SnackSwtCat) ~ BMI, family = binomial, full.df), exponentiate = TRUE, conf.int = TRUE)

# A tibble: 2 × 7
  term        estimate std.error statistic p.value conf.low conf.high
  <chr>          <dbl>     <dbl>     <dbl>   <dbl>    <dbl>     <dbl>
1 (Intercept)     1.54    0.997      0.431   0.666    0.192     10.4 
2 BMI             1.02    0.0349     0.470   0.638    0.953      1.10

    ## Snack calories
    tidy(glm(as.factor(SnackCalCat) ~ BMI, family = binomial, full.df), exponentiate = TRUE, conf.int = TRUE)

# A tibble: 2 × 7
  term        estimate std.error statistic p.value conf.low conf.high
  <chr>          <dbl>     <dbl>     <dbl>   <dbl>    <dbl>     <dbl>
1 (Intercept)    0.664    0.880     -0.465   0.642    0.114      3.79
2 BMI            1.01     0.0302     0.302   0.762    0.950      1.07

## The effect of restrictive eating on snack choice
    ## Sweet snack
    tidy(glm(as.factor(SnackSwtCat) ~ restraint, family = binomial, full.df), exponentiate = TRUE, conf.int = TRUE)

# A tibble: 2 × 7
  term        estimate std.error statistic p.value conf.low conf.high
  <chr>          <dbl>     <dbl>     <dbl>   <dbl>    <dbl>     <dbl>
1 (Intercept)    5.83     1.03       1.72   0.0860    0.819     48.3 
2 restraint      0.939    0.0704    -0.900  0.368     0.815      1.08

    ## Snack calories
    tidy(glm(as.factor(SnackCalCat) ~ restraint, family = binomial, full.df), exponentiate = TRUE, conf.int = TRUE)

# A tibble: 2 × 7
  term        estimate std.error statistic p.value conf.low conf.high
  <chr>          <dbl>     <dbl>     <dbl>   <dbl>    <dbl>     <dbl>
1 (Intercept)    0.537    0.921     -0.676   0.499   0.0840      3.24
2 restraint      1.03     0.0646     0.529   0.597   0.912       1.18

## The effect of uncontrolled eating on snack choice
    ## Sweet snack
    tidy(glm(as.factor(SnackSwtCat) ~ uncontrolled, family = binomial, full.df), exponentiate = TRUE, conf.int = TRUE)

# A tibble: 2 × 7
  term         estimate std.error statistic p.value conf.low conf.high
  <chr>           <dbl>     <dbl>     <dbl>   <dbl>    <dbl>     <dbl>
1 (Intercept)      1.42    1.06       0.333   0.739    0.175     11.7 
2 uncontrolled     1.03    0.0517     0.517   0.605    0.929      1.14

    ## Snack calories
    tidy(glm(as.factor(SnackCalCat) ~ uncontrolled, family = binomial, full.df), exponentiate = TRUE, conf.int = TRUE)

# A tibble: 2 × 7
  term         estimate std.error statistic p.value conf.low conf.high
  <chr>           <dbl>     <dbl>     <dbl>   <dbl>    <dbl>     <dbl>
1 (Intercept)     0.928    0.964    -0.0770   0.939    0.137      6.28
2 uncontrolled    0.996    0.0465   -0.0857   0.932    0.908      1.09

## The effect of emotional eating on snack choice
    ## Sweet snack
    tidy(glm(as.factor(SnackSwtCat) ~ emotional, family = binomial, full.df), exponentiate = TRUE, conf.int = TRUE)

# A tibble: 2 × 7
  term        estimate std.error statistic p.value conf.low conf.high
  <chr>          <dbl>     <dbl>     <dbl>   <dbl>    <dbl>     <dbl>
1 (Intercept)     1.36    0.689      0.450   0.653    0.348      5.35
2 emotional       1.05    0.0535     0.886   0.376    0.948      1.17

    ## Snack calories
    tidy(glm(as.factor(SnackCalCat) ~ emotional, family = binomial, full.df), exponentiate = TRUE, conf.int = TRUE)

# A tibble: 2 × 7
  term        estimate std.error statistic p.value conf.low conf.high
  <chr>          <dbl>     <dbl>     <dbl>   <dbl>    <dbl>     <dbl>
1 (Intercept)    1.43     0.630      0.562   0.574    0.416      5.03
2 emotional      0.960    0.0472    -0.874   0.382    0.872      1.05

Logistic regression - sugar intake

## The effect of sweet preference on sugar intake
full.df %>%
  ggplot(aes(pref.mean, SugarCat)) +
  geom_point() +
  geom_smooth(method = lm) +
  stat_cor() +
  theme_bw()

`geom_smooth()` using formula = 'y ~ x'

full.df$SugarCat <- as.factor(full.df$SugarCat)
summary(multinom(SugarCat ~ pref.mean, data = full.df, exponentiate = TRUE, conf.int = TRUE))

# weights:  9 (4 variable)
initial  value 71.409799 
final  value 64.743704 
converged

Call:
multinom(formula = SugarCat ~ pref.mean, data = full.df, exponentiate = TRUE, 
    conf.int = TRUE)

Coefficients:
  (Intercept) pref.mean
2   -1.132652  3.203846
3   -2.005743  4.857173

Std. Errors:
  (Intercept) pref.mean
2   0.5920126  1.469968
3   0.6825497  1.531867

Residual Deviance: 129.4874 
AIC: 137.4874 

tidy(multinom(SugarCat ~ pref.mean, data = full.df), exponentiate = TRUE, conf.int = TRUE)

# weights:  9 (4 variable)
initial  value 71.409799 
final  value 64.743704 
converged

# A tibble: 4 × 8
  y.level term        estimate std.error statistic p.value conf.low conf.high
  <chr>   <chr>          <dbl>     <dbl>     <dbl>   <dbl>    <dbl>     <dbl>
1 2       (Intercept)    0.322     0.592     -1.91 0.0557    0.101      1.03 
2 2       pref.mean     24.6       1.47       2.18 0.0293    1.38     439.   
3 3       (Intercept)    0.135     0.683     -2.94 0.00330   0.0353     0.513
4 3       pref.mean    129.        1.53       3.17 0.00152   6.39    2590.   

## The effect of intensity perception on sugar intake
full.df %>%
  ggplot(aes(intensity, SugarCat)) +
  geom_point() +
  geom_smooth(method = lm) +
  stat_cor() +
  theme_bw()

`geom_smooth()` using formula = 'y ~ x'

summary(multinom(SugarCat ~ intensity, data = full.df))

# weights:  9 (4 variable)
initial  value 71.409799 
iter  10 value 67.364423
iter  10 value 67.364423
final  value 67.364423 
converged

Call:
multinom(formula = SugarCat ~ intensity, data = full.df)

Coefficients:
  (Intercept)    intensity
2   0.6106706 -0.007915053
3   4.1239119 -0.057562602

Std. Errors:
  (Intercept)  intensity
2    1.859806 0.02378435
3    1.807014 0.02431695

Residual Deviance: 134.7288 
AIC: 142.7288 

tidy(multinom(SugarCat ~ intensity, data = full.df), exponentiate = TRUE, conf.int = TRUE)

# weights:  9 (4 variable)
initial  value 71.409799 
iter  10 value 67.364423
iter  10 value 67.364423
final  value 67.364423 
converged

# A tibble: 4 × 8
  y.level term        estimate std.error statistic p.value conf.low conf.high
  <chr>   <chr>          <dbl>     <dbl>     <dbl>   <dbl>    <dbl>     <dbl>
1 2       (Intercept)    1.84     1.86       0.328  0.743    0.0481    70.5  
2 2       intensity      0.992    0.0238    -0.333  0.739    0.947      1.04 
3 3       (Intercept)   61.8      1.81       2.28   0.0225   1.79    2134.   
4 3       intensity      0.944    0.0243    -2.37   0.0179   0.900      0.990

## The effect of BMI on sugar intake
full.df %>%
  ggplot(aes(BMI, SugarCat)) +
  geom_point() +
  geom_smooth(method = lm) +
  stat_cor() +
  theme_bw()

`geom_smooth()` using formula = 'y ~ x'

summary(multinom(SugarCat ~ BMI, data = full.df))

# weights:  9 (4 variable)
initial  value 71.409799 
final  value 70.457561 
converged

Call:
multinom(formula = SugarCat ~ BMI, data = full.df)

Coefficients:
  (Intercept)        BMI
2  -0.2907245 0.01078466
3  -1.3907493 0.04774725

Std. Errors:
  (Intercept)        BMI
2    1.151044 0.04122592
3    1.121589 0.03864272

Residual Deviance: 140.9151 
AIC: 148.9151 

tidy(multinom(SugarCat ~ BMI, data = full.df), exponentiate = TRUE, conf.int = TRUE)

# weights:  9 (4 variable)
initial  value 71.409799 
final  value 70.457561 
converged

# A tibble: 4 × 8
  y.level term        estimate std.error statistic p.value conf.low conf.high
  <chr>   <chr>          <dbl>     <dbl>     <dbl>   <dbl>    <dbl>     <dbl>
1 2       (Intercept)    0.748    1.15      -0.253   0.801   0.0783      7.14
2 2       BMI            1.01     0.0412     0.262   0.794   0.932       1.10
3 3       (Intercept)    0.249    1.12      -1.24    0.215   0.0276      2.24
4 3       BMI            1.05     0.0386     1.24    0.217   0.972       1.13

## The effect of BAS on sugar intake
full.df %>%
  ggplot(aes(BAS, SugarCat)) +
  geom_point() +
  geom_smooth(method = lm) +
  stat_cor() +
  theme_bw()

`geom_smooth()` using formula = 'y ~ x'

summary(multinom(SugarCat ~ BAS, data = full.df))

# weights:  9 (4 variable)
initial  value 71.409799 
final  value 67.573150 
converged

Call:
multinom(formula = SugarCat ~ BAS, data = full.df)

Coefficients:
  (Intercept)         BAS
2    2.318531 -0.06427259
3    3.860242 -0.11275701

Std. Errors:
  (Intercept)        BAS
2    1.503832 0.04078139
3    1.569354 0.04437967

Residual Deviance: 135.1463 
AIC: 143.1463 

tidy(multinom(SugarCat ~ BAS, data = full.df), exponentiate = TRUE, conf.int = TRUE)

# weights:  9 (4 variable)
initial  value 71.409799 
final  value 67.573150 
converged

# A tibble: 4 × 8
  y.level term        estimate std.error statistic p.value conf.low conf.high
  <chr>   <chr>          <dbl>     <dbl>     <dbl>   <dbl>    <dbl>     <dbl>
1 2       (Intercept)   10.2      1.50        1.54  0.123     0.533   194.   
2 2       BAS            0.938    0.0408     -1.58  0.115     0.866     1.02 
3 3       (Intercept)   47.5      1.57        2.46  0.0139    2.19   1029.   
4 3       BAS            0.893    0.0444     -2.54  0.0111    0.819     0.975

## Sugar log Mod 1 
sugar.logmod.1 <- multinom(SugarCat ~ pref.mean, full.df, exponentiate = TRUE, conf.int = TRUE)

# weights:  9 (4 variable)
initial  value 71.409799 
final  value 64.743704 
converged

tidy(sugar.logmod.1, exponentiate = TRUE, conf.int = TRUE)

# A tibble: 4 × 8
  y.level term        estimate std.error statistic p.value conf.low conf.high
  <chr>   <chr>          <dbl>     <dbl>     <dbl>   <dbl>    <dbl>     <dbl>
1 2       (Intercept)    0.322     0.592     -1.91 0.0557    0.101      1.03 
2 2       pref.mean     24.6       1.47       2.18 0.0293    1.38     439.   
3 3       (Intercept)    0.135     0.683     -2.94 0.00330   0.0353     0.513
4 3       pref.mean    129.        1.53       3.17 0.00152   6.39    2590.   

summary(sugar.logmod.1) # AIC = 137.5

Call:
multinom(formula = SugarCat ~ pref.mean, data = full.df, exponentiate = TRUE, 
    conf.int = TRUE)

Coefficients:
  (Intercept) pref.mean
2   -1.132652  3.203846
3   -2.005743  4.857173

Std. Errors:
  (Intercept) pref.mean
2   0.5920126  1.469968
3   0.6825497  1.531867

Residual Deviance: 129.4874 
AIC: 137.4874 

## Sugar log Mod 2 
sugar.logmod.2 <- multinom(SugarCat ~ pref.mean + intensity , full.df, exponentiate = TRUE, conf.int = TRUE)

# weights:  12 (6 variable)
initial  value 71.409799 
iter  10 value 62.770344
final  value 62.763444 
converged

tidy(sugar.logmod.2, exponentiate = TRUE, conf.int = TRUE)

# A tibble: 6 × 8
  y.level term        estimate std.error statistic p.value conf.low conf.high
  <chr>   <chr>          <dbl>     <dbl>     <dbl>   <dbl>    <dbl>     <dbl>
1 2       (Intercept)    0.134    2.28      -0.879 0.380    0.00153     11.8 
2 2       pref.mean     27.0      1.55       2.13  0.0335   1.29       564.  
3 2       intensity      1.01     0.0263     0.416 0.677    0.960        1.06
4 3       (Intercept)    1.94     2.26       0.292 0.770    0.0229     164.  
5 3       pref.mean     72.5      1.62       2.64  0.00833  3.01      1747.  
6 3       intensity      0.967    0.0272    -1.24  0.214    0.917        1.02

summary(sugar.logmod.2) # AIC = 137.5

Call:
multinom(formula = SugarCat ~ pref.mean + intensity, data = full.df, 
    exponentiate = TRUE, conf.int = TRUE)

Coefficients:
  (Intercept) pref.mean   intensity
2  -2.0065416  3.295694  0.01095777
3   0.6619808  4.283690 -0.03379538

Std. Errors:
  (Intercept) pref.mean  intensity
2    2.284022  1.550338 0.02632898
3    2.264425  1.623615 0.02719374

Residual Deviance: 125.5269 
AIC: 137.5269 

## Sugar log Mod 3
sugar.logmod.3 <- multinom(SugarCat ~ pref.mean + intensity + BAS, full.df, exponentiate = TRUE, conf.int = TRUE)

# weights:  15 (8 variable)
initial  value 71.409799 
iter  10 value 59.040765
final  value 58.994046 
converged

tidy(sugar.logmod.3)

# A tibble: 8 × 6
  y.level term        estimate std.error statistic p.value
  <chr>   <chr>          <dbl>     <dbl>     <dbl>   <dbl>
1 2       (Intercept)   0.906     2.80       0.323 0.747  
2 2       pref.mean     4.05      1.75       2.31  0.0207 
3 2       intensity     0.0118    0.0274     0.430 0.667  
4 2       BAS          -0.0897    0.0478    -1.88  0.0606 
5 3       (Intercept)   5.02      3.03       1.66  0.0979 
6 3       pref.mean     5.10      1.85       2.75  0.00593
7 3       intensity    -0.0337    0.0292    -1.15  0.249  
8 3       BAS          -0.133     0.0541    -2.45  0.0143 

summary(sugar.logmod.3) # AIC = 140.1

Call:
multinom(formula = SugarCat ~ pref.mean + intensity + BAS, data = full.df, 
    exponentiate = TRUE, conf.int = TRUE)

Coefficients:
  (Intercept) pref.mean   intensity         BAS
2   0.9057692  4.052565  0.01177329 -0.08966691
3   5.0207753  5.099192 -0.03365032 -0.13255490

Std. Errors:
  (Intercept) pref.mean  intensity        BAS
2    2.804485  1.751189 0.02738006 0.04779571
3    3.033200  1.853252 0.02921885 0.05411463

Residual Deviance: 117.9881 
AIC: 133.9881 

## Sugar log Mod 4 - with interactions 
sugar.logmod.4 <- multinom(SugarCat ~ pref.mean*intensity + BMI, full.df)

# weights:  18 (10 variable)
initial  value 71.409799 
iter  10 value 61.527396
iter  20 value 61.321607
iter  30 value 61.297906
iter  40 value 61.296459
final  value 61.296349 
converged

summary(sugar.logmod.4) # AIC = 142.6

Call:
multinom(formula = SugarCat ~ pref.mean * intensity + BMI, data = full.df)

Coefficients:
  (Intercept) pref.mean   intensity          BMI pref.mean:intensity
2    1.845888 -6.174594 -0.03727547 -0.002882568          0.12592011
3    1.922837 -1.640080 -0.06526975  0.040127527          0.08053323

Std. Errors:
  (Intercept) pref.mean  intensity        BMI pref.mean:intensity
2    4.248144  7.904272 0.04889250 0.04438127          0.10431813
3    4.394039  7.012175 0.05509319 0.04390515          0.09634732

Residual Deviance: 122.5927 
AIC: 142.5927 

## Visualization
ggcoef_multinom(sugar.logmod.1, variable_labels = c(pref.mean = "Preferred sucrose concentration, M"),
                y.level_label = c("3" = "High vs low sugar intake", "2" = "Medium vs low sugar intake")) +
labs(x = "Odds ratio)")

! `broom::tidy()` failed to tidy the model.
✔ `tidy_parameters()` used instead.
ℹ Add `tidy_fun = broom.helpers::tidy_parameters` to quiet these messages.

Linear regression - sugar intake

• Regression models were built using variables that are significantly correlated with added sugars intake (i.e., BMI, BAS, sweet intensity perception, and sweet preference)
• BMI is correlated with BAS
• Sweet intensity perception is correlated with sweet preference

## The effect of sweet preference on sugar intake
full.df %>%
  ggplot(aes(log10(pref.mean), sugar.intake)) +
  geom_point() +
  geom_smooth(method = lm) +
  stat_cor() +
  theme_bw()

tidy(glm(sugar.intake ~ log10(pref.mean), family = gaussian, full.df), exponentiate = TRUE, conf.int = TRUE)

# A tibble: 2 × 7
  term             estimate std.error statistic  p.value conf.low  conf.high
  <chr>               <dbl>     <dbl>     <dbl>    <dbl>    <dbl>      <dbl>
1 (Intercept)      4.25e300      55.5     12.5  1.13e-18 2.45e253 Inf       
2 log10(pref.mean) 6.88e183      98.9      4.28 6.50e- 5 4.33e 99   1.09e268

## The effect of intensity perception on sugar intake
full.df %>%
  ggplot(aes(intensity, sugar.intake)) +
  geom_point() +
  geom_smooth(method = lm) +
  stat_cor() +
  theme_bw()

tidy(glm(sugar.intake ~ intensity, family = gaussian, full.df), exponentiate = TRUE, conf.int = TRUE)

# A tibble: 2 × 7
  term         estimate std.error statistic     p.value  conf.low conf.high
  <chr>           <dbl>     <dbl>     <dbl>       <dbl>     <dbl>     <dbl>
1 (Intercept) Inf          165.        5.85 0.000000192 6.14e+277   Inf    
2 intensity     0.00172      2.19     -2.90 0.00507     2.34e-  5     0.126

## The effect of BMI on sugar intake
full.df %>%
  ggplot(aes(BMI, sugar.intake)) +
  geom_point() +
  geom_smooth(method = lm) +
  stat_cor() +
  theme_bw()

tidy(glm(sugar.intake ~ BMI, family = gaussian, full.df), exponentiate = TRUE, conf.int = TRUE)

# A tibble: 2 × 7
  term        estimate std.error statistic p.value conf.low conf.high
  <chr>          <dbl>     <dbl>     <dbl>   <dbl>    <dbl>     <dbl>
1 (Intercept)  1.07e87    116.        1.73  0.0890 1.92e-12  5.95e185
2 BMI          3.62e 4      3.98      2.64  0.0106 1.47e+ 1  8.92e  7

## The effect of BAS on sugar intake
full.df %>%
  ggplot(aes(BAS, sugar.intake)) +
  geom_point() +
  geom_smooth(method = lm) +
  stat_cor() +
  theme_bw()

tidy(glm(sugar.intake ~ BAS, family = gaussian, full.df), exponentiate = TRUE, conf.int = TRUE)

# A tibble: 2 × 7
  term             estimate std.error statistic  p.value   conf.low  conf.high
  <chr>               <dbl>     <dbl>     <dbl>    <dbl>      <dbl>      <dbl>
1 (Intercept) Inf              134.        7.78 8.55e-11 Inf        Inf       
2 BAS           0.000000135      3.77     -4.20 8.55e- 5   8.41e-11   0.000217

## Sugar lin Mod 1 
sugar.linmod.1 <- glm(sugar.intake ~ log10(pref.mean) , family = gaussian, full.df)
summary(sugar.linmod.1)

Call:
glm(formula = sugar.intake ~ log10(pref.mean), family = gaussian, 
    data = full.df)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-371.81  -178.10   -53.26   102.74   673.35  

Coefficients:
                 Estimate Std. Error t value Pr(>|t|)    
(Intercept)        692.22      55.50  12.473  < 2e-16 ***
log10(pref.mean)   423.30      98.92   4.279  6.5e-05 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for gaussian family taken to be 60149.45)

    Null deviance: 4890850  on 64  degrees of freedom
Residual deviance: 3789415  on 63  degrees of freedom
AIC: 903.73

Number of Fisher Scoring iterations: 2

## Sugar lin Mod 2 
sugar.linmod.2 <- glm(sugar.intake ~ log10(pref.mean) + intensity , family = gaussian, full.df)
summary(sugar.linmod.2)

Call:
glm(formula = sugar.intake ~ log10(pref.mean) + intensity, family = gaussian, 
    data = full.df)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-395.42  -147.77   -42.47   102.28   688.80  

Coefficients:
                 Estimate Std. Error t value Pr(>|t|)    
(Intercept)       901.924    153.923   5.860 1.92e-07 ***
log10(pref.mean)  357.667    107.879   3.315  0.00153 ** 
intensity          -3.268      2.241  -1.459  0.14969    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for gaussian family taken to be 59091.6)

    Null deviance: 4890850  on 64  degrees of freedom
Residual deviance: 3663679  on 62  degrees of freedom
AIC: 903.54

Number of Fisher Scoring iterations: 2

## Sugar lin Mod 3
sugar.linmod.3 <- glm(sugar.intake ~ log10(pref.mean) + intensity + BAS, family = gaussian, full.df)
summary(sugar.linmod.3)    

Call:
glm(formula = sugar.intake ~ log10(pref.mean) + intensity + BAS, 
    family = gaussian, data = full.df)

Deviance Residuals: 
   Min      1Q  Median      3Q     Max  
-371.4  -169.6   -26.1   148.8   627.8  

Coefficients:
                 Estimate Std. Error t value Pr(>|t|)    
(Intercept)      1271.988    163.148   7.797 9.85e-11 ***
log10(pref.mean)  359.385     95.944   3.746 0.000401 ***
intensity          -1.677      2.029  -0.826 0.411807    
BAS               -14.082      3.377  -4.170 9.80e-05 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for gaussian family taken to be 46739.57)

    Null deviance: 4890850  on 64  degrees of freedom
Residual deviance: 2851114  on 61  degrees of freedom
AIC: 889.24

Number of Fisher Scoring iterations: 2

## Sugar lin Mod 4 - with interactions 
sugar.linmod.4 <- glm(sugar.intake ~ log10(pref.mean)*intensity + BAS, family = gaussian, full.df)
summary(sugar.linmod.4)

Call:
glm(formula = sugar.intake ~ log10(pref.mean) * intensity + BAS, 
    family = gaussian, data = full.df)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-368.46  -162.05   -26.14   138.51   637.36  

Coefficients:
                           Estimate Std. Error t value Pr(>|t|)    
(Intercept)                1338.497    232.349   5.761 3.07e-07 ***
log10(pref.mean)            557.639    499.232   1.117 0.268451    
intensity                    -2.582      3.028  -0.852 0.397363    
BAS                         -14.065      3.401  -4.136 0.000112 ***
log10(pref.mean):intensity   -2.540      6.274  -0.405 0.687088    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for gaussian family taken to be 47389.16)

    Null deviance: 4890850  on 64  degrees of freedom
Residual deviance: 2843350  on 60  degrees of freedom
AIC: 891.06

Number of Fisher Scoring iterations: 2

## Compare models
anova(sugar.linmod.2, sugar.linmod.1, test = "Chisq")

Analysis of Deviance Table

Model 1: sugar.intake ~ log10(pref.mean) + intensity
Model 2: sugar.intake ~ log10(pref.mean)
  Resid. Df Resid. Dev Df Deviance Pr(>Chi)
1        62    3663679                     
2        63    3789415 -1  -125736   0.1446

anova(sugar.linmod.2, sugar.linmod.3, test = "Chisq")

Analysis of Deviance Table

Model 1: sugar.intake ~ log10(pref.mean) + intensity
Model 2: sugar.intake ~ log10(pref.mean) + intensity + BAS
  Resid. Df Resid. Dev Df Deviance  Pr(>Chi)    
1        62    3663679                          
2        61    2851114  1   812565 3.052e-05 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

anova(sugar.linmod.2, sugar.linmod.4, test = "Chisq")

Analysis of Deviance Table

Model 1: sugar.intake ~ log10(pref.mean) + intensity
Model 2: sugar.intake ~ log10(pref.mean) * intensity + BAS
  Resid. Df Resid. Dev Df Deviance  Pr(>Chi)    
1        62    3663679                          
2        60    2843350  2   820329 0.0001742 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Logistic regression - sHEI

## The effect of sweet preference on sugar intake
full.df %>%
  ggplot(aes(log10(pref.mean), sHEI)) +
  geom_point() +
  geom_smooth(method = lm) +
  stat_cor() +
  theme_bw()

`geom_smooth()` using formula = 'y ~ x'

sHEI.logmod.1 <- multinom(as.factor(DQ) ~ pref.mean, family = binomial, data = full.df)

# weights:  3 (2 variable)
initial  value 45.054567 
final  value 44.853743 
converged

tidy(multinom(as.factor(DQ) ~ pref.mean, family = binomial, data = full.df), exponentiate = TRUE, conf.int = TRUE)

# weights:  3 (2 variable)
initial  value 45.054567 
final  value 44.853743 
converged

# A tibble: 2 × 8
  y.level term        estimate std.error statistic p.value conf.low conf.high
  <chr>   <chr>          <dbl>     <dbl>     <dbl>   <dbl>    <dbl>     <dbl>
1 1       (Intercept)    1.24      0.471     0.462   0.644   0.494       3.13
2 1       pref.mean      0.558     0.943    -0.619   0.536   0.0879      3.54

ggcoef_multinom(sHEI.logmod.1)

! `broom::tidy()` failed to tidy the model.
✔ `tidy_parameters()` used instead.
ℹ Add `tidy_fun = broom.helpers::tidy_parameters` to quiet these messages.

Linear regression - sHEI

## The effect of sweet preference on sHEI
summary(glm(sHEI ~ log10(pref.mean), family = gaussian, full.df))

Call:
glm(formula = sHEI ~ log10(pref.mean), family = gaussian, data = full.df)

Deviance Residuals: 
     Min        1Q    Median        3Q       Max  
-19.5521   -6.0489   -0.8705    6.8385   16.4449  

Coefficients:
                 Estimate Std. Error t value Pr(>|t|)    
(Intercept)        41.847      1.951  21.448   <2e-16 ***
log10(pref.mean)   -5.888      3.478  -1.693   0.0954 .  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for gaussian family taken to be 74.34753)

    Null deviance: 4897.0  on 64  degrees of freedom
Residual deviance: 4683.9  on 63  degrees of freedom
AIC: 468.5

Number of Fisher Scoring iterations: 2