---
title: "R Notebook"
output: html_notebook
---

```{r}
library(readxl)
library(lubridate)
library(writexl)
library(ggplot2)
library(colorspace)
library(reshape2)
library(psych)
library(car)

setwd("/Users/Libby/Desktop/UVA/Research/Data Inventory/Additional Measurements/Sediment Accumulation")

all = read_excel("Sep_July.xlsx")
describeBy(all$`Mass(gpercm2)`, group = all$`Section/Location`)
```
```{r}
describeBy(all$`Mass(gpercm2)`, group = all$Location)
```


```{r}
describeBy(all$`Mass(9percm2)`, group = all$`Section`)
```



```{r}
library(lm)
sed.lmm<- lmer(all$`Mass(gpercm2)` ~ Section + Date + (1 | Section:Location), data = all)
summary(sed.lmm)
```
```{r}
all$sedlog = log10(all$`Mass(gpercm2)`)
sed.anova<- aov(all$sedlog ~ Section, data = all)
summary(sed.anova) 
```

```{r}
sed.anova2<- aov(all$sedlog ~ Location, data = all)
summary(sed.anova2)
TukeyHSD(sed.anova2)
```
```{r}
sed.anova_resids <- sed.anova$residuals
qqPlot(sed.anova_resids) 
hist(sed.anova_resids)
#Shapiro Wilk's Test for Normality
shapiro.test(sed.anova_resids) #p > 0.05, so points are normally distributed 
leveneTest(all$`Mass(gpercm2)`, all$Section) #p value > 0.05, so there is homogeneity of variation
```
```{r}
sed.anova2_resids <- sed.anova2$residuals
qqPlot(sed.anova2_resids) 
hist(sed.anova2_resids)
#Shapiro Wilk's Test for Normality
shapiro.test(sed.anova2_resids) #p > 0.05, so points are normally distributed 
leveneTest(all$`Mass(gpercm2)`, all$Location) #p value > 0.05, so there is homogeneity of variation
```
```{r}
sed.anova3<- aov(all$sedlog ~ Date, data = all)
summary(sed.anova3)
TukeyHSD(sed.anova3)
```
```{r}
mypal=c("red1", "blue3")
ggplot(all,aes(x=Section,y=all$`Mass(gpercm2)`, fill=Section)) + geom_boxplot(show.legend = FALSE) +
  scale_fill_manual(values=mypal) +
  ylab("Sediment mass (g per m2)") + 
  xlab("Section") + 
  theme_light() +
  theme(aspect.ratio=3/2)
ggsave("boxplot.png")
```

```{r}
#calculate 95% confidence intervals
S4ci = 0.0107 + c(-2.11,2.11)*0.0023
S7ci = 0.00603 + c(-2.13,2.13)*0.0019
```


```{r}
#First create vectors for your grouping variable, dependent variable, and 95% CIs
#independent variable
group = c("S4", "S7")
#depependent variable 
yvar = c(0.0107, 0.00603)
#95% ci of means for each group
ci = c(0.0049,0.0041)

#Combine into data frame
errorbars = data.frame(group, yvar, ci)

#create barplot w/error bars
```


```{r}
mypal=c("red1","blue3")
ggplot(errorbars, aes(x=group, y=yvar,fill=group)) +
  geom_bar(stat="identity", width = 0.3, color="black", lwd=0.3, show.legend = FALSE)+
  geom_errorbar(aes(ymin=yvar-ci, ymax=yvar+ci), width=0.1, color="black") +
  theme_light() +
  scale_fill_manual(values=mypal) +
  ylab("Mean sediment mass (g per m2)") +
  xlab("Section") + 
  theme(aspect.ratio=3/2) 
ggsave("barchart3.png")
```



