#set working directory to file location
setwd("C:/Users/elias/Documents/R/Internship_2")

#install pakages necessary for scripts to run and load the libraries
#install.packages("vegan")     #for PERMANOVA
#install.packages("readxl")    #to read excell files
#install.packages("devtools")  #to allow other packages to run
library(devtools)             
#install.packages('cluster')   #for pairwise comparisons
#install_github("pmartinezarbizu/pairwiseAdonis/pairwiseAdonis")  #for pairwise comparisons, this packages can give errors, make sure to disable errors from warnings otherwise installarion will fail
          #when errors occur use this piece of code: Sys.setenv("R_REMOTES_NO_ERRORS_FROM_WARNINGS"=TRUE)
library(vegan)
library(readxl)
library(cluster)
library(pairwiseAdonis)

#WARINGING: NEVER RUN BOTH FIRST AND SECOND PART OF THIS SCRIPT AT THE SAME TIME!!!


#FIRST PART:
#First part is the permanova based on the full database, this is including all timepoints

#create the factor file that permanova needs and check if names are imported correctly
factor <- read_xlsx("Edited_Data/factor.xlsx") 
names(factor)


#load the exell files and check if names are imported correctly
coral <- read_xlsx("Edited_Data/permanovacoral.xlsx")      #datasheet with square root transformed cover percentages for all coral genera for all timepoints
names(coral)
benthos <- read_xlsx("Edited_Data/permanovabenthos.xlsx")      #datasheet with square root transformed cover percentages for benthic cover for all timepoints
names(benthos)
lhs <- read_xlsx("Edited_Data/permanovalhs_updated.xlsx")      #datasheet with square root transformed cover percentages for all life history strategies for all timepoints
names(lhs)
lhs_ex_poc <- read_xlsx("Edited_Data/permanovalhs_poc_excluded.xlsx")      #datasheet with square root transformed cover percentages for all life history strategies for all timepoints with pocillopora excluded
names(lhs_ex_poc)
corallite <- read_xlsx("Edited_Data/permanovaCIS.xlsx")    #datasheet with square root transformed cover percentages for high and low Corallite Integration Scores for all timepoints
names(corallite)
growth <- read_xlsx("Edited_Data/permanovagrowthform.xlsx")    #datasheet with square root transformed cover percentages for all growth forms for all timepoints
names(growth)


#two way permanova non scaled data
adonis (coral ~ factor$Time * factor$Zone, method = "bray", perm=999)
adonis (benthos ~ factor$Time * factor$Zone, method = "bray", perm=999)
adonis (lhs ~ factor$Time * factor$Zone, method = "bray", perm=999)
adonis (lhs_ex_poc ~ factor$Time * factor$Zone, method = "bray", perm=999)
adonis (corallite ~ factor$Time * factor$Zone, method = "bray", perm=999)
adonis (growth ~ factor$Time * factor$Zone, method = "bray", perm=999)


#pairwise test on two way permanova
pairwise.adonis(x=coral, factors=factor$TimeZone, sim.method='bray',p.adjust.m='holm')
pairwise.adonis(x=benthos, factors=factor$TimeZone, sim.method='bray',p.adjust.m='holm')
pairwise.adonis(x=lhs, factors=factor$TimeZone, sim.method='bray',p.adjust.m='holm')
pairwise.adonis(x=lhs_ex_poc, factors=factor$TimeZone, sim.method='bray',p.adjust.m='holm')
pairwise.adonis(x=corallite, factors=factor$TimeZone, sim.method='bray',p.adjust.m='holm')
pairwise.adonis(x=growth, factors=factor$TimeZone, sim.method='bray',p.adjust.m='holm')


#betadisp
dst_coral <- dist(coral)
coral.bd <- betadisper(dst_coral, factor$TimeZone)
coral.bd
permutest(coral.bd)

dst_benthos <- dist(benthos)
benthos.bd <- betadisper(dst_benthos, factor$TimeZone)
benthos.bd
permutest(benthos.bd)

dst_lhs <- dist(lhs)
lhs.bd <- betadisper(dst_lhs, factor$TimeZone)
lhs.bd
permutest(lhs.bd)

dst_lhs_ex_poc <- dist(lhs_ex_poc)
lhs.bd_ex_poc <- betadisper(dst_lhs_ex_poc, factor$TimeZone)
lhs.bd_ex_poc
permutest(lhs.bd_ex_poc)

