R/cata.R

In cata: Analysis of Check-All-that-Apply (CATA) Data

#' Cochran's Q test
#'
#' Conduct Cochran's Q test assuming equal columns proportions for matched binary
#' responses versus the alternative hypothesis of unequal column proportions.
#' @name cochranQ
#' @usage cochranQ(X, quiet = FALSE, digits = getOption("digits"))
#' @param X matrix of \eqn{I} assessors (rows) and \eqn{J} products (columns)
#' where values are \code{0} (not checked) or \code{1} (checked)
#' @param quiet if \code{FALSE} (default) then it prints information related to 
#' the test; if \code{TRUE} it returns only the test statistic (\code{Q})
#' @param digits for rounding
#' @return Cochran's Q test results (statistic, degrees of freedom, p-value)
#' @export
#' @details
#' Method returns test statistic, degrees of freedom, and p value from Cochran's
#' Q test.
#' @encoding UTF-8
#' @author J.C. Castura
#' @seealso \code{\link[cata]{mcnemarQ}}
#' @references  
#' Cochran, W.G. (1950). The comparison of percentages in matched samples. 
#' \emph{Biometrika}, 37, 256-266, \doi{10.2307/2332378} 
#' 
#' Meyners, M., Castura, J.C., & Carr, B.T. (2013). Existing and  
#' new approaches for the analysis of CATA data. \emph{Food Quality and Preference}, 
#' 30, 309-319, \doi{10.1016/j.foodqual.2013.06.010}
#' @examples
#' # Cochran's Q test on the first 50 consumers on the first attribute ("Fresh")
#' cochranQ(bread$cata[1:50, , 1], digits=3)
#' 
#' # Same, returning only test statistics for the first 4 attributes
#' t(res <- apply(bread$cata[1:50, , 1:4], 3, cochranQ, quiet=TRUE, digits=3))
cochranQ <- function(X, quiet = FALSE, digits = getOption("digits")){
  if(length(X) == 1){
    X <- matrix(X, nrow = 1) # 1 assessor
  }
  X <- matrix(as.integer(X), nrow = nrow(X), ncol = ncol(X), dimnames=dimnames(X))
  Jn <- ncol(X)
  X <- X[which((rowSums(X) < Jn) & (rowSums(X) > 0)), , drop=FALSE] 
  if(nrow(X) == 1){
    if(!quiet){
      print("Only 1 assessor")
    }
    return(c(Q = NA, df = Jn - 1, p.value = NA))
  }
  if(sum(X) == 0){
    if(!quiet){
      print("No variability in X")
    }
    return(c(Q = NA, df = Jn - 1, p.value = NA))
  }
  Jn <- ncol(X)
  # Cochran's Q test
  numQ <- Jn * (Jn-1) * (sum(colSums(X)^2) - (sum(X)^2)/Jn)
  denomQ <-  Jn * sum(X) - sum(rowSums(X)^2)
  Q = numQ / denomQ
  # chance-corrected effect size
  # In <- nrow(X) # blocks
  # if(In == 1){ 
  #   udelta <- delta <- 0
  #   R = NA
  # } else { # this condition is needed to calculate delta
  #   Ipairs <- utils::combn(In, 2)
  #   dnum <- sum(apply(Ipairs, 2, function(pindx) {
  #     sum(abs(X[pindx[1], ] - X[pindx[2], ]))
  #   }))
  #   delta <- dnum/(Jn * ncol(Ipairs))
  #   p_i <- apply(X, 1, mean)
  #   udelta <- 1/ncol(Ipairs) * ((sum(p_i)*(In - sum(p_i))) - sum((p_i*(1-p_i))))
  #   R = 1 - delta/udelta
  # }
  # report output
  o <- c(Q  = Q, df = Jn - 1, 
         p.value = stats::pchisq(Q, Jn - 1, lower.tail = FALSE))
  if(!quiet){
    oprint <- o
    oprint[c("Q", "p.value")] <- .digf(o[c("Q", "p.value")], digits = digits)
    cat(paste(c("",
                "Cochran's Q test", 
                "----------------",
                "H0: Citation rates equal for all products (columns)",
                "H1: Citation rates not equal for all products",
                ""), 
              collapse = '\n'))
    print(oprint, quote = FALSE)
    # cat(paste(c("",
    #             "Chance-corrected effect size: ", o["effect.size"], ""), 
    #           collapse = '\n'))
  }
  invisible(o)
}

