Sweet Foods Liking and Diet Quality

Load packages and call data

Functions

Useful function that will be used throughout the analyses.

## 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 groups

## 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: 193
alpha: 0.687

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

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

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: 192
alpha: 0.642

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: 136
alpha: 0.749

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: 205
alpha: 0.73

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: 203
alpha: 0.573

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: 185
alpha: 0.496

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: 203
alpha: 0.537

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: 183
alpha: 0.598

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: 209
alpha: 0.799

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: 186
alpha: 0.614

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: 185
alpha: 0.808

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: 165
alpha: 0.469

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: 199
alpha: 0.849

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

Decoding survey

## 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), 
         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),
         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")) 

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

## sHEI score
HEI.df <- survey.df %>%
  select("ID", "sex", "Fruit":"Water") %>%
  mutate(fruit.1 = case_when(
    Fruit == 1 ~ 0,
    Fruit == 2 ~ 2, 
    Fruit == 3 ~ 3.5,
    Fruit >= 4 ~ 5)) %>%
  mutate(fruit.2 = case_when (
    Juice == 1 ~ 0,
    Juice == 2 ~ 2,
    Juice == 3 ~ 3.5,
    Juice >= 4 ~ 5)) %>% 
  mutate(tfruitHEI = case_when(
    fruit.1 + fruit.2 <5 ~ 0,
    fruit.1 + fruit.2 >= 5 ~ 5)) %>%
  mutate(wfruitHEI = case_when(
    Fruit == 1 ~ 0,
    Fruit == 2 ~ 2.5, 
    Fruit >= 3 ~ 5)) %>%
  mutate(vegHEI = case_when(
    `Green veg` == 1 ~ 1.6,
    `Green veg` == 2 & `Starchy veg` >= 2 ~ 2.46,
    `Green veg` >= 2 & `Starchy veg` >= 2 ~ 3.24,
    `Green veg` >= 2 & `Starchy veg` == 1 ~ 3.56)) %>%
  mutate(bean.1 = case_when(
    `Green veg` == 1 ~ 0,
    `Green veg` >= 2 ~5)) %>%
  mutate(bean.2 = case_when(
    Beans == 1 ~ 0,
    Beans >= 2 ~ 5)) %>%
  mutate(beanHEI = case_when(
    bean.1 + bean.2 < 5 ~ 0,
    bean.1 + bean.2 >= 5 ~ 5)) %>%
  mutate(wgHEI = case_when(
    `Whole grains.day` == 1 ~ 0.51,
    sex == 1 & `Whole grains.day` >= 2 ~ 2.97,
    sex == 2 & `Whole grains.day` >= 2 & `Whole grains.day` <= 3 ~ 5.20,
    sex == 2 & `Whole grains.day` >= 4 ~ 6.94)) %>%
  mutate(dairyHEI = case_when(
    sex == 1 & `Milk.day` <= 3 ~ 3.22,
    sex == 2 & `Milk.day` <= 3 & `Low-fat milk.day` == 1 ~ 3.32,
    sex == 2 & `Milk.day` <= 3 & `Low-fat milk.day` >= 2 ~ 4.81,
    `Milk.day` >= 4 ~ 6.51)) %>%
  mutate(tproteinHEI = case_when(
    sex == 1 & `Seafood.freq` <= 4 ~ 4.11,
    sex == 1 & `Seafood.freq` >= 5 ~ 4.98,
    sex == 1 & is.na(`Seafood.freq`) ~ 4.11,
    sex == 2 ~ 4.97)) %>%
  mutate(sfpproteinHEI = case_when(
    sex == 1 & Nuts <= 2 ~ 0.49,
    sex == 2 & Nuts <= 2 ~ 1.50,
    Nuts >= 3 ~ 4.20)) %>%
  mutate(fatHEI = case_when(
    `Milk.day` >= 4 ~ 2.56,
    `Saturated fats` >= 2 & `Saturated fats` <= 3 & `Milk.freq` >= 1 & `Low-fat milk.freq` <= 2 ~ 2.63,
    `Saturated fats` >= 2 & `Saturated fats` <= 3 & is.na(`Milk.freq`) | is.na(`Low-fat milk.freq`) ~ 2.63,
    `Saturated fats` >= 2 & `Saturated fats` <= 3 & `Milk.freq` >= 1 & `Low-fat milk.freq` >= 3 ~ 4.54,
    `Saturated fats` == 1 & `Milk.freq` >= 1 ~ 5.93,
    `Saturated fats` == 1 & is.na(`Milk.freq`) ~ 5.93)) %>%
  mutate(grainHEI = case_when(
    `Green veg` == 1 ~ 2.13,
    `Grains.day` >= 3 & `Seafood.day` >= 2 & `Green veg` >= 2 ~ 2.27,
    `Grains.day` >= 3 & Nuts >= 1 & Nuts <= 2 & `Seafood.day` == 1 & `Green veg` >= 2 ~ 4.73,
    `Grains.day` >= 3 & Nuts >= 3 & `Seafood.day` == 1 & `Green veg` >= 2 ~ 8.45,
    `Grains.day` >= 1 & `Grains.day` <= 2 & `Green veg` >= 2 ~ 9.25)) %>%
  mutate(sodiumHEI = case_when(
    Fruit >= 1 & Fruit <= 2 & `Grains.day` >= 3 & Water == 3 ~ 0.70,
    Fruit >= 3 & `Grains.day` >= 3 & Water == 3 ~ 2.30,
    `Grains.day` >= 3 & Water >= 1 & Water <= 2 ~ 4.94,
    `Grains.day` >= 1 & `Grains.day` <= 2 ~ 6.07)) %>%
  mutate(ssb.cal = case_when(
    `SSB.day` == 1  ~ 0,
    `SSB.day` == 2 ~ 156,
    `SSB.day` == 3 ~ 312,
    `SSB.day` == 4 ~ 468,
    `SSB.day` == 5 ~ 624,
    `SSB.day` == 6 ~ 780,
    `SSB.day` == 7 ~ 936),
    sugar.cal = case_when(
    `Added sugars` == 1 ~ 130,
    `Added sugars` == 2 ~ 260,
    `Added sugars` == 3 ~ 520)) %>%
  mutate(sugar.intake = (ssb.cal + sugar.cal)) %>%
  mutate(sugarHEI = case_when(
    sugar.intake <= 130 ~ 10,
    sugar.intake > 130 & sugar.cal < 520 ~ 5,
    sugar.intake >= 520 ~ 0)) %>%
  mutate(sugar.tsp = case_when(
    `SSB.freq` <= 4 ~ 13.26,
    `SSB.freq` >= 5 & `SSB.freq` <=6 ~ 16.00,
    `SSB.day` ==2 ~ 16.00,
    `SSB.day` >=3 ~ 26.87)) %>%
  mutate(satfatHEI = case_when(
    `SSB.day` >= 3 ~ 1.82,
    `SSB.day` <= 2 & `Grains.day` <= 2 ~ 3.20,
    `SSB.day` <= 2 & `Grains.day` >= 3 & Nuts <= 2 ~ 4.64,
    `SSB.day` <= 2 & `Grains.day` >= 3 & Nuts >= 3 ~ 6.56)) %>%
  select("ID", "sugar.intake", "sugar.tsp", ends_with("HEI")) %>%
  mutate(sHEI = rowSums(select(., "tfruitHEI":"satfatHEI"), na.rm = TRUE)) %>%
  mutate(SuFatHEI = rowSums(select(., c("sugarHEI", "satfatHEI")))) 

## Putting scores together

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

full.df <- survey.df  %>%
  select(c("ID", "sex", "age", "race", "usb", "education", "hinc", "exec", "ageCat")) %>%
  filter(race %in% c("African", "Asian")) %>%
  left_join(anthro.df, by = "ID") %>%
  left_join(SST.df, by = "ID") %>%
  left_join(redjade.df, by = "ID") %>%
  left_join(liking.df, by = "ID") %>%
  left_join(HEI.df, by = "ID") %>%
  left_join(TFQ.df, by = "ID") %>%
  left_join(BAS.df, by = "ID") %>%
  left_join(IES.df, by = "ID") %>%
  left_join(AEBQ.df, by = "ID") %>% 
  mutate(BMIcat = case_when(
      BMI < 18.5 ~ "0",
      BMI >= 18.5 & BMI  < 25 & race != "Asian" ~ "1",
      BMI >= 25 & BMI < 30 & race != "Asian" ~ "2",
      BMI > 30 & race != "Asian" ~ "3",
      BMI >= 18.5 & BMI  < 23 & race == "Asian" ~ "1",
      BMI >= 23 & BMI < 25 & race == "Asian" ~ "2",
      BMI > 25 & race == "Asian" ~ "3"))  %>%
  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(SugarCat = ntile(sugar.intake, 3)) %>%
  mutate(DQ = ntile(sHEI, 2)) 

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

African.df <- full.df %>%
  filter(race == "African")
  
Asian.df <- full.df %>%
  filter(race == "Asian")

USB.df <- full.df %>%
  filter(race == "U.S. Born")
  
Immigrant.df <- full.df %>%
  filter(race == "Immigrant")

Data distribution

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

Warning: The dot-dot notation (`..density..`) was deprecated in ggplot2 3.4.0.
ℹ Please use `after_stat(density)` instead.

`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

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

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

    Shapiro-Wilk normality test

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

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

`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

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

