rm(list=ls())
setwd("E:/Dropbox/Research/_DeLoach/Ag_Data_News/")

#############################################################
## Code to create plots in Ag Data News article            ##
##  Title:  Which Milk Substitute is the Best?      ##
##  Link:  https://asmith.ucdavis.edu/news/milk    ##
#############################################################

# If a package is installed, it will be loaded. If any are not, the missing package(s) will be installed 
# from CRAN and then loaded.
pacman::p_load(tidyverse,readxl,RColorBrewer,httr,ggpubr,rnassqs,data.table,readxl,ggimage,ggplot2,scales,lubridate,blscrapeR)

logos <- data.frame(name=c("Almond", "Coconut", "Dairy (NF)","Dairy (2%)","Dairy (Whole)", "Oat","Soy" ),
                       logo=c("https://files.asmith.ucdavis.edu/almond.jpg",
                              "https://files.asmith.ucdavis.edu/coconut.jpg",
                              "https://files.asmith.ucdavis.edu/cow.jpg",
                              "https://files.asmith.ucdavis.edu/cow.jpg",
                              "https://files.asmith.ucdavis.edu/cow.jpg",
                              "https://files.asmith.ucdavis.edu/Oat.jpg",
                              "https://files.asmith.ucdavis.edu/soy.jpg"))



#Load milk tasting data
GET("https://files.asmith.ucdavis.edu/Milkrubric.xlsx", write_disk(ttf1 <- tempfile(fileext = ".xlsx")))
milk <- read_excel(ttf1,sheet="Data") 


# Make latte ratings plot
bar_plot_latte <- milk %>%
  mutate(Criteria=ifelse(Criteria=="Taste","Non-Blind",Criteria)) %>%
  select(-Nonfat,-Whole) %>%
  filter(Item=="Latte") %>%
  filter(Criteria %in% c("Blind","Non-Blind")) %>%
  filter(!Rater=="myfooddiary") %>%
  pivot_longer(cols =c("Almond","Coconut","Dairy (2%)","Oat","Soy")) %>%
  group_by(Rater,name,Criteria) %>%
  summarise(rating=sum(value)) %>%
  mutate(name = factor(name, levels=c("Soy","Oat","Dairy (2%)","Coconut","Almond"))) %>%
  left_join(logos) %>%
  ggplot(aes(x=name,y=rating,label=rating,fill=Criteria)) +
    geom_bar(stat="identity",  show.legend = FALSE) + 
    #geom_text(size = 6, position = position_stack(vjust = 1.04))+
    geom_image(aes(x=name,image=logo,y=.9),size = 0.15, hjust = 1,inherit.aes = FALSE) +
    facet_grid(cols=vars(Rater),rows=vars(Criteria))+
    ggtitle(paste0("Latte Ratings by Milk Type"))+
    theme_minimal()+
    coord_flip()+
    scale_fill_brewer(palette="Dark2")+
    scale_y_continuous(limits=c(0,5),oob = rescale_none,labels=NULL)  +
    labs(x="",y="")+
    theme(plot.title = element_text(hjust = 0.5,size=20), text = element_text(size=18))

# Draw and save plot
bar_plot_latte
ggsave("latte_ratings.png")





# Make milk ratings plot
bar_plot_milk <- milk %>%
  mutate(Criteria=ifelse(Criteria=="Taste","Non-Blind",Criteria)) %>%
  select(-Nonfat,-Whole) %>%
  filter(Item=="Milk") %>%
  filter(Criteria %in% c("Blind","Non-Blind")) %>%
  filter(!Rater=="myfooddiary") %>%
  pivot_longer(cols =c("Almond","Coconut","Dairy (2%)","Oat","Soy")) %>%
  group_by(Rater,name,Criteria) %>%
  summarise(rating=sum(value)) %>%
  mutate(name = factor(name, levels=c("Soy","Oat","Dairy (2%)","Coconut","Almond"))) %>%
  left_join(logos) %>%
  ggplot(aes(x=name,y=rating,label=rating,fill=Criteria)) +
    geom_bar(stat="identity",show.legend=FALSE) + 
    #geom_text(size = 6, position = position_stack(vjust = 1.04))+
    geom_image(aes(x=name,image=logo,y=.9),size = 0.15, hjust = 1,inherit.aes = FALSE) +
    facet_grid(cols=vars(Rater),rows=vars(Criteria))+
    ggtitle(paste0("Milk Ratings"))+
    theme_minimal()+
    coord_flip()+
    scale_fill_brewer(palette="Dark2")+
    scale_y_continuous(limits=c(0,5),oob = rescale_none,labels=NULL)  +
    labs(x="",y="")+
    theme(plot.title = element_text(hjust = 0.5,size=20), text = element_text(size=18))