#' McNemar's test
#'
#' Pairwise tests are conducted using the two-tailed binomial test. These tests
#' can be conducted after Cochran's Q test. 
#' @name mcnemarQ
#' @usage mcnemarQ(X, quiet = FALSE, digits = getOption("digits"))
#' @param X matrix of \eqn{I} assessors (rows) and \eqn{J} products (columns)
#' where values are \code{0} (not checked) or \code{1} (checked)
#' @param quiet if \code{FALSE} (default) then it prints information related to 
#' the test; if \code{TRUE} it returns only the test statistic (\code{Q})
#' @param digits for rounding
#' @return Test results for all McNemar pairwise tests conducted via the 
#' binomial test
#' @export
#' @encoding UTF-8
#' @author J.C. Castura
#' @seealso \code{\link[cata]{cochranQ}}
#' @references  
#' Cochran, W.G. (1950). The comparison of percentages in matched samples. 
#' \emph{Biometrika}, 37, 256-266. \doi{10.2307/2332378}
#' 
#' McNemar, Q. (1947). Note on the sampling error of the difference between 
#' correlated proportions or percentages. \emph{Psychometrika}, 12(2), 153-157.
#' \doi{10.1007/BF02295996}
#' 
#' Meyners, M., Castura, J.C., & Carr, B.T. (2013). Existing and  
#' new approaches for the analysis of CATA data. \emph{Food Quality and 
#' Preference}, 30, 309-319, \doi{10.1016/j.foodqual.2013.06.010}
#' 
#' @examples
#' # McNemar's exact pairwise test for all product pairs
#' # on the first 50 consumers and the first attribute ("Fresh")
#' mcnemarQ(bread$cata[1:50, , 1])
#' 
#' # Same, returning only results for the first 4 attributes
#' (res <- apply(bread$cata[1:50, , 1:4], 3, mcnemarQ, quiet=TRUE, simplify=FALSE))
mcnemarQ <- function(X, quiet = FALSE, digits = getOption("digits")){
  if(is.vector(X)){
    # assume data come from only 1 assessor
    # proceed to return results structured in the usual way
    X <- matrix(X, nrow = 1)
  }
  if(length(X) == 1){
    return(NA) # not much can be done when 1 assessor evaluates only 1 product
  }
  X <- matrix(as.integer(X), nrow = nrow(X), ncol = ncol(X), dimnames=dimnames(X))
  Jn <- ncol(X)
  X <- X[which((rowSums(X) < Jn) & (rowSums(X) > 0)), , drop=FALSE] # also deals with missings!
  if(sum(X) == 0){
    if(!quiet){
      print("No variability in X")
    }
    return(NA)
  }
  res <- matrix(NA, nrow = choose(Jn, 2), ncol = 6,
                dimnames = list(NULL, c("index.1", "index.2", "b", "c", 
                                        "proportion", "p.value")))
  res[, 1:2] <- t(utils::combn(Jn, 2))
  if(Jn == 2){
    res[, 3:4] <- colSums(X)
    res[, 5] <- res[,3] / sum(res[,3:4])
  } else {
    res[,3:4] <- apply(barray(X, values = "bc"), 2:3, sum)
    res[,5] <- res[, 3] / rowSums(res[, 3:4])
  }
  res[, 6] <- mapply(function(x, y){ 
    2 * sum(stats::dbinom(0:min(x, sum(x, y) - x), size = sum(x, y), prob = .5))}, 
    res[, 3], res[, 4])
  res[, 6] <- mapply(min, res[, 6], 1)
  if(!quiet){
    cat(paste(c("",
                "McNemar's pairwise test", 
                "-----------------------", "",
                "H0: Citation rates equal for both products",
                "H1: Citation rates not equal for both products",
                "", ""), 
              collapse = '\n'))
    res.print <- res
    res.print[, 5:6] <- .digf(res.print[, 5:6], digits = digits)
    print(as.data.frame(res.print, row.names = ""))
  }
  invisible(res)
}