dst_corallite <- dist(corallite)
corallite.bd <- betadisper(dst_corallite, factor$TimeZone)
corallite.bd
permutest(corallite.bd)

dst_growth <- dist(growth)
growth.bd <- betadisper(dst_growth, factor$TimeZone)
growth.bd
permutest(growth.bd)
#WHEN RUNNING FIRST PART STOP HERE

#To continue on with the second part first clear the workspace to prevent any interference
rm(list=ls())

#SECOND PART:

#The second part is the permanova based on only the january 2016 (pre-bleaching) and the october 2019 (long-term recovery) timepoint

#create the factor file that permanova needs and check if names are imported correctly
factor <- read_xlsx("Edited_Data/factor_JAN2019.xlsx") 
names(factor)

#load the exell files and check if names are imported correctly
coral <- read_xlsx("Edited_Data/permanovacoral_JAN2019.xlsx")     #datasheet with square root transformed cover percentages for all coral genera for january 2016 versus october 2019 timepoints
names(coral)
benthos <- read_xlsx("Edited_Data/permanovabenthos_JAN2019.xlsx")     #datasheet with square root transformed cover percentages for benthic cover for january 2016 versus october 2019 timepoints
names(benthos)
lhs <- read_xlsx("Edited_Data/permanovalhs_JAN2019_updated.xlsx")      #datasheet with square root transformed cover percentages for all life history strategies for all timepoints
names(lhs)
lhs_ex_poc <- read_xlsx("Edited_Data/permanovalhs_JAN2019_poc_excluded.xlsx")      #datasheet with square root transformed cover percentages for all life history strategies for all timepoints with pocillopora excluded
names(lhs_ex_poc)
corallite <- read_xlsx("Edited_Data/permanovaCIS_JAN2019.xlsx")     #datasheet with square root transformed cover percentages for high and low Corallite Integration Scores for january 2016 versus october 2019 timepoints
names(corallite)
growth <- read_xlsx("Edited_Data/permanovagrowthform_JAN2019.xlsx")     #datasheet with square root transformed cover percentages for all growth forms for january 2016 versus october 2019 timepoints
names(growth)



#two way permanova non scaled data
adonis (coral ~ factor$Time * factor$Zone, method = "bray", perm=999)
adonis (benthos ~ factor$Time * factor$Zone, method = "bray", perm=999)
adonis (lhs ~ factor$Time * factor$Zone, method = "bray", perm=999)
adonis (lhs_ex_poc ~ factor$Time * factor$Zone, method = "bray", perm=999)
adonis (corallite ~ factor$Time * factor$Zone, method = "bray", perm=999)
adonis (growth ~ factor$Time * factor$Zone, method = "bray", perm=999)

#pairwise test on two way permanova
pairwise.adonis(x=coral, factors=factor$TimeZone, sim.method='bray',p.adjust.m='holm')
pairwise.adonis(x=benthos, factors=factor$TimeZone, sim.method='bray',p.adjust.m='holm')
pairwise.adonis(x=lhs, factors=factor$TimeZone, sim.method='bray',p.adjust.m='holm')
pairwise.adonis(x=lhs_ex_poc, factors=factor$TimeZone, sim.method='bray',p.adjust.m='holm')
pairwise.adonis(x=corallite, factors=factor$TimeZone, sim.method='bray',p.adjust.m='holm')
pairwise.adonis(x=growth, factors=factor$TimeZone, sim.method='bray',p.adjust.m='holm')

#betadisp
dst_coral <- dist(coral)
coral.bd <- betadisper(dst_coral, factor$TimeZone)
coral.bd
permutest(coral.bd)

dst_benthos <- dist(benthos)
benthos.bd <- betadisper(dst_benthos, factor$TimeZone)
benthos.bd
permutest(benthos.bd)

dst_lhs <- dist(lhs)
lhs.bd <- betadisper(dst_lhs, factor$TimeZone)
lhs.bd
permutest(lhs.bd)

dst_lhs_ex_poc <- dist(lhs_ex_poc)
lhs.bd_ex_poc <- betadisper(dst_lhs_ex_poc, factor$TimeZone)
lhs.bd_ex_poc
permutest(lhs.bd_ex_poc)

dst_corallite <- dist(corallite)
corallite.bd <- betadisper(dst_corallite, factor$TimeZone)
corallite.bd
permutest(corallite.bd)

dst_growth <- dist(growth)
growth.bd <- betadisper(dst_growth, factor$TimeZone)
growth.bd
permutest(growth.bd)