#' 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))
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.