#Developed by Elie Maza (AgroToulouse) for Dolet's Oeno One paper (2026)
#********************************************************
library(agricolae)

X=read.table("CultivarScore.csv",sep=";",header=TRUE,stringsAsFactors=TRUE)

X$NumPanel=as.factor(X$NumPanel)

summary(X)
boxplot(X$Liking~paste(X$Colour,X$Resistance,X$Phase))

anova.hedo.4=aov(Liking~Colour*Resistance*Phase*Profiles*Climat,data=X)
sa=summary(anova.hedo.4)
sa

tCouleur=HSD.test(anova.hedo.4,"Colour")
tCouleur

tResistance=HSD.test(anova.hedo.4,"Resistance")
tResistance

tPhase=HSD.test(anova.hedo.4,"Phase")
tPhase

tProfiles=HSD.test(anova.hedo.4,"Profiles")
tProfiles

tClimat=HSD.test(anova.hedo.4,"Climat")
tClimat

pdf(paste("fig_tuckey_4_facteurs_plus.pdf",sep=""),width=12,height=8)
par(mar=c(5.1,5.1,2.1,1.1))
layout(matrix(c(1,2,2,3,4,5),byrow=TRUE,nrow=2,ncol=3))
plot(tCouleur,variation="SE",ylab="Liking",cex=1.5,cex.names=1.5,cex.axis=1.5,cex.lab=1.5,ylim=c(4.5,5.6),main="",xlab="Color")
legend("topright",legend=paste("p-value=",round(sa[[1]][1,5],5)),cex=1.5)
plot(tProfiles,variation="SE",ylab="Liking",cex=1.5,cex.names=1.5,cex.axis=1.5,cex.lab=1.5,ylim=c(4.5,5.6),main="",xlab="Profiles")
legend("topright",legend=paste("p-value=",round(sa[[1]][4,5],4)),cex=1.5)
plot(tResistance,variation="SE",ylab="Liking",cex=1.5,cex.names=1.5,cex.axis=1.5,cex.lab=1.5,ylim=c(4.5,5.6),main="",xlab="Resistance")
legend("topright",legend=paste("p-value=",round(sa[[1]][2,5],3)),cex=1.5)
plot(tClimat,variation="SE",ylab="Liking",cex=1.5,cex.names=1.5,cex.axis=1.5,cex.lab=1.5,ylim=c(4.5,5.6),main="",xlab="Climat")
legend("topright",legend=paste("p-value=",round(sa[[1]][5,5],2)),cex=1.5)
plot(tPhase,variation="SE",ylab="Liking",cex=1.5,cex.names=1.5,cex.axis=1.5,cex.lab=1.5,ylim=c(4.5,5.6),main="",xlab="Phase")
legend("topright",legend=paste("p-value=",round(sa[[1]][3,5],1)),cex=1.5)
dev.off()

#---------
# The End
#---------