R/utils.R

In perbrock/sensR: Thurstonian Models for Sensory Discrimination

getPguess <-
    function(method = c("duotrio", "tetrad", "threeAFC",
             "twoAFC", "triangle", "hexad", "twofive", "twofiveF"),
             double = FALSE)
{
    ## get guessing probability for all protocols
    method <- match.arg(method)
    double <- as.logical(double[1L])
    pg <- switch(method,
                 duotrio = 1/2,
                 twoAFC = 1/2,
                 threeAFC = 1/3,
                 triangle = 1/3,
                 tetrad = 1/3,
                 hexad = 1/10,
                 twofive = 1/10,
                 twofiveF = 2/5)
    if(double) pg^2 else pg
}

getFamily <-
    function(method = c("duotrio", "tetrad", "threeAFC",
             "twoAFC", "triangle", "hexad", "twofive", "twofiveF"),
             double = FALSE)
{
  method <- match.arg(method)
  double <- as.logical(double[1L])
  if(method %in% c("hexad", "twofive", "twofiveF") && double)
      stop("'double' method for 'hexad', 'twofive' and 'twofiveF' is not yet implemented")
  if(!double) {
      fam <- switch(method,
                    duotrio = duotrio(),
                    tetrad = tetrad(),
                    triangle = triangle(),
                    twoAFC = twoAFC(),
                    threeAFC = threeAFC(),
                    hexad = hexad(),
                    twofive = twofive(),
                    twofiveF = twofiveF())
  } else {
      fam <- switch(method,
                    duotrio = doubleduotrio(),
                    tetrad = doubletetrad(),
                    triangle = doubletriangle(),
                    twoAFC = doubletwoAFC(),
                    threeAFC = doublethreeAFC())
  }
  fam
}

rescale <-
  function(pc, pd, d.prime, std.err,
           method = c("duotrio", "tetrad", "threeAFC", "twoAFC",
             "triangle", "hexad", "twofive", "twofiveF"),
           double = FALSE)
{
  m <- match.call(expand.dots = FALSE)
  m[[1]] <- as.name("list")
  m <- eval.parent(m) # evaluate the *list* of arguments
  arg <- c("pc", "pd", "d.prime")
  isPresent <- sapply(arg, function(arg) !is.null(m[[arg]]))
  if(sum(isPresent) != 1)
    stop("One and only one of pc, pd and d.prime should be given")
  method <- match.arg(method)
  double <- as.logical(double[1])
  if(method %in% c("hexad", "twofive", "twofiveF") && double)
      stop("'double' method for 'hexad', 'twofive' and 'twofiveF' is not yet implemented")
  Pguess <- getPguess(method=method, double=double)
  par <- arg[isPresent]
  if(!is.null(se <- m$std.err)) {
    stopifnot(is.numeric(se) && length(se) == length(m[[par]]))
    stopifnot(all(se[!is.na(se)] > 0))
  }
  if(par == "pc") {
    pc <- m[[par]]
    stopifnot(is.numeric(pc) && all(pc >= 0) && all(pc <= 1))
    tooSmall <- pc < Pguess
    pc[tooSmall] <- Pguess
    pd <- pc2pd(pc, Pguess)
    d.prime <- psyinv(pc, method = method, double=double)
    if(!is.null(se)) {
      se.pc <- se
      se.pc[tooSmall] <- NA
      se.pd <- se.pc / (1 - Pguess)
      se.d.prime <-
          se.pc / psyderiv(d.prime, method = method, double = double)
    }
  }
  if(par == "pd") {
    pd <- m[[par]]
    stopifnot(is.numeric(pd) && all(pd >= 0) && all(pd <= 1))
    pc <- pd2pc(pd, Pguess)
    d.prime <- psyinv(pc, method = method, double = double)
    if(!is.null(se)) {
      se.pd <- se
      se.pc <- se.pd * (1 - Pguess)
      se.d.prime <-
          se.pc / psyderiv(d.prime, method = method, double = double)
    }
  }
  if(par == "d.prime") {
    stopifnot(is.numeric(d.prime) && all(d.prime >= 0))
    d.prime <- m[[par]]
    pc <- psyfun(d.prime, method = method, double = double)
    pd <- pc2pd(pc, Pguess)
    if(!is.null(se)) {
      se.d.prime <- se
      se.pc <- se * psyderiv(d.prime, method = method, double = double)
      se.pd <- se.pc / (1 - Pguess)
    }
  }
  coef <- data.frame(pc = pc, pd = pd, d.prime = d.prime)
  res <- list(coefficients = coef)
  if(!is.null(se))
    res$std.err <- data.frame(pc = se.pc, pd = se.pd,
                              d.prime = se.d.prime)
  res <- c(res, list(method=method, double=double))
  class(res) <- "rescale"
  return(res)
}

