setwd("E:/SENSOMETRIE/classiqueF")  #Répertoire contenant le fichier expert.csv
library(SensoMineR)

######### ################################################################################
######################### CHARGEMENT DES DONNEES     #####################################
#########################   ANALYSE DU FICHIER       #####################################
##########################################################################################
#Importation des données avec détail des notes attribuées par les 12 juges par descripteur et par session
#Dans le jeu de données, il y a 16 produits évalués par 12 juges lors de 2 sessions.
#Il y a donc des répétitions par juge.

classique <- read.table("Experts.csv", header=TRUE, sep=";", na.strings="", dec=",")
names(classique) #nom des variables
#Supprime les variables non utilisées (commentaires) et renomme les variables (libellés anglais)
classique<-classique[,-23]
classique<-classique[,-21]
classique<-classique[,-18]
classique<-classique[,-8]
classique<-classique[,-6]
names(classique)[1]="Session"
names(classique)[2]="Panelist"
names(classique)[3]="Rank"
names(classique)[4]="Juice"
names(classique)[5]="color.intensity"
names(classique)[6]="odor.intensity"
names(classique)[7]="strongness"
names(classique)[8]="sweet"
names(classique)[9]="acidity"
names(classique)[10]="bitterness"
names(classique)[11]="odor.orange"
names(classique)[12]="odor.banana"
names(classique)[13]="odor.mango"
names(classique)[14]="odor.lemon"
names(classique)[15]="persistence"
names(classique)[16]="pulp"
names(classique)[17]="thickness"
names(classique)[18]="hedonic"
dim(classique)   #dimension du jeu de données
for(i in 1:3) classique[,i]<-as.factor(classique[,i]) #transforme en facteur

#Pour voir si données équilibrées
table(classique$Session,classique$Panelist) #xhaque juge a goûté les 16 produits
table(classique$Session,classique$Juice)  #chaque produit a été testé aux 2 séances
table(classique$Panelist,classique$Juice) #chaque juge a goûté 2 fois chaque produit

numSummary(classique, statistics=c("mean", "sd")) #statistiques descriptives
#On voit ci-dessous qu'il y a des données manquantes ==> pb pour les analyses factorielles

#odor.intensity  5.243316 2.0612981 374 10
#strongness      6.209115 2.0038196 373 11
#sweet           6.093333 2.0460823 375  9
#acidity         5.657068 2.6924976 382  2
#bitterness      1.957560 1.4634928 377  7
#odor.orange     5.054974 2.4792290 382  2
#odor.banana     3.937337 3.0853887 383  1
#odor.mango      4.109375 2.9073627 384  0
#odor.lemon      4.373684 2.9672517 380  4
#persistence     6.171271 1.8972069 362 22
#pulp            2.194149 1.7133415 376  8
#thickness       5.301837 2.0965440 381  3
#hedonic         4.949735 2.1483724 378  6


#Boîte à moustaches pour chacun des descripteurs
boxprod(classique[,c("Session", "Panelist", "Rank", "Juice", "color.intensity", "odor.intensity", "strongness", "sweet", "acidity", "bitterness", 
  "odor.orange", "odor.banana", "odor.mango", "odor.lemon", "persistence", "pulp", "thickness", "hedonic")],firstvar=5,col.p=4,numr=2,numc=1)

######### ################################################################################
############################  ANALYSE MULTIVARIEE    #####################################
##########################################################################################

##################
#1ère METHODE: ACP via panellipse Menu SensoMineR/CharacterizationProducts/ultidimensional Sensory Profile = panellipse
##################
#Les données ne sont pas moyennées par juge dans le fichier mais panellipse les agrège.
#level.search.desc=0.2 indique le seuil pour conserver les variables dans les anova Juice~descripteur
#On ne peut pas avoir de variable illustrative avec panellipse.
#On ne peut pas utiliser le fichier de sorties tel quel pour lancer une classification (HCPC).
#Il faut récupérer les coordonnées factorielles sur les 2 premiers axes dans r_panellipse$coord$moyen[1:16,]
#les lignes suivantes correspondent aux projections des couples produits-juges.

r_panellipse=panellipse(classique[,c("Session", "Panelist", "Rank", "Juice", "color.intensity", "odor.intensity", 
  "strongness", "sweet", "acidity", "bitterness", "odor.orange", "odor.banana", "odor.mango", "odor.lemon", 
  "persistence", "pulp", "thickness")],col.p=4,col.j=2,firstvar=5,alpha=0.05,coord = c(1,2), nbsimul =500,nbchoix =NULL,
  level.search.desc=0.2,scale.unit=1,variability.variable =TRUE,centerbypanelist =TRUE,scalebypanelist=FALSE,
  name.panelist=FALSE)

r_panellipse$eig   #valeurs propres
dev.new()
barplot(r_panellipse$eig[,1], main="Eigenvalues", xlab="Dimension", ylab="Eigenvalues",
  names.arg=1:nrow(r_panellipse$eig))
