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

##
## a: non-P Correlation
##

## math testing data, check association between score 
## and socio-economic status (SES)
nels <- read.csv("nels_math.csv")
s1 <- which(nels$school == 1)
y <- nels$score[s1]
x <- nels$ses[s1]

## 
## Spearman's rho
##

## Spearman correlation
rshat <- cor(rank(y), rank(x)) 
rshat
cor(y, x, method="spearman")

## simulation
n <- length(y)
N <- 1000000
Rs <- rep(NA, N)
for(i in 1:N) {
  RYs <- sample(1:n, n)      ## random Y ranks
  RXs <- sample(1:n, n)      ## random X ranks
  Rs[i] <- cor(RXs, RYs)     ## Pearson correlation on ranks
}

## visualize
hist(Rs, main="", xlim=c(-0.75, 1)) 
abline(v=c(rshat, -rshat), lty=1:2, col=2, lwd=2)
legend("topright", c("obs", "reflect"), lty=1:2, lwd=2, col=2, bty="n")

## p-value
2*mean(Rs >= rshat)

## the math way, first do install.packages("SuppDists")
library(SuppDists)
2*pSpearman(rshat, n, lower.tail=FALSE)
          
##
## Kendall's tau
## 

## paste from the book for visual of concordant and discordant pairs
plot(x, y, xlab="SES", ylab="math test score")
i <- 1
abline(h=y[i], lty=2, lwd=2)
abline(v=x[i], lty=2, lwd=2)
inf <- 1000
gx <- c(-inf, x[i], x[i], inf, inf, x[i], x[i], inf)
gy <- c(y[i], y[i], inf, inf, y[i], y[i], -inf, -inf)
polygon(gx, gy, col=3, density=10, angle=45)
point <- paste0("(x[", i, "], y[", i, "])")
legend(x=0.2, y=45, c(point, "concordant"), col=c(1, 3), 
       lty=2:1, lwd=2:1, bty="n")

## there are ties in these data, but the library 
## symmetrizes the calculation
tauhat <- cor(y, x, method="kendall")
tauhat
cor(x, y, method="kendall")

## simulate pairs of anything, like "ranks"
Taus <- rep(NA, N)
for(i in 1:N) {
  Ys <- sample(1:n, n)                      ## random Ys (or ranks)
  Xs <- sample(1:n, n)                      ## random Xs (or ranks)
  Taus[i] <- cor(Xs, Ys, method="kendall")  ## Kendall's tau calculation
}

## p-value
2*mean(Taus >= tauhat)

## math way
2*pKendall(tauhat, n, lower.tail=FALSE)

##
## b: non-P SLR
##

## data on precipitation over last 100 years
precip <- read.csv("usprecip.csv")
x <- precip$date
y <- precip$inches
n <- length(y)
plot(x, y, type="b", xlab="year", ylab="inches") 

## testing

## test for a slope of H0: beta1 = 0.02, which would mean
## a 2" rise in average precipitation in 100 years
beta1 <- 0.02
u <- y - beta1*x
rs.ux <- cor(u, x, method="spearman")

## via simulation
Rs <- rep(NA, N)
for(i in 1:N) {
  RYs <- sample(1:n, n)
  RXs <- sample(1:n, n)
  Rs[i] <- cor(RXs, RYs)
}

## p-value
2*mean(Rs <= rs.ux)    ## note left tail
2*pSpearman(rs.ux, n)
## cannot reject H0 that rize is 2" in 100 years

## CI

## shorthand for all slopes
s <- outer(y, y, '-')/outer(x, x, '-')
dim(s)  ## too big/has everything twice
s <- s[upper.tri(s)]

## some denominators could be zero
s <- s[is.finite(s)]    ## there aren't any in the precip example
ns <- length(s)
c(ns, n*(n - 1)/2)

## quantiles of sampling distribution for tau, 
## skipping simulation version (see book)
alpha <- 0.05
q <- qKendall(1 - alpha/2, n)

## many sensible ways to use q to map to indicies 1, ..., ns
qr <- ceiling(ns*(q + 1)/2)          ## ceiling for integer
ql <- ns - qr

## use those indices on s in sorted order
so <- sort(s)
CI <- c(so[ql], so[qr])
CI
## does not include 0, but does include 0.02

## prediction

## point estimate
b1 <- median(s)
medx <- median(x)
b0 <- median(y) - b1*medx
c(b0=b0, b1=b1)
abline(b0, b1)

## but what about uncertainty?
## see book for simulation that combines Kendall's tau
## with bootstrap median