---
title: "Marsh Grass"
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/Marsh Grass")
MarshGrassSep2020 = read_excel("MarshGrass.xlsx", sheet="Sep 2020")
MarshGrassSep2020$Sample = as.character(MarshGrassSep2020$Sample)
MarshGrassJul2021 = read_excel("MarshGrass.xlsx", sheet="July 2021")
MarshGrassJul2021$Sample = as.character(MarshGrassJul2021$Sample)
MarshGrassOct2021 = read_excel("MarshGrass.xlsx", sheet="Oct 2021")
MarshGrassOct2021$Sample = as.character(MarshGrassOct2021$Sample)
all = read_excel("MarshGrass.xlsx", sheet="all")
```

```{r}
#testing for homogeneity of variation 
levene.test(MarshGrassSep2020$Height, MarshGrassSep2020$Site) #p value > 0.05, so there is homogeneity of variation
levene.test(MarshGrassSep2020$Density, MarshGrassSep2020$Site) #p value > 0.05, so there is homogeneity of variation
levene.test(MarshGrassSep2020$Height, MarshGrassSep2020$Location) #p value < 0.05, so there is NO homogeneity of variation
levene.test(MarshGrassSep2020$Density, MarshGrassSep2020$Location) #p value > 0.05, so there is homogeneity of variation

levene.test(MarshGrassJul2021$Height, MarshGrassJul2021$Site) #p value > 0.05, so there is homogeneity of variation
levene.test(MarshGrassJul2021$Density, MarshGrassJul2021$Site) #p value > 0.05, so there is homogeneity of variation
levene.test(MarshGrassJul2021$Height, MarshGrassJul2021$Location) #p value > 0.05, so there is homogeneity of variation
levene.test(MarshGrassJul2021$Density, MarshGrassJul2021$Location) #p value > 0.05, so there is homogeneity of variation

levene.test(MarshGrassOct2021$Height, MarshGrassOct2021$Site) #p value > 0.05, so there is homogeneity of variation
levene.test(MarshGrassOct2021$Density, MarshGrassOct2021$Site) #p value < 0.05, so there is NO homogeneity of variation
levene.test(MarshGrassOct2021$Height, MarshGrassOct2021$Location) #p value > 0.05, so there is homogeneity of variation
levene.test(MarshGrassOct2021$Density, MarshGrassOct2021$Location) #p value > 0.05, so there is homogeneity of variation
```


```{r}
mypal=c("red1", "blue3")
ggplot(data = MarshGrassSep2020, aes(x=Site, y=Height, fill=Site)) +
  geom_boxplot() +
  ylab("Height (cm)") +
  xlab("Section") +
  theme_light() + 
  theme(aspect.ratio=1/1.5) +
  scale_fill_manual(values=mypal) +
  ggtitle("September 2020")
ggsave("Sep20Height.png")
#sig difference, p<0.05
```
```{r}
mypal=c("red1", "blue3")
ggplot(data = MarshGrassJul2021, aes(x=Site, y=Height, fill=Site)) +
  geom_boxplot() +
  ylab("Stem height (cm)") +
  xlab("Section") +
  theme_light() + 
  theme(aspect.ratio=1/1.5) +
  scale_fill_manual(values=mypal) +
  ggtitle("July 2021")
ggsave("Jul21HEIGHT.png")
#sig difference, p<0.05
```
```{r}
mypal=c("red1", "blue3")
ggplot(data = MarshGrassOct2021, aes(x=Site, y=Height, fill=Site)) +
  geom_boxplot() +
  ylab("Height (cm)") +
  xlab("Section") +
  theme_light() + 
  theme(aspect.ratio=1/1.5) +
  scale_fill_manual(values=mypal) +
  ggtitle("October 2021")
