setwd()
library(agricolae)
library(ggplot2)
library (tidyverse)
library(dplyr)

mydata =  read.table("ROOT.txt", header=TRUE, sep="\t")
head(mydata)

## check str for confirmation
str(mydata)

## model anova
mydata$year<-as.factor(mydata$year)
mydata$rep<-as.factor(mydata$rep)
mydata$A<-as.factor(mydata$A)
mydata$B<-as.factor(mydata$B)
mydata$C<-as.factor(mydata$C)
mod = aov(X1 ~year*A*B*C+Error(year/rep/A+year/rep),data = mydata)
summary(mod)
###Coefficient of Variance for MainPlot
cv.mainplot=sqrt(1019)*100/mean(mydata$X1)
cv.mainplot

#Coefficient of Variance for SubPlot
cv.subplot=sqrt(139)*100/mean(mydata$X1)
cv.subplot

attach(mydata)
pairwise.t.test(X1,year,p.adj="bonferroni")
pairwise.t.test(X1,A,p.adj="bonferroni")
pairwise.t.test(X1,year:A,p.adj="bonferroni")
pairwise.t.test(X1,B,p.adj="bonferroni")
pairwise.t.test(X1,C,p.adj="bonferroni")
pairwise.t.test(X1,year:B,p.adj="bonferroni")
pairwise.t.test(X1,A:B,p.adj="bonferroni")
pairwise.t.test(X1,year:C,p.adj="bonferroni")
pairwise.t.test(X1,A:C,p.adj="bonferroni")
pairwise.t.test(X1,B:C,p.adj="bonferroni")
pairwise.t.test(X1,year:A:B,p.adj="bonferroni")
pairwise.t.test(X1,year:A:C,p.adj="bonferroni")
pairwise.t.test(X1,year:B:C,p.adj="bonferroni")
pairwise.t.test(X1,A:B:C,p.adj="bonferroni")
pairwise.t.test(X1,year:A:B:C,p.adj="bonferroni")

Duncan_Y<-duncan.test(X1, year, 7, 12.3, 0.05)
Duncan_Y
Duncan_A<-duncan.test(X1, A, 7, 12.3, 0.05)
Duncan_A
Duncan_YA<-duncan.test(X1, year:A, 7, 12.3, 0.05)
Duncan_YA
Duncan_B<-duncan.test(X1, B, 276, 139, 0.05)
Duncan_B
Duncan_C<-duncan.test(X1, C, 276, 139, 0.05)
Duncan_C
Duncan_YB<-duncan.test(X1, year:B,276, 139, 0.05)
Duncan_YB
Duncan_AB<-duncan.test(X1, A:B, 276, 139, 0.05)
Duncan_AB
Duncan_YC<-duncan.test(X1, year:C,276, 139, 0.05 )
Duncan_YC
Duncan_AC<-duncan.test(X1, A:C, 276, 139, 0.05)
Duncan_AC
Duncan_BC<-duncan.test(X1, B:C, 276, 139, 0.05)
Duncan_BC
Duncan_YAB<-duncan.test(X1, year:A:B, 276, 139, 0.05)
Duncan_YAB
Duncan_YAC<-duncan.test(X1, year:A:C, 276, 139, 0.05)
Duncan_YAC
Duncan_YBC<-duncan.test(X1, year:B:C, 276, 139, 0.05)
Duncan_YBC
Duncan_ABC<-duncan.test(X1, A:B:C, 276, 139, 0.05)
Duncan_ABC
Duncan_YABC<-duncan.test(X1, year:A:B:C, 276, 139, 0.05)
Duncan_YABC

##duncan year###
mod1 = aov(X1 ~year*A*B*C,data = mydata)
out_Y= Duncan_Y$groups %>% group_by(row.names(Duncan_Y$groups))%>% arrange(row.names(Duncan_Y$groups))
out_Y
MeanSE_Y=mydata%>%group_by(year)%>%summarise(avg_Y = mean(X1),
                                             sde = sd(X1)/sqrt(length(X1)))