shapiro.test(full.df$BMI) ## Not normal

    Shapiro-Wilk normality test

data:  full.df$BMI
W = 0.86029, p-value = 1.285e-12

## SST
histfunc(full.df, SST, "Simple Sweet Test histogram") +
  xlim(-130, 130) +
  labs(y= "Density", x = "Sweet Taste Liking Score")

`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Warning: Removed 1 row containing non-finite outside the scale range
(`stat_bin()`).

Warning: Removed 1 row containing non-finite outside the scale range
(`stat_density()`).

qqfunc(full.df, SST, "Simple Sweet QQ plot") +
  labs(y= "Sample Quantiles", x = "Theoretical Quantiles")

Warning: Removed 1 row containing non-finite outside the scale range
(`stat_qq()`).

Warning: Removed 1 row containing non-finite outside the scale range
(`stat_qq_line()`).

shapiro.test(full.df$SST) ## Normal

    Shapiro-Wilk normality test

data:  full.df$SST
W = 0.99247, p-value = 0.3949

full.df %>%
  ggplot(aes(SST)) +
  geom_density(aes(fill= group), alpha = 0.4) + 
    xlim(-150, 180) +
    labs(title="Sweet Taste Liking by Ancestry Groups and Immigrantion Status") +
  theme_classic()

Warning: Removed 1 row containing non-finite outside the scale range
(`stat_density()`).

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

`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

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

shapiro.test(full.df$sfbl.liking) ## Not normal

    Shapiro-Wilk normality test

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

## Added sugars intake
histfunc(full.df, sugar.intake, "Added sugars intake histogram") ## transform

`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

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

shapiro.test(full.df$sugar.intake) ## Not normal

    Shapiro-Wilk normality test

data:  full.df$sugar.intake
W = 0.92699, p-value = 1.836e-08

## Total HEI
histfunc(full.df, sHEI, "Total HEI intake histogram")

`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

qqfunc(full.df, sHEI, "Overall HEI intake QQ plot")

shapiro.test(full.df$sHEI) ## Normal

    Shapiro-Wilk normality test

data:  full.df$sHEI
W = 0.99162, p-value = 0.3004

Sensory reliability - intensity & liking

## 0.09M intensity
full.df %>% rownames_to_column(var = "label") %>%
   ggplot(aes(x = lint_d1, y = lint_d2, label = label)) +
    geom_point() +
    geom_abline(intercept = 0, slope = 1, color = "red", linetype = "dashed") +
    geom_text_repel() +
    ylim(0, 30) +
    xlim(0, 30) +
    stat_cor() +
    ggtitle("Reliability - 0.09M Intensity") +
    theme_light() 

Warning: Removed 14 rows containing non-finite outside the scale range
(`stat_cor()`).

Warning: Removed 14 rows containing missing values or values outside the scale range
(`geom_point()`).

Warning: Removed 14 rows containing missing values or values outside the scale range
(`geom_text_repel()`).

Warning: ggrepel: 96 unlabeled data points (too many overlaps). Consider
increasing max.overlaps

    ## African
    full.df %>% rownames_to_column(var = "label") %>%
        filter(race == "African") %>%
       ggplot(aes(x = lint_d1, y = lint_d2, label = label)) +
        geom_point() +
        geom_abline(intercept = 0, slope = 1, color = "red", linetype = "dashed") +
        geom_text_repel() +
        ylim(0, 30) +
        xlim(0, 30) +
        stat_cor() +
        ggtitle("African: Reliability - 0.09M Intensity") +
        theme_light() 

Warning: Removed 9 rows containing non-finite outside the scale range
(`stat_cor()`).

Warning: Removed 9 rows containing missing values or values outside the scale range
(`geom_point()`).

Warning: Removed 9 rows containing missing values or values outside the scale range
(`geom_text_repel()`).

Warning: ggrepel: 43 unlabeled data points (too many overlaps). Consider
increasing max.overlaps

    ## Asian
    full.df %>% rownames_to_column(var = "label") %>%
        filter(race == "Asian") %>%
       ggplot(aes(x = lint_d1, y = lint_d2, label = label)) +
        geom_point() +
        geom_abline(intercept = 0, slope = 1, color = "red", linetype = "dashed") +
        geom_text_repel() +
        ylim(0, 30) +
        xlim(0, 30) +
        stat_cor() +
        ggtitle("Asian: Reliability - 0.09M Intensity") +
        theme_light() 

Warning: Removed 5 rows containing non-finite outside the scale range
(`stat_cor()`).

Warning: Removed 5 rows containing missing values or values outside the scale range
(`geom_point()`).

Warning: Removed 5 rows containing missing values or values outside the scale range
(`geom_text_repel()`).

## 1.05M intensity
full.df %>% rownames_to_column(var = "label") %>%
   ggplot(aes(x = hint_d1, y = hint_d2, label = label)) +
    geom_point() +
    geom_abline(intercept = 0, slope = 1, color = "red", linetype = "dashed") +
    geom_text_repel() +
    ylim(0, 100) +
    xlim(0, 100) +
    stat_cor() +
    ggtitle("Reliability - 1.05M Intensity") +
    theme_light() 

Warning: Removed 3 rows containing non-finite outside the scale range
(`stat_cor()`).

Warning: Removed 3 rows containing missing values or values outside the scale range
(`geom_point()`).

Warning: Removed 3 rows containing missing values or values outside the scale range
(`geom_text_repel()`).

Warning: ggrepel: 78 unlabeled data points (too many overlaps). Consider
increasing max.overlaps

    ## African
    full.df %>% rownames_to_column(var = "label") %>%
        filter(race == "African") %>%
       ggplot(aes(x = hint_d1, y = hint_d2, label = label)) +
        geom_point() +
        geom_abline(intercept = 0, slope = 1, color = "red", linetype = "dashed") +
        geom_text_repel() +
        ylim(0, 100) +
        xlim(0, 100) +
        stat_cor() +
        ggtitle("African: Reliability - 1.05M Intensity") +
        theme_light() 

Warning: Removed 3 rows containing non-finite outside the scale range
(`stat_cor()`).

Warning: Removed 3 rows containing missing values or values outside the scale range
(`geom_point()`).

Warning: Removed 3 rows containing missing values or values outside the scale range
(`geom_text_repel()`).

Warning: ggrepel: 15 unlabeled data points (too many overlaps). Consider
increasing max.overlaps

    ## Asian
    full.df %>% rownames_to_column(var = "label") %>%
        filter(race == "Asian") %>%
       ggplot(aes(x = hint_d1, y = hint_d2, label = label)) +
        geom_point() +
        geom_abline(intercept = 0, slope = 1, color = "red", linetype = "dashed") +
        geom_text_repel() +
        ylim(0, 100) +
        xlim(0, 100) +
        stat_cor() +
        ggtitle("Asian: Reliability - 1.05M Intensity") +
        theme_light() 

Warning: ggrepel: 1 unlabeled data points (too many overlaps). Consider
increasing max.overlaps

## 0.09M liking
full.df %>% rownames_to_column(var = "label") %>%
 ggplot(aes(x = llike_d1, y = llike_d2, label = label)) +
  geom_point() +
  geom_abline(intercept = 0, slope = 1, color = "red", linetype = "dashed") +
  geom_text_repel() +
  ylim(0, 100) +
  xlim(0, 100) +  
  stat_cor() +
  ggtitle("0.09M Liking Reliability") +
  theme_light() 

Warning: Removed 3 rows containing non-finite outside the scale range
(`stat_cor()`).

Warning: Removed 3 rows containing missing values or values outside the scale range
(`geom_point()`).

Warning: Removed 3 rows containing missing values or values outside the scale range
(`geom_text_repel()`).

Warning: ggrepel: 61 unlabeled data points (too many overlaps). Consider
increasing max.overlaps

  ## Africans
  full.df %>% rownames_to_column(var = "label") %>%
    filter(race == "African") %>%
 ggplot(aes(x = llike_d1, y = llike_d2, label = label)) +
  geom_point() +
  geom_abline(intercept = 0, slope = 1, color = "red", linetype = "dashed") +
  geom_text_repel() +
  ylim(0, 100) +
  xlim(0, 100) +  
  stat_cor() +
  ggtitle("0.09M Liking Reliability: African") +
  theme_light() 

Warning: Removed 3 rows containing non-finite outside the scale range
(`stat_cor()`).

Warning: Removed 3 rows containing missing values or values outside the scale range
(`geom_point()`).

Warning: Removed 3 rows containing missing values or values outside the scale range
(`geom_text_repel()`).

Warning: ggrepel: 9 unlabeled data points (too many overlaps). Consider
increasing max.overlaps

  ## Asians
    full.df %>% rownames_to_column(var = "label") %>%
    filter(race == "Asian") %>%
 ggplot(aes(x = llike_d1, y = llike_d2, label = label)) +
  geom_point() +
  geom_abline(intercept = 0, slope = 1, color = "red", linetype = "dashed") +
  geom_text_repel() +
  ylim(0, 100) +
  xlim(0, 100) +  
  stat_cor() +
  ggtitle("0.09M Liking Reliability: Asian") +
  theme_light() 

## 1.05M Liking
full.df %>% rownames_to_column(var = "label") %>%
 ggplot(aes(x = hlike_d1, y = hlike_d2, label = label)) +
  geom_point() +
  geom_abline(intercept = 0, slope = 1, color = "red", linetype = "dashed") +
  geom_text_repel() +
  ylim(0, 100) +
  xlim(0, 100) +
  stat_cor() +
  ggtitle("1.05M Liking Reliability") +
  theme_light() 

