R/kTest.R

In vcoscrato/ktest: Estimated kernel densities comparisons

#' kTest
#'
#' Performs a hypothesis test for equality of distributions based on the estimated kernel densities and the permutation test.
#'
#' @param data Either a list of numeric vectors, a numeric vector (with classes parameter defined), or a stacked data frame (first column with numeric values and second column with classes).
#' @param classes Classes relative to data parameter, should be used only when data is a numeric vector.
#' @param perm Boolean indicating weather to obtain the p-value trough the permutation test or just return return the common area between densities.
#' @param B Number of permutations.
#' @param pairwise Boolean indicating weather to obtain pairwise p-values (valid only if > 2 groups). If the number of groups is big, this will make the execution much slower.
#' @param bw The bandwidth used to estimate the kernel densities.
#' @param npoints The number of points used to estimate the kernel densities.
#' @param threads Number of cores to be used for parallel computing.
#'
#' @return A list containing:
#'
#' - densities: The estimated densities for each group.
#' 
#' - labels: The group labels.
#'
#' - ca: Common area between the kernel densities.
#'
#' - pvalue: The p-value generated by the permutation test (if perm = TRUE).
#'
#' - pairwiseca: Groups pairwise common area triangular matrix (just if >= 3 groups).
#'
#' - pairwisepvalue: Groups pairwise p-value triangular matrix generated by the permutation test (if perm = TRUE, pairwise = TRUE and >= 3 groups).
#'
#' @export
#'
#' @examples data = list(x = rnorm(30), y = rexp(50), z = rpois(70, 1))
#' test = kTest(data)
#' print(test)
#' plot(test)
#' pairs(test)

kTest = function(data, classes = NULL, perm = TRUE, B = 5000, pairwise = FALSE, bw = bw.nrd0, npoints = 512, threads = detectCores() - 1) {

  output = list()

  data = convertDataStacked(data, classes)

  densities = densitiesEval(data, bw = bw(data[,1]), npoints)

  output$densities = densities$densities

  output$labels = densities$labels

  if(length(output$labels) <= 2 & pairwise == TRUE) {

    warning('pairwise not relevant when number of groups is less then 3.')

    pairwise = FALSE

  }

  output$ca = commonArea(output$densities)

  if(length(output$labels) <= 2) {

    output$pairwiseca = NULL

  } else {

    output$pairwiseca = pairwiseCommonArea(output$densities, output$labels)

  }

  if(perm) {

    output$pvalue = permTest(data, B, threads, bw = bw(data[,1]), npoints, output$ca)

    if(pairwise) {

      output$pairwisepvalue = pairwisePermTest(data, output$labels, B, threads, bw, npoints, output$pairwiseca)

    } else {

      output$pairwisepvalue = NULL

    }

  } else {

    output$pairwisepvalue = NULL

    output$pvalue = NULL

  }

  class(output) = c('kTest', 'densities')

  return(output)
}


#' @export

print.kTest = function(output) {

  cat(paste('\n', length(output$labels), 'densities kTest results:\n\n'))

  cat(paste('- Common area between all densities:', round(output$ca, 4)))

  if(!is.null(output$pvalue)) {

    cat(paste('\n\n- p-value for H0 (All densities are equal):', round(output$pvalue, 4),'\n\n'))

  }

  pander(output$pairwiseca, plain.ascii = TRUE, caption = 'Pairwise Common Area')

  pander(output$pairwisepvalue, plain.ascii = TRUE, caption = 'Pairwise p-value')
}


#' @export

pairs.kTest = function(output) {

  par(mfrow = c(length(output$labels) - 1, length(output$labels) - 1))

  for(i in 1:length(output$labels)) {

    for(j in i:length(output$labels)) {

      if(i == j & i != 1 & i != length(output$labels)) {

        plot.new()

      }

      if(i != j) {

        xlabel = paste('Common area =', round(output$pairwiseca[i,j], 4))

        if(!is.null(output$pairwisepvalue)) {

          xlabel = paste(xlabel, '/ p-value =', round(output$pairwisepvalue[i,j], 4))

        }

        plot(output$densities[[i]], col = 1, main = "", xlab = xlabel, ylab = "", ylim = c(0, max(c(output$densities[[i]]$y, output$densities[[j]]$y))))

        lines(output$densities[[j]], col = 2)

        legend('topright', legend = c(output$labels[i], output$labels[j]), col = 1:2, lwd = 1, bty = 'n')

      }

    }

  }

  par(mfrow=c(1,1))
}