R Code for Wilson's confidence interval for a proportion
(Journal of the American Statistical
Association 1927). This is the score CI, based on inverting the
asymptotic normal test
using the null standard error.
----------------------------------------------------------------------------- scoreci <- function(x,n,conflev) { zalpha <- abs(qnorm((1-conflev)/2)) phat <- x/n bound <- (zalpha*((phat*(1-phat)+(zalpha**2)/(4*n))/n)**(1/2))/ (1+(zalpha**2)/n) midpnt <- (phat+(zalpha**2)/(2*n))/(1+(zalpha**2)/n) uplim <- round(midpnt + bound,digits=4) lowlim <- round(midpnt - bound,digits=4) results <- data.frame(lowlim,uplim) cat("\n") cat("With confidence level",conflev," and sample proportion", round(phat,digits=4), " \nthe lower and upper limits for the score confidence interval are: \n") cat("\n") print(results) cat("\n") # This function computes a confidence interval for a proportion. # It is based on inverting the large-sample normal score test for the # proportion. The input parameters are the number of successes, x, the # total sample size, n, and the confidence level, conflev (e.g. 0.90). # Returned are the endpoints for the 100*conflev % confidence # interval for the proportion. # binconf(x,n,conflev) } -----------------------------------------------------------------------------
R Code for Blaker's `exact' CI for a binomial proportion
(Canadian J. Statist., 2000).
Based on code presented by Blaker, and later correction by Blaker in same
journal.
----------------------------------------------------------------------------- blakerci <- function(x,n,level,tolerance=1e-05){ lower = 0 upper = 1 if (x!=0){lower = qbeta((1-level)/2, x, n-x+1) while (acceptbin(x, n, lower + tolerance) < (1 - level)) lower = lower+tolerance } if (x!=n){upper = qbeta(1 - (1-level)/2, x+1, n-x) while (acceptbin(x, n, upper - tolerance) < (1 - level)) upper = upper-tolerance } c(lower,upper) } # Computes the Blaker exact ci (Canadian J. Stat 2000) # for a binomial success probability # for x successes out of n trials with # confidence coefficient = level; uses acceptbin function acceptbin = function(x, n, p){ #computes the Blaker acceptability of p when x is observed # and X is bin(n, p) p1 = 1 - pbinom(x - 1, n, p) p2 = pbinom(x, n, p) a1 = p1 + pbinom(qbinom(p1, n, p) - 1, n, p) a2 = p2 + 1 - pbinom(qbinom(1 - p2, n, p), n, p) return(min(a1,a2)) } -----------------------------------------------------------------------------
R Code for Clopper-Pearson `exact' CI for a binomial
parameter.
This inverts two binomial tests.
----------------------------------------------------------------------------- exactci <- function(x,n,conflev){ alpha <- (1 - conflev) if (x == 0) { ll <- 0 ul <- 1 - (alpha/2)^(1/n) } else if (x == n) { ll <- (alpha/2)^(1/n) ul <- 1 } else { ll <- 1/(1 + (n - x + 1) / (x * qf(alpha/2, 2 * x, 2 * (n-x+1)))) ul <- 1/(1 + (n - x) / ((x + 1) * qf(1-alpha/2, 2 * (x+1), 2 * (n-x)))) } c(ll,ul) } # Computes the Clopper/Pearon exact ci # for a binomial success probability # for x successes out of n trials with # confidence coefficient conflev -----------------------------------------------------------------------------
R Code for Agresti-Coull CI for a binomial proportion, based on adding
successes
and failures before computing the Wald CI (American
Statistician 1998).
The CI is truncated when it overshoots the boundary.
----------------------------------------------------------------------------- addz2ci <- function(x,n,conflev){ z = abs(qnorm((1-conflev)/2)) tr = z^2 #the number of trials added suc = tr/2 #the number of successes added ptilde = (x+suc)/(n+tr) stderr = sqrt(ptilde * (1-ptilde)/(n+tr)) ul = ptilde + z * stderr ll = ptilde - z * stderr if(ll < 0) ll = 0 if(ul > 1) ul = 1 c(ll,ul) } # Computes the Agresti-Coull CI for x successes out of n trials # with confidence coefficient conflev. -----------------------------------------------------------------------------
R Code for Agresti-Coull add-4 CI for a binomial proportion,
based on adding 2 successes and
2 failures before computing the
Wald CI (American Statistician 1998; see also article by
Agresti
and Caffo in American Statistician 2000). The CI is
truncated when it overshoots the boundary.
----------------------------------------------------------------------------- add4ci <- function(x,n,conflev){ ptilde = (x+2)/(n+4) z = abs(qnorm((1-conflev)/2)) stderr = sqrt(ptilde * (1-ptilde)/(n+4)) ul = ptilde + z * stderr ll = ptilde - z * stderr if(ll < 0) ll = 0 if(ul > 1) ul = 1 c(ll,ul) } # Computes the Agresti-Coull `add 4' CI for x successes out of n trials # with confidence coefficient conflev. Adds 2 successes and # 4 trials. -----------------------------------------------------------------------------
R Code written by Anna Gottard of the Univ. of Firenze for the mid-P confidence interval adaptation of the Clopper-Pearson interval, with confidence coefficient (1 - alpha). See, e.g., comments by Agresti and Gottard in Statistical Science, vol. 20, November 2005, pp. 367-371.
----------------------------------------------------------------------------- midPci <- function(x,n,alpha){ pp<-seq(0.0001, 1 , 0.0005) uplim<-1 lowlim<-0 if (x == 0) uplim <- 1-alpha^(1/n) if (x == n) lowlim <- (alpha)^(1/n) if (x>0 & x<n){ a2 <- 0.5*pbinom(x-1, n , pp,lower.tail = T) + 0.5*pbinom(x, n , pp, lower.tail = T, log.p = FALSE) uplim=pp[ max(which(a2>(alpha/2))) ] lowlim=pp[ min(which(a2<(1-alpha/2))) ] } c(lowlim,uplim) } -----------------------------------------------------------------------------