### FMV17_CHLA_statistics.r

# Roeland C. van de Vijsel
# NIOZ Yerseke, The Netherlands
# www.nioz.nl
# V1.0 - 10/07/2020
# V2.0 - 22/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: Chlorophyll-a 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("~/Documents/20181029_Documents_OLD/Articles/Writing/Article3_Experiments_BGM_feedbacks/v11_Revision1/REVISION1/Fig5+6_two-way_ANOVA/Figure6_R1/rawdata")
# setwd("[WRITE THE PATH LOCATION OF THE CURRENT R-SCRIPT HERE]")
workdir = getwd()

### 1.1 Sort FMV17 data. OK!

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

##########################################################################
# Chlorophyll-a content (CHLA)
##########################################################################

### 2.1 CHLA: plot. OK!

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

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

## Original (untransformed) data ==> NOT NEEDED, SINCE CHLA DATA IS NOT TRANSFORMED, SO TRANSFORMED = ORIGINAL
# CHLAsumORIG <- summaryBy(CHLA_mug_per_g_ORIGINAL ~ Treatment, data = FMV, FUN=c(length, mean, sd, sem, min, max)) # Give data summary (one-way anova)
# CHLAsumORIG <- summaryBy(CHLA_mug_per_g_ORIGINAL ~ Time*Vaucheria, data = FMV, FUN=c(length, mean, sd, sem, min, max)) # Give data summary (two-way anova)
# CHLAsumORIG
# write.csv(CHLAsumORIG,"output/CHLA_ORIGINAL_summary.csv")

### 2.2 CHLA: ANOVA. OK!

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

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

attach(FMV)
# (CHLAthsd <- TukeyHSD(aov(CHLA_mug_per_g ~ Treatment, data=FMV))) # (one-way anova)
# write.csv(CHLAthsd$Treatment,"output/CHLA_tHSD.csv") # (one-way anova)
(CHLAthsd <- TukeyHSD(aov(CHLA_mug_per_g ~ Time*Vaucheria, data=FMV))) # (two-way anova)
write.csv(CHLAthsd$Time,"output/CHLA_tHSD_Time.csv") # (two-way anova)
write.csv(CHLAthsd$Vaucheria,"output/CHLA_tHSD_Vaucheria.csv") # (two-way anova)
write.csv(CHLAthsd$`Time:Vaucheria`,"output/CHLA_tHSD_Time+Vaucheria.csv") # (two-way anova)

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

### 2.4: CHLA Test assumption of equal variances.

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

### 2.5: CHLA Test assumption of normality.

png('output/CHLA_qq.png')
# lm.FMV <- lm(CHLA_mug_per_g ~ Treatment, data=FMV) # (one-way anova)
lm.FMV <- lm(CHLA_mug_per_g ~ Time*Vaucheria, data=FMV) # (one-way anova)
res <- resid(lm.FMV)
qqPlot(res)
dev.off()
### ==> Assumption of normality doesn't seem valid!

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

png('output/CHLA_hist.png')
hist(res)
dev.off()
### Histogram does not appear normally distributed.

## Because non-normal and heteroscedastic: try non-parametric test instead (Kruskal-Wallis Rank Sum Test)
## Since (to our knowledge), non-parametric tests are only suitable for cases with one fixed factor, we use Treatment instead of Time*Vaucheria
FMV17_CHLA_Kruskal <- kruskal.test(CHLA_mug_per_g ~ Treatment, data=FMV) # (non-paramatric alternative for one-way anova)
FMV17_CHLA_Kruskal
write.csv(FMV17_CHLA_Kruskal$statistic,"output/CHLA_kruskal_chi2.csv")
write.csv(FMV17_CHLA_Kruskal$parameter,"output/CHLA_kruskal_df.csv")
write.csv(FMV17_CHLA_Kruskal$p.value,"output/CHLA_kruskal_pvalue.csv")

## Pairwise comparison: Dunn's test (also only for cases with one fixed factor)
library("dunn.test")
DATA <- FMV$CHLA_mug_per_g
GROUPS <- FMV$Treatment

dunn <- dunn.test(DATA, g=GROUPS, method="bonferroni",kw = TRUE, label = TRUE, wrap = TRUE, alpha = 0.05 )
dunn
write.csv(dunn,"output/CHLA_dunn.csv")