## Sterre Witte June 2023

## SAND COVERAGE

##---------------------------------
## Load packages

library(tidyverse)
library(dplyr)
library(lme4)
library(car)
library(emmeans)

##---------------------------------
## Load data

setwd("~...folder containing data") 
data <- read.csv("Sediment data.csv", sep=";")

##---------------------------------
## Filter & summarise

data <- data %>%
  mutate_if(is.character, as.factor) %>%
  mutate(Block=as.factor(Block)) %>%
  filter(Season=="Fall 2020" | Season=="Fall 2021")

year <- data %>%
  group_by(Season) %>%
  summarise(mean = mean(Percentage_sanded_over, na.rm=T), SD=sd(Percentage_sanded_over, na.rm=T))

(year[year$Season=="Fall 2021",]$mean-year[year$Season=="Fall 2020",]$mean)/
  year[year$Season=="Fall 2020",]$mean*100

block <- data %>%
  filter(Season=="Fall 2021")%>%
  group_by(Block)%>%
  summarise(min=min(Percentage_sanded_over, na.rm = T), max=max(Percentage_sanded_over, na.rm=T))

ggplot(data, aes(x=Season, y=Percentage_sanded_over)) + 
  geom_boxplot()

peryear <- glmer(data=data, Percentage_sanded_over~Season+(1|Block), family=poisson)
peryearn <- glm(data=data, Percentage_sanded_over~Season, family=poisson)
AIC(peryear)-AIC(peryearn)
Anova(peryear)
summary(peryear)

substrate <- data %>%
  group_by(Substrate) %>%
  summarise(mean = mean(Percentage_sanded_over, na.rm=T), SD=sd(Percentage_sanded_over, na.rm=T))

(substrate[substrate$Substrate=="BESE",]$mean-mean(substrate[substrate$Substrate!="BESE",]$mean))/
  mean(substrate[substrate$Substrate!="BESE",]$mean)*100

(mean(substrate[substrate$Substrate=="Cobbles" | substrate$Substrate=="Pebbles" | substrate$Substrate=="Shell",]$mean)
  -mean(substrate[substrate$Substrate=="Reef" | substrate$Substrate=="Wood",]$mean))/
  mean(substrate[substrate$Substrate=="Reef" | substrate$Substrate=="Wood",]$mean)*100

ggplot(data, aes(x=Substrate, y=Percentage_sanded_over)) + 
  geom_boxplot()+
  geom_point()

persubstrate <- glmer(data=data, Percentage_sanded_over~Substrate+(1|Block), family=poisson)
persubstraten <- glm(data=data, Percentage_sanded_over~Substrate, family=poisson)
AIC(persubstrate)-AIC(persubstraten)
Anova(persubstrate)
emmeans(persubstrate, pairwise~Substrate, adjust='fdr')
summary(persubstrate)

block <- data %>%
  group_by(Block) %>%
  summarise(mean = mean(Percentage_sanded_over, na.rm=T), SD=sd(Percentage_sanded_over, na.rm=T),
            min=min(Percentage_sanded_over, na.rm=T), max=max(Percentage_sanded_over, na.rm=T))