library(limma)

#	Simulation of 10K microarray with 1% non-null effects with SNR 3


m		<- 10^4
beta.vec	<- c(rep(0,0.99*m),rnorm(0.01*m,mean = 3,sd = 0.1))
H.vec		<- c(rep(0,0.99*m),rep(1,0.01*m))


d0		<- 5
s0		<- 1

sig.g.vec	<- sqrt(d0*s0^2 / rchisq(m,5))
summary(sig.g.vec)


#	Experiment with 4 replicates


y.array	<- array(rnorm(m*4,mean = rep(beta.vec,4),sd =  rep(sig.g.vec,4)),dim=c(m,4))

beta.hat.vec	<- apply(y.array,1,mean)
sg.vec		<- apply(y.array,1,sd)

summary(sg.vec)

#	==>	v.g = 4, dg = 3

t.vec	<- beta.hat.vec / (sg.vec / 2)

summary(t.vec[H.vec == 0])
summary(t.vec[H.vec == 1])

p.t.vec	<- 2*(1 - pt(abs(t.vec),3))

par(mfcol=c(1,2))
plot(sort(t.vec))
plot(sort(p.t.vec))

sum(p.adjust(p.t.vec[H.vec == 0],method = "BH") < 0.05)
sum(p.adjust(p.t.vec,method = "BH") < 0.05)

qt(1 - 0.025/m,3)


#	Let's use LIMMA's ebayes function


#  ? squeezeVar

eBayes.est	<- squeezeVar(sg.vec^2,3)

str(eBayes.est)

mod.t.vec	<- beta.hat.vec / sqrt(eBayes.est$var.post  / 4)

par(mfcol=c(1,1))
plot(t.vec,mod.t.vec)
abline(0,1,col=2)



summary(mod.t.vec[H.vec == 0])
summary(mod.t.vec[H.vec == 1])

p.mod.t.vec	<- 2*(1 - pt(abs(mod.t.vec),3+eBayes.est$df.prior))


par(mfcol=c(1,2))
plot(sort(mod.t.vec))
plot(sort(p.mod.t.vec))

qt(1 - 0.025/m,3+5)
qt(1 - 50*0.025/m,3+5)

sum(p.adjust(p.mod.t.vec[H.vec == 0],method = "BH") < 0.05)

(R.BH	<- sum(p.adjust(p.mod.t.vec,method = "BH") < 0.05))
(V.BH	 <-sum(p.adjust(p.mod.t.vec,method = "BH") < 0.05 & H.vec == 0))