#' Convert CATA binary data to CATA differences array
#'
#' Converts a three-dimensional array (\eqn{I} assessors, \eqn{J} products, 
#' \eqn{M} attributes) to an output array of \eqn{J(J-1)/2} product comparisons.
#'  
#' @name barray
#' @usage barray(X, values = "bc", type.in = "binary", type.out = "binary")
#' @param X three-dimensional array (\eqn{I} assessors, \eqn{J}
#' products, \eqn{M} attributes) where values are \code{0} (not checked) 
#' or \code{1} (checked)
#' @param values \code{"bc"} returns a four-way array separating two outcomes: 
#' \code{b} and \code{c}; \code{"bc-diff"} (default) returns a three-way array 
#' of differences between outcomes \code{b-c}; \code{"abcd"} returns a four-way 
#' array separating four outcomes: \code{a}, \code{b}, \code{c}, \code{d}.
#' @param type.in type of data submitted; default (\code{binary}) may be set to
#' \code{ordinal} or \code{scale}.
#' @param type.out currently only \code{binary} is implemented
#' @return for \code{"values = bc-diff"}, a three-way array of \eqn{I} ]assessors, 
#' \eqn{J(J-1)/2} product comparisons, and \eqn{M} attributes; otherwise, a four-way 
#' array of \eqn{I} assessors, \eqn{J(J-1)/2} product comparisons, binary outcomes 
#' for the paired difference, and \eqn{M} attributes.
#' @export
#' @encoding UTF-8
#' @author J.C. Castura
#' @references Castura, J.C., Meyners, M., Varela, P., & Næs, T. (2022). 
#' Clustering consumers based on product discrimination in check-all-that-apply 
#' (CATA) data. \emph{Food Quality and Preference}, 104564. 
#' \doi{10.1016/j.foodqual.2022.104564}.
#' 
#' @examples
#' # Get the 4d array of CATA differences for the first 8 consumers
#' b <- barray(bread$cata[1:8,,])
barray <- function(X, values = "bc", type.in = "binary", type.out = "binary"){
  abcd_1attribute <- function(X, type.in = "binary"){
    .abcd_1pair <- function(x, y, type.in = "binary"){
      x <- c(x) # product 'x' (1 response per consumer)
      y <- c(y) # product 'y' (1 response per consumer)
      if(type.in == "binary"){
        return(list(a = sum(x+y == 2), b = sum(x-y == 1),
                    c = sum(x-y == -1), d = sum(x+y == 0)))
      }
      if(type.in %in% c("scale", "ordinal")){
        return(list(a = sum(all(x == y, x != 0)), 
                    b = sum(x > y),
                    c = sum(x < y), 
                    d = sum(all(x == y, x == 0))))
      }
    }
    if(is.vector(X)){
      X <- matrix(X, nrow = 1)
    } else {
      X <- as.matrix(X) # assessors x products
    }
    pair.items <- utils::combn(ncol(X), 2)
    this.pair <- matrix(0, nrow = ncol(pair.items), ncol = 4, 
                        dimnames = list(NULL, letters[1:4]))
    for(i in 1:ncol(pair.items)){
      this.pair[i, ] <- unlist(.abcd_1pair(X[,pair.items[1,i]], 
                                           X[,pair.items[2,i]], type.in = type.in))
    }
    return(this.pair)
  }
  bc_1attribute <- function(X, type.in = "binary"){
    .bc_1pair <- function(x, y, type.in = "binary"){
      x <- c(x) # product 'x' (1 response per consumer)
      y <- c(y) # product 'y' (1 response per consumer)
      if(type.in == "binary"){
        return(list(b = sum(x - y == 1), c = sum(x - y == -1)))
      }
      if(type.in %in% c("scale", "ordinal")){
        return(list(b = sum(x > y),
                    c = sum(x < y)))
      }
    }
    if(is.vector(X)){
      X <- matrix(X, nrow = 1)
    } else {
      X <- as.matrix(X) # assessors x products
    }
    pair.items <- utils::combn(ncol(X), 2)
    this.pair <- matrix(0, nrow = ncol(pair.items), ncol = 2, 
                        dimnames = list(NULL, letters[2:3]))
    for(i in 1:ncol(pair.items)){
      this.pair[i, ] <- unlist(.bc_1pair(X[,pair.items[1,i]], 
                                         X[,pair.items[2,i]], type.in = type.in))
    }
    return(this.pair)
  }
  dimX <- dim(X)
  namesX <- dimnames(X)
  if(length(dimX) == 2){
    oX <- X
    X <- array(NA, c(dimX[1], dimX[2], 1), dimnames = list(
      namesX[[1]], namesX[[2]], 1))
    X[,,1] <- oX
    dimX <- dim(X)
    namesX <- dimnames(X)
  }
  if (values %in% "bc-diff") {
    out <- array(0, c(dimX[1], choose(dimX[2], 2), 2, dimX[3]), dimnames = list(
      namesX[[1]], apply(t(utils::combn(dimX[2], 2)), 1, paste0, collapse = "_"), 
      letters[2:3], namesX[[3]]))
    for (i in 1:dimX[1]){
      for(m in 1:dimX[3]){
        this.data <- X[i,,m]
        out[i,,,m] <- bc_1attribute(this.data, type.in = type.in)
      }
    }
    return(array(out[,,1,, drop = FALSE] - out[,,2,, drop = FALSE], dim = dim(out)[c(1,2,4)]))
  }
  if (values %in% "bc") {
    out <- array(0, c(dimX[1], choose(dimX[2], 2), 2, dimX[3]), dimnames = list(
      namesX[[1]], apply(t(utils::combn(dimX[2], 2)), 1, paste0, collapse = "_"), 
      letters[2:3], namesX[[3]]))
    for (i in 1:dimX[1]){
      for(m in 1:dimX[3]){
        this.data <- X[i,,m]
        out[i,,,m] <- bc_1attribute(this.data, type.in = type.in)
      }
    }
    return(out)
  }
  if (values %in% "abcd") {
    out <- array(0, c(dimX[1], choose(dimX[2], 2), 4, dimX[3]), dimnames = list(
      namesX[[1]], apply(t(utils::combn(dimX[2], 2)), 1, paste0, collapse = "_"), 
      letters[1:4], namesX[[3]]))
    for (i in 1:dimX[1]){
      for(m in 1:dimX[3]){
        this.data <- X[i,,m]
        out[i,,1:4,m] <- abcd_1attribute(this.data, type.in = type.in)
      }
    }
    return(out)
  }
  if(!(values %in% c("bc-diff", "bc", "abcd"))){
    return(print(paste("Invalid option: values =", values), quote = FALSE))
  }
}


