###########################
## H0: Liver enzyme data ##
###########################

require(coin)

## Base R
boxplot(asat~group, asat, col='gray90')
t.test(asat~group, asat)
wilcox.test(asat~group, asat)

## Evaluation of ASL
wilcox_test(asat~group, asat)
wilcox_test(asat~group, asat, distribution="exact")
wilcox_test(asat~group, asat, distribution="approximate")
wilcox_test(asat~group, asat, distribution=approximate(B=10000))

## Choice of scores
wilcox_test(asat~group, asat)
normal_test(asat~group, asat)
median_test(asat~group, asat)

## Same as
a.wilcox <- rank(asat$asat)
independence_test(a.wilcox~group, asat)
n <- nrow(asat)
a.vdw <- qnorm(rank(asat$asat)/(n+1))
independence_test(a.vdw~group, asat)
a.median <- rank(asat$asat) > (n+1)/2
independence_test(a.median~group, asat)
pairs(data.frame(a.wilcox, a.vdw, a.median), pch=19)

#########################################
## H1: Cystic fibrosis crossover study ##
#########################################

## Base R
cf <- read.delim("http://web.as.uky.edu/statistics/users/pbreheny/621/data/cysticfibrosis.txt")
with(cf, hist(Placebo-Drug, col="gray70", border=FALSE))
with(cf, t.test(Placebo-Drug))
with(cf, wilcox.test(Drug, Placebo, paired=TRUE))

## coin
wilcoxsign_test(Drug~Placebo, cf, distribution="exact")
 
## Same as:
n <- nrow(cf)
ID <- factor(rep(1:n, rep(2, n)))
s <- factor(rep(c("-1","+1"), n))
calcA <- function(scores, signs) {
  pos <- scores * (signs > 0)
  neg <- scores * (signs < 0)
  drop(as.vector(t(cbind(pos, neg))))
}
diffs <- with(cf, Placebo-Drug)
r <- rank(abs(diffs))
a.wilcox <- calcA(r, sign(diffs))
independence_test(a.wilcox~s|ID, cf, distribution="exact")

## Choice of scores
a.vdw <- calcA(qnorm(.5+r/(2*(n+1))), sign(diffs))
independence_test(a.vdw~s|ID, cf)
a.sign <- calcA(1, sign(diffs))
independence_test(a.sign~s|ID, cf, distribution="exact")
binom.test(11,14) ## Same as above
pairs(data.frame(a.wilcox, a.vdw, a.sign), pch=19)

########################
## H2: US Arrest data ##
########################

## Base R
pairs(USArrests, pch=19)
cor.test(~ Murder + UrbanPop, data = USArrests)
cor.test(~ Murder + UrbanPop, data = USArrests, method="spearman")
cor.test(~ Assault + UrbanPop, data = USArrests, method="spearman")
cor.test(~ Rape + UrbanPop, data = USArrests, method="spearman")

## Coin
spearman_test(Assault ~ UrbanPop, data = USArrests)

## Same as
independence_test(rank(Assault)~rank(UrbanPop), USArrests)

## Choice of scores
n <- nrow(USArrests)
a1 <- with(USArrests, qnorm(rank(Assault)/(n+1)))
a2 <- with(USArrests, qnorm(rank(UrbanPop)/(n+1)))
independence_test(a1~a2, USArrests) ## Van der Waerden
a1 <- with(USArrests, rank(Assault) > (n+1)/2)
a2 <- with(USArrests, rank(UrbanPop) > (n+1)/2)
independence_test(a1~a2, USArrests) ## Quadrant

##############################
## Multivariate: Ozone data ##
##############################

## Base R
boxplot(Ozone ~ Month, data = airquality)
summary(lm(Ozone ~ factor(Month), data = airquality))
kruskal.test(Ozone ~ Month, data = airquality)

## coin
kruskal_test(Ozone ~ factor(Month), data = airquality)
independence_test(rank(Ozone) ~ factor(Month), data = subset(airquality, !is.na(Ozone)), teststat="quad")

## Other scores
n <- sum(!is.na(airquality$Ozone))
independence_test(qnorm(rank(Ozone)/(n+1)) ~ factor(Month), data = subset(airquality, !is.na(Ozone)), teststat="quad")
independence_test(rank(Ozone) > (n+1)/2 ~ factor(Month), data = subset(airquality, !is.na(Ozone)), teststat="quad")

## Max vs. quadratic
independence_test(rank(Ozone) ~ factor(Month), data = subset(airquality, !is.na(Ozone)), teststat="quad")
independence_test(rank(Ozone) ~ factor(Month), data = subset(airquality, !is.na(Ozone)), teststat="max")

## By hand
Z <- with(subset(airquality, !is.na(Ozone)), model.matrix(~0+factor(Month)))
r <- with(subset(airquality, !is.na(Ozone)), rank(Ozone))
u <- crossprod(Z, r)
u0 <- mean(r) * apply(Z, 2, sum)
V <- as.numeric(var(r)) * crossprod(scale(Z, scale=FALSE))
Vinv <- MASS:::ginv(V)
teststat <- crossprod(u-u0, Vinv %*% (u-u0))
1-pchisq(teststat, df=4)

###################################
## Friedman test: Rounding times ##
###################################
rounding <- read.delim("http://web.as.uky.edu/statistics/users/pbreheny/621/data/rounding.txt")
boxplot(times~methods, rounding, col="gray")
boxplot(times~block, rounding, col="gray")
friedman.test(times~methods|block, rounding)
rounding$block <- factor(rounding$block)
friedman_test(times~methods|block, rounding)
kruskal.test(times~methods, rounding)