library(tidyverse)

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

ggplot(met, aes(x=paste(Strain, Diet, sep=':'), y=Cholesterol, color=Diet)) + 
  geom_boxplot() +
  xlab("Group") + ylab("Cholesterol (mg/dl)")

chol.strain.diet <- aov(Cholesterol ~ Strain + Diet, data=met)
summary(chol.strain.diet)

chol.strain.diet$coefficients

met %>% 
  group_by(Strain, Diet) %>%
  summarize(Mean_Cholesterol = mean(Cholesterol))   

confint(chol.strain.diet)


gluc.strain.diet <- aov(Glucose ~ Strain + Diet, data=met)
summary(gluc.strain.diet)

gluc.strain.diet$coefficients

met %>% 
  group_by(Strain, Diet) %>%
  summarize(Mean_Glucose = mean(Glucose))   

gluc.strain.diet.int <- aov(Glucose ~ Strain + Diet + Strain:Diet, data=met)
summary(gluc.strain.diet.int)
confint(chol.strain.diet)

gluc.strain.diet.int$coefficients

confint(gluc.strain.diet.int)

params <- as_tibble(confint(gluc.strain.diet.int), rownames="Parameter")
params <- mutate(params, estimate = gluc.strain.diet.int$coefficients)
ggplot(params, aes(x=Parameter, y=estimate)) +
  geom_point() +
  geom_errorbar(aes(ymin=`2.5 %`, ymax=`97.5 %`), width=0.3) +
  geom_hline(yintercept=0, color='#00b140') +
  ylab("Estimated effect (mg/dl)") +
  ggtitle("95% Confidence Intervals for effect on Glucose level")