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: #################################