library(SensoMineR)
# 6.1
library(ggplot2)
liking <- read.csv("data/sensory/perfumes_liking.csv")
length(unique(liking$consumer))
liking$consumer <- as.factor(liking$consumer)
liking$product <- as.factor(liking$product)
consumer.names <- levels(liking$consumer)
product.names <- levels(liking$product)

iLike=NULL
for (i in consumer.names) {
  tmp=subset(liking, consumer==i)[,"liking"]
  iLike=rbind(iLike,tmp)
}
rownames(iLike)=consumer.names
colnames(iLike)=product.names

#write.csv(iLike,file="iLike.csv")


rawPCA <- PCA(iLike,graph=FALSE)
dev.new();plot(rawPCA,choix="ind",title="Liking is not centered")#Fig 6.1-4
iLike[c("11174","13311"),]
iLike[c("4640","1755","8840"),]

dev.new()
plot(rawPCA,choix="var",title="Liking is not centered") #Fig 6.1-5

iLike2 <- cbind(iLike,apply(iLike,1,mean))
colnames(iLike2)[13] <- "Average"
rawPCA2 <- PCA(iLike2,quanti.sup=13,graph=FALSE)
dev.new();plot(rawPCA2,choix="var", col.quanti.sup="red") #Fig 6.1-6

iLike_centered <- t(scale(t(iLike),scale=FALSE))
head(iLike_centered)
 
center.pca <- PCA(iLike_centered,graph=FALSE)
dev.new()
plot(center.pca,choix="var",title="Liking is centered")#Fig 6.1-7

dev.new()
plot(rawPCA,choix="var",axes=c(2,3),title="Liking is NOT centered") #Fig 6.1-8



options(contrasts=c("contr.sum", "contr.sum"))
output1 <- lm(liking~product+consumer,data=liking)
anova(output1)
coef(summary(output1))[1:12,]

model <- aov(liking~product+consumer,data=liking)
#library(agricolae)
LSD.compare <- agricolae::LSD.test(model,"product",p.adj="none",group=FALSE,
                    main="Results of the Fisher LSD")
names(LSD.compare)
LSD.compare$statistics
LSD.compare$comparison

LSD.group <- agricolae::LSD.test(model,"product",p.adj="none",group=TRUE,
                                   main="Results of the Fisher LSD")
LSD.group$groups #Fig 6.1-11

data.mdpref <- t(iLike)
iLike.mdpref <- PCA(data.mdpref,ncp=4,graph=FALSE)

dev.new();plot(iLike.mdpref,choix="ind") #Fig 6.1-12
dev.new();plot(iLike.mdpref,choix="var") #Fig 6.1-13

iLike.HCPC<- HCPC(iLike.mdpref, graph=FALSE)

dev.new();plot(iLike.HCPC,choice="tree")#Fig 6.1-14
dev.new();plot(iLike.HCPC,choice="3D.map")#Fig 6.1-15



# 6.2

liking <- read.csv("data/sensory/perfumes_liking.csv")
liking$consumer <- as.factor(liking$consumer)
liking$product <- as.factor(liking$product)
consumer.names <- levels(liking$consumer)
product.names <- levels(liking$product)

hedonic.c=NULL
for (i in consumer.names) {
  tmp=subset(liking, consumer==i)[,"liking"]
  hedonic.c=rbind(hedonic.c,tmp)
}
rownames(hedonic.c)=paste0("C",consumer.names)
colnames(hedonic.c)=product.names
Hedonic=t(hedonic.c)


experts<- read.csv("data/sensory/perfumes_qda_experts.csv")
experts$Session <- as.factor(experts$Session)
experts$Rank <- as.factor(experts$Rank)

#library(SensoMineR)
sensory <- SensoMineR::averagetable(experts,formul="~Product+Panelist",firstvar=5)
sensory=sensory[rownames(Hedonic),]

dat <- cbind(Hedonic,sensory)
head(round(dat,2))

mdpref <- PCA(dat,
              quanti.sup=seq(to=ncol(dat),length=12),
              graph=FALSE)
dev.new();plot(mdpref,choix="ind")
dev.new();plot(mdpref,choix="var",label="quanti.sup")

res.dimdesc <- dimdesc(mdpref)
CORR1=res.dimdesc$'Dim.1'$quanti
select.supp <- which(rownames(CORR1) %in% colnames(sensory))
res.dimdesc$'Dim.1'$quanti[select.supp,]


PCA.hedonic  <- PCA(Hedonic,graph=FALSE)
dev.new();plot(PCA.hedonic, choix="ind")
dev.new();plot(PCA.hedonic,choix="var", col.var="blue")

PCA.sensory <- PCA(sensory,graph=FALSE)
dev.new();plot(PCA.sensory,choix="ind")
dev.new();plot(PCA.sensory,choix="var", col.var="blue")




Hedonic.avg <- as.matrix(apply(Hedonic,1,mean))
rownames(Hedonic.avg) <- rownames(Hedonic)
data4lm <- cbind(Hedonic.avg[rownames(sensory),],sensory)
colnames(data4lm)[1] <- "Liking"

coef(summary(lm(Liking~Spicy,data=data4lm)))

coef(summary(lm(Liking~Marine,data=data4lm)))



sensoryPCA <- PCA(sensory,graph=FALSE)
Hedonic=as.data.frame(Hedonic)
dev.new();carto(sensoryPCA$ind$coord[,1:2],
                 Hedonic,
                 level=-0.5,
                 regmod=1) #6.2-7