Warning: Removed 3 rows containing non-finite outside the scale range
(`stat_cor()`).

Warning: Removed 3 rows containing missing values or values outside the scale range
(`geom_point()`).

Warning: Removed 3 rows containing missing values or values outside the scale range
(`geom_text_repel()`).

Warning: ggrepel: 28 unlabeled data points (too many overlaps). Consider
increasing max.overlaps

  ## African
  full.df %>% rownames_to_column(var = "label") %>%
    filter(race == "African") %>%
   ggplot(aes(x = hlike_d1, y = hlike_d2, label = label)) +
    geom_point() +
    geom_abline(intercept = 0, slope = 1, color = "red", linetype = "dashed") +
    geom_text_repel() +
    ylim(0, 100) +
    xlim(0, 100) +
    stat_cor() +  
    ggtitle("1.05M Liking Reliability - African") +
    theme_light() 

Warning: Removed 3 rows containing non-finite outside the scale range
(`stat_cor()`).

Warning: Removed 3 rows containing missing values or values outside the scale range
(`geom_point()`).

Warning: Removed 3 rows containing missing values or values outside the scale range
(`geom_text_repel()`).

  ## Asian
  full.df %>% rownames_to_column(var = "label") %>%
    filter(race == "Asian") %>%
   ggplot(aes(x = hlike_d1, y = hlike_d2, label = label)) +
    geom_point() +
    geom_abline(intercept = 0, slope = 1, color = "red", linetype = "dashed") +
    geom_text_repel() +
    ylim(0, 100) +
    xlim(0, 100) +
    stat_cor() +
    ggtitle("1.05M Liking Reliability - Asian") +
    theme_light() 

Warning: ggrepel: 3 unlabeled data points (too many overlaps). Consider
increasing max.overlaps

Sensory reliability - slopes

## SST
full.df %>% rownames_to_column(var = "label") %>%
 ggplot(aes(x = SST.d1, y = SST.d2, label = label)) +
  geom_point() +
  geom_abline(intercept = 0, slope = 1, color = "red", linetype = "dashed") +
  geom_text_repel() +
  ylim(-100, 100) +
  xlim(-100, 100) +
  stat_cor() +
  ggtitle("Simple Sweet Test Reliability") +
  theme_light() 

Warning: Removed 3 rows containing non-finite outside the scale range
(`stat_cor()`).

Warning: Removed 3 rows containing missing values or values outside the scale range
(`geom_point()`).

Warning: Removed 3 rows containing missing values or values outside the scale range
(`geom_text_repel()`).

Warning: ggrepel: 105 unlabeled data points (too many overlaps). Consider
increasing max.overlaps

  ## African
  full.df %>% rownames_to_column(var = "label") %>%
    filter(race == "African") %>%
     ggplot(aes(x = SST.d1, y = SST.d2, label = label)) +
    geom_point() +
    geom_abline(intercept = 0, slope = 1, color = "red", linetype = "dashed") +
    geom_text_repel(segment.linetype = 3, max.overlaps = Inf) +
    ylim(-100, 100) +
    xlim(-100, 100) +
    stat_cor() +
    ggtitle("Simple Sweet Test Reliability - African") +
    theme_light() 

Warning: Removed 3 rows containing non-finite outside the scale range
(`stat_cor()`).

Warning: Removed 3 rows containing missing values or values outside the scale range
(`geom_point()`).

Warning: Removed 3 rows containing missing values or values outside the scale range
(`geom_text_repel()`).

  ## Asian
  full.df %>% rownames_to_column(var = "label") %>%
    filter(race == "Asian") %>%
   ggplot(aes(x = SST.d1, y = SST.d2, label = label)) +
    geom_point() +
    geom_abline(intercept = 0, slope = 1, color = "red", linetype = "dashed") +
    geom_text_repel(segment.linetype = 3, max.overlaps = Inf) +
    ylim(-100, 100) +
    xlim(-100, 100) +
    stat_cor() +
    ggtitle("Simple Sweet Test Reliability - Asian") +
    theme_light() 

  ## By ancestry and USB

full.df %>% rownames_to_column(var = "label") %>%
   ggplot(aes(x = SST.d1, y = SST.d2, label = label)) +
    geom_point() +
    geom_abline(intercept = 0, slope = 1, color = "red", linetype = "dashed") +
    stat_cor(method = "pearson", p.accuracy = 0.001, r.accuracy = 0.01)+
    facet_wrap(vars(race, usb)) +
    ggtitle("Sweet Liking Test-Retest Reliability") +
    theme_bw()+
    scale_x_continuous(name = "Sweet Taste Liking, Day 1", limits = c(-100, 100)) +
    scale_y_continuous(name = "Sweet Taste Liking, Day 2", limits = c(-100, 100)) +
  theme(strip.text.x = element_text(size = 12, face = "bold"))

Warning: Removed 3 rows containing non-finite outside the scale range
(`stat_cor()`).

Warning: Removed 3 rows containing missing values or values outside the scale range
(`geom_point()`).

Demographics and counts

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

# A tibble: 2 × 2
# Groups:   race [2]
  race        n
  <chr>   <int>
1 African   121
2 Asian      80

      ## Graph
      full.df %>%
        group_by(race, ageCat) %>%
        count() %>%
        group_by(race) %>%
        mutate(percent = (n / sum(n)) * 100) %>%
        ggplot(aes(x = race, y = percent, fill = ageCat)) +
        geom_bar(stat = "identity", position = position_dodge()) +
        labs(title = "Bar Graph of Percentages", x = "Ancestry", y = "Percentage") +
        theme_bw() +
        theme(strip.text.x = element_text(size = 12, face = "bold")) 

      ## Chi-square
      chisq.test(full.df$race, full.df$ageCat)

Warning in chisq.test(full.df$race, full.df$ageCat): Chi-squared approximation
may be incorrect

    Pearson's Chi-squared test

data:  full.df$race and full.df$ageCat
X-squared = 1.3155, df = 3, p-value = 0.7255

## Counts - by usb and age
full.df %>%
  group_by(usb, ageCat) %>%
  count()

# A tibble: 8 × 3
# Groups:   usb, ageCat [8]
  usb       ageCat         n
  <chr>     <chr>      <int>
1 Immigrant AgeGroup.1    56
2 Immigrant AgeGroup.2     7
3 Immigrant AgeGroup.3     4
4 Immigrant AgeGroup.4     4
5 U.S. Born AgeGroup.1   106
6 U.S. Born AgeGroup.2    22
7 U.S. Born AgeGroup.3     1
8 U.S. Born AgeGroup.4     1

      ## Graph
      full.df %>%
        group_by(usb, ageCat) %>%
        count() %>%
        group_by(usb) %>%
        mutate(percent = (n / sum(n)) * 100) %>%
        ggplot(aes(x = usb, y = percent, fill = ageCat)) +
        geom_bar(stat = "identity", position = position_dodge()) +
        labs(title = "Bar Graph of Percentages", x = "Birth Origin", y = "Percentage") +
        theme_bw() +
        theme(strip.text.x = element_text(size = 12, face = "bold")) 

      ## Chi-square
      chisq.test(full.df$usb, full.df$ageCat)

Warning in chisq.test(full.df$usb, full.df$ageCat): Chi-squared approximation
may be incorrect

    Pearson's Chi-squared test

data:  full.df$usb and full.df$ageCat
X-squared = 10.365, df = 3, p-value = 0.0157

## Counts - by usb and age
full.df %>%
  group_by(usb, ageCat) %>%
  count()

# A tibble: 8 × 3
# Groups:   usb, ageCat [8]
  usb       ageCat         n
  <chr>     <chr>      <int>
1 Immigrant AgeGroup.1    56
2 Immigrant AgeGroup.2     7
3 Immigrant AgeGroup.3     4
4 Immigrant AgeGroup.4     4
5 U.S. Born AgeGroup.1   106
6 U.S. Born AgeGroup.2    22
7 U.S. Born AgeGroup.3     1
8 U.S. Born AgeGroup.4     1

      ## Graph
      full.df %>%
        group_by(usb, ageCat) %>%
        count() %>%
        group_by(usb) %>%
        mutate(percent = (n / sum(n)) * 100) %>%
        ggplot(aes(x = usb, y = percent, fill = ageCat)) +
        geom_bar(stat = "identity", position = position_dodge()) +
        labs(title = "Bar Graph of Percentages", x = "Birth Origin", y = "Percentage") +
        theme_bw() +
        theme(strip.text.x = element_text(size = 12, face = "bold")) 

      ## Chi-square
      chisq.test(full.df$usb, full.df$ageCat)      

Warning in chisq.test(full.df$usb, full.df$ageCat): Chi-squared approximation
may be incorrect

    Pearson's Chi-squared test

data:  full.df$usb and full.df$ageCat
X-squared = 10.365, df = 3, p-value = 0.0157

## Summary table
demo.race.df <- full.df %>% 
  group_by(race) %>%
  select(-ID) %>%
  get_summary_stats()

write_xlsx(demo.race.df,"/Users/maycheung/Documents/Analyses/SweetLiking//Demographics.Race.xlsx")

demo.usb.df <- full.df %>%
  group_by(usb) %>%
  select(-ID) %>%
  get_summary_stats()

