# 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 ###

# See Potts et al_PeerJ_Fig_S5.R for further post-hoc comparison




#############
# Main figure: site vs root mass

# Version 1: facet_wrapped(). .:. constant y-axis
ind_uRootMass%>%
  group_by(SamplingEvent,Site)%>%
  summarise(seRootMass=sd(uRootMass),
            uRootMass=mean(uRootMass),
            plusbar=uRootMass+seRootMass,
            minusbar=uRootMass-seRootMass,
            n=n()
  )->
  ind_uRootMass_stats

ind_uRootMass%>%
  ggplot(aes(Site,uRootMass))+
  geom_point(alpha=0.5,shape=4)+
  geom_linerange(aes(Site,ymin=minusbar,ymax=plusbar),data=ind_uRootMass_stats,linewidth=1)+
  geom_point(aes(Site,uRootMass),data=ind_uRootMass_stats)+
  facet_wrap(~SamplingEvent)

# Version 2
plotter <- function(i) {
  ind_uRootMass%>%
    filter(SamplingEvent==i)%>%
    group_by(Site)%>%
    summarise(sdRootMass=sd(uRootMass),
              seRootMass=sd(uRootMass)/sqrt(n()),
              uRootMass=mean(uRootMass),
              plusbar_sd=uRootMass+sdRootMass,
              minusbar_sd=uRootMass-sdRootMass,
              plusbar_se=uRootMass+seRootMass,
              minusbar_se=uRootMass-seRootMass,
              n=n()
    )->
    ind_uRootMass_stats
  
  ind_uRootMass%>%
      filter(SamplingEvent==i)%>%
      group_by(DaysSinceStart,SamplingEvent,Site)%>%
      shapiro_test(uRootMass)%>%
    mutate(not_norm=ifelse(p<=0.05,"*",""))->
    shap_norm_result
    
  ind_uRootMass%>%
    filter(SamplingEvent==i)%>%
    group_by(SamplingEvent,DaysSinceStart)%>%
    tukey_hsd(uRootMass~Site)-> post_hoc
  post_hoc$p.adj->ph_v
  names(ph_v) <- paste(post_hoc$group1,post_hoc$group2,sep="-")  
  cld<- multcompLetters(ph_v)  
  
  bind_cols(Letters=cld$Letters,
            Site=cld$Letters%>%names())->
    Tk
  Tk%>%
    left_join(select(shap_norm_result,Site,not_norm))%>%
    mutate(Letters=paste0(Letters,not_norm))->
    Tk
  letters_y <- ind_uRootMass%>%filter(SamplingEvent==i)%>%.$uRootMass%>%max()+0.05
  
  
  ind_uRootMass%>%
    filter(SamplingEvent==i)%>%
    ggplot(aes(Site,uRootMass))+
    geom_linerange(aes(Site,ymin=minusbar_sd,ymax=plusbar_sd),
                 data=ind_uRootMass_stats,linewidth=1,colour="black")+
    # geom_linerange(aes(Site,ymin=minusbar_se,ymax=plusbar_se),
    #                data=ind_uRootMass_stats,linewidth=2,colour="black")+
    geom_point(aes(Site,uRootMass),data=ind_uRootMass_stats,size=3)+
    geom_text(data = Tk, aes(x = Site, y = letters_y, label = Letters),fontface = "italic")+
    geom_point(shape=4,colour="grey",size=3)+
    theme_bw()+
    ylab("")+xlab("")+
    theme(axis.title.x = element_blank(),axis.text.x = element_blank(),axis.ticks.x = element_blank())
}

plotter(4)->A
plotter(5)->B
plotter(6)->C
plotter(7)->D
plotter(8)->E

(A+ggtitle("35 days"))/
  (B+ggtitle("42 days"))/
  (C+ggtitle("48 days")+
     ylab("Mean root mass (g)"))/
  (D+ggtitle("56 days"))/(E+
           theme(axis.text.x = element_text(),axis.ticks.x = element_line(),
                 axis.title.x=element_text())+
           #scale_y_continuous(limits = c(NA, 0.6))+
           #xlab("Site")+
             ggtitle("103 days")+xlab("Site"))

ggsave("FinalFiguresV2/Fig_4.jpg",width=5,height=8)