R/consistency_cata.R

In ClustBlock: Clustering of Datasets

## =============================================================================


##' @title Test the consistency of each attribute in a CATA experiment
##'
##' @description
##' Permutation test on the agreement between subjects for each attribute in a CATA experiment
##'
##' @usage
##' consistency_cata(Data,nblo, nperm=100, alpha=0.05, printAttrTest=FALSE)
##'
##'
##' @param Data data frame or matrix. Correspond to all the blocks of variables merged horizontally
##'
##' @param nblo  numerical. Number of blocks (subjects).
##'
##' @param nperm numerical. How many permutations are required? Default: 100
##'
##' @param alpha numerical between 0 and 1. What is the threshold? Default: 0.05
##'
##' @param printAttrTest logical. Print the number of remaining attributes to be tested? Default: FALSE
##'
##'
##'
##' @return a list with:
##'         \itemize{
##'          \item consist: the consistent attributes
##'          \item no_consist: the inconsistent attributes
##'          \item pval: pvalue for each test
##'          }
##'
##'
##'
##' @keywords CATA RATA
##'
##' @references
##' Llobell, F., Giacalone, D., Labenne, A.,  Qannari, E.M. (2019).	Assessment of the agreement and cluster analysis of the respondents in a CATA experiment.	Food Quality and Preference, 77, 184-190.
##'
##'
##' @examples
##' \donttest{
##'  data(straw)
##' #with only 40 subjects
##' consistency_cata(Data=straw[,1:(16*40)], nblo=40)
##' #with all subjects
##' consistency_cata(Data=straw, nblo=114, printAttrTest=TRUE)
##' }
##'
##' @seealso   \code{\link{consistency_cata_panel}}, \code{\link{change_cata_format}}, \code{\link{change_cata_format2}}
##'
##' @export


## =============================================================================

consistency_cata <- function(Data, nblo, nperm = 100, alpha = 0.05, printAttrTest = FALSE) {
  nprod <- nrow(Data)
  nattr <- ncol(Data) / nblo
  n <- nprod
  Blocks <- rep(nattr, nblo)
  J <- rep(1:nblo, times = Blocks) # indicates which block each variable belongs to
  if (is.null(colnames(Data))) colnames(Data) <- rep(paste0("A", 1:Blocks[1]), nblo)

  # parapet for numerical Data
  for (i in 1:ncol(Data))
  {
    if (is.numeric(Data[, i]) == FALSE) {
      stop(paste("The data must be numeric (column", i, ")"))
    }
  }

  # no NA
  if (sum(is.na(Data)) > 0) {
    print("NA detected:")
    tabna <- which(is.na(Data), arr.ind = TRUE)
    print(tabna)
    stop(paste("NA are not accepted"))
  }


  # prepare data
  Xi <- array(0, dim = c(n, Blocks[1], nblo)) # array with all subjects matrices
  for (j in 1:nblo)
  {
    Xi[, , j] <- as.matrix(Data[, J == j])
  }

  Tabi <- array(0, dim = c(nprod, nblo, nattr))
  for (j in 1:nattr)
  {
    for (k in 1:nblo)
    {
      normk <- sqrt(sum(diag(tcrossprod(Xi[, j, k], Xi[, j, k]))))
      if (normk > 0) {
        Tabi[, k, j] <- Xi[, j, k] / normk
      } else {
        Tabi[, k, j] <- Xi[, j, k]
      }
    }
  }

  # S matrices:
  Sarray <- array(0, dim = c(nblo, nblo, nattr))
  for (k in 1:nattr)
  {
    S <- matrix(0, nblo, nblo)
    diag(S) <- rep(1, nblo)
    for (i in 1:(nblo - 1))
    {
      for (j in (i + 1):nblo)
      {
        S[i, j] <- sum(diag(tcrossprod(Tabi[, i, k], Tabi[, j, k])))
        S[j, i] <- S[i, j]
      }
      Sarray[, , k] <- S
    }
  }

  l1 <- NULL
  u1 <- NULL
  for (i in 1:nattr)
  {
    ei <- eigen(Sarray[, , i])
    lambda1 <- ei$values[1]
    l1 <- c(l1, lambda1)
    u1 <- cbind(u1, abs(ei$vectors[, 1]))
  }
  names(l1) <- colnames(Data[, 1:nattr])
  hom <- (l1 / nblo) * 100

  # Permutations:
  res_perm2 <- matrix(0, nperm, nattr)
  for (k in 1:nattr)
  {
    lambdaun_perm <- NULL
    eig <- NULL
    Tabi_per <- Tabi
    for (perm in 1:nperm)
    {
      for (j in 1:nblo)
      {
        Tabi_per[, j, k] <- Tabi[sample(1:n), j, k]
      }

      S_perm <- matrix(0, nblo, nblo)
      diag(S_perm) <- rep(1, nblo)
      for (i in 1:(nblo - 1))
      {
        for (j in (i + 1):nblo)
        {
          S_perm[i, j] <- sum(diag(tcrossprod(Tabi_per[, i, k], Tabi_per[, j, k])))
          S_perm[j, i] <- S_perm[i, j]
        }
      }
      eig <- c(eig, eigen(S_perm)$values[1])
    }
    res_perm2[, k] <- eig
    if (printAttrTest == TRUE) {
      print(nattr - k)
    }
  }

  # results
  quant <- NULL
  consist <- NULL
  no_consist <- NULL
  pval <- NULL
  for (k in 1:nattr)
  {
    pval[k] <- sum(l1[k] < res_perm2[, k]) / nperm
    if (pval[k] > alpha) {
      no_consist <- c(no_consist, colnames(Data)[k])
    } else {
      consist <- c(consist, colnames(Data)[k])
    }
  }

  names(pval) <- names(l1)

  return(list(consist = consist, no_consist = no_consist, pval = pval))
}