ggsave("Oct21Height.png")
#no sig difference, p>0.05
```
```{r}
mypal=c("red1", "blue3")
ggplot(data = MarshGrassSep2020, aes(x=Site, y=Density, fill=Site)) +
  geom_boxplot() +
  ylab("Density") +
  xlab("Section") +
  theme_light() + 
  theme(aspect.ratio=1/1.5) +
  scale_fill_manual(values=mypal) +
  ggtitle("September 2020")
ggsave("Sep20Density.png")
#sig difference, p<0.05
```
```{r}
mypal=c("red1", "blue3")
ggplot(data = MarshGrassJul2021, aes(x=Site, y=MarshGrassJul2021$`Density(percm2)`, fill=Site)) +
  geom_boxplot() +
  ylab("Stem density (per cm2)") +
  xlab("Section") +
  theme_light() + 
  theme(aspect.ratio=1/1.5) +
  scale_fill_manual(values=mypal) +
  ggtitle("July 2021")
ggsave("Jul21DENSITY.png")
#sig difference, p<0.05
```
```{r}
mypal=c("red1", "blue3")
ggplot(data = MarshGrassOct2021, aes(x=Site, y=Density, fill=Site)) +
  geom_boxplot() +
  ylab("Density") +
  xlab("Section") +
  theme_light() + 
  theme(aspect.ratio=1/1.5) +
  scale_fill_manual(values=mypal) +
  ggtitle("October 2021")
ggsave("Oct21Density.png")
#no sig difference, p>0.05
```
```{r}
mypal=c("red1", "blue3")
ggplot(data = MarshGrassSep2020, aes(x=Location, y=Height, fill=Site)) +
  geom_boxplot() +
  ylab("Height (cm)") +
  theme_light() + 
  theme(aspect.ratio=1/1.5) +
  scale_fill_manual(values=mypal) +
  ggtitle("September 2020")
ggsave("Sep20Height2.png")
```
```{r}
mypal=c("red1", "blue3")
ggplot(data = MarshGrassJul2021, aes(x=Location, y=Height, fill=Site)) +
  geom_boxplot() +
  ylab("Height (cm)") +
  theme_light() + 
  theme(aspect.ratio=1/1.5) +
  scale_fill_manual(values=mypal) +
  ggtitle("July 2021")
ggsave("Jul21Height2.png")
```
```{r}
mypal=c("red1", "blue3")
ggplot(data = MarshGrassOct2021, aes(x=Location, y=Height, fill=Site)) +
  geom_boxplot() +
  ylab("Height (cm)") +
  theme_light() + 
  theme(aspect.ratio=1/1.5) +
  scale_fill_manual(values=mypal) +
  ggtitle("October 2021")
ggsave("Oct21Height2.png")
```
```{r}
mypal=c("red1", "blue3")
ggplot(data = MarshGrassSep2020, aes(x=Location, y=Density, fill=Site)) +
  geom_boxplot() +
  ylab("Density") +
  theme_light() + 
  theme(aspect.ratio=1/1.5) +
  scale_fill_manual(values=mypal) +
  ggtitle("September 2020")
ggsave("Sep20Density2.png")
```
```{r}
mypal=c("red1", "blue3")
ggplot(data = MarshGrassJul2021, aes(x=Location, y=Density, fill=Site)) +
  geom_boxplot() +
  ylab("Density") +
  theme_light() + 
  theme(aspect.ratio=1/1.5) +
  scale_fill_manual(values=mypal) +
  ggtitle("July 2021")
ggsave("July21Density2.png")
```
```{r}
mypal=c("red1", "blue3")
ggplot(data = MarshGrassOct2021, aes(x=Location, y=Density, fill=Site)) +
  geom_boxplot() +
  ylab("Density") +
  theme_light() + 
  theme(aspect.ratio=1/1.5) +
  scale_fill_manual(values=mypal) +
  ggtitle("October 2021")