write_xlsx(demo.usb.df,"/Users/maycheung/Documents/Analyses/SweetLiking//Demographics.USB.xlsx")

demo.race.usb.df <- full.df %>% 
  group_by(race, usb) %>%
  select(-ID) %>%
  get_summary_stats()

write_xlsx(demo.race.usb.df,"/Users/maycheung/Documents/Analyses/SweetLiking//Demographics.RaceUSB.xlsx")

Differences: by ancestry

## Age and sex chi-square
chisq.test(full.df$race, full.df$sex)

Warning in chisq.test(full.df$race, full.df$sex): Chi-squared approximation may
be incorrect

    Pearson's Chi-squared test

data:  full.df$race and full.df$sex
X-squared = 7.356, df = 2, p-value = 0.02527

## Age
summary(full.df %>%
  aov(age ~ as.factor(race), .))

                 Df Sum Sq Mean Sq F value Pr(>F)
as.factor(race)   1     53    53.3   1.263  0.262
Residuals       199   8397    42.2               

## BMI
summary(full.df %>%
  aov(BMI ~ as.factor(race), .))

                 Df Sum Sq Mean Sq F value   Pr(>F)    
as.factor(race)   1   1543  1542.7   36.22 8.33e-09 ***
Residuals       199   8475    42.6                     
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## hinc
summary(full.df %>%
  aov(hinc ~ as.factor(race), .))

                 Df Sum Sq Mean Sq F value Pr(>F)  
as.factor(race)   1   16.2  16.233   5.985 0.0153 *
Residuals       199  539.8   2.712                 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## low intensity
summary(full.df %>%
  aov(lint ~ as.factor(race), .))

                 Df Sum Sq Mean Sq F value Pr(>F)
as.factor(race)   1     52   52.21   1.029  0.312
Residuals       196   9948   50.76               
3 observations deleted due to missingness

## high intensity
summary(full.df %>%
  aov(hint ~ as.factor(race), .))

                 Df Sum Sq Mean Sq F value Pr(>F)
as.factor(race)   1     75   74.55   0.257  0.613
Residuals       196  56829  289.94               
3 observations deleted due to missingness

## SST
summary(full.df %>%
  aov(SST ~ as.factor(race), .))

                 Df Sum Sq Mean Sq F value Pr(>F)  
as.factor(race)   1   8176    8176   6.154 0.0139 *
Residuals       198 263072    1329                 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
1 observation deleted due to missingness

## Sweet food liking
summary(full.df %>%
  aov(sweets.liking ~ as.factor(race), .))

                 Df Sum Sq Mean Sq F value   Pr(>F)    
as.factor(race)   1  24112   24112   27.94 3.28e-07 ***
Residuals       199 171758     863                     
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## Sugar sweetened beverages liking
summary(full.df %>%
  aov(ssb.liking ~ as.factor(race), .))

                 Df Sum Sq Mean Sq F value Pr(>F)
as.factor(race)   1   1804    1804   1.412  0.236
Residuals       199 254252    1278               

## Sweet foods and beverages liking
summary(full.df %>%
  aov(sfbl.liking ~ as.factor(race), .))

                 Df Sum Sq Mean Sq F value   Pr(>F)    
as.factor(race)   1  11143   11143   14.53 0.000184 ***
Residuals       199 152666     767                     
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## Sugar intake
summary(full.df %>%
  aov(sugar.intake ~ as.factor(race), .))

                 Df   Sum Sq Mean Sq F value   Pr(>F)    
as.factor(race)   1  1616863 1616863   27.25 4.47e-07 ***
Residuals       199 11805921   59326                     
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## sHEI
summary(full.df %>%
  aov(sHEI ~ as.factor(race), .))

                 Df Sum Sq Mean Sq F value Pr(>F)
as.factor(race)   1      2    2.08   0.024  0.878
Residuals       199  17609   88.49               

## Emotional
summary(full.df %>%
  aov(emotional ~ as.factor(race), .))

                 Df Sum Sq Mean Sq F value Pr(>F)
as.factor(race)   1     15   14.88   0.741   0.39
Residuals       199   3994   20.07               

## Uncontrolled
summary(full.df %>%
  aov(uncontrolled ~ as.factor(race), .))

                 Df Sum Sq Mean Sq F value Pr(>F)
as.factor(race)   1     12   11.66   0.482  0.488
Residuals       199   4818   24.21               

## Restraint
summary(full.df %>%
  aov(restraint ~ as.factor(race), .))

                 Df Sum Sq Mean Sq F value Pr(>F)
as.factor(race)   1    6.2   6.235   0.512  0.475
Residuals       199 2425.3  12.188               

## BAS
summary(full.df %>%
  aov(BAS ~ as.factor(race), .))

                 Df Sum Sq Mean Sq F value Pr(>F)
as.factor(race)   1      1    0.91   0.015  0.902
Residuals       199  11987   60.24               

## IES
summary(full.df %>%
  aov(IES ~ as.factor(race), .))

                 Df Sum Sq Mean Sq F value Pr(>F)  
as.factor(race)   1   1.03  1.0257   4.601 0.0332 *
Residuals       199  44.36  0.2229                 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Differences: by birth origin

## Age and sex chi-square
chisq.test(full.df$usb, full.df$sex)

Warning in chisq.test(full.df$usb, full.df$sex): Chi-squared approximation may
be incorrect

    Pearson's Chi-squared test

data:  full.df$usb and full.df$sex
X-squared = 2.8019, df = 2, p-value = 0.2464

## Age
summary(full.df %>%
  aov(age ~ as.factor(usb), .))

                Df Sum Sq Mean Sq F value  Pr(>F)   
as.factor(usb)   1    317   317.1   7.758 0.00586 **
Residuals      199   8134    40.9                   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## BMI
summary(full.df %>%
  aov(BMI ~ as.factor(usb), .))

                Df Sum Sq Mean Sq F value Pr(>F)  
as.factor(usb)   1    226   226.3     4.6 0.0332 *
Residuals      199   9791    49.2                 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## hinc
summary(full.df %>%
  aov(hinc ~ as.factor(usb), .))

                Df Sum Sq Mean Sq F value Pr(>F)
as.factor(usb)   1      0  0.0247   0.009  0.925
Residuals      199    556  2.7939               

## low intensity
summary(full.df %>%
  aov(lint ~ as.factor(usb), .))

                Df Sum Sq Mean Sq F value Pr(>F)  
as.factor(usb)   1    310  310.47    6.28  0.013 *
Residuals      196   9690   49.44                 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
3 observations deleted due to missingness

## high intensity
summary(full.df %>%
  aov(hint ~ as.factor(usb), .))

                Df Sum Sq Mean Sq F value Pr(>F)
as.factor(usb)   1    762   761.8   2.659  0.105
Residuals      196  56142   286.4               
3 observations deleted due to missingness

## SST
summary(full.df %>%
  aov(SST ~ as.factor(usb), .))

                Df Sum Sq Mean Sq F value Pr(>F)
as.factor(usb)   1    109     109    0.08  0.778
Residuals      198 271139    1369               
1 observation deleted due to missingness

## Sweet food liking
summary(full.df %>%
  aov(sweets.liking ~ as.factor(usb), .))

                Df Sum Sq Mean Sq F value Pr(>F)  
as.factor(usb)   1   3702    3702   3.833 0.0516 .
Residuals      199 192169     966                 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## Sugar sweetened beverages liking
summary(full.df %>%
  aov(ssb.liking ~ as.factor(usb), .))

                Df Sum Sq Mean Sq F value Pr(>F)
as.factor(usb)   1   2841    2841   2.233  0.137
Residuals      199 253215    1272               

## Sweet foods and beverages liking
summary(full.df %>%
  aov(sfbl.liking ~ as.factor(usb), .))

                Df Sum Sq Mean Sq F value Pr(>F)  
as.factor(usb)   1   3180    3180   3.939 0.0486 *
Residuals      199 160630     807                 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## Sugar intake
summary(full.df %>%
  aov(sugar.intake ~ as.factor(usb), .))

                Df   Sum Sq Mean Sq F value Pr(>F)
as.factor(usb)   1   129347  129347   1.936  0.166
Residuals      199 13293437   66801               

## sHEI
summary(full.df %>%
  aov(sHEI ~ as.factor(usb), .))

                Df Sum Sq Mean Sq F value Pr(>F)
as.factor(usb)   1      3    2.51   0.028  0.866
Residuals      199  17609   88.49               

## Emotional
summary(full.df %>%
  aov(emotional ~ as.factor(usb), .))

                Df Sum Sq Mean Sq F value Pr(>F)
as.factor(usb)   1      0   0.286   0.014  0.905
Residuals      199   4009  20.145               

## Uncontrolled
summary(full.df %>%
  aov(uncontrolled ~ as.factor(usb), .))

                Df Sum Sq Mean Sq F value Pr(>F)
as.factor(usb)   1      0   0.194   0.008  0.929
Residuals      199   4830  24.271               

## Restraint
summary(full.df %>%
  aov(restraint ~ as.factor(usb), .))

                Df Sum Sq Mean Sq F value Pr(>F)
as.factor(usb)   1    0.3   0.301   0.025  0.875
Residuals      199 2431.3  12.217               

## BAS
summary(full.df %>%
  aov(BAS ~ as.factor(usb), .))

                Df Sum Sq Mean Sq F value Pr(>F)   
