R code for Bayesian tail confidence interval for p1 - p2 with confidence level = conflev, using independent beta(a, b) and beta(c, d) priors, based on binomial variates (x1, x2) with sample sizes (n1, n2).

-----------------------------------------------------------------------------
# approximate tail method using simulation
diff.app<- function(a1,b1,c1,d1,conflev,nsim=100000)
{ z1 <- rbeta(nsim, a1,b1)
  z2 <- rbeta(nsim, c1,d1)
  z <- z1 - z2
  z <- sort(z)
  lq <- nsim * (1-conflev)/2
  uq <- nsim * (1 - (1-conflev)/2)
  ci <- array(0,2)
  ci[1] <- z[lq]
  ci[2] <- z[uq]
  return(ci) }
# Bayes tail interval with beta priors
fct.F1<- function(x,t,a1,b1,c1,d1){dbeta(x,c1,d1)*pbeta(x+t,a1,b1)}
fct.F2<- function(x,t,a1,b1,c1,d1){dbeta(x,c1,d1)*(1-pbeta(x+t,a1,b1))}
diff.F <- function(t,a1,b1,c1,d1)
{ if(t < 0)
  Fvalue <- integrate(fct.F1,-t,1,t=t,a1=a1,b1=b1,c1=c1,d1=d1)$value
  else
  Fvalue <- 1-integrate(fct.F2,0,1-t,t=t,a1=a1,b1=b1,c1=c1,d1=d1)$value
  return(Fvalue) }
diff.fct <- function(ab,a1,b1,c1,d1,conflev)
{ abs(diff.F(ab[2],a1,b1,c1,d1) - (1 - (1-conflev)/2))+
  abs(diff.F(ab[1],a1,b1,c1,d1) - (1-conflev)/2) }
diffCI <- function(x1,n1,x2,n2,a,b,c,d,conflev=.95)
{ a1 <- a + x1
  b1 <- b + n1 - x1
  c1 <- c + x2
  d1 <- d + n2 - x2
  start <- diff.app(a1,b1,c1,d1,conflev)
  tailci <- optim(start,diff.fct,a1=a1,b1=b1,c1=c1,d1=d1,
          conflev=conflev,control=list(maxit=20000))$par
  if(tailci[1] < -1) tailci[1]  <- -1
  if(tailci[2] >  1) tailci[2]  <- 1
  return(tailci) }
-----------------------------------------------------------------------------

-----------------------------------------------------------------------------
# R code for Bayesian HPD confidence interval for p_1 - p_2
# with confidence level = conflev, using independent beta(a, b) and
# beta(c, d) priors, based on
# binomial variates (x_1, x_2) with sample sizes (n_1, n_2).

# approximate HPD confidence interval
diff.app<- function(a1,b1,c1,d1,conflev)
{
  mu1 <- a1/(a1+b1)
  mu2 <- c1/(c1+d1)
  v1 <- mu1*(1-mu1)/(a1 + b1)
  v2 <- mu2*(1-mu2)/(c1 + d1)
  z  <- qnorm(1-(1-conflev)/2)
  ci <- array(0,2)
  ci[1] <- mu1-mu2 - z*sqrt(v1+v2)
  ci[2] <- mu1-mu2 + z*sqrt(v1+v2)
  return(ci)
}

# Bayes HPD interval with beta priors
fct.F1<- function(x,t,a1,b1,c1,d1){
  dbeta(x,c1,d1)*pbeta(x+t,a1,b1)}

fct.F2<- function(x,t,a1,b1,c1,d1){
  dbeta(x,c1,d1)*(1-pbeta(x+t,a1,b1))}

diff.F <- function(t,a1,b1,c1,d1)
{ if(t < 0)
  Fvalue <- integrate(fct.F1,-t,1,t=t,a1=a1,b1=b1,c1=c1,d1=d1)$value
  else
  Fvalue <- 1-integrate(fct.F2,0,1-t,t=t,a1=a1,b1=b1,c1=c1,d1=d1)$value
  return(Fvalue)
}

fct.f<- function(x,t,a1,b1,c1,d1){
 dbeta(x,c1,d1)*dbeta(x+t,a1,b1)
}

diff.f <- function(t,a1,b1,c1,d1)
{
  if(t < -1) fvalue <- 100
  else if(t > 1)  fvalue <- 100
  else if((t >= -1) && (t <= 0))
    fvalue <- integrate(fct.f,-t,1,t=t,a1=a1,b1=b1,c1=c1,d1=d1)$value
  else
    fvalue <- integrate(fct.f,0,1-t,t=t,a1=a1,b1=b1,c1=c1,d1=d1)$value
  return(fvalue)
}


diff <- function(ab,a1,b1,c1,d1,conflev)
{
  1000*abs(diff.F(ab[2],a1,b1,c1,d1) -
      diff.F(ab[1],a1,b1,c1,d1) - conflev)+
  abs(diff.f(ab[1],a1,b1,c1,d1) -
        diff.f(ab[2],a1,b1,c1,d1))
}

diffCIhpd <- function(x1,n1,x2,n2,a,b,c,d,conflev=.95)
{
  y <- x1
  if( y > n1/2 ) {
   x2 <- n2-x2
   x1 <- n1-x1
  }
  a1 <- a + x1
  b1 <- b + n1 - x1
  c1 <- c + x2
  d1 <- d + n2 - x2
  start <- diff.app(a1,b1,c1,d1,conflev)
  hdrci <- optim(start,diff,a1=a1,b1=b1,c1=c1,d1=d1,conflev=conflev)$par
  if(hdrci[1] < -1) hdrci[1]  <- -1
  if(hdrci[2] >  1) hdrci[1]  <- 1
   if( y > n1/2 ) {
   ci <- hdrci
   hdrci[1] <- -ci[2]
   hdrci[2] <- -ci[1]
  }
  return(hdrci)
}

