#-------------------
#	EXAMPLE 1
#-------------------

z.vec	<-  seq(10,-6,length = 10000)
pi.0	<- 0.80
mu.0	<- 0
mu.1	<- 4

f.0		<- dnorm(z.vec,mean=mu.0,sd=1)
f.1		<- dnorm(z.vec,mean=mu.1,sd=1)
f		<- pi.0*f.0 + (1-pi.0)*f.1

fdr.vec	<- pi.0*f.0 / f
Fdr.vec <- cumsum(fdr.vec * f) / cumsum(f)
(z.05	<- approx(Fdr.vec,z.vec,xout=0.05)$y)
(Fdr.1	<- approx(z.vec,Fdr.vec,xout=1)$y)

plot(z.vec,f.0,type="l",col=2,ylim=c(0,1),xlab="Z",ylab="",lwd=2)
lines(z.vec,f.1,type="l",col=3,lwd=2)
lines(z.vec,f,type="l",col=4,lwd=2)
lines(z.vec,fdr.vec,type="l",col=5,lwd=2)
legend(5,0.8,lwd=2,c("f.0","f.1","f","fdr"),col=c(2,3,4,5))

plot(z.vec,f,type="l",col=4,ylim=c(0,1),xlab="Z",ylab="Fdr/fdr",lwd=2)
lines(z.vec,fdr.vec,type="l",col=5,lwd=2)
lines(z.vec,Fdr.vec,type="l",col=6,lwd=2)
legend(5,0.8,lwd=2,c("f","fdr","Fdr"),col=c(4,5,6))
arrows(1,-.02,8,-.02,length=.05,lwd=3,col=4,code=2,angle=45)
abline(v=1,lty=2,col=4)
arrows(1,Fdr.1,-4,Fdr.1,length=.05,lwd=1,col=4,code=2,angle=45)


plot(z.vec,f,type="l",col=4,ylim=c(0,1),xlab="Z",ylab="Fdr/fdr",lwd=2)
lines(z.vec,Fdr.vec,type="l",col=6,lwd=2)
abline(h = 0.05,col=6,lty=2)
arrows(z.05,0.05,z.05,-.02,length=.05,lwd=1,col=6,code=2,angle=45)
arrows(z.05,-.02,8,-.02,length=.05,lwd=3,col=4,code=2,angle=45)


#	1.a Controlling 1D Bayes FDR

library(locfdr)
library(ggplot2)
library(qvalue)

m		<- 500000
H.vec	<- rbinom(m,1,1 - pi.0)
Z.vec	<- H.vec * mu.1 + rnorm(m)
P.vec	<- 1 - pnorm(Z.vec,mean=mu.0)

par(mfcol=c(1,1))
hist(Z.vec)


#	1.a Bayes classifier results

t.test(1 - H.vec[1 <= Z.vec])
t.test(1 - H.vec[z.05 <= Z.vec])

z.05
sum(z.05 <= Z.vec)
mean(1 - H.vec[z.05 <= Z.vec])
sum(1 - H.vec[z.05 <= Z.vec])

Z.fdr	<- approx(z.vec,fdr.vec,xout = Z.vec)$y
mean(Z.fdr[z.05 <= Z.vec])

#	1.b locfdr results

a 		<- locfdr(Z.vec,df=18)
z.locfdr.05	<- approx(a$mat[,4],a$mat[,1],xout = 0.05)$y

z.locfdr.05
sum(z.locfdr.05 <= Z.vec)
mean(1 - H.vec[z.locfdr.05 <= Z.vec])
sum(1 - H.vec[z.locfdr.05 <= Z.vec])

plot(a$mat[,1],a$mat[,4],xlab="z",ylab = "Fdrright",type="l",col=4,lwd=2)
lines(z.vec,Fdr.vec,type="l",col=6,lwd=2)
legend(3.5,0.8,lwd=2,c("Fdr","locfdr Fdr estimate"),col=c(6,4))


