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

##
## a: Squared Ranks test
##

## numbers of bicycles on Green (y) and Orange (x) line 
## Metro trains
y <- c(46, 51, 61, 70, 66, 69, 50, 44, 58, 59, 49, 54)
x <- c(74, 72, 90, 51, 65, 63, 100, 57, 82)

## dimensions
ny <- length(y)
nx <- length(x)

## absolute deviations
v <- abs(y - mean(y))
u <- abs(x - mean(x))

## ranking the combined samplke
r <- rank(c(v, u))

## calculate the test stastic
rv <- r[1:ny]
r2bar <- sum(rv^2)

## simulation
N <- 1000000
n <- length(r)
R2bars <- rep(NA, N)
for(i in 1:N) {
  Rvus <- sample(1:n, n)
  Rvs <- Rvus[1:ny]
  R2bars[i] <- sum(Rvs^2)
}

## skipping the visual (but always check what side), p-value
2*mean(R2bars < r2bar)

## math way requires SSP, download code from book webpage
## requires install.packages(c("FLSSS", "RcppAlgos"))
source("ranks.R")

## closed form calculation, can be very slow if ny or nx is large
2*psqranks(r2bar, ny, nx)

##
## B: Kruskal-Wallis
##

## combining all of the farecard data
ylist <- list(
  VA=c(80, 76, 56, 67, 73, 58, 51, 65, 68, 61),
  DC=c(83, 66, 71, 82, 81, 89, 97, 59, 74),
  MD=c(67, 94, 83, 98, 35, 73, 29, 36, 60, 105, 34, 84, 89, 76, 79, 92, 
       49, 97, 46, 32, 104, 60, 61, 59, 98, 101),
  DC2=c(107, 114, 87, 102, 85, 94, 109, 91, 99, 59, 97, 92, 103, 90, 43, 
        108, 61, 111, 87, 112, 109, 115, 89, 105, 109, 83, 35, 114))

## don't forget your ANOVA shortcuts
nj <- sapply(ylist, length)
m <- length(nj)
n <- sum(nj)

## ranking the combined sample
r <- rank(unlist(ylist))

## here is a tidy way to extract the rank sums for each group
g <- factor(rep.int(seq_len(m), nj))
rg <- tapply(r, g, sum)

## and then aggregating their size-adjusted squares
r2bar <- sum(rg^2/nj)

## simulating
for(i in 1:N) {
  Rs <- sample(1:n, n)                      ## random permutation
  Rgs <- tapply(Rs, g, sum)                 ## break into groups and sum             
  R2bars[i] <- sum(Rgs^2/nj)                ## calculate statistic
}

## visualizing
hist(R2bars, main="") 
abline(v=r2bar, col=2, lwd=2)
legend("top", "obs", lwd=2, col=2, bty="n")

## like ANOVA, the p-value is right-tailed only
mean(R2bars > r2bar)

## CLT-likek approximation
denom <- (sum(r^2) - n*(n + 1)^2/4)/(n - 1)
x2 <- (r2bar - n*(n + 1)^2/4)/denom

## p-value
pchisq(x2, m - 1, lower.tail=FALSE)