as.factor(usb)   1    441   440.7   7.595 0.0064 **
Residuals      199  11547    58.0                  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## IES
summary(full.df %>%
  aov(IES ~ as.factor(usb), .))

                Df Sum Sq Mean Sq F value Pr(>F)
as.factor(usb)   1   0.02 0.01784   0.078   0.78
Residuals      199  45.37 0.22798               

Differences: by ancestry X birth origin

## Age
summary(full.df %>%
  aov(age ~ as.factor(usb) * as.factor(race) , .))

                                Df Sum Sq Mean Sq F value  Pr(>F)   
as.factor(usb)                   1    317   317.1   7.707 0.00603 **
as.factor(race)                  1      7     6.6   0.161 0.68878   
as.factor(usb):as.factor(race)   1     22    22.1   0.538 0.46420   
Residuals                      197   8105    41.1                   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## BMI
summary(full.df %>%
  aov(BMI ~  as.factor(usb) * as.factor(race) , .))

                                Df Sum Sq Mean Sq F value   Pr(>F)    
as.factor(usb)                   1    226   226.3   5.275   0.0227 *  
as.factor(race)                  1   1337  1337.3  31.164 7.78e-08 ***
as.factor(usb):as.factor(race)   1      1     0.6   0.014   0.9057    
Residuals                      197   8454    42.9                     
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## hinc
summary(full.df %>%
  aov(hinc ~ as.factor(usb) * as.factor(race) , .))

                                Df Sum Sq Mean Sq F value Pr(>F)  
as.factor(usb)                   1    0.0   0.025   0.009 0.9239  
as.factor(race)                  1   17.2  17.151   6.355 0.0125 *
as.factor(usb):as.factor(race)   1    7.2   7.182   2.661 0.1044  
Residuals                      197  531.7   2.699                 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## low intensity
summary(full.df %>%
  aov(lint ~ as.factor(usb) * as.factor(race) , .))

                                Df Sum Sq Mean Sq F value Pr(>F)  
as.factor(usb)                   1    310  310.47   6.318 0.0128 *
as.factor(race)                  1    157  156.97   3.195 0.0754 .
as.factor(usb):as.factor(race)   1      0    0.31   0.006 0.9365  
Residuals                      194   9533   49.14                 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
3 observations deleted due to missingness

## high intensity
summary(full.df %>%
  aov(hint ~ as.factor(usb) * as.factor(race) , .))

                                Df Sum Sq Mean Sq F value Pr(>F)
as.factor(usb)                   1    762   761.8   2.638  0.106
as.factor(race)                  1      1     1.3   0.004  0.948
as.factor(usb):as.factor(race)   1    129   128.7   0.446  0.505
Residuals                      194  56012   288.7               
3 observations deleted due to missingness

## SST
summary(full.df %>%
  aov(SST ~ as.factor(usb) * as.factor(race) , .))

                                Df Sum Sq Mean Sq F value Pr(>F)  
as.factor(usb)                   1    109     109   0.081 0.7757  
as.factor(race)                  1   8274    8274   6.184 0.0137 *
as.factor(usb):as.factor(race)   1    612     612   0.457 0.4997  
Residuals                      196 262254    1338                 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
1 observation deleted due to missingness

## Sweet food liking
summary(full.df %>%
  aov(sweets.liking ~ as.factor(usb) * as.factor(race) , .))

                                Df Sum Sq Mean Sq F value  Pr(>F)    
as.factor(usb)                   1   3702    3702   4.256  0.0404 *  
as.factor(race)                  1  20791   20791  23.900 2.1e-06 ***
as.factor(usb):as.factor(race)   1      8       8   0.009  0.9249    
Residuals                      197 171370     870                    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## Sugar sweetened beverages liking
summary(full.df %>%
  aov(ssb.liking ~ as.factor(race) * as.factor(race) , .))

                 Df Sum Sq Mean Sq F value Pr(>F)
as.factor(race)   1   1804    1804   1.412  0.236
Residuals       199 254252    1278               

## Sweet foods and beverages liking
summary(full.df %>%
  aov(sfbl.liking ~ as.factor(usb) * as.factor(race) , .))

                                Df Sum Sq Mean Sq F value   Pr(>F)    
as.factor(usb)                   1   3180    3180   4.135 0.043358 *  
as.factor(race)                  1   8797    8797  11.439 0.000867 ***
as.factor(usb):as.factor(race)   1    336     336   0.436 0.509672    
Residuals                      197 151498     769                     
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## Sugar intake
summary(full.df %>%
  aov(sugar.intake ~ as.factor(usb) * as.factor(race) , .))

                                Df   Sum Sq Mean Sq F value   Pr(>F)    
as.factor(usb)                   1   129347  129347   2.159    0.143    
as.factor(race)                  1  1487763 1487763  24.834 1.36e-06 ***
as.factor(usb):as.factor(race)   1     3514    3514   0.059    0.809    
Residuals                      197 11802160   59909                     
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## sHEI
summary(full.df %>%
  aov(sHEI ~ as.factor(usb) * as.factor(race) , .))

                                Df Sum Sq Mean Sq F value Pr(>F)
as.factor(usb)                   1      3    2.51   0.028  0.866
as.factor(race)                  1      4    3.78   0.043  0.836
as.factor(usb):as.factor(race)   1    180  180.48   2.040  0.155
Residuals                      197  17424   88.45               

## Emotional
summary(full.df %>%
  aov(emotional ~ as.factor(race) * as.factor(usb), .))

                                Df Sum Sq Mean Sq F value Pr(>F)
as.factor(race)                  1     15  14.882   0.735  0.392
as.factor(usb)                   1      0   0.281   0.014  0.906
as.factor(race):as.factor(usb)   1      6   5.503   0.272  0.603
Residuals                      197   3989  20.246               

## Uncontrolled
summary(full.df %>%
  aov(uncontrolled ~ as.factor(race) * as.factor(usb), .))

                                Df Sum Sq Mean Sq F value Pr(>F)
as.factor(race)                  1     12  11.665   0.478  0.490
as.factor(usb)                   1      0   0.254   0.010  0.919
as.factor(race):as.factor(usb)   1     11  11.403   0.467  0.495
Residuals                      197   4807  24.400               

## Restraint
summary(full.df %>%
  aov(restraint ~ as.factor(race) * as.factor(usb), .))

                                Df Sum Sq Mean Sq F value Pr(>F)
as.factor(race)                  1    6.2   6.235   0.506  0.478
as.factor(usb)                   1    0.0   0.018   0.001  0.970
as.factor(race):as.factor(usb)   1    0.1   0.071   0.006  0.939
Residuals                      197 2425.2  12.311               

## BAS
summary(full.df %>%
  aov(BAS ~ as.factor(race) * as.factor(usb), .))

                                Df Sum Sq Mean Sq F value  Pr(>F)   
as.factor(race)                  1      1     0.9   0.016 0.90091   
as.factor(usb)                   1    464   464.0   7.958 0.00528 **
as.factor(race):as.factor(usb)   1     35    35.5   0.609 0.43628   
Residuals                      197  11488    58.3                   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## IES
summary(full.df %>%
  aov(IES ~ as.factor(race) * as.factor(usb), .))

                                Df Sum Sq Mean Sq F value Pr(>F)  
as.factor(race)                  1   1.03  1.0257   4.563 0.0339 *
as.factor(usb)                   1   0.02  0.0214   0.095 0.7579  
as.factor(race):as.factor(usb)   1   0.06  0.0588   0.262 0.6095  
Residuals                      197  44.28  0.2248                 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Differences: by sweet liking phenotype

## Sweet liking phenotype by race only

      ## Graph
      full.df %>%
        filter(!is.na(sPhen)) %>%
        group_by(race, sPhen) %>%
        count() %>%
        group_by(race) %>%
        mutate(percent = (n / sum(n)) * 100) %>%
        ggplot(aes(x = race, y = percent, fill = as.factor(sPhen))) +
        geom_bar(stat = "identity", position = position_dodge()) +
        labs(title = "Bar Graph of Percentages", x = "Ancestry", y = "Percentage") +
        theme_bw() +
        theme(strip.text.x = element_text(size = 12, face = "bold")) +
        scale_fill_manual(values = c("#689F38","#FBC02D", "#C62828" ))

      ## Chi-square
      chisq.test(full.df$race, full.df$sPhen)

    Pearson's Chi-squared test

data:  full.df$race and full.df$sPhen
X-squared = 4.7593, df = 2, p-value = 0.09258

## Sweet liking phenotype by race and USB

      ## Graph
      full.df %>%
        filter(!is.na(sPhen)) %>%
        group_by(race, sPhen, usb) %>%
        count() %>%
        group_by(race, usb) %>%
        mutate(percent = (n / sum(n)) * 100) %>%
        ggplot(aes(x = usb, y = percent, fill = as.factor(sPhen))) +
        geom_bar(stat = "identity", position = position_dodge()) +
        labs(title = "Bar Graph of Percentages", x = "Birth Origin", y = "Percentage") +
        facet_wrap(vars(race)) +
        theme_bw() +
        theme(strip.text.x = element_text(size = 12, face = "bold")) +
        scale_fill_manual(values = c("#689F38","#FBC02D", "#C62828" ))

      ## Chi-square
      
          ## Overall
          chisq.test(full.df$usb, full.df$sPhen)

    Pearson's Chi-squared test