# 6.3 JAR
jar <- read.csv("data/sensory/perfumes_jar.csv")
jar$consumer <- as.factor(jar$consumer)

summary(jar)
intensity.freq <- table(jar$product,jar$intensity)
intensity.freq

intensity.pct <- round(100*prop.table(intensity.freq,margin=1),2)
intensity.pct

percentage <- list()
for (i in seq(21)) {
attribute.freq <- table(jar$product,jar[,i+2])
percentage[[i]] <- as.matrix(100*prop.table(attribute.freq,margin=1))
}
names(percentage) <- colnames(jar)[3:23]

dev.new();barplot(t(percentage$vanilla),main='Vanilla',beside=TRUE,
        legend.text=c('Much Too Little','Too Little','JAR','Too Much',
                        'Much Too Much'))


library(car)
jar2 <- jar
for (i in 1:21){
  jar2[,i+2] <- recode(jar[,i+2], recodes="c(-2,-1)='low'; 0='jar'; c(1,2)='high'",as.factor=TRUE)
  jar2[,i+2] <- relevel(jar2[,i+2],ref="jar")

  }
  summary(jar2)

 


options(contrasts=c("contr.treatment","contr.poly"))
jar2$product <- as.factor(jar2$product)
product.names <- levels(jar2[,2])

penalty <- list()
for (j in 1:length(product.names)) {

  penalty[[j]] <- list()

jar2.p <- jar2[jar2[,2]==product.names[j],]
  
for (i in 1:21) {
penalty[[j]][[i]] <- coef(summary(lm(jar2.p$liking~jar2.p[,i+2])))
                             
}

names(penalty[[j]]) <- colnames(jar2)[3:23]

  }

names(penalty) <- product.names  
round(penalty$Shalimar$intensity,3)
  
attribute <- colnames(jar)[3:23]
ID=sample(product.names,1);ID
ID.no=which(ID==names(penalty))
dev.new();plot(0,0,type="n",xlim=c(0,65),ylim=c(-0.3,3),
       xlab="Proportion (%)",ylab="Penalty",main=paste0("Product: ",ID))  #Fig 6.3-4
for (i in 1:21){
 points(sum(percentage[[i]][rownames(percentage[[i]])==ID,4:5]),
         abs(penalty[[ID.no]][[i]][2,1]),pch=20,col="blue")
 if (penalty[[ID.no]][[i]][2,4]<=0.05){
 text(sum(percentage[[i]][rownames(percentage[[i]])==ID,4:5]),
       abs(penalty[[ID.no]][[i]][2,1])+0.05,paste(attribute[i],"*",sep=""),
       col="blue",cex=0.8)
 } else {
 text(sum(percentage[[i]][rownames(percentage[[i]])==ID,4:5]),
       abs(penalty[[ID.no]][[i]][2,1])+0.05,attribute[i],col="blue",cex=0.8)
}
 points(sum(percentage[[i]][rownames(percentage[[i]])==ID,1:2]),
         abs(penalty[[ID.no]][[i]][3,1]),pch=20,col="lightblue")
 if (penalty[[ID.no]][[i]][3,4]<=0.05) {
 text(sum(percentage[[i]][rownames(percentage[[i]])==ID,1:2]),
       abs(penalty[[ID.no]][[i]][3,1])+0.05,paste(attribute[i],"*",sep=""),
       col="lightblue",cex=0.8)
 } else {
 text(sum(percentage[[i]][rownames(percentage[[i]])==ID,1:2]),
       abs(penalty[[ID.no]][[i]][3,1])+0.05,attribute[i],col="lightblue",
       cex=0.8)
}
}
abline(v=20,lwd=2,lty=2,col="red")
legend("topleft",bty="n",legend=c("Low 不足","High 偏多"),col=
           c("lightblue","blue"),pch=20,cex=0.8)  
  
  
percentage=lapply(percentage,round,2)
  

# 6.4 IPM: Ideal Profile Method
ideal <- read.csv("data/sensory/perfumes_ideal.csv")
ideal[,1] <- as.factor(ideal[,1])
summary(ideal[,1:6])



library(SensoMineR)
dev.new();res.MultiId <- MultiIdeal(ideal,col.p=2,col.j=1,id.recogn="id_") #6.4-2
round(res.MultiId,3)


dev.new(); res.IdMap <- IdMap(ideal,col.p=2,col.j=1,col.lik=ncol(ideal),id.recogn='id_')
names(res.IdMap)
dev.new();plot(res.IdMap,choix="IdMap")#6.4-3


dev.new();plot(res.IdMap,
               xlim=c(-10,10),
               ylim=c(-10,10),
               color=TRUE,
               inverse=TRUE)#6.4-3
t(res.IdMap$ideal$profiles)

res.IdMap$ideal$pct.conso


ideal.ref <- round(res.IdMap$ideal$profiles,2)

att.no <- seq(from=3,to=ncol(ideal)-1,by=2)
senso.data <- averagetable(ideal[,c(1,2,att.no)],
                                formul="~product+consumer",firstvar=3)

id=sample(product.names,1) #id="Aromatics Elixir"
difference <- t(round(ideal.ref-senso.data[id,],2))
difference=data.frame(value=difference[order(difference),])

library(ggplot2)
gg.dat=data.frame(attributes=rownames(difference),
                  value=difference$value)
dev.new();ggplot(gg.dat, aes(x=value, y=reorder(attributes,-value),fill=attributes)) + 
  geom_bar(stat = "identity", show.legend = FALSE) +
  xlab(id)
