### VMT18_statistics.r

# Roeland C. van de Vijsel
# NIOZ Yerseke, The Netherlands
# www.nioz.nl
# V1.0 - 30/06/2020

# In this script, statistical analyses are applied to the data from the "Vaucheria Micro-Topography 2018" experiment (VMT18).
# Goal of experiment: test whether Vaucheria algae grow better on better-drained sedimentary substrates
# Four treatments, i.e. different types of substrate: 1) Hollow, 2) Flat, 3) Hummock, 4) Control. 
# The Control is a positive control, to test whether Vaucheria grows well enough at all.
# Algal growth is quantified via measurement of chlorophyll-a content and organic carbon content.
# Setup of script: i) Plot data; ii) ANOVA; iii) Pairwise-comparison; iv) Test homoscedasticity assumption; v) Test normality assumption

### 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 FULL FILE PATH OF THE CURRENT R-SCRIPT]") # Set working directory
workdir = getwd()

### 1.1 Sort VMT18 data. OK!

VMT <- read.csv("VMT18.csv", header=T) # Read raw data set
head(VMT) # Show header of data set
levels(VMT$Substrate) # Show the different groups (treatments) that are defined

##########################################################################
# Chlorophyl-a
##########################################################################

### 2.1 CHLA: plot. OK!

pdf('output/CHLA_boxplot.pdf')
par(cex=0.8) # Prepare plot
plot(CHLA_total_mug ~ Substrate, data=VMT) # Simple boxplot
dev.off()

sem <- function(x) { sd(x)/sqrt(length(x)) } # Define standard error of the mean
library(doBy)
options(digits=7)
CHLAsum <- summaryBy(CHLA_total_mug ~ Substrate, data = VMT, FUN=c(length, mean, sd, sem, min, max)) # Give data summary
CHLAsum
write.csv(CHLAsum,"output/CHLA_summary.csv")     

### 2.2 CHLA: ANOVA. OK!

lm.VMT <- lm(CHLA_total_mug ~ Substrate, data=VMT)
CHLAanova <- anova(lm.VMT)
CHLAanova
write.csv(CHLAanova,"output/CHLA_anova.csv")

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

attach(VMT)
(CHLAthsd <- TukeyHSD(aov(CHLA_total_mug ~ Substrate, data=VMT)))
write.csv(CHLAthsd$Substrate,"output/CHLA_tHSD.csv")

pdf('output/CHLA_tHSDplot.pdf')
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. OK!

library(car)
CHLAlevene <- leveneTest(CHLA_total_mug ~ Substrate, data=VMT)
CHLAlevene
write.csv(CHLAlevene,"output/CHLA_levene.csv")
### ==> H0 (equal variances) not rejected ==> Assumption not invalid.

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

pdf('output/CHLA_qq.pdf')
lm.VMT <- lm(CHLA_total_mug ~ Substrate, data=VMT)
res <- resid(lm.VMT)
qqPlot(res)
dev.off()
### ==> Assumption of normality seems valid

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

pdf('output/CHLA_hist.pdf')
hist(res)
dev.off()
### Histogram appears quite normal.

##########################################################################
# Organic Carbon
##########################################################################

### 3.1 orgC: plot. OK!

pdf('output/orgC_boxplot.pdf')
par(cex=0.8) # Prepare plot
plot(orgC_total_g ~ Substrate, data=VMT) # Simple boxplot
dev.off()

sem <- function(x) { sd(x)/sqrt(length(x)) } # Define standard error of the mean
library(doBy)
options(digits=7)
orgCsum <- summaryBy(orgC_total_g ~ Substrate, data = VMT, FUN=c(length, mean, sd, sem, min, max)) # Give data summary
orgCsum
write.csv(orgCsum,"output/orgC_summary.csv")     

### 3.2 orgC: ANOVA. OK!

lm.VMT <- lm(orgC_total_g ~ Substrate, data=VMT)
orgCanova <- anova(lm.VMT)
orgCanova
write.csv(orgCanova,"output/orgC_anova.csv")

### 3.3: orgC Post-hoc tests. OK!

attach(VMT)
(orgCthsd <- TukeyHSD(aov(orgC_total_g ~ Substrate, data=VMT)))
write.csv(orgCthsd$Substrate,"output/orgC_tHSD.csv")

pdf('output/orgC_tHSDplot.pdf')
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(orgCthsd)
dev.off()

### 3.4: orgC Test assumption of equal variances. OK!

library(car)
orgClevene <- leveneTest(orgC_total_g ~ Substrate, data=VMT)
orgClevene
write.csv(orgClevene,"output/orgC_levene.csv")
### ==> H0 (equal variances) not rejected ==> Assumption not invalid.

### 3.5: orgC Test assumption of normality. OK!

pdf('output/orgC_qq.pdf')
lm.VMT <- lm(orgC_total_g ~ Substrate, data=VMT)
res <- resid(lm.VMT)
qqPlot(res)
dev.off()
### ==> Assumption of normality seems valid

orgCshapiro <- print(shapiro.test(res))
write.csv(orgCshapiro$p.value,"output/orgC_shapiro_pvalue.csv")
### ==> H0 (normal distribution) rejected when tested formally ...

pdf('output/orgC_hist.pdf')
hist(res)
dev.off()
### Histogram appears quite normal; porbably sample size is small.


##########################################################################
# 4. Organic Carbon - Logtransformed analysis ( orgC_total_g_LOG = log10(orgC_total_g) )
##########################################################################

### 4.1 orgC: plot. OK!

pdf('output/orgClog_boxplot.pdf')
par(cex=0.8) # Prepare plot
plot(orgC_total_g_LOG ~ Substrate, data=VMT) # Simple boxplot
dev.off()

sem <- function(x) { sd(x)/sqrt(length(x)) } # Define standard error of the mean
library(doBy)
options(digits=7)
orgCsum <- summaryBy(orgC_total_g_LOG ~ Substrate, data = VMT, FUN=c(length, mean, sd, sem, min, max)) # Give data summary
orgCsum
write.csv(orgCsum,"output/orgClog_summary.csv")     

### 4.2 orgC: ANOVA. OK!

lm.VMT <- lm(orgC_total_g_LOG ~ Substrate, data=VMT)
orgCanova <- anova(lm.VMT)
orgCanova
write.csv(orgCanova,"output/orgClog_anova.csv")

### 4.3: orgC Post-hoc tests. OK!

attach(VMT)
(orgCthsd <- TukeyHSD(aov(orgC_total_g_LOG ~ Substrate, data=VMT)))
write.csv(orgCthsd$Substrate,"output/orgClog_tHSD.csv")

pdf('output/orgClog_tHSDplot.pdf')
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(orgCthsd)
dev.off()

### 4.4: orgC Test assumption of equal variances. OK!

library(car)
orgClevene <- leveneTest(orgC_total_g_LOG ~ Substrate, data=VMT)
orgClevene
write.csv(orgClevene,"output/orgClog_levene.csv")
### ==> H0 (equal variances) not rejected ==> Assumption not invalid.

### 4.5: orgC Test assumption of normality. OK!

pdf('output/orgClog_qq.pdf')
lm.VMT <- lm(orgC_total_g_LOG ~ Substrate, data=VMT)
res <- resid(lm.VMT)
qqPlot(res)
dev.off()
### ==> Assumption of normality seems valid

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

pdf('output/orgClog_hist.pdf')
hist(res)
dev.off()
### Histogram appears quite normal; porbably sample size is small.