Sweet Foods Liking and Diet Quality

Functions

Useful function that will be used throughout the analyses.

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

Load packages and call data

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: 388
alpha: 0.708

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: 412
alpha: 0.557

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: 409
alpha: 0.557

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: 382
alpha: 0.686

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: 250
alpha: 0.762

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: 402
alpha: 0.701

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: 401
alpha: 0.567

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: 377
alpha: 0.504

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: 399
alpha: 0.563

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: 369
alpha: 0.612

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: 413
alpha: 0.763

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: 365
alpha: 0.589

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: 362
alpha: 0.778

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: 308
alpha: 0.511

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: 395
alpha: 0.842

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)/2, 
         fruit.liking = rowMeans(select(., "FL-07_3":"FL-11_3"), na.rm = TRUE)/2, 
         saltyfat.liking = rowMeans(select(., "FL-12_3":"FL-15_3"), na.rm = TRUE)/2, 
         hfprotein.liking = rowMeans(select(., "FL-16_3":"FL-20_3"), na.rm = TRUE)/2,   
         alcohol.liking = rowMeans(select(., "FL-22_3":"FL-25_3"), na.rm = TRUE)/2, 
         wg.liking = rowMeans(select(., "FL-26_3":"FL-29_3"), na.rm = TRUE)/2, 
         carbs.liking = rowMeans(select(., "FL-30_3":"FL-34_3"), na.rm = TRUE)/2, 
         healthyfat.liking = rowMeans(select(., "FL-35_3":"FL-38_3"), na.rm = TRUE)/2, 
        sweets.liking = rowMeans(select(., "FL-39_3":"FL-43_3"), na.rm = TRUE)/2,
        ssb.liking = rowMeans(select(., "FL-44_3":"FL-47_3"), na.rm = TRUE)/2,
        sfbl.liking = rowMeans(select(., "FL-39_3":"FL-47_3"), na.rm = TRUE)/2,
         unhealthyfat.liking = rowMeans(select(., "FL-48_3":"FL-52_3"), na.rm = TRUE)/2,
         protein.liking = rowMeans(select(., "FL-53_3":"FL-57_3"), na.rm = TRUE)/2, 
         spicy.liking = rowMeans(select(., "FL-58_3":"FL-60_3"), na.rm = TRUE)/2) %>%
  select(c("ID", ends_with(".liking"))) 

## Three factor questionnaire
TFQ.df <- survey.df %>% 
  select("ID", "TFQ-01":"TFQ-21") %>%
  mutate(uncontrolled = rowSums(select(., "TFQ-03", "TFQ-06", "TFQ-08", "TFQ-09", "TFQ-12", "TFQ-13", "TFQ-15", "TFQ-19", "TFQ-20"), na.rm = TRUE)) %>%  
  mutate(uncontrolled = ((uncontrolled-1)/36)*100, na.rm = TRUE) %>%
  mutate(restraint = rowSums(select(., "TFQ-01", "TFQ-05", "TFQ-11"), na.rm = TRUE)) %>%
  mutate(restraint = ((restraint-1)/12)*100, na.rm = TRUE) %>%
  mutate(emotional = rowSums(select(., "TFQ-02", "TFQ-04", "TFQ-07", "TFQ-10", "TFQ-14", "TFQ-16"), na.rm = TRUE)) %>%
  mutate(emotional = ((emotional-1)/24)*100, na.rm = TRUE) %>%
    select(c("ID", "uncontrolled", "restraint", "emotional")) 

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

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

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


## Combining AFR and ASI genotype files
AFR.geno <- AFR.geno.gwas.df %>%
  select(-c("CHR:POS", genetic.ancestry, phenotype.cat, effect_allele, other_allele, major_allele, minor_allele, effect, source)) %>%
  pivot_longer(cols = c(starts_with("1")), names_to = "ID", values_to = "genotype") %>%
  pivot_wider(names_from = rsID, values_from = genotype)

ASI.geno <- ASI.geno.gwas.df %>%
  select(-c("CHR:POS", genetic.ancestry, phenotype.cat, effect_allele, other_allele, major_allele, minor_allele, effect, source, "0_1001a-1-0424948735_1001a-1-0424948735")) %>%
  pivot_longer(cols = c(starts_with("1")), names_to = "ID", values_to = "genotype") %>%
  pivot_wider(names_from = rsID, values_from = genotype)

genotypes.df <- bind_rows(AFR.geno, ASI.geno)

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

survey.cleaned <- survey.df %>%
  select(c("ID", "sex", "age", "race", "ethn", "usb", "education", "hinc", "exec", "atus", "tius", "dcus", "scus")) 

full.predropout <- anthro.df %>%
  left_join(SST.df, by = "ID") %>%
  left_join(survey.cleaned, 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") %>% 
  left_join(genetic.ancestry, by = "ID") %>%
  left_join(genotypes.df, by = "ID") %>%
  mutate(ageCat = case_when(
      age >= 18 & age < 21 ~ "1",
      age >= 21 & age < 25 ~ "2",
      age >= 25 & age < 35  ~ "3",
      age >= 35 & age < 45  ~ "4",
      age >= 45 & age < 65  ~ "5",
      age >= 65 ~ "5")) %>%
  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)) %>%  
  mutate(atusCat = case_when(
      atus == 2 ~ "< 4 years",
      atus >= 3 & atus <= 4 ~ "4 to 12 years",
      atus == 5 ~ "13 to 18 years",
      atus >= 6 & atus <= 7 ~ "> 18 years")) %>%
  mutate(tiusCat = case_when(
      tius <= 4 ~ "< 5 years",
      tius == 5 ~ "5 to 10 years",
      tius >= 6 & tius <= 7 ~ "10 to 20 years", 
      tius == 8 ~ "> 20 years")) %>%  
    mutate(dcusCat = case_when(
      dcus <= 3 ~ "little to no change",
      dcus == 4 ~ "changed moderately",
      dcus >= 5 ~ "changed a lot")) 
  

full.df <- full.predropout %>%
  filter(!is.na(age) & !is.na(SST))

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(usb == "U.S. Born")
  
Immigrant.df <- full.df %>%
  filter(usb == "Immigrant")

Data distribution

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

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

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

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 `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.78883, p-value < 2.2e-16

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

`stat_bin()` using `bins = 30`. Pick better value `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.85032, p-value < 2.2e-16

## 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 `binwidth`.

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

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

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

    Shapiro-Wilk normality test

data:  full.df$SST
W = 0.99533, p-value = 0.2519

full.df %>%
  filter(!is.na(genetic.ancestry)) %>%
  ggplot(aes(SST)) +
  geom_density(aes(fill= genetic.ancestry), alpha = 0.4) + 
    xlim(-150, 180) +
    labs(title="Sweet Taste Liking by Ancestry Groups") +
  theme_classic()

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

`stat_bin()` using `bins = 30`. Pick better value `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.9836, p-value = 0.0001258

full.df %>%
  filter(!is.na(genetic.ancestry)) %>%
  ggplot(aes(sfbl.liking)) +
  geom_density(aes(fill= genetic.ancestry), alpha = 0.4) + 
    xlim(-150, 180) +
    labs(title="Sweet Food and Beverage Liking by Ancestry Groups") +
  theme_classic()

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

`stat_bin()` using `bins = 30`. Pick better value `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.90818, p-value = 4.078e-15

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

