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

# Leuk example
summ <- summary(lm(X~y))
tstat <- sapply(summ, function(s) s$coef[2,3])
p <- sapply(summ, function(s) s$coef[2,4])
z <- qnorm(p/2) * sign(-tstat)
h <- length(p)

# Slide 3
q <- p.adjust(p, method='BH')
sum(q < .01)
range(-qnorm(p/2)[q < .01])

# Slide 8
(h/100)/sum(p < .02 & p > .01)

# Note: the function lfdrPlot is defined in
# myweb.uiowa.edu/pbreheny/7600/s16/notes/fun.R

# Slide 10
lfdr1 <- lfdrPlot(z)

# Slide 12
lfdr2 <- lfdrPlot(z, pi0=.56)

# z Slide13
f1 <- approxfun(z, lfdr1)
f2 <- approxfun(z, lfdr2)
zz <- seq(-5, 5, len=99)
plot(zz, f1(zz), type='l', lwd=3, bty='n', las=1, ylim=c(0, 1), col=pal(2)[1],
     xlab='z', ylab='Local FDR')
lines(zz, f2(zz), lwd=3, col=pal(2)[2])
toplegend(legend=c(expression(hat(pi)[0]*'='*1), expression(hat(pi)[0]*'='*0.56)),
          lwd=3, col=pal(2))
lines(zz, rep(0.2, length(zz)), lty=2)

uniroot(function(x) f1(x)-0.2, lower=2, upper=4)$root
uniroot(function(x) f2(x)-0.2, lower=2, upper=4)$root
sum(lfdr1 < .2)
sum(lfdr2 < .2)

# Slide 14
require(locfdr)
res <- locfdr(z)
head(res$fdr)

# Slide 15
sum(q < .1)
sum(lfdr1 < .1)
min(abs(z[q < .1]))
min(abs(z[lfdr1 < .1]))

# Slide 16
F <- ecdf(z)
zz <- seq(-6, -2, len=99)
plot(F(zz), pnorm(zz), type='l', lwd=3, bty='n', las=1, ylim=c(0, 0.025),
     xlab='F(z)', ylab=expression(F[0](z)))
zz <- quantile(F, .08)
lines(c(0, 0.08), c(0, pnorm(zz)), lwd=3, col=pal(2)[1])
dens <- density(z, bw="nrd")
f <- approxfun(dens$x, dens$y)
slope <- dnorm(zz)/f(zz)
lines(.08+c(-1,1)*pnorm(zz)/slope, c(0, 2*pnorm(zz)), lwd=3, col=pal(2)[2])

# Slide 18
p.left <- pnorm(z)
FDR <- h*p.left/rank(p.left)
mean(lfdr1[FDR < .1])
p.right <- pnorm(-z)
FDR <- h*p.right/rank(p.right)
mean(lfdr1[FDR < .1])

# Slide 20
x <- c(1, 3, 5, 6, 6.2, 6.6, 7, 8, 10)
n <- length(x)
xx <- seq(-2, 10, len=101)
fn <- paste("kernel", 1:3, ".pdf", sep="")
hist(x, border="white", main="", freq=FALSE, las=1)
rug(x, col="red", lwd=3)
for (i in 1:n) lines(xx, dnorm(xx, x[i], 0.75)/10, col="#008DFFB2", lwd=3)
d <- density(x, bw=0.75)
lines(d, lwd=3)

# Slide 21
brk <- seq(min(z), max(z), length = 99)
h <- hist(z, breaks=brk, main="", border=FALSE, col="lightgray", freq=FALSE, las=1)
dens <- density(z, bw=.05)
lines(dens, col='slateblue', lwd=3)
h <- hist(z, breaks=brk, main="", border=FALSE, col="lightgray", freq=FALSE, las=1)
dens <- density(z)
lines(dens, col='slateblue', lwd=3)
h <- hist(z, breaks=brk, main="", border=FALSE, col="lightgray", freq=FALSE, las=1)
dens <- density(z, bw=1)
lines(dens, col='slateblue', lwd=3)

# Slide 22
res <- locfdr(z, nulltype=0)
lfdr <- lfdrPlot(z, pi0=res$fp0[1,3])
cor(res$fdr, lfdr)

# Slide 25
res$fp0[c(1,3,5),]

# Slide 26
res <- locfdr(z, plot=FALSE)
lfdr <- lfdrPlot(z, pi0=0.56)
lfdr <- lfdrPlot(z, sigma=res$fp0[5,2], pi0=res$fp0[5,3])
lfdr <- lfdrPlot(z, delta=res$fp0[5,1], sigma=res$fp0[5,2], pi0=res$fp0[5,3])