ggsave("Oct21Density2.png")
```
```{r}
mypal=c("red1", "blue3")
ggplot(data = MarshGrassSep2020, aes(x=Location2, y=Height, fill = Location2)) +
  geom_boxplot() +
  ylab("Height (cm)") +
  xlab("Location") +
  theme_light() + 
  theme(aspect.ratio=1/1.5) +
  scale_fill_manual(values=mypal) +
  ggtitle("September 2020")
ggsave("Sep20Height3.png")
#no sig difference, p > 0.05
```
```{r}
mypal=c("red1", "blue3")
ggplot(data = MarshGrassJul2021, aes(x=Location2, y=Height, fill = Location2)) +
  geom_boxplot() +
  ylab("Stem height (cm)") +
  xlab("Location") +
  theme_light() + 
  theme(aspect.ratio=1/1.5) +
  scale_fill_manual(values=mypal) +
  ggtitle("July 2021")
ggsave("Jul21Height3.png")
#no sig difference, p > 0.05
```
```{r}
mypal=c("red1", "blue3")
ggplot(data = MarshGrassOct2021, aes(x=Location2, y=Height, fill = Location2)) +
  geom_boxplot() +
  ylab("Height (cm)") +
  xlab("Location") +
  theme_light() + 
  theme(aspect.ratio=1/1.5) +
  scale_fill_manual(values=mypal) +
  ggtitle("October 2021")
ggsave("Oct21Height3.png")
#no sig difference, p > 0.05
```
```{r}
mypal=c("red1", "blue3")
ggplot(data = MarshGrassSep2020, aes(x=Location2, y=Density, fill = Location2)) +
  geom_boxplot() +
  ylab("Density") +
  xlab("Location") +
  theme_light() + 
  theme(aspect.ratio=1/1.5) +
  scale_fill_manual(values=mypal) +
  ggtitle("September 2020")
ggsave("Sep20Density3.png")
#no sig difference, p > 0.05
```
```{r}
mypal=c("red1", "blue3")
ggplot(data = MarshGrassJul2021, aes(x=Location2, y=Density, fill = Location2)) +
  geom_boxplot() +
  ylab("Stem density (per cm2") +
  xlab("Location") +
  theme_light() + 
  theme(aspect.ratio=1/1.5) +
  scale_fill_manual(values=mypal) +
  ggtitle("July 2021")
ggsave("Jul21Density3.png")
#no sig difference, p > 0.05
```


```{r}
mypal=c("red1", "blue3")
ggplot(data = MarshGrassOct2021, aes(x=Location2, y=Density, fill = Location2)) +
  geom_boxplot() +
  ylab("Density") +
  xlab("Location") +
  theme_light() + 
  theme(aspect.ratio=1/1.5) +
  scale_fill_manual(values=mypal) +
  ggtitle("October 2021")
