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))
}

diffCI <- 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)
}
-----------------------------------------------------------------------------
```

Alan Agresti 2003-12-09