r_panellipse$coordinates
r_panellipse$correl
r_panellipse$hotelling
coltable(r_panellipse$hotelling, main.title = "P-values for the Hotelling T2 tests")


##################
#2ème METHODE 
##################
#Une ACP classique (menu FactoMineR/PCA) sur le fichier ne peut convenir ici car données enregistrées sur 2 séances 
#Il faut moyenner les données par juge et faire une ACP sur un tableau 16 produits*13 descripteurs.
#On peut garder des variables illustratives.
#Avec le calcul des données moyennées par juge avec Menu SensoMineR/Tools Functions/average by product and by descriptor 
#puis PCA, on n'a pas les mêmes inerties qu'avec panellipse (problème d'estimation des données manquantes ?).

#Pour s'en affranchir : le tableau moyennes récupère les moyennes par produit (col.p=4) tous juges (col.j=2) et sessions confondus  soit 16 lignes
#ainsi que les moyennes par produit et par juge toutes sessions confondues soit 16*12 = 192 lignes et 16 variables

moyennes=scalebypanelist(classique,firstvar=5,center=FALSE,scale=FALSE,col.p=4,col.j=2)
dim(moyennes) #208 14
names(moyennes)
#Panelist Juice color.intensity odor.intensity strongness sweet acidity bitterness odor.orange odor.banana odor.mango
#odor.lemon persistence pulp thickness hedonic

average=moyennes[1:16,] #seul les moyennes par produit nous intéressent
#rownames a été modifiée par la fonction scalebypanelist ; on lui réattribue le nom du jus de fruit (variable 2)
row.names(average)=average[,2] 

#On lance une ACP sur le tableau des moyennes, avec "note hédonique" (variable 14) en supplémentaire
res_pcamoy=PCA(average[1:16,-(1:2)],quanti.sup=14:14)
#Représentation des jus de fruits sur le plan (1,2)
plot.PCA(res_pcamoy, axes=c(1, 2), choix="ind", habillage="none", col.ind="black", col.ind.sup="blue", col.quali="magenta", 
  label=c("ind", "ind.sup", "quali"), title="")

#Représentation des 13 descripteurs + note hédonique (en supplémentaire) sur le cercle des corrélations.
plot.PCA(res_pcamoy, axes=c(1, 2), choix="var", col.var="black", col.quanti.sup="blue", label=c("var", "quanti.sup"), 
  lim.cos2.var=0, title="")

#sauvegarde des résultats
write.infile(res_pcamoy$eig, file ="PCA_moy.csv",append=FALSE)
write.infile(res_pcamoy$var, file ="PCA_moy.csv",append=TRUE)
write.infile(res_pcamoy$ind, file ="PCA_moy.csv",append=TRUE)
write.infile(res_pcamoy$quanti.sup, file ="PCA_moy.csv",append=TRUE)
write.infile(dimdesc(res_pcamoy, axes=c(1, 2)), file ="PCA_moy.csv",append=TRUE)
  
#Sélection des variables les plus représentatives
var_mieux<-dimdesc(res_pcamoy, axes=c(1, 2))
write.infile(var_representative, file ="var_mieux_acp_moy.csv",append=FALSE)

#######################################################
#######################             CAH                
#######################################################
#On peut ensuite lancer une CAH automatique sur les résultats de l'ACP.

res.HCPC<-HCPC(res_pcamoy,nb.clust=-1)
res.HCPC      

#Répartition des classes
res.HCPC$desc.ind


################################################################################
##############              Analyses univariées                    #############
################################################################################

#Donne les produits en fonction des descripteurs les plus importants 
#Modèle note_descripteur~Juice+Panelist+Juice:Panelist
resdecat=decat(classique[,c("Session", "Panelist", "Rank", "Juice", "color.intensity", "odor.intensity", "strongness", 
  "sweet", "acidity", "bitterness", "odor.orange", "odor.banana", "odor.mango", "odor.lemon", "persistence", "pulp", 
  "thickness","hedonic")],firstvar=5,formul=~Juice+Panelist+Juice:Panelist,proba=0.05,graph=TRUE,col.lower="mistyrose",
  col.upper="lightblue")
resdecat

#
respanelperf=panelperf(classique[,c("Session", "Panelist", "Rank", "Juice", "color.intensity", "odor.intensity", 
  "strongness", "sweet", "acidity", "bitterness", "odor.orange", "odor.banana", "odor.mango", "odor.lemon", 
  "persistence", "pulp", "thickness","hedonic")],firstvar=5,formul=~"Juice+Panelist+Juice:Panelist",random=1)
respanelperf$p.value
coltable(magicsort(respanelperf$p.value, sort.mat = respanelperf$p.value[,1], bycol = FALSE,method = "median"),
  main.title = "Panel performance (sorted by product P-value)")
respanelperf$complete
respanelperf$variability
respanelperf$res