#' Convert 3d array of CATA differences
#'
#' products, \eqn{M} attributes) to a four-dimensional array of product
#' comparisons (\eqn{I} assessors, \eqn{J(J-1)/2}
#' product comparisons, two outcomes (of type \code{b} or \code{c}), \eqn{M} 
#' attributes)
#'  
#' @name bcdiff
#' @usage bcdiff(X)
#' @param X three-dimensional array (\eqn{I} assessors, \eqn{J}
#' products, \eqn{M} attributes) where values are \code{0} (not checked) 
#' or \code{1} (checked)
#' @return three-way array of \eqn{I} ]assessors, \eqn{J(J-1)/2} product 
#' comparisons, and \eqn{M} attributes.
#' @export
#' @encoding UTF-8
#' @author J.C. Castura
#' @references Castura, J.C., Meyners, M., Varela, P., & Næs, T. (2022). 
#' Clustering consumers based on product discrimination in check-all-that-apply 
#' (CATA) data. \emph{Food Quality and Preference}, 104564. 
#' \doi{10.1016/j.foodqual.2022.104564}.
#' 
#' @examples
#' # Get the 3d array of CATA differences for the first 8 consumers
#' head(bcdiff <- bcdiff(bread$cata[1:8,,]))
bcdiff <- function(X){
  return(barray(X, values = "bc-diff", type.in = "binary", type.out = "binary"))
}