print(MeanSE_Y)
attach(MeanSE_Y)


ggplot(MeanSE_Y, aes(x=year, y=avg_Y))+
  geom_bar(stat = "identity", color = "black", position = position_dodge(width=0.5),width = 0.4)+theme_bw()+
  geom_text(aes(x=year,y=avg_Y+sde,label=as.matrix(out_Y$groups)),position = position_dodge(width = 0.5),angle = 90,vjust=-(0.1))+
  geom_errorbar( aes(ymax = avg_Y + sde,ymin = avg_Y - sde), position = position_dodge(width=0.5),width = 0.1)+
  theme(axis.text.x = element_text(angle = 90, vjust =-0.1)) +
  ggtitle("duncan_Y")+labs(x="factor Y")


###duncan A ###

mod1 = aov(X1 ~year*A*B*C,data = mydata)

out_A= Duncan_A$groups %>% group_by(row.names(Duncan_A$groups))%>% arrange(row.names(Duncan_A$groups))
MeanSE_A=mydata%>%group_by(A)%>%summarise(avg_A = mean(X1),
                                          sde = sd(X1)/sqrt(length(X1)))

print(MeanSE_A)
attach(MeanSE_A)


ggplot(MeanSE_A, aes(x=A, y=avg_A))+
  geom_bar(stat = "identity", color = "black", position = position_dodge(width=0.5),width = 0.4)+theme_bw()+
  geom_text(aes(x=A,y=avg_A+sde,label=as.matrix(out_A$groups)),position = position_dodge(width = 0.5),vjust=-(0.3))+
  geom_errorbar( aes(ymax = avg_A + sde,ymin = avg_A - sde), position = position_dodge(width=0.5),width = 0.1)+
  theme(axis.text.x = element_text(angle = 90, vjust = 1)) +
  ggtitle("duncan_A")+labs(x="factor A")


####duncan YEAR:A ###
mod1 = aov(X1 ~year*A*B*C,data = mydata)

out_YA= Duncan_YA$groups %>% group_by(row.names(Duncan_YA$groups))%>% arrange(row.names(Duncan_YA$groups))
MeanSE_YA=mydata%>%group_by(year,A)%>%summarise(avg_YA = mean(X1),
                                                sde = sd(X1)/sqrt(length(X1)))
print(MeanSE_YA)
attach(MeanSE_YA)

ggplot(MeanSE_YA, aes(x =A,  y = avg_YA, fill = factor(year)))+
  geom_bar(stat = "identity", color = "black", position = position_dodge(width=0.5),width = 0.4)+theme_bw()+
  geom_text(aes(x=A,y=avg_YA+sde,label=as.matrix(out_YA$groups)),position = position_dodge(width = 0.5),vjust=-(0.3))+
  geom_errorbar( aes(ymax = avg_YA + sde,ymin = avg_YA - sde), position = position_dodge(width=0.5),width = 0.1)+
  theme(axis.text.x = element_text(angle = 90, vjust = 1)) +
  ggtitle("duncan_YA")+labs(x="factor year:A")

###DUNCAN B###


mod1 = aov(X1 ~year*A*B*C,data = mydata)
out_B= Duncan_B$groups %>% group_by(row.names(Duncan_B$groups))%>% arrange(row.names(Duncan_B$groups))
MeanSE_B=mydata%>%group_by(B)%>%summarise(avg_B = mean(X1),
                                          sde = sd(X1)/sqrt(length(X1)))

print(MeanSE_B)
attach(MeanSE_B)
ggplot(MeanSE_B, aes(x=B, y=avg_B))+
  geom_bar(stat = "identity", color = "black", position = position_dodge(width=0.5),width = 0.5)+theme_bw()+
  geom_text(aes(x=B,y=avg_B+sde,label=as.matrix(out_B$groups)),position = position_dodge(width = 0.5),vjust=-(0.3))+
  geom_errorbar( aes(ymax = avg_B + sde,ymin = avg_B - sde), position = position_dodge(width=0.5),width = 0.1)+
  theme(axis.text.x = element_text(angle = 90, vjust = 1)) +
  ggtitle("duncan_B")+labs(x="factor B")