# Draw and save plot
bar_plot_milk
ggsave("milk_ratings.png")






# Make calories plot
bar_plot_cal <- milk %>%
  rename(`Dairy (NF)`=Nonfat,`Dairy (Whole)`=Whole) %>%
  filter(Criteria %in% c("Carbohydrate","Fat","Protein")) %>%
  pivot_longer(cols =c("Almond","Coconut","Dairy (NF)","Dairy (2%)","Dairy (Whole)","Oat","Soy")) %>%
  mutate(value=ifelse(Criteria=="Protein",4*value,value)) %>%
  mutate(value=ifelse(Criteria=="Carbohydrate",4*value,value)) %>%
  mutate(value=ifelse(Criteria=="Fat",9*value,value)) %>%
  mutate(name = factor(name, levels=c("Almond","Coconut","Dairy (NF)","Dairy (2%)","Dairy (Whole)","Oat","Soy"))) %>%
  mutate(Criteria = factor(Criteria, levels=c( "Fat","Protein", "Carbohydrate"))) %>%
  left_join(logos) %>%
  ggplot(aes(x=name,y=value,fill=Criteria,label=value)) +
    geom_bar(stat="identity", position = "stack") + 
    geom_text(size = 4, position = position_stack(vjust = .9))+
    geom_image(aes(x=name,image=logo,y=-10),size = 0.12, hjust = .5,inherit.aes = FALSE) +
    ggtitle(paste0("Calories in a 12oz Iced Latte"))+
    scale_fill_brewer(palette="Set2")+
    theme_minimal()+
    labs(x="Milk Type",y="Calories",fill="Nutrient",caption="source: myfooddiary.com")+
    scale_x_discrete()+
    scale_y_continuous(breaks = seq(0,120,20),limits=c(-10,120))+
    theme(plot.title = element_text(hjust = 0.5,size=20), text = element_text(size=16))

# Draw and save plot
bar_plot_cal
ggsave("latte_calories.png")




# Make sugar scatter plot
scatter_plot <- milk %>%
  select(-Nonfat,-Whole) %>%
  filter(Criteria %in% c("Sugar")) %>%
  pivot_longer(cols =c("Almond","Coconut","Dairy (2%)","Oat","Soy")) %>%
  select(-Item,-Rater,-Criteria) %>%
  rename(Sugar=value) %>%
  right_join(milk %>% select(-Nonfat,-Whole) %>% 
               filter(Criteria %in% c("Blind")) %>%
               pivot_longer(cols =c("Almond","Coconut","Dairy (2%)","Oat","Soy"))) %>%
  left_join(logos) %>%
  ggplot(aes(x=Sugar,y=value)) +
    geom_image(aes(image=logo),asp = 1.1,size=.15,position=position_jitter(width=0, height=0.3)) + 
    ggtitle(paste0("Blind Milk Ratings vs Sugar Content"))+
    facet_grid(rows=vars(Item),cols=vars(Rater)) +
    xlim(0,10)+
    ylim(0,5)+
    labs(x="Sugar (g)",y="Rating (0-5)")+
    theme(plot.title = element_text(hjust = 0.5,size=20), text = element_text(size=18),panel.background = element_rect(fill="white",color="gray"))

scatter_plot
ggsave("milk_scatter.png")



