# Coded by Alastair Potts
# For review

library(groundhog)
set.groundhog.folder("C:\\Groundhog_R\\R4.2.3_2023-04-01")
GroundhogDay <- '2023-04-01'
# set.groundhog.folder("C:\\Groundhog_R\\R4.2.3_2024-02-01")
# GroundhogDay <- '2024-02-01'

groundhog.library("dplyr",GroundhogDay)
groundhog.library("ggplot2",GroundhogDay)
groundhog.library("tibble",GroundhogDay)
groundhog.library("tidyr",GroundhogDay)
groundhog.library("lubridate",GroundhogDay)
groundhog.library("rstatix",GroundhogDay)
groundhog.library("patchwork",GroundhogDay)
groundhog.library("ggpmisc",GroundhogDay)
groundhog.library("multcompView",GroundhogDay)

setwd("D:\\Google Drive\\0_SpekboomResearchGroup2.0\\3_Experiments\\RegionalTransect21(10)_ParentPlantMaterialExperiment\\manuscript/")

read.csv("data_manuscript/harvestdata_2023-08-29.csv")%>%
  as_tibble()->
  dat_all


read.csv("data_manuscript/harvest_days.csv")%>%
  as_tibble()%>%
  mutate(Date=dmy(Date))%>%
  mutate(DaysSincePlanting=as.factor(DaysSincePlanting))->
  harvestdates


dat_all%>%
  mutate(Site=factor(Site,labels=paste(1:10)))->
  dat_all

dat_all%>%
  ungroup()%>%
  filter(SamplingEvent!=2)%>%
  filter(SamplingEvent!=3)%>%
  #filter(SamplingEvent!=4)%>% # Removing week 4 until sorted out
  filter(!is.na(RootMass))%>%
  left_join(harvestdates)%>%
  mutate(Individual=as.factor(Individual))->
  dat_rootmass

dat_rootmass%>%
    group_by(SamplingEvent,Site,Individual,DaysSinceStart)%>%
    summarise(uRootMass=mean(RootMass,na.rm=T))%>%
    left_join(harvestdates)->
    ind_uRootMass

### ANOVA ASSUMPTIONS ###
## Testing if normal distributions within each site is violated
# Harvest event 7 has normally distributed data for all sites
(ind_uRootMass%>%
  group_by(DaysSinceStart,SamplingEvent,Site)%>%
  shapiro_test(uRootMass)%>%
  filter(p<0.1)->
  non_normal_sites)
# Result: some sites have non-normal distributions within each harvest event (
# DaysSinceStart) 

##Testing for equal variances
# All events have equal variances
ind_uRootMass%>%
  group_by(SamplingEvent,DaysSinceStart)%>%
  levene_test(uRootMass~Site)


### ANOVA (on mean values) TEST ###

## Removing non-normal sites from each harvest event, and conducting an anova.
ind_uRootMass%>%
  anti_join(select(non_normal_sites,SamplingEvent,Site))%>%
  group_by(SamplingEvent)%>%
  anova_test(uRootMass~Site)

## Keeping the non-normal sites, and conducting an anova.
ind_uRootMass%>%
  #anti_join(select(non_normal_sites,SamplingEvent,Site))%>%
  df_group_by(SamplingEvent)%>%
  anova_test(uRootMass~Site)


### POST-HOC TUKEY TEST ###

# post-hoc with the non-normal distribution sites within each harvest event 
# removed.
(ind_uRootMass%>%
  anti_join(select(non_normal_sites,SamplingEvent,Site))%>%
  group_by(SamplingEvent,DaysSinceStart)%>%
  tukey_hsd(uRootMass~Site)->tukey_1)%>%
  filter(p.adj.signif!="ns")

# post-hoc with all sites (including those with non-normal distributions)
(ind_uRootMass%>%
#  anti_join(select(non_normal_sites,SamplingEvent,Site))%>%
  group_by(SamplingEvent,DaysSinceStart)%>%
  tukey_hsd(uRootMass~Site)->tukey_2)%>%
  filter(p.adj.signif!="ns")

left_join(tukey_1,tukey_2,by=c("SamplingEvent","group1","group2","DaysSinceStart"))%>%
  mutate(DaysSinceStart=factor(DaysSinceStart,labels=paste(unique(DaysSinceStart)," days")))%>%
  ggplot(aes(p.adj.x,p.adj.y))+
  geom_vline(xintercept=0.05,colour="darkgrey")+
  geom_hline(yintercept=0.05,colour="darkgrey")+
  geom_smooth(method="lm")+
  stat_poly_eq(aes(label =  sprintf("%s*\"; \"*%s",
                                    after_stat(rr.label),
                                    after_stat(p.value.label))),
               small.r=TRUE,
               small.p=TRUE)+
  geom_point()+
  facet_wrap(~DaysSinceStart)+
  ylab("Non-normal sites included")+
  xlab("Non-normal sites removed")+
  theme_bw()+
  coord_equal()
ggsave("FinalFiguresV2/Fig_s5.jpg",width=8,height=5)