###DUNCAN C###

mod1 = aov(X1 ~year*A*B*C,data = mydata)


out_C= Duncan_C$groups %>% group_by(row.names(Duncan_C$groups))%>% arrange(row.names(Duncan_C$groups))
print(out_C)
MeanSE_C= mydata%>%group_by(C)%>% summarise(avg_C = mean(X1),
                                            sde = sd(X1)/sqrt(length(X1)))


print(MeanSE_C)
attach(MeanSE_C)
for (i in which(nchar(out_C$groups)>2)) {
  out_C$groups[i]=paste0(str_sub(out_C$groups[i],1,1),"-",
                         str_sub(out_C$groups[i],nchar(out_C$groups[i])))
}
MeanSE_C=MeanSE_C[out_C$`row.names(Duncan_C$groups)`,]

ggplot(MeanSE_C, aes(x=C, y=avg_C,))+
  geom_bar(stat = "identity", color = "black", position = position_dodge(width=0.9),width = 0.7)+theme_bw()+
  geom_text(aes(x=C,y=avg_C+sde,label=as.matrix(out_C$groups)),position = position_dodge(width = 0.9),vjust=-(0.3),angle = 0,size =4, vjust = -(0.3))+
  geom_errorbar( aes(ymax = avg_C + sde,ymin = avg_C - sde), position = position_dodge(width=0.9),width = 0.1)+
  theme(axis.text.x = element_text(angle = 0, vjust = 1)) +
  ggtitle("duncan_C")+labs(x="factor C")+
  scale_y_continuous(name="avg C", limits=c(0, 200))  

#####YB###
mod1 = aov(X1 ~year*A*B*C,data = mydata)
out_YB= Duncan_YB$groups %>% group_by(row.names(Duncan_YB$groups))%>% arrange(row.names(Duncan_YB$groups))
MeanSE_YB=mydata%>%group_by(year,B)%>%summarise(avg_YB = mean(X1),
                                                sde = sd(X1)/sqrt(length(X1)))
print(MeanSE_YB)
attach(MeanSE_YB)
ggplot(MeanSE_YB, aes(x =B,  y = avg_YB, fill = factor(year)))+
  geom_bar(stat = "identity", color = "black", position = position_dodge(width=0.5),width = 0.4)+theme_bw()+
  geom_text(aes(x=B,y=avg_YB+sde,label=as.matrix(out_YB$groups)),position = position_dodge(width = 0.5),vjust=-(0.3))+
  geom_errorbar( aes(ymax = avg_YB + sde,ymin = avg_YB - sde), position = position_dodge(width=0.5),width = 0.1)+
  theme(axis.text.x = element_text(angle = 90, vjust = 1)) +
  ggtitle("duncan_YB")+labs(x="factor year:B")

####AB###
mod1 = aov(X1 ~year*A*B*C,data = mydata)

out_AB= Duncan_AB$groups %>% group_by(row.names(Duncan_AB$groups))%>% arrange(row.names(Duncan_AB$groups))
MeanSE_AB=mydata%>%group_by(A:B)%>%summarise(avg_AB = mean(X1),
                                             sde = sd(X1)/sqrt(length(X1)))
print(MeanSE_AB)
attach(MeanSE_AB)
ggplot(MeanSE_AB, aes(x =A,  y = avg_AB, fill = factor(B)))+
  geom_bar(stat = "identity", color = "black", position = position_dodge(width=0.5),width = 0.4)+theme_bw()+
  geom_text(aes(x=A, y=avg_AB+sde,label=as.matrix(out_AB$groups)),position = position_dodge(width = 0.5),vjust=-(0.3))+
  geom_errorbar( aes(ymax = avg_AB + sde,ymin = avg_AB - sde), position = position_dodge(width=0.5),width = 0.1)+
  theme(axis.text.x = element_text(angle = 90, vjust = 1)) +
  ggtitle("duncan_AB")+labs(x="factor A:B")


