library(dplyr)
library(lme4)
library(ggplot2)
library(multcomp)
library(tidyverse)

setwd("/Users/harliecope/Desktop/GenX:PFOA Study/GenX PFOA Mouse Study/GenXPFOA Study")
## Set working directory

## Read in data file containing animal IDs, treatment groups, etc
x<- read_csv("~/Desktop/GenX:PFOA Study/GenX PFOA Mouse Study/QC Data Files/csv for R/R format PND 5.5 pooled serum and urine.csv", na="Na")
x$group <-as.factor(x$group)
str(x)

names(x)
levels(x$group)


summary.stats<-x%>%
  dplyr::group_by(group, matrix) %>%
  dplyr::summarise(mean.matrix = mean(ug.ml, na.rm = TRUE),
                   sd.matrix = sd(ug.ml, na.rm = TRUE))
summary.stats

write.csv(summary.stats, "PND 5.5 pooled serum and urine summary stats.csv",row.names = F)

summary.stats$ug.ml<-summary.stats$mean.matrix
summary.stats$ug.ml.lo<-summary.stats$mean.matrix-summary.stats$sd.matrix
summary.stats$ug.ml.hi<-summary.stats$mean.matrix+summary.stats$sd.matrix

##statistical analysis

serum<-subset(x, matrix=="pooled serum")
serumPFOA<-subset(serum, chem=="PFOA")
serumGenX<-subset(serum, chem=="GenX")

urine<-subset(x, matrix=="pooled urine")
urinePFOA<-subset(urine, chem=="PFOA")
urineGenX<-subset(urine, chem=="GenX")

#serum PFOA
summary(glht(aov(ng.ml~group,data=serumPFOA), linfct=mcp(group="Dunnett")))


#serum GenX
summary(glht(aov(ng.ml~group,data=serumGenX), linfct=mcp(group="Dunnett")))



p1<-ggplot(summary.stats, aes(x=group, y=ug.ml))+
  theme_bw()+
  theme(legend.position="none",
        axis.text=element_text(size=12),
        axis.text.x=element_text(angle =-45, hjust=0),
        axis.title=element_text(size=14),
        plot.margin=unit(c(2,2,2,2),"cm"))+
  geom_jitter(data=x, aes(x=group, y=ug.ml), position = position_jitter(width = 0.2))+
  geom_crossbar(data=summary.stats, aes(ymin = ug.ml.lo, ymax = ug.ml.hi), width = 0.4)+
  facet_wrap(matrix~., scales="free")+
  ggtitle("PND 5.5 Pooled Pup Urine and Serum Dosimetry")+
  ylab("Analyte concentration (ug/mL)\n")+
  xlab("\nGroup")
p1

tiff("PND 5.5 pooled serum urine dosimetry.tiff", units="in", width=9, height=7, res=300)
# insert ggplot code
p1
dev.off()