### FMV17_SAMPLING_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: Water content
# 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 THE PATH LOCATION OF THE CURRENT R-SCRIPT HERE]")
workdir = getwd()

### 1.1 Sort FMV17 data. OK!

FMV <- read.csv("FMV17_SAMPLING_rawdata_duplicateT0.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

##########################################################################
# Water Content (WC)
##########################################################################

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

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

# png('output/WC_ORIGINAL_boxplot.png')
# par(cex=0.8) # Prepare plot
# plot(WC_perc_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
# WCsum <- summaryBy(WC_perc ~ Treatment, data = FMV, FUN=c(length, mean, sd, sem, min, max)) # Give data summary (one-way anova)
WCsum <- summaryBy(WC_perc ~ Time*Vaucheria, data = FMV, FUN=c(length, mean, sd, sem, min, max)) # Give data summary
WCsum
write.csv(WCsum,"output/WC_summary.csv")     

## Original (untransformed) data
# WCsumORIG <- summaryBy(WC_perc_ORIGINAL ~ Treatment, data = FMV, FUN=c(length, mean, sd, sem, min, max)) # Give data summary (one-way anova)
WCsumORIG <- summaryBy(WC_perc_ORIGINAL ~ Time*Vaucheria, data = FMV, FUN=c(length, mean, sd, sem, min, max)) # Give data summary
WCsumORIG
write.csv(WCsumORIG,"output/WC_ORIGINAL_summary.csv")

### 2.2 WC: ANOVA. OK!

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

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

attach(FMV)
# (WCthsd <- TukeyHSD(aov(WC_perc ~ Treatment, data=FMV))) # (one-way anova)
# write.csv(WCthsd$Treatment,"output/WC_tHSD.csv") # (one-way anova)
(WCthsd <- TukeyHSD(aov(WC_perc ~ Time*Vaucheria, data=FMV))) # (two-way anova)
write.csv(WCthsd$Time,"output/WC_tHSD_Time.csv") # (two-way anova)
write.csv(WCthsd$Vaucheria,"output/WC_tHSD_Vaucheria.csv") # (two-way anova)
write.csv(WCthsd$`Time:Vaucheria`,"output/WC_tHSD_Time+Vaucheria.csv") # (two-way anova)

png('output/WC_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(WCthsd)
dev.off()

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

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

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

png('output/WC_qq.png')
# lm.FMV <- lm(WC_perc ~ Treatment, data=FMV) # (one-way anova)
lm.FMV <- lm(WC_perc ~ 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.

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

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