`stat_bin()` using `bins = 30`. Pick better value `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.99236, p-value = 0.03296

Demographics and counts

## Counts - by ancestry and USB

full.df %>%
  filter(race %in% c("African", "Asian")) %>%
  group_by(race, usb) %>%
  count()

# A tibble: 4 × 3
# Groups:   race, usb [4]
  race    usb           n
  <chr>   <chr>     <int>
1 African Immigrant    86
2 African U.S. Born   116
3 Asian   Immigrant    95
4 Asian   U.S. Born   112

full.df %>%
  filter(race %in% c("African", "Asian")) %>%
  group_by(race, atus) %>%
  count()

# A tibble: 14 × 3
# Groups:   race, atus [14]
   race     atus     n
   <chr>   <dbl> <int>
 1 African     2    11
 2 African     3    12
 3 African     4    14
 4 African     5    28
 5 African     6    15
 6 African     7     6
 7 African    NA   116
 8 Asian       2    12
 9 Asian       3    13
10 Asian       4    11
11 Asian       5    27
12 Asian       6    26
13 Asian       7     6
14 Asian      NA   112

full.df %>%
  filter(race %in% c("African", "Asian")) %>%
  group_by(race, tius) %>%
  count()

# A tibble: 16 × 3
# Groups:   race, tius [16]
   race     tius     n
   <chr>   <dbl> <int>
 1 African     2     6
 2 African     3     7
 3 African     4    13
 4 African     5    20
 5 African     6    12
 6 African     7    16
 7 African     8    12
 8 African    NA   116
 9 Asian       2    12
10 Asian       3     8
11 Asian       4    10
12 Asian       5    22
13 Asian       6    14
14 Asian       7    18
15 Asian       8    11
16 Asian      NA   112

full.df %>%
  filter(race %in% c("African", "Asian")) %>%
  group_by(race, dcus) %>%
  count()

# A tibble: 12 × 3
# Groups:   race, dcus [12]
   race     dcus     n
   <chr>   <dbl> <int>
 1 African     2     8
 2 African     3    21
 3 African     4    18
 4 African     7    28
 5 African     8    10
 6 African    NA   117
 7 Asian       2    14
 8 Asian       3    20
 9 Asian       4    36
10 Asian       7    23
11 Asian       8     2
12 Asian      NA   112

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

# A tibble: 14 × 3
# Groups:   race, scus [14]
   race     scus     n
   <chr>   <dbl> <int>
 1 African     2     4
 2 African     3     3
 3 African     4    14
 4 African     5    35
 5 African     6    29
 6 African    NA   117
 7 Asian       2     3
 8 Asian       3    11
 9 Asian       4    21
10 Asian       5    48
11 Asian       6    12
12 Asian      NA   112
13 <NA>        6     1
14 <NA>       NA     3

## Graph, time living in US
      full.df %>%
        group_by(race, usb, tius) %>%
        filter(!is.na(tius) & !is.na(usb) & !is.na(race)) %>% 
        count() %>%
        group_by(race, usb) %>%
        filter(race %in% c("African", "Asian")) %>%
        mutate(percent = (n / sum(n)) * 100) %>%
        ggplot(aes(x = usb, y = percent, fill = as.factor(tius))) +
        geom_bar(stat = "identity", position = position_dodge()) +
        labs(title = "Duration living in the US (percentages)", y = "Percentage", fill = "ageCat") +
        scale_fill_discrete(labels = c("< 1 year", "1 to 2 years", "3 to 5 years", "5 to 10 years", "10 to 15 years", "15 to 20 years", "> 20 years")) +
        theme_classic() +
        facet_wrap(~race) +
        theme(strip.text.x = element_text(size = 12, face = "bold")) 

## Graph, age when moved to the US
      full.df %>%
        group_by(race, usb, atus) %>%
        filter(!is.na(atus) & !is.na(usb) & !is.na(race)) %>% 
        count() %>%
        group_by(race, usb) %>%
        filter(race %in% c("African", "Asian")) %>%
        mutate(percent = (n / sum(n)) * 100) %>%
        ggplot(aes(x = usb, y = percent, fill = as.factor(atus))) +
        geom_bar(stat = "identity", position = position_dodge()) +
        labs(title = "Age when moved to the US (percentages)", y = "Percentage", fill = "ageCat") +
        scale_fill_discrete(labels = c("< 1 year", "1 to 2 years", "3 to 5 years", "5 to 10 years", "10 to 15 years", "15 to 20 years", "> 20 years")) +
        theme_classic() +
        facet_wrap(~race) +
        theme(strip.text.x = element_text(size = 12, face = "bold")) 

## Graph, diet change after moving to the US
      full.df %>%
        group_by(race, usb, dcus) %>%
        filter(!is.na(dcus) & !is.na(usb) & !is.na(race)) %>% 
        count() %>%
        group_by(race, usb) %>%
        filter(race %in% c("African", "Asian")) %>%
        mutate(percent = (n / sum(n)) * 100) %>%
        ggplot(aes(x = usb, y = percent, fill = as.factor(dcus))) +
        geom_bar(stat = "identity", position = position_dodge()) +
        labs(title = "Diet change after moving to the US (percentages)", y = "Percentage", fill = "ageCat") +
        scale_fill_discrete(labels = c("No change", "Changed a little", "Changed moderately", "Changed a lot", "Changed completely")) +
        theme_classic() +
        facet_wrap(~race) +
        theme(strip.text.x = element_text(size = 12, face = "bold"))       

Correlations - Corrplot

adjust.cor.pmat <- function(p.mat, method = "fdr") {
  p.mat[upper.tri(p.mat)] <- p.adjust(p.mat[upper.tri(p.mat)], method = method)
  p.mat[lower.tri(p.mat)] <- t(p.mat)[lower.tri(p.mat)]
  diag(p.mat) <- 0
  return(p.mat)
}

## All corr
corr.df <- full.df %>%
  select(age, sex, BMI, hinc, restraint, emotional, uncontrolled, sugar.intake, sHEI, sfbl.liking, SST) %>%
  rename("household income" = hinc, "sugar intake" = sugar.intake, 
         "healthy eating index" = sHEI, "sweet foods liking" = sfbl.liking, "sweet taste liking" = SST) 

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

p.mat <- cor_pmat(corr.df) %>%
  column_to_rownames()
p.mat.adj <- adjust.cor.pmat(p.mat, method = "fdr")

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

ggsave("All.cor.png", heigh = 5, width = 5)


## Corr - diet change
dc.corr.df <- full.df %>%
  filter(usb == "Immigrant") %>%
  select(age, sex, BMI, sHEI, sugar.intake, hinc, tius, atus, dcus, scus) %>%
  rename("healthy eating index" = sHEI, "household income" = hinc, "sugar intake" = sugar.intake, "duration in US" = tius, "age when moved to US" = atus, 
         "diet change after moving to US" = dcus, "sugar intake change after moving to US" = scus) 

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

p.mat <- cor_pmat(dc.corr.df) %>%
  column_to_rownames()

p.mat.adj <- adjust.cor.pmat(p.mat, method = "fdr")

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

## African diet change corr
afr.dc.corr.df <- full.df %>%
  filter(usb == "Immigrant", race == "African") %>%
  select(age, sex, BMI, sHEI, sugar.intake, hinc, tius, atus, dcus, scus) %>%
  rename("healthy eating index" = sHEI, "household income" = hinc, "sugar intake" = sugar.intake, "duration in US" = tius, "age when moved to US" = atus, "diet change after moving to US" = dcus, "sugar intake change after moving to US" = scus) 

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

p.mat <- cor_pmat(afr.dc.corr.df) %>%
  column_to_rownames()

p.mat.adj <- adjust.cor.pmat(p.mat, method = "fdr")

ggcorrplot(afr.dc.corr, type = "lower", p.mat = p.mat.adj, insig = "blank", lab = TRUE, lab_size = 2.5, colors = c("#E46726", "white", "#6D9EC1")) + ggtitle("African")

## Asian corr
asi.dc.corr.df <- full.df %>%
  filter(usb == "Immigrant", race == "Asian") %>%
  select(age, sex, BMI, sHEI, sugar.intake, hinc, tius, atus, dcus, scus) %>%
  rename("healthy eating index" = sHEI, "household income" = hinc, "sugar intake" = sugar.intake, "duration in US" = tius, "age when moved to US" = atus, "diet change after moving to US" = dcus, "sugar intake change after moving to US" = scus) 

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

p.mat <- cor_pmat(asi.dc.corr.df) %>%
  column_to_rownames()

p.mat.adj <- adjust.cor.pmat(p.mat, method = "fdr")

ggcorrplot(asi.dc.corr, type = "lower", p.mat = p.mat.adj, insig = "blank", lab = TRUE, lab_size = 2.5, colors = c("#E46726", "white", "#6D9EC1")) + ggtitle("Asian")

Differences - BMI

## BMI, ANCOVA, African vs Asian
bmi.ancova <- full.df %>%
  aov(BMI ~ race + ageCat, data = .)
summary(bmi.ancova)

             Df Sum Sq Mean Sq F value   Pr(>F)    
race          1   2609    2609  66.924 3.73e-15 ***
ageCat        4    740     185   4.746 0.000943 ***
Residuals   403  15710      39                     
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
4 observations deleted due to missingness

emmeans(bmi.ancova, pairwise ~ race, adjust = "tukey")

$emmeans
 race    emmean    SE  df lower.CL upper.CL
 African   29.5 0.622 403     28.3     30.7
 Asian     24.7 0.623 403     23.4     25.9

Results are averaged over the levels of: ageCat 
Confidence level used: 0.95 

$contrasts
 contrast        estimate    SE  df t.ratio p.value
 African - Asian     4.83 0.622 403   7.776  <.0001

Results are averaged over the levels of: ageCat 

emmeans(bmi.ancova, ~race)

 race    emmean    SE  df lower.CL upper.CL
 African   29.5 0.622 403     28.3     30.7
 Asian     24.7 0.623 403     23.4     25.9

Results are averaged over the levels of: ageCat 
Confidence level used: 0.95 

## BMI, ANCOVA, African, USB vs Immigrant
bmi.afr.ancova <- full.df %>%
  filter(race == "African") %>%
  aov(BMI ~ usb + ageCat, data = .)
summary(bmi.afr.ancova)

             Df Sum Sq Mean Sq F value  Pr(>F)   
usb           1     16   15.60   0.246 0.62022   
ageCat        4    937  234.34   3.701 0.00627 **
Residuals   196  12410   63.32                   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

emmeans(bmi.afr.ancova, ~usb)

 usb       emmean   SE  df lower.CL upper.CL
 Immigrant   29.6 1.15 196     27.4     31.9
 U.S. Born   30.6 1.14 196     28.4     32.8

Results are averaged over the levels of: ageCat 
Confidence level used: 0.95 

## BMI, ANCOVA, African, immigrant, atus
BMI.afr.atus.ancova <- full.df %>%
  filter(race == "African", usb == "Immigrant") %>%
  filter(!is.na(BMI) & !is.na(atusCat) & !is.na(ageCat)) %>%
  aov(BMI ~ atusCat + ageCat, data = .)
summary(BMI.afr.atus.ancova)

            Df Sum Sq Mean Sq F value Pr(>F)
atusCat      3    286   95.28   1.578  0.201
ageCat       4     74   18.41   0.305  0.874
Residuals   78   4709   60.37               

emmeans(BMI.afr.atus.ancova, ~atusCat)

 atusCat        emmean   SE df lower.CL upper.CL
 < 4 years        28.8 2.72 78     23.4     34.2
 > 18 years       31.4 1.80 78     27.8     35.0
 13 to 18 years   28.7 1.94 78     24.8     32.6
 4 to 12 years    27.6 2.05 78     23.5     31.7

Results are averaged over the levels of: ageCat 
Confidence level used: 0.95 

## BMI, ANCOVA, African, tius
BMI.afr.tius.ancova <- full.df %>%
  filter(race == "African", usb == "Immigrant") %>%
  filter(!is.na(BMI) & !is.na(tiusCat) & !is.na(ageCat)) %>%
  aov(BMI ~ tiusCat + ageCat, data = .)
summary(BMI.afr.tius.ancova)

            Df Sum Sq Mean Sq F value Pr(>F)
tiusCat      3    318  105.96   1.796  0.155
ageCat       4    149   37.21   0.631  0.642
Residuals   78   4602   59.00               

emmeans(BMI.afr.tius.ancova, ~tiusCat)

 tiusCat        emmean   SE df lower.CL upper.CL
 < 5 years        30.7 1.93 78     26.9     34.6
 > 20 years       31.6 2.63 78     26.4     36.9
 10 to 20 years   26.9 1.94 78     23.1     30.8
 5 to 10 years    28.6 2.23 78     24.1     33.0

Results are averaged over the levels of: ageCat 
Confidence level used: 0.95 

## BMI, ANCOVA, African, dcus
BMI.afr.dcus.ancova <- full.df %>%
  filter(race == "African", usb == "Immigrant") %>%
    filter(!is.na(BMI) & !is.na(dcusCat) & !is.na(ageCat)) %>%
  aov(BMI ~ dcusCat + ageCat, data = .)
summary(BMI.afr.dcus.ancova)

            Df Sum Sq Mean Sq F value Pr(>F)
dcusCat      2      4    1.80   0.029  0.971
ageCat       4    213   53.22   0.865  0.489
Residuals   78   4798   61.52               

emmeans(BMI.afr.dcus.ancova, pairwise ~ dcusCat, adjust = "tukey")

$emmeans
 dcusCat             emmean   SE df lower.CL upper.CL
 changed a lot         30.0 1.69 78     26.6     33.3
 changed moderately    30.1 2.11 78     25.9     34.3
 little to no change   29.2 1.80 78     25.6     32.8

Results are averaged over the levels of: ageCat 
Confidence level used: 0.95 

$contrasts
 contrast                                 estimate   SE df t.ratio p.value
 changed a lot - changed moderately         -0.155 2.29 78  -0.068  0.9975
 changed a lot - little to no change         0.775 2.08 78   0.374  0.9260
 changed moderately - little to no change    0.931 2.44 78   0.381  0.9232

Results are averaged over the levels of: ageCat 
P value adjustment: tukey method for comparing a family of 3 estimates 

emmeans(BMI.afr.dcus.ancova, ~dcusCat)

 dcusCat             emmean   SE df lower.CL upper.CL
 changed a lot         30.0 1.69 78     26.6     33.3
 changed moderately    30.1 2.11 78     25.9     34.3
 little to no change   29.2 1.80 78     25.6     32.8

Results are averaged over the levels of: ageCat 
Confidence level used: 0.95 

## BMI, ANCOVA, Asian, USB vs Immigrant
bmi.asi.ancova <- full.df %>%
  filter(race == "Asian") %>%
  aov(BMI ~ usb + ageCat, data = .)
summary(bmi.asi.ancova)

             Df Sum Sq Mean Sq F value Pr(>F)
usb           1   20.1   20.08   1.348  0.247
ageCat        4   72.8   18.19   1.221  0.303
Residuals   201 2994.4   14.90               

emmeans(bmi.asi.ancova, ~usb)

 usb       emmean    SE  df lower.CL upper.CL
 Immigrant   24.1 0.528 201     23.0     25.1
 U.S. Born   23.7 0.577 201     22.6     24.9

Results are averaged over the levels of: ageCat 
Confidence level used: 0.95 

## BMI, ANCOVA, Asian, atus
BMI.asi.atus.ancova <- full.df %>%
  filter(race == "Asian", usb == "Immigrant") %>%
  filter(!is.na(BMI) & !is.na(atus) & !is.na(ageCat)) %>%
  aov(BMI ~ atusCat + ageCat, data = .)
summary(BMI.asi.atus.ancova)

            Df Sum Sq Mean Sq F value Pr(>F)
atusCat      3   12.6    4.19   0.219  0.883
ageCat       4   95.6   23.90   1.252  0.295
Residuals   87 1661.4   19.10               

emmeans(BMI.asi.atus.ancova, ~atusCat)

 atusCat        emmean    SE df lower.CL upper.CL
 < 4 years        25.4 1.470 87     22.5     28.3
 > 18 years       23.5 0.846 87     21.9     25.2
 13 to 18 years   24.6 1.030 87     22.5     26.6
 4 to 12 years    25.3 1.130 87     23.0     27.5

Results are averaged over the levels of: ageCat 
Confidence level used: 0.95 

## BMI, ANCOVA, Asian, tius
BMI.asi.tius.ancova <- full.df %>%
  filter(race == "Asian", usb == "Immigrant") %>%
  filter(!is.na(BMI) & !is.na(tiusCat) & !is.na(ageCat)) %>%
  aov(BMI ~ tiusCat + ageCat, data = .)
summary(BMI.asi.tius.ancova)

            Df Sum Sq Mean Sq F value Pr(>F)
tiusCat      3   83.2   27.74   1.450  0.234
ageCat       4   21.6    5.41   0.283  0.888
Residuals   87 1664.7   19.13               

emmeans(BMI.asi.tius.ancova, ~tiusCat)

 tiusCat        emmean   SE df lower.CL upper.CL
 < 5 years        23.6 1.15 87     21.3     25.9
 > 20 years       25.8 1.50 87     22.8     28.8
 10 to 20 years   23.6 1.16 87     21.3     25.9
 5 to 10 years    22.7 1.33 87     20.0     25.3

Results are averaged over the levels of: ageCat 
Confidence level used: 0.95 

## BMI, ANCOVA, Asian, dcus
BMI.asi.dcus.ancova <- full.df %>%
  filter(race == "Asian", usb == "Immigrant") %>%
    filter(!is.na(BMI) & !is.na(dcusCat) & !is.na(ageCat)) %>%
  aov(BMI ~ dcusCat + ageCat, data = .)
summary(BMI.asi.dcus.ancova)

            Df Sum Sq Mean Sq F value Pr(>F)
dcusCat      2   13.1   6.561   0.343  0.711
ageCat       4   70.8  17.708   0.924  0.454
Residuals   88 1685.6  19.154               

emmeans(BMI.asi.dcus.ancova, ~dcusCat)

 dcusCat             emmean    SE df lower.CL upper.CL
 changed a lot         23.9 0.968 88     22.0     25.8
 changed moderately    24.8 0.903 88     23.0     26.6
 little to no change   23.9 0.904 88     22.1     25.7

Results are averaged over the levels of: ageCat 
Confidence level used: 0.95 

Differences - Healthy Eating Index

## HEI, ANCOVA, African vs Asian
HEI.ancova <- full.df %>%
  aov(sHEI ~ race + ageCat, data = .)
summary(HEI.ancova)

             Df Sum Sq Mean Sq F value Pr(>F)  
race          1     28   27.90   0.320 0.5719  
ageCat        4    775  193.82   2.224 0.0657 .
Residuals   403  35122   87.15                 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
4 observations deleted due to missingness

emmeans(HEI.ancova, pairwise ~ race, adjust = "tukey")

$emmeans
 race    emmean    SE  df lower.CL upper.CL
 African   46.7 0.930 403     44.9     48.6
 Asian     46.4 0.932 403     44.5     48.2

Results are averaged over the levels of: ageCat 
Confidence level used: 0.95 

$contrasts
 contrast        estimate   SE  df t.ratio p.value
 African - Asian    0.375 0.93 403   0.403  0.6868

Results are averaged over the levels of: ageCat 

emmeans(HEI.ancova, ~race)

 race    emmean    SE  df lower.CL upper.CL
 African   46.7 0.930 403     44.9     48.6
 Asian     46.4 0.932 403     44.5     48.2

Results are averaged over the levels of: ageCat 
Confidence level used: 0.95 

## HEI, ANCOVA, African, usb vs immigrant
HEI.afr.ancova <- full.df %>%
  filter(!is.na(sHEI) & !is.na(usb) & !is.na(ageCat)) %>%
  filter(race == "African") %>%
  aov(sHEI ~ usb + ageCat + sex, data = .)
summary(HEI.afr.ancova)

             Df Sum Sq Mean Sq F value   Pr(>F)    
usb           1     12    11.9   0.133   0.7155    
ageCat        4    760   190.0   2.128   0.0788 .  
sex           1   2075  2074.8  23.237 2.88e-06 ***
Residuals   195  17412    89.3                     
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

emmeans(HEI.afr.ancova, ~usb)

 usb       emmean   SE  df lower.CL upper.CL
 Immigrant   45.8 1.37 195     43.1     48.5
 U.S. Born   47.5 1.35 195     44.8     50.1

Results are averaged over the levels of: ageCat 
Confidence level used: 0.95 

## HEI, ANCOVA, African, atus
HEI.afr.atus.ancova <- full.df %>%
  filter(race == "African", usb == "Immigrant") %>%
  filter(!is.na(sHEI) & !is.na(atusCat) & !is.na(ageCat)) %>%
  aov(sHEI ~ atusCat + ageCat + sex, data = .)
summary(HEI.afr.atus.ancova)

            Df Sum Sq Mean Sq F value Pr(>F)
atusCat      3    104   34.60   0.332  0.802
ageCat       4    170   42.38   0.406  0.804
sex          1      6    5.60   0.054  0.817
Residuals   77   8030  104.28               

emmeans(HEI.afr.atus.ancova, ~atusCat)

 atusCat        emmean   SE df lower.CL upper.CL
 < 4 years        45.4 3.76 77     37.9     52.9
 > 18 years       43.0 2.74 77     37.6     48.5
 13 to 18 years   46.7 2.72 77     41.2     52.1
 4 to 12 years    43.6 2.91 77     37.8     49.4

Results are averaged over the levels of: ageCat, sex 
Confidence level used: 0.95 

## HEI, ANCOVA, African, tius
HEI.afr.tius.ancova <- full.df %>%
  filter(race == "African", usb == "Immigrant") %>%
  filter(!is.na(sHEI) & !is.na(tiusCat) & !is.na(ageCat)) %>%
  aov(sHEI ~ tiusCat + ageCat + sex, data = .)
summary(HEI.afr.tius.ancova)

            Df Sum Sq Mean Sq F value Pr(>F)
tiusCat      3    280   93.37   0.911  0.440
ageCat       4    138   34.39   0.336  0.853
sex          1      1    0.81   0.008  0.929
Residuals   77   7890  102.47               

emmeans(HEI.afr.tius.ancova, ~tiusCat)

 tiusCat        emmean   SE df lower.CL upper.CL
 < 5 years        46.9 2.91 77     41.1     52.7
 > 20 years       43.2 3.63 77     35.9     50.4
 10 to 20 years   45.5 2.80 77     40.0     51.1
 5 to 10 years    41.7 3.12 77     35.5     47.9

Results are averaged over the levels of: ageCat, sex 
Confidence level used: 0.95 

## HEI, ANCOVA, African, dcus
HEI.afr.dcus.ancova <- full.df %>%
  filter(race == "African", usb == "Immigrant") %>%
    filter(!is.na(sHEI) & !is.na(dcusCat) & !is.na(ageCat)) %>%
  aov(sHEI ~ dcusCat + ageCat + sex, data = .)
summary(HEI.afr.dcus.ancova)

            Df Sum Sq Mean Sq F value Pr(>F)
dcusCat      2    381  190.36   1.866  0.162
ageCat       4     14    3.54   0.035  0.998
sex          1     11   11.02   0.108  0.743
Residuals   77   7857  102.04               

emmeans(HEI.afr.dcus.ancova, ~dcusCat)

 dcusCat             emmean   SE df lower.CL upper.CL
 changed a lot         41.9 2.55 77     36.8     47.0
 changed moderately    45.6 2.90 77     39.8     51.4
 little to no change   46.2 2.57 77     41.1     51.3

Results are averaged over the levels of: ageCat, sex 
Confidence level used: 0.95 

## HEI, ANCOVA, Asian, usb vs immigrant
HEI.asi.ancova <- full.df %>%
  filter(race == "Asian") %>%
  aov(sHEI ~ usb + ageCat + sex, data = .)
summary(HEI.asi.ancova)

             Df Sum Sq Mean Sq F value   Pr(>F)    
usb           1    247   246.6   3.534 0.061590 .  
ageCat        4    545   136.2   1.952 0.103344    
sex           1    889   889.5  12.745 0.000447 ***
Residuals   200  13958    69.8                     
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

emmeans(HEI.asi.ancova, ~usb) 

 usb       emmean   SE  df lower.CL upper.CL
 Immigrant   47.1 1.15 200     44.9     49.4
 U.S. Born   45.4 1.25 200     43.0     47.9

Results are averaged over the levels of: ageCat 
Confidence level used: 0.95 

## Graph here? it is borderline significant

## HEI, ANCOVA, Asian, atus
HEI.asi.atus.ancova <- full.df %>%
  filter(race == "Asian", usb == "Immigrant") %>%
  filter(!is.na(sHEI) & !is.na(atusCat) & !is.na(ageCat)) %>%
  aov(sHEI ~ atusCat + ageCat + sex, data = .)
summary(HEI.asi.atus.ancova)

            Df Sum Sq Mean Sq F value Pr(>F)  
atusCat      3     95   31.63   0.482 0.6954  
ageCat       4    705  176.24   2.688 0.0365 *
sex          1     32   31.69   0.483 0.4888  
Residuals   86   5639   65.57                 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

emmeans(HEI.asi.atus.ancova, ~atusCat)

 atusCat        emmean   SE df lower.CL upper.CL
 < 4 years        49.6 2.75 86     44.2     55.1
 > 18 years       47.0 1.57 86     43.9     50.1
 13 to 18 years   49.0 1.90 86     45.2     52.8
 4 to 12 years    47.0 2.10 86     42.9     51.2

Results are averaged over the levels of: ageCat 
Confidence level used: 0.95 

## ANCOVA, Asian, tius
HEI.asi.tius.ancova <- full.df %>%
  filter(race == "Asian") %>%
  filter(!is.na(sHEI) & !is.na(tiusCat) & !is.na(ageCat)) %>%
  aov(sHEI ~ tiusCat + ageCat + sex, data = .)
summary(HEI.asi.tius.ancova)

            Df Sum Sq Mean Sq F value Pr(>F)  
tiusCat      3    236   78.71   1.247 0.2976  
ageCat       4    776  194.09   3.076 0.0203 *
sex          1     32   32.26   0.511 0.4765  
Residuals   86   5426   63.09                 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

emmeans(HEI.asi.tius.ancova, ~tiusCat)

 tiusCat        emmean   SE df lower.CL upper.CL
 < 5 years        48.1 2.09 86     44.0     52.3
 > 20 years       44.4 2.74 86     39.0     49.9
 10 to 20 years   49.3 2.13 86     45.0     53.5
 5 to 10 years    52.3 2.41 86     47.5     57.1

Results are averaged over the levels of: ageCat 
Confidence level used: 0.95 

## ANCOVA, Asian, dcus
HEI.asi.dcus.ancova <- full.df %>%
  filter(race == "Asian", usb == "Immigrant") %>%
  filter(!is.na(sHEI) & !is.na(dcusCat) & !is.na(ageCat)) %>%
  aov(sHEI ~ dcusCat + ageCat + sex, data = .)
summary(HEI.asi.dcus.ancova)

            Df Sum Sq Mean Sq F value Pr(>F)  
dcusCat      2     44   21.83   0.336 0.7155  
ageCat       4    746  186.40   2.870 0.0276 *
sex          1     31   31.13   0.479 0.4905  
Residuals   87   5650   64.95                 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

emmeans(HEI.asi.dcus.ancova, ~dcusCat)

 dcusCat             emmean   SE df lower.CL upper.CL
 changed a lot         46.2 1.79 87     42.7     49.8
 changed moderately    48.7 1.67 87     45.4     52.1
 little to no change   48.0 1.67 87     44.7     51.3

Results are averaged over the levels of: ageCat 
Confidence level used: 0.95 

Differences - sugar intake

## sugar, ANCOVA, African vs Asian
sugar.ancova <- full.df %>%
  aov(sugar.intake ~ race + ageCat, data = .)
summary(sugar.ancova)

             Df   Sum Sq Mean Sq F value   Pr(>F)    
race          1  2578198 2578198  42.340 2.28e-10 ***
ageCat        4   468211  117053   1.922    0.106    
Residuals   403 24539881   60893                     
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
4 observations deleted due to missingness

emmeans(sugar.ancova, pairwise ~ race, adjust = "tukey")

$emmeans
 race    emmean   SE  df lower.CL upper.CL
 African    528 24.6 403      480      576
 Asian      365 24.6 403      316      413

Results are averaged over the levels of: ageCat 
Confidence level used: 0.95 

$contrasts
 contrast        estimate   SE  df t.ratio p.value
 African - Asian      163 24.6 403   6.638  <.0001

Results are averaged over the levels of: ageCat 

emmeans.sugar.race <- emmeans(sugar.ancova, ~race) %>% as.data.frame()

 p.values <- pairs(emmeans(sugar.ancova, ~race), adjust = "tukey") %>%
    as.data.frame() %>%
    mutate(group1 = sub(" - .*", "", contrast), group2 = sub(".* - ", "", contrast),
    group1 = gsub(" ", " ", group1), group2 = gsub(" ", " ", group2),
    p.adj = p.value,
    p.adj.signif = case_when(p.adj < 0.001 ~ "***", p.adj < 0.01 ~ "**", p.adj < 0.05 ~ "*", TRUE ~ "ns"))

## Graph
emmeans.sugar.race %>%
    filter(race %in% c("African", "Asian")) %>%
  ggplot(aes(x = race, y = emmean, fill = race)) +
  geom_bar(stat = "identity", position = position_dodge(), color = "black") +
  geom_errorbar(aes(ymin = emmean-SE, ymax = emmean+SE), 
                width = 0.2, position = position_dodge(0.9)) +
    stat_pvalue_manual(p.values, hide.ns = T, label = "p.adj.signif", label.size = 5, bracket.size = 0.5, tip.length = 0.05, step.increase = 0.2, y.position = 650) +
  labs(x = " ", 
       y = "Sugar Intake in Calories (Adjusted Mean)",
       title = " ",
       fill = "Ancestry") +
  ylim(0, 800) +
  theme_classic() +
        theme(axis.text.x = element_text(angle = 45, hjust = 1, size = 14), axis.text.y = element_text(size = 14), axis.title.y = element_text(size = 16), plot.title = element_text(size = 18), legend.text = element_text(size = 14), legend.title = element_text(size = 16))  

## sugar, ANCOVA, African, usb vs immigrant
sugar.afr.usb.ancova <- full.df %>%
  filter(!is.na(sHEI) & !is.na(usb) & !is.na(ageCat)) %>%
  filter(race == "African") %>%
  aov(sugar.intake ~ usb + ageCat, data = .)
summary(sugar.afr.usb.ancova)

             Df   Sum Sq Mean Sq F value Pr(>F)
usb           1      680     680   0.009  0.926
ageCat        4   548000  137000   1.725  0.146
Residuals   196 15569262   79435               

emmeans(sugar.afr.usb.ancova, ~usb)

 usb       emmean   SE  df lower.CL upper.CL
 Immigrant    504 40.7 196      424      584
 U.S. Born    498 40.3 196      418      577

Results are averaged over the levels of: ageCat 
Confidence level used: 0.95 

## sugar, ANCOVA, African, immigrant, atus
sugar.afr.atus.ancova <- full.df %>%
  filter(race == "African", usb == "Immigrant") %>%
  filter(!is.na(sugar.intake) & !is.na(atusCat) & !is.na(ageCat)) %>%
  aov(sugar.intake ~ atusCat + ageCat, data = .)
summary(sugar.afr.atus.ancova)

            Df  Sum Sq Mean Sq F value Pr(>F)
atusCat      3  150787   50262   0.652  0.584
ageCat       4  279441   69860   0.907  0.464
Residuals   78 6008626   77034               

emmeans(sugar.afr.atus.ancova, ~atusCat)

 atusCat        emmean   SE df lower.CL upper.CL
 < 4 years         531 97.2 78      338      725
 > 18 years        511 64.2 78      383      639
 13 to 18 years    480 69.3 78      342      618
 4 to 12 years     561 73.1 78      415      706

Results are averaged over the levels of: ageCat 
Confidence level used: 0.95 

## sugar, ANCOVA, African, tius
sugar.afr.tius.ancova <- full.df %>%
  filter(race == "African", usb == "Immigrant") %>%
  filter(!is.na(sugar.intake) & !is.na(tiusCat) & !is.na(ageCat)) %>%
  aov(sugar.intake ~ tiusCat + ageCat, data = .)
summary(sugar.afr.tius.ancova)

            Df  Sum Sq Mean Sq F value Pr(>F)
tiusCat      3   59689   19896   0.257  0.856
ageCat       4  333720   83430   1.076  0.374
Residuals   78 6045444   77506               

emmeans(sugar.afr.tius.ancova, ~tiusCat)

 tiusCat        emmean   SE df lower.CL upper.CL
 < 5 years         479 70.0 78      339      618
 > 20 years        553 95.2 78      363      742
 10 to 20 years    526 70.5 78      386      666
 5 to 10 years     493 80.8 78      332      654

Results are averaged over the levels of: ageCat 
Confidence level used: 0.95 

## sugar, ANCOVA, African, dcus
sugar.afr.dcus.ancova <- full.df %>%
  filter(race == "African", usb == "Immigrant") %>%
    filter(!is.na(sugar.intake) & !is.na(dcusCat) & !is.na(ageCat)) %>%
  aov(sugar.intake ~ dcusCat + ageCat, data = .)
summary(sugar.afr.dcus.ancova)

            Df  Sum Sq Mean Sq F value Pr(>F)  
dcusCat      2  497585  248792   3.584 0.0324 *
ageCat       4  182537   45634   0.657 0.6235  
Residuals   78 5414981   69423                 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

emmeans(sugar.afr.dcus.ancova, pairwise ~ dcusCat, adjust = "tukey")

$emmeans
 dcusCat             emmean   SE df lower.CL upper.CL
 changed a lot          573 56.8 78      460      686
 changed moderately     387 70.9 78      245      528
 little to no change    516 60.5 78      396      637

Results are averaged over the levels of: ageCat 
Confidence level used: 0.95 

$contrasts
 contrast                                 estimate   SE df t.ratio p.value
 changed a lot - changed moderately          186.1 76.9 78   2.419  0.0465
 changed a lot - little to no change          56.5 69.7 78   0.810  0.6979
 changed moderately - little to no change   -129.6 82.1 78  -1.579  0.2606

Results are averaged over the levels of: ageCat 
P value adjustment: tukey method for comparing a family of 3 estimates 

emmeans(sugar.afr.dcus.ancova, ~dcusCat)

 dcusCat             emmean   SE df lower.CL upper.CL
 changed a lot          573 56.8 78      460      686
 changed moderately     387 70.9 78      245      528
 little to no change    516 60.5 78      396      637

Results are averaged over the levels of: ageCat 
Confidence level used: 0.95 

## sugar, ANCOVA, Asian, usb vs immigrant
sugar.asi.ancova <- full.df %>%
  filter(!is.na(sugar.intake) & !is.na(usb) & !is.na(ageCat)) %>%
  filter(race == "Asian") %>%
  aov(sugar.intake ~ usb + ageCat, data = .)
summary(sugar.asi.ancova)

             Df  Sum Sq Mean Sq F value Pr(>F)
usb           1   38035   38035   0.880  0.349
ageCat        4  161913   40478   0.936  0.444
Residuals   201 8690203   43235               

emmeans(sugar.asi.ancova, ~usb)

 usb       emmean   SE  df lower.CL upper.CL
 Immigrant    377 28.5 201      321      433
 U.S. Born    395 31.1 201      334      456

Results are averaged over the levels of: ageCat 
Confidence level used: 0.95 

## sugar, ANCOVA, Asian, atus
sugar.asi.atus.ancova <- full.df %>%
  filter(race == "Asian", usb == "Immigrant") %>%
  filter(!is.na(sugar.intake) & !is.na(atus) & !is.na(ageCat)) %>%
  aov(sugar.intake ~ atusCat + ageCat, data = .)
summary(sugar.asi.atus.ancova)

            Df  Sum Sq Mean Sq F value Pr(>F)
atusCat      3  161444   53815   1.704  0.172
ageCat       4  107845   26961   0.854  0.495
Residuals   87 2747692   31583               

emmeans(sugar.asi.atus.ancova, ~atusCat)

 atusCat        emmean   SE df lower.CL upper.CL
 < 4 years         414 59.7 87      295      532
 > 18 years        364 34.4 87      295      432
 13 to 18 years    347 41.8 87      264      430
 4 to 12 years     441 45.9 87      350      533

Results are averaged over the levels of: ageCat 
Confidence level used: 0.95 

## sugar, ANCOVA, Asian, tius
sugar.asi.tius.ancova <- full.df %>%
  filter(race == "Asian", usb == "Immigrant") %>%
  filter(!is.na(sugar.intake) & !is.na(tiusCat) & !is.na(ageCat)) %>%
  aov(sugar.intake ~ tiusCat + ageCat, data = .)
summary(sugar.asi.tius.ancova)

            Df  Sum Sq Mean Sq F value Pr(>F)
tiusCat      3  116533   38844   1.222  0.307
ageCat       4  134856   33714   1.061  0.381
Residuals   87 2765591   31788               

emmeans(sugar.asi.tius.ancova, ~tiusCat)

 tiusCat        emmean   SE df lower.CL upper.CL
 < 5 years         393 47.0 87      300      487
 > 20 years        309 61.3 87      187      430
 10 to 20 years    448 47.2 87      354      542
 5 to 10 years     390 54.1 87      282      497

Results are averaged over the levels of: ageCat 
Confidence level used: 0.95 

## sugar, ANCOVA, Asian, dcus
sugar.asi.dcus.ancova <- full.df %>%
  filter(race == "Asian", usb == "Immigrant") %>%
    filter(!is.na(sugar.intake) & !is.na(dcusCat) & !is.na(ageCat)) %>%
  aov(sugar.intake ~ dcusCat + ageCat, data = .)
summary(sugar.asi.dcus.ancova)

            Df  Sum Sq Mean Sq F value Pr(>F)
dcusCat      2   83680   41840   1.310  0.275
ageCat       4  122501   30625   0.959  0.434
Residuals   88 2810801   31941               

emmeans(sugar.asi.dcus.ancova, ~dcusCat)

 dcusCat             emmean   SE df lower.CL upper.CL
 changed a lot          396 39.5 88      318      475
 changed moderately     394 36.9 88      321      468
 little to no change    339 36.9 88      266      413

Results are averaged over the levels of: ageCat 
Confidence level used: 0.95 

Differences - sweet food and beverage liking

## sfbl, ANCOVA, African vs Asian
sfbl.ancova <- full.df %>%
  aov(sfbl.liking ~ race + ageCat, data = .)
summary(sfbl.ancova)

             Df Sum Sq Mean Sq F value   Pr(>F)    
race          1   2676  2675.7  16.659 5.40e-05 ***
ageCat        4   4570  1142.5   7.114 1.53e-05 ***
Residuals   403  64728   160.6                     
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
4 observations deleted due to missingness

emmeans(sfbl.ancova, pairwise ~ race, adjust = "tukey")

$emmeans
 race    emmean   SE  df lower.CL upper.CL
 African   18.6 1.26 403     16.1     21.0
 Asian     13.1 1.26 403     10.6     15.6

Results are averaged over the levels of: ageCat 
Confidence level used: 0.95 

$contrasts
 contrast        estimate   SE  df t.ratio p.value
 African - Asian     5.47 1.26 403   4.332  <.0001

Results are averaged over the levels of: ageCat 

emmeans.sfbl.race <- emmeans(sfbl.ancova, ~race) %>% as.data.frame()

## Graph
p.values <- pairs(emmeans(sfbl.ancova, ~race), adjust = "tukey") %>%
    as.data.frame() %>%
    mutate(group1 = sub(" - .*", "", contrast), group2 = sub(".* - ", "", contrast),
    group1 = gsub(" ", " ", group1), group2 = gsub(" ", " ", group2),
    p.adj = p.value,
    p.adj.signif = case_when(p.adj < 0.001 ~ "***", p.adj < 0.01 ~ "**", p.adj < 0.05 ~ "*", TRUE ~ "ns"))

  full.df %>%
  filter(race %in% c("African","Asian")) %>%
  ggplot(aes(x = race, y = sfbl.liking, fill = race)) +
  geom_boxplot(width = 0.5, alpha = 0.5, outlier.shape = NA) +
  geom_jitter(width = 0.2, alpha = 0.3, size = 1) + 
  geom_hline(yintercept = 0, linetype = "dashed", color = "red") +  
  stat_pvalue_manual(p.values, hide.ns = T, label = "p.adj.signif", label.size = 5, bracket.size = 0.5, tip.length = 0.05, step.increase = 0.1, y.position = 70) +
  labs(x = " ", 
       y = "Sweet Food and Beverage Liking (0-100)",
       title = " ",
       fill = "Genotype") +
  ylim(-50, 85) +
    scale_y_continuous(breaks = c(-50, -25, 0, 25, 50)) +
  theme_classic() +
        theme(axis.text.x = element_text(angle = 45, hjust = 1, size = 14), axis.text.y = element_text(size = 14), axis.title.y = element_text(size = 16), plot.title = element_text(size = 18), legend.text = element_text(size = 14), legend.title = element_text(size = 16))

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

## sfbl, ANCOVA, African, usb vs immigrant
SFBL.afr.ancova <- full.df %>%
  filter(!is.na(sHEI) & !is.na(usb) & !is.na(ageCat)) %>%
  filter(race == "African") %>%
  aov(sfbl.liking ~ usb + ageCat, data = .)
summary(SFBL.afr.ancova)

             Df Sum Sq Mean Sq F value Pr(>F)  
usb           1     97    96.8   0.524 0.4699  
ageCat        4   2330   582.5   3.156 0.0153 *
Residuals   196  36172   184.5                 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

emmeans(SFBL.afr.ancova, ~usb)

 usb       emmean   SE  df lower.CL upper.CL
 Immigrant   19.5 1.96 196     15.7     23.4
 U.S. Born   20.3 1.94 196     16.5     24.1

Results are averaged over the levels of: ageCat 
Confidence level used: 0.95 

## sfbl, ANCOVA, African, atus
SFBL.afr.atus.ancova <- full.df %>%
  filter(race == "African", usb == "Immigrant") %>%
  filter(!is.na(sfbl.liking) & !is.na(atusCat) & !is.na(ageCat)) %>%
  aov(sfbl.liking ~ atusCat + ageCat, data = .)
summary(SFBL.afr.atus.ancova)

            Df Sum Sq Mean Sq F value  Pr(>F)   
atusCat      3    183    61.1   0.363 0.77956   
ageCat       4   2488   622.1   3.701 0.00822 **
Residuals   78  13110   168.1                   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

emmeans(SFBL.afr.atus.ancova, ~atusCat)

 atusCat        emmean   SE df lower.CL upper.CL
 < 4 years        22.3 4.54 78     13.3     31.4
 > 18 years       19.6 3.00 78     13.6     25.6
 13 to 18 years   19.8 3.23 78     13.4     26.3
 4 to 12 years    21.7 3.42 78     14.9     28.5

Results are averaged over the levels of: ageCat 
Confidence level used: 0.95 

## sfbl, ANCOVA, African, tius
SFBL.afr.tius.ancova <- full.df %>%
  filter(race == "African", usb == "Immigrant") %>%
  filter(!is.na(sfbl.liking) & !is.na(tiusCat) & !is.na(ageCat)) %>%
  aov(sfbl.liking ~ tiusCat + ageCat, data = .)
summary(SFBL.afr.tius.ancova)

            Df Sum Sq Mean Sq F value Pr(>F)  
tiusCat      3   1013   337.8   2.037 0.1155  
ageCat       4   1835   458.7   2.766 0.0331 *
Residuals   78  12933   165.8                 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

emmeans(SFBL.afr.tius.ancova, ~tiusCat)

 tiusCat        emmean   SE df lower.CL upper.CL
 < 5 years        18.9 3.24 78     12.5     25.4
 > 20 years       20.8 4.40 78     12.0     29.6
 10 to 20 years   22.5 3.26 78     16.0     29.0
 5 to 10 years    18.3 3.74 78     10.8     25.7

Results are averaged over the levels of: ageCat 
Confidence level used: 0.95 

## sfbl, ANCOVA, African, dcus
SFBL.afr.dcus.ancova <- full.df %>%
  filter(race == "African", usb == "Immigrant") %>%
    filter(!is.na(sfbl.liking) & !is.na(dcusCat) & !is.na(ageCat)) %>%
  aov(sfbl.liking ~ dcusCat + ageCat, data = .)
summary(SFBL.afr.dcus.ancova)

            Df Sum Sq Mean Sq F value  Pr(>F)   
dcusCat      2    324   162.0   0.991 0.37574   
ageCat       4   2428   607.1   3.716 0.00805 **
Residuals   78  12744   163.4                   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

emmeans(SFBL.afr.dcus.ancova, ~dcusCat)

 dcusCat             emmean   SE df lower.CL upper.CL
 changed a lot         22.3 2.75 78     16.8     27.8
 changed moderately    19.7 3.44 78     12.9     26.6
 little to no change   18.1 2.93 78     12.3     24.0

Results are averaged over the levels of: ageCat 
Confidence level used: 0.95 

## sfbl, ANCOVA, Asian, usb vs immigrant
SFBL.asi.ancova <- full.df %>%
  filter(!is.na(sfbl.liking) & !is.na(usb) & !is.na(ageCat)) %>%
  filter(race == "Asian") %>%
  aov(sfbl.liking ~ usb + ageCat, data = .)
summary(SFBL.asi.ancova)

             Df Sum Sq Mean Sq F value  Pr(>F)   
usb           1    646   646.0   4.720 0.03099 * 
ageCat        4   2541   635.2   4.641 0.00132 **
Residuals   201  27513   136.9                   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

emmeans(SFBL.asi.ancova, pairwise ~ usb, adjust = "tukey")

$emmeans
 usb       emmean   SE  df lower.CL upper.CL
 Immigrant   11.2 1.60 201     8.05     14.4
 U.S. Born   13.5 1.75 201    10.02     16.9

Results are averaged over the levels of: ageCat 
Confidence level used: 0.95 

$contrasts
 contrast              estimate  SE  df t.ratio p.value
 Immigrant - U.S. Born    -2.26 1.7 201  -1.331  0.1846

Results are averaged over the levels of: ageCat 

emmeans(SFBL.asi.ancova, ~usb) 

 usb       emmean   SE  df lower.CL upper.CL
 Immigrant   11.2 1.60 201     8.05     14.4
 U.S. Born   13.5 1.75 201    10.02     16.9

Results are averaged over the levels of: ageCat 
Confidence level used: 0.95 

## sfbl, ANCOVA, Asian, atus
SFBL.asi.atus.ancova <- full.df %>%
  filter(race == "Asian", usb == "Immigrant") %>%
  filter(!is.na(sfbl.liking) & !is.na(atus) & !is.na(ageCat)) %>%
  aov(sfbl.liking ~ atusCat + ageCat, data = .)
summary(SFBL.asi.atus.ancova)

            Df Sum Sq Mean Sq F value Pr(>F)  
atusCat      3   1083   361.0   2.525 0.0628 .
ageCat       4   1077   269.3   1.884 0.1205  
Residuals   87  12440   143.0                 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

emmeans(SFBL.asi.atus.ancova, pairwise ~ atusCat, adjust = "tukey")

$emmeans
 atusCat        emmean   SE df lower.CL upper.CL
 < 4 years        18.7 4.01 87    10.69     26.6
 > 18 years       10.6 2.31 87     5.99     15.2
 13 to 18 years   10.1 2.81 87     4.49     15.7
 4 to 12 years    12.3 3.09 87     6.21     18.5

Results are averaged over the levels of: ageCat 
Confidence level used: 0.95 

$contrasts
 contrast                       estimate   SE df t.ratio p.value
 < 4 years - > 18 years            8.077 4.49 87   1.798  0.2813
 < 4 years - 13 to 18 years        8.588 4.20 87   2.046  0.1792
 < 4 years - 4 to 12 years         6.324 4.26 87   1.485  0.4510
 > 18 years - 13 to 18 years       0.511 3.44 87   0.148  0.9988
 > 18 years - 4 to 12 years       -1.753 3.60 87  -0.488  0.9617
 13 to 18 years - 4 to 12 years   -2.264 3.39 87  -0.668  0.9089

Results are averaged over the levels of: ageCat 
P value adjustment: tukey method for comparing a family of 4 estimates 

emmeans(SFBL.asi.atus.ancova, ~atusCat)

 atusCat        emmean   SE df lower.CL upper.CL
 < 4 years        18.7 4.01 87    10.69     26.6
 > 18 years       10.6 2.31 87     5.99     15.2
 13 to 18 years   10.1 2.81 87     4.49     15.7
 4 to 12 years    12.3 3.09 87     6.21     18.5

Results are averaged over the levels of: ageCat 
Confidence level used: 0.95 

## sfbl, ANCOVA, Asian, tius
SFBL.asi.tius.ancova <- full.df %>%
  filter(race == "Asian", usb == "Immigrant") %>%
  filter(!is.na(sfbl.liking) & !is.na(tiusCat) & !is.na(ageCat)) %>%
  aov(sfbl.liking ~ tiusCat + ageCat, data = .)
summary(SFBL.asi.tius.ancova)

            Df Sum Sq Mean Sq F value Pr(>F)
tiusCat      3    627   209.0   1.410  0.245
ageCat       4   1074   268.6   1.812  0.134
Residuals   87  12898   148.3               

emmeans(SFBL.asi.tius.ancova, ~tiusCat)

 tiusCat        emmean   SE df lower.CL upper.CL
 < 5 years        9.54 3.21 87     3.16     15.9
 > 20 years      12.48 4.18 87     4.16     20.8
 10 to 20 years  12.64 3.22 87     6.23     19.1
 5 to 10 years   10.01 3.70 87     2.66     17.4

Results are averaged over the levels of: ageCat 
Confidence level used: 0.95 

## sfbl, ANCOVA, Asian, dcus
SFBL.asi.dcus.ancova <- full.df %>%
  filter(race == "Asian", usb == "Immigrant") %>%
    filter(!is.na(sfbl.liking) & !is.na(dcusCat) & !is.na(ageCat)) %>%
  aov(sfbl.liking ~ dcusCat + ageCat, data = .)
summary(SFBL.asi.dcus.ancova)

            Df Sum Sq Mean Sq F value Pr(>F)  
dcusCat      2    726   363.2   2.565 0.0827 .
ageCat       4   1412   352.9   2.492 0.0488 *
Residuals   88  12462   141.6                 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

emmeans(SFBL.asi.dcus.ancova, ~dcusCat)

 dcusCat             emmean   SE df lower.CL upper.CL
 changed a lot         7.35 2.63 88     2.12     12.6
 changed moderately   13.00 2.46 88     8.12     17.9
 little to no change  13.40 2.46 88     8.51     18.3

Results are averaged over the levels of: ageCat 
Confidence level used: 0.95