options(help_type = "html")
library(sensR)
#sessionInfo()

#################################
## Same-diff:
#################################
## nsamesame = 8, ndiffsame = 5, nsamediff = 4, ndiffdiff = 9
Anna <- samediff(8, 5, 4, 9)
Anna

## summary(Anna, profile = FALSE)
summary(Anna)

## Distributions of sensory intensity:
plot(Anna, ylim = c(0, .45))
text(-1, .35, "X")
text(3, .35, "Y")
arrows(0, .405, coef(Anna)[2], .405, length=.1, lwd = 2)
text(1, .42, expression(delta))
segments(4, .3, coef(Anna)[1] + 4, .3, lwd = 2)
text(4.6, .32, expression(tau))

## Similarity testing:
Mark <- samediff(8, 5, 11, 2)
summary(Mark)
plot(profile(Mark))

#################################
## A-not A
#################################

## sensR function to compute delta:
(ana <- AnotA(39, 100, 19, 100))
ana$coef

## Obtain standard error, CI and P-value:
confint(ana)

## Similarity testing:

## Illustration of the distributions of sensory intensity:
plot(AnotA(57, 100, 42, 100), xlim = c(-3, 3), ylim = c(0, .45))
arrows(0, .405, coef(res)[2], .405, length=.1, lwd = 2)
text(.2, .42, expression(delta))
abline(v = -coef(res)[1], lty=3)

#################################
## ROC curves and AUC
#################################
source("funs.R")
## ROC curve for delta-estimate
(delta <- ana$coef)
ROC(delta)
AUC(delta)

## computation of Sensitivity (AUC):
pnorm(AnotA(57, 100, 42, 100)$coef / sqrt(2))
pnorm(AnotA(39, 100, 19, 100)$coef / sqrt(2))
## Computation of distribution overlap:
overlap(delta, 1)

#################################
## A-not A with sureness
#################################
## read in data:
wts <- c(132, 161, 65, 41, 121, 219, 96, 99, 50, 57, 156, 650)
dat <- data.frame(sureness = factor(rep(1:6, 2), ordered=TRUE),
                  prod = factor(rep(c("ref", "test"), each = 6)),
                  freq = wts)
dat

library(ordinal)
## str(dat)
## see also ?soup and data(soup)

## Table of data:
tab <- t(xtabs(freq ~ sureness + prod, data = dat))
tab

sure_bi <- unclass(dat$sureness) >= 3
tab_bi <- xtabs(freq ~ prod + sure_bi, data=dat)
rowSums(tab_bi)
rowSums(tab)
tab_bi/10


## Conventional test for "any difference":
chisq.test(tab)

## cumulative link models:

## Equal-variances model:
fm1 <- clm(sureness ~ prod, data = dat, weights=freq,
           link = "probit")
fm1
summary(fm1) ## read off d-prime
coef(fm1) ## d-prime
fm1$beta
## Sensitivity:
pnorm(fm1$beta / sqrt(2))
## Overlap:
overlap(d=fm1$beta, sigma=1)
2 * pnorm(-fm1$beta / 2)
## Test of H_0: d' = 0:
## 1) Estimate null model:
fm0 <- clm(sureness ~ 1, data = dat, weights=freq,
           link = "probit")
## or: fm0 <- update(fm1, ~.-prod)
## 2) Test hypothesis in a likelihood ratio test:
anova(fm0, fm1)


## Unequal-variances model:
fm2 <- clm(sureness ~ prod, scale = ~prod, data = dat,
            weights = freq, link = "probit")
fm2
summary(fm2)

## Test if we need unequal variances:
anova(fm1, fm2) ## LR test

## Are products different (yes, but here is a formal test:)
anova(fm0, fm2) ## note 2 df

## Scale estimate:
exp(0.2167) # 1.24197
exp(fm2$zeta) # 1.24197
## and d-prime:
fm2$beta

## Sensitivity:
pnorm(fm2$beta, sqrt(1 + exp(fm2$zeta)^2))
## or: pnorm(fm2$beta, sqrt(1 + exp(2 * fm2$zeta)))

## Distribution overlap:
overlap(fm2$beta, exp(fm2$zeta))

#################################
## Extras:
#################################

### Figure 1:
## Normalized Profile Likelihood for delta:
x <- plot(profile(Anna, max=3), which=2)
lim <- sapply(c(0.95, .99), function(x) exp(-qchisq(x, df=1)/2) )
abline(h = lim)
text(4.5, .17, "95% limit")
text(4.5, .06, "99% limit")

## Gaussian approximation:
dlta <- coef(Anna)[2]
lines(x[,3], exp(-(x[,3] - dlta)^2/(2 * vcov(Anna)[2, 2])), lty = 2)


### Section 5.1. Example: Anna's Experiment:

## Pearson's test (without Yates' continuity correction):
(mat <- matrix(c(8, 4, 5, 9), 2))
chisq.test(mat, correct = FALSE)$p.value/2

## Likelihood Ratio (G^2) test:
names(chisq.test(mat, correct = FALSE))
Expected <- chisq.test(mat, correct = FALSE)$expected
with(chisq.test(mat, correct = FALSE),
     pchisq(2* sum(observed * log(observed/expected)), df=1,
            lower=FALSE)/2) ## or:
pchisq(2 * sum(mat * log(mat/Expected)), df=1, lower=FALSE)/2

## Pearson's test with Yates' continuity correction:
chisq.test(mat, correct = TRUE)$p.value/2

## Fisher's 'exact' test:
fisher.test(mat, alternative = "greater")$p.value

## Likelihood root statistic for delta > 0.5
sp <- with(x, splinefun(delta, nplDelta,  method = "natural"))
pnorm( sqrt(-2 * log(sp(1/2))), lower = FALSE)

## Corresponding Wald p-value:
pnorm((coef(Anna)[2] - 1/2)/sqrt(vcov(Anna)[2,2]), lower = FALSE)


#################################
## Simulation and Power:

## Finding the power of a discrimination test with a sensory delta of 2
## (alternative hypothesis) versus a null hypothesis of delta = 0 with
## a sample of size 2 x 30 and a type I level of .05. n should be higher
## for a reasonable precision:
n <- 1e3
samediffPwr(n=n, tau = 1, delta = 2, Ns = 30, Nd = 30)

## Approximate uncertainty of Power, if Power = 0.8 and n = n
## standard error of estimated Power:
sqrt(.8*.2/n)
## Approximate half 95% confidence interval width:
2*sqrt(.8*.2/n)
## So CI for estimated Power:

#################################