Body Appreciation, Eating Behaviors, and Food Liking on Diet Quality in Adult Women (< 40 years old), Qualtric Survey Results

*** Caren thesis hypotheses

Aim 1: Determine the influence of sweet foods and beverages liking on dietary patterns in young and early adult women Hypothesis 1: The degree of liking for sweet foods and beverages will be negatively associated with added sugars intake and the healthy eating index

Aim 2: Determine the influence of self-body appreciation on dietary patterns in young and early adult women Hypothesis 2: Self-body appreciation will be positively associated with added sugars intake and the healthy eating index

Aim 3: Determine the influence of eating behaviors on dietary patterns in young and early adult women Hypothesis 3.a: Uncontrolled eating score will be negatively associated with added sugars intake and the healthy eating index Hypothesis 3.b: Emotional eating score will be negatively associated with added sugars intake and the healthy eating index score
Hypothesis 3.c: Cognitive restraint score will be positively associated with added sugars intake and the healthy eating index

Aim 4: Determine the interactions across sweet foods and beverages liking, self-body appreciation, and eating behaviors on dietary patterns in young and early adult women Hypothesis 4.a: Self-body appreciation and eating behaviors (uncontrolled eating, emotional eating, and cognitive restraint) will interact with the relationship between sweet foods and beverages liking and added sugars intake Hypothesis 4.b: Self-body appreciation and eating behaviors (uncontrolled eating, emotional eating, and cognitive restraint) will interact with the relationship between sweet foods and beverages liking and the healthy eating index

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

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

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

## survey df
survey.df <- read_excel("EatBehav_23OCT24.xlsx")

survey.df <- survey.df %>% 
  filter(!row_number() %in% 1) %>%
  mutate(across(-c("DM-05", "DM-05_7_TEXT", "DM-06", "DM-08.1"), as.numeric)) %>%
    filter(Duration > 900) %>%
  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", pob = "DM-08.1", 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, age <= 40)

Functions

Useful functions used in this script

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

## Correlation function
corrfunc <- function(data, Cor1, Cor2, title) {
  ggplot(data = data, aes(x = {{Cor1}}, y = {{Cor2}})) +
  geom_point()+
  geom_smooth(method = lm)+
  ggtitle(title) +
  stat_cor(method = "pearson") +
  theme_light() 
}

## Interaction visualization function
intfunc <- function(data, Cor1, Cor2, Group, title) {
  ggplot(data = data, aes(x = {{Cor1}}, y = {{Cor2}}, color = {{Group}})) +
  geom_point()+
  geom_smooth(method = lm, se = FALSE)+
  ggtitle(title) +
  stat_cor(method = "pearson") +
  theme_light() 
}

Cronbach’s alpha for food liking groups

Agreement across food items within the groups from the Food Liking Survey

## Vegetable 
veg.cor <- survey.df %>%
  select("FL-01_3":"FL-06_3") %>%
  mutate(overall = rowMeans(select(., "FL-01_3":"FL-06_3"), na.rm = TRUE)) %>%
  cor(., use = "complete.obs")

survey.df %>% 
  select("FL-01_3":"FL-06_3")%>%
  na.omit() %>%
  cronbach.alpha(., standardized = TRUE)

Standardized Cronbach's alpha for the '.' data-set

Items: 6
Sample units: 256
alpha: 0.694

corrplot(veg.cor, method = "number", title = "Corrplot: liking for vegetables")

## Fruit 
fruit.cor <- survey.df %>%
  select("FL-07_3":"FL-11_3") %>%
  mutate(overall = rowMeans(select(., "FL-07_3":"FL-11_3"), na.rm = TRUE)) %>%
  cor(., use = "complete.obs")

survey.df %>% 
  select("FL-07_3":"FL-11_3")%>%
  na.omit() %>%
  cronbach.alpha(., standardized = TRUE)

Standardized Cronbach's alpha for the '.' data-set

Items: 5
Sample units: 267
alpha: 0.552

corrplot(fruit.cor, method = "number", title = "Corrplot: liking for fruits")

## Salty/fat 
saltyfat.cor <- survey.df %>%
  select("FL-12_3":"FL-15_3") %>%
  mutate(overall = rowMeans(select(., "FL-12_3":"FL-15_3"), na.rm = TRUE)) %>%
  cor(., use = "complete.obs")

survey.df %>% 
  select("FL-12_3":"FL-15_3")%>%
  na.omit() %>%
  cronbach.alpha(., standardized = TRUE)

Standardized Cronbach's alpha for the '.' data-set

Items: 4
Sample units: 267
alpha: 0.487

corrplot(saltyfat.cor, method = "number", title = "Corrplot: liking for salty & fat foods")

## High-fat proteins
saltyfat.cor <- survey.df %>%
  select("FL-16_3":"FL-20_3") %>%
  mutate(overall = rowMeans(select(., "FL-16_3":"FL-20_3"), na.rm = TRUE)) %>%
  cor(., use = "complete.obs")

survey.df %>% 
  select("FL-16_3":"FL-20_3")%>%
  na.omit() %>%
  cronbach.alpha(., standardized = TRUE)

Standardized Cronbach's alpha for the '.' data-set

Items: 5
Sample units: 242
alpha: 0.625

corrplot(saltyfat.cor, method = "number", title = "Corrplot: liking for salty & fat foods")

## Alcohol 
alcohol.cor <- survey.df %>%
  select("FL-21_3":"FL-25_3") %>%
  mutate(overall = rowMeans(select(., "FL-21_3":"FL-25_3"), na.rm = TRUE)) %>%
  cor(., use = "complete.obs")

survey.df %>% 
  select("FL-21_3":"FL-25_3")%>%
  na.omit() %>%
  cronbach.alpha(., standardized = TRUE)

Standardized Cronbach's alpha for the '.' data-set

Items: 5
Sample units: 185
alpha: 0.755

corrplot(alcohol.cor, method = "number", title = "Corrplot: liking for alcohol")

## Whole grain 
wg.cor <- survey.df %>%
  select("FL-26_3":"FL-29_3") %>%
  mutate(overall = rowMeans(select(., "FL-26_3":"FL-29_3"), na.rm = TRUE)) %>%
  cor(., use = "complete.obs")

survey.df %>% 
  select("FL-26_3":"FL-29_3")%>%
  na.omit() %>%
  cronbach.alpha(., standardized = TRUE)

Standardized Cronbach's alpha for the '.' data-set

Items: 4
Sample units: 266
alpha: 0.751

corrplot(wg.cor, method = "number", title = "Corrplot: liking for whole grains")

## Carbs
carbs.cor <- survey.df %>%
  select("FL-30_3":"FL-34_3") %>%
  mutate(overall = rowMeans(select(., "FL-30_3":"FL-34_3"), na.rm = TRUE)) %>%
  cor(., use = "complete.obs")

survey.df %>% 
  select("FL-30_3":"FL-34_3")%>%
  na.omit() %>%
  cronbach.alpha(., standardized = TRUE)

Standardized Cronbach's alpha for the '.' data-set

Items: 5
Sample units: 264
alpha: 0.574

corrplot(carbs.cor, method = "number", title = "Corrplot: liking for carbohydrates")

## Healthy fats
healthyfat.cor <- survey.df %>%
  select("FL-35_3":"FL-38_3") %>%
  mutate(overall = rowMeans(select(., "FL-35_3":"FL-38_3"), na.rm = TRUE)) %>%
  cor(., use = "complete.obs")

survey.df %>% 
  select("FL-35_3":"FL-38_3")%>%
  na.omit() %>%
  cronbach.alpha(., standardized = TRUE)

Standardized Cronbach's alpha for the '.' data-set

Items: 4
Sample units: 247
alpha: 0.591

corrplot(healthyfat.cor, method = "number", title = "Corrplot: liking for healthy fats")

## Unhealthy fats
unhealthyfat.cor <- survey.df %>%
  select("FL-48_3":"FL-52_3") %>%
  mutate(overall = rowMeans(select(., "FL-48_3":"FL-52_3"), na.rm = TRUE)) %>%
  cor(., use = "complete.obs")

survey.df %>% 
  select("FL-48_3":"FL-52_3")%>%
  na.omit() %>%
  cronbach.alpha(., standardized = TRUE)

Standardized Cronbach's alpha for the '.' data-set

Items: 5
Sample units: 265
alpha: 0.547

corrplot(unhealthyfat.cor, method = "number", title = "Corrplot: liking for unhealthy fats")

## All fat
allfat.cor <- survey.df %>%
  select("FL-35_3":"FL-38_3", "FL-48_3":"FL-52_3") %>%
  mutate(overall = rowMeans(select(.,"FL-35_3":"FL-38_3", "FL-48_3":"FL-52_3"), na.rm = TRUE)) %>%
  cor(., use = "complete.obs")

survey.df %>% 
  select("FL-35_3":"FL-38_3", "FL-48_3":"FL-52_3")%>%
  na.omit() %>%
  cronbach.alpha(., standardized = TRUE)

Standardized Cronbach's alpha for the '.' data-set

Items: 9
Sample units: 244
alpha: 0.603

corrplot(allfat.cor, method = "number", title = "Corrplot: liking for all fats")

## Sweets
sweets.cor <- survey.df %>%
  select("FL-39_3":"FL-43_3") %>%
  mutate(overall = rowMeans(select(., "FL-39_3":"FL-43_3"), na.rm = TRUE)) %>%
  cor(., use = "complete.obs")