#' Converts 3d array of CATA data to a wide 2d matrix format
#'
#' Converts a three-dimensional array (\eqn{I} assessors, \eqn{J}
#' products, \eqn{M} attributes) to a two-dimensional matrix
#' (\eqn{J} products, (\eqn{I} assessors, \eqn{M} attributes))
#'  
#' @name toWideMatrix
#' @usage toWideMatrix(X)
#' @param X three-dimensional array (\eqn{I} assessors, \eqn{J}
#' products, \eqn{M} attributes) where values are \code{0} (not checked) 
#' or \code{1} (checked)
#' @return matrix with \code{J} products in rows and \eqn{I} 
#' assessors \eqn{\times M} attributes in columns (data from each assessor are
#' in adjacent columns; attribute order is the same within each assessor)
#' @export
#' @encoding UTF-8
#' @author J.C. Castura
#' @examples
#' # convert CATA results from the first 8 consumers and the first 4 attributes
#' # to a wide matrix
#' toWideMatrix(bread$cata[1:8,,1:4])
toWideMatrix <- function(X){
  dimX <- dim(X)
  if(length(dimX) != 3){
    return("Function requires a 3d array")
  } 
  namesI <- dimnames(X)[[1]]
  namesJ <- dimnames(X)[[2]]
  namesM <- dimnames(X)[[3]]
  if(is.null(namesI)[1]) namesI <- 1:dimX[1]
  if(is.null(namesJ)[1]) namesJ <- 1:dimX[2]
  if(is.null(namesM)[1]) namesM <- 1:dimX[3]
  out <- array(NA, dim = c(dimX[2], prod(dimX[c(1,3)])), dimnames = list(
    namesJ, paste0(rep(namesM, times = dimX[1]), "_", rep(namesI, each = dimX[3]))))
  for(i in 1:dimX[1]){
    out[, (1+(i-1)*dimX[3]):(i*dimX[3])] <- array(X[i,,], dim = dimX[2:3])
  }
  return(out)
}

