### FMV17_statistics.r

# Roeland C. van de Vijsel
# NIOZ Yerseke, The Netherlands
# www.nioz.nl
# V1.0 - 01/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: Relative Area Increase of loadstep 4 (fully loaded), i.e. deformation ratio.
# Changes in V2.0: Apply two-way anova instead of one-way anova

### 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.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

##########################################################################
# Relative Area Increase 4 (RAI4)
##########################################################################

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

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

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

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

### 2.2 RAI4: ANOVA. OK!

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

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

attach(FMV)
# (RAI4thsd <- TukeyHSD(aov(RelAreaIncr4 ~ Treatment, data=FMV))) # (one-way anova)
# write.csv(RAI4thsd$Treatment,"output/RAI4_tHSD.csv") # (one-way anova)
(RAI4thsd <- TukeyHSD(aov(RelAreaIncr4 ~ Time*Vaucheria, data=FMV))) # (two-way anova)
write.csv(RAI4thsd$Time,"output/RAI4_tHSD_Time.csv") # (two-way anova)
write.csv(RAI4thsd$Vaucheria,"output/RAI4_tHSD_Vaucheria.csv") # (two-way anova)
write.csv(RAI4thsd$`Time:Vaucheria`,"output/RAI4_tHSD_Time+Vaucheria.csv") # (two-way anova)

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

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

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

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

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

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

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