#Problem 4.8: data with ties
#Using Charlie Geyer's Computer Code
exdata <- read.table("http://edoras.sdsu.edu/~babailey/stat672/4.8.data", header=T)
attach(exdata)
exdata

#Wilcoxon Rank Sum Test
mu <- 0 # hypothesized value of median difference
x <- x[!is.na(x)]
y <- y[!is.na(y)]
print(nx <- length(x))
print(ny <- length(y))
y <- y - mu
data <- c(x, y)
names(data) <- c(rep("x", nx), rep("y", ny))
data <- sort(data)
r <- rank(data)
rbind(data, r)
print(w <- sum(r[names(data) == "y"]))
print(u <- w - ny * (ny + 1) / 2)
pwilcox(u, nx, ny)

#Hodges-Lehmann estimator: median of the pairwise difference 
x <- x[!is.na(x)]
y <- y[!is.na(y)]
print(diffs <- sort(as.vector(outer(y, x, "-"))))
median(diffs)

#Associated CI: counts in k from each end of the list of 
#sorted pairwise differences
conf.level <- 0.95
x <- x[!is.na(x)]
y <- y[!is.na(y)]
print(nx <- length(x))
print(ny <- length(y))
print(diffs <- sort(as.vector(outer(y, x, "-"))))
print(m <- length(diffs))
alpha <- 1 - conf.level
k <- qwilcox(alpha / 2, nx, ny)
if (k == 0) k <- k + 1
print(k)
cat("achieved confidence level:",
    1 - 2 * pwilcox(k - 1, nx, ny), "\n")
c(diffs[k], diffs[m + 1 - k])

#Or use the wilcox.exact test
library("exactRankTests")
wilcox.exact(y, x, alternative="less")
wilcox.exact(y, x, conf.int=T)