#' Converts 3d array of CATA data to a tall 2d matrix format
#'
#' Converts a three-dimensional array (\eqn{I} assessors, \eqn{J}
#' products, \eqn{M} attributes) to a two-dimensional matrix with
#' (\eqn{I} assessors, \eqn{J} products) rows and (\eqn{M} 
#' attributes) columns, optionally preceded by two columns of row headers.
#'  
#' @name toTallMatrix
#' @usage toTallMatrix(X, header.rows = TRUE)
#' @param X three-dimensional array (\eqn{I} assessors, \eqn{J}
#' products, \eqn{M} attributes) where values are \code{0} (not checked) 
#' or \code{1} (checked)
#' @param header.rows \code{TRUE} (default) includes row headers; set to
#' \code{FALSE} to exclude these headers
#' @return matrix with \eqn{I} assessors \eqn{\times J} products in rows
#' and \eqn{M} attributes in columns (optionally preceded by 2 columns of 
#' headers if \code{header.rows = TRUE})
#' @export
#' @encoding UTF-8
#' @author J.C. Castura
#' @examples
#' # convert CATA results from the first 8 consumers and the first 4 attributes
#' # to a tall matrix
#' toTallMatrix(bread$cata[1:8,,1:4])
toTallMatrix <- function(X, header.rows = TRUE){
  dimX <- dim(X)
  if(length(dimX) != 3){
    return("Function requires a 3d array")
  } 
  namesI <- dimnames(X)[[1]]
  namesJ <- dimnames(X)[[2]]
  namesM <- dimnames(X)[[3]]
  if(is.null(namesI)[1]) namesI <- 1:dimX[1]
  if(is.null(namesJ)[1]) namesJ <- 1:dimX[2]
  if(is.null(namesM)[1]) namesM <- 1:dimX[3]
  out <- matrix(NA, ncol = dimX[3], nrow = prod(dimX[1:2]),
                 dimnames = list(paste0(rep(namesI, each = dimX[2]), "_",
                                        rep(namesJ, times = dimX[1])), namesM))
  for(cc in 1:dimX[3]){
    out[, cc] <- c(aperm(array(X[,, cc], dim = dimX[1:2]), 2:1))
  }
  if(header.rows){
    out <- data.frame(assessor = as.factor(rep(namesI, each = dimX[2])),
                      product = as.factor(rep(namesJ, times = dimX[1])), out)
  }
  return(out)
}

#' Calculate \eqn{RV} Coefficient
#'
#' Calculate \eqn{RV} coefficient
#' @name rv.coef
#' @usage rv.coef(X, Y, method = 1)
#' @param X input matrix (same dimensions as \code{Y})
#' @param Y input matrix (same dimensions as \code{X})
#' @param method \code{1} (default) and \code{2} give identical \eqn{RV} coefficients
#' @return \eqn{RV} coefficient
#' @export
#' @encoding UTF-8
#' @author J.C. Castura
#' @references Robert, P., & Escoufier, Y. (1976). A unifying tool for linear 
#' multivariate statistical methods: the RV-coefficient. \emph{Journal of the Royal 
#' Statistical Society: Series C (Applied Statistics)}, 25, 257-265.
#' \doi{10.2307/2347233}
#' @examples
#' # Generate some data
#' set.seed(123)
#' X <- matrix(rnorm(8), nrow = 4)
#' Y <- matrix(rnorm(8), nrow = 4)
#' 
#' # get the RV coefficient
#' rv.coef(X, Y)
rv.coef <- function(X, Y, method = 1){
  tr <- function(X){
    sum(diag(X), na.rm = TRUE)
  }
  X <- sweep(X, 2, colMeans(X))#apply(X, 2, mean), "-")
  Y <- sweep(Y, 2, colMeans(Y))#apply(Y, 2, mean), "-")
  XX <- tcrossprod(X,X)
  YY <- tcrossprod(Y,Y)
  if(method==1){
    out <- tr(XX %*% YY) /
      sqrt(tr(XX %*% XX) %*% tr(YY %*% YY) )
  }
  if(method==2){
    out <- sum(c(XX)*c(YY), na.rm=TRUE) /
      sqrt(sum(XX^2, na.rm=TRUE)*sum(YY^2, na.rm=TRUE))
  }
  return(c(out))
}

#' Salton's cosine measure
#'
#' Calculate Salton's cosine measure
#' @name salton
#' @usage salton(X, Y)
#' @param X input matrix (same dimensions as \code{Y})
#' @param Y input matrix (same dimensions as \code{X})
#' @return Salton's cosine measure
#' @export
#' @encoding UTF-8
#' @author J.C. Castura
#' @references Salton, G., & McGill, M.J. (1983). \emph{Introduction to Modern 
#' Information Retrieval}. Toronto: McGraw-Hill.
#' @examples
#' # Generate some data
#' set.seed(123)
#' X <- matrix(rnorm(8), nrow = 4)
#' Y <- matrix(rnorm(8), nrow = 4)
#' 
#' # get Salton's cosine measure
#' salton(X, Y)
salton <- function(X, Y){
  out <- sum(c(X)*c(Y)) / sqrt(sum(X^2)*sum(Y^2))
  return(out)
}


