next up previous
Next: About this document ...

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 $z^2/2$ successes
and $z^2/2$ 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)
  }
-----------------------------------------------------------------------------




next up previous
Next: About this document ...
Alan Agresti 2006-03-06