# Coded by Alastair Potts
# For review

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

groundhog.library("dplyr",GroundhogDay)
groundhog.library("ggplot2",GroundhogDay)
groundhog.library("lubridate",GroundhogDay)
groundhog.library("tibble",GroundhogDay)
groundhog.library("tidyr",GroundhogDay)
groundhog.library("rstatix",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

dat_all%>%
  mutate(Site=factor(Site,labels=paste("Site:",1:10)))->
  dat_all
  
read.csv("data_manuscript/harvest_days.csv")%>%
  as_tibble()->
  harvestdates


dat_all%>%
  group_by(SamplingEvent,Site,Individual)%>%
  summarise(n=n(),
            Y=sum(RootedYN=="Y",na.rm=T),
            N=sum(RootedYN=="N",na.rm=T))%>%
  mutate(RootedPerc= round(Y/(Y+N)*100,1))%>%
  left_join(harvestdates,by="SamplingEvent")%>%
  mutate(DaysSincePlanting=factor(DaysSincePlanting))%>%
  mutate(Ind2=LETTERS[Individual])%>%
  filter(!is.na(RootedPerc))->dat_perc

### Testing for normality



## Transforming percentage data using arcsin
dat_perc%>%
  group_by(SamplingEvent,Site)%>%
  shapiro_test(RootedPerc)

dat_perc%>%
  filter(SamplingEvent==2)%>%
  mutate(arcsin_RootedPerc=asin(RootedPerc/100))%>%
  group_by(SamplingEvent,Site)%>%
  shapiro_test(arcsin_RootedPerc)
# Result: not all sites have normally distributed data within each harvest event

# .:. switching to the non-parametric kruskal test
dat_perc%>%
  group_by(SamplingEvent)%>%
  kruskal_test(RootedPerc~Site)%>%
  get_test_label(type="text",detailed=TRUE)

dat_perc%>%
  group_by(SamplingEvent)%>%
  kruskal_test(RootedPerc~Site)%>%
  mutate(text=paste0("H=",round(statistic,1),", df=",df,", p=",p))%>%
  select(text)


# All harvest events contain significant differences amongst sites for the 
# percentage of rooted cuttings (parent plant).

# conducting a post-hoc dunn_test.
# https://www.rdocumentation.org/packages/stats/versions/3.6.2/topics/p.adjust
dat_perc%>%
  group_by(SamplingEvent)%>%
  dunn_test(RootedPerc~Site,p.adjust.method = "hochberg")->
  post_hoc

cld<-list()
for (i in 2:8) {
  post_hoc%>%
    filter(SamplingEvent==i)->
    post_hoc_f
  post_hoc_f$p.adj->ph_v
  names(ph_v) <- paste(post_hoc_f$group1,post_hoc_f$group2,sep="-")  
  cld[[i]]<- multcompLetters(ph_v)$Letters
}
data.frame(SamplingEvent=sort(rep(2:8,10)),letters=unlist(cld),Site=paste0("Site: ",1:10))->
  post_hoc_letters

dat_perc%>%
  group_by(SamplingEvent,Site)%>%
  summarise(uP=mean(RootedPerc),sdP=sd(RootedPerc))%>%
  left_join(post_hoc_letters,by=c("SamplingEvent","Site"))%>%
  mutate(text=paste0(round(uP,0),"\U00B1",round(sdP,0),letters))%>%
  select(SamplingEvent,Site,text)%>%
  pivot_wider(names_from = Site,values_from = text)%>%
  write.csv("results/Table_percRooting.csv")

(dat_all%>%
    group_by(SamplingEvent,Site,Individual)%>%
    summarise(n=n(),
              Y=sum(RootedYN=="Y",na.rm=T),
              N=sum(RootedYN=="N",na.rm=T))%>%
    mutate(RootedPerc= round(Y/(Y+N)*100,1))%>%
    left_join(harvestdates,by="SamplingEvent")%>%
    mutate(DaysSincePlanting=factor(DaysSincePlanting))%>%
    mutate(Ind2=LETTERS[Individual])%>%
    filter(!is.na(RootedPerc))->tata1)%>%
  ungroup()%>%
  group_by(SamplingEvent,Site)%>%
  summarise(uRootingPerc=mean(RootedPerc),
            sdRootingPerc=sd(RootedPerc),
            plusbar=uRootingPerc+sdRootingPerc,
            minusbar=uRootingPerc-sdRootingPerc)%>%
  left_join(harvestdates,by="SamplingEvent")%>%
  mutate(DaysSincePlanting=factor(DaysSincePlanting))%>%
  ggplot(aes(DaysSincePlanting,uRootingPerc,group=Site))+
  geom_point(data=tata1,aes(DaysSincePlanting,RootedPerc),colour="grey",alpha=0.5)+
  geom_line(data=tata1,
            aes(DaysSincePlanting,RootedPerc,group=Individual,linetype=as.factor(Individual)),
            colour="grey",alpha=0.5)+
  geom_linerange(aes(DaysSincePlanting,ymin=minusbar,ymax=plusbar),linewidth=0.2,colour="black")+
  geom_line(linewidth=1.2)+
  geom_point(size=3)+
  #geom_text_repel(aes(label=Pop))+
  facet_wrap(~Site,nrow=2)+
  ylab("Rooted Percentage (n=6 per plant)")+
  xlab("Days since start of experiment")+
  theme_bw()+
  coord_cartesian(ylim = c(0, 100))+
  theme(legend.position="none") 

ggsave("FinalFiguresV2/Fig_3.jpg",width=10,height=5)