# Make heat map
bar_plot_heat <- milk %>%
  filter(Criteria %in% c("Almond","Coconut","Dairy (2%)","Oat","Soy")) %>%
  select(-Nonfat,-Whole) %>%
  pivot_longer(cols =c("Almond","Coconut","Dairy (2%)","Oat","Soy")) %>%
  left_join(logos) %>%
  ggplot(aes(x=name,y=factor(Criteria,levels=c("Soy","Oat","Dairy (2%)","Coconut","Almond")),fill=factor(value))) +
    geom_tile(color="gray",show.legend=FALSE) + 
    ggtitle(paste0("Did We Know What We Were Drinking?"))+
    facet_grid(cols=vars(Rater),rows=vars(Item))+
    scale_fill_manual(values=c("red","green"),na.value="white") +
    labs(x="Guess",y="Truth",fill="Nutrient")+
    theme(plot.title = element_text(hjust = 0.5,size=20), text = element_text(size=16),axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))

# Draw and save plot
bar_plot_heat
ggsave("milk_heat.png")






#Load milk consumption data
GET("https://www.ers.usda.gov/webdocs/DataFiles/48685/pcconsp_1_.xlsx", write_disk(tf2 <- tempfile(fileext = ".xlsx")))
cons <- read_excel(tf2,sheet="per capita",range = cell_rows(2:50)) %>%
  rename(`Fluid Milk`=`Fluid \r\nbeverage \r\nmilk 2`) %>%
  mutate(`American Cheese`=as.numeric(Cheese)) %>%
  mutate(`Other Cheese`=as.numeric(`...4`)) %>%
  mutate(`Ice Cream`=as.numeric(`Frozen products`)+as.numeric(`...12`)) %>%
  mutate(Yogurt=as.numeric(`Yogurt, other than frozen`)) %>%
  rename(Total = `All products, milk-fat \r\nmilk-equivalent basis`) %>%
  select(Year,`Fluid Milk`,`American Cheese`,`Other Cheese`,`Ice Cream`,`Yogurt`,Butter,Total) %>%
  filter(!is.na(Year))


# Line graph of fluid milk consumption
plot_milk_cons <- cons  %>%
  select(Year,`Fluid Milk`) %>%
  ggplot(aes(x=Year,y=`Fluid Milk`))+
    geom_line(size=1.5,color="blue") +
    ggtitle(paste0("Annual US Fluid Milk Consumption Per Capita"))+
    theme_minimal()+
    labs(x="Year",y="lbs/person",caption="Data source: USDA ERS")+
    theme(plot.title = element_text(hjust = 0.5,size=16), text = element_text(size=14))

# Draw and save plot
plot_milk_cons
ggsave("milk_cons.png")



# Line graph of total milk consumption
plot_total_milk_cons <- cons  %>%
  select(Year,Total) %>%
  ggplot(aes(x=Year,y=Total))+
    geom_line(size=1.5,color="blue") +
    ggtitle(paste0("Annual US Total Milk Consumption Per Capita"))+
    theme_minimal()+
    labs(x="Year",y="lbs/person",caption="Data source: USDA ERS")+
    theme(plot.title = element_text(hjust = 0.5,size=16), text = element_text(size=14))

# Draw and save plot
plot_total_milk_cons
ggsave("total_milk_cons.png")





# Line graph of dairy consumption
plot_dairy_product_cons <- cons  %>%
  select(-`Fluid Milk`,-Total) %>%
  pivot_longer(-Year) %>%
  ggplot(aes(x=Year,y=value,color=factor(name,levels=c("Butter","American Cheese","Other Cheese","Ice Cream","Yogurt"))))+
    geom_line(size=1.5) +
    ggtitle(paste0("Annual US Dairy Product Consumption Per Capita"))+
    theme_minimal()+
    labs(x="Year",y="lbs/person",color="Product",caption="Data source: USDA ERS")+
    theme(plot.title = element_text(hjust = 0.5,size=16), text = element_text(size=14))

# Draw and save plot
plot_dairy_product_cons
ggsave("dairy_product_cons.png")