data:  full.df$usb and full.df$sPhen
X-squared = 0.75298, df = 2, p-value = 0.6863

          ## African
          full.df %>%
            filter(race == "African") %>%
            count(usb = factor(usb), sPhen = as.factor(sPhen)) %>%
            pivot_wider(names_from = sPhen, values_from = n, values_fill = list(n = 0)) %>%
            column_to_rownames(., "usb") %>%
            chisq.test(usb.)

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

    Pearson's Chi-squared test

data:  .
X-squared = 0.39348, df = 3, p-value = 0.9416

          ## Asian
          full.df %>%
            filter(race == "Asian") %>%
            count(usb = factor(usb), sPhen = factor(sPhen)) %>%
            pivot_wider(names_from = sPhen, values_from = n, values_fill = list(n = 0)) %>%
            column_to_rownames(., "usb") %>%
            chisq.test(usb.)

    Pearson's Chi-squared test

data:  .
X-squared = 0.37602, df = 2, p-value = 0.8286

Correlations: Corrplot, SST x diet

## Corr DF
corr.df <- full.df %>%
  select(age, race, usb, education, hinc, BMI, sfbl.liking, restraint, uncontrolled, emotional, sugar.intake, sHEI, SST) 

    ## All Participants
    all.corr <- corr.df %>%
      select(-race, -usb) %>%
      cor(use = "pairwise.complete.obs")
    
    all.res <- corr.df %>%
      select(-c("race", "usb")) %>% 
      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("All Participants")       

    ## All Africans
    Afr.corr <- corr.df %>%
      filter(race == "African") %>%
      select(-race, -usb) %>%
      cor(use = "pairwise.complete.obs")
    
    Afr.res <- corr.df %>%
      filter(race == "African") %>%
      select(-c("race", "usb")) %>% 
      cor.mtest(., conf.level = 0.95)
    
    corrplot(Afr.corr, p.mat = Afr.res$p, method = 'circle',tl.col = 'black', type = 'lower', insig='blank', addCoef.col ='black', number.cex = 0.7, order = 'original', diag=FALSE) 
      title("All Africans")  

    ## All East Asians
    Asi.corr <- corr.df %>%
      filter(race == "Asian") %>%
      select(-race, -usb) %>%
      cor(use = "pairwise.complete.obs")
    
    Asi.res <- corr.df %>%
      filter(race == "Asian") %>%
      select(-c("race", "usb")) %>% 
      cor.mtest(., conf.level = 0.95)
    
    corrplot(Asi.corr, p.mat = Asi.res$p, method = 'circle',tl.col = 'black', type = 'lower', insig='blank', addCoef.col ='black', number.cex = 0.7, order = 'original', diag=FALSE) 
      title("All East Asians")     

    ## All USB
    USB.corr <- corr.df %>%
      filter(usb == "U.S. Born") %>%
      select(-race, -usb) %>%
      cor(use = "pairwise.complete.obs")
    
    USB.res <- corr.df %>%
      filter(race == "African") %>%
      select(-c("race", "usb")) %>% 
      cor.mtest(., conf.level = 0.95)
    
    corrplot(USB.corr, p.mat = USB.res$p, method = 'circle',tl.col = 'black', type = 'lower', insig='blank', addCoef.col ='black', number.cex = 0.7, order = 'original', diag=FALSE) 
      title("All USB")  

    ## All Immigrants
    IMG.corr <- corr.df %>%
      filter(usb == "Immigrant") %>%
      select(-race, -usb) %>%
      cor(use = "pairwise.complete.obs")
    
    IMG.res <- corr.df %>%
      filter(usb == "Immigrant") %>%
      select(-c("race", "usb")) %>% 
      cor.mtest(., conf.level = 0.95)
    
    corrplot(IMG.corr, p.mat = IMG.res$p, method = 'circle',tl.col = 'black', type = 'lower', insig='blank', addCoef.col ='black', number.cex = 0.7, order = 'original', diag=FALSE) 
      title("All Immigrants")    

    ## African USB
    AfrUSB.corr <- corr.df %>%
      filter(race == "African", usb == "U.S. Born") %>%
      select(-race, -usb) %>%
      cor(use = "pairwise.complete.obs")
    
    AfrUSB.res <- corr.df %>%
      filter(race == "African", usb == "U.S. Born") %>%
      select(-c("race", "usb")) %>% 
      cor.mtest(., conf.level = 0.95)
    
    corrplot(AfrUSB.corr, p.mat = AfrUSB.res$p, method = 'circle',tl.col = 'black', type = 'lower', insig='blank', addCoef.col ='black', number.cex = 0.7, order = 'original', diag=FALSE)

  title("African, U.S. Born")    

    ## African immigrant
    AfrIm.corr <- corr.df %>%
      filter(race == "African", usb == "Immigrant") %>%
      select(-race, -usb) %>%
      cor(use = "pairwise.complete.obs")
    
    AfrIm.res <- corr.df %>%
      filter(race == "African", usb == "Immigrant") %>%
      select(-c("race", "usb")) %>% 
      cor.mtest(., conf.level = 0.95)
    
    corrplot(AfrIm.corr, p.mat = AfrIm.res$p, method = 'circle',tl.col = 'black', type = 'lower', insig='blank', addCoef.col ='black', number.cex = 0.7, order = 'original', diag=FALSE)

    title("African, Immigrants")       

    ## Asian USB
    AsiUSB.corr <- corr.df %>%
      filter(race == "Asian", usb == "U.S. Born") %>%
      select(-race, -usb) %>%
      cor(use = "pairwise.complete.obs")
    
    AsiUSB.res <- corr.df %>%
      filter(race == "Asian", usb == "U.S. Born") %>%
      select(-c("race", "usb")) %>% 
      cor.mtest(., conf.level = 0.95)
    
      corrplot(AsiUSB.corr, p.mat = AsiUSB.res$p, method = 'circle',tl.col = 'black', type = 'lower', insig='blank', addCoef.col ='black', number.cex = 0.7, order = 'original', diag=FALSE)
      
    title("Asian, U.S. Born")       

    ## Asian immigrant
    AsiIm.corr <- corr.df %>%
      filter(race == "Asian", usb == "Immigrant") %>%
      select(-race, -usb) %>%
      cor(use = "pairwise.complete.obs")
    
    AsiIm.res <- corr.df %>%
      filter(race == "Asian", usb == "Immigrant") %>%
      select(-c("race", "usb")) %>% 
      cor.mtest(., conf.level = 0.95)
    
    corrplot(AsiIm.corr, p.mat = AsiIm.res$p, method = 'circle',tl.col = 'black', type = 'lower', insig='blank', addCoef.col ='black', number.cex = 0.7, order = 'original', diag=FALSE)

      title("Asian, Immigrants")       

Linear regression - sugar intake, African

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

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

Warning: Removed 1 row containing non-finite outside the scale range
(`stat_smooth()`).

Warning: Removed 1 row containing non-finite outside the scale range
(`stat_cor()`).

Warning: Removed 1 row containing missing values or values outside the scale range
(`geom_point()`).

summary(lm(sugar.intake ~ SST, data = African.df))

Call:
lm(formula = sugar.intake ~ SST, data = African.df)

Residuals:
    Min      1Q  Median      3Q     Max 
-556.40 -184.13  -28.16  184.73  670.37 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 546.9573    26.7205  20.470  < 2e-16 ***
SST           1.9124     0.6499   2.943  0.00392 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 267 on 118 degrees of freedom
  (1 observation deleted due to missingness)
Multiple R-squared:  0.06836,   Adjusted R-squared:  0.06047 
F-statistic: 8.659 on 1 and 118 DF,  p-value: 0.00392

tidy(lm(sugar.intake ~ SST, data = African.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.47e237    26.7       20.5  1.23e-40 3.63e214  3.32e260
2 SST         6.77e  0     0.650      2.94 3.92e- 3 1.87e  0  2.45e  1

## Sweet foods and beverages liking on sugar intake 
African.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 = African.df))

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

Residuals:
    Min      1Q  Median      3Q     Max 
-534.22 -151.55  -37.37  157.67  657.19 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 390.2278    41.9722   9.297 8.64e-16 ***
sfbl.liking   4.3382     0.8083   5.367 4.01e-07 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 247.5 on 119 degrees of freedom
Multiple R-squared:  0.1949,    Adjusted R-squared:  0.1881 
F-statistic:  28.8 on 1 and 119 DF,  p-value: 4.014e-07

tidy(lm(sugar.intake ~ sfbl.liking, data = African.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.98e169    42.0        9.30 8.64e-16 2.40e133  3.70e205
2 sfbl.liking 7.66e  1     0.808      5.37 4.01e- 7 1.55e  1  3.79e  2

## BMI on sugar intake 
African.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 = African.df))

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

Residuals:
    Min      1Q  Median      3Q     Max 
-494.98 -223.30  -10.25  179.05  723.12 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  402.402     91.860   4.381 2.56e-05 ***
BMI            6.253      3.109   2.011   0.0465 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 271.3 on 119 degrees of freedom
Multiple R-squared:  0.03288,   Adjusted R-squared:  0.02475 
F-statistic: 4.046 on 1 and 119 DF,  p-value: 0.04655

tidy(lm(sugar.intake ~ BMI, data = African.df), exponentiate = TRUE, conf.int = TRUE)

# A tibble: 2 × 7
  term        estimate std.error statistic   p.value conf.low conf.high
  <chr>          <dbl>     <dbl>     <dbl>     <dbl>    <dbl>     <dbl>