#	1.c qvalue results

b		<- qvalue(P.vec)
b$pi0
ord.z	<- order(Z.vec)
sum(b$qvalues <= 0.05)
z.qvalue.05	<- min(Z.vec[b$qvalues <= 0.05])
qvalue.z.Fdr	<- approx(Z.vec[ord.z],b$qvalues[ord.z],xout=z.vec)


z.qvalue.05
sum(z.qvalue.05 <= Z.vec)
mean(1 - H.vec[z.qvalue.05 <= Z.vec])
sum(1 - H.vec[z.qvalue.05 <= Z.vec])



plot(qvalue.z.Fdr,xlab="z",ylab = "Fdr",type="l",col=3,lwd=2)
lines(z.vec,Fdr.vec,type="l",col=6,lwd=2)
lines(qvalue.z.Fdr,col=3,lwd=2)
legend(3.5,0.8,lwd=2,c("Fdr","locfdr Fdr estimate"),col=c(6,3))
#lines(a$mat[,1],a$mat[,4],type="l",col=4,lwd=2)#

plot(qvalue.z.Fdr,xlab="z",ylab = "Fdr",type="l",col=3,lwd=2,xlim=c(2,3),ylim=c(0,0.08))
lines(z.vec,Fdr.vec,type="l",col=6,lwd=2)
lines(a$mat[,1],a$mat[,4],type="l",col=4,lwd=2)
abline(h = 0.05,col=6,lty=2)
arrows(z.05,0.05,z.05,0,length=.05,lwd=1,col=6,code=2,angle=45)
arrows(z.qvalue.05,0.05,z.qvalue.05,0,length=.05,lwd=1,col=3,code=2,angle=45)
arrows(z.locfdr.05,0.05,z.locfdr.05,0,length=.05,lwd=1,col=4,code=2,angle=45)




#	1.d  BH procedure results


P.adj	<- p.adjust(P.vec,method="BH")
sum(P.adj <= 0.05)
min(Z.vec[P.adj <= 0.05])
sum(1 - H.vec[P.adj <= 0.05])
mean(1 - H.vec[P.adj <= 0.05])


#	1.e  adaptive BH procedure results

pi.0.hat	<- 2 * mean(0.5 <= P.vec)

pi.0.hat
sum(P.adj <= 0.05/pi.0.hat)
min(Z.vec[P.adj <= 0.05/pi.0.hat])
sum(1 - H.vec[P.adj <= 0.05/pi.0.hat])
sum(1 - H.vec[P.adj <= 0.05/b$pi0])



#-------------------
#	EXAMPLE 2
#-------------------

x.vec	<- seq(-3,6,length=500)
y.vec	<- seq(-3,6,length=500)

mu.00	<- c(0,0)
mu.01	<- c(0,3)
mu.10	<- c(3,0)
mu.11	<- c(3,3)

f.00	<- outer(dnorm(x.vec,mean=mu.00[1]),dnorm(y.vec,mean=mu.00[2])) 
f.01	<- outer(dnorm(x.vec,mean=mu.01[1]),dnorm(y.vec,mean=mu.01[2])) 
f.10	<- outer(dnorm(x.vec,mean=mu.10[1]),dnorm(y.vec,mean=mu.10[2])) 
f.11	<- outer(dnorm(x.vec,mean=mu.11[1]),dnorm(y.vec,mean=mu.11[2])) 

Fdr.fdr		<- f.00 * 0
Fdr.p			<- f.00 * 0
Fdr.BH		<- f.00 * 0

px.outr		<- 1 - pnorm(outer(x.vec,rep(1,500)))
py.outr		<- 1 - pnorm(outer(rep(1,500),y.vec))
p.no.assoc	<- 1 - pchisq(-2*(log(px.outr) + log(py.outr)),2*2)
p.no.rep	<- pmax(px.outr,py.outr)



# 2.1 Plots for testing no association

