Nothing
# version 1.1.5 October 28, 2017 by alexis.dinno@pdx.edu
# perform Conover-Iman test of multiple comparisons using rank sums
p.adjustment.methods <- c("none","bonferroni","sidak","holm","hs","hochberg","bh","by")
conover.test <- function(x=NA, g=NA, method=p.adjustment.methods, kw=TRUE, label=TRUE, wrap=FALSE, table=TRUE, list=FALSE, rmc=FALSE, alpha=0.05, altp=FALSE) {
# FUNCTIONS
# kwallis.test: a custom Kruskal Wallis test function to support conover.test
# Note: does not currently play nicely with missing values.
kwallis.test <- function(x=NA,g=NA,pass=0) {
# FUNCTIONS
# tiedranks: enumerates tied values from a vector of ranks
# ranks: a vector of rank values
# returns: ties, a vector of values that are tied
tiedranks <- function(ranks) {
ranks <- sort(ranks)
#enumerate tied values
ties <- c()
for (i in 2:length(ranks)) {
if (ranks[i-1] == ranks[i]) {
if (length(ties) > 0) {
if (ranks[i-1] != tail(ties,n=1)) {
ties <- c(ties,ranks[i-1])
}
}
else {
ties <- c(ranks[i-1])
}
}
}
return(ties)
}
#set up data by lists
if (is.list(x)) {
N <- 0
for (i in 1:length(x)) {
N <- N + length(x[[i]])
}
Data <- matrix(NA,N,4)
for (i in 1:N) {
Data[i,1] <- i
}
obs <- c()
group <- c()
for (i in 1:length(x)) {
obs <- c(obs,x[[i]])
group <- c(group,rep(i,length(x[[i]])))
}
Data[,2] <- obs
if (length(g) > 1) {
Data[,3] <- g
}
else {
Data[,3] <- group
}
Data[,4] <- rank(Data[,2], ties.method="average", na.last=NA)
}
#set up data by groups
if (!is.list(x)) {
N <- length(x)
Data <- matrix(NA,length(x),4)
Data[,1] <- 1:length(x)
Data[,2] <- x
Data[,3] <- g
Data[,4] <- rank(Data[,2], ties.method="average", na.last=NA)
}
k <- length(unique(Data[,3]))
#calculate ties adjustment
ranks <- Data[,4]
ranks <- ranks[!is.na(ranks)]
ties <- tiedranks(ranks)
r <- length(ties)
tiesadj <- 0
if (r > 0) {
for (s in 1:r) {
tau <- sum(ranks==ties[s])
tiesadj <- tiesadj + (tau^{3} - tau)
}
}
tiesadj <- 1-(tiesadj/((N^3) - N))
#calculate H
ranksum <- 0
for (i in 1:k) {
ranksum <- ranksum + ((sum(Data[,4][Data[,3]==i]))^2)/(sum(Data[,3]==i))
}
H <- ((12/(N*(N+1)))*ranksum - 3*(N+1))/tiesadj
df <- k-1
p <- pchisq(H,k-1,lower.tail=FALSE)
#present output
output <- paste("Kruskal-Wallis chi-squared = ",round(H,digits=4),", df = ",df,", p-value = ",round(p,digits=2),"\n\n" ,sep="")
invisible(list(output=output,H=H,df=df,p=p,N=N,Data=Data,data.name=data))
}
# all.integers robustly tests whether all elements of a vector are integers
all.integers <- function(x) {
for (i in length(x)) {
if (is.na(x[i]) | is.list(x[i]) | length(x[i]) > 1 | !is.numeric(x[i])) {
return(FALSE)
}
if (x[i]%%1!=0) {
return(FALSE)
}
}
return(TRUE)
}
# tpad: returns the pad string, n times, defaulting to spaces
# n: a number of replications; pad: string to replicate
tpad <- function(n=1, pad=" ") {
if (n == 0) {
return("")
}
return( paste0(rep(x=pad,times=n),collapse="") )
}
# tformat: formats t values for display in table
# t: a real t-value
# returns: a formatted string
tformat <- function(t) {
if (t < 0) {
sign <- "-"
}
else {
sign <- " "
}
leftspaces <- max(2,floor(log10(abs(t)))+2)
leftdigits <- floor(abs(t))
rightspaces <- 8 - leftspaces
rightdigits <- substr(paste0(abs(t) - floor(abs(t)),"00000000"),3,rightspaces+2)
return(paste0(sign,leftdigits,".",rightdigits))
}
# pformat: formats p values for display in table
# p: a real p-value
# returns: a formatted string
pformat <- function(p) {
if (p == 1) {
return("1.0000")
}
else {
rightspaces <- 4
rightdigits <- substr(paste0(sprintf("%1.4f",p),"000000000"),3,rightspaces+2)
return(paste0("0.",rightdigits))
}
}
# centertext: centers a string within a specific width
centertext <- function(text,width=80,lmargin=2,rmargin=2) {
textwidth <- nchar(text)
if (textwidth <= width-lmargin-rmargin) {
text <- substr(text,1,width-lmargin-rmargin)
}
buff <- (width-lmargin-rmargin-textwidth)
if (buff%%2 == 0) {
return(paste(paste(rep(" ",buff/2),collapse=""),text,paste(rep(" ",buff/2),collapse=""),"\n",collapse=""))
}
else {
return(paste(paste(rep(" ",1+(buff-1)/2),collapse=""),text,paste(rep(" ",(buff-1)/2),collapse=""),"\n",collapse=""))
}
}
# conovertestheader: displays the Conover-Iman test table headers.
conovertestheader <- function(groupvar,colstart,colstop,rmc) {
if (rmc==FALSE) {
cat("Col Mean-|\nRow Mean |")
}
else {
cat("Row Mean-|\nCol Mean |")
}
groupvalues <- levels(factor(groupvar))
for (col in colstart:colstop) {
vallab <- substr(groupvalues[col],1,8)
pad <- 8-nchar(vallab)
colhead <- paste0(tpad(n=pad),vallab)
cat(paste0(" ",substr(colhead,1,8)),sep="")
}
cat("\n")
separatorlength = 10 + 11*(colstop-colstart)+1
cat(tpad(n=9,pad="-"),"+",tpad(n=separatorlength,pad="-"),"\n",sep="")
}
# conovertestttable displays the Conover-Iman test t values
conovertestttable <- function(groupvar, index, Tvalues, colstart, colstop) {
groupvalues <- levels(factor(groupvar))
# Row headers
vallab <- substr(groupvalues[index],1,8)
pad <- 8-nchar(vallab)
rowhead <- paste0(tpad(n=pad),vallab)
cat(paste0(rowhead," |"))
# Table t entries
for (i in colstart:colstop) {
t <- Tvalues[(((index-2)*(index-1))/2) + i]
cat(paste0(" ",tformat(t)))
}
cat("\n")
}
# conovertestptable: displays the Conover-Iman test p values
conovertestptable <- function(index, P, colstart, colstop, Reject,last) {
# Blank row header
cat(" |")
# Table p entries
for (i in colstart:colstop) {
p <- P[(((index-2)*(index-1))/2) + i]
if ( Reject[(((index-2)*(index-1))/2) + i] == 0) {
cat(paste0(" ",pformat(p)),sep="")
}
else {
cat(paste0(" ",pformat(p)),"*",sep="")
}
}
cat("\n")
# Close out with another blank row header
if (last == 0) {
cat(" |\n")
}
}
# alphabetize.factor: alphabetizes a factor
# this is very quick and dirty with no checking.
alphabetize.factor <- function(x) {
return(factor(sort(c(as.character(x)))))
}
# VALIDATIONS & PREPARATIONS
# names for output
xname <- paste(if(is.name(substitute(x))) {deparse(substitute(x))} else {"x"})
gname <- paste(if(is.name(substitute(g))) {deparse(substitute(g))} else {"group"})
xaslist <- is.list(x)
tempg <- 1:length(x)
# casewise deletion of missing data
if (length(g) > 1) {
if (label==TRUE) {
if (is.factor(g)) {
glevels <- levels(g)
}
else {
g <- factor(g)
glevels <- levels(g)
}
}
else {
glevels <- names(table(addNA(g, ifany = TRUE)))
}
Data <- data.frame(x,g)[order(levels(g)[c(g)]),]
Data <- Data[!is.na(unlist(Data$x)),]
Data <- Data[!is.na(Data$g),]
x <- Data[,1]
g <- alphabetize.factor(Data$g)
}
else {
g <- c()
for (i in 1:length(x)) {
for (j in 1:length(x[[i]])) {
g <- c(g,i)
}
}
x <- as.numeric(unlist(x)[!is.na(unlist(x))])
}
# validate method
if (length(method) > 1) {
method <- "none"
}
if (tolower(method) != "none" & tolower(method) != "bonferroni" & tolower(method) != "sidak" & tolower(method) != "holm" & tolower(method) != "hs" & tolower(method) != "hochberg" & tolower(method) != "bh" & tolower(method) != "by") {
stop("method must be one of: none, bonferroni, sidak, hs, bh, or by")
}
if (tolower(method)=="none") {
Name <- "(No adjustment)"
}
if (tolower(method)=="bonferroni") {
Name <- "(Bonferroni)"
}
if (tolower(method)=="sidak") {
Name <- "(\u0160id\u00E1k)"
}
if (tolower(method)=="holm") {
Name <- "(Holm)"
}
if (tolower(method)=="hs") {
Name <- "(Holm-\u0160id\u00E1k)"
}
if (tolower(method)=="hochberg") {
Name <- "(Hochberg)"
}
if (tolower(method)=="bh") {
Name <- "(Benjamini-Hochberg)"
}
if (tolower(method)=="by") {
Name <- "(Benjamini-Yekuteili)"
}
# validate that x is longer than 1
if (length(unlist(x))==1) {
stop("too few observations in x.")
}
# validate that x is a numeric vector of data values, or a list of numeric
# data vectors, and is not NA
if (TRUE %in% is.na(unlist(x)) | !is.vector(unlist(x), mode="numeric") ) {
stop("x must contain a numeric vector of data values, or a list of numeric data vectors.")
}
# validate that g is not missing if x is a vector
if (!xaslist & TRUE %in% is.na(g) ) {
stop("when specifying x as a vector, you must include g.")
}
# validate that g is not NA
if (length(g) > 1 & TRUE %in% is.na(g)) {
stop("g must have no missing values.")
}
# validate that g is factor or vector.
if (length(g) > 1 & ( !is.factor(g) & !is.vector(g) ) ) {
stop("g must be a factor, character vector, or integer vector.")
}
# validate that g is a vector of mode = character or mode = integer.
if (length(g) > 1 & is.vector(g) ) {
if ( !is.vector(g, mode="character") & !all.integers(g) ) {
stop("g must be a factor, character vector, or integer vector.")
}
}
# CALCULATIONS
out <- NULL
if (xaslist & length(g)==1) {
kwallis.test(x, 1:length(x), pass=1) -> out
}
if (length(g)>1 & xaslist) {
kwallis.test(x, g, pass=1) -> out
}
if (length(g)>1 & !xaslist) {
kwallis.test(x, g, pass=0) -> out
}
if (kw==TRUE) {
cat(" Kruskal-Wallis rank sum test\n\ndata: ",xname," and ",gname,"\n",sep="")
cat(out$output)
}
chi2 <- out$H
k <- out$df + 1
N <- out$N
Data <- out$Data
m <- k*(k-1)/2
Tvalues <- rep(0,m)
P <- rep(0,m)
nu <- N - k
# calculate pooled variance to be used in generating t statistics
ranks <- Data[,4]
varA <- 1/(N-1)
varB <- 0
for (i in 1:N) {
varB <- varB + ranks[i]^2
}
varB <- varB - (N*(((N+1)^2)/4))
S2 <- varA*varB
# Calculate t test statistics
for (i in 2:k) {
for (j in 1:(i-1)) {
ni <- sum(Data[,3]==i)
nj <- sum(Data[,3]==j)
meanranki <- mean(Data[,4][Data[,3]==i])
meanrankj <- mean(Data[,4][Data[,3]==j])
if (rmc == FALSE) {
t <- (meanrankj-meanranki)/(sqrt(S2*((N-1-chi2)/(N-k)))*sqrt(1/nj + 1/ni))
}
if (rmc == TRUE) {
t <- (meanranki-meanrankj)/(sqrt(S2*((N-1-chi2)/(N-k)))*sqrt(1/nj + 1/ni))
}
index <- ((i-2)*(i-1)/2) + j
Tvalues[index] <- t
}
}
# Calculate p-values for t statistics, and adjust as needed
# If p = P(T >= |t|)
if (altp==FALSE) {
P <- pt(q=abs(Tvalues), df=nu, lower.tail=FALSE)
}
# Otherwise, if p = P(|T| >= |t|)
else {
P <- 2*pt(q=abs(Tvalues), df=nu, lower.tail=FALSE)
}
#calculate adjusted p-values based on method argument
Reject <- rep(0,m)
# No adjustment for multiple comparisons
if (tolower(c(method))=="none") {
P.adjust <- P
# If p = P(T >= |t|)
if (altp==FALSE) {
Reject <- P.adjust <= alpha/2
}
# Otherwise, if p= P(|T| >= |t|)
else {
Reject <- P.adjust <= alpha
}
}
# Control FWER using (Dunn's) Bonferroni
if (tolower(c(method))=="bonferroni") {
P.adjust <- pmin(1,P*m)
# If p = P(T >= |t|)
if (altp==FALSE) {
Reject <- P.adjust <= alpha/2
}
# Otherwise, if p= P(|T| >= |t|)
else {
Reject <- P.adjust <= alpha
}
}
# Control FWER using Šidák
if (tolower(c(method))=="sidak") {
P.adjust <- pmin(1,1 - (1-P)^m)
# If p = P(T >= |t|)
if (altp==FALSE) {
Reject <- P.adjust <= alpha/2
}
# Otherwise, if p= P(|T| >= |t|)
else {
Reject <- P.adjust <= alpha
}
}
# Control FWER using Holm(-Bonferroni)
if (tolower(c(method))=="holm") {
Psort <- matrix(c(P,1:m,rep(0,m)),3,m,byrow=TRUE)
Psort <- Psort[,order(Psort[1,])]
for (i in 1:m) {
adjust <- m+1-i
Psort[1,i] <- pmin(1,Psort[1,i]*adjust)
# If p = P(T >= |t|)
if (altp==FALSE) {
if (i==1) {
Psort[3,i] <- Psort[1,i] <= alpha/2
}
else {
Psort[3,i] <- ((Psort[1,i] <= alpha/2) & Psort[3,i-1] != 0)
}
}
# Otherwise, if p = P(|T| >= |t|)
else {
if (i==1) {
Psort[3,i] <- Psort[1,i] <= alpha
}
else {
Psort[3,i] <- ((Psort[1,i] <= alpha) & Psort[3,i-1] != 0)
}
}
}
Psort <- Psort[,order(Psort[2,])]
P.adjust <- Psort[1,]
Reject <- Psort[3,]
}
# Control FWER using Holm-Šidák
if (tolower(c(method))=="hs") {
Psort <- matrix(c(P,1:m,rep(0,m)),3,m,byrow=TRUE)
Psort <- Psort[,order(Psort[1,])]
for (i in 1:m) {
adjust <- m+1-i
Psort[1,i] <- pmin(1,(1 - ((1 - Psort[1,i])^adjust)))
# If p = P(T >= |t|)
if (altp==FALSE) {
if (i==1) {
Psort[3,i] <- Psort[1,i] <= alpha/2
}
else {
Psort[3,i] <- ((Psort[1,i] <= alpha/2) & Psort[3,i-1] != 0)
}
}
# Otherwise, if p = P(|T| >= |t|)
else {
if (i==1) {
Psort[3,i] <- Psort[1,i] <= alpha
}
else {
Psort[3,i] <- ((Psort[1,i] <= alpha) & Psort[3,i-1] != 0)
}
}
}
Psort <- Psort[,order(Psort[2,])]
P.adjust <- Psort[1,]
Reject <- Psort[3,]
}
# Control FWER using Hochberg
if (tolower(c(method))=="hochberg") {
Psort <- matrix(c(P,1:m,rep(0,m)),3,m,byrow=TRUE)
Psort <- Psort[,order(Psort[1,], decreasing=TRUE)]
for (i in 1:m) {
adjust <- i
Psort[1,i] <- min(1,Psort[1,i]*adjust)
# If p = P(T >= |t|)
if (altp==FALSE) {
if (i==1) {
Psort[3,i] <- Psort[1,i] <= alpha/2
}
else {
Psort[3,i] <- ((Psort[1,i] <= alpha/2) | Psort[3,i-1] == 1)
}
}
# Otherwise, if p = P(|T| >= |t|)
else {
if (i==1) {
Psort[3,i] <- Psort[1,i] <= alpha
}
else {
Psort[3,i] <- ((Psort[1,i] <= alpha) | Psort[3,i-1] == 1)
}
}
}
Psort <- Psort[,order(Psort[2,])]
P.adjust <- Psort[1,]
Reject <- Psort[3,]
}
# Control FDR using Benjamini-Hochberg
if (tolower(c(method))=="bh") {
Psort <- matrix(c(P,1:m,rep(0,m)),3,m,byrow=TRUE)
Psort <- Psort[,order(Psort[1,], decreasing=TRUE)]
for (i in 1:m) {
adjust <- (m/(m+1-i))
Psort[1,i] <- min(1,Psort[1,i]*adjust)
# If p = P(T >= |t|)
if (altp==FALSE) {
if (i==1) {
Psort[3,i] <- Psort[1,i] <= alpha/2
}
else {
Psort[3,i] <- ((Psort[1,i] <= alpha/2) | Psort[3,i-1] == 1)
}
}
# Otherwise, if p = P(|T| >= |t|)
else {
if (i==1) {
Psort[3,i] <- Psort[1,i] <= alpha
}
else {
Psort[3,i] <- ((Psort[1,i] <= alpha) | Psort[3,i-1] == 1)
}
}
}
Psort <- Psort[,order(Psort[2,])]
P.adjust <- Psort[1,]
Reject <- Psort[3,]
}
# Control FDR using Benjamini-Yekuteili
if (tolower(c(method))=="by") {
Psort <- matrix(c(P,1:m,rep(0,m)),3,m,byrow=TRUE)
Psort <- Psort[,order(Psort[1,], decreasing=TRUE)]
for (i in 1:m) {
adjust <- (m/(m+1-i))*sum(1/(1:m))
Psort[1,i] <- min(1,Psort[1,i]*adjust)
# If p = P(T >= |t|)
if (altp==FALSE) {
if (i==1) {
Psort[3,i] <- Psort[1,i] <= alpha/2 # reverse sorted, so m-i+1, rather than i
}
else {
Psort[3,i] <- ((Psort[1,i] <= alpha/2) | Psort[3,i-1] == 1)
}
}
# Otherwise, if p = P(|T| >= |t|)
else {
if (i==1) {
Psort[3,i] <- Psort[1,i] <= alpha
}
else {
Psort[3,i] <- ((Psort[1,i] <= alpha) | Psort[3,i-1] == 1)
}
}
}
Psort <- Psort[,order(Psort[2,])]
P.adjust <- Psort[1,]
Reject <- Psort[3,]
}
# OUTPUT
cat("\n")
if (table==TRUE | list==TRUE) {
if ((TRUE %in% is.na(g))) {
title <- paste("Comparison of ",xname," across ",k," groups",sep="")
}
else {
title <- paste("Comparison of ",xname," by ",gname,sep="")
}
cat(centertext(title))
cat(centertext(Name))
}
if (table==TRUE) {
# Need to determine how many tables (reps) to output
reps <- floor((k-1)/6)
laststart <- (reps*6) + 1
kminusone <- k - 1
if (label==FALSE) {
g <- as.numeric(g)
}
if (length(g)==1) {
g <- 1:k
}
# Replication loop for >7 groups, no wrap
if (wrap==FALSE) {
if (k > 7) {
for (rep in 1:reps) {
colstart <- (6*rep)-5
colstop <- 6*rep
conovertestheader(g,colstart,colstop,rmc)
# Table body
for (i in (colstart+1):k) {
colstop <- min(i-1,6*rep)
conovertestttable(g,i,Tvalues,colstart,colstop)
if (i < k) {
conovertestptable(i,P.adjust,colstart,colstop,Reject,0)
}
else {
conovertestptable(i,P.adjust,colstart,colstop,Reject,1)
}
}
}
# End of table
if (laststart < k) {
conovertestheader(g,laststart,kminusone,rmc)
# Table body
for (i in (laststart+1):k) {
colstop <- i-1
conovertestttable(g,i,Tvalues,laststart, colstop)
if (i < k) {
conovertestptable(i,P.adjust,laststart, colstop,Reject,0)
}
else {
conovertestptable(i,P.adjust,laststart, colstop,Reject,1)
}
}
}
}
# Replication loop for <=7 groups
if (k <= 7) {
conovertestheader(g,1,kminusone,rmc)
# Table body
for (i in 2:k) {
colstop <- i-1
conovertestttable(g,i,Tvalues,1,colstop)
if (i < k) {
conovertestptable(i,P.adjust,1,colstop,Reject,0)
}
else {
conovertestptable(i,P.adjust,1,colstop,Reject,1)
}
}
}
}
# Replication loop for >7 groups, with wrap
if (wrap==TRUE) {
conovertestheader(g,1,kminusone,rmc)
# Table body
for (i in 2:k) {
colstop <- i-1
conovertestttable(g,i,Tvalues,1,colstop)
if (i < k) {
conovertestptable(i,P.adjust,1,colstop,Reject,0)
}
else {
conovertestptable(i,P.adjust,1,colstop,Reject,1)
}
}
}
cat("\n")
}
# Output pairwise comparisons as a list if requested.
if (list==TRUE) {
groupvalues <- levels(factor(g))
# get the lengths of each group name (whether explicitly labeled or not)
lengths <- sort(unlist(lapply(groupvalues,nchar)))
# stringlength will be the sum of the two largest values in lengths plus 6
stringlength <- lengths[length(lengths)] + lengths[length(lengths)-1] + 6
# Output list header
if (tolower(method)=="none") {
cat("\nList of pairwise comparisons: t statistic (p-value)\n")
}
else {
cat("\nList of pairwise comparisons: t statistic (adjusted p-value)\n")
}
cat(paste0(rep("-",stringlength+18+max(Reject)),collapse=""))
cat("\n")
index <- 0
for (i in 2:k) {
for (j in 1:(i-1)) {
index <- index + 1
buffer <- max(stringlength - (nchar(groupvalues[i]) + nchar(groupvalues[j]) + 4) - 2,0)
if ( Reject[index] == 0) {
pformatted <- paste0("(",pformat(P.adjust[index]),")",sep="")
}
else {
pformatted <- paste0("(",pformat(P.adjust[index]),")","*",sep="")
}
if (rmc==FALSE) {
cat(groupvalues[j]," - ",groupvalues[i],paste(rep(" ",buffer))," : ",tformat(Tvalues[index]), " ",pformatted,"\n", sep="")
}
if (rmc==TRUE) {
cat(groupvalues[i]," - ",groupvalues[j],paste(rep(" ",buffer))," : ",tformat(Tvalues[index]), " ",pformatted,"\n", sep="")
}
}
}
cat("\n")
}
cat("alpha = ",alpha,"\n",sep="")
# If p = P(T <= |t|)
if (altp==FALSE) {
cat("Reject Ho if p <= alpha/2\n")
}
# Otherwise, if p = P(|T| <= |t|)
else {
cat("Reject Ho if p <= alpha\n")
}
# Create comparisons variable for returned values (whether the list option
# is TRUE or FALSE
comparisons <- rep(NA,(k*(k-1)/2))
groupvalues <- levels(factor(g))
index <- 0
for (i in 2:k) {
for (j in 1:(i-1)) {
index <- index + 1
if (rmc==FALSE) {
comparisons[index] <- paste0(groupvalues[j]," - ",groupvalues[i])
}
if (rmc==TRUE) {
comparisons[index] <- paste0(groupvalues[i]," - ",groupvalues[j])
}
}
}
# If p = P(T <= |t|)
if (altp==FALSE) {
invisible(list(chi2=chi2, T=Tvalues, P=P, P.adjusted=P.adjust, comparisons=comparisons))
}
# Otherwise, if p = P(|T| <= |t|)
else {
invisible(list(chi2=chi2, T=Tvalues, altP=P, altP.adjusted=P.adjust, comparisons=comparisons))
}
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.