---
title: "Sed OM"
output: html_notebook
---

```{r}
library(readxl)
library(lubridate)
library(writexl)
library(ggplot2)
library(colorspace)
library(reshape2)
library(lawstat)
library(car)
setwd("/Users/Libby/Desktop/UVA/Research/Data Inventory/Additional Measurements/Infauna and OM")
SedOM = read_excel("Infauna_OM.xlsx", sheet="Sed OM")
```

```{r}
SedOM$Sample = as.factor(SedOM$Sample)
SedOM$Date = as.factor(SedOM$Date)
SedOM$Site = as.factor(SedOM$Site)
SedOM$`N/O` = as.factor(SedOM$`N/O`)
SedOM$OMPLog = log10(SedOM$`OM percent`)
SedOM$OMPSqrt = sqrt(SedOM$`OM percent`)

model1 = aov(SedOM$`OM percent`~Date, data=SedOM)
summary(model1)
TukeyHSD(model1)

model2 = aov(SedOM$OMPLog~Date, data=SedOM)
summary(model2)
TukeyHSD(model2)

model3 = aov(SedOM$OMPSqrt~Date, data=SedOM)
summary(model3)
TukeyHSD(model3)
```

```{r}
model1_resids <- model1$residuals
qqPlot(model1_resids) 
hist(model1_resids)
#Shapiro Wilk's Test for Normality
shapiro.test(model1_resids) #p > 0.05, so points are normally distributed 
leveneTest(SedOM$`OM percent`, SedOM$Date) #p value > 0.05, so there is homogeneity of variation
```
```{r}
model2_resids <- model2$residuals
qqPlot(model2_resids) 
hist(model2_resids)
#Shapiro Wilk's Test for Normality
shapiro.test(model2_resids) #p > 0.05, so points are normally distributed 
leveneTest(SedOM$`OM percent`, SedOM$Date) #p value > 0.05, so there is homogeneity of variation
```
```{r}
model3_resids <- model3$residuals
qqPlot(model3_resids) 
hist(model3_resids)
#Shapiro Wilk's Test for Normality
shapiro.test(model3_resids) #p > 0.05, so points are normally distributed 
leveneTest(SedOM$`OM percent`, SedOM$Date) #p value > 0.05, so there is homogeneity of variation
```


```{r}
SedOM$Date = factor(SedOM$Date,levels = c("Sep 2020", "Jul 2021", "Oct 2021", "Jul 2022"))
ggplot(SedOM,aes(x=Date,y=SedOM$`OM percent`)) + geom_boxplot(fill="sienna",show.legend = FALSE) +
  ylab("% sediment organic matter content") + 
  xlab("Date") + 
  theme_light() +
  theme(aspect.ratio=2/3)
ggsave("sedOMboxplot.png")
```
```{r}
library(psych)
describeBy(SedOM$`OM percent`, group = SedOM$Date)
```