pi.00	<- 0.85
pi.01	<- 0.05
pi.10	<- 0.05
pi.11	<- 0.05

f		<- pi.00*f.00 + pi.01*f.01 + pi.10*f.10 + pi.11*f.11
f.0		<- pi.00*f.00 
f.1		<- pi.01*f.01 + pi.10*f.10 + pi.11*f.11
fdr		<- f.0 / f

ord.fdr		<- order(fdr)
p.outr		<- p.no.assoc
ord.p		<- order(p.outr)

Fdr.fdr[ord.fdr]	<- cumsum((fdr*f)[ord.fdr]) / cumsum((f)[ord.fdr])
Fdr.p[ord.p]	<- cumsum((fdr*f)[ord.p])   / cumsum((f)[ord.p])  
Fdr.BH[ord.p]	<- p.outr[ord.p]   / (cumsum((f)[ord.p])  / sum(f)) 

Fdr.fdr	<- array(dim=c(500,500),Fdr.fdr)
Fdr.p		<- array(dim=c(500,500),Fdr.p)
Fdr.BH	<- array(dim=c(500,500),Fdr.BH)

contour(x.vec,y.vec,f.0,xlim=c(-2,5),ylim=c(-2,5),col=2,xlab="Z1",ylab="Z2")
contour(x.vec,y.vec,f.1,col=3,add=T)
points(mu.00[1],mu.00[2],col=2,pch=3,lwd=2)
points(mu.01[1],mu.01[2],col=3,pch=3,lwd=2)
points(mu.10[1],mu.10[2],col=3,pch=3,lwd=2)
points(mu.11[1],mu.11[2],col=3,pch=3,lwd=2)
contour(x.vec,y.vec,Fdr.fdr,levels=0.05,col=6,add=T,lwd = 3)

contour(x.vec,y.vec,Fdr.fdr,levels=0.05,col=6,xlim=c(-2,5),ylim=c(-2,5),lwd=2,xlab="Z1",ylab="Z2")
contour(x.vec,y.vec,Fdr.p,levels=0.05,col=3,add=T,lwd = 2)
contour(x.vec,y.vec,Fdr.BH,levels=0.05,col=4,add=T,lwd = 2)
legend(-1,5,c("Fdr=0.05 Bayes classifier","Fdr=0.05 p-value based classifier",
"BH estimated Fdr=0.05 p-value based classifier"),lty=1,col=c(6,3,4),lwd=2)

sum(f[Fdr.fdr <= 0.05]) / sum(f)
sum(f[Fdr.p <= 0.05]) / sum(f)
sum(f[Fdr.BH <= 0.05]) / sum(f)


# 2.2 Plots for testing no association DIFFERENT HYPERPARAMETER values

pi.00	<- 0.85
pi.01	<- 0.00
pi.10	<- 0.15
pi.11	<- 0.00

f		<- pi.00*f.00 + pi.01*f.01 + pi.10*f.10 + pi.11*f.11
f.0		<- pi.00*f.00 
f.1		<- pi.01*f.01 + pi.10*f.10 + pi.11*f.11
fdr		<- f.0 / f

ord.fdr		<- order(fdr)
p.outr		<- p.no.assoc
ord.p		<- order(p.outr)

Fdr.fdr[ord.fdr]	<- cumsum((fdr*f)[ord.fdr]) / cumsum((f)[ord.fdr])
Fdr.p[ord.p]	<- cumsum((fdr*f)[ord.p])   / cumsum((f)[ord.p])  
Fdr.BH[ord.p]	<- p.outr[ord.p]   / (cumsum((f)[ord.p])  / sum(f)) 

Fdr.fdr	<- array(dim=c(500,500),Fdr.fdr)
Fdr.p	<- array(dim=c(500,500),Fdr.p)
Fdr.BH	<- array(dim=c(500,500),Fdr.BH)