####Duncan_YABC##
Duncan_YABC<-duncan.test(X1, year:A:B:C, 276, 139, 0.05,p.adj="bonferroni",group = TRUE, console = TRUE)
Duncan_YABC

out_YABC= Duncan_YABC$groups %>% group_by(row.names(Duncan_YABC$groups))%>% arrange(row.names(Duncan_YABC$groups))
MeanSE_YABC=mydata%>%group_by(year,A,B,C)%>%summarise(avg_YABC = mean(X1),
                                                      sde = sd(X1)/sqrt(length(X1)))
print(MeanSE_YABC)
attach(MeanSE_YABC)



for (i in which(nchar(out_YABC$groups)>2)) {
  out_YABC$groups[i]=paste0(str_sub(out_YABC$groups[i],1,1),"-",
                            str_sub(out_YABC$groups[i],nchar(out_YABC$groups[i])))
}
MeanSE_YABC$X=paste0(MeanSE_YABC$year,":",MeanSE_YABC$A,":",MeanSE_YABC$B,":",MeanSE_YABC$C)
rownames(MeanSE_YABC)=MeanSE_YABC$X
MeanSE_YABC=MeanSE_YABC[out_YABC$`row.names(Duncan_YABC$groups)`,]
ggplot(MeanSE_YABC[MeanSE_YABC$A==1,], aes(x =year:A:B:C,  y= avg_YABC))+
  geom_bar(stat = "identity", color = "black", position = position_dodge(width=0.5),width = 0.4)+theme_bw()+
  geom_text(aes(x=year:A:B:C,y=avg_YABC+sde,label=as.matrix(out_YABC$groups[which(MeanSE_YABC$A==1)])),
            position = position_dodge(width = 0.5),vjust=-(0.3),angle = 90, vjust = 1,size =2.5)+
  geom_errorbar( aes(ymax = avg_YABC + sde,ymin = avg_YABC - sde), position = position_dodge(width=0.5),width = 0.1)+
  theme(axis.text.x = element_text(angle = 90, vjust = 1)) +
  ggtitle("duncan_YABC for Factor 1")+labs(x="factorY:A:B:C")

ggplot(MeanSE_YABC[MeanSE_YABC$A==2,], aes(x =year:A:B:C,  y= avg_YABC))+
  geom_bar(stat = "identity", color = "black", position = position_dodge(width=0.5),width = 0.4)+theme_bw()+
  geom_text(aes(x=year:A:B:C,y=avg_YABC+sde,label=as.matrix(out_YABC$groups[which(MeanSE_YABC$A==2)])),
            position = position_dodge(width = 0.5),vjust=-(0.3),angle = 90, vjust = 1,size =2)+
  geom_errorbar( aes(ymax = avg_YABC + sde,ymin = avg_YABC - sde), position = position_dodge(width=0.5),width = 0.1)+
  theme(axis.text.x = element_text(angle = 90, vjust = 1)) +
  ggtitle("duncan_YABC for Factor 2")+labs(x="factorY:A:B:C")

ggplot(MeanSE_YABC[MeanSE_YABC$A==3,], aes(x =year:A:B:C,  y= avg_YABC))+
  geom_bar(stat = "identity", color = "black", position = position_dodge(width=0.5),width = 0.4)+theme_bw()+
  geom_text(aes(x=year:A:B:C,y=avg_YABC+sde,label=as.matrix(out_YABC$groups[which(MeanSE_YABC$A==3)])),
            position = position_dodge(width = 0.5),vjust=-(0.3),angle = 90, vjust = 1,size =2)+
  geom_errorbar( aes(ymax = avg_YABC + sde,ymin = avg_YABC - sde), position = position_dodge(width=0.5),width = 0.1)+
  theme(axis.text.x = element_text(angle = 90, vjust = 1)) +
  ggtitle("duncan_YABC for Factor 3")+labs(x="factorY:A:B:C")