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

##
## a: bootstrap
##

## lengths of fish caught in the Chesapeake Bay
y <- c(9.0, 46.2, 9.0, 13.8, 11.5, 18.3, 15.2, 21.0, 11.2, 30.9, 7.1, 
       15.2, 18.1, 25.3, 42.6, 32.8, 11.9, 14.3, 8.3, 11.0, 10.5, 10.7, 
       36.4, 11.6, 21.7, 26.5, 25.6, 18.8, 55.7, 24.8, 18.5)
n <- length(y)
n

## calculating medians
ymed <- median(y)
ymed

## single Bootstrap sample median
Ystar <- sample(y, n, replace=TRUE)    ## always take n = length(y) samples
Ymed <- median(Ystar)                  ## this is s(Ystar)
Ymed

## bootstrap sample has duplicates/missing some
setdiff(y, Ystar)
table(round(Ystar, 2))

## full generic Bootstrap resample
B <- 10000
Ystar <- matrix(NA, nrow=B, ncol=n)
for(b in 1:B) {
  Ystar[b,] <- sample(y, n, replace=TRUE)
}

## calculating Bootstrap median
Ymeds <- apply(Ystar, 1, median)
table(Ymeds)

## visualizing the empirical/sampling distribution
hist(Ymeds)
## clearly lower resolution than a simulation from a model

## this is because the median itself is low resolution

## two options when using a Bootstrap to form a CI

## 1: direct quantiles
p <- c(0.025, 0.975)
quantile(Ymeds, p)
## less desirable given the low resoultion sampling distribution

## 2: smooth over with the CLT (technically just a Gaussian approx.)
q <- qnorm(p)
mean(Ymeds) + q*sd(Ymeds)

## with sd instead of median
Ysds <- apply(Ystar, 1, sd)
hist(Ysds)
quantile(Ysds, p)
mean(Ysds) + q*sd(Ysds)

## 
## b: permutation tests
##

## more fish, caught on another day
x <- c(23.4, 18.1, 37.8, 32.5, 26.1, 22.0, 15.0, 32.2, 37.5, 19.2, 15.9, 
       31.4, 13.0, 14.8, 24.1, 21.6, 34.6, 28.0, 27.1, 25.2, 26.1, 37.6, 
       33.1, 14.7, 25.0)

## form the combined sample and calculat ybar - xbar stat
nx <- length(x)
ny <- length(y)     ## formerly n, but now n is for combined sample
yx <- c(y, x)
n <- length(yx)                                   ## same as nx + ny
bdiff <- mean(yx[1:ny]) - mean(yx[(ny + 1):n])    ## this is s(.)
c(bdiff, mean(y) - mean(x))

## lots of random permutations
P <- 10000
Bds <- rep(NA, P)
for(p in 1:P) {
  YXs <- sample(yx, n)                               ## rand permutation
  Bds[p] <- mean(YXs[1:ny]) - mean(YXs[(ny + 1):n])  ## s(YXs)
}

## familiar visual
hist(Bds, main="", xlab="differences in means", xlim=c(-11, 13))
abline(v=bdiff, col=2, lwd=2)
abline(v=-bdiff, col=2, lty=2, lwd=2)
legend("topright", c("obs", "reflect"), lty=1:2, lwd=2, col=2, bty="n")

## p-value calculation
pval.bd <- 2*mean(Bds <= bdiff)
pval.bd
## fail to reject H0 with this stat

## what about another stat, like ratio of standard deviation
sdiff <- sd(yx[1:ny]) - sd(yx[(ny + 1):n])

### another loop
Sds <- rep(NA, P)
for(p in 1:P) {
  YXs <- sample(yx, n)
  Sds[p] <- sd(YXs[1:ny]) - sd(YXs[(ny + 1):n])
}

## p-value
pval.sd <- 2*mean(Sds >= sdiff)
pval.sd
## so the distributions are indeed different