## wilcox.test: Lead IQ data
lead <- read.delim("http://myweb.uiowa.edu/pbreheny/data/lead-iq.txt")
t.test(IQ~Smelter, lead)
wilcox.test(IQ~Smelter, lead)
x <- with(lead, IQ[Smelter=="Near"])
y <- with(lead, IQ[Smelter=="Far"])
wilcox.test(x, y) ## Same thing

## Adding an outlier
x <- c(x, 210)
t.test(x, y)$p.value
wilcox.test(x, y)$p.value

## Exact method: Problem with ties
wilcox.test(IQ~Smelter, lead, exact=TRUE)

## Exact vs approximate
lead2 <- lead
lead2$IQ <- lead2$IQ + runif(nrow(lead2), -0.1, 0.1)
wilcox.test(IQ~Smelter, lead2, exact=TRUE)     ## Exact result
wilcox.test(IQ~Smelter, lead2)                 ## Approximation 1
wilcox.test(IQ~Smelter, lead2, correct=FALSE)  ## Approximation 2

## Exact vs approximate: Time
x <- runif(200)
y <- runif(200)
system.time(wilcox.test(x, y, exact=TRUE))
system.time(wilcox.test(x, y))

## Semiparametric interval
x <- with(lead, IQ[Smelter=="Near"])
y <- with(lead, IQ[Smelter=="Far"])
wilcox.test(x, y, mu=3) ## This is the test
wilcox.test(x, y+3)     ## we are inverting
wilcox.test(x, y, conf.int=TRUE)

## Do-it-yourself signed rank test
cf <- read.delim("http://myweb.uiowa.edu/pbreheny/data/cystic-fibrosis.txt")
Diff <- apply(cf, 1, diff)
r <- rank(abs(Diff))  ## Absolute ranks
Obs <- sum(r[Diff>0]) ## Sum of positive ranks
n <- nrow(cf)
N <- 10000
ts <- numeric(N)
for (i in 1:N) {
  s <- rbinom(n, 1, 0.5) ## Random signs
  ts[i] <- sum(r[s==1])
}
hist(ts, col="gray", border="white", xlim=c(0,105), breaks=0:105, main="", xlab="Test statistic")
v <- c(mean(ts) - (Obs-mean(ts)), Obs)
abline(v=v, col="slateblue", lwd=3)
p <- sum(ts <= v[1] | ts >= v[2])/N
mtext(paste("p =",p), line=1)

## Do-it-yourself signed rank test, exact
l <- rep(list(0:1), n)
X <- expand.grid(l)
N <- nrow(X)
ts <- numeric(N)
source("http://myweb.uiowa.edu/pbreheny/571/f14/labs/displayProgressBar.R")
for (i in 1:N) {
  ts[i] <- sum(r[X[i,]==1])
  displayProgressBar(i, N)
}
tab <- table(ts)
obs <- as.character(Obs)
v <- as.numeric(names(tab)[tab==tab[obs]])
hist(ts, col="gray", border="white", xlim=c(0,105), breaks=0:105, main="", xlab="Test statistic")
abline(v=v, col="slateblue", lwd=3)
p <- sum(tab[tab <= tab[obs]])/N
mtext(paste("p =", round(p, 5)), line=1)

## wilcox.test: Signed
wilcox.test(Diff)

## wilcox.test: Confidence interval
wilcox.test(Diff, mu=4) ## Test of location
wilcox.test(Diff-4)     ## Same thing
wilcox.test(Diff, conf.int=TRUE)

## Bootstrapping the median
B <- 10000
b <- numeric(B)
for (i in 1:B) {
  x <- sample(Diff, replace=TRUE)
  b[i] <- median(x)
}
quantile(b, c(.025, .975))                ## Percentile interval
median(Diff) + qnorm(c(.025, .975))*sd(b) ## "z" interval