ggsave("Oct21Density3.png")
#no sig difference, p > 0.05
```


```{r}
t.test(Density~Site,data=MarshGrassSep2020,alternative="two.sided",var.equal=TRUE)
```
```{r}
t.test(Density~Site,data=MarshGrassJul2021,alternative="two.sided",var.equal=TRUE)
```
```{r}
t.test(Density~Site,data=MarshGrassOct2021,alternative="two.sided",var.equal=TRUE)
```
```{r}
t.test(Height~Site,data=MarshGrassSep2020,alternative="two.sided",var.equal=TRUE)
```
```{r}
t.test(Height~Site,data=MarshGrassJul2021,alternative="two.sided",var.equal=TRUE)
```
```{r}
t.test(Height~Site,data=MarshGrassOct2021,alternative="two.sided",var.equal=TRUE)
```
```{r}
model1 = aov(Density~Location, data=MarshGrassSep2020)
summary(model1)
```
```{r}
model2 = aov(Density~Location, data=MarshGrassJul2021)
summary(model2)
```
```{r}
model3 = aov(Density~Location, data=MarshGrassOct2021)
summary(model3)
```
```{r}
model4 = aov(Height~Location, data=MarshGrassSep2020)
summary(model4)
```
```{r}
model5 = aov(Height~Location, data=MarshGrassJul2021)
summary(model5)
```
```{r}
model6 = aov(Height~Location, data=MarshGrassOct2021)
summary(model6)
```
```{r}
t.test(Density~Location2,data=MarshGrassSep2020,alternative="two.sided",var.equal=TRUE)
```
```{r}
t.test(Density~Location2,data=MarshGrassJul2021,alternative="two.sided",var.equal=TRUE)
```
```{r}
t.test(Density~Location2,data=MarshGrassOct2021,alternative="two.sided",var.equal=TRUE)
```
```{r}
t.test(Height~Location2,data=MarshGrassSep2020,alternative="two.sided",var.equal=TRUE)
```
```{r}
t.test(Height~Location2,data=MarshGrassJul2021,alternative="two.sided",var.equal=TRUE)
```
```{r}
t.test(Height~Location2,data=MarshGrassOct2021,alternative="two.sided",var.equal=TRUE)
```
#no significant difference between edge and interior density and height for both sites 
#significant difference between heights at S4 and S7 in Sep and July, but not Oct 
#significant difference between densities at S4 and S7 in Sep and July, but not Oct 
```{r}
#Conducting Tukey HSD test - all possible pairwise comparisons for locations (S4E, S4I, S7E, S7I)
tukey.test <- TukeyHSD(model1)
tukey.test
#sig difference in density between S4I/S4E, S7E/S4E, S7I/S4E, S7E/S4I
```
```{r}
tukey.test <- TukeyHSD(model3)
tukey.test
#sig difference in density between S4I/S4E, S7E/S4E
```
```{r}
tukey.test <- TukeyHSD(model4)
tukey.test
#sig difference in height between S4I/S4E, S7E/S4E, S7I/S4E, S7E/S4I, S7I/S7E
```
```{r}
tukey.test <- TukeyHSD(model5)
tukey.test
#sig difference in height between S7I/S4E
```
```{r}
mean(MarshGrassSep2020$Height)
mean(MarshGrassJul2021$Height)
mean(MarshGrassOct2021$Height)
mean(MarshGrassSep2020$Density)
mean(MarshGrassJul2021$Density)
mean(MarshGrassOct2021$Density)
```
```{r}
Sep20HeightS4 = subset(MarshGrassSep2020$Height,MarshGrassSep2020$Site == "S4")
Sep20HeightS4 = as.data.frame(Sep20HeightS4)
mean(Sep20HeightS4$Sep20HeightS4)

Sep20HeightS7 = subset(MarshGrassSep2020$Height,MarshGrassSep2020$Site == "S7")
Sep20HeightS7 = as.data.frame(Sep20HeightS7)
mean(Sep20HeightS7$Sep20HeightS7)

Sep20DensityS4 = subset(MarshGrassSep2020$Density,MarshGrassSep2020$Site == "S4")
Sep20DensityS4 = as.data.frame(Sep20DensityS4)
mean(Sep20DensityS4$Sep20DensityS4)

Sep20DensityS7 = subset(MarshGrassSep2020$Density,MarshGrassSep2020$Site == "S7")
Sep20DensityS7 = as.data.frame(Sep20DensityS7)
mean(Sep20DensityS7$Sep20DensityS7)
```
```{r}
Jul21HeightS4 = subset(MarshGrassJul2021$Height,MarshGrassJul2021$Site == "S4")
Jul21HeightS4 = as.data.frame(Jul21HeightS4)
mean(Jul21HeightS4$Jul21HeightS4)

Jul21HeightS7 = subset(MarshGrassJul2021$Height,MarshGrassJul2021$Site == "S7")
Jul21HeightS7 = as.data.frame(Jul21HeightS7)
mean(Jul21HeightS7$Jul21HeightS7)

Jul21DensityS4 = subset(MarshGrassJul2021$Density,MarshGrassJul2021$Site == "S4")
Jul21DensityS4 = as.data.frame(Jul21DensityS4)
mean(Jul21DensityS4$Jul21DensityS4)

