.statisticsHolmsAdjustemnt <- function(p){
if (all(nna <- !is.na(p)))
nna <- TRUE
p2 <- p[nna]
n <- length(p2)
i <- n:1
o <- order(p2)
ro <- order(o)
p[nna] <- pmin(1, cummax(i* p2[o]))[ro]
p
}
.statisticsShafferAdjustemnt <- function(p){
n <- length(p)
o <- order(p)
ro <- order(o)
# n = k * (k-1) / 2 --> solve the quadratic equation
k <- (1+sqrt(1+8*n))/2
t <- .shafferS(k)
lT <- length(t)
extT <- c(t[lT])
# At least we are comparing 3 methods, so at least lT will be 3 ({0,1,3})
for(i in lT:3){
extT <- c(extT, rep(t[i-1],t[i]-t[i-1]))
}
pmin(1, cummax(n:1 * p[o]))[ro]
}
.shafferS <- function(k){
prev <- list()
for(i in 0:k){
tmp <- .shafferSIterative(i,prev)
#The plus 1 is because indexes of lists (must be >0)
prev[[i+1]] <- tmp
}
tmp
}
.shafferSIterative <- function(k, previous){
if(k<=1){
mySet <- c(0)
}
else{
mySet <- c()
for(j in 1:k){
#The plus 1 is because indexes of lists (must be >0)
tmp <- previous[[k-j+1]]
tmp <- tmp + choose(j,2)
mySet <- union(mySet,tmp)
}
}
mySet <- mySet[order(mySet)]
mySet
}
#' @importFrom stats psignrank
.statisticsWilcoxon <- function(x,y, alternative = "less"){
# Performs the paired wilcoxon signed rank test. Ties are allowed.
#
# Args:
# x: a numeric array (results for the first method).
# y: a numeric array (results for the second method).
#
# Returns:
# a statisticsWilcoxon object to be used in the toolbox
# PARAMETER VALIDATION:
#Check that x is a numeric array
if( !is.numeric(x) )
stop("The parameter x must be numeric.")
#Check that y is a numeric array
if( !is.numeric(y) )
stop("The parameter y must be numeric.")
#Check that x and y have the same length
if( length(x) != length(y) )
stop("The length of x and y must be the same.")
#Check that alternative is equal to "less" or "greater"
if( alternative != "less" && alternative != "greater")
stop("The alternative parameter must be \"less\" or \"greater\"")
N <- length(x)
d <- x-y
r <- rank(abs(d),ties.method = "average")
#R+
rP <- sum( r[x<y] )
#R-
rM <- sum( r[x>y] )
#In case of ties
rEq = r[x==y]
#If it is not odd, we discard one case
if(length(r[x==y]) %% 2==1){
rEq <- rEq[-1]
N <- N - 1
}
rEq <- sum(rEq)
rP <- rP+rEq/2.0
rM <- rM+rEq/2.0
if(alternative=="less")
t = min(rP,rM)
else
t = max(rP,rM)
#If there is only one data in x and y, an it is the same it would be removed and N would be equal to 0.
#In that case, psignrank(0,0) would rise a warning message.
#Instead of that, we undo the removal of that case (only 1 draw). Therefore,
#we set the pvalue equal to 1 and the statistic to 0.5
if(N==0){
t <- 0.5
pvalue <- 1
}
else
pvalue <- psignrank(t,N)
results <- list(
"statistic" = t,
"pvalue" = pvalue
)
class(results) <- c("statisticsWilcoxon")
results
}
.statisticsFriedman <- function(data, rankOrder){
# Performs the friedman test.
#
# Args:
# data: a data.frame containing the results for each method (rows) for each problem (columns)
# rankOrder: tell us if we want to maximize ("max") or minimize ("min") the output variable.
#
# Returns:
# a .statisticsFriedman object to be used in the toolbox
# PARAMETER VALIDATION:
#Check that variable x is a data.frame
if( !is.data.frame(data) || all(dim(data)==dim(data.frame())) )
stop("The parameter data must be a non-empty data.frame")
#Check that rankOrder is one of the two possible values
if(rankOrder!="max" && rankOrder!="min")
stop("The parameter rankOrder must be \"max\" or \"min\".")
#Number of instantiations
k <- dim(data)[1]
#Number of problems
n <- dim(data)[2]
if (rankOrder == "min" )
rankModifier <- 1
else
rankModifier <- -1
ranks <- sapply(data*rankModifier, FUN=rank)
#If k is equal to 1, sapply simplifies it to an array. In that case, we use the t function
#which traspose x. If it is a numeric array, it just makes what we want to get (one row matrix)
if(!is.matrix(ranks))
ranks <- t(ranks)
rownames(ranks) <- rownames(data)
friedmanRanks <- rowMeans(ranks)
#We obtain the chi-square statistic. At least k is equal to 1, so diviion by zero is avoided.
chisqStatistic <- 12 * n / ( k * (k+1) ) * ( sum(friedmanRanks^2) - ( k * (k+1)^2 ) / 4 )
friedmanDegreeOfFreedom <- k-1
friedmanPValue <- pchisq(chisqStatistic, friedmanDegreeOfFreedom, lower.tail=FALSE)
results <- list(
"statistic" = chisqStatistic,
"degreesOfFreedom" = friedmanDegreeOfFreedom,
"pvalue" = friedmanPValue,
"n" = n,
"k" = k,
"ranks" = ranks
)
class(results) <- c(".statisticsFriedman")
results
}
#' @importFrom stats pf
.statisticsImanDavenport <- function(f){
# Performs the iman davenport test.
#
# Args:
# f: a .statisticsFriedman object
#
# Returns:
# a .statisticsImanDavenport object to be used in the toolbox
#Number of instantiations
k <- f$k
#Number of problems
n <- f$n
if( k<2 || n < 2)
stop("To compute the Iman Davenport test there must be at least 3 methods and at least 3 problems")
chisqStatistic <- f$statistic
imanDavenStatistic <- ( (n-1) * chisqStatistic ) / ( n * (k-1) - chisqStatistic )
imanDavenDegreeOfFreedom_1 <- k-1
imanDavenDegreeOfFreedom_2 <- (k-1) * (n-1)
imanDavenPValue <- pf(imanDavenStatistic, imanDavenDegreeOfFreedom_1, imanDavenDegreeOfFreedom_2, lower.tail=FALSE)
results <- list(
"statistic" = imanDavenStatistic,
"degreesOfFreedom1" = imanDavenDegreeOfFreedom_1,
"degreesOfFreedom2" = imanDavenDegreeOfFreedom_2,
"pvalue" = imanDavenPValue,
"n" = n,
"k" = k,
"ranks" = f$ranks
)
class(results) <- c(".statisticsImanDavenport")
results
}
#' @importFrom stats pf
.statisticsGeneralControl <- function(f){
# Performs the general control test.
#
# Args:
# f: a testFriedman object
#
# Returns:
# a .statisticsGeneralControl object to be used in the toolbox
# PARAMETER VALIDATION:
#Check that variable f is a testFriedman object
if( !is.testFriedman(f) )
stop("The parameter f must be a testFriedman object")
friedmanRanks <- rowMeans(f$ranks)
#Number of instantiations
k <- dim(f$ranks)[1]
#Number of problems
n <- dim(f$ranks)[2]
controlIdx <- which.min(friedmanRanks)
criticalDiff <- sqrt( (k*(k+1))/(6*n) )
# Generate matrix with statisctics and p-values for post-hoc test:
p <- c()
z <- c()
for (i in 1:k){
z[i] <- abs( ( friedmanRanks[i] - friedmanRanks[controlIdx] ) / criticalDiff )
p[i] <- 2 * pnorm( z[i], lower.tail=FALSE)
}
#The p-value in the position controlIdx is equal to pnorm( 0, lower.tail=FALSE ) --> It must be NA
p[controlIdx] <- NA
results <- list(
"statistic" = z,
"criticalDifference" = criticalDiff,
"controlIndex" = controlIdx,
"pvalue" = p
)
class(results) <- c(".statisticsGeneralControl")
results
}
#' @importFrom stats pf
.statisticsGeneralPairwise <- function(f){
# Performs the general pairwise test.
#
# Args:
# f: a testFriedman object
#
# Returns:
# a .statisticsGeneralPairwise object to be used in the toolbox
# PARAMETER VALIDATION:
#Check that variable f is a testFriedman object
if( !is.testFriedman(f) )
stop("The parameter f must be a testFriedman object")
friedmanRanks <- rowMeans(f$ranks)
# Order data according to ranks. It is necessary because we apply the test comparisson only in one direction (greater versus worse)
correct_data_order <- order(friedmanRanks)
friedmanRanks <- friedmanRanks[correct_data_order]
algorithmNames <- names(friedmanRanks)
#Number of instantiations
k <- dim(f$ranks)[1]
#Number of problems
n <- dim(f$ranks)[2]
criticalDiff <- sqrt( (k*(k+1))/(6*n) )
met1 <- c()
met2 <- c()
p <- c()
z <- c()
counter <- 1
for (i in 1:(k-1))
{
for (j in (i+1):k)
{
met1[counter] <- algorithmNames[i]
met2[counter] <- algorithmNames[j]
z[counter] <- abs( ( friedmanRanks[j] - friedmanRanks[i] ) / criticalDiff )
p[counter] <- 2 * pnorm( z[counter], lower.tail=FALSE)
counter <- counter + 1
}
}
names <- data.frame( method1 = met1, method2 = met2)
results <- list(
"statistic" = z,
"criticalDifference" = criticalDiff,
"names" = names,
"pvalue" = p
)
class(results) <- c(".statisticsGeneralPairwise")
results
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.