contour(x.vec,y.vec,f.0,xlim=c(-2,5),ylim=c(-2,5),col=2,xlab="Z1",ylab="Z2")
contour(x.vec,y.vec,f.1,col=3,add=T)
points(mu.00[1],mu.00[2],col=2,pch=3,lwd=2)
points(mu.01[1],mu.01[2],col=3,pch=3,lwd=2)
points(mu.10[1],mu.10[2],col=3,pch=3,lwd=2)
points(mu.11[1],mu.11[2],col=3,pch=3,lwd=2)
contour(x.vec,y.vec,Fdr.fdr,levels=0.05,col=6,add=T,lwd = 3)

contour(x.vec,y.vec,Fdr.fdr,levels=0.05,col=6,xlim=c(-2,5),ylim=c(-2,5),lwd=2,xlab="Z1",ylab="Z2")
contour(x.vec,y.vec,Fdr.p,levels=0.05,col=3,add=T,lwd = 2)
contour(x.vec,y.vec,Fdr.BH,levels=0.05,col=4,add=T,lwd = 2)
legend(-1,5,c("Fdr=0.05 Bayes classifier","Fdr=0.05 p-value based classifier",
"BH estimated Fdr=0.05 p-value based classifier"),lty=1,col=c(6,3,4),lwd=2)

sum(f[Fdr.fdr <= 0.05]) / sum(f)
sum(f[Fdr.p <= 0.05]) / sum(f)
sum(f[Fdr.BH <= 0.05]) / sum(f)




# 2.3 Testing no replicability


pi.00	<- 0.85
pi.01	<- 0.05
pi.10	<- 0.05
pi.11	<- 0.05

f		<- pi.00*f.00 + pi.01*f.01 + pi.10*f.10 + pi.11*f.11
f.0		<- pi.00*f.00 + pi.01*f.01 + pi.10*f.10 
f.1		<- pi.11*f.11
fdr		<- f.0 / f

ord.fdr		<- order(fdr)
p.outr		<- p.no.rep
ord.p		<- order(p.outr)

Fdr.fdr[ord.fdr]	<- cumsum((fdr*f)[ord.fdr]) / cumsum((f)[ord.fdr])
Fdr.p[ord.p]	<- cumsum((fdr*f)[ord.p])   / cumsum((f)[ord.p])  
Fdr.BH[ord.p]	<- p.outr[ord.p]   / (cumsum((f)[ord.p])  / sum(f)) 

Fdr.fdr	<- array(dim=c(500,500),Fdr.fdr)
Fdr.p	<- array(dim=c(500,500),Fdr.p)
Fdr.BH	<- array(dim=c(500,500),Fdr.BH)


contour(x.vec,y.vec,f.0,levels=c(.001,.005,seq(.01,.15,length=15)),xlim=c(-2,5),ylim=c(-2,5),col=2,xlab="Z1",ylab="Z2")
contour(x.vec,y.vec,f.1,col=3,add=T)
points(mu.00[1],mu.00[2],col=2,pch=3,lwd=2)
points(mu.01[1],mu.01[2],col=2,pch=3,lwd=2)
points(mu.10[1],mu.10[2],col=2,pch=3,lwd=2)
points(mu.11[1],mu.11[2],col=3,pch=3,lwd=2)
contour(x.vec,y.vec,Fdr.fdr,levels=0.05,col=6,add=T,lwd = 3)

contour(x.vec,y.vec,Fdr.fdr,levels=0.05,col=6,xlim=c(-2,5),ylim=c(-2,5),lwd=2,xlab="Z1",ylab="Z2")
contour(x.vec,y.vec,Fdr.p,levels=0.05,col=3,add=T,lwd = 2)
contour(x.vec,y.vec,Fdr.BH,levels=0.05,col=4,add=T,lwd = 2)
legend(-1,0,c("Fdr=0.05 Bayes classifier","Fdr=0.05 p-value based classifier",
"BH estimated Fdr=0.05 p-value based classifier"),lty=1,col=c(6,3,4),lwd=2)


