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