print.rescale <- function(x, digits = getOption("digits"), ...)
{
    txt <- if(x$double) {
        paste("\nEstimates for the double", x$method, "protocol:\n", sep = " ")
    } else {
        paste("\nEstimates for the", x$method, "protocol:\n", sep = " ")
    }
  cat(txt)
  print(coef(x))
  if(!is.null(x$std.err)) {
    cat("\nStandard errors:\n")
    print(x$std.err)
  }
  return(invisible(x))
}

pc2pd <- function(pc, Pguess)
### Maps pc to pd

### arg: pc: numeric vector; 0 <= pc <= 1
###      Pguess: the guessing probability; numeric scalar,
###              0 <= pc <= 1
### res: pd: numeric vector; 0 <= pc <= 1
{
  stopifnot(is.numeric(Pguess) && length(Pguess) == 1 &&
            Pguess >= 0 && Pguess <= 1)
  stopifnot(is.numeric(pc) && all(pc >= 0) && all(pc <= 1))
  pd <- (pc - Pguess) / (1 - Pguess)
  pd[pc <= Pguess] <- 0
  names(pd) <- names(pc)
  return(pd)
}

pd2pc <- function(pd, Pguess) {
### Maps pd to pc

### arg: pd: numeric vector; 0 <= pc <= 1
###      Pguess: the guessing probability; numeric scalar,
###              0 <= pc <= 1
### res: pc: numeric vector; 0 <= pc <= 1
  stopifnot(is.numeric(Pguess) && length(Pguess) == 1 &&
            Pguess >= 0 && Pguess <= 1)
  stopifnot(is.numeric(pd) && all(pd >= 0) && all(pd <= 1))
  pc <- Pguess + pd * (1 - Pguess)
  names(pc) <- names(pd)
  return(pc)
}

psyfun <-
  function(d.prime,
           method = c("duotrio", "tetrad", "threeAFC", "twoAFC",
             "triangle", "hexad", "twofive", "twofiveF"),
           double = FALSE)
### Maps d.prime to pc for sensory discrimination protocols

### arg: d.prime: non-negative numeric vector
### res: pc: numeric vector
{
  method <- match.arg(method)
  double <- as.logical(double[1L])
  stopifnot(all(is.numeric(d.prime)) && all(d.prime >= 0),
            length(double) == 1L && is.logical(double))
  psyFun <- getFamily(method=method, double=double)$linkinv
  pc <- numeric(length(d.prime))
### Extreme cases are not handled well in the links, so we need:
  OK <- d.prime < Inf
  if(sum(OK) > 0)
    pc[OK] <- psyFun(d.prime[OK])
  pc[!OK] <- 1
  names(pc) <- names(d.prime)
  return(pc)
}

psyinv <- function(pc,
                   method = c("duotrio", "tetrad", "threeAFC", "twoAFC",
                   "triangle", "hexad", "twofive", "twofiveF"),
                   double = FALSE)
### Maps pc to d.prime for sensory discrimination protocols

