Nonparametric R code
Professor: Jeng-Tung, Chiang ; Student: Chih-Hui, Wang
Tuesday, January 27, 2015
Chapter2 Two-Sample Method
2.2 Permutation Tests Based on the Median and Trimmed Means
Ref: 2-26~2-28
#Table 2.1.1
New <- c(37, 49, 55, 57)
Tradition <- c(23, 31, 46)
#Combine all and find out the index for new method
all <- c(New, Tradition)
index <- seq(along=all)
indexInNew <- combn(index, 4)
#Median difference
medianD <- rep(NA, dim(indexInNew)[2])
for(i in 1:dim(indexInNew)[2]){
medianD[i] <- median(all[indexInNew[, i]]) - median(all[-indexInNew[, i]])
}
#Exact p-value
(p_value <- sum(medianD >= medianD[1])/length(medianD))## [1] 0.057142862.6 Mann-Whitney Test and a Confidence Interval
Ref: 2-64~2-65
#Table2.6.2
Granite <- c(33.63, 39.86, 69.32, 42.13, 58.36, 74.11)
Basalt <- c(26.15, 18.56, 17.55, 9.84, 28.29, 34.15)
#Function to compute pair difference
Pairdiff <- function(x, y){
n1 <- length(x);n2 <- length(y)
paird <- c()
for(i in 1:n1){
paird <- c(paird, rep(x[i], n2) - y)
}
paird
}
d <- Pairdiff(Granite, Basalt)
sort(d)## [1] -0.52 5.34 5.71 7.48 7.98 11.57 13.71 13.84 15.07 15.98 16.08
## [12] 21.30 22.31 23.57 23.79 24.21 24.58 30.02 30.07 32.21 32.29 35.17
## [23] 39.80 39.96 40.81 41.03 43.17 45.82 47.96 48.52 50.76 51.77 55.55
## [34] 56.56 59.48 64.27#Hodes-Lehmann estimator of delta
median(d)## [1] 30.045#CI for delta
deltaCI <- function(x, y, q=0.05){
#Above function
Pairdiff <- function(x, y){
n1 <- length(x);n2 <- length(y)
paird <- c()
for(i in 1:n1){
paird <- c(paird, rep(x[i], n2) - y)
}
paird
}
#Save paired difference in d
d <- sort(Pairdiff(x, y))
nlower <- qwilcox(q, length(x), length(y))
nupper <- qwilcox(q, length(x), length(y), lower.tail=FALSE)
cat("Confidence Interval:", "\n")
cat(d[nlower], "<", "delta", "<=", d[nupper + 1])
}
deltaCI(Granite, Basalt)## Confidence Interval:
## 13.84 < delta <= 47.962.8 Test for Equality of Scale Parameters and an Ominibus Test
#A package that have function to help us find the all possible permutation outcomes
library(combinat)2.8.2 Test on Deviances(RMD)
Ref: 2-81~2-83 & 2-86~2-87
#Table 2.8.1
#Test on Deviances (Ratio of the absolute Mean Deviance)
trt1 <- c(16.55, 15.36, 15.94, 16.43, 16.01)
trt2 <- c(16.05, 15.98, 16.10, 15.88, 15.91)
dev1 <- trt1 - median(trt1);dev2 <- trt2 - median(trt2)
all <- c(dev1, dev2)
index <- seq(along=all)
indexIntrt1 <- combn(index, 5)
RMD <- rep(NA, dim(indexIntrt1)[2])
for(i in 1:dim(indexIntrt1)[2]){
RMD[i] <- mean(abs(all[indexIntrt1[, i]]))/mean(abs(all[-indexIntrt1[, i]]))
}
(p_value <- sum(RMD >= RMD[1])/length(RMD))## [1] 0.06349206#Histogram
hist(RMD, main="Distribution of RMD", xlab="RMD", ylab="Frequency")
Chapter4 Paired Comparisons and Blocked Designs
4.1 Paired-Comparison Permutation Test
Ref: 4-4
#Table 4.1.1
Hour <- c(1530, 2130, 2940, 1960, 2270)
Survey <- c(1290, 2250, 2430, 1900, 2120)#Paired-Comparison Permutation Test
library(exactRankTests)#Note that this test use sum as statistic, not mean.
perm.test(Hour, Survey, paired=TRUE, alternative="greater")##
## 1-sample Permutation Test
##
## data: Hour and Survey
## T = 960, p-value = 0.09375
## alternative hypothesis: true mu is greater than 0Ref: 4-9
#Table 4.1.3
drug1 <- c(74, 55, 61, 41, 53, 74, 52, 31, 50, 58, 54, 53, 69, 60, 61, 54, 57)
drug2 <- c(63, 58, 49, 47, 50, 69, 67, 40, 44, 38, 56, 38, 47, 41, 46, 47, 44)#Exact p-value p
perm.test(drug1, drug2, paired=TRUE, alternative="two.sided")##
## 1-sample Permutation Test
##
## data: drug1 and drug2
## T = 148, p-value = 0.02502
## alternative hypothesis: true mu is not equal to 04.4 A Permutation Test for a Randomized Complete Block Design
4.4.2 Permutation F-test for Randomized Complete Block Design
Ref 4-44
#table 4.4.3
trt1 <- c(120, 208, 199, 194, 177, 195)
trt2 <- c(207, 188, 181, 164, 155, 175)
trt3 <- c(122, 137, 177, 177, 160, 138)
trt4 <- c(128, 128, 160, 142, 157, 179)
m <- rbind(trt1, trt2, trt3, trt4)#Approximate pR
RCBDperR <- function(m, pernum=5000){
nt <- dim(m)[1]
nb <- dim(m)[2]
#Save each permutation of block in blki
blk1 <- permn(m[, 1]);blk2 <- permn(m[, 2])
blk3 <- permn(m[, 3]);blk4 <- permn(m[, 4])
blk5 <- permn(m[, 5]);blk6 <- permn(m[, 6])
#Create an array to put each simulation matrix
w=array(0, dim=c(nt, nb, pernum))
w[, , 1] <- m
for(i in 2:pernum){
#Store simulate number
simN <- c()
for(j in 1:6){
simN[j] <- sample(1:factorial(nt), 1)
}
w[, , i] <- as.array(cbind(blk1[[simN[1]]], blk2[[simN[2]]],
blk3[[simN[3]]], blk4[[simN[4]]],
blk5[[simN[5]]], blk6[[simN[6]]]))
}
#Now we have the array w which its
#third dimension represent each simulation
#Here I compute the F statistics
#However, as you know, you can use SST, SSX
SST <- rep(0, pernum)
for(i in 1:pernum){
rowmean <- apply(w[, , i], 1, mean)
SST[i] <- nb*(sum((rowmean - rep(mean(rowmean), nt))^2))
}
colmean <- apply(w[, , 1], 2, mean)
SSBL <- rep(nt*(sum((colmean - rep(mean(colmean), nb))^2)), pernum)
#Each SSTO will be the same
SSTO <- rep(var(as.numeric(w[, , 1]))*(nt*nb - 1), pernum)
SSE <- SSTO - SST - SSBL
#Compute F statistic
F <- (SST/(nt - 1))/(SSE/((nt - 1)*(nb - 1)))
(pR <- mean(F >= F[1]))
}
RCBDperR(m)## [1] 0.0594#Simulate 10000 times
RCBDperR(m, pernum=10000)## [1] 0.05674.5 Friedman’s Test for a Randomized Complete Block Design
4.5.2 Adjusted for Ties
Ref 4-49~4-51
#Table 4.5.3
trt1 <- c(100, 80, 50)
trt2 <- c(100, 80, 60)
trt3 <- c(150, 80, 80)
trt4 <- c(150, 90, 90)
m <- rbind(trt1, trt2, trt3, trt4)#Exact p-value p
FriedmanTieper <- function(m){
nt <- dim(m)[1]
nb <- dim(m)[2]
#Save each permutation of block in blki
blk1 <- permn(rank(m[, 1]))
blk2 <- permn(rank(m[, 2]))
blk3 <- permn(rank(m[, 3]))
#Create an array to put each simulation matrix
w=array(0,dim=c(nt, nb, factorial(nt)^nb))
for(i in 1: factorial(nt)){
for(j in 1: factorial(nt)){
for(k in 1: factorial(nt)){
w[, , (i - 1)*(factorial(nt)^(nb - 1)) +
(j - 1)*(factorial(nt)^(nb - 2)) +
(k - 1) + 1] <-
as.array(cbind(blk1[[i]], blk2[[j]], blk3[[k]]))
}
}
}
#Now we have the array w which its
#third dimension represent each simulation
#Here I use sum(Ribar^2)
Ribarsq <- c()
for(i in 1:factorial(nt)^nb){
Ribarsq[i] <- sum(apply(w[, , i], 1, mean)^2)
}
(p <- mean(Ribarsq >= Ribarsq[1]))
}
FriedmanTieper(m)## [1] 0.04166667#Rank
m <- rbind(trt1, trt2, trt3, trt4)
(r <- apply(m, 2, rank))## [,1] [,2] [,3]
## trt1 1.5 2 1
## trt2 1.5 2 2
## trt3 3.5 2 3
## trt4 3.5 4 4Asypmtotic p-value pA
#Asypmtotic p-value pA
k <- dim(m)[1]
b <- dim(m)[2]
FM <- ((b^2)/sum(apply(r, 2, var)))*(sum((apply(r, 1, mean) - (k + 1)/2)^2))
(pA <- pchisq(FM, k - 1, lower.tail=F))## [1] 0.057558454.6 Ordered Alternative for a Randomized Complete Block Design
Page’s Test
Ref 4-62
#Table 4.4.3
trt1 <- c(120, 208, 199, 194, 177, 195)
trt2 <- c(207, 188, 181, 164, 155, 175)
trt3 <- c(122, 137, 177, 177, 160, 138)
trt4 <- c(128, 128, 160, 142, 157, 179)
m <- rbind(trt1, trt2, trt3, trt4)#Approximate p-value pR
myPagetestR <- function(m, pernum=5000){
nt <- dim(m)[1]
nb <- dim(m)[2]
r <- apply(m, 2, rank)
#Save each permutation of block in blki
blk1 <- permn(r[, 1]);blk2 <- permn(r[, 2])
blk3 <- permn(r[, 3]);blk4 <- permn(r[, 4])
blk5 <- permn(r[, 5]);blk6 <- permn(r[, 6])
#Create an array to put each simulation matrix
w=array(0, dim=c(nt, nb, pernum))
w[, , 1] <- r
for(i in 2:pernum){
#Store simulate number
simN <- c()
for(j in 1:6){
simN[j] <- sample(1:factorial(nt), 1)
}
w[, , i] <- as.array(cbind(blk1[[simN[1]]], blk2[[simN[2]]],
blk3[[simN[3]]], blk4[[simN[4]]],
blk5[[simN[5]]], blk6[[simN[6]]]))
}
#Compute sum(i*Ri)
sumRi <- rep(0, pernum)
for(i in 1:pernum){
Rankrowmean <- apply(w[, , i], 1, sum)
sumRi[i] <- sum((1:nt)*Rankrowmean)
}
(pR <- mean(sumRi <= sumRi[1]))
}
myPagetestR(m)## [1] 0.013#Asypmtotic p-value pA
myPagetestA <- function(m){
k <- dim(m)[1]
b <- dim(m)[2]
r <- apply(m, 2, rank)
PGobs <- sum((1:k)*apply(r, 1, sum))
EPG <- (b*k*(k + 1)^2)/4
varPG <- (b*(k - 1)*(k^2)*(k + 1)^2)/144
(pA <- pnorm((PGobs - EPG)/sqrt(varPG), lower.tail=T))
}
myPagetestA(m)## [1] 0.01182581