load(url('http://myweb.uiowa.edu/pbreheny/data/Golub1999.RData'))
source('http://myweb.uiowa.edu/pbreheny/7600/s16/notes/fun.R')

# Slide 6
summ <- summary(lm(X~y))
tstat <- sapply(summ, function(s) s$coef[2,3])
p <- sapply(summ, function(s) s$coef[2,4])
h <- length(p)
n <- length(y)

# Slide 7
z <- qnorm(p/2) * sign(-tstat)

# Slide 8
hist(p, main='', col='gray', border='white', las=1, breaks=seq(0, 1, .025), freq=FALSE)
lines(0:1, c(1,1), lwd=3, col='slateblue')

# Slide 9
hist(z, main='', col='gray', border='white', las=1, breaks=seq(-6, 8.5, 0.5),
     freq=FALSE, ylim=c(0, 0.4), xlab='z')
zz <- seq(-6, 8.5, .01)
lines(zz, dnorm(zz), lwd=3, col='slateblue')

# Slide 13 (adds to the figure from Slide 9)
lines(rep(qnorm(.05/h), 2), c(0, 0.4), lwd=2, lty=2)
lines(rep(-qnorm(.05/h), 2), c(0, 0.4), lwd=2, lty=2)
sum(p < .05)
sum(p < .05/h)

# Slide 18
sum(p.adjust(p, method='bonf') < .05)
sum(p.adjust(p, method='holm') < .05)

# Westfall-Young
# This one takes several minutes; you may want to set B = 25 if you're in a hurry
# There are faster ways to carry out the calculations here, but I felt they
# make it harder to see what's happening
B <- 500
minP <- matrix(NA, B, h)
pb <- txtProgressBar(1, B, style=3)
ind <- order(p)
XX <- X[,ind]
for (i in 1:B) {
  summ.i <- summary(lm(XX~sample(y)))
  p.i <- sapply(summ.i, function(s) s$coef[2,4])
  minP[i,] <- rev(cummin(rev(p.i)))
  setTxtProgressBar(pb, i)
}
pi.hat <- apply(t(minP) <= p[ind], 1, mean)
p.adj.sorted <- cummax(pi.hat)
p.adj <- p.adj.sorted[match(1:h, ind)]
sum(p.adj < .05)

# Slide 23
P <- cbind(p.adjust(p, method='bonf'), p.adjust(p), p.adj)
S <- apply(P, 2, sort)
matplot(S[1:500,], col=pal(3), type='l', bty='n', lwd=3, lty=1, xlab="Ordered index", ylab='Adjusted p', las=1)
toplegend(legend=c('Bonferroni', 'Holm', 'Westfall'), lwd=3, col=pal(3))