### FMV17_ElevChange_statistics.r

# Roeland C. van de Vijsel
# NIOZ Yerseke, The Netherlands
# www.nioz.nl
# V1.0 - 10/07/2020
# V2.0 - 21/09/2021

# In this script, statistical analyses are applied to the data from the "Fluid Mud Vaucheria 2017" experiment (FMV17).
# Goal of experiment: test whether Vaucheria algae increase the cohesion of fluid mud and whether this in turn increases sediment elevation
# Focus of this script: elevation change
# V2.0: two-way ANOVA (instead of one-way)

### 1.0 Set up script. OK!

remove(list=ls()) # Remove all variables from memory
graphics.off()    # Closes all figures 
cat("\014")       # Clears the console 

setwd("[WRITE HERE THE PATH LOCATION OF THE CURRENT R-SCRIPT]")
workdir = getwd()

### 1.1 Sort FMV17 data. OK!

FMV <- read.csv("FMV17_ElevChange_outliers.csv", header=T) # Read raw data set
head(FMV) # Show header of data set
levels(FMV$Time)      # Show the different Time treatments
levels(FMV$Vaucheria) # Show the different Vaucheria treatments
levels(FMV$Treatment) # Show all combinations of Time and Vaucheria treatments

##########################################################################
# Elevation Change (ElevChange)
##########################################################################

### 2.1 ElevChange: overview plot and summary table. OK!

png('output/ElevChange_boxplot.png')
par(cex=0.8) # Prepare plot
plot(ElevChange_mL ~ Treatment, data=FMV) # Simple boxplot
dev.off()

png('output/ElevChange_ORIGINAL_boxplot.png')
par(cex=0.8) # Prepare plot
plot(ElevChange_mL_ORIGINAL ~ Treatment, data=FMV) # Simple boxplot
dev.off()

sem <- function(x) { sd(x)/sqrt(length(x)) } # Define standard error of the mean
library(doBy)
options(digits=7)

## Transformed data
# ElevChangesum <- summaryBy(ElevChange_mL ~ Treatment, data = FMV, FUN=c(length, mean, sd, sem, min, max)) # Give data summary (one-way anova)
ElevChangesum <- summaryBy(ElevChange_mL ~ Time*Vaucheria, data = FMV, FUN=c(length, mean, sd, sem, min, max)) # Give data summary (two-way anova)
ElevChangesum
write.csv(ElevChangesum,"output/ElevChange_summary.csv")     

# ElevChangesumORIG <- summaryBy(ElevChange_mL_ORIGINAL ~ Treatment, data = FMV, FUN=c(length, mean, sd, sem, min, max)) # Give data summary (one-way anova)
ElevChangesumORIG <- summaryBy(ElevChange_mL_ORIGINAL ~ Time*Vaucheria, data = FMV, FUN=c(length, mean, sd, sem, min, max)) # Give data summary (two-way anova)
ElevChangesumORIG
write.csv(ElevChangesumORIG,"output/ElevChange_ORIGINAL_summary.csv")

### 2.2 ElevChange: ANOVA. OK!

# lm.FMV <- lm(ElevChange_mL ~ Treatment, data=FMV) # (one-way anova)
lm.FMV <- lm(ElevChange_mL ~ Time*Vaucheria, data=FMV) # (two-way anova)
ElevChangeanova <- anova(lm.FMV)
ElevChangeanova
write.csv(ElevChangeanova,"output/ElevChange_anova.csv")

### 2.3: ElevChange: Post-hoc tests. OK!

attach(FMV)
# (ElevChangethsd <- TukeyHSD(aov(ElevChange_mL ~ Treatment, data=FMV))) # (one-way anova)
# write.csv(ElevChangethsd$Treatment,"output/ElevChange_tHSD.csv") # (one-way anova)
(ElevChangethsd <- TukeyHSD(aov(ElevChange_mL ~ Time*Vaucheria, data=FMV))) # (two-way anova)
write.csv(ElevChangethsd$Time,"output/ElevChange_tHSD_Time.csv")
write.csv(ElevChangethsd$Vaucheria,"output/ElevChange_tHSD_Vaucheria.csv")
write.csv(ElevChangethsd$`Time:Vaucheria`,"output/ElevChange_tHSD_Time+Vaucheria.csv")

png('output/ElevChange_tHSDplot.png')
par(mfrow=c(1,1),mar=c(3.5, 9, 2.5, 1) + 0.1, mgp=c(2.5,1.5,0),las=1,cex=0.8)
plot(ElevChangethsd)
dev.off()

### 2.4: ElevChange: Test assumption of equal variances. OK!

library(car)
# ElevChangelevene <- leveneTest(ElevChange_mL ~ Treatment, data=FMV) # (one-way anova)
ElevChangelevene <- leveneTest(ElevChange_mL ~ Time*Vaucheria, data=FMV) # (two-way anova)
ElevChangelevene
write.csv(ElevChangelevene,"output/ElevChange_levene.csv")
### ==> H0 (equal variances) not rejected ==> Assumption not invalid. OK!

### 2.5: ElevChange: Test assumption of normality. OK!

png('output/ElevChange_qq.png')
# lm.FMV <- lm(ElevChange_mL ~ Treatment, data=FMV) # (one-way anova)
lm.FMV <- lm(ElevChange_mL ~ Time*Vaucheria, data=FMV) # (two-way anova)
res <- resid(lm.FMV)
qqPlot(res)
dev.off()
### ==> Some outliers, but assumption of normality does not seem invalid ==> Perform Shapiro-Wilk test to confirm.

ElevChangeshapiro <- print(shapiro.test(res))
write.csv(ElevChangeshapiro$p.value,"output/ElevChange_shapiro_pvalue.csv")
### ==> H0 (normal distribution) not rejected ==> Assumption not invalid

png('output/ElevChange_hist.png')
hist(res)
dev.off()
### Histogram appears quite normal.