# Chapter 3
library(SensoMineR)

experts <- read.csv("data/sensory/perfumes_qda_experts.csv")
experts$Session <- as.factor(experts$Session)
#summary(experts)
names(experts)
dev.new();hist(experts$Greedy,col="lightgray",border="steelblue",proba=TRUE)
lines(density(experts$Greedy),lwd=2,col="blue")


{dev.new()
par(mar = c(5, 9, 2, 5))
boxplot(Floral~Product,ylab="",main="",xlab="Floral",
                  data=experts,
                  col="lightgray",
                  main="Wrapping",
                  horizontal = TRUE,las = 1)
}


out.aov <- aov(Floral~(Product+Panelist+Session)^2,data=experts)
summary(out.aov)

dev.new();SensoMineR::graphinter(experts,col.p=4,col.j=2,firstvar=16,lastvar=16,numr=1,numc=1)

dev.new();SensoMineR::graphinter(experts,col.p=4,col.j=2,firstvar=5,lastvar=8,numr=2,numc=2)



out.panelperf <- SensoMineR::panelperf(experts,firstvar=5,formul="~(Product+Panelist+Session)^2")
names(out.panelperf)

out.panelperf$p.value[,c(3,5,6)]=round(out.panelperf$p.value[,c(3,5,6)],4)

out.panelperf$variability
out.panelperf$res
out.panelperf$r2

dev.new();SensoMineR::coltable(out.panelperf$p.value[order(out.panelperf$p.value[,1]),])


out.paneliperf <- SensoMineR::paneliperf(experts,
                             formul="~(Product+Panelist+Session)^2",
                             formul.j="~Product+Session",
                             col.j=1,
                             firstvar=5,
                             synthesis=TRUE)
names(out.paneliperf)
out.paneliperf$prob.ind=round(out.paneliperf$prob.ind,4)
out.paneliperf$prob.ind


summary(aov(Spicy~Product+Session,data=subset(experts,Panelist=="CM")))

out.lm=lm(Spicy~Product+Session,data=subset(experts,Panelist=="CM"))
lmtest::waldtest(out.lm,"Product")
summary.aov(out.lm)
coef(summary(out.lm))
summary(out.lm)$sigma
sum(out.lm$residuals^2)/11

out.paneliperf$agree.ind=round(out.paneliperf$agree.ind,2)
out.paneliperf$agree.ind


aa <- averagetable(experts, formul = "~(Product+Panelist+Session)^2",
                   firstvar = 5,
                   lastvar = ncol(experts),method = "coeff")
options(contrasts=c("contr.sum","contr.sum"))
out.a <- lm(Spicy~(Product+Panelist+Session)^2-1, data=experts)
coef.a <- coef(summary(out.a))
a <- coef.a[1:12,1,drop=FALSE];a

bb <- averagetable(experts, formul = "~Product+Session", 
                   subset = (experts[,1] == "CM"), 
                   firstvar = 5, 
                   lastvar = ncol(experts))
out.b <- lm(Spicy~(Product+Session)-1,data=subset(experts,Panelist=="CM"))
coef.b <- coef(summary(out.b))
b <- coef.b[1:12,1,drop=FALSE]
diag(cor(aa, bb, use = "pairwise.complete.obs"))

cor(a,b, use = "pairwise.complete.obs")

out.paneliperf$res.ind=round(out.paneliperf$res.ind,2)
out.paneliperf$res.ind

summary(lm(Spicy~Product+Session, 
            data = subset(experts, Panelist == "CM")))$sigma


out.paneliperf$r2.ind=round(out.paneliperf$r2.ind,2)
out.paneliperf$r2.ind


out.magicsort <- magicsort(out.paneliperf$res.ind,method="median")
round(out.magicsort,2)
dev.new()
coltable(round(out.magicsort,2),level.lower=0.9328,level.upper=2.033,
         col.lower="lightgray",col.upper="darkgray")


out.PO <- panelperf(experts,firstvar=5,
                       formul="~(Product+Panelist+Session)^2+Rank")

dev.new();coltable(round(out.PO$p.value,3),col.lower="lightblue")




