## Chapter 9
##
## Code to support Video Vignette
##
## Videos and supporting code are not a complete portrayal
## of Chapter content.

##
## a: Wilcoxon rank sum test
##

## some "data" values
y <- c(73, 14, 6, 9, 85, 3, 12, 37, 36, 23, 11)

## ranks
rank(y)
order(order(y))

## tied ranks for duplicates
y <- c(y, y[2], y[7], y[7]) 
rank(y)

## number of VA (y) and MD (x) Washington Metro
## riders who said that they used a prepaid paper fare card
y <- c(80, 76, 56, 67, 73, 58, 51, 65, 68, 61)
x <- c(83, 66, 71, 82, 81, 89, 97, 59, 74)

## extract sizes
ny <- length(y)
nx <- length(x)

## rank the combined sample
yx <- c(y, x)
n <- ny + nx
ryx <- rank(yx) 
ryx

## extract the parts that go with each original sample
ry <- ryx[1:ny]
rx <- ryx[(ny + 1):n]

## test statistic (rank sum)
rbar <- sum(ry)
rbar

## simulating ranks involves permuting {1,...,n}
N <- 1000000
Rbars <- rep(NA, N)
for(i in 1:N) {
  Ryxs <- sample(1:n, n)
  Rys <- Ryxs[1:ny]
  Rbars[i] <- sum(Rys)
}

## visualize
hist(Rbars, main="") 
abline(v=c(rbar, 2*mean(Rbars) - rbar), col=2, lty=1:2, lwd=2)
legend("topright", c("obs", "reflect"), lty=1:2, lwd=2, col=2, bty="n")

## p-value calculation
2*mean(Rbars <= rbar)

## the math way
2*pwilcox(rbar - ny*(ny + 1)/2, ny, nx)


##
## b: Wilcoxon signed ranks test
##

## 12 mice given a drug to improve their blood circulation,
## measured, in milliliters per minute, before (x) and after (y)
## administering a drug
circ <- data.frame(
  after=c(0.0541, 0.0196, 0.0234, 0.0217, 0.0165, 0.0305, 0.0351, 0.0155, 
          0.0272, 0.02425, 0.02605, 0.0513),
  before=c(0.0335, 0.0240, 0.0375, 0.0133, 0.0305, 0.0173, 0.0138, 0.0259, 
           0.0134, 0.0365, 0.0205, 0.0352))
y <- circ$after
x <- circ$before

## calculate differences and remove zeros
nprime <- length(y)
d <- y - x
d <- d[d != 0]
n <- length(d)

## signed ranks statistic
r <- rank(abs(d))
sr <- r*sign(d)               ## sr for r^{+/-} 
rbarp <- sum(sr[sr > 0])
rbarp

## signed ranks simulation. Sampling ranks isn't
## really necessary. The random sign locations is
## what's important.  
Rbarps <- rep(NA, length=N)
for(i in 1:N) {
  Rs <- sample(1:n, n)
  Signs <- sample(c(-1 ,1), n, replace=TRUE)
  SRs <- Rs*Signs
  Rbarps[i] <- sum(SRs[Signs > 0])
}

## skipping the visual, but it helps to check
## what tail you're in
c(rbarp, mean(Rbarps))

## right tail
2*mean(Rbarps >= rbarp)

## right-tail calculations must subtract 1 due to how psignrank
## counts from the left
2*psignrank(rbarp - 1, n, lower.tail=FALSE)

##
## one-sample signed ranks test
##

## IQ scores for kids in Blacksburg
y <- c(95, 101, 100, 115, 105, 99, 101, 107, 120, 100, 98, 77, 89, 121,
       102, 97, 106, 111, 101, 107, 102, 88, 112, 128)

## follow the signed ranks test for H0: E(Y) = mu = 100
x <- 100 

## cutting and pasting from above
nprime <- length(y)
d <- y - x
d <- d[d != 0] 
n <- length(d)
r <- rank(abs(d))
sr <- r*sign(d)
rbarp <- sum(sr[sr > 0])
rbarp

## more cutting and pasting
Rbarps <- rep(NA, length=N)
for(i in 1:N) {
  Rs <- sample(1:n, n)
  Signs <- sample(c(-1 ,1), n, replace=TRUE)
  SRs <- Rs*Signs
  Rbarps[i] <- sum(SRs[Signs > 0])
}

## skipping the visual, but it helps to check
## what tail you're in
c(rbarp, mean(Rbarps))

## right tail
2*mean(Rbarps >= rbarp)

## right-tail calculations must subtract 1 due to how pwilcox
## counts from the left
2*psignrank(rbarp - 1, n, lower.tail=FALSE)