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

# Leukemia example
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)

# Slide 9
minP <- readRDS('minP.rds')   # From previous lecture
ind <- order(p)
pi.hat <- apply(t(minP) <= p[ind], 1, mean)
p.adj.sorted <- cummax(pi.hat)
p.wy <- p.adj.sorted[match(1:h, ind)]
P <- cbind(p.adjust(p, method='bonf'), p.wy, p.adjust(p, method='BH'))
colnames(P) <- c('Bonferroni', 'Westfall', 'Benjamini')
S <- apply(P, 2, sort)
col <- c("#FF4E37FF", "#00B500FF", "#008DFFFF")
matplot(S[1:2000,], col=col, type='l', bty='n', lwd=3, lty=1, xlab="Ordered index", ylab='FWER or FDR', las=1)
hdrm:::toplegend(legend=colnames(P), lwd=3, col=col)
apply(P <= 0.1, 2, sum)

# Slide 14
ph <- function(tt) sapply(tt, function(t) sum(p>t)/(h*(1-t)))
tt <- seq(.01, .99, .01)
plot(tt, ph(tt), pch=19, xlab='t', ylab=expression(hat(pi)(t)), bty='n', las=1, ylim=c(0.4, 1), col='gray')
fit <- smooth.spline(tt, ph(tt), df=6)
lines(fit$x, fit$y, col='red', lwd=3)
pi.hat <- tail(fit$y, 1)
h0.hat <- pi.hat*h

# Slide 15
H <- H0 <- hist(p, breaks=seq(0,1,.025), plot=FALSE)
H0$counts <- pmin(h0.hat*.025, H$counts)
plot(H, main="", xlab="p", ylim=c(0, max(H$counts)), border="white", col="gray", las=1)
plot(H0, add=TRUE, border="white", col="gray50")
lines(0:1, rep(h0.hat*.025, 2), lwd=3, col='red')