Golub1999 <- readRDS(url('https://s3.amazonaws.com/pbreheny-data-sets/Golub1999.rds'))

# Slide 6
summ <- summary(lm(Golub1999$X~Golub1999$y))
tstat <- sapply(summ, function(s) s$coef[2,3])
p <- sapply(summ, function(s) s$coef[2,4])
h <- length(p)
n <- length(Golub1999$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 (zooming in here on the z-values between -6 and 6)
z6 <- z[z > -6 & z < 6]
zz <- seq(-6, 6, .01)
hist(z6, main='', col='gray', border='white', las=1, breaks=seq(-6, 6, 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 <- Golub1999$X[,ind]
set.seed(6)  # Feel free to use a diffferent seed but this is the one
             # that will give the results that appear in the notes
for (i in 1:B) {
  summ.i <- summary(lm(XX~sample(Golub1999$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)
col <- c("#FF4E37", "#00B500", "#008DFF")
matplot(S[1:500,], col=col, type='l', bty='n', lwd=3, lty=1, xlab="Ordered index", ylab='Adjusted p', las=1)
hdrm:::toplegend(legend=c('Bonferroni', 'Holm', 'Westfall'), lwd=3, col=col)