1 (Intercept) 5.76e174     91.9       4.38 0.0000256  5.84e95  5.69e253
2 BMI         5.20e  2      3.11      2.01 0.0465     1.10e 0  2.45e  5

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

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

Residuals:
    Min      1Q  Median      3Q     Max 
-534.22 -151.55  -37.37  157.67  657.19 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 390.2278    41.9722   9.297 8.64e-16 ***
sfbl.liking   4.3382     0.8083   5.367 4.01e-07 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 247.5 on 119 degrees of freedom
Multiple R-squared:  0.1949,    Adjusted R-squared:  0.1881 
F-statistic:  28.8 on 1 and 119 DF,  p-value: 4.014e-07

Afr.sugar.lm.M2 <- lm(sugar.intake ~ sfbl.liking + SST, data = African.df) 
summary(Afr.sugar.lm.M2) # r2 = 0.22

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

Residuals:
    Min      1Q  Median      3Q     Max 
-508.85 -157.46  -36.48  187.75  608.39 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 378.2230    41.5895   9.094 2.95e-15 ***
sfbl.liking   4.0371     0.8066   5.005 1.99e-06 ***
SST           1.4077     0.6009   2.343   0.0208 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 243.4 on 117 degrees of freedom
  (1 observation deleted due to missingness)
Multiple R-squared:  0.2327,    Adjusted R-squared:  0.2195 
F-statistic: 17.74 on 2 and 117 DF,  p-value: 1.871e-07

Afr.sugar.lm.M3 <- lm(sugar.intake ~ sfbl.liking + SST + BMI, data = African.df) 
summary(Afr.sugar.lm.M3) # r2 = 0.22

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

Residuals:
    Min      1Q  Median      3Q     Max 
-485.01 -163.08  -43.18  185.58  642.41 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 285.5262    84.8455   3.365  0.00104 ** 
sfbl.liking   3.8190     0.8233   4.639 9.28e-06 ***
SST           1.4348     0.5998   2.392  0.01836 *  
BMI           3.5722     2.8521   1.252  0.21291    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 242.8 on 116 degrees of freedom
  (1 observation deleted due to missingness)
Multiple R-squared:  0.2429,    Adjusted R-squared:  0.2233 
F-statistic:  12.4 on 3 and 116 DF,  p-value: 4.285e-07

Afr.sugar.lm.M4 <- lm(sugar.intake ~ sfbl.liking*SST, data = African.df) 
summary(Afr.sugar.lm.M4) # r2 = 0.22

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

Residuals:
    Min      1Q  Median      3Q     Max 
-497.92 -155.03  -33.96  181.84  612.18 

Coefficients:
                 Estimate Std. Error t value Pr(>|t|)    
(Intercept)     388.20472   42.61178   9.110 2.88e-15 ***
sfbl.liking       3.67422    0.87538   4.197 5.32e-05 ***
SST               0.52176    1.02664   0.508    0.612    
sfbl.liking:SST   0.02280    0.02143   1.064    0.290    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 243.3 on 116 degrees of freedom
  (1 observation deleted due to missingness)
Multiple R-squared:  0.2401,    Adjusted R-squared:  0.2204 
F-statistic: 12.21 on 3 and 116 DF,  p-value: 5.288e-07

Afr.sugar.lm.M5 <- lm(sugar.intake ~ sfbl.liking*SST + BMI, data = African.df) 
summary(Afr.sugar.lm.M5) # r2 = 0.23

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

Residuals:
    Min      1Q  Median      3Q     Max 
-470.96 -157.64  -36.74  185.62  649.38 

Coefficients:
                 Estimate Std. Error t value Pr(>|t|)    
(Intercept)     289.05790   84.75501   3.411 0.000895 ***
sfbl.liking       3.39857    0.89579   3.794 0.000238 ***
SST               0.45405    1.02423   0.443 0.658374    
BMI               3.86288    2.85788   1.352 0.179137    
sfbl.liking:SST   0.02530    0.02144   1.180 0.240352    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 242.4 on 115 degrees of freedom
  (1 observation deleted due to missingness)
Multiple R-squared:  0.252, Adjusted R-squared:  0.2259 
F-statistic: 9.683 on 4 and 115 DF,  p-value: 8.729e-07

## Model comparison
anova(Afr.sugar.lm.M5, Afr.sugar.lm.M3)

Analysis of Variance Table

Model 1: sugar.intake ~ sfbl.liking * SST + BMI
Model 2: sugar.intake ~ sfbl.liking + SST + BMI
  Res.Df     RSS Df Sum of Sq      F Pr(>F)
1    115 6756874                           
2    116 6838714 -1    -81840 1.3929 0.2404

Logistic regression - sugar intake, African

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

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

Warning: Removed 1 row containing non-finite outside the scale range
(`stat_smooth()`).

Warning: Removed 1 row containing non-finite outside the scale range
(`stat_cor()`).

Warning: Removed 1 row containing missing values or values outside the scale range
(`geom_point()`).

summary(multinom(SugarCat ~ SST, data = African.df)) # AIC = 256.0

# weights:  9 (4 variable)
initial  value 131.833475 
final  value 124.006228 
converged

Call:
multinom(formula = SugarCat ~ SST, data = African.df)

Coefficients:
  (Intercept)          SST
2   0.2580705 0.0001314602
3   0.6141680 0.0088287878

Std. Errors:
  (Intercept)         SST
2   0.2668980 0.006825261
3   0.2541097 0.006340889

Residual Deviance: 248.0125 
AIC: 256.0125 

tidy(multinom(SugarCat ~ SST, data = African.df), exponentiate = TRUE, conf.int = TRUE)

# weights:  9 (4 variable)
initial  value 131.833475 
final  value 124.006228 
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.29   0.267      0.967   0.334     0.767      2.18
2 2       SST             1.00   0.00683    0.0193  0.985     0.987      1.01
3 3       (Intercept)     1.85   0.254      2.42    0.0157    1.12       3.04
4 3       SST             1.01   0.00634    1.39    0.164     0.996      1.02

## Sweet foods and beverages liking on sugar intake 
African.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 = African.df)) # AIC = 239.4

# weights:  9 (4 variable)
initial  value 132.932087 
final  value 115.685242 
converged

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

Coefficients:
  (Intercept) sfbl.liking
2  -0.5022691  0.02243138
3  -0.9394810  0.04235171

Std. Errors:
  (Intercept) sfbl.liking
2   0.4283089 0.009854636
3   0.4821049 0.010472970

Residual Deviance: 231.3705 
AIC: 239.3705 

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

# weights:  9 (4 variable)
initial  value 132.932087 
final  value 115.685242 
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.605   0.428       -1.17 0.241        0.261      1.40
2 2       sfbl.liking    1.02    0.00985      2.28 0.0228       1.00       1.04
3 3       (Intercept)    0.391   0.482       -1.95 0.0513       0.152      1.01
4 3       sfbl.liking    1.04    0.0105       4.04 0.0000526    1.02       1.06

## BMI on sugar intake 
African.df %>%
  ggplot(aes(as.numeric(BMIcat), SugarCat)) +
  geom_point() +
  geom_smooth(method = lm) +
  stat_cor() +
  theme_bw()

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

summary(multinom(SugarCat ~ as.numeric(BMIcat), data = African.df)) # AIC = 256.6

# weights:  9 (4 variable)
initial  value 132.932087 
iter  10 value 121.210118
iter  10 value 121.210118
final  value 121.210118 
converged

Call:
multinom(formula = SugarCat ~ as.numeric(BMIcat), data = African.df)

Coefficients:
  (Intercept) as.numeric(BMIcat)
2   -1.170748          0.8044177
3   -0.692988          0.8255954

Std. Errors:
  (Intercept) as.numeric(BMIcat)
2   0.6053199          0.3127138
3   0.5389737          0.2875646

Residual Deviance: 242.4202 
AIC: 250.4202 

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

