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