survey.df %>% 
  select("FL-39_3":"FL-43_3")%>%
  na.omit() %>%
  cronbach.alpha(., standardized = TRUE)

Standardized Cronbach's alpha for the '.' data-set

Items: 5
Sample units: 270
alpha: 0.767

corrplot(sweets.cor, method = "number", title = "Corrplot: liking for sweets")

## Sugar-sweetened beverages 
ssb.cor <- survey.df %>%
  select("FL-44_3":"FL-47_3") %>%
  mutate(ssb.overall = rowMeans(select(., "FL-44_3":"FL-47_3"), na.rm = TRUE)) %>% 
  cor(., use = "complete.obs")

survey.df %>% 
  select("FL-44_3":"FL-47_3")%>%
  na.omit() %>%
  cronbach.alpha(., standardized = TRUE)

Standardized Cronbach's alpha for the '.' data-set

Items: 4
Sample units: 233
alpha: 0.621

corrplot(ssb.cor, method = "number", title = "Corrplot: liking for SSB")

## Sweet and fatty foods 
sufat.cor <- survey.df %>%
  select("FL-39_3":"FL-40_3", "FL-43_3", "FL-45_3") %>%
  mutate(sufat.overall = rowMeans(select(., "FL-39_3":"FL-40_3", "FL-43_3", "FL-45_3"), na.rm = TRUE)) %>% 
  cor(., use = "complete.obs")

survey.df %>% 
  select("FL-39_3":"FL-40_3", "FL-43_3", "FL-45_3")%>%
  na.omit() %>%
  cronbach.alpha(., standardized = TRUE)

Standardized Cronbach's alpha for the '.' data-set

Items: 4
Sample units: 265
alpha: 0.659

corrplot(sufat.cor, method = "number", title = "Corrplot: sweet + fatty foods")

## Sweet foods and sugar-sweetened beverages 
sfbl.cor <- survey.df %>%
  select("FL-39_3":"FL-47_3") %>%
  mutate(ssb.overall = rowMeans(select(., "FL-39_3":"FL-47_3"), na.rm = TRUE)) %>% 
  cor(., use = "complete.obs")

survey.df %>% 
  select("FL-39_3":"FL-47_3")%>%
  na.omit() %>%
  cronbach.alpha(., standardized = TRUE)

Standardized Cronbach's alpha for the '.' data-set

Items: 9
Sample units: 233
alpha: 0.78

corrplot(sfbl.cor, method = "number", title = "Corrplot: liking for SSB")

## Protein
protein.cor <- survey.df %>%
  select("FL-53_3":"FL-57_3") %>%
  mutate(overall = rowMeans(select(., "FL-53_3":"FL-57_3"), na.rm = TRUE)) %>%
  cor(., use = "complete.obs") ## soy protein does not fit

survey.df %>% 
  select("FL-53_3":"FL-57_3")%>%
  na.omit() %>%
  cronbach.alpha(., standardized = TRUE)

Standardized Cronbach's alpha for the '.' data-set

Items: 5
Sample units: 203
alpha: 0.423

corrplot(protein.cor, method = "number", title = "Corrplot: liking for protein")

## Spicy
spicy.cor <- survey.df %>%
  select("FL-58_3":"FL-59_3") %>%
  mutate(overall = rowMeans(select(., "FL-58_3":"FL-59_3"), na.rm = TRUE)) %>%
  cor(., use = "complete.obs")

survey.df %>% 
  select("FL-58_3":"FL-59_3")%>%
  na.omit() %>%
  cronbach.alpha(., standardized = TRUE)

Standardized Cronbach's alpha for the '.' data-set

Items: 2
Sample units: 258
alpha: 0.875