# weights:  9 (4 variable)
initial  value 132.932087 
iter  10 value 120.801777
iter  10 value 120.801777
final  value 120.801777 
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.0505    1.20       -2.48 0.0130   0.00479     0.532
2 2       as.numeric(BM…   1.13      0.0454      2.68 0.00737  1.03        1.23 
3 3       (Intercept)      0.0952    1.13       -2.08 0.0374   0.0104      0.872
4 3       as.numeric(BM…   1.12      0.0434      2.71 0.00667  1.03        1.22 

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

# weights:  9 (4 variable)
initial  value 132.932087 
final  value 115.685242 
converged

summary(Afr.sugar.logit.M1) # AIC = 239.4

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

Coefficients:
  (Intercept) sfbl.liking
2  -0.5022691  0.02243138
3  -0.9394810  0.04235171

Std. Errors:
  (Intercept) sfbl.liking
2   0.4283089 0.009854636
3   0.4821049 0.010472970

Residual Deviance: 231.3705 
AIC: 239.3705 

Afr.sugar.logit.M2 <- multinom(SugarCat ~ sfbl.liking + SST, data = African.df)

# weights:  12 (6 variable)
initial  value 131.833475 
iter  10 value 114.140363
final  value 114.140357 
converged

summary(Afr.sugar.logit.M2) # AIC = 240.3

Call:
multinom(formula = SugarCat ~ sfbl.liking + SST, data = African.df)

Coefficients:
  (Intercept) sfbl.liking          SST
2  -0.4881552  0.02274737 -0.002330884
3  -1.0048809  0.04137997  0.004610894

Std. Errors:
  (Intercept) sfbl.liking         SST
2   0.4283979 0.009979187 0.006978600
3   0.4887997 0.010580605 0.006836799

Residual Deviance: 228.2807 
AIC: 240.2807 

Afr.sugar.logit.M3 <- multinom(SugarCat ~ sfbl.liking + SST + BMI, data = African.df)

# weights:  15 (8 variable)
initial  value 131.833475 
iter  10 value 111.066898
final  value 110.942270 
converged

summary(Afr.sugar.logit.M3) # AIC = 237.9

Call:
multinom(formula = SugarCat ~ sfbl.liking + SST + BMI, data = African.df)

Coefficients:
  (Intercept) sfbl.liking          SST        BMI
2   -2.981779  0.01712012 -0.001779845 0.10002191
3   -3.295255  0.03653191  0.005228128 0.09190839

Std. Errors:
  (Intercept) sfbl.liking         SST        BMI
2    1.208405  0.01023240 0.007151602 0.04556570
3    1.212543  0.01070755 0.006962590 0.04435979

Residual Deviance: 221.8845 
AIC: 237.8845 

Afr.sugar.logit.M4 <- multinom(SugarCat ~ sfbl.liking*SST, data = African.df)

# weights:  15 (8 variable)
initial  value 131.833475 
iter  10 value 114.083380
final  value 113.868426 
converged

summary(Afr.sugar.logit.M4) # AIC = 243.7

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

Coefficients:
  (Intercept) sfbl.liking          SST sfbl.liking:SST
2  -0.4558896  0.02052427 -0.008941733    0.0002027997
3  -0.9866636  0.04003002  0.001931134    0.0001051114

Std. Errors:
  (Intercept) sfbl.liking        SST sfbl.liking:SST
2   0.4335414  0.01040347 0.01184906    0.0002828817
3   0.5051038  0.01125520 0.01274526    0.0002903209

Residual Deviance: 227.7369 
AIC: 243.7369 

Afr.sugar.logit.M5 <- multinom(SugarCat ~ sfbl.liking*SST + BMI, data = African.df)

# weights:  18 (10 variable)
initial  value 131.833475 
iter  10 value 111.751337
final  value 110.565280 
converged

summary(Afr.sugar.logit.M5) # AIC = 241.1

Call:
multinom(formula = SugarCat ~ sfbl.liking * SST + BMI, data = African.df)

Coefficients:
  (Intercept) sfbl.liking          SST        BMI sfbl.liking:SST
2   -3.010557  0.01443035 -0.009945624 0.10238590    0.0002451509
3   -3.295205  0.03468735  0.001392988 0.09299244    0.0001369054

Std. Errors:
   (Intercept) sfbl.liking        SST        BMI sfbl.liking:SST
2 0.0006244147  0.01063575 0.01251034 0.01642007    0.0002914043
3 0.0005062265  0.01093835 0.01320953 0.01789327    0.0002959124

Residual Deviance: 221.1306 
AIC: 241.1306 

## Model comparison
anova(Afr.sugar.logit.M3, Afr.sugar.logit.M2)

Likelihood ratio tests of Multinomial Models

Response: SugarCat
                    Model Resid. df Resid. Dev   Test    Df LR stat.    Pr(Chi)
1       sfbl.liking + SST       234   228.2807                                 
2 sfbl.liking + SST + BMI       232   221.8845 1 vs 2     2 6.396173 0.04084028

Linear regression - sugar intake, Asian

## Sweet foods and beverages liking on sugar intake 
Asian.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 = Asian.df))

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

Residuals:
    Min      1Q  Median      3Q     Max 
-329.63 -100.61  -14.05   87.97  776.82 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 348.6603    29.5657   11.79   <2e-16 ***
sfbl.liking   1.6942     0.7496    2.26   0.0266 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 181.9 on 78 degrees of freedom
Multiple R-squared:  0.06146,   Adjusted R-squared:  0.04943 
F-statistic: 5.108 on 1 and 78 DF,  p-value: 0.02661

tidy(lm(sugar.intake ~ sfbl.liking, data = Asian.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.64e151    29.6       11.8  5.16e-19 7.22e125  9.64e176
2 sfbl.liking 5.44e  0     0.750      2.26 2.66e- 2 1.22e  0  2.42e  1

## BMI on sugar intake 
Asian.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 = Asian.df))

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

Residuals:
    Min      1Q  Median      3Q     Max 
-340.10 -127.07    4.83   66.00  722.95 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)  
(Intercept)   95.901    143.246   0.669   0.5052  
BMI           13.212      6.218   2.125   0.0368 *
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 182.6 on 78 degrees of freedom
Multiple R-squared:  0.05471,   Adjusted R-squared:  0.04259 
F-statistic: 4.514 on 1 and 78 DF,  p-value: 0.03678

tidy(lm(sugar.intake ~ BMI, data = Asian.df), exponentiate = TRUE, conf.int = TRUE)

# A tibble: 2 × 7
  term        estimate std.error statistic p.value conf.low conf.high
  <chr>          <dbl>     <dbl>     <dbl>   <dbl>    <dbl>     <dbl>
1 (Intercept)  4.46e41    143.       0.669  0.505  6.26e-83  3.18e165
2 BMI          5.47e 5      6.22     2.12   0.0368 2.30e+ 0  1.30e 11

## Models
Asi.sugar.lm.M1 <- lm(sugar.intake ~ sfbl.liking + BMI, data = African.df) 
summary(Asi.sugar.lm.M1) # r2 = 0.2

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

Residuals:
    Min      1Q  Median      3Q     Max 
-513.35 -160.68  -32.52  173.88  633.86 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  307.499     85.844   3.582 0.000496 ***
sfbl.liking    4.147      0.826   5.021 1.84e-06 ***
BMI            3.201      2.899   1.104 0.271641    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 247.3 on 118 degrees of freedom
Multiple R-squared:  0.2031,    Adjusted R-squared:  0.1896 
F-statistic: 15.04 on 2 and 118 DF,  p-value: 1.522e-06

Logistic regression - sugar intake, Asian

## Sweet foods and beverages liking on sugar intake 
Asian.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 = Asian.df)) # AIC = 153.9

# weights:  9 (4 variable)
initial  value 87.888983 
iter  10 value 72.969989
final  value 72.969977 
converged

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

Coefficients:
  (Intercept) sfbl.liking
2  -0.8032714  0.01988917
3  -2.2036881  0.02028945

Std. Errors:
  (Intercept) sfbl.liking
2   0.3963950  0.01019039
3   0.6936231  0.01694458

Residual Deviance: 145.94 
AIC: 153.94 

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

# weights:  9 (4 variable)
initial  value 87.888983 
iter  10 value 72.969989
final  value 72.969977 
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.448    0.396      -2.03 0.0427    0.206      0.974
2 2       sfbl.liking    1.02     0.0102      1.95 0.0510    1.00       1.04 
3 3       (Intercept)    0.110    0.694      -3.18 0.00149   0.0283     0.430
4 3       sfbl.liking    1.02     0.0169      1.20 0.231     0.987      1.05 

## BMI on sugar intake 
Asian.df %>%
  ggplot(aes(as.numeric(BMIcat), SugarCat)) +
  geom_point() +
  geom_smooth(method = lm) +
  stat_cor() +
  theme_bw()

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

summary(multinom(SugarCat ~ as.numeric(BMIcat), data = Asian.df)) # AIC = 159.2

# weights:  9 (4 variable)
initial  value 87.888983 
iter  10 value 74.299094
iter  10 value 74.299094
final  value 74.299094 
converged

Call:
multinom(formula = SugarCat ~ as.numeric(BMIcat), data = Asian.df)

Coefficients:
  (Intercept) as.numeric(BMIcat)
2  -0.8620788          0.4032553
3  -1.9954063          0.2533913

Std. Errors:
  (Intercept) as.numeric(BMIcat)
2   0.4920259          0.2688414
3   0.7881725          0.4315884

Residual Deviance: 148.5982 
AIC: 156.5982 

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

# weights:  9 (4 variable)
initial  value 87.888983 
iter  10 value 74.686453
iter  10 value 74.686452
final  value 74.686452 
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.140     1.71       -1.15   0.252  4.86e-3      4.03
2 2       as.numeric(BM…   1.08      0.0746      1.03   0.304  9.33e-1      1.25
3 3       (Intercept)      0.0147    2.66       -1.59   0.112  8.00e-5      2.69
4 3       as.numeric(BM…   1.12      0.113       1.01   0.313  8.98e-1      1.40

## Models
Asi.sugar.logit.M1 <- multinom(SugarCat ~ sfbl.liking + BMI, data = Asian.df)

# weights:  12 (6 variable)
initial  value 87.888983 
iter  10 value 72.353709
final  value 72.329147 
converged

summary(Asi.sugar.logit.M1) # AIC = 156.7

Call:
multinom(formula = SugarCat ~ sfbl.liking + BMI, data = Asian.df)

Coefficients:
  (Intercept) sfbl.liking        BMI
2   -2.380473  0.01949257 0.06997633
3   -4.601081  0.01944413 0.10574328

Std. Errors:
  (Intercept) sfbl.liking        BMI
2    1.770689  0.01018217 0.07639325
3    2.700675  0.01683093 0.11353218

Residual Deviance: 144.6583 
AIC: 156.6583