Article originally posted on Data Science Central. Visit Data Science Central

Below is an R code for Friedman Test that includes post-hoc tests as well in case the alternative hypothesis is rejected.

Feel free to use the code after copying and pasting it into R workspace.

friedman.test <- function(data, alpha=0.05) {
#—————————————————————————–
# Author : Okan OYMAK, MS in Operations Research at the NPS Monterey, CA, USA
# Date : 12 March 2015
#
# Data:
# The data consist of b mutually independent k-variate random variables
# (X_i1, X_i2, …, X_ik), called b blocks, i=1,2,…,b.
# The random variable is X_ij in block i and is associated with treatment j .
# The b blocks are arranged as follows.
#
# Treatments
# Block 1 2 … k
# 1 X_11 X_12 … X_1k
# 2 X_21 X_22 … X_2k
# . …. …. … ….
# . …. …. … ….
# . …. …. … ….
# b X_b1 X_b2 … X_bk
#
# Let R(Xij) be the rank, from 1 to k, assigned to Xij within block (row) i.
# That is, for block i the random variables X_i1, X_i2, …, X_ik are
# compared with each other and rank 1 is assigned to the smallest observed
# value, rank 2 to the second smallest, and so on to the rank k, which is
# assigned to the largest observation in block i.
#
# Assumptions:
# 1. The b k-variate random variables are mutually independent. (The results
# within one block do not influence the results within the other blocks)
# 2. Within each block the observations may be ranked acoording to some
# criteria of interest.
#
# Hypotheses:
# Ho: Each ranking of the random variables within a block is equally likely.
# (i.e. the treatments have identical effectes.)
# Ha: At least one of the treatments tends to yield larger observed values
# than at least one of the other treatment.
#—————————————————————————–
# LET THE FUN BEGIN!
#—————————————————————————–
if(!is.matrix(data)) {
data <- as.matrix(data)
}
tie.check <- numeric()

for (i in 1: dim(data)[1]) {
if(length(unique(data[i,])) == dim(data)[2]) {
tie.check[i] <- 1
} else {
tie.check[i] <- 0
}
}

b <- dim(data)[1]
k <- dim(data)[2]

ranks <- t(apply(data, 1, rank))
Rj <- apply(ranks, 2, sum)

if (sum(tie.check) < b ) { # means there are ties

A1 <- sum(ranks^2)
C1 <- (b*k*(k+1)^2)/4
T1 <- ((k-1)*sum((Rj-(b*(k+1))/2)^2))/(A1-C1)
T2 <- ((b-1)*T1)/(b*(k-1)-T1)
k1 <- k-1
k2 <- (b-1)*k1
cr1 <- qchisq(1-alpha, k1)
p.value.1.T1 <- 1-pchisq(T1, k1)
cr <- qf(1-alpha, k1, k2)
p.value.1 <- 1-pf(T2, k1, k2)

cat(“Ties exist.”, “n”, “n”)
cat(“Friedman Test Statistic based on the Chi-Square Distribution, T1 =”, round(T1, 4), “n”)
cat(“Critical region of size”, alpha,
“corresponds to all values of T1 greater than”, round(cr1, 4),”n”)
cat(“p-value = “, p.value.1.T1, ” df = “, k1, “n”, “n”)

cat(“Friedman Test Statistic based on the F Distribution, T2 =”, round(T2, 4), “n”)
cat(“Critical region of size”, alpha,
“corresponds to all values of T2 greater than”, round(cr, 4),”n”)

if(p.value.1 <= alpha) {
t.stat <- qt(1 – alpha/2, k2)
cat(“p-value = “, p.value.1, “<“, “alpha = “, alpha, “n”, “n”)
cat(“Therefore, reject the null hypothesis”, “n”, “n”)
cat(“INFERENCE: At least one of the treatments tends to yield larger observed values than at least one of the other treatments.”, “n”, “n”)
cat(“Multiple comparisons are as follows”, “n”)
combinations <- combn(c(1:dim(data)[2]), 2)
for(i in 1:dim(combinations)[2]) {
test.1 <- abs(Rj[combinations[1,i]]-Rj[combinations[2,i]])
test.2 <- t.stat*sqrt( (2*(b*A1-sum(Rj^2))) / ((b-1)*(k-1)) )
cat(“Treatments:”, combinations[1,i], “-“, combinations[2,i], “n”)
cat(“|Rj – Ri| =”, test.1, ” compared to”, round(test.2, 4))
if (test.1 <= test.2) {
cat(” No Difference”, “n”, “n”)
}
else {
cat(” ***Different***”, “n”, “n”)

}
}

}
else {
cat(“p-value = “, p.value.1, “>”, “alpha = “, alpha, “n”, “n”)
cat(“Therefore, do not reject the null hypothesis”, “n”, “n”)
cat(“INFERENCE: Each ranking of the random variables within a block is equally likely.”, “n”)
cat(” (i.e. the treatments have identical effects.)”, “n”)
}

} else if (sum(tie.check) == b) { # means there are no ties

A1 <- b*k*(k+1)*(2*k+1)/6
T1 <- (12/(b*k*(k+1)))*sum((Rj-b*(k+1)/2)^2)
T2 <- ((b-1)*T1)/(b*(k-1)-T1)
k1 <- k-1
k2 <- (b-1)*k1
cr1 <- qchisq(1-alpha, k1)
p.value.0.T1 <- 1-pchisq(T1, k1)
cr <- qf(1-alpha, k1, k2)
p.value.0 <- 1-pf(T2, k1, k2)

cat(“No ties exist.”, “n”, “n”)
cat(“Friedman Test Statistic based on the Chi-Square Distribution, T1 =”, round(T1, 4), “n”)
cat(“Critical region of size”, alpha,
“corresponds to all values of T1 greater than”, round(cr1, 4),”n”)
cat(“p-value = “, p.value.0.T1, ” df= “, k1, “n”, “n”)

cat(“Friedman Test Statistic based on the F Distribution, T2 =”, round(T2, 4), “n”)
cat(“Critical region of size”, alpha,
“corresponds to all values of T2 greater than”, round(cr, 4),”n”)

if(p.value.0 <= alpha) {
t.stat <- qt(1 – alpha/2, k2)
cat(“p-value = “, p.value.0, “<“, “alpha = “, alpha, “n”, “n”)
cat(“Therefore, reject the null hypothesis”, “n”, “n”)
cat(“INFERENCE: At least one of the treatments tends to yield larger observed values than at least one of the other treatments.”, “n”, “n”)
cat(“Multiple comparisons are as follows”, “n”)
combinations <- combn(c(1:dim(data)[2]), 2)
for(i in 1:dim(combinations)[2]) {
test.1 <- abs(Rj[combinations[1,i]]-Rj[combinations[2,i]])
test.2 <- t.stat*sqrt( (2*(b*A1-sum(Rj^2))) / ((b-1)*(k-1)) )
cat(“Treatments:”, combinations[1,i], “-“, combinations[2,i], “n”)
cat(“|Rj – Ri| =”, test.1, ” compared to”, round(test.2, 4))
if (test.1 <= test.2) {
cat(” No Difference”, “n”, “n”)
}
else {
cat(” ***Different***”, “n”, “n”)

}
}

}
else {
cat(“p-value = “, p.value.0, “>”, “alpha = “, alpha, “n”, “n”)
cat(“Therefore, do not reject the null hypothesis”, “n”, “n”)
cat(“INFERENCE: Each ranking of the random variables within a block is equally likely.”, “n”)
cat(” (i.e. the treatments have identical effects.)”, “n”)

}

}
}