corrplot(spicy.cor, method = "number", title = "Corrplot: liking for spicy foods")

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(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")) 

## 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(
    `WG.day` == 1 ~ 0.51,
    sex == 1 & `WG.day` >= 2 ~ 2.97,
    sex == 2 & `WG.day` >= 2 & `WG.day` <= 3 ~ 5.20,
    sex == 2 & `WG.day` >= 4 ~ 6.94)) %>%
  mutate(dairyHEI = case_when(
    sex == 1 & `Milk.day` <= 3 ~ 3.22,
    sex == 2 & `Milk.day` <= 3 & `LFmilk.day` == 1 ~ 3.32,
    sex == 2 & `Milk.day` <= 3 & `LFmilk.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.day` >= 1 & `LFmilk.day` <= 2 ~ 2.63,
    `Saturated Fats` >= 2 & `Saturated Fats` <= 3 & is.na(`Milk.day`) | is.na(`LFmilk.day`) ~ 2.63,
    `Saturated Fats` >= 2 & `Saturated Fats` <= 3 & `Milk.day` >= 1 & `LFmilk.day` >= 3 ~ 4.54,
    `Saturated Fats` == 1 & `Milk.day` >= 1 ~ 5.93,
    `Saturated Fats` == 1 & is.na(`Milk.day`) ~ 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(
    `SS beverages` == 1  ~ 0,
    `SS beverages` == 2 ~ 156,
    `SS beverages` == 3 ~ 312,
    `SS beverages` == 4 ~ 468,
    `SS beverages` == 5 ~ 624,
    `SS beverages` == 6 ~ 780,
    `SS beverages` == 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 frequency` <= 4 ~ 13.26,
    `SSB frequency` >= 5 & `SSB frequency` <=6 ~ 16.00,
    `SS beverages` ==2 ~ 16.00,
    `SS beverages` >=3 ~ 26.87)) %>%
  mutate(satfatHEI = case_when(
    `SS beverages` >= 3 ~ 1.82,
    `SS beverages` <= 2 & `Grains.day` <= 2 ~ 3.20,
    `SS beverages` <= 2 & `Grains.day` >= 3 & Nuts <= 2 ~ 4.64,
    `SS beverages` <= 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, 3))

## 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  %>%
  mutate(BMI = (weight/(height^2)*703)) %>%
  select(c("ID", "sex", "age", "BMI", "race", "usb", "education", "hinc", "exec")) %>%   
  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, 3)) %>%
  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)) 

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

Data distribution

## 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.90471, p-value = 4.763e-12

## 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.82501, p-value < 2.2e-16

## 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.97872, p-value = 0.0004628

## 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.9851, p-value = 0.006561

## 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.98244, p-value = 0.002094

## 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.94537, p-value = 1.752e-08

## 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.9809, p-value = 0.001106

## Sugar intake
histfunc(full.df, sugar.intake , "Added Sugars intake index histogram")

qqfunc(full.df, sugar.intake, "Added Sugars intake QQ plot") 

shapiro.test(full.df$sugar.intake) 

    Shapiro-Wilk normality test

data:  full.df$sugar.intake
W = 0.89448, p-value = 8.846e-13

## sHEI eating
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.99193, p-value = 0.1472

Demographics and counts

## Summary table - by migration
demo.DQ.df <- full.df %>% 
  select(-ID) %>%
  mutate(n = n()) %>%
  group_by(DQ) %>%
  get_summary_stats()

write_xlsx(demo.DQ.df,"/Users/maycheung/Documents/Analyses/EatBehav//Demographics.DQ.xlsx")

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

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

## Counts - by ancestry and age
full.df %>%
  group_by(race) %>%
  count()

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

## Graph
      full.df %>%
        group_by(race) %>%
        count() %>%
        ggplot(aes(x = race, y = n, fill = as.factor(race))) +
        geom_bar(stat = "identity", position = position_dodge()) +
        labs(x = "Ancestry", y = "n") +
        theme_bw() +
        theme(strip.text.x = element_text(size = 12, face = "bold")) 

Correlations - overall corrplot

## Corr DF
corr.df <- full.df %>%
  select(age, BMI, BAS, hinc, uncontrolled, restraint, emotional, sfbl.liking, 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"))

Linear regression - added sugars intake

## Models
sugar.lm.M1 <- lm(sugar.intake ~ sfbl.liking, data = full.df)
summary(sugar.lm.M1) # r2 = 0.188

Call:
lm(formula = sugar.intake ~ sfbl.liking, data = full.df)

Residuals:
    Min      1Q  Median      3Q     Max 
-495.46 -153.76  -41.99  125.33  854.99 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 331.2003    22.9435  14.435  < 2e-16 ***
sfbl.liking   4.0005     0.5031   7.951 5.18e-14 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 232.4 on 268 degrees of freedom
Multiple R-squared:  0.1909,    Adjusted R-squared:  0.1879 
F-statistic: 63.22 on 1 and 268 DF,  p-value: 5.182e-14

sugar.lm.M2 <- lm(sugar.intake ~ sfbl.liking*BAS , data = full.df)
summary(sugar.lm.M2) # r2 = 0.197

Call:
lm(formula = sugar.intake ~ sfbl.liking * BAS, data = full.df)

Residuals:
    Min      1Q  Median      3Q     Max 
-483.35 -156.20  -44.45  130.27  803.57 

Coefficients:
                 Estimate Std. Error t value Pr(>|t|)    
(Intercept)     419.97238  106.64470   3.938 0.000105 ***
sfbl.liking       5.23913    2.32345   2.255 0.024953 *  
BAS              -2.46317    2.77648  -0.887 0.375796    
sfbl.liking:BAS  -0.03518    0.06043  -0.582 0.560939    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 231.2 on 266 degrees of freedom
Multiple R-squared:  0.2055,    Adjusted R-squared:  0.1966 
F-statistic: 22.94 on 3 and 266 DF,  p-value: 3.07e-13

sugar.lm.M3 <- lm(sugar.intake ~ sfbl.liking*uncontrolled , data = full.df)
summary(sugar.lm.M3) # r2 = 0.206

Call:
lm(formula = sugar.intake ~ sfbl.liking * uncontrolled, data = full.df)

Residuals:
    Min      1Q  Median      3Q     Max 
-492.60 -147.89  -38.82  124.69  793.37 

Coefficients:
                          Estimate Std. Error t value Pr(>|t|)    
(Intercept)              276.27047   82.25001   3.359 0.000897 ***
sfbl.liking                1.42541    1.86735   0.763 0.445942    
uncontrolled               3.01750    4.38652   0.688 0.492113    
sfbl.liking:uncontrolled   0.12801    0.09658   1.325 0.186172    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 229.9 on 266 degrees of freedom
Multiple R-squared:  0.2146,    Adjusted R-squared:  0.2057 
F-statistic: 24.22 on 3 and 266 DF,  p-value: 6.865e-14

sugar.lm.M4 <- lm(sugar.intake ~ sfbl.liking*BAS + sfbl.liking*uncontrolled, data = full.df)
summary(sugar.lm.M4) # r2 = 0.207

Call:
lm(formula = sugar.intake ~ sfbl.liking * BAS + sfbl.liking * 
    uncontrolled, data = full.df)

Residuals:
    Min      1Q  Median      3Q     Max 
-465.83 -147.26  -39.03  124.96  807.68 

Coefficients:
                           Estimate Std. Error t value Pr(>|t|)  
(Intercept)              385.649363 154.968803   2.489   0.0134 *
sfbl.liking                1.771675   3.562022   0.497   0.6193  
BAS                       -2.497744   2.881033  -0.867   0.3868  
uncontrolled               2.023059   4.569931   0.443   0.6584  
sfbl.liking:BAS           -0.006712   0.063750  -0.105   0.9162  
sfbl.liking:uncontrolled   0.121411   0.102459   1.185   0.2371  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 229.7 on 264 degrees of freedom
Multiple R-squared:  0.2215,    Adjusted R-squared:  0.2067 
F-statistic: 15.02 on 5 and 264 DF,  p-value: 5.58e-13

sugar.lm.M5 <- lm(sugar.intake ~ sfbl.liking*BAS + sfbl.liking*uncontrolled + age + hinc, data = full.df)
summary(sugar.lm.M5) # r2 = 0.211

Call:
lm(formula = sugar.intake ~ sfbl.liking * BAS + sfbl.liking * 
    uncontrolled + age + hinc, data = full.df)

Residuals:
    Min      1Q  Median      3Q     Max 
-465.12 -146.57  -44.76  104.52  830.89 

Coefficients:
                           Estimate Std. Error t value Pr(>|t|)   
(Intercept)              484.640759 174.983512   2.770  0.00601 **
sfbl.liking                0.944575   3.580184   0.264  0.79212   
BAS                       -2.669014   2.878043  -0.927  0.35459   
uncontrolled               0.553664   4.628055   0.120  0.90487   
age                        0.045718   3.082230   0.015  0.98818   
hinc                     -11.185241   6.283200  -1.780  0.07620 . 
sfbl.liking:BAS           -0.001727   0.063699  -0.027  0.97839   
sfbl.liking:uncontrolled   0.146958   0.103103   1.425  0.15525   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 229.1 on 262 degrees of freedom
Multiple R-squared:  0.2316,    Adjusted R-squared:  0.211 
F-statistic: 11.28 on 7 and 262 DF,  p-value: 1.734e-12

## Model comparison
anova(sugar.lm.M5, sugar.lm.M1) ## significant

Analysis of Variance Table

Model 1: sugar.intake ~ sfbl.liking * BAS + sfbl.liking * uncontrolled + 
    age + hinc
Model 2: sugar.intake ~ sfbl.liking
  Res.Df      RSS Df Sum of Sq      F  Pr(>F)  
1    262 13749134                              
2    268 14477308 -6   -728174 2.3126 0.03414 *
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

anova(sugar.lm.M5, sugar.lm.M2)

Analysis of Variance Table

Model 1: sugar.intake ~ sfbl.liking * BAS + sfbl.liking * uncontrolled + 
    age + hinc
Model 2: sugar.intake ~ sfbl.liking * BAS
  Res.Df      RSS Df Sum of Sq      F  Pr(>F)  
1    262 13749134                              
2    266 14214795 -4   -465661 2.2184 0.06739 .
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

anova(sugar.lm.M5, sugar.lm.M3)

Analysis of Variance Table

Model 1: sugar.intake ~ sfbl.liking * BAS + sfbl.liking * uncontrolled + 
    age + hinc
Model 2: sugar.intake ~ sfbl.liking * uncontrolled
  Res.Df      RSS Df Sum of Sq      F Pr(>F)
1    262 13749134                           
2    266 14053410 -4   -304276 1.4496 0.2181

anova(sugar.lm.M5, sugar.lm.M4)

Analysis of Variance Table

Model 1: sugar.intake ~ sfbl.liking * BAS + sfbl.liking * uncontrolled + 
    age + hinc
Model 2: sugar.intake ~ sfbl.liking * BAS + sfbl.liking * uncontrolled
  Res.Df      RSS Df Sum of Sq      F Pr(>F)
1    262 13749134                           
2    264 13929881 -2   -180747 1.7221 0.1807

Logistic regression - added sugars intake

## Models
sugar.logit.M1 <- multinom(SugarCat ~ sfbl.liking, data = full.df)

# weights:  9 (4 variable)
initial  value 296.625318 
iter  10 value 264.991749
iter  10 value 264.991749
final  value 264.991749 
converged

tidy(sugar.logit.M1, 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.482   0.242       -3.02 2.53e- 3   0.300      0.774
2 2       sfbl.liking    1.03    0.00624      4.06 4.99e- 5   1.01       1.04 
3 3       (Intercept)    0.154   0.339       -5.52 3.38e- 8   0.0790     0.299
4 3       sfbl.liking    1.05    0.00759      6.74 1.63e-11   1.04       1.07 

summary(sugar.logit.M1) # AIC = 538.0

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

Coefficients:
  (Intercept) sfbl.liking
2  -0.7294329  0.02532107
3  -1.8730620  0.05113017

Std. Errors:
  (Intercept) sfbl.liking
2   0.2415343 0.006242891
3   0.3392963 0.007590864

Residual Deviance: 529.9835 
AIC: 537.9835 

sugar.logit.M2 <- multinom(SugarCat ~ sfbl.liking*BAS , data = full.df)

# weights:  15 (8 variable)
initial  value 296.625318 
iter  10 value 260.301344
final  value 260.291714 
converged

tidy(sugar.logit.M2, exponentiate = TRUE, conf.int = TRUE)

# A tibble: 8 × 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)      5.50   0.000424  4027.    0           5.50      5.51 
2 2       sfbl.liking      1.00   0.0199       0.201 8.40e- 1    0.966     1.04 
3 2       BAS              0.936  0.00667     -9.93  3.18e-23    0.924     0.948
4 2       sfbl.liking:…    1.00   0.000546     1.03  3.01e- 1    0.999     1.00 
5 3       (Intercept)      1.77   0.000438  1301.    0           1.77      1.77 
6 3       sfbl.liking      1.04   0.0194       2.02  4.38e- 2    1.00      1.08 
7 3       BAS              0.935  0.00941     -7.13  9.71e-13    0.918     0.952
8 3       sfbl.liking:…    1.00   0.000539     0.589 5.56e- 1    0.999     1.00 

summary(sugar.logit.M2) # AIC = 536.58

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

Coefficients:
  (Intercept) sfbl.liking         BAS sfbl.liking:BAS
2   1.7055250 0.004007762 -0.06616344    0.0005643821
3   0.5701559 0.039023569 -0.06710934    0.0003170021

Std. Errors:
   (Intercept) sfbl.liking         BAS sfbl.liking:BAS
2 0.0004235405  0.01990180 0.006665093    0.0005460977
3 0.0004383018  0.01936088 0.009406249    0.0005386382

Residual Deviance: 520.5834 
AIC: 536.5834 

sugar.logit.M3 <- multinom(SugarCat ~ sfbl.liking*uncontrolled , data = full.df)

# weights:  15 (8 variable)
initial  value 296.625318 
iter  10 value 262.078374
final  value 261.967235 
converged

tidy(sugar.logit.M3, exponentiate = TRUE, conf.int = TRUE)

# A tibble: 8 × 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.216   0.824       -1.86  6.31e- 2   0.0429    1.09  
2 2       sfbl.liking     1.03    0.0234       1.32  1.86e- 1   0.985     1.08  
3 2       uncontrolled    1.05    0.0436       1.02  3.09e- 1   0.960     1.14  
4 2       sfbl.liking:…   1.00    0.00122     -0.269 7.88e- 1   0.997     1.00  
5 3       (Intercept)     0.0277  0.266      -13.5   2.74e-41   0.0165    0.0467
6 3       sfbl.liking     1.06    0.0179       3.00  2.66e- 3   1.02      1.09  
7 3       uncontrolled    1.10    0.0206       4.47  7.94e- 6   1.05      1.14  
8 3       sfbl.liking:…   1.00    0.000965    -0.208 8.35e- 1   0.998     1.00  

summary(sugar.logit.M3) # AIC = 523.93

Call:
multinom(formula = SugarCat ~ sfbl.liking * uncontrolled, data = full.df)

Coefficients:
  (Intercept) sfbl.liking uncontrolled sfbl.liking:uncontrolled
2   -1.532351  0.03094270   0.04432119            -0.0003278974
3   -3.585215  0.05376868   0.09191967            -0.0002006819

Std. Errors:
  (Intercept) sfbl.liking uncontrolled sfbl.liking:uncontrolled
2   0.8243882  0.02337448   0.04359835             0.0012205889
3   0.2663901  0.01789514   0.02057903             0.0009651083

Residual Deviance: 523.9345 
AIC: 539.9345 

sugar.logit.M4 <- multinom(SugarCat ~ sfbl.liking*BAS + sfbl.liking*uncontrolled, data = full.df)

# weights:  21 (12 variable)
initial  value 296.625318 
iter  10 value 270.545821
iter  20 value 258.327530
iter  20 value 258.327530
iter  20 value 258.327530
final  value 258.327530 
converged

tidy(sugar.logit.M4, exponentiate = TRUE, conf.int = TRUE)

# A tibble: 12 × 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)      2.83   0.000861  1210.    0          2.83      2.84 
 2 2       sfbl.liking      1.01   0.0306       0.376 0.707      0.953     1.07 
 3 2       BAS              0.938  0.0198      -3.25  0.00114    0.902     0.975
 4 2       uncontrolled     1.03   0.0380       0.857 0.392      0.959     1.11 
 5 2       sfbl.liking:…    1.00   0.000671     0.830 0.407      0.999     1.00 
 6 2       sfbl.liking:…    1.00   0.00121     -0.322 0.748      0.997     1.00 
 7 3       (Intercept)      0.296  0.00100  -1214.    0          0.295     0.296
 8 3       sfbl.liking      1.04   0.0300       1.43  0.151      0.984     1.11 
 9 3       BAS              0.941  0.0284      -2.12  0.0336     0.890     0.995
10 3       uncontrolled     1.09   0.0523       1.59  0.113      0.981     1.20 
11 3       sfbl.liking:…    1.00   0.000753     0.503 0.615      0.999     1.00 
12 3       sfbl.liking:…    1.00   0.00138     -0.275 0.784      0.997     1.00 

summary(sugar.logit.M4) # AIC = 516.66

Call:
multinom(formula = SugarCat ~ sfbl.liking * BAS + sfbl.liking * 
    uncontrolled, data = full.df)

Coefficients:
  (Intercept) sfbl.liking         BAS uncontrolled sfbl.liking:BAS
2    1.041807  0.01151822 -0.06441351   0.03258162    0.0005568012
3   -1.217735  0.04303711 -0.06036361   0.08306175    0.0003790823
  sfbl.liking:uncontrolled
2            -0.0003904770
3            -0.0003785954

Std. Errors:
   (Intercept) sfbl.liking        BAS uncontrolled sfbl.liking:BAS
2 0.0008608408  0.03064461 0.01979883   0.03803042    0.0006710809
3 0.0010030760  0.02999624 0.02841231   0.05234780    0.0007533132
  sfbl.liking:uncontrolled
2              0.001213571
3              0.001378089

Residual Deviance: 516.6551 
AIC: 540.6551 

sugar.logit.M5 <- multinom(SugarCat ~ sfbl.liking*BAS + sfbl.liking*uncontrolled + age + hinc, data = full.df)

# weights:  27 (16 variable)
initial  value 296.625318 
iter  10 value 262.170219
iter  20 value 254.681716
final  value 254.671484 
converged

tidy(sugar.logit.M5, exponentiate = TRUE, conf.int = TRUE)

# A tibble: 16 × 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)      4.10   0.000832 1696.     0          4.09      4.10 
 2 2       sfbl.liking      1.01   0.0336      0.336  0.737      0.947     1.08 
 3 2       BAS              0.938  0.0233     -2.76   0.00581    0.896     0.982
 4 2       uncontrolled     1.03   0.0408      0.755  0.450      0.952     1.12 
 5 2       age              0.990  0.0326     -0.317  0.751      0.928     1.06 
 6 2       hinc             0.986  0.0695     -0.205  0.838      0.860     1.13 
 7 2       sfbl.liking:…    1.00   0.000719    0.742  0.458      0.999     1.00 
 8 2       sfbl.liking:…    1.00   0.00126    -0.288  0.773      0.997     1.00 
 9 3       (Intercept)      2.16   0.00109   707.     0          2.16      2.17 
10 3       sfbl.liking      1.03   0.0345      0.881  0.378      0.963     1.10 
11 3       BAS              0.936  0.0314     -2.09   0.0366     0.881     0.996
12 3       uncontrolled     1.06   0.0556      1.06   0.288      0.951     1.18 
13 3       age              0.984  0.0380     -0.423  0.672      0.913     1.06 
14 3       hinc             0.834  0.0808     -2.25   0.0246     0.712     0.977
15 3       sfbl.liking:…    1.00   0.000804    0.604  0.546      0.999     1.00 
16 3       sfbl.liking:…    1.00   0.00144     0.0224 0.982      0.997     1.00 

summary(sugar.logit.M5) # AIC = 509.34

Call:
multinom(formula = SugarCat ~ sfbl.liking * BAS + sfbl.liking * 
    uncontrolled + age + hinc, data = full.df)

Coefficients:
  (Intercept) sfbl.liking         BAS uncontrolled         age        hinc
2   1.4103657  0.01131023 -0.06435534   0.03081459 -0.01034661 -0.01420854
3   0.7722927  0.03040546 -0.06566710   0.05906263 -0.01606100 -0.18159067
  sfbl.liking:BAS sfbl.liking:uncontrolled
2    0.0005333484            -3.624714e-04
3    0.0004856087             3.236772e-05

Std. Errors:
   (Intercept) sfbl.liking        BAS uncontrolled        age       hinc
2 0.0008317089  0.03364017 0.02333354   0.04081206 0.03263007 0.06945721
3 0.0010917661  0.03450672 0.03141103   0.05557827 0.03797068 0.08080109
  sfbl.liking:BAS sfbl.liking:uncontrolled
2    0.0007188063              0.001259064
3    0.0008036372              0.001444313

Residual Deviance: 509.343 
AIC: 541.343 

## Model comparison
anova(sugar.logit.M5, sugar.logit.M1)

Likelihood ratio tests of Multinomial Models

Response: SugarCat
                                                        Model Resid. df
1                                                 sfbl.liking       536
2 sfbl.liking * BAS + sfbl.liking * uncontrolled + age + hinc       524
  Resid. Dev   Test    Df LR stat.    Pr(Chi)
1   529.9835                                 
2   509.3430 1 vs 2    12 20.64053 0.05589867

anova(sugar.logit.M5, sugar.logit.M2)

Likelihood ratio tests of Multinomial Models

Response: SugarCat
                                                        Model Resid. df
1                                           sfbl.liking * BAS       532
2 sfbl.liking * BAS + sfbl.liking * uncontrolled + age + hinc       524
  Resid. Dev   Test    Df LR stat.  Pr(Chi)
1   520.5834                               
2   509.3430 1 vs 2     8 11.24046 0.188443

anova(sugar.logit.M5, sugar.logit.M3)

Likelihood ratio tests of Multinomial Models

Response: SugarCat
                                                        Model Resid. df
1                                  sfbl.liking * uncontrolled       532
2 sfbl.liking * BAS + sfbl.liking * uncontrolled + age + hinc       524
  Resid. Dev   Test    Df LR stat.    Pr(Chi)
1   523.9345                                 
2   509.3430 1 vs 2     8  14.5915 0.06759237

anova(sugar.logit.M5, sugar.logit.M4)

Likelihood ratio tests of Multinomial Models

Response: SugarCat
                                                        Model Resid. df
1              sfbl.liking * BAS + sfbl.liking * uncontrolled       528
2 sfbl.liking * BAS + sfbl.liking * uncontrolled + age + hinc       524
  Resid. Dev   Test    Df LR stat.   Pr(Chi)
1   516.6551                                
2   509.3430 1 vs 2     4 7.312093 0.1202864

Linear regression - sHEI

## Models
sHEI.lm.M1 <- lm(sHEI ~ sfbl.liking, data = full.df)
summary(sHEI.lm.M1) # r2 = 0.027

Call:
lm(formula = sHEI ~ sfbl.liking, data = full.df)

Residuals:
     Min       1Q   Median       3Q      Max 
-23.7366  -6.4683   0.4186   5.9963  21.3516 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 48.69103    0.85702  56.815  < 2e-16 ***
sfbl.liking -0.05453    0.01879  -2.901  0.00402 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 8.682 on 268 degrees of freedom
Multiple R-squared:  0.03045,   Adjusted R-squared:  0.02684 
F-statistic: 8.418 on 1 and 268 DF,  p-value: 0.004024

sHEI.lm.M2 <- lm(sHEI ~ sugar.intake*sfbl.liking, data = full.df)
summary(sHEI.lm.M2) # r2 = 0.125

Call:
lm(formula = sHEI ~ sugar.intake * sfbl.liking, data = full.df)

Residuals:
     Min       1Q   Median       3Q      Max 
-20.6058  -6.2422   0.3397   6.0467  21.8734 

Coefficients:
                           Estimate Std. Error t value Pr(>|t|)    
(Intercept)               5.562e+01  1.709e+00  32.553  < 2e-16 ***
sugar.intake             -2.045e-02  4.375e-03  -4.675 4.68e-06 ***
sfbl.liking              -8.594e-02  3.613e-02  -2.379   0.0181 *  
sugar.intake:sfbl.liking  1.936e-04  7.765e-05   2.493   0.0133 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 8.234 on 266 degrees of freedom
Multiple R-squared:  0.1344,    Adjusted R-squared:  0.1246 
F-statistic: 13.76 on 3 and 266 DF,  p-value: 2.271e-08

sHEI.lm.M3 <- lm(sHEI ~ sugar.intake*restraint + sfbl.liking, data = full.df)
summary(sHEI.lm.M3) # r2 = 0.129

Call:
lm(formula = sHEI ~ sugar.intake * restraint + sfbl.liking, data = full.df)

Residuals:
    Min      1Q  Median      3Q     Max 
-24.323  -5.458   0.404   5.970  23.047 

Coefficients:
                         Estimate Std. Error t value Pr(>|t|)    
(Intercept)            48.6841904  4.3334122  11.235   <2e-16 ***
sugar.intake           -0.0154677  0.0082017  -1.886   0.0604 .  
restraint               0.2410752  0.3047423   0.791   0.4296    
sfbl.liking            -0.0127922  0.0198380  -0.645   0.5196    
sugar.intake:restraint  0.0003753  0.0005819   0.645   0.5195    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 8.213 on 265 degrees of freedom
Multiple R-squared:  0.142, Adjusted R-squared:  0.129 
F-statistic: 10.96 on 4 and 265 DF,  p-value: 3.055e-08

sHEI.lm.M4 <- lm(sHEI ~ sugar.intake + sfbl.liking*restraint, data = full.df)
summary(sHEI.lm.M4) # r2 = 0.128

Call:
lm(formula = sHEI ~ sugar.intake + sfbl.liking * restraint, data = full.df)

Residuals:
     Min       1Q   Median       3Q      Max 
-24.0706  -5.5309   0.3728   5.8164  22.5758 

Coefficients:
                       Estimate Std. Error t value Pr(>|t|)    
(Intercept)           46.817930   3.210654  14.582  < 2e-16 ***
sugar.intake          -0.010307   0.002187  -4.713 3.95e-06 ***
sfbl.liking           -0.028680   0.071143  -0.403   0.6872    
restraint              0.378056   0.218454   1.731   0.0847 .  
sfbl.liking:restraint  0.001063   0.004833   0.220   0.8260    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 8.219 on 265 degrees of freedom
Multiple R-squared:  0.1408,    Adjusted R-squared:  0.1278 
F-statistic: 10.86 on 4 and 265 DF,  p-value: 3.642e-08

sHEI.lm.M5 <- lm(sHEI ~ sugar.intake*restraint + sfbl.liking*restraint, data = full.df)
summary(sHEI.lm.M5) # r2 = 0.126

Call:
lm(formula = sHEI ~ sugar.intake * restraint + sfbl.liking * 
    restraint, data = full.df)

Residuals:
     Min       1Q   Median       3Q      Max 
-24.3596  -5.4574   0.4031   5.9336  23.0212 

Coefficients:
                         Estimate Std. Error t value Pr(>|t|)    
(Intercept)            48.6468432  4.4054397  11.042   <2e-16 ***
sugar.intake           -0.0156616  0.0090872  -1.723    0.086 .  
restraint               0.2439834  0.3108137   0.785    0.433    
sfbl.liking            -0.0090099  0.0782501  -0.115    0.908    
sugar.intake:restraint  0.0003885  0.0006399   0.607    0.544    
restraint:sfbl.liking  -0.0002654  0.0053110  -0.050    0.960    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 8.229 on 264 degrees of freedom
Multiple R-squared:  0.142, Adjusted R-squared:  0.1257 
F-statistic: 8.738 on 5 and 264 DF,  p-value: 1.101e-07

## Model comparison
anova(sHEI.lm.M3, sHEI.lm.M1) ## significant

Analysis of Variance Table

Model 1: sHEI ~ sugar.intake * restraint + sfbl.liking
Model 2: sHEI ~ sfbl.liking
  Res.Df   RSS Df Sum of Sq      F    Pr(>F)    
1    265 17876                                  
2    268 20200 -3   -2323.7 11.482 4.217e-07 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

anova(sHEI.lm.M3, sHEI.lm.M2)

Analysis of Variance Table

Model 1: sHEI ~ sugar.intake * restraint + sfbl.liking
Model 2: sHEI ~ sugar.intake * sfbl.liking
  Res.Df   RSS Df Sum of Sq      F Pr(>F)
1    265 17876                           
2    266 18035 -1   -158.99 2.3568 0.1259

anova(sHEI.lm.M3, sHEI.lm.M4)

Analysis of Variance Table

Model 1: sHEI ~ sugar.intake * restraint + sfbl.liking
Model 2: sHEI ~ sugar.intake + sfbl.liking * restraint
  Res.Df   RSS Df Sum of Sq F Pr(>F)
1    265 17876                      
2    265 17901  0   -24.789         

anova(sHEI.lm.M3, sHEI.lm.M5)

Analysis of Variance Table

Model 1: sHEI ~ sugar.intake * restraint + sfbl.liking
Model 2: sHEI ~ sugar.intake * restraint + sfbl.liking * restraint
  Res.Df   RSS Df Sum of Sq      F Pr(>F)
1    265 17876                           
2    264 17876  1   0.16911 0.0025 0.9602

Logistic regression - sHEI

## Models
sHEI.logit.M1 <- multinom(DQ ~ sfbl.liking, data = full.df)

# weights:  9 (4 variable)
initial  value 296.625318 
final  value 292.943295 
converged

tidy(sHEI.logit.M1, 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)    1.29    0.264       0.972 0.331      0.770     2.17 
2 2       sfbl.liking    0.993   0.00557    -1.18  0.237      0.983     1.00 
3 3       (Intercept)    1.69    0.252       2.08  0.0373     1.03      2.77 
4 3       sfbl.liking    0.985   0.00555    -2.64  0.00829    0.975     0.996

summary(sHEI.logit.M1) # AIC = 593.89

Call:
multinom(formula = DQ ~ sfbl.liking, data = full.df)

Coefficients:
  (Intercept)  sfbl.liking
2   0.2563379 -0.006589123
3   0.5247267 -0.014649056

Std. Errors:
  (Intercept) sfbl.liking
2   0.2638221 0.005574394
3   0.2519171 0.005548693

Residual Deviance: 585.8866 
AIC: 593.8866 

sHEI.logit.M2 <- multinom(DQ ~ sugar.intake*sfbl.liking, data = full.df)

# weights:  15 (8 variable)
initial  value 296.625318 
iter  10 value 283.549908
final  value 283.529972 
converged

tidy(sHEI.logit.M2, exponentiate = TRUE, conf.int = TRUE)

# A tibble: 8 × 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)      4.07  0.000125   11242.   0           4.07      4.07 
2 2       sugar.intake     0.997 0.000637      -4.90 9.66e- 7    0.996     0.998
3 2       sfbl.liking      0.990 0.00766       -1.25 2.12e- 1    0.976     1.01 
4 2       sugar.intake…    1.00  0.0000146      1.68 9.36e- 2    1.00      1.00 
5 3       (Intercept)      7.40  0.000126   15852.   0           7.40      7.41 
6 3       sugar.intake     0.996 0.000652      -6.45 1.13e-10    0.995     0.997
7 3       sfbl.liking      0.983 0.00804       -2.17 3.02e- 2    0.967     0.998
8 3       sugar.intake…    1.00  0.0000161      1.98 4.81e- 2    1.00      1.00 

summary(sHEI.logit.M2) # AIC = 583.06

Call:
multinom(formula = DQ ~ sugar.intake * sfbl.liking, data = full.df)

Coefficients:
  (Intercept) sugar.intake  sfbl.liking sugar.intake:sfbl.liking
2    1.403228 -0.003120136 -0.009558481             2.455887e-05
3    2.001990 -0.004202985 -0.017438546             3.175223e-05

Std. Errors:
   (Intercept) sugar.intake sfbl.liking sugar.intake:sfbl.liking
2 0.0001248255 0.0006369709 0.007661142             1.464557e-05
3 0.0001262933 0.0006517841 0.008043924             1.606437e-05

Residual Deviance: 567.0599 
AIC: 583.0599 

sHEI.logit.M3 <- multinom(DQ ~ sugar.intake*restraint + sfbl.liking, data = full.df)

# weights:  18 (10 variable)
initial  value 296.625318 
iter  10 value 279.957554
final  value 279.818896 
converged

tidy(sHEI.logit.M3, exponentiate = TRUE, conf.int = TRUE)

# A tibble: 10 × 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.02  0.00179       9.74  2.06e- 22    1.01      1.02 
 2 2       sugar.inta…    0.998 0.00124      -1.81  7.02e-  2    0.995     1.00 
 3 2       restraint      1.07  0.0253        2.52  1.17e-  2    1.01      1.12 
 4 2       sfbl.liking    1.00  0.00629       0.159 8.73e-  1    0.989     1.01 
 5 2       sugar.inta…    1.00  0.0000944     0.332 7.40e-  1    1.00      1.00 
 6 3       (Intercept)    0.942 0.00173     -34.4   1.76e-259    0.939     0.945
 7 3       sugar.inta…    0.996 0.00143      -2.51  1.19e-  2    0.994     0.999
 8 3       restraint      1.11  0.0255        4.16  3.21e-  5    1.06      1.17 
 9 3       sfbl.liking    0.995 0.00633      -0.831 4.06e-  1    0.982     1.01 
10 3       sugar.inta…    1.00  0.000104      0.681 4.96e-  1    1.00      1.00 

summary(sHEI.logit.M3) # AIC = 579.64

Call:
multinom(formula = DQ ~ sugar.intake * restraint + sfbl.liking, 
    data = full.df)

Coefficients:
  (Intercept) sugar.intake  restraint  sfbl.liking sugar.intake:restraint
2  0.01742751 -0.002251036 0.06386268  0.001002945           3.131278e-05
3 -0.05950834 -0.003584389 0.10594164 -0.005257391           7.050097e-05

Std. Errors:
  (Intercept) sugar.intake  restraint sfbl.liking sugar.intake:restraint
2 0.001789458  0.001243459 0.02532161 0.006294143           9.439122e-05
3 0.001729371  0.001425818 0.02548002 0.006326969           1.035855e-04

Residual Deviance: 559.6378 
AIC: 579.6378 

sHEI.logit.M4 <- multinom(DQ ~ sugar.intake + sfbl.liking*restraint, data = full.df)

# weights:  18 (10 variable)
initial  value 296.625318 
iter  10 value 279.936069
final  value 279.802818 
converged

tidy(sHEI.logit.M4, exponentiate = TRUE, conf.int = TRUE)

# A tibble: 10 × 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.869  0.472      -0.297  7.67e-1    0.345     2.19 
 2 2       sugar.intake     0.998  0.000663   -2.77   5.58e-3    0.997     0.999
 3 2       sfbl.liking      0.999  0.0173     -0.0407 9.68e-1    0.966     1.03 
 4 2       restraint        1.08   0.0392      1.93   5.40e-2    0.999     1.16 
 5 2       sfbl.liking:…    1.00   0.00121     0.110  9.12e-1    0.998     1.00 
 6 3       (Intercept)      0.774  0.489      -0.525  5.99e-1    0.297     2.02 
 7 3       sugar.intake     0.997  0.000738   -3.49   4.78e-4    0.996     0.999
 8 3       sfbl.liking      0.987  0.0186     -0.719  4.72e-1    0.951     1.02 
 9 3       restraint        1.13   0.0389      3.08   2.08e-3    1.04      1.22 
10 3       sfbl.liking:…    1.00   0.00126     0.448  6.54e-1    0.998     1.00 

summary(sHEI.logit.M4) # AIC = 579.61

Call:
multinom(formula = DQ ~ sugar.intake + sfbl.liking * restraint, 
    data = full.df)

Coefficients:
  (Intercept) sugar.intake   sfbl.liking  restraint sfbl.liking:restraint
2  -0.1400477 -0.001836388 -0.0007060572 0.07550574          0.0001330216
3  -0.2566899 -0.002577185 -0.0133848630 0.11976390          0.0005651518

Std. Errors:
  (Intercept) sugar.intake sfbl.liking  restraint sfbl.liking:restraint
2   0.4722156 0.0006625792  0.01734411 0.03918077           0.001205067
3   0.4888042 0.0007378550  0.01862703 0.03889997           0.001260308

Residual Deviance: 559.6056 
AIC: 579.6056 

sHEI.logit.M5 <- multinom(DQ ~ sugar.intake*restraint + sfbl.liking*restraint, data = full.df)

# weights:  21 (12 variable)
initial  value 296.625318 
iter  10 value 281.914319
final  value 279.778155 
converged

tidy(sHEI.logit.M5, exponentiate = TRUE, conf.int = TRUE)

# A tibble: 12 × 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.01   0.00181    4.17    2.98e- 5    1.00      1.01 
 2 2       sugar.intake    0.998  0.00203   -1.10    2.70e- 1    0.994     1.00 
 3 2       restraint       1.07   0.0255     2.53    1.15e- 2    1.01      1.12 
 4 2       sfbl.liking     1.00   0.0242     0.0383  9.69e- 1    0.954     1.05 
 5 2       sugar.intak…    1.00   0.000146   0.203   8.39e- 1    1.00      1.00 
 6 2       restraint:s…    1.00   0.00170    0.00989 9.92e- 1    0.997     1.00 
 7 3       (Intercept)     0.984  0.00175   -9.40    5.49e-21    0.980     0.987
 8 3       sugar.intake    0.997  0.00215   -1.52    1.28e- 1    0.993     1.00 
 9 3       restraint       1.11   0.0256     4.00    6.39e- 5    1.05      1.16 
10 3       sfbl.liking     0.989  0.0248    -0.449   6.54e- 1    0.942     1.04 
11 3       sugar.intak…    1.00   0.000151   0.325   7.45e- 1    1.00      1.00 
12 3       restraint:s…    1.00   0.00171    0.240   8.11e- 1    0.997     1.00 

summary(sHEI.logit.M5) # AIC = 583.56

Call:
multinom(formula = DQ ~ sugar.intake * restraint + sfbl.liking * 
    restraint, data = full.df)

Coefficients:
   (Intercept) sugar.intake  restraint  sfbl.liking sugar.intake:restraint
2  0.007564414 -0.002236875 0.06438084  0.000929529           2.978888e-05
3 -0.016401919 -0.003266371 0.10239668 -0.011132947           4.909948e-05
  restraint:sfbl.liking
2          1.684975e-05
3          4.099066e-04

Std. Errors:
  (Intercept) sugar.intake  restraint sfbl.liking sugar.intake:restraint
2 0.001811969  0.002028217 0.02547127  0.02423965           0.0001464117
3 0.001745013  0.002146930 0.02561221  0.02480117           0.0001512278
  restraint:sfbl.liking
2           0.001703574
3           0.001710787

Residual Deviance: 559.5563 
AIC: 583.5563 

## Model comparison
anova(sHEI.logit.M3, sHEI.logit.M1) ## significant

Likelihood ratio tests of Multinomial Models

Response: DQ
                                   Model Resid. df Resid. Dev   Test    Df
1                            sfbl.liking       536   585.8866             
2 sugar.intake * restraint + sfbl.liking       530   559.6378 1 vs 2     6
  LR stat.      Pr(Chi)
1                      
2  26.2488 0.0002000909

anova(sHEI.logit.M3, sHEI.logit.M2) ## significant

Likelihood ratio tests of Multinomial Models

Response: DQ
                                   Model Resid. df Resid. Dev   Test    Df
1             sugar.intake * sfbl.liking       532   567.0599             
2 sugar.intake * restraint + sfbl.liking       530   559.6378 1 vs 2     2
  LR stat.    Pr(Chi)
1                    
2 7.422151 0.02445121

anova(sHEI.logit.M3, sHEI.logit.M4)

Likelihood ratio tests of Multinomial Models

Response: DQ
                                   Model Resid. df Resid. Dev   Test    Df
1 sugar.intake * restraint + sfbl.liking       530   559.6378             
2 sugar.intake + sfbl.liking * restraint       530   559.6056 1 vs 2     0
    LR stat. Pr(Chi)
1                   
2 0.03215723       0

anova(sHEI.logit.M3, sHEI.logit.M5)

Likelihood ratio tests of Multinomial Models

Response: DQ
                                               Model Resid. df Resid. Dev
1             sugar.intake * restraint + sfbl.liking       530   559.6378
2 sugar.intake * restraint + sfbl.liking * restraint       528   559.5563
    Test    Df   LR stat.   Pr(Chi)
1                                  
2 1 vs 2     2 0.08148264 0.9600775

Differences - by sweet foods and beverages liking

## BMI by sweet foods and beverages liking
SFBL.BMI.aov <- full.df %>%
  aov(BMI ~ as.factor(SFBLcat), .)
summary(SFBL.BMI.aov)

                    Df Sum Sq Mean Sq F value  Pr(>F)   
as.factor(SFBLcat)   2    634   316.9    6.21 0.00231 **
Residuals          267  13627    51.0                   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

TukeyHSD(SFBL.BMI.aov)

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

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

$`as.factor(SFBLcat)`
        diff        lwr      upr     p adj
2-1 1.558401 -0.9515580 4.068361 0.3102843
3-1 3.736040  1.2260806 6.245999 0.0015309
3-2 2.177639 -0.3323207 4.687598 0.1036452

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

full.df %>%
  group_by(SFBLcat) %>%
  summarize(mean = mean(BMI), sd =sd(BMI)) %>%
  ggplot(aes(x = as.factor(SFBLcat), y = mean, fill = as.factor(SFBLcat))) +
  geom_bar(stat = "identity", position = position_dodge()) +
  labs(x = "Sweet Foods and Beverages Liking", y = "Body Mass Index") +
  geom_errorbar(aes(ymin = mean - sd/sqrt(sd), ymax = mean + sd/sqrt(sd), width = 0.25)) +
  stat_pvalue_manual(p.sig, hide.ns = T, label = "p.adj.signif", label.size = 10, bracket.size = 1, tip.length = 0.05, y.position = 40) +
  ylim(0, 45)+
  scale_fill_npg() +
  theme_bw() 

## Sugar Intake by sweet foods and beverages liking
SFBL.sugar.aov <- full.df %>%
  aov(sugar.intake ~ as.factor(SFBLcat), .)
summary(SFBL.sugar.aov)

                    Df   Sum Sq Mean Sq F value   Pr(>F)    
as.factor(SFBLcat)   2  2710389 1355195   23.83 2.99e-10 ***
Residuals          267 15182201   56862                     
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

TukeyHSD(SFBL.sugar.aov)

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

Fit: aov(formula = sugar.intake ~ as.factor(SFBLcat), data = .)

$`as.factor(SFBLcat)`
        diff       lwr      upr     p adj
2-1 102.8444  19.06467 186.6242 0.0114685
3-1 244.4000 160.62022 328.1798 0.0000000
3-2 141.5556  57.77578 225.3353 0.0002591

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

full.df %>%
  group_by(SFBLcat) %>%
  summarize(mean = mean(sugar.intake), sd =sd(sugar.intake)) %>%
  ggplot(aes(x = as.factor(SFBLcat), y = mean, fill = as.factor(SFBLcat))) +
  geom_bar(stat = "identity", position = position_dodge()) +
  labs(x = "Sweet Foods and Beverages Liking", y = "Sugar Intake (kcal/day)") +
  geom_errorbar(aes(ymin = mean - sd/sqrt(sd), ymax = mean + sd/sqrt(sd), width = 0.25)) +
    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(550, 645, 700)) +
  ylim(0, 700)+
  scale_fill_npg() +
  theme_bw() 

## Diet quality by sweet foods and beverages liking
SFBL.sHEI.aov <- full.df %>%
  aov(sHEI ~ as.factor(SFBLcat), .)
summary(SFBL.sHEI.aov)

                    Df Sum Sq Mean Sq F value Pr(>F)  
as.factor(SFBLcat)   2    607  303.43   4.005 0.0193 *
Residuals          267  20227   75.76                 
---
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 = sHEI ~ as.factor(SFBLcat), data = .)

$`as.factor(SFBLcat)`
         diff       lwr        upr     p adj
2-1 -1.005111 -4.063144  2.0529218 0.7189120
3-1 -3.561444 -6.619477 -0.5034116 0.0177075
3-2 -2.556333 -5.614366  0.5016995 0.1216253

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

full.df %>%
  group_by(SFBLcat) %>%
  summarize(mean = mean(sHEI), sd =sd(sHEI)) %>%
  ggplot(aes(x = as.factor(SFBLcat), y = mean, fill = as.factor(SFBLcat))) +
  geom_bar(stat = "identity", position = position_dodge()) +
  labs(x = "Sweet Foods and Beverages Liking", y = "Short Healthy Eating Index") +
  geom_errorbar(aes(ymin = mean - sd/sqrt(sd), ymax = mean + sd/sqrt(sd), width = 0.25)) +
  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)) +
  ylim(0, 65) +
  scale_fill_npg() +
  theme_bw() 

Differences - by BAS

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

                   Df Sum Sq Mean Sq F value  Pr(>F)   
as.factor(BAScat)   2    692   346.0   6.808 0.00131 **
Residuals         267  13569    50.8                   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

TukeyHSD(BAS.BMI.aov)

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

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

$`as.factor(BAScat)`
         diff       lwr         upr     p adj
2-1 -1.290725 -3.795331  1.21388122 0.4455769
3-1 -3.852017 -6.356623 -1.34741061 0.0010065
3-2 -2.561292 -5.065898 -0.05668582 0.0437226

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

full.df %>%
  group_by(BAScat) %>%
  summarize(mean = mean(BMI), sd =sd(BMI)) %>%
  ggplot(aes(x = as.factor(BAScat), y = mean, fill = as.factor(BAScat))) +
  geom_bar(stat = "identity", position = position_dodge()) +
  labs(x = "Body Appreciation", y = "Body Mass Index") +
  geom_errorbar(aes(ymin = mean - sd/sqrt(sd), ymax = mean + sd/sqrt(sd), width = 0.25)) +
    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(35,40)) +
  ylim(0, 45) +
  scale_fill_npg() +
  theme_bw() 

## Sugar Intake by BAS
BAS.sugar.aov <- full.df %>%
  aov(sugar.intake ~ as.factor(BAScat), .)
summary(BAS.sugar.aov)

                   Df   Sum Sq Mean Sq F value Pr(>F)  
as.factor(BAScat)   2   403953  201976   3.084 0.0474 *
Residuals         267 17488638   65501                 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

TukeyHSD(BAS.sugar.aov)

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

Fit: aov(formula = sugar.intake ~ as.factor(BAScat), data = .)

$`as.factor(BAScat)`
         diff       lwr        upr     p adj
2-1 -69.04444 -158.9631 20.8742143 0.1683335
3-1 -90.71111 -180.6298 -0.7924523 0.0474706
3-2 -21.66667 -111.5853 68.2519921 0.8372998

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

full.df %>%
  group_by(BAScat) %>%
  summarize(mean = mean(sugar.intake), sd =sd(sugar.intake)) %>%
  ggplot(aes(x = as.factor(BAScat), y = mean, fill = as.factor(BAScat))) +
  geom_bar(stat = "identity", position = position_dodge()) +
  labs(x = "Body Appreciation", y = "Sugar Intake (kcal/day)") +
  geom_errorbar(aes(ymin = mean - sd/sqrt(sd), ymax = mean + sd/sqrt(sd), width = 0.25)) +
  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(600)) +
  ylim(0, 700) +
  scale_fill_npg() +
  theme_bw() 

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

                   Df Sum Sq Mean Sq F value Pr(>F)
as.factor(BAScat)   2     38   18.85   0.242  0.785
Residuals         267  20797   77.89               

TukeyHSD(BAS.sHEI.aov)

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

Fit: aov(formula = sHEI ~ as.factor(BAScat), data = .)

$`as.factor(BAScat)`
          diff       lwr      upr     p adj
2-1 -0.8512222 -3.951980 2.249536 0.7942215
3-1 -0.1340000 -3.234758 2.966758 0.9942972
3-2  0.7172222 -2.383536 3.817980 0.8490410

full.df %>%
  group_by(BAScat) %>%
  summarize(mean = mean(sHEI), sd =sd(sHEI)) %>%
  ggplot(aes(x = as.factor(BAScat), y = mean, fill = as.factor(BAScat))) +
  geom_bar(stat = "identity", position = position_dodge()) +
  labs(x = "Body Appreciation", y = "Short Healthy Eating Index") +
  geom_errorbar(aes(ymin = mean - sd/sqrt(sd), ymax = mean + sd/sqrt(sd), width = 0.25)) +
  ylim(0, 65) +
  scale_fill_npg() +
  theme_bw() 

Differences - by uncontrolled eating

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

                  Df Sum Sq Mean Sq F value Pr(>F)  
as.factor(UEcat)   2    330  165.21   3.167 0.0437 *
Residuals        267  13930   52.17                 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

TukeyHSD(UE.BMI.aov)

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

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

$`as.factor(UEcat)`
         diff         lwr      upr     p adj
2-1 -0.379705 -2.91745712 2.158047 0.9337705
3-1  2.133674 -0.40407840 4.671426 0.1187341
3-2  2.513379 -0.02437342 5.051131 0.0528850

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

full.df %>%
  group_by(UEcat) %>%
  summarize(mean = mean(BMI), sd =sd(BMI)) %>%
  ggplot(aes(x = as.factor(UEcat), y = mean, fill = as.factor(UEcat))) +
  geom_bar(stat = "identity", position = position_dodge()) +
  labs(x = "Uncontrolled Eating", y = "Body Mass Index") +
  geom_errorbar(aes(ymin = mean - sd/sqrt(sd), ymax = mean + sd/sqrt(sd), width = 0.25)) +
  scale_fill_npg() +
  theme_bw() 

## Sugar Intake by uncontrolled eating
UE.sugar.aov <- full.df %>%
  aov(sugar.intake ~ as.factor(UEcat), .)
summary(UE.sugar.aov)

                  Df   Sum Sq Mean Sq F value Pr(>F)  
as.factor(UEcat)   2   525332  262666   4.038 0.0187 *
Residuals        267 17367259   65046                 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

TukeyHSD(UE.sugar.aov)

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

Fit: aov(formula = sugar.intake ~ as.factor(UEcat), data = .)

$`as.factor(UEcat)`
         diff       lwr      upr     p adj
2-1  54.60000 -35.00608 144.2061 0.3238493
3-1 108.04444  18.43837 197.6505 0.0133601
3-2  53.44444 -36.16163 143.0505 0.3393802

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

full.df %>%
  group_by(UEcat) %>%
  summarize(mean = mean(sugar.intake), sd =sd(sugar.intake)) %>%
  ggplot(aes(x = as.factor(UEcat), y = mean, fill = as.factor(UEcat))) +
  geom_bar(stat = "identity", position = position_dodge()) +
  labs(x = "Uncontrolled Eating", y = "Sugar Intake (kcal/day)") +
  geom_errorbar(aes(ymin = mean - sd/sqrt(sd), ymax = mean + sd/sqrt(sd), width = 0.25)) +
  stat_pvalue_manual(p.sig, hide.ns = T, label = "p.adj.signif", label.size = 10, bracket.size = 1, tip.length = 0.05, y.position = 600) +
  ylim(0, 700)+
  scale_fill_npg() +
  theme_bw() 

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

                  Df Sum Sq Mean Sq F value Pr(>F)
as.factor(UEcat)   2     49   24.45   0.314  0.731
Residuals        267  20785   77.85               

TukeyHSD(UE.sHEI.aov)

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

Fit: aov(formula = sHEI ~ as.factor(UEcat), data = .)

$`as.factor(UEcat)`
          diff       lwr      upr     p adj
2-1 -0.6180000 -3.717923 2.481923 0.8854979
3-1  0.4181111 -2.681812 3.518034 0.9458325
3-2  1.0361111 -2.063812 4.136034 0.7108851

full.df %>%
  group_by(UEcat) %>%
  summarize(mean = mean(sHEI), sd =sd(sHEI)) %>%
  ggplot(aes(x = as.factor(UEcat), y = mean, fill = as.factor(UEcat))) +
  geom_bar(stat = "identity", position = position_dodge()) +
  labs(x = "Uncontrolled Eating", y = "Short Healthy Eating Index") +
  geom_errorbar(aes(ymin = mean - sd/sqrt(sd), ymax = mean + sd/sqrt(sd), width = 0.25)) +
  ylim(0, 65) +
  scale_fill_npg() +
  theme_bw() 

Differences - by emotional eating

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

                  Df Sum Sq Mean Sq F value  Pr(>F)   
as.factor(EMcat)   2    497  248.43   4.819 0.00879 **
Residuals        267  13764   51.55                   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

TukeyHSD(EM.BMI.aov)

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

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

$`as.factor(EMcat)`
         diff        lwr      upr     p adj
2-1 0.7594191 -1.7631255 3.281964 0.7580684
3-1 3.1812442  0.6586997 5.703789 0.0090215
3-2 2.4218252 -0.1007193 4.944370 0.0629252

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

full.df %>%
  group_by(EMcat) %>%
  summarize(mean = mean(BMI), sd =sd(BMI)) %>%
  ggplot(aes(x = as.factor(EMcat), y = mean, fill = as.factor(EMcat))) +
  geom_bar(stat = "identity", position = position_dodge()) +
  labs(x = "Emotional Eating", y = "Body Mass Index") +
  geom_errorbar(aes(ymin = mean - sd/sqrt(sd), ymax = mean + sd/sqrt(sd), width = 0.25)) +
  stat_pvalue_manual(p.sig, hide.ns = T, label = "p.adj.signif", label.size = 10, bracket.size = 1, tip.length = 0.05, y.position = 40) +
  ylim(0, 45)+
  scale_fill_npg() +
  theme_bw() 

## Sugar Intake by emotional eating
EM.sugar.aov <- full.df %>%
  aov(sugar.intake ~ as.factor(EMcat), .)
summary(EM.sugar.aov)

                  Df   Sum Sq Mean Sq F value Pr(>F)
as.factor(EMcat)   2   180827   90414   1.363  0.258
Residuals        267 17711763   66336               

TukeyHSD(EM.sugar.aov)

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

Fit: aov(formula = sugar.intake ~ as.factor(EMcat), data = .)

$`as.factor(EMcat)`
         diff       lwr      upr     p adj
2-1 -9.244444 -99.73489  81.2460 0.9685509
3-1 49.688889 -40.80156 140.1793 0.3996984
3-2 58.933333 -31.55711 149.4238 0.2762422

full.df %>%
  group_by(EMcat) %>%
  summarize(mean = mean(sugar.intake), sd =sd(sugar.intake)) %>%
  ggplot(aes(x = as.factor(EMcat), y = mean, fill = as.factor(EMcat))) +
  geom_bar(stat = "identity", position = position_dodge()) +
  labs(x = "Emotional Eating", y = "Sugar Intake (kcal/day)") +
  geom_errorbar(aes(ymin = mean - sd/sqrt(sd), ymax = mean + sd/sqrt(sd), width = 0.25)) +
  ylim(0, 700)+
  scale_fill_npg() +
  theme_bw() 

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

                  Df Sum Sq Mean Sq F value Pr(>F)
as.factor(EMcat)   2     74   37.19   0.478   0.62
Residuals        267  20760   77.75               

TukeyHSD(EM.sHEI.aov)

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

Fit: aov(formula = sHEI ~ as.factor(EMcat), data = .)

$`as.factor(EMcat)`
          diff       lwr      upr     p adj
2-1  1.2855556 -1.812468 4.383579 0.5914085
3-1  0.6532222 -2.444801 3.751245 0.8728309
3-2 -0.6323333 -3.730357 2.465690 0.8803269

full.df %>%
  group_by(EMcat) %>%
  summarize(mean = mean(sHEI), sd =sd(sHEI)) %>%
  ggplot(aes(x = as.factor(EMcat), y = mean, fill = as.factor(EMcat))) +
  geom_bar(stat = "identity", position = position_dodge()) +
  labs(x = "Emotional Eating", y = "Short Healthy Eating Index") +
  geom_errorbar(aes(ymin = mean - sd/sqrt(sd), ymax = mean + sd/sqrt(sd), width = 0.25)) +
  ylim(0, 65) +
  scale_fill_npg() +
  theme_bw() 

Differences - by cognitive restraint

## BMI by cognitive restraint
CR.BMI.aov <- full.df %>%
  aov(BMI ~ as.factor(CRcat), .)
summary(CR.BMI.aov)

                  Df Sum Sq Mean Sq F value Pr(>F)  
as.factor(CRcat)   2    450  225.09   4.352 0.0138 *
Residuals        267  13810   51.72                 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

TukeyHSD(CR.BMI.aov)

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

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

$`as.factor(CRcat)`
        diff        lwr      upr     p adj
2-1 1.890073 -0.6367470 4.416893 0.1841530
3-1 3.141299  0.6144794 5.668120 0.0102636
3-2 1.251226 -1.2755937 3.778046 0.4739219

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

full.df %>%
  group_by(CRcat) %>%
  summarize(mean = mean(BMI), sd =sd(BMI)) %>%
  ggplot(aes(x = as.factor(CRcat), y = mean, fill = as.factor(CRcat))) +
  geom_bar(stat = "identity", position = position_dodge()) +
  labs(x = "Cognitive Restraint", y = "Body Mass Index") +
  geom_errorbar(aes(ymin = mean - sd/sqrt(sd), ymax = mean + sd/sqrt(sd), width = 0.25)) +
  stat_pvalue_manual(p.sig, hide.ns = T, label = "p.adj.signif", label.size = 10, bracket.size = 1, tip.length = 0.05, y.position = 40) +
  ylim(0, 45)+
  scale_fill_npg() +
  theme_bw() 

## Sugar Intake by cognitive restraint
CR.sugar.aov <- full.df %>%
  aov(sugar.intake ~ as.factor(CRcat), .)
summary(CR.sugar.aov)

                  Df   Sum Sq Mean Sq F value Pr(>F)
as.factor(CRcat)   2   274010  137005   2.076  0.127
Residuals        267 17618580   65987               

TukeyHSD(CR.sugar.aov)

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

Fit: aov(formula = sugar.intake ~ as.factor(CRcat), data = .)

$`as.factor(CRcat)`
          diff        lwr      upr     p adj
2-1 -63.844444 -154.09654 26.40765 0.2197534
3-1 -70.777778 -161.02987 19.47431 0.1560535
3-2  -6.933333  -97.18543 83.31876 0.9820911

full.df %>%
  group_by(CRcat) %>%
  summarize(mean = mean(sugar.intake), sd =sd(sugar.intake)) %>%
  ggplot(aes(x = as.factor(CRcat), y = mean, fill = as.factor(CRcat))) +
  geom_bar(stat = "identity", position = position_dodge()) +
  labs(x = "Cognitive Restraint", y = "Sugar Intake (kcal/day)") +
  geom_errorbar(aes(ymin = mean - sd/sqrt(sd), ymax = mean + sd/sqrt(sd), width = 0.25)) +
  ylim(0, 700)+
  scale_fill_npg() +
  theme_bw() 

## Diet quality by cognitive restraint
CR.sHEI.aov <- full.df %>%
  aov(sHEI ~ as.factor(CRcat), .)
summary(CR.sHEI.aov)

                  Df Sum Sq Mean Sq F value  Pr(>F)   
as.factor(CRcat)   2    808   404.2   5.389 0.00508 **
Residuals        267  20026    75.0                   
---
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 = sHEI ~ as.factor(CRcat), data = .)

$`as.factor(CRcat)`
         diff        lwr      upr     p adj
2-1 3.4662222  0.4234620 6.508982 0.0209850
3-1 3.8455556  0.8027953 6.888316 0.0088460
3-2 0.3793333 -2.6634269 3.422094 0.9535339

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

full.df %>%
  group_by(CRcat) %>%
  summarize(mean = mean(sHEI), sd =sd(sHEI)) %>%
  ggplot(aes(x = as.factor(CRcat), y = mean, fill = as.factor(CRcat))) +
  geom_bar(stat = "identity", position = position_dodge()) +
  labs(x = "Cognitive Restraint", y = "Short Healthy Eating Index") +
  geom_errorbar(aes(ymin = mean - sd/sqrt(sd), ymax = mean + sd/sqrt(sd), width = 0.25)) +
  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(55,62)) +
  ylim(0, 65) +
  scale_fill_npg() +
  theme_bw()