#
# Functions written by Jon Wakefield 17th April, 2007.
#
# For more details see "A Bayesian Measure of the Probability of False
# Discovery in Genetic Epidemiology Studies", under revision for American
# Journal of Human Genetics.
#
# This function calculates BFDP, the approximate Pr( H0 | thetahat ),
# given an estiamte of the log relative risk, thetahat, the variance of
# this estimate, V, the prior variance, W, and the prior probability of
# a non-null association.
#
BFDPfunV <- function(thetahat,V,W,pi1){
  pH0 <- dnorm(thetahat,m=0,s=sqrt(V))
  postvar <- V + W
  pH1 <- dnorm(thetahat,m=0,s=sqrt(postvar))
  BF <- pH0/pH1; PO <- (1-pi1)/pi1
  BFDP <- BF*PO/(BF*PO+1)
  list(BF=BF,pH0=pH0,pH1=pH1,BFDP=BFDP)
}
#
# This function calculates BFDP, the approximate Pr( H0 | thetahat ),
# given an estiamte of the relative risk, RRhat, and the upper point of
# a 95% confidence interval RRhi, the prior variance, W, and the prior
# probability of a non-null association.
#
BFDPfunCI <- function(RRhat,RRhi,W,pi1){
  thetahat <- log(RRhat)
  V <- ((log(RRhi)-log(RRhat))/1.96)^2
  pH0 <- dnorm(thetahat,m=0,s=sqrt(V))
  postvar <- V + W
  pH1 <- dnorm(thetahat,m=0,s=sqrt(postvar))
  BF <- pH0/pH1; pi0 <- 1-pi1; PO <- pi0/(1-pi0)
  BFDP <- BF*PO/(BF*PO+1)
  list(BF=BF,pH0=pH0,pH1=pH1,BFDP=BFDP)
}
BFDPfunV2 <- function(thetahat1,thetahat2,V1,V2,W,pi1){
  BF1 <-dnorm(thetahat1,m=0,s=sqrt(V1))/dnorm(thetahat1,m=0,s=sqrt(V1+W))
  BF2 <-dnorm(thetahat2,m=0,s=sqrt(V2))/dnorm(thetahat2,m=0,s=sqrt(V2+W))
  Rfac <- W/(V1*W + V2*W + V1*V2)
  term1 <- thetahat1^2/V1^2
  term2 <- 2*thetahat1*thetahat2/(V1*V2)
  term3 <- thetahat2^2/V2^2
  BF12 <- sqrt(W/(Rfac*V1*V2))*exp(-.5*(term1+term2+term3)*V1*V2*Rfac)
  BF2given1 <- BF12/BF1
  PO <- (1-pi1)/pi1
  BFDP1 <- BF1*PO/(BF1*PO+1)
  BFDP2 <- BF2*PO/(BF2*PO+1)
  BFDP12 <- BF12*PO/(BF12*PO+1)
  smallr1 <- W/(V1+W)
  smallr2 <- W/(V2+W)
# pm1 and pv1 are the posterior median and posterior variance from the first study only
  pm1 <- smallr1*thetahat1; pv1 <- smallr1*V1
# pm2 and pv2 are the posterior median and posterior variance from the second study only
  pm2 <- smallr2*thetahat2; pv2 <- smallr2*V2
# pm12 and pv12 are the posterior median and posterior variance from the first two studies combined
  pm12 <- thetahat1*V2*Rfac + thetahat2*V1*Rfac; pv12 <- Rfac*V1*V2
  list(BF1=BF1,BF2=BF2,BF12=BF12,BF2given1=BF2given1,BFDP1=BFDP1,BFDP2=BFDP2,
       pm1=pm1,pv1=pv1,pm2=pm2,pv2=pv2,pm12=pm12,pv12=pv12)
}
BFDPfunCI2 <- function(RRhat1,RRhat2,RRhi1,RRhi2,W,pi1){
  thetahat1 <- log(RRhat1); thetahat2 <- log(RRhat2) 
  V1 <- ((log(RRhi1)-log(RRhat1))/1.96)^2
  V2 <- ((log(RRhi2)-log(RRhat2))/1.96)^2
  BF1 <-dnorm(thetahat1,m=0,s=sqrt(V1))/dnorm(thetahat1,m=0,s=sqrt(V1+W))
  BF2 <-dnorm(thetahat2,m=0,s=sqrt(V2))/dnorm(thetahat2,m=0,s=sqrt(V2+W))
  Rfac <- W/(V1*W + V2*W + V1*V2)
  term1 <- thetahat1^2/V1^2
  term2 <- 2*thetahat1*thetahat1/(V1*V2)
  term3 <- thetahat2^2/V2^2
  BF12 <- sqrt(W/(Rfac*V1*V2))*exp(-.5*(term1+term2+term3)*V1*V2*Rfac)
  BF2given1 <- BF12/BF1
  PO <- (1-pi1)/pi1
  BFDP1 <- BF1*PO/(BF1*PO+1)
  BFDP2 <- BF2*PO/(BF2*PO+1)
  BFDP12 <- BF12*PO/(BF12*PO+1)
  smallr1 <- W/(V1+W)
  smallr2 <- W/(V2+W)
# pm1 and pv1 are the posterior median and posterior variance from the first study only
  pm1 <- smallr1*thetahat1; pv1 <- smallr1*V1
# pm2 and pv2 are the posterior median and posterior variance from the second study only
  pm2 <- smallr2*thetahat2; pv2 <- smallr2*V2
# pm12 and pv12 are the posterior median and posterior variance from the first two studies combined
  pm12 <- thetahat1*V2*Rfac + thetahat2*V1*Rfac; pv12 <- Rfac*V1*V2
  list(BF1=BF1,BF2=BF2,BF12=BF12,BF2given1=BF2given1,BFDP1=BFDP1,BFDP2=BFDP2,BFDP12=BFDP12,
       pm1=pm1,pv1=pv1,pm2=pm2,pv2=pv2,pm12=pm12,pv12=pv12,V1=V1,V2=V2)
}