#################################
### How can sensR help us analyzing sensory similarity test data?
#################################
library(sensR)
discrim(40, 100, pd0=0.25, conf.level=0.90, method="triangle", test="similarity")

#Compare with the difference test data analysis: (only the test info changes - NOT CIs)
discrim(40, 100, conf.level=0.90, method = "triangle")

# SAME data may provide non-significant results for BOTH difference and similarity tests:
discrim(40, 100, pd0=0.2, conf.level=0.90, method="triangle", test="similarity")
discrim(40, 100, conf.level=0.90, method = "triangle")

# SAME data may provide significant results for BOTH difference and similarity tests:
discrim(80, 200, pd0=0.25, conf.level=0.90, method="triangle", test="similarity")
discrim(80, 200, conf.level=0.90, method = "triangle")

# What if similarity is defined as dprime < 0.5:
mypd0=pc2pd(psyfun(0.5,method="threeAFC"),Pguess=1/3)
#Or the same:
mypd0=rescale(d.prime=0.5,method="threeAFC")$coefficients[1,2]
mypd0
discrim(40,100,pd0=mypd0,conf.level=0.90,method="threeAFC", test="similarity")


#################################
### How can sensR help us planning sensory similarity tests?
#################################

# Assuming true pd to be zero:
discrimPwr(pdA=0,pd0=0.25, sample.size=100, alpha=0.05,pGuess =1/3,test="similarity")

# Assuming true pd to be 0.2:
discrimPwr(pdA=0.2,pd0=0.25, sample.size=100, alpha=0.05,pGuess =1/3,test="similarity")

# The critical value for the similarity test can be found by findcr: 
cr=findcr(sample.size=100, alpha = .05, p0 = 1/3, pd0 = 0.25, test = c("similarity"))
cr
discrim(41, 100, pd0=0.25, conf.level=0.90, method="triangle", test="similarity")
discrim(42, 100, pd0=0.25, conf.level=0.90, method="triangle", test="similarity")

# The power can be now explained as follows:
# First for pdA=0:
pbinom(cr,100,1/3)

# Nextfor pdA=0.2:
pdA=0.2
pcA=pd2pc(pdA,Pguess=1/3)
pcA
pbinom(cr,100,pcA)

# Power of similarity tests when similarity is defined on d.prime scale:
d.primePwr(d.primeA=0,d.prime0=0.5, sample.size=100, alpha=0.05, method="triangle", test="similarity")

# Sample size computations to meet required power:

# First on pd-scale with pd0=0.25, best case:
discrimSS(pdA=0,pd0=0.25, target.power=0.9, alpha=0.05,pGuess =1/3,
          test="similarity")

# Still on pd-scale with pd0=0.25, "bad" case:
discrimSS(pdA=.2,pd0=0.25 , target.power=0.9, alpha=0.05,pGuess=1/3,
          test="similarity")

# Then on d.prime-scale with d.prime0=0.5, best case:
d.primeSS(d.primeA=0, d.prime0=0.5, target.power=0.9, alpha=0.05,method="triangle",
          test="similarity")

# Then on d.prime-scale with d.prime0=0.5, "bad" case:
d.primeSS(d.primeA=0.4, d.prime0=0.5, target.power=0.9, alpha=0.05,method="triangle",
          test="similarity")

#################################
### How can sensR help us analyze replicated sensory difference tests?
#################################
discrim(11,12,method="triangle")
discrim(70,180,method="triangle")

# First we look at possible outcomes of replicated triangles:
# First a case where individuals are "clones" - they are NOT really different,
# AND there is NO product difference:
rs
windows()
par(mfrow=c(1,1))
rs=discrimSim(20, replicates = 12, d.prime = 0,method = "triangle",sd.indiv=0)
barplot(rs[order(rs)],col="blue",ylim=c(0,12),main="dprime=0, sd.indiv=0", names.arg=1:20,cex.names=0.8)
abline(h=mean(rs))

# Next a case where individuals are "clones" - they are NOT really different,
# but there IS a product difference, d.prime=2:
rs=discrimSim(20, replicates = 12, d.prime = 2,method = "triangle",sd.indiv=0)
barplot(rs[order(rs)],col="blue",ylim=c(0,12),main="dprime=2, sd.indiv=0", names.arg=1:20,cex.names=0.8)
abline(h=mean(rs))

# Finally a case where individuals are really different, (sd.indiv=2)
# AND there IS a product difference, d.prime=2:
rs=discrimSim(20, replicates = 12, d.prime = 2,method = "triangle",sd.indiv=2)
barplot(rs[order(rs)],col="blue",ylim=c(0,12),main="dprime=2, sd.indiv=2", names.arg=1:20,cex.names=0.8)
abline(h=mean(rs))

#  Analysis for the Table 7.4 data from the book:
x <- c(2, 2, 3, 3, 4, 4, 4, 4, 4, 4, 4, 5, 6, 10, 11)
X2 <- cbind(x, 12)
X2

summary(betabin(X2,method="triangle"))

# Naive analysis:
sum(x)
discrim(70,180,method="triangle")


# Exercise 5
# Set the situation: (make your own choices)
myN=15 #Number of individuals
myK=10 #Number of replications
myd.prime=1 # Level of product difference
mysd.indiv=1 # individual heterogeneity 

# Simulate and plot some data:
rs=discrimSim(myN, replicates = myK, d.prime = myd.prime,method = "triangle",sd.indiv=mysd.indiv)
barplot(rs[order(rs)],col="blue",ylim=c(0,myK), names.arg=1:myN,cex.names=0.8)
abline(h=mean(rs))

# Put the simulated data on the form needed for the betabin analysis:
X=matrix(c(rs,rep(myK,myN)),ncol=2)
X

# Do the betabin analysis:
summary(betabin(X,method="triangle"))