### arg: pc: numeric vector; 0 <= pc <= 1
### res: d.prime: numeric vector
{
  method <- match.arg(method)
  double <- as.logical(double[1L])
  stopifnot(all(is.numeric(pc)) && all(pc >= 0) && all(pc <= 1),
            length(double) == 1L && is.logical(double))
  psyInv <- getFamily(method=method, double=double)$linkfun
  d.prime <- numeric(length(pc))
### Extreme cases are not handled well in the links, so we need:
  OK <- pc < 1
  if(sum(OK) > 0)
    d.prime[OK] <- psyInv(pc[OK])
  d.prime[!OK] <- Inf
  names(d.prime) <- names(pc)
  return(d.prime)
}


psyderiv <-
  function(d.prime,
           method = c("duotrio", "tetrad", "threeAFC", "twoAFC",
             "triangle", "hexad", "twofive", "twofiveF"),
           double = FALSE)
### Computes the derivative of the psychometric functions at some
### d.prime for sensory discrimination protocols.

### arg: d.prime: non-negative numeric vector
### res: pc: numeric vector
{
  method <- match.arg(method)
  double <- as.logical(double[1L])
  stopifnot(all(is.numeric(d.prime)) && all(d.prime >= 0),
            length(double) == 1L && is.logical(double))
  psyDeriv <- getFamily(method=method, double=double)$mu.eta
  Deriv <- numeric(length(d.prime))
### Extreme cases are not handled well in the links, so we need:
  OK <- (d.prime > 0 & d.prime < Inf)
  if(sum(OK) > 0)
    Deriv[OK] <- psyDeriv(d.prime[OK])
  Deriv[d.prime == 0] <- NA
  Deriv[d.prime == Inf] <- 0
  names(Deriv) <- names(d.prime)
  return(Deriv)
}

## findcr <-
##   function (sample.size, alpha = .05, p0 = .5, pd0 = 0,
##             type = c("difference", "similarity"))
## {
##   ## Find the critical value of a one-tailed binomial test.
##   type <- match.arg(type)
##   ss <- sample.size
##   if(ss != trunc(ss) | ss <= 0)
##     stop("'sample.size' has to be a positive integer")
##   if(alpha <= 0 | alpha >= 1)
##     stop("'alpha' has to be between zero and one")
##   if(p0 <= 0 | p0 >= 1)
##     stop("'p0' has to be between zero and one")
##   if(pd0 < 0 | pd0 > 1)
##     stop("'pd0' has to be between zero and one")
##   ## Core function:
##   i <- 0
##   if(type == "difference") {
##     while (1 - pbinom(i, ss, pd0 + p0*(1-pd0)) > alpha) i <- i + 1
##     i + 1
##   }
##   else {
##     while(pbinom(i, ss, pd0 + p0*(1-pd0)) < alpha) i <- i + 1
##     i - 1
##   }
## }

test.crit <-
  function(xcr, sample.size, p.correct = 0.5, alpha = 0.05, test)
### Is xcr the critical value of a one-tailed binomial test?
### Result: boolean

### OBS: there is deliberately no requirement that xcr should be
### positive or less than sample.size.
{
  if(test %in% c("difference", "greater")) ## alternative is "greater"
    ((1 - pbinom(q = xcr - 1, size = sample.size, prob = p.correct) <= alpha) &&
     (1 - pbinom(q = xcr - 2, size = sample.size, prob = p.correct) > alpha))
  else if(test %in% c("similarity", "less")) ## alternative is "less"
    ((pbinom(q = xcr, size = sample.size, prob = p.correct) <= alpha) &&
     (pbinom(q = xcr + 1, size = sample.size, prob = p.correct) > alpha))
  else
    stop("unknown 'test' argument")
}

findcr <-
  function(sample.size, alpha = 0.05, p0 = 0.5, pd0 = 0,
           test = c("difference", "similarity"))
### Find the critical value of a one-tailed binomial
### test. "difference" means a "greater" alternative hypothesis and
### "similarity" means a "less" alternative hypothesis.

