library(tidyverse)
met <- read_csv("https://denvirlab.marshall.edu/BMR617-2022/data/TH-B6-metabolic.csv") %>%
  separate(MouseID, into=c("Strain", "Diet", "ID"), sep='-') %>%
  mutate(Strain = factor(Strain), Diet=factor(Diet))

aov.full <- aov(Cholesterol ~ Strain * Diet, data = met)
aov.full$coefficients
confint(aov.full)
summary(aov.full)

png(filename = "figure1.png", width=4*960, height = 4*480, res=4*72)
ggplot(met, aes(x=Strain, y=Cholesterol, fill=Diet)) +
  geom_boxplot(outlier.shape = NA) + 
  geom_point(position=position_jitterdodge(jitter.width = 0.1, dodge.width = 0.75)) +
  ylab("Cholesterol Level (mg/dl)") +
  ggtitle("Cholesterol Level by Mouse Strain and Diet")
dev.off()

sessionInfo()

citation()
citation("tidyverse")

resultsTable <- as_tibble(confint(aov.full), rownames="Parameter")
resultsTable <- cbind(resultsTable, tibble(Estimate = aov.full$coefficients))
resultsTable <- resultsTable %>%
  mutate(`2.5 %` = round(`2.5 %`, digits = 2),
         `97.5 %` = round(`97.5 %`, digits = 2),
         Estimate = round(Estimate, digits = 2)) %>%
  mutate(`95% Confidence Interval`=paste0("[", `2.5 %`, ", ", `97.5 %`, "]"))
resultsTable <- resultsTable %>%
  select(Parameter, Estimate, `95% Confidence Interval`)

png("lowResChol.png", width=480, height=480)
ggplot(met, aes(x=Strain, y=Cholesterol, fill=Diet)) +
  geom_boxplot(outlier.shape = NA) + 
  geom_point(position=position_jitterdodge(jitter.width = 0.1, dodge.width = 0.75)) +
  ylab("Cholesterol Level (mg/dl)") +
  ggtitle("Cholesterol Level by Mouse Strain and Diet")
dev.off()

png("hiResChol.png", width=4*480, height=4*480)
ggplot(met, aes(x=Strain, y=Cholesterol, fill=Diet)) +
  geom_boxplot(outlier.shape = NA) + 
  geom_point(position=position_jitterdodge(jitter.width = 0.1, dodge.width = 0.75)) +
  ylab("Cholesterol Level (mg/dl)") +
  ggtitle("Cholesterol Level by Mouse Strain and Diet")
dev.off()


png("hiResChol.png", width=4*480, height=4*480, res=4*72)
ggplot(met, aes(x=Strain, y=Cholesterol, fill=Diet)) +
  geom_boxplot(outlier.shape = NA) + 
  geom_point(position=position_jitterdodge(jitter.width = 0.1, dodge.width = 0.75)) +
  ylab("Cholesterol Level (mg/dl)") +
  ggtitle("Cholesterol Level by Mouse Strain and Diet")
dev.off()