### PaardenschorSeedlings_statistics.r

# Roeland C. van de Vijsel
# NIOZ Yerseke, The Netherlands
# www.nioz.nl
# V1 - 23 July 2020

# This R-script performs statistical analysis of seedling density data, for 
# the boxplot in Figure 4 in the PhD Thesis of Roeland C. van de Vijsel 
# (https://doi.org/10.33612/diss.160081233).

### 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("[FILL IN CURRENT WORKING DIRECTORY HERE]")
workdir = getwd()

### 1.1 Sort PaardenschorSeedlings data. OK!

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

##########################################################################
# SeedlingDensity [1/m2]
##########################################################################

### 2.1 SeedDens: plot. OK!

pdf('output/SeedDens_boxplot.pdf')
par(cex=0.8) # Prepare plot
plot(SeedlingDensity_per_m2 ~ Substrate, data=PSS) # Simple boxplot
dev.off()

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

### 2.2 SeedDens: ANOVA. OK!

lm.PSS <- lm(SeedlingDensity_per_m2 ~ Substrate, data=PSS)
SeedDensanova <- anova(lm.PSS)
SeedDensanova
write.csv(SeedDensanova,"output/SeedDens_anova.csv")

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

attach(PSS)
(SeedDensthsd <- TukeyHSD(aov(SeedlingDensity_per_m2 ~ Substrate, data=PSS)))
write.csv(SeedDensthsd$Substrate,"output/SeedDens_tHSD.csv")

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

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

library(car)
SeedDenslevene <- leveneTest(SeedlingDensity_per_m2 ~ Substrate, data=PSS)
SeedDenslevene
write.csv(SeedDenslevene,"output/SeedDens_levene.csv")
### ==> H0 (equal variances) not rejected ==> Assumption not invalid.

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

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

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

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

##

kruskal.test(SeedlingDensity_per_m2 ~ Substrate, data=PSS)
  
##
library("dunn.test")
DATA <- PSS$SeedlingDensity_per_m2
GROUPS <- PSS$Substrate

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