### FIXME: What should this function do/return if the critical value
### is larger than the sample size? Or when it is negative as can
### happen with similarity? Examples:
### (xcr <- findcr(sample.size = 1, p0 = psyfun(1, "twoAFC"))) ## 2
### (xcr <- findcr(sample.size = 1, test = "similarity")) ## -1
### This means that there is no number large/small enough for this
### sample size that could give a significant p-value. Maybe this
### should just be a deliberate feature.
{
  ## match and test arguments:
  test <- match.arg(test)
  ss <- sample.size
### FIXME: Does this test work as intented?
  if(ss != trunc(ss) | ss <= 0)
    stop("'sample.size' has to be a positive integer")
  if(alpha <= 0 | alpha >= 1)
    stop("'alpha' has to be between zero and one")
  if(p0 <= 0 | p0 >= 1)
    stop("'p0' has to be between zero and one")
  if(pd0 < 0 | pd0 > 1)
    stop("'pd0' has to be between zero and one")
  ## core function:
  pc <- pd2pc(pd0, p0)
  if(test == "difference") {
    crdiff <-  function(cr)
      1 - pbinom(q = cr - 1, size = ss, prob = pc) - alpha
    interval <- c(0, ss + 2) ## deliberately outside allowed range
  }
  else if(test == "similarity") {
    crdiff <- function(cr)
      pbinom(q = cr + 1, size = ss, prob = pc) - alpha
    interval <- c(-2, ss) ## deliberately outside allowed range
  }
  else ## should never occur
    stop("'test' not recognized")
  xcr <- round(uniroot(crdiff, interval = interval)$root)
  ## is xcr the critical value?:
  is.crit <- test.crit(xcr = xcr, sample.size = ss, p.correct = pc,
                       alpha = alpha, test = test)
  if(is.crit) return(xcr)
  ## if uniroot fails, then do a simple search around the vicinity of
  ## the result from uniroot:
  max.iter <- 20 ## avoid infinite loop
  xcr <- delimit(xcr - 10, lower = -1)
  i <- 0
  if(test == "difference") {
    while(1 - pbinom(q = xcr + i, size = ss, prob = pc) > alpha) {
      if(i > max.iter || xcr + i > ss) break
      i <- i + 1
    }
    xcr <- xcr + i + 1
  }
  if(test == "similarity") {
    while(pbinom(q = xcr + i, size = ss, prob = pc) < alpha) {
      if(i > max.iter || xcr + i > ss) break
      i <- i + 1
    }
    xcr <- xcr + i - 1
  }
  ## is xcr now the critical value?:
  is.crit <- test.crit(xcr = xcr, sample.size = ss, p.correct = pc,
                       alpha = alpha, test = test)
  if(is.crit) return(xcr)
  else stop("Failed to find critical value")
}

delimit <- function(x, lower, upper, set.na = FALSE)
### Sets the values of x < lower to lower and values of x > upper to
### upper. If set.na is TRUE values are set to NA. If both lower and
### upper are supplied, the (lower < upper) has to be TRUE.
{
  m <- match.call()
  m[[1]] <- as.name("list")
  m <- eval.parent(m)
  if(!is.null(m$lower) && !is.null(m$upper))
    stopifnot(m$lower < m$upper)
  if(!is.null(m$lower))
    x[x < m$lower] <- ifelse(set.na, NA, m$lower)
  if(!is.null(m$upper))
    x[x > m$upper] <- ifelse(set.na, NA, m$upper)
  return(x)
}

normalPvalue <-
### Computes the p-value for a statistic that follows a standard
### normal distribution under the null hypothesis.

### Arguments:
### statistic - a numerical (vector?)
### alternative - the type of alternative hypothesis

### Value:
### the p-value, possibly a vector.
  function(statistic, alternative = c("two.sided", "less", "greater"))
{
  alternative <- match.arg(alternative)
  stopifnot(all(is.finite(statistic)))
  p.value <-
    switch(alternative,
           "greater" = pnorm(statistic, lower.tail = FALSE),
           "less" = pnorm(statistic, lower.tail = TRUE),
           "two.sided" = 2 * pnorm(abs(statistic), lower.tail = FALSE))
  return(p.value)
}


## Do not partially match arguments.
## If possible give functions explicitly named  arguments - preferably
##   with default values.
## Value readability over speed.
## Value accuracy over speed.
## Use small functions with conceptual - easy-to-understand tasks.
##