hist(MPdata$MPConc)
MPdata$log.MPConc <- log(MPdata$MPConc)
hist(MPdata$log.MPConc)

#log data looks pretty normal

shapiro.test(MPdata$log.MPConc)

#log MP conc data are statistically normal

#some summary stats for tables

aggregate (MPConc~Sex+Age, MPdata, mean)
aggregate (MPConc~Sex+Age, MPdata, sd)
aggregate (MPConc~Sex+Age, MPdata, length)


aggregate (MPConc~Sex, MPdata, mean)
aggregate (MPConc~Sex, MPdata, sd)
aggregate (MPConc~Sex, MPdata, length)

aggregate (MPConc~Age, MPdata, mean)
aggregate (MPConc~Age, MPdata, sd)
aggregate (MPConc~Age, MPdata, length)




########## Stats ############

#Age and sex differences?

male.dat <-MPdata[which(MPdata$Sex=="M"),]
female.dat <-MPdata[which(MPdata$Sex=="F"),]

t.test(male.dat$log.MPConc,female.dat$log.MPConc)


adult.dat <-MPdata[which(MPdata$Sex=="A"),]
juv.dat <-MPdata[which(MPdata$Sex=="J"),]

t.test(adult.dat$log.MPConc,juv.dat$log.MPConc)


#Effect of MP conc on body mass with a graph
hist(MPdata$BatMass)

plot(MPdata$MPConc~MPdata$BatMass)

plot(MPdata$BatMass~MPdata$log.MPConc)

cor.test(MPdata$BatMass,MPdata$log.MPConc)

mass.mod <- glm(MPdata$BatMass~MPdata$log.MPConc)
summary(mass.mod)
anova(mass.mod)

hist(mass.mod$residuals)
shapiro.test(mass.mod$residuals) #residuals are normal

predicted_df <- df(mpg_pred=predict(mass.mod, df))


library(ggplot2)


custom.cb.palette <- c("#7b3294", "#fdae61", "#a6dba0", "#008837")

df <- MPdata


ggplot(df, aes(log.MPConc, BatMass)) + geom_point(aes(color=Sex), size = 3) + scale_color_manual(values=c("#7b3294", "#fdae61", "#a6dba0", "#008837")) + ylim (5, 25) +xlab("Microplastic Concentration (log(n/g))") + ylab("Bat Mass (g)") +theme_minimal() + geom_smooth(method="lm", colour="black", size=1.0)


#saving as high res

ggsave("Plot_BatMass_MP_600dpi.jpg", units="in", width=5, height=4, dpi=600)



###### Are Bat concentrations sig diff than controls?

MP.EPFU.BLANK.dat <-MPdata #reload all data and assign 

blank.mod <-aov(MP.EPFU.BLANK.dat$log.MPConc ~ MP.EPFU.BLANK.dat$Species)
anova(blank.mod)
summary(lm(blank.mod))
TukeyHSD(blank.mod)

age.sex.mod <- aov(df$log.MPConc~ df$Sex + df$Age)
anova(age.sex.mod)
summary(age.sex.mod)




########## Make a figure showing body mass per month with mp conc overlayed
library(ggplot2)
library(tidyr)
library(dplyr)

par(mfrow = c(2, 1)) 

mids <-barplot(MPBodyMassdata$MeanBatMass, names.arg = MPBodyMassdata$CollectionMonth,
            
               ylab= "Bat Body Mass (g)",
               ylim=c(0,30),
               beside=TRUE  
)


arrows (mids, MPBodyMassdata$MeanBatMass, mids, MPBodyMassdata$MeanBatMass  + MPBodyMassdata$SDBatMass, code=3, angle =90, length=0.1)
arrows (mids, MPBodyMassdata$MeanBatMass, mids, MPBodyMassdata$MeanBatMass  - MPBodyMassdata$SDBatMass, code=3, angle =90, length=0.1)

mids2 <-barplot(MPBodyMassdata$MeanMPConc, names.arg = MPBodyMassdata$CollectionMonth,
                     xlab="Collection Month", 
                     ylab= "MP Concentration (n/g)",
                     ylim=c(0,30),
                     beside=TRUE  
)


arrows (mids2, MPBodyMassdata$MeanMPConc, mids2, MPBodyMassdata$MeanMPConc  + MPBodyMassdata$SDMPConcentration, code=3, angle =90, length=0.1)
arrows (mids2, MPBodyMassdata$MeanMPConc, mids2, MPBodyMassdata$MeanMPConc  - MPBodyMassdata$SDMPConcentration, code=3, angle =90, length=0.1)

ggsave("Plot_month.tiff", units="in", width=5, height=8, dpi=600)