#' Apply top-k box coding to scale data
#'
#' Apply top-k box coding to scale data. Using defaults give top-2 box (T2B) coding.
#' @name code.topk
#' @usage code.topk(X, zero.below = 8, one.above = 7)
#' @param X input matrix
#' @param zero.below default is \code{8}; values below this numeric threshold 
#' will be coded \code{0}; use \code{NULL} if there is no such threshold
#' @param one.above default is \code{7}; values above this numeric threshold 
#' will be coded \code{1}; use \code{NULL} if there is no such threshold
#' @return matrix \code{X} with top-k coding applied
#' @export
#' @encoding UTF-8
#' @author J.C. Castura
#' @references Castura, J.C., Meyners, M., Pohjanheimo, T., Varela, P., & Næs, T. (2023). 
#' An approach for clustering consumers by their top-box and top-choice responses. 
#' \emph{Journal of Sensory Studies}, e12860. \doi{10.1111/joss.12860}
#' @examples
#' # Generate some data
#' set.seed(123)
#' X <- matrix(sample(1:9, 100, replace = TRUE), nrow = 5)
#' 
#' # apply top-2 box (T2B) coding
#' code.topk(X, zero.below = 8, one.above = 7)
code.topk <- function(X, zero.below = 8, one.above = 7){
  Y <- X
  if(!is.null(zero.below) && !is.null(one.above)){
    if(zero.below != one.above + 1){
      return(
        print("zero.below should be equal to one.above+1 if both specified"))
    }
  }
  if(!is.null(zero.below)){
    Y[Y < zero.below] <- 0
  }
  if(!is.null(one.above)){
    Y[Y > one.above] <- 1
  }
  return(matrix(as.integer(Y), nrow = nrow(Y), ncol = ncol(Y), dimnames = dimnames(Y)) )
  #return(Y)
}

#' Apply top-c choices coding to a vector of scale data from a respondent
#'
#' Apply top-c choices coding to a vector of scale data from a respondent
#' @name topc
#' @usage topc(x, c = 2, coding = "B")
#' @param x input matrix
#' @param c number of top choices considered to be 'success'; other choices are 
#' considered to be 'failure' and are coded \code{0}
#' @param coding \code{"B"} (default) codes all successes as \code{1}; 
#' \code{"N"} codes all successes with their numeric coding
#' @return matrix \code{X} with top-k coding applied
#' @export
#' @encoding UTF-8
#' @author J.C. Castura
#' @references Castura, J.C., Meyners, M., Pohjanheimo, T., Varela, P., & Næs, T. (2023). 
#' An approach for clustering consumers by their top-box and top-choice responses. 
#' \emph{Journal of Sensory Studies}, e12860. \doi{10.1111/joss.12860}
#' @examples
#' # Generate some data
#' set.seed(123)
#' X <- matrix(sample(1:9, 100, replace = TRUE), nrow = 5)
#' 
#' # apply top-2 choice (T2C) coding
#' apply(X, 1, topc)
topc <- function(x, c = 2, coding = "B"){
  #anything with a value of 0 is special, so give in the max value
  y <- max(x)-x+1
  ranky <- rank(y, ties.method = "average")
  ranky.u <- sort(unique(ranky))
  it <- 0
  indx <- c()
  nc <- 0
  while(nc < c){
    it <- it + 1
    indx <- c(indx, which(ranky == ranky.u[it]))
    nc <- length(indx)
  }
  out <- x*0
  if(coding %in% "B"){
    out[indx] <- 1
  }
  if(coding %in% "N"){
    out[indx] <- x[indx]
  }
  return(out)
}