R code for Bayesian tail confidence interval for relative risk p1/p2 with confidence level = conflev, using independent beta(a, b) and beta(c, d) priors, based on binomial variates (x1, x2) with sample sizes (n1, n2).

-----------------------------------------------------------------------------
# riskCI generates a Bayesian tail confidence interval for p_1/ p_2
# with confidence level = conflev, using independent beta(a, b) and
# beta(c, d) priors, based on
# binomial variates (x1, x2) with sample sizes (n1, n2).


# approximate tail method using simulation
risk.app<- function(a1,b1,c1,d1,conflev,nsim=100000)
{
  z1 <- rbeta(nsim, a1,b1)
  z2 <- rbeta(nsim, c1,d1)
  z <- z1/z2
  z <- sort(z)
  lq <- nsim * (1-conflev)/2
  uq <- nsim * (1 - (1-conflev)/2)
  ci <- array(0,2)
  ci[1] <- z[lq]
  ci[2] <- z[uq]
  return(ci)
}

# Bayes tail interval with beta priors

fct.F1<- function(x,t,a1,b1,a2,b2){
 dbeta(x,a2,b2)*pbeta(x*t,a1,b1)}

fct.F2<- function(x,t,a1,b1,a2,b2){
 dbeta(x,a2,b2)*(1-pbeta(x*t,a1,b1))}

risk.F <- function(t,a1,b1,a2,b2)
{
  if((0<t) && (t<=1)){
    return(integrate(fct.F1,0,1,t=t,a1=a1,b1=b1,a2=a2,b2=b2)$value)
  }
  else{
    return(1-integrate(fct.F2,0,1/t,t=t,a1=a1,b1=b1,a2=a2,b2=b2)$value)
  }
}

risk.fct <- function(ab,a1,b1,c1,d1,conflev)
{
  abs(risk.F(ab[2],a1,b1,c1,d1) - (1 - (1-conflev)/2)) +
  abs(risk.F(ab[1],a1,b1,c1,d1) -(1-conflev)/2)
}

riskCI <- function(x1,n1,x2,n2,a,b,c,d,conflev=.95)
{
  a1 <- a + x1
  b1 <- b + n1 - x1
  c1 <- c + x2
  d1 <- d + n2 - x2
  start <- risk.app(a1,b1,c1,d1,conflev)
  tailci <- optim(start,risk.fct,a1=a1,b1=b1,c1=c1,d1=d1,
            conflev=conflev,control=list(maxit=20000))$par
  if(tailci[1] < 0) tailci[1]  <- 0
  return(tailci)
}
-----------------------------------------------------------------------------

R code for Bayesian tail confidence interval for odds ratio with confidence level = conflev, using independent beta(a, b) and beta(c, d) priors, based on binomial variates (x1, x2) with sample sizes (n1, n2).

-----------------------------------------------------------------------------
# orCI generates a Bayesian tail confidence interval for
# p1(1-p2)/(p2(1-p1))  with confidence level = conflev, using
# independent beta(a, b) and  beta(c, d) priors, based on
# binomial variates (x1, x2) with sample sizes (n1, n2).


# approximate tail method using simulation

or.app<- function(a1,b1,c1,d1,conflev,nsim=1000000)
{
  z1 <- rf(nsim, 2*a1,2*b1)
  z2 <- rf(nsim, 2*c1,2*d1)
  a <- (d1/c1)/(b1/a1)
  z <- a*z1/z2
  z <- sort(z)
  lq <- nsim * (1-conflev)/2
  uq <- nsim * (1 - (1-conflev)/2)
  ci <- array(0,2)
  ci[1] <- z[lq]
  ci[2] <- z[uq]
  return(ci)
}


# Bayes tail interval with beta priors

fct.F<- function(x,t,a1,b1,a2,b2){
  c <- (b2/a2)/(b1/a1)
  df(x,2*a2,2*b2)*pf(x*t/c,2*a1,2*b1)
}

or.F <- function(t,a1,b1,a2,b2)
{
  return(integrate(fct.F,0,Inf,t=t,a1=a1,b1=b1,a2=a2,b2=b2)$value)
}

or.fct <- function(ab,a1,b1,c1,d1,conflev)
{
 abs(or.F(ab[2],a1,b1,c1,d1) - (1 - (1-conflev)/2))+
 abs(or.F(ab[1],a1,b1,c1,d1) - (1-conflev)/2)
}

orCI <- function(x1,n1,x2,n2,a,b,c,d,conflev=.95)
{
  if(x2!=n2){
  a1 <- a + x1
  b1 <- b + n1 - x1
  c1 <- c + x2
  d1 <- d + n2 - x2
  start <- or.app(a1,b1,c1,d1,conflev)
  tailci <- optim(start,or.fct,a1=a1,b1=b1,c1=c1,d1=d1,
           conflev=conflev,control=list(maxit=20000))$par
  if(tailci[1] < 0) tailci[1]  <- 0 }
  else{
  a1 <- a + n1 - x1
  b1 <- b +  x1
  c1 <- c + n2 - x2
  d1 <- d + x2
  start <- or.app(a1,b1,c1,d1,conflev)
  tailci1 <- optim(start,or.fct,a1=a1,b1=b1,c1=c1,d1=d1,
           conflev=conflev,control=list(maxit=20000))$par
  if(tailci[1] < 0) tailci[1]  <- 0
  tailci <- array(0,2)
  tailci[1] <- 1/ tailci1[2]
  tailci[2] <- 1/ tailci1[1]
 }
  return(tailci)
}
-----------------------------------------------------------------------------