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

# Slide 9
p0 <- runif(1000)
A <- function(tt) sapply(tt, function(t) sum(p0<t))
tt <- seq(1, .01, len=999)
M <- A(tt)/tt
dev <- max(abs(M-length(p0)))
ylim <- c(length(p0)-dev, length(p0)+dev)
plot(tt, M, type='l', lwd=3, bty='n', las=1, xlab='t', ylab='A(t)/t', xlim=rev(range(tt)), ylim=ylim)

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

# Slide 10 (requires running the code in 1-25 if you want to inclue Westfall-Young)
P <- cbind(p.adjust(p, method='bonf'), p.wy, p.adjust(p, method='BH'))
S <- apply(P, 2, sort)
matplot(S[1:2000,], col=pal(3)[c(1,2,3)], type='l', bty='n', lwd=3, lty=1, xlab="Ordered index", ylab='FWER or FDR', las=1)
toplegend(legend=c('Bonferroni', 'Westfall', 'Benjamini'), lwd=3, col=pal(3))

# Slide 12
q <- p.adjust(p, method='BH')
sum(q < .1)

# Slide 15
ph <- function(tt) sapply(tt, function(t) sum(p>t)/(h*(1-t)))
tt <- seq(.01, .99, .01)
par(mar=c(5, 5, 0.5, 0.5))
plot(tt, ph(tt), pch=19, xlab='t', ylab=expression(hat(pi)(t)), bty='n', las=1, ylim=c(0.4, 1))
pi.hat <- ph(.6)
h0.hat <- pi.hat*h

# Slide 16
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')