sum(f[Fdr.fdr <= 0.05]) / sum(f)
sum(f[Fdr.p <= 0.05]) / sum(f)
sum(f[Fdr.BH <= 0.05]) / sum(f)
max(p.outr[Fdr.p <= 0.05])  / (sum((f)[Fdr.p <= 0.05]) / sum(f))




#-------------------------------
#	Analyze 2D data
#-------------------------------


library(repfdr)

m <- 500000

pi.00	<- 0.85
pi.01	<- 0.05
pi.10	<- 0.05
pi.11	<- 0.05

n.vec	<- sample(1:4,m,replace = T,prob=c(pi.00,pi.01,pi.10,pi.11))
table(n.vec)

H.mat	<- array(0,dim=c(m,2))
H.mat[n.vec == 3 | n.vec == 4,1]	<- 1
H.mat[n.vec == 2 | n.vec == 4,2]	<- 1

Z.mat	<- H.mat*3  + cbind(rnorm(m),rnorm(m))

P.mat		<- 1 - pnorm(Z.mat,mean=mu.0)
P.no.assoc	<- 1 - pchisq(-2*log(P.mat[,1]) -2*log(P.mat[,2]),4)
P.no.replic	<- apply(P.mat,1,max)

round(cbind(P.mat,P.no.assoc,P.no.replic)[1:10,],5)


par(mfcol=c(1,1))
plot(Z.mat[1:2000,1],Z.mat[1:2000,2],xlab = "1st study Z-scores",ylab = "2nd study Z-scores")
abline(h = 0,lty=2,lwd = 0.5)
abline(v = 0,lty=2,lwd = 0.5)
abline(h = 3,lty=2,lwd = 0.5)
abline(v = 3,lty=2,lwd = 0.5)

#cbind(n.vec,H.mat,Z.mat,P.mat)

input.to.repfdr <- ztobins(Z.mat, 2)


str(input.to.repfdr)
str(input.to.repfdr$pdf.binned.z)

par(mfcol=c(2,1))
plot(1:119,input.to.repfdr$pdf.binned.z[1,,1],type="l",main = "Study 1 pdf's",xlab = "Bin number")
lines(1:119,input.to.repfdr$pdf.binned.z[1,,2],type="l")
plot(1:119,input.to.repfdr$pdf.binned.z[2,,1],type="l",main = "Study 2 pdf's",xlab = "Bin number")
lines(1:119,input.to.repfdr$pdf.binned.z[2,,2],type="l")

pbz_sim    <- input.to.repfdr$pdf.binned.z
bz_sim     <- input.to.repfdr$binned.z.mat


output.rep <- repfdr(pbz_sim, bz_sim, "meta-analysis")
output.rep$Pi
ord.fdr	<- order(output.rep$mat[,1])
cbind(ord.fdr[1:100],output.rep$mat[ord.fdr[1:100],])

table(n.vec[output.rep$mat[,2] <= 0.05])
sum(n.vec == 1 & output.rep$mat[,2] <= 0.05) / sum(output.rep$mat[,2] <= 0.05) 

table(p.adjust(P.no.assoc,method = "BH") <= 0.05)
table(output.rep$mat[,2] <= 0.05)
table(p.adjust(P.no.assoc,method = "BH") <= 0.05,output.rep$mat[,2] <= 0.05)

sum(n.vec == 1 & p.adjust(P.no.assoc,method = "BH") <= 0.05) / sum(p.adjust(P.no.assoc,method = "BH") <= 0.05)
0.05 * 0.85


output.rep <- repfdr(pbz_sim, bz_sim, "replication")
table(n.vec[output.rep$mat[,2] <= 0.05])
sum(n.vec <= 3 & output.rep$mat[,2] <= 0.05) / sum(output.rep$mat[,2] <= 0.05) 

table(n.vec[p.adjust(P.no.replic,method = "BH") <= 0.05])