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

Load packages and call data

• Load packages and call Qualtrics data file
• I 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_31MAY24.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: 225
alpha: 0.702

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: 234
alpha: 0.588

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: 234
alpha: 0.476

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: 211
alpha: 0.618

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: 167
alpha: 0.78

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: 234
alpha: 0.769

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: 233
alpha: 0.553

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: 218
alpha: 0.625

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: 233
alpha: 0.523

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: 215
alpha: 0.594

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: 237
alpha: 0.779

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: 206
alpha: 0.622

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

## 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: 206
alpha: 0.785

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: 179
alpha: 0.4

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: 226
alpha: 0.891

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-47_3", "FL-48_3":"FL-52_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, 2)) %>%
  mutate(UEcat = ntile(uncontrolled, 3)) %>%
  mutate(CRcat = ntile(restraint, 3)) %>%
  mutate(EMcat = ntile(emotional, 3)) %>%
  mutate(SFBLcat = ntile(sfbl.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)

test <- full.df %>% 
  select(c(SugarCat, sugar.intake))

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.8975, p-value = 1.254e-11

## 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.83442, p-value = 3.46e-15

## 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.98145, p-value = 0.003401

## 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.98926, p-value = 0.07528

## 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.98247, p-value = 0.005017

## 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.94634, p-value = 1.162e-07

## 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.97721, p-value = 0.0007262

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

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

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

    Shapiro-Wilk normality test

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

## 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.99364, p-value = 0.4115

Demographics

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

anova_test(full.df, BMI ~ BMIcat)

ANOVA Table (type II tests)

  Effect DFn DFd       F        p p<.05   ges
1 BMIcat   3 233 159.893 2.73e-56     * 0.673

Differences - by BMI

full.df$BMIcat <- fct_relevel(full.df$BMIcat, "under", "healthy", "overweight", "obese")

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

                   Df Sum Sq Mean Sq F value Pr(>F)
as.factor(BMIcat)   3    276   92.08   1.205  0.309
Residuals         233  17812   76.45               

tidy(TukeyHSD(sHEI.BMI.aov))  

# A tibble: 6 × 7
  term              contrast  null.value estimate conf.low conf.high adj.p.value
  <chr>             <chr>          <dbl>    <dbl>    <dbl>     <dbl>       <dbl>
1 as.factor(BMIcat) healthy-…          0  -2.44      -8.49      3.61       0.723
2 as.factor(BMIcat) overweig…          0  -3.99     -10.4       2.41       0.373
3 as.factor(BMIcat) obese-un…          0  -3.89     -10.3       2.56       0.404
4 as.factor(BMIcat) overweig…          0  -1.55      -5.23      2.14       0.698
5 as.factor(BMIcat) obese-he…          0  -1.45      -5.23      2.33       0.754
6 as.factor(BMIcat) obese-ov…          0   0.0988    -4.22      4.42       1.00 

full.df %>%
  group_by(BMIcat) %>%
  summarize(mean = mean(sHEI), sd=sd(sHEI), n = n()) %>%
  ggplot(aes(x = BMIcat, y = mean, fill = BMIcat)) +
  geom_bar(stat = "identity", position = position_dodge()) +
  labs(title = "BMI and sHEI", x = "BMI Category", y = "Healthy Eating Index") +
  geom_errorbar(aes(ymin = mean - sd/sqrt(n), ymax = mean + sd/sqrt(n), width = 0.25)) +
     scale_fill_npg() +
  theme_bw() 

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

             Df Sum Sq Mean Sq F value   Pr(>F)    
BMIcat        3   1094   364.7   6.454 0.000325 ***
Residuals   233  13165    56.5                     
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

p.sig <- full.df %>%
  tukey_hsd(BAS ~ BMIcat) %>% 
  add_significance() %>%
  rename(., BMIcat = term) 

full.df %>%
  group_by(BMIcat) %>%
  summarize(mean = mean(BAS), sd=sd(BAS), n = n()) %>%
  ggplot(aes(x = BMIcat, y = mean, fill = BMIcat)) +
  geom_bar(stat = "identity", position = position_dodge()) +
  labs(title = "BAS and BMI", x = "BMI Category", y = "Body Appreciation") +
  geom_errorbar(aes(ymin = mean - sd/sqrt(n), ymax = mean + sd/sqrt(n), 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(45,52)) +
    ylim(0,60) +
   scale_fill_npg() +
  theme_bw() 

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

                   Df Sum Sq Mean Sq F value  Pr(>F)   
as.factor(BMIcat)   3    256   85.37   4.054 0.00781 **
Residuals         233   4907   21.06                   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

p.sig <- full.df %>%
  tukey_hsd(uncontrolled ~ BMIcat) %>% 
  add_significance() %>%
  rename(., BMIcat = term) 

full.df %>%
  group_by(BMIcat) %>%
  summarize(mean = mean(uncontrolled), sd=sd(uncontrolled), n = n()) %>%
  ggplot(aes(x = BMIcat, y = mean, fill = BMIcat)) +
  geom_bar(stat = "identity", position = position_dodge()) +
  labs(title = "Uncontrolled eating and BMI", x = "BMI category", y = "Uncontrolled eating") +
  geom_errorbar(aes(ymin = mean - sd/sqrt(n), ymax = mean + sd/sqrt(n), 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(25)) +
    ylim(0, 30) +
    scale_fill_npg() +
  theme_bw() 

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

                   Df Sum Sq Mean Sq F value   Pr(>F)    
as.factor(BMIcat)   3  224.6   74.86   6.385 0.000355 ***
Residuals         233 2731.7   11.72                     
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

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

full.df %>%
  group_by(BMIcat) %>%
  summarize(mean = mean(restraint), sd=sd(restraint), n = n()) %>%
  ggplot(aes(x = BMIcat, y = mean, fill = BMIcat)) +
  geom_bar(stat = "identity", position = position_dodge()) +
  labs(title = "Restrictive eating and BMI", x = "BMI category", y = "Restrictive eating") +
      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(20, 24, 28)) +  
    scale_fill_npg() +
  ylim(0, 30) +
  geom_errorbar(aes(ymin = mean - sd/sqrt(n), ymax = mean + sd/sqrt(n), width = 0.25)) +
  theme_bw()

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

                   Df Sum Sq Mean Sq F value   Pr(>F)    
as.factor(BMIcat)   3    357  118.97   5.784 0.000788 ***
Residuals         233   4793   20.57                     
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

p.sig <- full.df %>%
  tukey_hsd(emotional ~ BMIcat) %>% 
  add_significance() %>%
  rename(., BMIcat = term)

full.df %>%
  group_by(BMIcat) %>%
  summarize(mean = mean(emotional), sd=sd(emotional), n = n()) %>%
  ggplot(aes(x = BMIcat, y = mean, fill = BMIcat)) +
  geom_bar(stat = "identity", position = position_dodge()) +
  labs(title = "Emotional eating and BMI", x = "BMI category", y = "Emotional eating") +
  geom_errorbar(aes(ymin = mean - sd/sqrt(n), ymax = mean + sd/sqrt(n), 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(20)) +
     scale_fill_npg() +
  ylim(0, 30) +
  theme_bw() 

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

                   Df Sum Sq Mean Sq F value  Pr(>F)   
as.factor(BMIcat)   3   9680    3227   4.118 0.00717 **
Residuals         233 182553     783                   
---
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 = sfbl.liking ~ as.factor(BMIcat), data = .)

$`as.factor(BMIcat)`
                        diff        lwr      upr     p adj
healthy-under      -5.146334 -24.515276 14.22261 0.9017813
overweight-under    7.435551 -13.056694 27.92780 0.7839165
obese-under         8.690345 -11.970693 29.35138 0.6969730
overweight-healthy 12.581885   0.779188 24.38458 0.0316205
obese-healthy      13.836679   1.743290 25.93007 0.0176995
obese-overweight    1.254794 -12.566429 15.07602 0.9954247

full.df %>%
  group_by(BMIcat) %>%
  summarize(mean = mean(sfbl.liking), sd=sd(sfbl.liking), n = n()) %>%
  ggplot(aes(x = BMIcat, y = mean, fill = BMIcat)) +
  geom_bar(stat = "identity", position = position_dodge()) +
  labs(title = "Sweet foods and beverages liking and BMI", x = "BMI category", y = "Sweet foods and beverages liking") +
  geom_errorbar(aes(ymin = mean - sd/sqrt(n), ymax = mean + sd/sqrt(n), width = 0.25)) +
     scale_fill_npg() +
  theme_bw() 

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

                   Df Sum Sq Mean Sq F value Pr(>F)
as.factor(BMIcat)   3   1851   617.2   0.805  0.492
Residuals         233 178614   766.6               

TukeyHSD(UHF.BMI.aov)

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

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

$`as.factor(BMIcat)`
                        diff        lwr      upr     p adj
healthy-under       8.123086 -11.035734 27.28191 0.6916422
overweight-under   11.927851  -8.342085 32.19779 0.4255488
obese-under         9.127123 -11.309776 29.56402 0.6553833
overweight-healthy  3.804765  -7.869891 15.47942 0.8336926
obese-healthy       1.004037 -10.958158 12.96623 0.9963736
obese-overweight   -2.800728 -16.472014 10.87056 0.9516931

full.df %>%
  group_by(BMIcat) %>%
  summarize(mean = mean(unhealthyfat.liking), sd=sd(unhealthyfat.liking), n = n()) %>%
  ggplot(aes(x = BMIcat, y = mean, fill = BMIcat)) +
  geom_bar(stat = "identity", position = position_dodge()) +
  labs(title = "Unhealthy fat liking and BMI", x = "BMI category", y = "Unhealthy fat liking") +
  geom_errorbar(aes(ymin = mean - sd/sqrt(n), ymax = mean + sd/sqrt(n), width = 0.25)) +
     scale_fill_npg() +
  theme_bw() 

Differences - by sHEI

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

               Df Sum Sq Mean Sq F value Pr(>F)
as.factor(DQ)   2    211  105.52   2.077  0.128
Residuals     234  11890   50.81               

TukeyHSD(BMI.sHEI.aov)

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

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

$`as.factor(DQ)`
          diff       lwr       upr     p adj
2-1 -2.2975655 -4.972792 0.3776609 0.1083051
3-1 -1.3679823 -4.043209 1.3072442 0.4507185
3-2  0.9295833 -1.745643 3.6048097 0.6912526

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 = "BMI and sHEI", x = "Diet Quality", y = "Body Mass Index") +
  geom_errorbar(aes(ymin = mean - sd/sqrt(n), ymax = mean + sd/sqrt(n), width = 0.25)) +
  theme_bw() +
  scale_fill_manual(values = c("#F44336", "#FFD54F", "#66BB6A"))

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

               Df Sum Sq Mean Sq F value Pr(>F)
as.factor(DQ)   2    127   63.57   1.053  0.351
Residuals     234  14132   60.39               

TukeyHSD(BAS.sHEI.aov)

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

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

$`as.factor(DQ)`
          diff       lwr      upr     p adj
2-1  1.0506329 -1.865950 3.967215 0.6724990
3-1 -0.7341772 -3.650760 2.182405 0.8236092
3-2 -1.7848101 -4.701393 1.131772 0.3203912

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 = "BAS and sHEI", 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", "#FFD54F", "#66BB6A"))

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

               Df Sum Sq Mean Sq F value Pr(>F)
as.factor(DQ)   2      6   3.245   0.147  0.863
Residuals     234   5156  22.035               

TukeyHSD(UE.sHEI.aov)

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

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

$`as.factor(DQ)`
          diff       lwr      upr     p adj
2-1  0.4050633 -1.356653 2.166779 0.8504878
3-1  0.1898734 -1.571842 1.951589 0.9650067
3-2 -0.2151899 -1.976906 1.546526 0.9552823

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 = "Uncontrolled eating and sHEI", 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", "#FFD54F", "#66BB6A"))

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

               Df Sum Sq Mean Sq F value  Pr(>F)   
as.factor(DQ)   2  159.9   79.95    6.69 0.00149 **
Residuals     234 2796.4   11.95                   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

TukeyHSD(CR.sHEI.aov)

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

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

$`as.factor(DQ)`
         diff         lwr      upr     p adj
2-1 1.3670886  0.06969748 2.664480 0.0362213
3-1 1.9620253  0.66463418 3.259416 0.0012686
3-2 0.5949367 -0.70245442 1.892328 0.5264129

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 = "Restrictive eating and sHEI", 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", "#FFD54F", "#66BB6A"))

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

               Df Sum Sq Mean Sq F value Pr(>F)
as.factor(DQ)   2     37   18.38   0.841  0.433
Residuals     234   5113   21.85               

TukeyHSD(EM.sHEI.aov)

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

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

$`as.factor(DQ)`
         diff       lwr      upr     p adj
2-1  0.000000 -1.754309 1.754309 1.0000000
3-1 -0.835443 -2.589752 0.918866 0.5006917
3-2 -0.835443 -2.589752 0.918866 0.5006917

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 = "Emotional eating and sHEI", 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", "#FFD54F",  "#66BB6A"))

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

               Df Sum Sq Mean Sq F value Pr(>F)  
as.factor(DQ)   2   5145  2572.5   3.218 0.0418 *
Residuals     234 187089   799.5                 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

TukeyHSD(SFBL.sHEI.aov)

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

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

$`as.factor(DQ)`
          diff       lwr       upr     p adj
2-1  -3.052567 -13.66450  7.559366 0.7762447
3-1 -11.049955 -21.66189 -0.438022 0.0390832
3-2  -7.997388 -18.60932  2.614545 0.1793671

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

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

               Df Sum Sq Mean Sq F value Pr(>F)
as.factor(DQ)   2   1465   732.3   0.957  0.385
Residuals     234 179001   765.0               

TukeyHSD(UHF.sHEI.aov)

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

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

$`as.factor(DQ)`
         diff        lwr       upr     p adj
2-1  2.635443  -7.744589 13.015475 0.8208504
3-1 -3.436076 -13.816108  6.943956 0.7151903
3-2 -6.071519 -16.451551  4.308513 0.3532101

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

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

               Df Sum Sq Mean Sq F value   Pr(>F)    
as.factor(DQ)   2  16689    8345   8.976 0.000175 ***
Residuals     234 217546     930                     
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

TukeyHSD(HF.sHEI.aov)

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

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

$`as.factor(DQ)`
         diff       lwr      upr     p adj
2-1  6.701477 -4.741699 18.14465 0.3523283
3-1 20.179325  8.736149 31.62250 0.0001323
3-2 13.477848  2.034672 24.92102 0.0162328

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

Differences - by BAS

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

                   Df Sum Sq Mean Sq F value  Pr(>F)   
as.factor(BAScat)   1    359   358.8    7.18 0.00789 **
Residuals         235  11742    50.0                   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

TukeyHSD(BMI.BAS.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 -2.460794 -4.270012 -0.6515766 0.0078916

full.df %>%
  group_by(BAScat) %>%
  summarize(mean = mean(BMI), sd=sd(BMI), n = n()) %>%
  ggplot(aes(x = as.factor(BAScat), y = mean, fill = as.factor(BAScat))) +
  geom_bar(stat = "identity", position = position_dodge()) +
  labs(title = "BMI on BAS", x = "Body Appreciation", y = "Body Mass Index") +
  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"))

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

                   Df Sum Sq Mean Sq F value Pr(>F)   
as.factor(BAScat)   1    162  162.29   7.627 0.0062 **
Residuals         235   5000   21.28                  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

TukeyHSD(UE.BAS.aov)

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

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

$`as.factor(BAScat)`
         diff       lwr        upr     p adj
2-1 -1.655035 -2.835678 -0.4743917 0.0062042

full.df %>%
  group_by(BAScat) %>%
  summarize(mean = mean(uncontrolled), sd=sd(uncontrolled), n = n()) %>%
  ggplot(aes(x = as.factor(BAScat), y = mean, fill = as.factor(BAScat))) +
  geom_bar(stat = "identity", position = position_dodge()) +
  labs(title = "BAS on Uncontrolled Eating Behavior", x = "Body Appreciation", 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"))

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

                   Df Sum Sq Mean Sq F value   Pr(>F)    
as.factor(BAScat)   1  155.4  155.41   13.04 0.000373 ***
Residuals         235 2800.9   11.92                     
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

TukeyHSD(CR.BAS.aov)

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

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

$`as.factor(BAScat)`
        diff       lwr        upr     p adj
2-1 -1.61957 -2.503189 -0.7359506 0.0003729

full.df %>%
  group_by(BAScat) %>%
  summarize(mean = mean(restraint), sd=sd(restraint), n = n()) %>%
  ggplot(aes(x = as.factor(BAScat), y = mean, fill = as.factor(BAScat))) +
  geom_bar(stat = "identity", position = position_dodge()) +
  labs(title = "BAS on Restrictive Eating Behavior", x = "Body Appreciation", 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"))

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

                   Df Sum Sq Mean Sq F value   Pr(>F)    
as.factor(BAScat)   1    504   504.5   25.52 8.79e-07 ***
Residuals         235   4645    19.8                     
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

TukeyHSD(EM.BAS.aov)

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

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

$`as.factor(BAScat)`
        diff       lwr      upr p adj
2-1 -2.91796 -4.055901 -1.78002 9e-07

full.df %>%
  group_by(BAScat) %>%
  summarize(mean = mean(emotional), sd=sd(emotional), n = n()) %>%
  ggplot(aes(x = as.factor(BAScat), y = mean, fill = as.factor(BAScat))) +
  geom_bar(stat = "identity", position = position_dodge()) +
  labs(title = "BAS on Emotional Eating Behavior", x = "Body Appreciation", 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"))

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

                   Df Sum Sq Mean Sq F value Pr(>F)
as.factor(BAScat)   1    109  109.25   1.428  0.233
Residuals         235  17979   76.51               

TukeyHSD(DQ.BAS.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 -1.357911 -3.596608 0.880785 0.2332929

full.df %>%
  group_by(BAScat) %>%
  summarize(mean = mean(sHEI), sd=sd(sHEI), n = n()) %>%
  ggplot(aes(x = as.factor(BAScat), y = mean, fill = as.factor(BAScat))) +
  geom_bar(stat = "identity", position = position_dodge()) +
  labs(title = "BAS on Healthy Eating Index", x = "Body Appreciation", y = "Healthy Eating Index") +
  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 EBs

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

                  Df Sum Sq Mean Sq F value Pr(>F)
as.factor(UEcat)   2    127   63.28   1.237  0.292
Residuals        234  11974   51.17               

TukeyHSD(BMI.UE.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.2574195 -2.427295 2.942134 0.9722000
3-1 1.6627532 -1.021961 4.347467 0.3117330
3-2 1.4053337 -1.279381 4.090048 0.4339452

full.df %>%
  group_by(UEcat) %>%
  summarize(mean = mean(BMI), sd=sd(BMI), n = n()) %>%
  ggplot(aes(x = as.factor(UEcat), y = mean, fill = as.factor(UEcat))) +
  geom_bar(stat = "identity", position = position_dodge()) +
  labs(title = "Uncontrolled Eating on BMI", x = "Uncontrolled Eating", 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", "#FFD54F", "#F44336"))

    full.df %>%
      group_by(UEcat) %>%
      summarize(mean = mean(BMI), sd=sd(BMI), n = n()) %>%
      ggplot(aes(x = as.factor(UEcat), y = mean, fill = as.factor(UEcat))) +
      geom_bar(stat = "identity", position = position_dodge()) +
      labs(title = "Uncontrolled Eating on BMI", x = "Uncontrolled Eating", 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", "#FFD54F", "#F44336"))

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

                  Df Sum Sq Mean Sq F value Pr(>F)
as.factor(UEcat)   2     37   18.48    0.24  0.787
Residuals        234  18051   77.14               

TukeyHSD(sHEI.UE.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.3256962 -2.970561 3.621953 0.9705044
3-1 0.9517722 -2.344485 4.248029 0.7747575
3-2 0.6260759 -2.670181 3.922333 0.8953351

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

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

                  Df Sum Sq Mean Sq F value Pr(>F)  
as.factor(CRcat)   2    446  223.07   4.479 0.0123 *
Residuals        234  11655   49.81                 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

TukeyHSD(BMI.CR.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 2.5264459 -0.1221991 5.175091 0.0651407
3-1 3.1825961  0.5339511 5.831241 0.0137919
3-2 0.6561502 -1.9924948 3.304795 0.8286493

full.df %>%
  group_by(CRcat) %>%
  summarize(mean = mean(BMI), sd=sd(BMI), n = n()) %>%
  ggplot(aes(x = as.factor(CRcat), y = mean, fill = as.factor(CRcat))) +
  geom_bar(stat = "identity", position = position_dodge()) +
  labs(title = "Restrictive Eating on BMI", x = "Restrictive Eating", 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", "#FFD54F", "#F44336"))

## sHEI by restrictive eating
sHEI.CR.aov <- full.df %>%
  aov(sHEI ~ as.factor(CRcat), .)
summary(sHEI.CR.aov)

                  Df Sum Sq Mean Sq F value  Pr(>F)   
as.factor(CRcat)   2    912   456.0   6.212 0.00235 **
Residuals        234  17176    73.4                   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

TukeyHSD(sHEI.CR.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.401646  0.1862746 6.617017 0.0352945
3-1 4.639873  1.4245025 7.855244 0.0022486
3-2 1.238228 -1.9771431 4.453599 0.6355780

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

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

                  Df Sum Sq Mean Sq F value Pr(>F)  
as.factor(EMcat)   2    279  139.35   2.758 0.0655 .
Residuals        234  11822   50.52                 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

TukeyHSD(BMI.EM.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.749185 -1.91842018 3.416790 0.7854786
3-1 2.581548 -0.08605734 5.249153 0.0602315
3-2 1.832363 -0.83524236 4.499968 0.2391320

full.df %>%
  group_by(EMcat) %>%
  summarize(mean = mean(BMI), sd=sd(BMI), n = n()) %>%
  ggplot(aes(x = as.factor(EMcat), y = mean, fill = as.factor(EMcat))) +
  geom_bar(stat = "identity", position = position_dodge()) +
  labs(title = "Emotional Eating on BMI", x = "Emotional Eating", 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", "#FFD54F", "#F44336"))

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

                  Df Sum Sq Mean Sq F value Pr(>F)
as.factor(EMcat)   2    128   63.86   0.832  0.436
Residuals        234  17960   76.75               

TukeyHSD(sHEI.EM.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.7981013 -1.489859 5.086062 0.4022633
3-1  0.8821519 -2.405808 4.170112 0.8021938
3-2 -0.9159494 -4.203910 2.372011 0.7885306

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

Differences - by Liking

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

                    Df Sum Sq Mean Sq F value  Pr(>F)   
as.factor(SFBLcat)   2    498  249.04   5.022 0.00732 **
Residuals          234  11603   49.59                   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

TukeyHSD(BMI.SFBL.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 2.126944 -0.5157936 4.769682 0.1414081
3-1 3.526054  0.8833168 6.168792 0.0052756
3-2 1.399110 -1.2436272 4.041848 0.4257833

full.df %>%
  group_by(SFBLcat) %>%
  summarize(mean = mean(BMI), sd=sd(BMI), n = n()) %>%
  ggplot(aes(x = as.factor(SFBLcat), y = mean, fill = as.factor(SFBLcat))) +
  geom_bar(stat = "identity", position = position_dodge()) +
  labs(title = "Sweet foods and beverages liking on BMI", x = "Sweet Foods and Beverages Liking", 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", "#FFD54F", "#F44336"))

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

                    Df Sum Sq Mean Sq F value Pr(>F)  
as.factor(SFBLcat)   2    369  184.44   2.436 0.0897 .
Residuals          234  17719   75.72                 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

TukeyHSD(sHEI.SFBL.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 -0.6673418 -3.933153 2.5984695 0.8799005
3-1 -2.9163291 -6.182140 0.3494822 0.0907329
3-2 -2.2489873 -5.514799 1.0168240 0.2374108

full.df %>%
  group_by(SFBLcat) %>%
  summarize(mean = mean(sHEI), sd=sd(sHEI), n = n()) %>%
  ggplot(aes(x = as.factor(SFBLcat), y = mean, fill = as.factor(SFBLcat))) +
  geom_bar(stat = "identity", position = position_dodge()) +
  labs(title = "Sweet foods and beverages liking on sHEI", x = "Sweet Foods and Beverages Liking", y = "Healthy Eating Index") +
  geom_errorbar(aes(ymin = mean - sd/sqrt(n), ymax = mean + sd/sqrt(n), width = 0.25)) +
  theme_bw() +
  scale_fill_manual(values = c("#66BB6A", "#FFD54F", "#F44336"))

## BMI by unhealthy fats liking
BMI.UHF.aov <- full.df %>%
  aov(BMI ~ as.factor(UHFcat), .)
summary(BMI.UHF.aov)

                   Df Sum Sq Mean Sq F value Pr(>F)
as.factor(UHFcat)   2     18    8.86   0.172  0.842
Residuals         234  12083   51.64               

TukeyHSD(BMI.UHF.aov)

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

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

$`as.factor(UHFcat)`
         diff       lwr      upr     p adj
2-1 0.1038248 -2.593063 2.800713 0.9954644
3-1 0.6248719 -2.072016 3.321760 0.8483578
3-2 0.5210471 -2.175841 3.217935 0.8919097

full.df %>%
  group_by(UHFcat) %>%
  summarize(mean = mean(BMI), sd=sd(BMI), n = n()) %>%
  ggplot(aes(x = as.factor(UHFcat), y = mean, fill = as.factor(UHFcat))) +
  geom_bar(stat = "identity", position = position_dodge()) +
  labs(title = "Unhealthy fats liking on BMI", x = "Unhealthy Fats Liking", 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", "#FFD54F", "#F44336"))

## sHEI by unhealthy fats liking
sHEI.UHF.aov <- full.df %>%
  aov(sHEI ~ as.factor(UHFcat), .)
summary(sHEI.UHF.aov)

                   Df Sum Sq Mean Sq F value Pr(>F)
as.factor(UHFcat)   2    159   79.61   1.039  0.355
Residuals         234  17929   76.62               

TukeyHSD(sHEI.UHF.aov)

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

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

$`as.factor(UHFcat)`
          diff       lwr      upr     p adj
2-1 -2.0077215 -5.292797 1.277354 0.3213229
3-1 -0.9994937 -4.284570 2.285582 0.7532989
3-2  1.0082278 -2.276848 4.293304 0.7495709

full.df %>%
  group_by(UHFcat) %>%
  summarize(mean = mean(sHEI), sd=sd(sHEI), n = n()) %>%
  ggplot(aes(x = as.factor(UHFcat), y = mean, fill = as.factor(UHFcat))) +
  geom_bar(stat = "identity", position = position_dodge()) +
  labs(title = "Unhealthy fats liking on sHEI", x = "Unhealthy Fats Liking", y = "Healthy Eating Index") +
  geom_errorbar(aes(ymin = mean - sd/sqrt(n), ymax = mean + sd/sqrt(n), width = 0.25)) +
  theme_bw() +
  scale_fill_manual(values = c("#66BB6A", "#FFD54F", "#F44336"))

## BMI by healthy fats liking
BMI.HF.aov <- full.df %>%
  aov(BMI ~ as.factor(HFcat), .)
summary(BMI.HF.aov)

                  Df Sum Sq Mean Sq F value Pr(>F)
as.factor(HFcat)   2    166   83.19   1.631  0.198
Residuals        234  11935   51.00               

TukeyHSD(BMI.HF.aov)

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

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

$`as.factor(HFcat)`
          diff        lwr      upr     p adj
2-1 -1.2977455 -3.9779922 1.382501 0.4892086
3-1  0.7280626 -1.9521841 3.408309 0.7977812
3-2  2.0258080 -0.6544386 4.706055 0.1775925

full.df %>%
  group_by(HFcat) %>%
  summarize(mean = mean(BMI), sd=sd(BMI), n = n()) %>%
  ggplot(aes(x = as.factor(HFcat), y = mean, fill = as.factor(HFcat))) +
  geom_bar(stat = "identity", position = position_dodge()) +
  labs(title = "Healthy fats liking on BMI", x = "Healthy Fats Liking", y = "BMI") +
  geom_errorbar(aes(ymin = mean - sd/sqrt(n), ymax = mean + sd/sqrt(n), width = 0.25)) +
  theme_bw() +
  scale_fill_npg() 

## sHEI by healthy fats liking
sHEI.HF.aov <- full.df %>%
  aov(sHEI ~ as.factor(HFcat), .)
summary(sHEI.HF.aov)

                  Df Sum Sq Mean Sq F value  Pr(>F)    
as.factor(HFcat)   2   1147   573.3   7.919 0.00047 ***
Residuals        234  16941    72.4                    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

TukeyHSD(sHEI.HF.aov)

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

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

$`as.factor(HFcat)`
        diff         lwr      upr     p adj
2-1 3.107468 -0.08586557 6.300802 0.0584159
3-1 5.365443  2.17210911 8.558777 0.0002888
3-2 2.257975 -0.93535924 5.451309 0.2197865

full.df %>%
  group_by(HFcat) %>%
  summarize(mean = mean(sHEI), sd=sd(sHEI), n = n()) %>%
  ggplot(aes(x = as.factor(HFcat), y = mean, fill = as.factor(HFcat))) +
  geom_bar(stat = "identity", position = position_dodge()) +
  labs(title = "Unhealthy fats liking on sHEI", x = "Healthy Fats Liking", y = "Healthy Eating Index") +
  geom_errorbar(aes(ymin = mean - sd/sqrt(n), ymax = mean + sd/sqrt(n), width = 0.25)) +
  theme_bw() +
  scale_fill_npg() 

Correlations - overall corrplot

## Corr DF
corr.df <- full.df %>%
  select(age, BMI, BAS, hinc, uncontrolled, restraint, emotional, sfbl.liking, healthyfat.liking, unhealthyfat.liking, sugar.intake, sHEI) 

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

    ## Women <= 35
    mwom.corr <- corr.df %>%
      filter(age <= 35) %>%
      cor(use = "pairwise.complete.obs")
    
    mwom.res <- corr.df %>%
      filter(age <= 35) %>%
      cor.mtest(., conf.level = 0.95)
    
    corrplot(mwom.corr, p.mat = mwom.res$p, method = 'circle',tl.col = 'black', type = 'lower', insig='blank', addCoef.col ='black', number.cex = 0.7, order = 'original', diag=FALSE) 
      title("Women <= 35 years old")      

    ## Women <= 25  
    ywom.corr <- corr.df %>%
      filter(age <= 25) %>%
      cor(use = "pairwise.complete.obs")
    
    ywom.res <- corr.df %>%
      filter(age <= 25) %>%
      cor.mtest(., conf.level = 0.95)
    
    corrplot(ywom.corr, p.mat = ywom.res$p, method = 'circle',tl.col = 'black', type = 'lower', insig='blank', addCoef.col ='black', number.cex = 0.7, order = 'original', diag=FALSE) 
      title("Women <= 25 years old")      

Partial correlations - BAS x EBs

## BAS x uncontrolled
  corrfunc(full.df, BAS, uncontrolled, "BAS x uncontrolled eating, overall")

  pcor.df <- full.df%>%select(c("BAS", "uncontrolled", "age")) %>% na.omit()
  pcor.test(pcor.df$BAS,pcor.df$uncontrolled,pcor.df[,c("age")],method="pearson")

   estimate     p.value statistic   n gp  Method
1 -0.203771 0.001650385   -3.1839 237  1 pearson

  ## Low BAS
  full.df %>% filter(BAScat == "1") %>%
  corrfunc(., BAS, uncontrolled, "BAS x uncontrolled eating, low BAS")

  pcor.df <- full.df%>% filter(BAScat == "1") %>%
  select(c("BAS", "uncontrolled", "age")) %>% na.omit()
  pcor.test(pcor.df$BAS,pcor.df$uncontrolled,pcor.df[,c("age")],method="pearson")

    estimate   p.value statistic   n gp  Method
1 -0.1202421 0.1946416 -1.304511 119  1 pearson

  ## High BAS
  full.df %>% filter(BAScat == "2") %>%
  corrfunc(., BAS, uncontrolled, "BAS x uncontrolled eating, high BAS")

  pcor.df <- full.df%>% filter(BAScat == "2") %>%
  select(c("BAS", "uncontrolled", "age")) %>% na.omit()
  pcor.test(pcor.df$BAS,pcor.df$uncontrolled,pcor.df[,c("age")],method="pearson")

     estimate p.value  statistic   n gp  Method
1 -0.07503973 0.42134 -0.8069867 118  1 pearson

## BAS x restrictive
  corrfunc(full.df, BAS, restraint, "BAS x restrictive eating, overall")

  pcor.df <- full.df%>%select(c("BAS", "restraint", "age")) %>% na.omit()
  pcor.test(pcor.df$BAS,pcor.df$restraint,pcor.df[,c("age")],method="pearson")

    estimate      p.value statistic   n gp  Method
1 -0.2388758 0.0002122597 -3.763036 237  1 pearson

  ## Low BAS
  full.df %>% filter(BAScat == "1") %>%
  corrfunc(., BAS, restraint, "BAS x restrictive eating, low BAS")

  pcor.df <- full.df%>% filter(BAScat == "1") %>%
  select(c("BAS", "restraint", "age")) %>% na.omit()
  pcor.test(pcor.df$BAS,pcor.df$restraint,pcor.df[,c("age")],method="pearson")

     estimate   p.value statistic   n gp  Method
1 -0.09951174 0.2836613 -1.077121 119  1 pearson

  ## High BAS
  full.df %>% filter(BAScat == "2") %>%
  corrfunc(., BAS, restraint, "BAS x restrictive eating, high BAS")

  pcor.df <- full.df%>% filter(BAScat == "2") %>%
  select(c("BAS", "restraint", "age")) %>% na.omit()
  pcor.test(pcor.df$BAS,pcor.df$restraint,pcor.df[,c("age")],method="pearson")

     estimate   p.value  statistic   n gp  Method
1 -0.06662462 0.4754055 -0.7160604 118  1 pearson

## BAS x uncontrolled
  corrfunc(full.df, BAS, emotional, "BAS x emotional eating, overall")  

  pcor.df <- full.df%>%select(c("BAS", "emotional", "age")) %>% na.omit()
  pcor.test(pcor.df$BAS,pcor.df$emotional,pcor.df[,c("age")],method="pearson")

    estimate      p.value statistic   n gp  Method
1 -0.3699795 4.547943e-09 -6.091878 237  1 pearson

  ## Low BAS
  full.df %>% filter(BAScat == "1") %>%
  corrfunc(., BAS, emotional, "BAS x emotional eating, low BAS")

  pcor.df <- full.df%>% filter(BAScat == "1") %>%
  select(c("BAS", "emotional", "age")) %>% na.omit()
  pcor.test(pcor.df$BAS,pcor.df$emotional,pcor.df[,c("age")],method="pearson")

    estimate   p.value statistic   n gp  Method
1 -0.2089502 0.0231625 -2.301259 119  1 pearson

  ## High BAS
  full.df %>% filter(BAScat == "2") %>%
  corrfunc(., BAS, emotional, "BAS x emotional eating, high BAS")

  pcor.df <- full.df%>% filter(BAScat == "2") %>%
  select(c("BAS", "emotional", "age")) %>% na.omit()
  pcor.test(pcor.df$BAS,pcor.df$emotional,pcor.df[,c("age")],method="pearson")

    estimate    p.value statistic   n gp  Method
1 -0.2068261 0.02526042  -2.26698 118  1 pearson

Partial correlations - Sweet liking x EBs

## Sweet foods and beverages liking x uncontrolled
  corrfunc(full.df, sfbl.liking, uncontrolled, "Sweet foods and beverages liking x uncontrolled eating, overall")

  pcor.df <- full.df%>% 
  select(c("sfbl.liking", "uncontrolled", "age")) %>% na.omit()
  pcor.test(pcor.df$sfbl.liking,pcor.df$uncontrolled,pcor.df[,c("age")],method="pearson")

   estimate    p.value statistic   n gp  Method
1 0.1312968 0.04390021  2.025994 237  1 pearson

        ## Low SFBL
        full.df %>% filter(SFBLcat == "1") %>%
        corrfunc(., sfbl.liking, uncontrolled, "Sweet foods and beverages liking x uncontrolled eating, low SFBL")

        pcor.df <- full.df%>% filter(SFBLcat == "1") %>%
  select(c("sfbl.liking", "uncontrolled", "age")) %>% na.omit()
  pcor.test(pcor.df$sfbl.liking,pcor.df$uncontrolled,pcor.df[,c("age")],method="pearson")

    estimate   p.value statistic  n gp  Method
1 0.07716861 0.5018806 0.6747524 79  1 pearson

        ## Medium SFBL
        full.df %>% filter(SFBLcat == "2") %>%
        corrfunc(., sfbl.liking, uncontrolled, "Sweet foods and beverages liking x uncontrolled eating, medium SFBL")

        pcor.df <- full.df%>% filter(SFBLcat == "2") %>%
  select(c("sfbl.liking", "uncontrolled", "age")) %>% na.omit()
  pcor.test(pcor.df$sfbl.liking,pcor.df$uncontrolled,pcor.df[,c("age")],method="pearson")

     estimate   p.value  statistic  n gp  Method
1 -0.03032326 0.7921309 -0.2644737 79  1 pearson

        ## High SFBL
        full.df %>% filter(SFBLcat == "3") %>%
        corrfunc(., sfbl.liking, uncontrolled, "Sweet foods and beverages liking x uncontrolled eating, high SFBL")

        pcor.df <- full.df%>% filter(SFBLcat == "3") %>%
  select(c("sfbl.liking", "uncontrolled", "age")) %>% na.omit()
  pcor.test(pcor.df$sfbl.liking,pcor.df$uncontrolled,pcor.df[,c("age")],method="pearson")

     estimate   p.value  statistic  n gp  Method
1 -0.09949783 0.3861025 -0.8717277 79  1 pearson

## Sweet foods and beverages liking x restrictive
  corrfunc(full.df, sfbl.liking, restraint, "Sweet foods and beverages liking x restrictive eating, overall")

        pcor.df <- full.df%>% 
        select(c("sfbl.liking", "restraint", "age")) %>% na.omit()
        pcor.test(pcor.df$sfbl.liking,pcor.df$restraint,pcor.df[,c("age")],method="pearson")

     estimate   p.value  statistic   n gp  Method
1 0.006277204 0.9235832 0.09602465 237  1 pearson

        ## Low SFBL
        full.df %>% filter(SFBLcat == "1") %>%
        corrfunc(., sfbl.liking, restraint, "Sweet foods and beverages liking x restrictive eating, low SFBL")

        pcor.df <- full.df%>% filter(SFBLcat == "1") %>%
        select(c("sfbl.liking", "restraint", "age")) %>% na.omit()
        pcor.test(pcor.df$sfbl.liking,pcor.df$restraint,pcor.df[,c("age")],method="pearson")        

     estimate   p.value  statistic  n gp  Method
1 0.008594863 0.9404665 0.07493105 79  1 pearson

        ## Medium SFBL
        full.df %>% filter(SFBLcat == "2") %>%
        corrfunc(., sfbl.liking, restraint, "Sweet foods and beverages liking x restrictive eating, medium SFBL")

         pcor.df <- full.df%>% filter(SFBLcat == "2") %>%
        select(c("sfbl.liking", "restraint", "age")) %>% na.omit()
        pcor.test(pcor.df$sfbl.liking,pcor.df$restraint,pcor.df[,c("age")],method="pearson")        

    estimate   p.value statistic  n gp  Method
1 0.02513713 0.8270744 0.2192097 79  1 pearson

        ## High SFBL
        full.df %>% filter(SFBLcat == "3") %>%
        corrfunc(., sfbl.liking, restraint, "Sweet foods and beverages liking x restrictive eating, high SFBL")

        pcor.df <- full.df%>% filter(SFBLcat == "3") %>%
        select(c("sfbl.liking", "restraint", "age")) %>% na.omit()
        pcor.test(pcor.df$sfbl.liking,pcor.df$restraint,pcor.df[,c("age")],method="pearson")        

    estimate   p.value statistic  n gp  Method
1 0.09853482 0.3907394 0.8632073 79  1 pearson

## Sweet foods and beverages liking x emotional
  corrfunc(full.df, sfbl.liking, emotional, "Sweet foods and beverages liking x emotional eating, overall")  

        pcor.df <- full.df%>% 
        select(c("sfbl.liking", "emotional", "age")) %>% na.omit()
        pcor.test(pcor.df$sfbl.liking,pcor.df$emotional,pcor.df[,c("age")],method="pearson")

    estimate   p.value statistic   n gp  Method
1 0.09125631 0.1622994  1.401802 237  1 pearson

      ## Low SFBL
      full.df %>% filter(SFBLcat == "1") %>%
      corrfunc(., sfbl.liking, emotional, "Sweet foods and beverages liking x emotional eating, low SFBL")

        pcor.df <- full.df%>% filter(SFBLcat == "1") %>%
        select(c("sfbl.liking", "emotional", "age")) %>% na.omit()
        pcor.test(pcor.df$sfbl.liking,pcor.df$emotional,pcor.df[,c("age")],method="pearson")

      estimate   p.value   statistic  n gp  Method
1 -0.008345793 0.9421887 -0.07275947 79  1 pearson

      ## Medium SFBL
      full.df %>% filter(SFBLcat == "2") %>%
      corrfunc(., sfbl.liking, emotional, "Sweet foods and beverages liking x emotional eating, medium SFBL")

        pcor.df <- full.df%>% filter(SFBLcat == "2") %>%
        select(c("sfbl.liking", "emotional", "age")) %>% na.omit()
        pcor.test(pcor.df$sfbl.liking,pcor.df$emotional,pcor.df[,c("age")],method="pearson")        

    estimate p.value statistic  n gp  Method
1 0.06381328 0.57886 0.5574474 79  1 pearson

      ## High SFBL
      full.df %>% filter(SFBLcat == "3") %>%
      corrfunc(., sfbl.liking, emotional, "Sweet foods and beverages liking x emotional eating, high SFBL")

        pcor.df <- full.df%>% filter(SFBLcat == "3") %>%
        select(c("sfbl.liking", "emotional", "age")) %>% na.omit()
        pcor.test(pcor.df$sfbl.liking,pcor.df$emotional,pcor.df[,c("age")],method="pearson")        

    estimate   p.value  statistic  n gp  Method
1 -0.0577299 0.6156366 -0.5041184 79  1 pearson

Partial correlations - Unhealthy fat liking X EBs

## Unhealthy fat liking x uncontrolled
  corrfunc(full.df, unhealthyfat.liking, uncontrolled, "Unhealthy fat liking x uncontrolled eating, overall")

        pcor.df <- full.df%>% 
        select(c("unhealthyfat.liking", "uncontrolled", "age")) %>% na.omit()
        pcor.test(pcor.df$unhealthyfat.liking,pcor.df$uncontrolled,pcor.df[,c("age")],method="pearson")      

    estimate   p.value statistic   n gp  Method
1 0.02178643 0.7391706 0.3333474 237  1 pearson

        ## Low UHF
        full.df %>% filter(UHFcat == "1") %>%
        corrfunc(., unhealthyfat.liking, uncontrolled, "Unhealthy fat liking x uncontrolled eating, low UHF")

        pcor.df <- full.df%>% filter(UHFcat == "1") %>%
        select(c("unhealthyfat.liking", "uncontrolled", "age")) %>% na.omit()
        pcor.test(pcor.df$unhealthyfat.liking,pcor.df$uncontrolled,pcor.df[,c("age")],method="pearson")      

   estimate    p.value statistic  n gp  Method
1 0.1913587 0.09328927  1.699635 79  1 pearson

        ## Medium UHF
        full.df %>% filter(UHFcat == "2") %>%
        corrfunc(., unhealthyfat.liking, uncontrolled, "Unhealthy fat liking x uncontrolled eating, medium UHF")

        pcor.df <- full.df%>% filter(UHFcat == "2") %>%
        select(c("unhealthyfat.liking", "uncontrolled", "age")) %>% na.omit()
        pcor.test(pcor.df$unhealthyfat.liking,pcor.df$uncontrolled,pcor.df[,c("age")],method="pearson")  

    estimate   p.value statistic  n gp  Method
1 0.03499055 0.7610273 0.3052275 79  1 pearson

        ## High UHF
        full.df %>% filter(UHFcat == "3") %>%
        corrfunc(., unhealthyfat.liking, uncontrolled, "Unhealthy fat liking x uncontrolled eating, high UHF")

        pcor.df <- full.df%>% filter(UHFcat == "3") %>%
        select(c("unhealthyfat.liking", "uncontrolled", "age")) %>% na.omit()
        pcor.test(pcor.df$unhealthyfat.liking,pcor.df$uncontrolled,pcor.df[,c("age")],method="pearson")  

   estimate   p.value statistic  n gp  Method
1 -0.140271 0.2206109 -1.235065 79  1 pearson

## Unhealthy fat liking x restrictive
  corrfunc(full.df, unhealthyfat.liking, restraint, "Unhealthy fat liking x restrictive eating, overall")

        pcor.df <- full.df%>% 
        select(c("unhealthyfat.liking", "restraint", "age")) %>% na.omit()
        pcor.test(pcor.df$unhealthyfat.liking,pcor.df$restraint,pcor.df[,c("age")],method="pearson") 

     estimate   p.value statistic   n gp  Method
1 -0.09098121 0.1635745 -1.397541 237  1 pearson

        ## Low UHF
        full.df %>% filter(UHFcat == "1") %>%
        corrfunc(., unhealthyfat.liking, restraint, "Unhealthy fat liking x restrictive eating, low UHF")

        pcor.df <- full.df %>% filter(UHFcat == "1") %>%
        select(c("unhealthyfat.liking", "restraint", "age")) %>% na.omit()
        pcor.test(pcor.df$unhealthyfat.liking,pcor.df$restraint,pcor.df[,c("age")],method="pearson") 

     estimate   p.value  statistic  n gp  Method
1 -0.05411163 0.6379762 -0.4724264 79  1 pearson

        ## Medium UHF
        full.df %>% filter(UHFcat == "2") %>%
        corrfunc(., unhealthyfat.liking, restraint, "Unhealthy fat liking x restrictive eating, medium UHF")

        pcor.df <- full.df%>% filter(UHFcat == "2") %>%
        select(c("unhealthyfat.liking", "restraint", "age")) %>% na.omit()
        pcor.test(pcor.df$unhealthyfat.liking,pcor.df$restraint,pcor.df[,c("age")],method="pearson") 

     estimate  p.value  statistic  n gp  Method
1 -0.03989028 0.728778 -0.3480325 79  1 pearson

        ## High UHF
        full.df %>% filter(UHFcat == "3") %>%
        corrfunc(., unhealthyfat.liking, restraint, "Unhealthy fat liking x restrictive eating, high UHF")

        pcor.df <- full.df%>% filter(UHFcat == "3") %>%
        select(c("unhealthyfat.liking", "restraint", "age")) %>% na.omit()
        pcor.test(pcor.df$unhealthyfat.liking,pcor.df$restraint,pcor.df[,c("age")],method="pearson") 

   estimate   p.value statistic  n gp  Method
1 0.1845948 0.1056814  1.637399 79  1 pearson

## Unhealthy fat liking x emotional
  corrfunc(full.df, unhealthyfat.liking, emotional, "Unhealthy fat liking x emotional eating, overall")  

        pcor.df <- full.df%>% 
        select(c("unhealthyfat.liking", "emotional", "age")) %>% na.omit()
        pcor.test(pcor.df$unhealthyfat.liking,pcor.df$emotional,pcor.df[,c("age")],method="pearson") 

     estimate   p.value  statistic   n gp  Method
1 -0.06384448 0.3287738 -0.9786293 237  1 pearson

      ## Low UHF
      full.df %>% filter(UHFcat == "1") %>%
      corrfunc(., unhealthyfat.liking, emotional, "Unhealthy fat liking x emotional eating, low UHF")

        pcor.df <- full.df%>% filter(UHFcat == "1") %>%
        select(c("unhealthyfat.liking", "emotional", "age")) %>% na.omit()
        pcor.test(pcor.df$unhealthyfat.liking,pcor.df$emotional,pcor.df[,c("age")],method="pearson") 

   estimate    p.value statistic  n gp  Method
1 0.2066116 0.06953618  1.840919 79  1 pearson

      ## Medium UHF
      full.df %>% filter(UHFcat == "2") %>%
      corrfunc(., unhealthyfat.liking, emotional, "Unhealthy fat liking x emotional eating, medium UHF")

        pcor.df <- full.df%>% filter(UHFcat == "2") %>%
        select(c("unhealthyfat.liking", "emotional", "age")) %>% na.omit()
        pcor.test(pcor.df$unhealthyfat.liking,pcor.df$emotional,pcor.df[,c("age")],method="pearson") 

     estimate   p.value  statistic  n gp  Method
1 -0.03183228 0.7820364 -0.2776481 79  1 pearson

      ## High UHF
      full.df %>% filter(UHFcat == "3") %>%
      corrfunc(., unhealthyfat.liking, emotional, "Unhealthy fat liking x emotional eating, high UHF")

        pcor.df <- full.df%>% filter(UHFcat == "3") %>%
        select(c("unhealthyfat.liking", "emotional", "age")) %>% na.omit()
        pcor.test(pcor.df$unhealthyfat.liking,pcor.df$emotional,pcor.df[,c("age")],method="pearson") 

     estimate   p.value  statistic  n gp  Method
1 -0.03412445 0.7667722 -0.2976635 79  1 pearson

Partial correlations - BAS/EBs x sHEI

## BAS x sHEI
  corrfunc(full.df, BAS, sHEI, "BAS x sHEI, overall")

        pcor.df <- full.df%>% 
        select(c("BAS", "sHEI", "age")) %>% na.omit()
        pcor.test(pcor.df$BAS,pcor.df$sHEI,pcor.df[,c("age")],method="pearson") 

     estimate   p.value statistic   n gp  Method
1 -0.09364102 0.1515564 -1.438754 237  1 pearson

  ## low BAS 
    full.df %>%
      filter(BAScat == "1") %>%
      corrfunc(., BAS, sHEI, "BAS x sHEI, low")

        pcor.df <- full.df%>% filter(BAScat == "1") %>%
        select(c("BAS", "sHEI", "age")) %>% na.omit()
        pcor.test(pcor.df$BAS,pcor.df$sHEI,pcor.df[,c("age")],method="pearson") 

     estimate   p.value  statistic   n gp  Method
1 -0.06504639 0.4840494 -0.7020579 119  1 pearson

  ## high BAS 
    full.df %>%
      filter(BAScat == "2") %>%
      corrfunc(., BAS, sHEI, "BAS x sHEI, high")

        pcor.df <- full.df%>% filter(BAScat == "2") %>%
        select(c("BAS", "sHEI", "age")) %>% na.omit()
        pcor.test(pcor.df$BAS,pcor.df$sHEI,pcor.df[,c("age")],method="pearson") 

     estimate   p.value  statistic   n gp  Method
1 -0.02773346 0.7666043 -0.2975227 118  1 pearson

## Uncontrolled x sHEI
  corrfunc(full.df, uncontrolled, sHEI, "Uncontrolled eating x sHEI, overall")  

        pcor.df <- full.df%>%
        select(c("uncontrolled", "sHEI", "age")) %>% na.omit()
        pcor.test(pcor.df$uncontrolled,pcor.df$sHEI,pcor.df[,c("age")],method="pearson")

    estimate   p.value statistic   n gp  Method
1 0.04886814 0.4549506  0.748433 237  1 pearson

  ## Low uncontrolled
  full.df %>%
    filter(UEcat == "1") %>%
    corrfunc(., uncontrolled, sHEI, "Uncontrolled eating x sHEI, low")  

        pcor.df <- full.df%>% filter(UEcat == "1") %>%
        select(c("uncontrolled", "sHEI", "age")) %>% na.omit()
        pcor.test(pcor.df$uncontrolled,pcor.df$sHEI,pcor.df[,c("age")],method="pearson")

    estimate   p.value statistic  n gp  Method
1 0.03181443 0.7821556 0.2774923 79  1 pearson

  ## Medium uncontrolled
  full.df %>%
    filter(UEcat == "2") %>%
    corrfunc(., uncontrolled, sHEI, "Uncontrolled eating x sHEI, medium")  

        pcor.df <- full.df%>% filter(UEcat == "2") %>%
        select(c("uncontrolled", "sHEI", "age")) %>% na.omit()
        pcor.test(pcor.df$uncontrolled,pcor.df$sHEI,pcor.df[,c("age")],method="pearson")

    estimate   p.value statistic  n gp  Method
1 0.06773795 0.5556841 0.5918852 79  1 pearson

  ## High uncontrolled
  full.df %>%
    filter(UEcat == "3") %>%
    corrfunc(., uncontrolled, sHEI, "Uncontrolled eating x sHEI, high")  

        pcor.df <- full.df%>% filter(UEcat == "3") %>%
        select(c("uncontrolled", "sHEI", "age")) %>% na.omit()
        pcor.test(pcor.df$uncontrolled,pcor.df$sHEI,pcor.df[,c("age")],method="pearson")

      estimate   p.value   statistic  n gp  Method
1 0.0008380469 0.9941899 0.007305926 79  1 pearson

## Restrictive x sHEI
  corrfunc(full.df, restraint, sHEI, "Restrictive eating x sHEI, overall")  

        pcor.df <- full.df%>%
        select(c("restraint", "sHEI", "age")) %>% na.omit()
        pcor.test(pcor.df$restraint,pcor.df$sHEI,pcor.df[,c("age")],method="pearson")

  estimate      p.value statistic   n gp  Method
1 0.221675 0.0006032956  3.477493 237  1 pearson

  ## Low restrictive
  full.df %>%
    filter(CRcat == "1") %>%
    corrfunc(., restraint, sHEI, "Restrictive eating x sHEI, low")  

        pcor.df <- full.df%>% filter(CRcat == "1") %>%
        select(c("restraint", "sHEI", "age")) %>% na.omit()
        pcor.test(pcor.df$restraint,pcor.df$sHEI,pcor.df[,c("age")],method="pearson")

    estimate   p.value statistic  n gp  Method
1 0.07370727 0.5213094 0.6443177 79  1 pearson

  ## Medium restrictive
  full.df %>%
    filter(CRcat == "2") %>%
    corrfunc(., restraint, sHEI, "Restrictive eating x sHEI, medium")  

        pcor.df <- full.df%>% filter(CRcat == "2") %>%
        select(c("restraint", "sHEI", "age")) %>% na.omit()
        pcor.test(pcor.df$restraint,pcor.df$sHEI,pcor.df[,c("age")],method="pearson")

     estimate   p.value  statistic  n gp  Method
1 -0.06924653 0.5468955 -0.6051298 79  1 pearson

  ## High restrictive
  full.df %>%
    filter(CRcat == "3") %>%
    corrfunc(., restraint, sHEI, "Restrictive eating x sHEI, high") 

        pcor.df <- full.df%>% filter(CRcat == "3") %>%
        select(c("restraint", "sHEI", "age")) %>% na.omit()
        pcor.test(pcor.df$restraint,pcor.df$sHEI,pcor.df[,c("age")],method="pearson")

    estimate   p.value statistic  n gp  Method
1 0.08588271 0.4546819 0.7514847 79  1 pearson

## Emotional x sHEI
  corrfunc(full.df, emotional, sHEI, "Emotional eating x sHEI, overall")  

        pcor.df <- full.df%>% 
        select(c("emotional", "sHEI", "age")) %>% na.omit()
        pcor.test(pcor.df$emotional,pcor.df$sHEI,pcor.df[,c("age")],method="pearson")

    estimate   p.value  statistic   n gp  Method
1 -0.0417432 0.5233796 -0.6391053 237  1 pearson

  ## Low emotional
  full.df %>%
    filter(EMcat == "1") %>%
    corrfunc(., emotional, sHEI, "Emotional eating x sHEI, low")  

        pcor.df <- full.df%>% filter(EMcat == "1") %>%
        select(c("emotional", "sHEI", "age")) %>% na.omit()
        pcor.test(pcor.df$emotional,pcor.df$sHEI,pcor.df[,c("age")],method="pearson")

    estimate    p.value statistic  n gp  Method
1 -0.2293216 0.04342337 -2.053914 79  1 pearson

  ## Medium emotional
  full.df %>%
    filter(EMcat == "2") %>%
    corrfunc(., emotional, sHEI, "Emotional eating x sHEI, medium")  

        pcor.df <- full.df%>% filter(EMcat == "2") %>%
        select(c("emotional", "sHEI", "age")) %>% na.omit()
        pcor.test(pcor.df$emotional,pcor.df$sHEI,pcor.df[,c("age")],method="pearson")

     estimate   p.value  statistic  n gp  Method
1 -0.05154173 0.6540428 -0.4499284 79  1 pearson

  ## High emotional
  full.df %>%
    filter(EMcat == "3") %>%
    corrfunc(., emotional, sHEI, "Emotional eating x sHEI, high") 

        pcor.df <- full.df%>% filter(EMcat == "3") %>%
        select(c("emotional", "sHEI", "age")) %>% na.omit()
        pcor.test(pcor.df$emotional,pcor.df$sHEI,pcor.df[,c("age")],method="pearson")

    estimate     p.value statistic  n gp  Method
1 -0.2949195 0.008762513 -2.690727 79  1 pearson

Partial correlations - Sweet liking x sHEI

## Sweet foods and beverages liking x sHEI
  corrfunc(full.df, sfbl.liking, sHEI, "Sweet foods and beverages liking x sHEI, overall")

  pcor.df <- full.df%>% 
  select(c("sfbl.liking", "sHEI", "age")) %>% na.omit()
  pcor.test(pcor.df$sfbl.liking,pcor.df$sHEI,pcor.df[,c("age")],method="pearson")

    estimate     p.value statistic   n gp  Method
1 -0.1769469 0.006422048 -2.750164 237  1 pearson

        ## Low SFBL
        full.df %>% filter(SFBLcat == "1") %>%
        corrfunc(., sfbl.liking, sHEI, "Sweet foods and beverages liking x sHEI, low SFBL")

        pcor.df <- full.df%>% filter(SFBLcat == "1") %>%
  select(c("sfbl.liking", "sHEI", "age")) %>% na.omit()
  pcor.test(pcor.df$sfbl.liking,pcor.df$sHEI,pcor.df[,c("age")],method="pearson")

    estimate   p.value statistic  n gp  Method
1 -0.1657273 0.1470359 -1.465036 79  1 pearson

        ## Medium SFBL
        full.df %>% filter(SFBLcat == "2") %>%
        corrfunc(., sfbl.liking, sHEI, "Sweet foods and beverages liking x sHEI, medium SFBL")

        pcor.df <- full.df%>% filter(SFBLcat == "2") %>%
  select(c("sfbl.liking", "sHEI", "age")) %>% na.omit()
  pcor.test(pcor.df$sfbl.liking,pcor.df$sHEI,pcor.df[,c("age")],method="pearson")

    estimate   p.value statistic  n gp  Method
1 -0.1217889 0.2881427 -1.069694 79  1 pearson

        ## High SFBL
        full.df %>% filter(SFBLcat == "3") %>%
        corrfunc(., sfbl.liking, sHEI, "Sweet foods and beverages liking x sHEI, high SFBL")

        pcor.df <- full.df%>% filter(SFBLcat == "3") %>%
  select(c("sfbl.liking", "sHEI", "age")) %>% na.omit()
  pcor.test(pcor.df$sfbl.liking,pcor.df$sHEI,pcor.df[,c("age")],method="pearson")

     estimate   p.value  statistic  n gp  Method
1 -0.03644525 0.7514073 -0.3179335 79  1 pearson

Partial correlations - BAS/EBs x BMI

## BAS x BMI
  corrfunc(full.df, BAS, BMI, "BAS x BMI, overall")

        pcor.df <- full.df%>% 
        select(c("BAS", "age", "BMI")) %>% na.omit()
        pcor.test(pcor.df$BAS,pcor.df$BMI,pcor.df[,c("age")],method="pearson") 

    estimate     p.value statistic   n gp  Method
1 -0.3057223 1.69494e-06 -4.911827 237  1 pearson

  ## low BAS 
    full.df %>%
      filter(BAScat == "1") %>%
      corrfunc(., BAS, BMI, "BAS x BMI, low")

        pcor.df <- full.df%>% filter(BAScat == "1") %>%
        select(c("BAS", "age", "BMI")) %>% na.omit()
        pcor.test(pcor.df$BAS,pcor.df$BMI,pcor.df[,c("age")],method="pearson") 

    estimate      p.value statistic   n gp  Method
1 -0.4022132 6.336776e-06 -4.731568 119  1 pearson

  ## high BAS 
    full.df %>%
      filter(BAScat == "2") %>%
      corrfunc(., BAS, BMI, "BAS x BMI, high")

        pcor.df <- full.df%>% filter(BAScat == "2") %>%
        select(c("BAS", "age", "BMI")) %>% na.omit()
        pcor.test(pcor.df$BAS,pcor.df$BMI,pcor.df[,c("age")],method="pearson") 

    estimate   p.value statistic   n gp  Method
1 -0.1272413 0.1715903 -1.375693 118  1 pearson

## Uncontrolled x BMI
  corrfunc(full.df, uncontrolled, BMI, "Uncontrolled eating x BMI, overall")  

        pcor.df <- full.df%>%
        select(c("uncontrolled", "age", "BMI")) %>% na.omit()
        pcor.test(pcor.df$uncontrolled,pcor.df$BMI,pcor.df[,c("age")],method="pearson")

   estimate    p.value statistic   n gp  Method
1 0.1264219 0.05242741  1.949525 237  1 pearson

  ## Low uncontrolled
  full.df %>%
    filter(UEcat == "1") %>%
    corrfunc(., uncontrolled, BMI, "Uncontrolled eating x BMI, overall")  

        pcor.df <- full.df%>% filter(UEcat == "1") %>%
        select(c("uncontrolled", "age", "BMI")) %>% na.omit()
        pcor.test(pcor.df$uncontrolled,pcor.df$BMI,pcor.df[,c("age")],method="pearson")

   estimate   p.value statistic  n gp  Method
1 0.0313948 0.7849593 0.2738285 79  1 pearson

  ## Medium uncontrolled
  full.df %>%
    filter(UEcat == "2") %>%
    corrfunc(., uncontrolled, BMI, "Uncontrolled eating x BMI, overall")  

        pcor.df <- full.df%>% filter(UEcat == "2") %>%
        select(c("uncontrolled", "age", "BMI")) %>% na.omit()
        pcor.test(pcor.df$uncontrolled,pcor.df$BMI,pcor.df[,c("age")],method="pearson")

   estimate   p.value statistic  n gp  Method
1 0.0303502 0.7919503 0.2647088 79  1 pearson

  ## High uncontrolled
  full.df %>%
    filter(UEcat == "3") %>%
    corrfunc(., uncontrolled, BMI, "Uncontrolled eating x BMI, overall")  

        pcor.df <- full.df%>% filter(UEcat == "3") %>%
        select(c("uncontrolled", "age", "BMI")) %>% na.omit()
        pcor.test(pcor.df$uncontrolled,pcor.df$BMI,pcor.df[,c("age")],method="pearson")

    estimate   p.value statistic  n gp  Method
1 0.07905412 0.4914574 0.6913415 79  1 pearson

## Restrictive x BMI
  corrfunc(full.df, restraint, BMI, "Restrictive eating x BMI, overall")  

        pcor.df <- full.df%>%
        select(c("restraint", "age", "BMI")) %>% na.omit()
        pcor.test(pcor.df$restraint,pcor.df$BMI,pcor.df[,c("age")],method="pearson")

   estimate      p.value statistic   n gp  Method
1 0.2326978 0.0003116607  3.660065 237  1 pearson

  ## Low restrictive
  full.df %>%
    filter(CRcat == "1") %>%
    corrfunc(., restraint, BMI, "Restrictive eating x BMI, overall")  

        pcor.df <- full.df%>% filter(CRcat == "1") %>%
        select(c("restraint", "age", "BMI")) %>% na.omit()
        pcor.test(pcor.df$restraint,pcor.df$BMI,pcor.df[,c("age")],method="pearson")

   estimate    p.value statistic  n gp  Method
1 0.2093856 0.06579185  1.866762 79  1 pearson

  ## Medium restrictive
  full.df %>%
    filter(CRcat == "2") %>%
    corrfunc(., restraint, BMI, "Restrictive eating x BMI, overall")  

        pcor.df <- full.df%>% filter(CRcat == "2") %>%
        select(c("restraint", "age", "BMI")) %>% na.omit()
        pcor.test(pcor.df$restraint,pcor.df$BMI,pcor.df[,c("age")],method="pearson")

   estimate   p.value statistic  n gp  Method
1 0.1113047 0.3319647 0.9763986 79  1 pearson

  ## High restrictive
  full.df %>%
    filter(CRcat == "3") %>%
    corrfunc(., restraint, BMI, "Restrictive eating x BMI, overall") 

        pcor.df <- full.df%>% filter(CRcat == "3") %>%
        select(c("restraint", "age", "BMI")) %>% na.omit()
        pcor.test(pcor.df$restraint,pcor.df$BMI,pcor.df[,c("age")],method="pearson")

    estimate   p.value statistic  n gp  Method
1 0.08659505 0.4509352 0.7577646 79  1 pearson

## Emotional x BMI
  corrfunc(full.df, emotional, BMI, "Emotional eating x BMI, overall")  

        pcor.df <- full.df%>% 
        select(c("emotional", "age", "BMI")) %>% na.omit()
        pcor.test(pcor.df$emotional,pcor.df$BMI,pcor.df[,c("age")],method="pearson")

   estimate     p.value statistic   n gp  Method
1 0.2196427 0.000679085  3.443988 237  1 pearson

  ## Low emotional
  full.df %>%
    filter(EMcat == "1") %>%
    corrfunc(., emotional, BMI, "Emotional eating x BMI, overall")  

        pcor.df <- full.df%>% filter(EMcat == "1") %>%
        select(c("emotional", "age", "BMI")) %>% na.omit()
        pcor.test(pcor.df$emotional,pcor.df$BMI,pcor.df[,c("age")],method="pearson")

    estimate   p.value statistic  n gp  Method
1 0.08476021 0.4606207 0.7415911 79  1 pearson

  ## Medium emotional
  full.df %>%
    filter(EMcat == "2") %>%
    corrfunc(., emotional, BMI, "Emotional eating x BMI, overall")  

        pcor.df <- full.df%>% filter(EMcat == "2") %>%
        select(c("emotional", "age", "BMI")) %>% na.omit()
        pcor.test(pcor.df$emotional,pcor.df$BMI,pcor.df[,c("age")],method="pearson")

   estimate    p.value statistic  n gp  Method
1 0.1940106 0.08875462  1.724104 79  1 pearson

  ## High emotional
  full.df %>%
    filter(EMcat == "3") %>%
    corrfunc(., emotional, BMI, "Emotional eating x BMI, overall") 

        pcor.df <- full.df%>% filter(EMcat == "3") %>%
        select(c("emotional", "BMI", "age", "BMI")) %>% na.omit()
        pcor.test(pcor.df$emotional,pcor.df$BMI,pcor.df[,c("age")],method="pearson")

   estimate    p.value statistic  n gp  Method
1 0.2720534 0.01597199  2.464668 79  1 pearson

Linear regression - sugar intake

## Sweet liking on sugar intake 
full.df %>%
  ggplot(aes(sfbl.liking, sugar.intake)) +
  geom_point() +
  geom_smooth(method = lm) +
  stat_cor() +
  theme_bw()

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

summary(lm(sugar.intake ~ sfbl.liking, data = full.df))

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

Residuals:
    Min      1Q  Median      3Q     Max 
-488.41 -148.05  -39.66  123.78  741.97 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 328.7873    23.2077  14.167  < 2e-16 ***
sfbl.liking   3.9374     0.5119   7.692 3.97e-13 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 224.4 on 235 degrees of freedom
Multiple R-squared:  0.2011,    Adjusted R-squared:  0.1977 
F-statistic: 59.16 on 1 and 235 DF,  p-value: 3.973e-13

tidy(lm(sugar.intake ~ sfbl.liking, data = 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) 6.17e142    23.2       14.2  2.38e-33 8.59e122  4.44e162
2 sfbl.liking 5.13e  1     0.512      7.69 3.97e-13 1.87e  1  1.41e  2

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

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

summary(lm(sugar.intake ~ BMI, data = full.df))

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

Residuals:
    Min      1Q  Median      3Q     Max 
-406.66 -178.02  -44.01  127.72  824.86 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  278.460     59.359   4.691 4.61e-06 ***
BMI            7.384      2.231   3.309  0.00108 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 245.5 on 235 degrees of freedom
Multiple R-squared:  0.04452,   Adjusted R-squared:  0.04046 
F-statistic: 10.95 on 1 and 235 DF,  p-value: 0.001083

tidy(lm(sugar.intake ~ BMI, data = 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) 8.58e120     59.4       4.69 0.00000461  1.40e70  5.27e171
2 BMI         1.61e  3      2.23      3.31 0.00108     1.98e 1  1.31e  5

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

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

summary(lm(sugar.intake ~ BAS, data = full.df))

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

Residuals:
    Min      1Q  Median      3Q     Max 
-408.00 -186.66  -56.75  129.69  866.39 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  624.993     77.078   8.109 2.83e-14 ***
BAS           -4.350      2.084  -2.087   0.0379 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 248.8 on 235 degrees of freedom
Multiple R-squared:  0.01821,   Adjusted R-squared:  0.01403 
F-statistic: 4.358 on 1 and 235 DF,  p-value: 0.03792

tidy(lm(sugar.intake ~ BAS, data = 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.70e+271     77.1       8.11 2.83e-14 3.04e+205   Inf    
2 BAS         1.29e-  2      2.08     -2.09 3.79e- 2 2.13e-  4     0.783

## Uncontrolled eating on sugar intake 
full.df %>%
  ggplot(aes(uncontrolled, sugar.intake)) +
  geom_point() +
  geom_smooth(method = lm) +
  stat_cor() +
  theme_bw()

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

summary(lm(sugar.intake ~ uncontrolled, data = full.df))

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

Residuals:
    Min      1Q  Median      3Q     Max 
-386.96 -196.24  -46.32  109.68  874.11 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)   353.024     67.538   5.227  3.8e-07 ***
uncontrolled    6.072      3.472   1.749   0.0817 .  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 249.5 on 235 degrees of freedom
Multiple R-squared:  0.01284,   Adjusted R-squared:  0.008643 
F-statistic: 3.058 on 1 and 235 DF,  p-value: 0.08166

tidy(lm(sugar.intake ~ uncontrolled, data = 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.07e153     67.5       5.23 0.000000380 3.39e+95  1.27e211
2 uncontrolled 4.33e  2      3.47      1.75 0.0817      4.63e- 1  4.05e  5

## Health fat liking on sugar intake 
full.df %>%
  ggplot(aes(healthyfat.liking, sugar.intake)) +
  geom_point() +
  geom_smooth(method = lm) +
  stat_cor() +
  theme_bw()

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

summary(lm(sugar.intake ~ healthyfat.liking, data = full.df))

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

Residuals:
    Min      1Q  Median      3Q     Max 
-445.58 -193.84  -44.58  120.97  849.54 

Coefficients:
                  Estimate Std. Error t value Pr(>|t|)    
(Intercept)       512.2727    26.0068  19.698   <2e-16 ***
healthyfat.liking  -1.1238     0.5136  -2.188   0.0297 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 248.6 on 235 degrees of freedom
Multiple R-squared:  0.01996,   Adjusted R-squared:  0.01579 
F-statistic: 4.787 on 1 and 235 DF,  p-value: 0.02967

tidy(lm(sugar.intake ~ healthyfat.liking, data = 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)       3.00e+222    26.0       19.7  1.16e-51 1.68e+200 5.36e+244
2 healthyfat.liking 3.25e-  1     0.514     -2.19 2.97e- 2 1.18e-  1 8.94e-  1

## HINC on sugar intake 
full.df %>%
  ggplot(aes(hinc, sugar.intake)) +
  geom_point() +
  geom_smooth(method = lm) +
  stat_cor() +
  theme_bw()

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

summary(lm(sugar.intake ~ hinc, data = full.df))

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

Residuals:
    Min      1Q  Median      3Q     Max 
-405.04 -192.15  -43.19  117.70  859.77 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  572.970     39.559  14.484  < 2e-16 ***
hinc         -18.964      6.513  -2.911  0.00394 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 246.7 on 235 degrees of freedom
Multiple R-squared:  0.03481,   Adjusted R-squared:  0.03071 
F-statistic: 8.476 on 1 and 235 DF,  p-value: 0.003944

tidy(lm(sugar.intake ~ hinc, data = 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) 6.88e+248     39.6      14.5  2.09e-34 9.80e+214 4.84e+282
2 hinc        5.81e-  9      6.51     -2.91 3.94e- 3 1.55e- 14 2.17e-  3

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

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

Residuals:
    Min      1Q  Median      3Q     Max 
-488.41 -148.05  -39.66  123.78  741.97 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 328.7873    23.2077  14.167  < 2e-16 ***
sfbl.liking   3.9374     0.5119   7.692 3.97e-13 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 224.4 on 235 degrees of freedom
Multiple R-squared:  0.2011,    Adjusted R-squared:  0.1977 
F-statistic: 59.16 on 1 and 235 DF,  p-value: 3.973e-13

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

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

Residuals:
    Min      1Q  Median      3Q     Max 
-460.72 -145.92  -39.28  116.79  740.16 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 446.3337    73.1786   6.099 4.37e-09 ***
sfbl.liking   3.8661     0.5116   7.556 9.31e-13 ***
BAS          -3.1804     1.8786  -1.693   0.0918 .  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 223.6 on 234 degrees of freedom
Multiple R-squared:  0.2108,    Adjusted R-squared:  0.204 
F-statistic: 31.25 on 2 and 234 DF,  p-value: 9.376e-13

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

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

Residuals:
    Min      1Q  Median      3Q     Max 
-483.68 -147.99  -39.25  124.74  733.10 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)  281.3150    61.5471   4.571 7.87e-06 ***
sfbl.liking    3.8747     0.5178   7.484 1.46e-12 ***
uncontrolled   2.6314     3.1594   0.833    0.406    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 224.6 on 234 degrees of freedom
Multiple R-squared:  0.2035,    Adjusted R-squared:  0.1967 
F-statistic: 29.89 on 2 and 234 DF,  p-value: 2.756e-12

sugar.lm.M4 <- lm(sugar.intake ~ sfbl.liking + healthyfat.liking, data = full.df)
summary(sugar.lm.M4) # r2 = 0.20

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

Residuals:
    Min      1Q  Median      3Q     Max 
-450.10 -150.88  -31.11  120.55  755.54 

Coefficients:
                  Estimate Std. Error t value Pr(>|t|)    
(Intercept)       361.8356    30.8149  11.742  < 2e-16 ***
sfbl.liking         3.8492     0.5130   7.503 1.29e-12 ***
healthyfat.liking  -0.7543     0.4648  -1.623    0.106    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 223.7 on 234 degrees of freedom
Multiple R-squared:   0.21, Adjusted R-squared:  0.2033 
F-statistic:  31.1 on 2 and 234 DF,  p-value: 1.052e-12

sugar.lm.M5 <- lm(sugar.intake ~ sfbl.liking + BMI, data = full.df)
summary(sugar.lm.M5) # r2 = 0.21

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

Residuals:
    Min      1Q  Median      3Q     Max 
-483.63 -141.32  -47.53  120.50  715.56 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 226.6241    54.3870   4.167 4.35e-05 ***
sfbl.liking   3.7125     0.5198   7.142 1.15e-11 ***
BMI           4.2965     2.0718   2.074   0.0392 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 222.9 on 234 degrees of freedom
Multiple R-squared:  0.2155,    Adjusted R-squared:  0.2088 
F-statistic: 32.15 on 2 and 234 DF,  p-value: 4.626e-13

sugar.lm.M6 <- lm(sugar.intake ~ sfbl.liking*BMI, data = full.df)
summary(sugar.lm.M6) # r2 = 0.21

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

Residuals:
    Min      1Q  Median      3Q     Max 
-485.37 -142.48  -44.33  125.98  714.56 

Coefficients:
                 Estimate Std. Error t value Pr(>|t|)  
(Intercept)     274.70749  115.69244   2.374   0.0184 *
sfbl.liking       2.62501    2.36631   1.109   0.2684  
BMI               2.25507    4.80450   0.469   0.6392  
sfbl.liking:BMI   0.04500    0.09552   0.471   0.6380  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 223.3 on 233 degrees of freedom
Multiple R-squared:  0.2163,    Adjusted R-squared:  0.2062 
F-statistic: 21.43 on 3 and 233 DF,  p-value: 2.697e-12

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

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

Residuals:
    Min      1Q  Median      3Q     Max 
-482.99 -152.07  -32.93  125.34  737.08 

Coefficients:
                          Estimate Std. Error t value Pr(>|t|)    
(Intercept)              290.27129   84.20165   3.447 0.000672 ***
sfbl.liking                3.57492    1.98737   1.799 0.073341 .  
uncontrolled               2.12922    4.51158   0.472 0.637408    
sfbl.liking:uncontrolled   0.01619    0.10365   0.156 0.875978    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 225.1 on 233 degrees of freedom
Multiple R-squared:  0.2036,    Adjusted R-squared:  0.1933 
F-statistic: 19.85 on 3 and 233 DF,  p-value: 1.709e-11

sugar.lm.M8 <- lm(sugar.intake ~ sfbl.liking*hinc, data = full.df)
summary(sugar.lm.M8) # r2 = 0.23

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

Residuals:
    Min      1Q  Median      3Q     Max 
-438.82 -139.25  -36.49  116.08  740.03 

Coefficients:
                 Estimate Std. Error t value Pr(>|t|)    
(Intercept)      365.9714    54.9205   6.664 1.91e-10 ***
sfbl.liking        4.6013     1.1873   3.875 0.000138 ***
hinc              -6.4561     7.7407  -0.834 0.405107    
sfbl.liking:hinc  -0.1361     0.1726  -0.788 0.431392    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 223.7 on 233 degrees of freedom
Multiple R-squared:  0.213, Adjusted R-squared:  0.2029 
F-statistic: 21.03 on 3 and 233 DF,  p-value: 4.326e-12

sugar.lm.M9 <- lm(sugar.intake ~ sfbl.liking*BMI + sfbl.liking*uncontrolled, data = full.df)
summary(sugar.lm.M9) # r2 = 0.20

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

Residuals:
    Min      1Q  Median      3Q     Max 
-481.89 -146.77  -46.97  115.83  729.02 

Coefficients:
                           Estimate Std. Error t value Pr(>|t|)  
(Intercept)              235.122911 140.496214   1.674   0.0956 .
sfbl.liking                2.730580   2.943020   0.928   0.3545  
BMI                        2.292302   4.824041   0.475   0.6351  
uncontrolled               2.148263   4.491208   0.478   0.6329  
sfbl.liking:BMI            0.041862   0.096273   0.435   0.6641  
sfbl.liking:uncontrolled  -0.003918   0.103802  -0.038   0.9699  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 224 on 231 degrees of freedom
Multiple R-squared:  0.2177,    Adjusted R-squared:  0.2007 
F-statistic: 12.85 on 5 and 231 DF,  p-value: 4.9e-11

## Model comparison
anova(sugar.lm.M8, sugar.lm.M5)

Analysis of Variance Table

Model 1: sugar.intake ~ sfbl.liking * hinc
Model 2: sugar.intake ~ sfbl.liking + BMI
  Res.Df      RSS Df Sum of Sq F Pr(>F)
1    233 11661547                      
2    234 11624672 -1     36875         

anova(sugar.lm.M8, sugar.lm.M1)

Analysis of Variance Table

Model 1: sugar.intake ~ sfbl.liking * hinc
Model 2: sugar.intake ~ sfbl.liking
  Res.Df      RSS Df Sum of Sq      F Pr(>F)
1    233 11661547                           
2    235 11838330 -2   -176783 1.7661 0.1733

anova(sugar.lm.M8, sugar.lm.M6)

Analysis of Variance Table

Model 1: sugar.intake ~ sfbl.liking * hinc
Model 2: sugar.intake ~ sfbl.liking * BMI
  Res.Df      RSS Df Sum of Sq F Pr(>F)
1    233 11661547                      
2    233 11613609  0     47937         

Logistic regression - sugar intake

## Sweet liking on sugar intake 
full.df %>%
  ggplot(aes(sfbl.liking, SugarCat)) +
  geom_point() +
  geom_smooth(method = lm) +
  stat_cor() +
  theme_bw()

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

summary(multinom(SugarCat ~ sfbl.liking, data = full.df)) # AIC = 466.4

# weights:  9 (4 variable)
initial  value 260.371112 
final  value 229.181914 
converged

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

Coefficients:
  (Intercept) sfbl.liking
2  -0.6690451  0.02431124
3  -2.0403853  0.05540583

Std. Errors:
  (Intercept) sfbl.liking
2   0.2494740 0.006567151
3   0.3805525 0.008458095

Residual Deviance: 458.3638 
AIC: 466.3638 

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

# weights:  9 (4 variable)
initial  value 260.371112 
final  value 229.181914 
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.512   0.249       -2.68 7.32e- 3   0.314      0.835
2 2       sfbl.liking    1.02    0.00657      3.70 2.14e- 4   1.01       1.04 
3 3       (Intercept)    0.130   0.381       -5.36 8.25e- 8   0.0617     0.274
4 3       sfbl.liking    1.06    0.00846      6.55 5.73e-11   1.04       1.07 

## 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)) # AIC = 515.6

# weights:  9 (4 variable)
initial  value 260.371112 
final  value 253.804893 
converged

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

Coefficients:
  (Intercept)        BMI
2   -1.710664 0.06950201
3   -2.341029 0.09318593

Std. Errors:
  (Intercept)        BMI
2   0.7230230 0.02890935
3   0.7249557 0.02864371

Residual Deviance: 507.6098 
AIC: 515.6098 

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

# weights:  9 (4 variable)
initial  value 260.371112 
final  value 253.804893 
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.181     0.723      -2.37 0.0180    0.0438     0.746
2 2       BMI           1.07      0.0289      2.40 0.0162    1.01       1.13 
3 3       (Intercept)   0.0962    0.725      -3.23 0.00124   0.0232     0.398
4 3       BMI           1.10      0.0286      3.25 0.00114   1.04       1.16 

## 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)) # AIC = 520.1

# weights:  9 (4 variable)
initial  value 260.371112 
final  value 256.028871 
converged

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

Coefficients:
  (Intercept)         BAS
2    1.843504 -0.05004382
3    2.128210 -0.05814438

Std. Errors:
  (Intercept)        BAS
2   0.8153929 0.02165589
3   0.8129093 0.02170883

Residual Deviance: 512.0577 
AIC: 520.0577 

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

# weights:  9 (4 variable)
initial  value 260.371112 
final  value 256.028871 
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)    6.32     0.815       2.26 0.0238     1.28     31.2  
2 2       BAS            0.951    0.0217     -2.31 0.0208     0.912     0.992
3 3       (Intercept)    8.40     0.813       2.62 0.00884    1.71     41.3  
4 3       BAS            0.944    0.0217     -2.68 0.00740    0.904     0.985

## Uncontrolled eating on sugar intake 
full.df %>%
  ggplot(aes(uncontrolled, SugarCat)) +
  geom_point() +
  geom_smooth(method = lm) +
  stat_cor() +
  theme_bw()

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

summary(multinom(SugarCat ~ uncontrolled, data = full.df)) # AIC = 522.7

# weights:  9 (4 variable)
initial  value 260.371112 
final  value 257.373195 
converged

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

Coefficients:
  (Intercept) uncontrolled
2  -0.4450276   0.02421409
3  -1.5638772   0.08231878

Std. Errors:
  (Intercept) uncontrolled
2   0.6601550   0.03486078
3   0.6871609   0.03518088

Residual Deviance: 514.7464 
AIC: 522.7464 

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

# weights:  9 (4 variable)
initial  value 260.371112 
final  value 257.373195 
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.641    0.660     -0.674  0.500    0.176      2.34 
2 2       uncontrolled    1.02     0.0349     0.695  0.487    0.957      1.10 
3 3       (Intercept)     0.209    0.687     -2.28   0.0229   0.0544     0.805
4 3       uncontrolled    1.09     0.0352     2.34   0.0193   1.01       1.16 

## Health fat liking on sugar intake 
full.df %>%
  ggplot(aes(healthyfat.liking, SugarCat)) +
  geom_point() +
  geom_smooth(method = lm) +
  stat_cor() +
  theme_bw()

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

summary(multinom(SugarCat ~ healthyfat.liking, data = full.df)) # AIC = 521.5

# weights:  9 (4 variable)
initial  value 260.371112 
final  value 256.757979 
converged

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

Coefficients:
  (Intercept) healthyfat.liking
2  0.04896833       -0.00112373
3  0.47357955       -0.01236937

Std. Errors:
  (Intercept) healthyfat.liking
2   0.2846460       0.005415522
3   0.2611889       0.005278698

Residual Deviance: 513.516 
AIC: 521.516 

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

# weights:  9 (4 variable)
initial  value 260.371112 
final  value 256.757979 
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.05    0.285       0.172  0.863     0.601     1.83 
2 2       healthyfat.li…    0.999   0.00542    -0.208  0.836     0.988     1.01 
3 3       (Intercept)       1.61    0.261       1.81   0.0698    0.962     2.68 
4 3       healthyfat.li…    0.988   0.00528    -2.34   0.0191    0.978     0.998

## HINC on sugar intake 
full.df %>%
  ggplot(aes(hinc, SugarCat)) +
  geom_point() +
  geom_smooth(method = lm) +
  stat_cor() +
  theme_bw()

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

summary(multinom(SugarCat ~ hinc, data = full.df)) # AIC = 517.3

# weights:  9 (4 variable)
initial  value 260.371112 
final  value 255.109343 
converged

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

Coefficients:
  (Intercept)        hinc
2   0.5520953 -0.09407147
3   1.1960762 -0.21701847

Std. Errors:
  (Intercept)       hinc
2   0.4106713 0.06452836
3   0.4119965 0.06916000

Residual Deviance: 510.2187 
AIC: 518.2187 

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

# weights:  9 (4 variable)
initial  value 260.371112 
final  value 255.109343 
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.74     0.411       1.34 0.179      0.777     3.88 
2 2       hinc           0.910    0.0645     -1.46 0.145      0.802     1.03 
3 3       (Intercept)    3.31     0.412       2.90 0.00369    1.47      7.42 
4 3       hinc           0.805    0.0692     -3.14 0.00170    0.703     0.922

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

# weights:  45 (28 variable)
initial  value 641.807898 
iter  10 value 563.440008
iter  20 value 484.393656
iter  30 value 477.496400
final  value 477.486345 
converged

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

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

Coefficients:
     (Intercept)   sfbl.liking
260    1.2703411 -0.0026185863
286   -0.8140608  0.0001620344
416    0.3014165  0.0297920619
442   -2.6429984  0.0284696381
520   -0.8608318  0.0112613689
572   -1.3300270  0.0462566254
676   -2.1686615  0.0524088793
728   -2.9475688  0.0541118329
832   -2.2005160  0.0530442902
884   -7.7993580  0.1224271722
988   -3.0097783  0.0553317279
1144  -8.0006144  0.1156272878
1196  -4.6407414  0.0698065682
1300  -5.8963461  0.0841240406

Std. Errors:
     (Intercept) sfbl.liking
260    0.3271742  0.00954765
286    0.5237114  0.01519173
416    0.3967368  0.01056109
442    1.1580028  0.02635824
520    0.5513528  0.01489620
572    0.5987281  0.01356242
676    0.7822200  0.01627566
728    1.0706596  0.02098414
832    0.7873970  0.01632286
884    2.4710296  0.03488778
988    1.0844796  0.02109944
1144   3.2272982  0.04538504
1196   1.8530094  0.03198911
1300   2.5460740  0.04050312

Residual Deviance: 954.9727 
AIC: 1010.973 

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

# weights:  60 (42 variable)
initial  value 641.807898 
iter  10 value 499.354722
iter  20 value 495.982750
iter  30 value 483.718057
iter  40 value 468.384885
iter  50 value 468.188931
final  value 468.188724 
converged

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

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

Coefficients:
     (Intercept)   sfbl.liking         BAS
260   -0.5349149 -1.969469e-03  0.04776625
286   -0.2059971  3.426398e-05 -0.01694688
416    0.9762416  2.982751e-02 -0.01897358
442   -7.7330117  2.815185e-02  0.12986594
520    0.2302823  1.116759e-02 -0.03083735
572   -0.5228702  4.650097e-02 -0.02301963
676   -0.7266515  5.343035e-02 -0.04257131
728   -2.4367329  5.414731e-02 -0.01432395
832   -2.8783710  5.184044e-02  0.01994355
884   -4.8580693  1.475268e-01 -0.13769646
988   -3.7550860  5.392076e-02  0.02202203
1144  -5.4302935  1.326357e-01 -0.10781052
1196  -2.5410767  7.375662e-02 -0.06709997
1300  -2.4810598  9.996728e-02 -0.13296571

Std. Errors:
     (Intercept) sfbl.liking        BAS
260     1.392822 0.009649426 0.03649601
286     2.065900 0.015280763 0.05597435
416     1.376403 0.010598320 0.03660794
442     4.447734 0.024871716 0.10359800
520     1.911512 0.014971745 0.05200501
572     1.621407 0.013644597 0.04264543
676     1.820968 0.016622691 0.04798306
728     2.334600 0.021083744 0.06012617
832     1.964884 0.016096881 0.04999879
884     3.034583 0.043698005 0.08453502
988     2.483592 0.020712479 0.06252450
1144    3.954384 0.053462909 0.10575001
1196    3.074947 0.034338117 0.07991625
1300    3.567728 0.050041553 0.09717506

Residual Deviance: 936.3774 
AIC: 1020.377 

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

# weights:  60 (42 variable)
initial  value 641.807898 
iter  10 value 528.815810
iter  20 value 521.935654
iter  30 value 486.433743
iter  40 value 473.788502
iter  50 value 473.379334
iter  60 value 473.372425
final  value 473.372405 
converged

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

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

Coefficients:
      (Intercept)   sfbl.liking uncontrolled
260    0.87141031 -0.0030904468   0.02266057
286   -1.29598550 -0.0003980417   0.02727730
416   -0.09924768  0.0293077551   0.02278754
442   -1.89451515  0.0295883475  -0.04458120
520   -2.85500504  0.0095698845   0.10688081
572   -3.00960658  0.0446313012   0.09118302
676   -3.72661917  0.0508698712   0.08492679
728   -5.29643324  0.0522620289   0.12433970
832   -2.58158314  0.0525027698   0.02187864
884   -7.53861948  0.1256418545  -0.02583825
988   -3.92018940  0.0542560593   0.05094846
1144 -14.36563506  0.1141316260   0.30316624
1196  -3.75103140  0.0710156477  -0.05289084
1300  -5.21618751  0.0853838110  -0.04138633

Std. Errors:
     (Intercept) sfbl.liking uncontrolled
260     1.060577 0.009628152   0.05780639
286     1.696496 0.015310409   0.09071443
416     1.118849 0.010636090   0.05977379
442     2.647460 0.026643848   0.14456737
520     1.697140 0.014997688   0.08410870
572     1.436439 0.013640822   0.07056594
676     1.707787 0.016350679   0.08117123
728     2.279354 0.021112372   0.10123576
832     1.647765 0.016392686   0.08201110
884     3.345673 0.036946250   0.14316154
988     2.159814 0.021165520   0.10208160
1144    6.247567 0.043972251   0.20933508
1196    3.077166 0.032365139   0.15071797
1300    3.926170 0.041434765   0.18163690

Residual Deviance: 946.7448 
AIC: 1030.745 

sugar.logit.M4 <- multinom(sugar.intake ~ sfbl.liking + healthyfat.liking, data = full.df)

# weights:  60 (42 variable)
initial  value 641.807898 
iter  10 value 525.989496
iter  20 value 518.078922
iter  30 value 489.909419
iter  40 value 470.309538
iter  50 value 470.177758
final  value 470.177656 
converged

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

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

Coefficients:
     (Intercept)  sfbl.liking healthyfat.liking
260    1.0942885 -0.002024747       0.003812436
286   -0.5110087 -0.001111860      -0.007187758
416    0.1896072  0.030021610       0.002564354
442   -2.9716129  0.028538683       0.007589368
520   -1.5571192  0.012239160       0.014505465
572   -1.1268409  0.047016724      -0.006201021
676   -1.6575720  0.054006978      -0.019279138
728   -2.8555940  0.054733991      -0.002944576
832   -2.0316043  0.053857974      -0.005298603
884   -8.9196478  0.126685615       0.016669823
988   -2.7235967  0.056439235      -0.009241731
1144  -7.4129638  0.114943166      -0.015664290
1196  -4.1324525  0.073253526      -0.027542048
1300  -5.9862167  0.084461952       0.001860922

Std. Errors:
     (Intercept) sfbl.liking healthyfat.liking
260    0.5479983 0.009720194       0.009668644
286    0.7965548 0.015805523       0.014403057
416    0.5807171 0.010667831       0.009678999
442    1.5265809 0.026030087       0.022563916
520    0.9206599 0.014621559       0.014839229
572    0.7404879 0.013775547       0.010784249
676    0.8951740 0.016920198       0.011347436
728    1.2083663 0.021142222       0.014853660
832    0.9166633 0.016515225       0.012078752
884    3.0436069 0.038321337       0.020770214
988    1.1979914 0.021379455       0.014194933
1144   3.1679664 0.043826999       0.022469827
1196   1.9140931 0.032748757       0.016992981
1300   2.7343526 0.040823156       0.024515692

Residual Deviance: 940.3553 
AIC: 1024.355 

sugar.logit.M5 <- multinom(sugar.intake ~ sfbl.liking + BMI, data = full.df)

# weights:  60 (42 variable)
initial  value 641.807898 
iter  10 value 527.446666
iter  20 value 523.364217
iter  30 value 499.896121
iter  40 value 472.189422
iter  50 value 472.002983
final  value 472.000723 
converged

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

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

Coefficients:
     (Intercept)   sfbl.liking         BMI
260     1.214756 -0.0026654114 0.002396548
286    -1.362085 -0.0004710984 0.023513928
416    -1.061468  0.0279850302 0.057469651
442    -4.053187  0.0265873889 0.059406564
520    -2.137859  0.0095731858 0.053979619
572    -3.248231  0.0435081873 0.079646399
676    -3.968040  0.0498696226 0.074981835
728    -5.061983  0.0510840007 0.087102593
832    -3.579561  0.0511271617 0.058272067
884    -8.086220  0.1255440641 0.003086398
988    -4.113013  0.0537874764 0.047047493
1144   -8.636989  0.1164185537 0.023118398
1196   -5.502547  0.0685564997 0.037095178
1300  -11.578854  0.0893789574 0.175985454

Std. Errors:
     (Intercept) sfbl.liking        BMI
260     1.322033  0.00962239 0.05548999
286     2.007110  0.01536698 0.08286043
416     1.309799  0.01064918 0.05336940
442     2.585027  0.02658546 0.09546611
520     1.744078  0.01510077 0.06973146
572     1.472855  0.01369076 0.05647949
676     1.649876  0.01637787 0.06045093
728     1.950747  0.02117203 0.06606113
832     1.681561  0.01638227 0.06243285
884     3.270974  0.03731201 0.10509563
988     2.085222  0.02115960 0.07481730
1144    4.207668  0.04762817 0.12967630
1196    3.036905  0.03220074 0.10231474
1300    4.396216  0.04772595 0.07499999

Residual Deviance: 944.0014 
AIC: 1028.001 

sugar.logit.M6 <- multinom(sugar.intake ~ sfbl.liking*BMI, data = full.df)

# weights:  75 (56 variable)
initial  value 641.807898 
iter  10 value 584.882554
iter  20 value 502.913820
iter  30 value 494.499175
iter  40 value 476.631208
iter  50 value 469.837925
iter  60 value 468.287957
iter  70 value 468.009095
iter  80 value 467.971699
iter  90 value 467.966239
iter 100 value 467.965757
final  value 467.965757 
stopped after 100 iterations

summary(sugar.logit.M6) # AIC = 1047.9

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

Coefficients:
     (Intercept)   sfbl.liking          BMI sfbl.liking:BMI
260    1.4197943 -0.0106011114 -0.006782621    3.470222e-04
286   -2.0103449  0.0238983453  0.051572223   -1.042896e-03
416   -0.2184915  0.0064972637  0.022195842    8.914539e-04
442   -4.3502342  0.0305505081  0.070316432   -1.252984e-04
520   -2.3214981  0.0125189359  0.061174441   -1.052021e-04
572   -3.5022214  0.0456237030  0.087351255   -2.022067e-05
676   -6.6384216  0.0995047149  0.172305959   -1.785180e-03
728   -1.0464423 -0.0304493840 -0.071065001    3.154683e-03
832   -3.3272458  0.0431302661  0.046092685    3.731517e-04
884   -2.6091173  0.0466121981 -0.228199004    3.323640e-03
988  -10.4040260  0.1750098203  0.281447914   -4.565477e-03
1144  12.1421573 -0.1582627536 -0.882138662    1.164039e-02
1196  -8.8164496  0.1251269086  0.160977781   -2.091850e-03
1300  -6.9926872 -0.0007714451  0.021869991    3.038531e-03

Std. Errors:
     (Intercept) sfbl.liking        BMI sfbl.liking:BMI
260  0.001190816  0.03787984 0.01423315     0.001608115
286  0.001629269  0.06070251 0.02232090     0.002574154
416  0.001664550  0.03550012 0.01669499     0.001498940
442  0.002121218  0.05851214 0.04610817     0.002403493
520  0.001386858  0.04639827 0.02296971     0.001907707
572  0.001344261  0.03684875 0.02380607     0.001566596
676  0.001378717  0.03850758 0.02979279     0.001669334
728  0.001783572  0.04211377 0.04021676     0.001744448
832  0.001473960  0.03894281 0.03116313     0.001664592
884  0.004244299  0.04935631 0.11018453     0.002103658
988  0.001640391  0.04449629 0.03847516     0.002023861
1144 0.006873898  0.06674157 0.14592958     0.002150049
1196 0.002755112  0.05264266 0.07064687     0.002474156
1300 0.003226464  0.05687144 0.09321242     0.002675494

Residual Deviance: 935.9315 
AIC: 1047.932 

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

# weights:  75 (56 variable)
initial  value 641.807898 
iter  10 value 605.037958
iter  20 value 541.546397
iter  30 value 508.190499
iter  40 value 485.904910
iter  50 value 470.478653
iter  60 value 468.330185
iter  70 value 467.649904
iter  80 value 467.428814
iter  90 value 467.341435
iter 100 value 467.297583
final  value 467.297583 
stopped after 100 iterations

summary(sugar.logit.M7) # AIC = 1046.6

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

Coefficients:
     (Intercept)  sfbl.liking uncontrolled sfbl.liking:uncontrolled
260    1.6381929 -0.052070578  -0.02016199             0.0027017333
286   -1.4246360  0.004484248   0.03258245            -0.0001853628
416    0.0472990  0.006749204   0.01214353             0.0013459894
442    1.6886407 -0.085046618  -0.27748672             0.0069406190
520   -1.3376743 -0.067527794   0.02453757             0.0040867230
572   -0.7754152 -0.028690003  -0.02849630             0.0039534860
676   -1.8492764 -0.012045192  -0.01685539             0.0034442096
728  -11.9570955  0.163958012   0.41439532            -0.0047196416
832   -0.9094421 -0.003788809  -0.07358490             0.0031963236
884  -19.0553538  0.268038180   0.51181910            -0.0063758147
988   -9.1325885  0.140480547   0.29861118            -0.0039908839
1144 -24.2205104  0.246099554   0.70444400            -0.0050635400
1196  -5.8572494  0.088775441   0.05641329            -0.0007836317
1300 -10.9462085  0.161110208   0.25287574            -0.0037647753

Std. Errors:
     (Intercept) sfbl.liking uncontrolled sfbl.liking:uncontrolled
260  0.798061000  0.03656844   0.04461023              0.001951590
286  0.245370711  0.05424156   0.02923204              0.002891617
416  1.126630031  0.04029724   0.06083240              0.002154405
442  0.053134441  0.05966946   0.07124675              0.003279536
520  0.431769393  0.04832154   0.03663153              0.002482252
572  0.199951132  0.03735250   0.03301716              0.002016224
676  0.069427572  0.04032934   0.04061014              0.002175777
728  0.021647497  0.04436895   0.04194689              0.002481989
832  0.077124807  0.03907552   0.04312978              0.002152009
884  0.005743853  0.04628569   0.12238336              0.002977082
988  0.029788188  0.04360758   0.04578449              0.002527306
1144 0.005115033  0.07636551   0.11793474              0.004226457
1196 0.010333526  0.05356556   0.10268555              0.003394853
1300 0.007453464  0.05913142   0.12885879              0.003984153

Residual Deviance: 934.5952 
AIC: 1046.595 

sugar.logit.M8 <- multinom(sugar.intake ~ sfbl.liking*hinc, data = full.df)

# weights:  75 (56 variable)
initial  value 641.807898 
iter  10 value 572.490019
iter  20 value 526.318429
iter  30 value 509.203209
iter  40 value 484.450328
iter  50 value 462.877044
iter  60 value 461.822705
iter  70 value 461.779129
iter  80 value 461.777743
final  value 461.777702 
converged

summary(sugar.logit.M8) # AIC = 1040.1

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

Coefficients:
     (Intercept)  sfbl.liking        hinc sfbl.liking:hinc
260    1.0990438 -0.016648438  0.03734950     0.0022801225
286   -2.0718125 -0.035129446  0.19574636     0.0054824106
416    0.9946835  0.012090046 -0.11361915     0.0028982243
442   -2.4599184 -0.028729226 -0.02670267     0.0091726485
520   -0.3762020 -0.005449790 -0.07127568     0.0026033869
572    0.1613359  0.022276869 -0.27760521     0.0042579254
676   -0.3621219  0.024106562 -0.34798307     0.0052418619
728   -0.9259092  0.049756525 -0.36814487    -0.0013062241
832   -0.4784571  0.007044848 -0.32009158     0.0084150808
884   -7.1211505  0.067174528 -0.16135171     0.0099000512
988   -3.9503477  0.078461138  0.18109608    -0.0049583489
1144   0.9764943 -0.037564559 -2.07430502     0.0328900669
1196  -1.7292543  0.060466689 -0.58007841    -0.0008791537
1300  -2.8356552  0.030930756 -0.61641793     0.0104552202

Std. Errors:
     (Intercept) sfbl.liking      hinc sfbl.liking:hinc
260    0.8096938  0.02240082 0.1102213      0.002961412
286    1.5187653  0.04559008 0.1903608      0.005605942
416    0.9393192  0.02442240 0.1441232      0.003780269
442    2.8624972  0.07337016 0.4265107      0.010564792
520    1.2974729  0.03544547 0.1956604      0.005465269
572    1.3828561  0.03099423 0.2557737      0.005518922
676    1.8053225  0.03720347 0.3570243      0.007028789
728    2.4175829  0.04573132 0.5439856      0.010422512
832    1.7852345  0.03790599 0.3259044      0.006555896
884    0.1614402  0.03124823 0.3955217      0.007355112
988    2.4739361  0.04438950 0.3391275      0.006524933
1144   3.3247532  0.05519251 1.1593218      0.014687988
1196   0.6506522  0.02920455 0.5260837      0.011536637
1300   0.2363936  0.03473632 0.5068826      0.009169387

Residual Deviance: 923.5554 
AIC: 1035.555 

sugar.logit.M9 <- multinom(sugar.intake ~ sfbl.liking*BMI + sfbl.liking*uncontrolled, data = full.df)

# weights:  105 (84 variable)
initial  value 641.807898 
iter  10 value 580.367033
iter  20 value 535.163630
iter  30 value 505.365644
iter  40 value 492.183660
iter  50 value 485.112377
iter  60 value 473.090310
iter  70 value 466.217070
iter  80 value 459.911488
iter  90 value 458.110448
iter 100 value 457.780060
final  value 457.780060 
stopped after 100 iterations

summary(sugar.logit.M9) # AIC = 1083.6

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

Coefficients:
     (Intercept) sfbl.liking         BMI uncontrolled sfbl.liking:BMI
260    1.9206705 -0.06470216 -0.01177286 -0.021272718    4.782972e-04
286   -2.8734182  0.03319133  0.05484089  0.041249290   -1.082400e-03
416   -0.3417529 -0.01977877  0.01774690  0.011360492    1.055760e-03
442   -0.9377181 -0.06904547  0.12615102 -0.303175378   -1.223612e-03
520   -2.8304987 -0.06282068  0.06461605  0.022358621   -4.148299e-04
572   -2.7027332 -0.02859587  0.08406938 -0.037695323   -4.055282e-05
676   -5.8846116  0.03043397  0.16541158 -0.029892442   -1.687652e-03
728  -10.4282691  0.09298475 -0.01326908  0.367479586    2.074279e-03
832   -2.1779567 -0.01192888  0.05080785 -0.072706341    3.387852e-04
884   -4.2453726  0.06962597 -0.43899368  0.308420490    5.968065e-03
988  -20.5873673  0.34098575  0.35172313  0.410632609   -6.058002e-03
1144  -0.2869188 -0.08524066 -0.71910227  0.407877084    9.600238e-03
1196  -8.1400814  0.11400586  0.13169501  0.001697476   -1.595813e-03
1300  -5.1828854 -0.03293880 -0.05418907  0.050341682    4.626531e-03
     sfbl.liking:uncontrolled
260              2.798563e-03
286             -3.742415e-04
416              1.299117e-03
442              7.656522e-03
520              4.306158e-03
572              3.983574e-03
676              3.565269e-03
728             -4.098748e-03
832              3.162534e-03
884             -3.816069e-03
988             -6.184832e-03
1144            -8.506333e-04
1196             4.744253e-05
1300            -1.231335e-03

Std. Errors:
     (Intercept) sfbl.liking        BMI uncontrolled sfbl.liking:BMI
260  0.003576089  0.05442541 0.04549285   0.05523493     0.001924991
286  0.001980610  0.08738033 0.07073865   0.08512026     0.003049082
416  0.005385615  0.05211414 0.05228796   0.06484124     0.001910523
442  0.002082266  0.08299192 0.14213583   0.24562830     0.003537347
520  0.001579228  0.06648137 0.07073215   0.09058341     0.002396308
572  0.003335451  0.05509728 0.07523272   0.10180596     0.002186528
676  0.002256959  0.05791125 0.09288199   0.13163425     0.002425599
728  0.001652140  0.06284011 0.11430009   0.11467539     0.002650458
832  0.002560709  0.05768977 0.09873994   0.13801168     0.002489358
884  0.003617475  0.07436859 0.38516559   0.40561175     0.005155482
988  0.002126191  0.06693356 0.11726095   0.14976342     0.003032674
1144 0.002835376  0.12594327 0.48676203   0.37743479     0.005488477
1196 0.004520557  0.07996025 0.22232150   0.32887931     0.004423314
1300 0.008460877  0.08420497 0.23878428   0.38403203     0.004809092
     sfbl.liking:uncontrolled
260               0.002083202
286               0.003338462
416               0.002202329
442               0.004982300
520               0.002952171
572               0.002762845
676               0.003197004
728               0.002911129
832               0.003259302
884               0.006561060
988               0.003528439
1144              0.007703255
1196              0.006157724
1300              0.006557072

Residual Deviance: 915.5601 
AIC: 1083.56 

Linear regression - sHEI

## Sweet liking on sHEI
full.df %>%
  ggplot(aes(sfbl.liking, sHEI)) +
  geom_point() +
  geom_smooth(method = lm) +
  stat_cor() +
  theme_bw()

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

summary(lm(sHEI ~ sfbl.liking, data = full.df)) # r2 = 0.02

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

Residuals:
     Min       1Q   Median       3Q      Max 
-23.7123  -6.4696   0.3182   5.7154  21.3145 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 48.63972    0.89438   54.38  < 2e-16 ***
sfbl.liking -0.05130    0.01973   -2.60  0.00991 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 8.65 on 235 degrees of freedom
Multiple R-squared:  0.02796,   Adjusted R-squared:  0.02383 
F-statistic: 6.761 on 1 and 235 DF,  p-value: 0.00991

tidy(lm(sHEI ~ sfbl.liking, data = 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.33e+21    0.894      54.4  3.95e-135 2.28e+20  7.75e+21
2 sfbl.liking 9.50e- 1    0.0197     -2.60 9.91e-  3 9.14e- 1  9.88e- 1

## Healthy fat liking on sHEI
full.df %>%
  ggplot(aes(healthyfat.liking, sHEI)) +
  geom_point() +
  geom_smooth(method = lm) +
  stat_cor() +
  theme_bw()

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

summary(lm(sHEI ~ healthyfat.liking, data = full.df)) # r2 = 0.07

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

Residuals:
     Min       1Q   Median       3Q      Max 
-23.7814  -5.6514   0.4074   5.5169  21.7526 

Coefficients:
                  Estimate Std. Error t value Pr(>|t|)    
(Intercept)       43.84820    0.88363  49.623  < 2e-16 ***
healthyfat.liking  0.07514    0.01745   4.305 2.45e-05 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 8.446 on 235 degrees of freedom
Multiple R-squared:  0.07311,   Adjusted R-squared:  0.06917 
F-statistic: 18.54 on 1 and 235 DF,  p-value: 2.449e-05

tidy(lm(sHEI ~ healthyfat.liking, data = 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.10e19    0.884      49.6  1.58e-126  1.94e18   6.30e19
2 healthyfat.liking  1.08e 0    0.0175      4.31 2.45e-  5  1.04e 0   1.12e 0

## Restraint on sHEI
full.df %>%
  ggplot(aes(restraint, sHEI)) +
  geom_point() +
  geom_smooth(method = lm) +
  stat_cor() +
  theme_bw()

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

summary(lm(sHEI ~ restraint, data = full.df)) # r2 = 0.05

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

Residuals:
     Min       1Q   Median       3Q      Max 
-25.0333  -5.5541  -0.0946   6.0067  23.9654 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  39.0860     2.2843  17.110  < 2e-16 ***
restraint     0.5499     0.1573   3.495 0.000566 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 8.554 on 235 degrees of freedom
Multiple R-squared:  0.04941,   Adjusted R-squared:  0.04537 
F-statistic: 12.22 on 1 and 235 DF,  p-value: 0.0005662

tidy(lm(sHEI ~ restraint, data = 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)  9.44e16     2.28      17.1  3.60e-43  1.05e15   8.50e18
2 restraint    1.73e 0     0.157      3.50 5.66e- 4  1.27e 0   2.36e 0

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

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

Residuals:
     Min       1Q   Median       3Q      Max 
-23.7814  -5.6514   0.4074   5.5169  21.7526 

Coefficients:
                  Estimate Std. Error t value Pr(>|t|)    
(Intercept)       43.84820    0.88363  49.623  < 2e-16 ***
healthyfat.liking  0.07514    0.01745   4.305 2.45e-05 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 8.446 on 235 degrees of freedom
Multiple R-squared:  0.07311,   Adjusted R-squared:  0.06917 
F-statistic: 18.54 on 1 and 235 DF,  p-value: 2.449e-05

sHEI.lm.M2 <- lm(sHEI ~ healthyfat.liking + sfbl.liking , data = full.df)
summary(sHEI.lm.M2) # r2 = 0.08

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

Residuals:
     Min       1Q   Median       3Q      Max 
-24.8598  -5.2330   0.3897   5.5929  21.4113 

Coefficients:
                  Estimate Std. Error t value Pr(>|t|)    
(Intercept)       45.52834    1.15387  39.457  < 2e-16 ***
healthyfat.liking  0.07101    0.01740   4.080 6.17e-05 ***
sfbl.liking       -0.04299    0.01921  -2.238   0.0262 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 8.375 on 234 degrees of freedom
Multiple R-squared:  0.09253,   Adjusted R-squared:  0.08478 
F-statistic: 11.93 on 2 and 234 DF,  p-value: 1.165e-05

sHEI.lm.M3 <- lm(sHEI ~ healthyfat.liking + restraint, data = full.df)
summary(sHEI.lm.M3) # r2 = 0.11

Call:
lm(formula = sHEI ~ healthyfat.liking + restraint, data = full.df)

Residuals:
     Min       1Q   Median       3Q      Max 
-26.2066  -5.1040   0.2085   5.8150  23.8263 

Coefficients:
                  Estimate Std. Error t value Pr(>|t|)    
(Intercept)       36.84058    2.27341  16.205  < 2e-16 ***
healthyfat.liking  0.07139    0.01713   4.169 4.32e-05 ***
restraint          0.50811    0.15244   3.333 0.000998 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 8.27 on 234 degrees of freedom
Multiple R-squared:  0.1151,    Adjusted R-squared:  0.1076 
F-statistic: 15.22 on 2 and 234 DF,  p-value: 6.097e-07

sHEI.lm.M4 <- lm(sHEI ~ healthyfat.liking*sfbl.liking , data = full.df)
summary(sHEI.lm.M4) # r2 = 0.08

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

Residuals:
     Min       1Q   Median       3Q      Max 
-24.8654  -5.2494   0.3936   5.5869  21.4146 

Coefficients:
                                Estimate Std. Error t value Pr(>|t|)    
(Intercept)                    4.549e+01  1.693e+00  26.861   <2e-16 ***
healthyfat.liking              7.180e-02  3.046e-02   2.357   0.0192 *  
sfbl.liking                   -4.204e-02  3.575e-02  -1.176   0.2408    
healthyfat.liking:sfbl.liking -1.969e-05  6.256e-04  -0.031   0.9749    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 8.393 on 233 degrees of freedom
Multiple R-squared:  0.09254,   Adjusted R-squared:  0.08085 
F-statistic:  7.92 on 3 and 233 DF,  p-value: 4.725e-05

# Model comparison
anova(sHEI.lm.M3, sHEI.lm.M2)

Analysis of Variance Table

Model 1: sHEI ~ healthyfat.liking + restraint
Model 2: sHEI ~ healthyfat.liking + sfbl.liking
  Res.Df   RSS Df Sum of Sq F Pr(>F)
1    234 16006                      
2    234 16414  0   -408.69         

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

Analysis of Variance Table

Model 1: sHEI ~ healthyfat.liking + restraint
Model 2: sHEI ~ healthyfat.liking * sfbl.liking
  Res.Df   RSS Df Sum of Sq F Pr(>F)
1    234 16006                      
2    233 16414  1   -408.62         

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

Analysis of Variance Table

Model 1: sHEI ~ healthyfat.liking + restraint
Model 2: sHEI ~ healthyfat.liking
  Res.Df   RSS Df Sum of Sq      F    Pr(>F)    
1    234 16006                                  
2    235 16766 -1   -759.96 11.111 0.0009977 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Logistic regression - sHEI

## Sweet liking on sHEI
full.df %>%
  ggplot(aes(sfbl.liking, DQ)) +
  geom_point() +
  geom_smooth(method = lm) +
  stat_cor() +
  theme_bw()

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

summary(multinom(DQ ~ sfbl.liking, data = full.df)) # AIC = 522.4

# weights:  9 (4 variable)
initial  value 260.371112 
final  value 257.198432 
converged

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

Coefficients:
  (Intercept)  sfbl.liking
2   0.1596656 -0.004151015
3   0.4796598 -0.013838294

Std. Errors:
  (Intercept) sfbl.liking
2   0.2768875 0.005881625
3   0.2601893 0.005805992

Residual Deviance: 514.3969 
AIC: 522.3969 

tidy(multinom(DQ ~ sfbl.liking, data = full.df), exponentiate = TRUE, conf.int = TRUE)

# weights:  9 (4 variable)
initial  value 260.371112 
final  value 257.198432 
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.17    0.277       0.577  0.564     0.682     2.02 
2 2       sfbl.liking    0.996   0.00588    -0.706  0.480     0.984     1.01 
3 3       (Intercept)    1.62    0.260       1.84   0.0653    0.970     2.69 
4 3       sfbl.liking    0.986   0.00581    -2.38   0.0172    0.975     0.998

## Healthy fat liking on sHEI
full.df %>%
  ggplot(aes(healthyfat.liking, DQ)) +
  geom_point() +
  geom_smooth(method = lm) +
  stat_cor() +
  theme_bw()

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

summary(multinom(DQ ~ healthyfat.liking, data = full.df)) # AIC = 510.5

# weights:  9 (4 variable)
initial  value 260.371112 
final  value 251.225077 
converged

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

Coefficients:
  (Intercept) healthyfat.liking
2  -0.2255535       0.006600251
3  -0.9692820       0.023322771

Std. Errors:
  (Intercept) healthyfat.liking
2   0.2364089       0.005066117
3   0.3035517       0.005926746

Residual Deviance: 502.4502 
AIC: 510.4502 

tidy(multinom(DQ ~ healthyfat.liking, data = full.df), exponentiate = TRUE, conf.int = TRUE)

# weights:  9 (4 variable)
initial  value 260.371112 
final  value 251.225077 
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.798   0.236      -0.954 3.40e-1    0.502     1.27 
2 2       healthyfat.li…    1.01    0.00507     1.30  1.93e-1    0.997     1.02 
3 3       (Intercept)       0.379   0.304      -3.19  1.41e-3    0.209     0.688
4 3       healthyfat.li…    1.02    0.00593     3.94  8.31e-5    1.01      1.04 

## Restraint on sHEI
full.df %>%
  ggplot(aes(restraint, DQ)) +
  geom_point() +
  geom_smooth(method = lm) +
  stat_cor() +
  theme_bw()

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

summary(multinom(DQ ~ restraint, data = full.df)) # AIC = 515.7

# weights:  9 (4 variable)
initial  value 260.371112 
final  value 253.825525 
converged

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

Coefficients:
  (Intercept) restraint
2   -1.553936 0.1137632
3   -2.297396 0.1645437

Std. Errors:
  (Intercept)  restraint
2   0.6614349 0.04693810
3   0.6938453 0.04821722

Residual Deviance: 507.6511 
AIC: 515.6511 

tidy(multinom(DQ ~ restraint, data = full.df), exponentiate = TRUE, conf.int = TRUE)

# weights:  9 (4 variable)
initial  value 260.371112 
final  value 253.825525 
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.211    0.661      -2.35 0.0188     0.0578     0.773
2 2       restraint      1.12     0.0469      2.42 0.0154     1.02       1.23 
3 3       (Intercept)    0.101    0.694      -3.31 0.000929   0.0258     0.392
4 3       restraint      1.18     0.0482      3.41 0.000644   1.07       1.30 

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

# weights:  9 (4 variable)
initial  value 260.371112 
final  value 251.225077 
converged

summary(DQ.logit.M1) # AIC = 510.5

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

Coefficients:
  (Intercept) healthyfat.liking
2  -0.2255535       0.006600251
3  -0.9692820       0.023322771

Std. Errors:
  (Intercept) healthyfat.liking
2   0.2364089       0.005066117
3   0.3035517       0.005926746

Residual Deviance: 502.4502 
AIC: 510.4502 

DQ.logit.M2 <- multinom(DQ ~ healthyfat.liking + sfbl.liking, data = full.df)

# weights:  12 (6 variable)
initial  value 260.371112 
iter  10 value 248.844932
iter  10 value 248.844932
iter  10 value 248.844932
final  value 248.844932 
converged

summary(DQ.logit.M2) # AIC = 509.7

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

Coefficients:
  (Intercept) healthyfat.liking  sfbl.liking
2 -0.06994517       0.006333468 -0.003809689
3 -0.51541215       0.022774803 -0.012458403

Std. Errors:
  (Intercept) healthyfat.liking sfbl.liking
2   0.3394437       0.005057540 0.006043198
3   0.3769700       0.006029911 0.006047557

Residual Deviance: 497.6899 
AIC: 509.6899 

DQ.logit.M3 <- multinom(DQ ~ healthyfat.liking + restraint, data = full.df)

# weights:  12 (6 variable)
initial  value 260.371112 
iter  10 value 245.177096
final  value 245.176433 
converged

summary(DQ.logit.M3) # AIC = 502.4

Call:
multinom(formula = DQ ~ healthyfat.liking + restraint, data = full.df)

Coefficients:
  (Intercept) healthyfat.liking restraint
2   -1.750571       0.006353406 0.1122357
3   -3.248097       0.023130144 0.1635274

Std. Errors:
  (Intercept) healthyfat.liking  restraint
2   0.6849041       0.005129672 0.04698814
3   0.7779846       0.006037818 0.05002966

Residual Deviance: 490.3529 
AIC: 502.3529 

DQ.logit.M4 <- multinom(DQ ~ healthyfat.liking + sfbl.liking + restraint, data = full.df)

# weights:  15 (8 variable)
initial  value 260.371112 
iter  10 value 242.725222
final  value 242.562530 
converged

summary(DQ.logit.M4) # AIC = 501.1

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

Coefficients:
  (Intercept) healthyfat.liking  sfbl.liking restraint
2    -1.58984       0.005902361 -0.004112507 0.1131906
3    -2.83542       0.022440523 -0.013334665 0.1687064

Std. Errors:
  (Intercept) healthyfat.liking sfbl.liking  restraint
2   0.7220845       0.005128117 0.006125133 0.04720895
3   0.8014702       0.006147039 0.006205341 0.05039013

Residual Deviance: 485.1251 
AIC: 501.1251 

# Model comparison
anova(DQ.logit.M4, DQ.logit.M3)

Likelihood ratio tests of Multinomial Models

Response: DQ
                                        Model Resid. df Resid. Dev   Test    Df
1               healthyfat.liking + restraint       468   490.3529             
2 healthyfat.liking + sfbl.liking + restraint       466   485.1251 1 vs 2     2
  LR stat.    Pr(Chi)
1                    
2 5.227806 0.07324811