Jul21DensityS7 = subset(MarshGrassJul2021$Density,MarshGrassJul2021$Site == "S7")
Jul21DensityS7 = as.data.frame(Jul21DensityS7)
mean(Jul21DensityS7$Jul21DensityS7)
```
```{r}
Oct21HeightS4 = subset(MarshGrassOct2021$Height,MarshGrassOct2021$Site == "S4")
Oct21HeightS4 = as.data.frame(Oct21HeightS4)
mean(Oct21HeightS4$Oct21HeightS4)

Oct21HeightS7 = subset(MarshGrassOct2021$Height,MarshGrassOct2021$Site == "S7")
Oct21HeightS7 = as.data.frame(Oct21HeightS7)
mean(Oct21HeightS7$Oct21HeightS7)

Oct21DensityS4 = subset(MarshGrassOct2021$Density,MarshGrassOct2021$Site == "S4")
Oct21DensityS4 = as.data.frame(Oct21DensityS4)
mean(Oct21DensityS4$Oct21DensityS4)

Oct21DensityS7 = subset(MarshGrassOct2021$Density,MarshGrassOct2021$Site == "S7")
Oct21DensityS7 = as.data.frame(Oct21DensityS7)
mean(Oct21DensityS7$Oct21DensityS7)
```
```{r}
model1_resids = model1$residuals
model2_resids = model2$residuals
model3_resids = model3$residuals
model4_resids = model4$residuals
model5_resids = model5$residuals
model6_resids = model6$residuals

#Shapiro-Wilk's Test for normality
shapiro.test(model1_resids) #p value > 0.05, so points are normally distributed 
shapiro.test(model2_resids) #p value > 0.05, so points are normally distributed 
shapiro.test(model3_resids) #p value > 0.05, so points are normally distributed 
shapiro.test(model4_resids) #p value > 0.05, so points are normally distributed 
shapiro.test(model5_resids) #p value > 0.05, so points are normally distributed 
shapiro.test(model6_resids) #p value > 0.05, so points are normally distributed 
```
```{r}
model1 = aov(Height~Site, data=all)
summary(model1)
```


```{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 NOT normally distributed 
leveneTest(all$Height, all$Location) #p value < 0.05, so there is NO homogeneity of variation
```
```{r}
model2 = aov(Density~Site, data=all)
summary(model2)
```
```{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(all$Height, all$Location) #p value < 0.05, so there is NO homogeneity of variation
```
```{r}
anova1 = aov(Density~Site, data=MarshGrassJul2021)
summary(anova1)
```
```{r}
anova2 = aov(Height~Site, data=MarshGrassJul2021)
summary(anova2)
```
```{r}
anova1_resids <- anova1$residuals
qqPlot(anova1_resids) 
hist(anova1_resids)
#Shapiro Wilk's Test for Normality
shapiro.test(anova1_resids) #p > 0.05, so points are normally distributed 
leveneTest(MarshGrassJul2021$Density, MarshGrassJul2021$Site) #p value > 0.05, so there is homogeneity of variation
```
```{r}
anova2_resids <- anova1$residuals
qqPlot(anova2_resids) 
hist(anova2_resids)
#Shapiro Wilk's Test for Normality
shapiro.test(anova2_resids) #p > 0.05, so points are normally distributed 
leveneTest(MarshGrassJul2021$Height, MarshGrassJul2021$Site) #p value > 0.05, so there is homogeneity of variation
```
```{r}
anova3 = aov(Height~Location2, data=MarshGrassJul2021)
summary(anova3)
```
```{r}
library(psych)
describeBy(MarshGrassJul2021$Height, group= MarshGrassJul2021$Site)
```
```{r}
describeBy(MarshGrassJul2021$`Density(percm2)`, group= MarshGrassJul2021$Site)
```

