#Example 3.2 Government vs Private Salaries
gov.data <- read.table("http://www.rohan.sdsu.edu/~babailey/stat672/t3-2.txt", header=T)
gov.data
attach(gov.data)
#From Charlie Geyer's code: Wilcoxon Sign Rank Test

mu <- 0 # hypothesized value of median
z <- sort(private - government)
z <- z - mu
z <- z[z != 0]
n <- length(z)
r <- rank(abs(z))
print(rbind(z, r))
print(tplus <- sum(r[z > 0]))
psignrank(floor(tplus) - 1, n, lower.tail=FALSE)

#Note there are ties, so wrong p-value!

#Sample median of Walsh averages

print(z <- sort(private - government))
walsh <- outer(z, z, "+") / 2
print(walsh <- sort(walsh[lower.tri(walsh, diag = TRUE)]))
median(walsh)

#Confidence Intervals
conf.level <- 0.95
print(z <- sort(private - government))
print(n <- length(z))
walsh <- outer(z, z, "+") / 2
print(walsh <- sort(walsh[lower.tri(walsh, diag = TRUE)]))
print(m <- length(walsh))
alpha <- 1 - conf.level
k <- qsignrank(alpha / 2, n)
if (k == 0) k <- k + 1
print(k)
cat("achieved confidence level:",1 - 2 * psignrank(k - 1, n), "\n")
c(walsh[k], walsh[m + 1 - k])


#or can use: 
wilcox.test(private, government, paired=TRUE, alternative="greater")
#Note there are ties, so wrong p-value!
wilcox.test(private, government, paired=TRUE, conf.int=TRUE)
#or the exact test:
library(exactRankTests)
wilcox.exact(private, government, paired = TRUE, alternative = "greater")
wilcox.exact(private, government, paired = TRUE, conf.int=TRUE)