binom.blaker.acc {BlakerCI}

Blaker's binomial acceptability function, optionally unimodalized.

Description

Calculates values of the acceptability function for the binomial distribution (see function acceptbin in Blaker (2000)) in a sequence of points (for, e.g., plotting purposes). The acceptability function may optionally be “unimodalized”, i.e. replaced with the smallest greater or equal unimodal function.

Usage

binom.blaker.acc(x, n, p, type = c("orig", "unimod"),
      acc.tol = 1e-10, ...)

Arguments

x

number of successes.

n

number of trials.

p

vector (length 1 allowed) of hypothesized binomial parameters (between 0 and 1). In case of more than one point, an increasing sequence required.

type

for type = "orig", original acceptability function calculated.
For type = "unimod", smallest unimodal function greater or equal to the acceptability function calculated instead.

acc.tol

numerical tolerance (relevant only for type = "unimod").

...

additional arguments to be passed to binom.blaker.acc.single.p; in fact, just maxiter (see BlakerCI-internal).

Details

For type = "orig", essentially the same is calculated as – for single points – by acceptbin function from Blaker (2000).

Single values of the “unimodalized” acceptability function (for type = "unimod") are computed by an iterative numerical algorithm implemented in internal function
binom.blaker.acc.single.p. The function cited is called just once in each of the intervals where the acceptability function is continuous (namely in the leftmost one of those points of p that fall into the interval when dealing with points below x/n, and the rightmost one when above x/n). The rest is done by function cummax. This is considerably faster than calling binom.blaker.acc.single.p for every point of p. Note that applying cummax directly to a vector of unmodified acceptability values is even faster and provides a unimodal output; it may, nevertheless, lack accuracy (see Examples).

Value

Vector of acceptability values (with or without unimodalization) in points of p.

Note

Inspired by M.P. Fay (2010), mentioning “unavoidable inconsistencies” between tests with non-unimodal acceptability functions and confidence intervals derived from them. When the acceptability functions are unimodalized and the test modified accordingly (i.e. p-values slightly increased in some cases), a perfectly matching test-CI pair is obtained.

Author(s)

Jan Klaschka [e-mail]

References

Blaker, H. (2000) Confidence curves and improved exact confidence intervals for discrete distributions. Canadian Journal of Statistics 28: 783-798.
(Corrigenda: Canadian Journal of Statistics 29: 681.)

Fay, M.P. (2010). Two-sided Exact Tests and Matching Confidence Intervals for Discrete Data. R Journal 2(1): 53-58.

Examples

p <- seq(0,1,length=1001)
acc <- binom.blaker.acc(3,10,p)
acc1 <- binom.blaker.acc(3,10,p,type="unimod")
## The two functions look the same at first glance.
plot(p,acc,type="l")
lines(p,acc1,col="red")
legend(x=.7,y=.8,c("orig","unimod"),col=c("black","red"),lwd=1)

## There is, nevertheless, a difference.
plot(p,acc1-acc,type="l")

## Focussing on the difference about p=0.4:
p <- seq(.395,.405,length=1001)
acc <- binom.blaker.acc(3,10,p)
acc1 <- binom.blaker.acc(3,10,p,type="unimod")
plot(p,acc,type="l",ylim=c(.749,.7495))
lines(p,acc1,col="red")
legend(x=.402,y=.7494,c("orig","unimod"),col=c("black","red"),lwd=1)

## Difference between type="unimod" and mere applying
## cummax to values obtained via type="orig":
p <- seq(0,1,length=1001)
x <- 59
n <- 355
## Upper confidence limit (at 0.95 level) is slightly above 0.209:
binom.blaker.limits(x,n) ## [1] 0.1300807 0.2090809
## Unmodified acceptability value fall below 0.05 at p = .209
## left to the limit (so that the null hypothesis p = .209
## would be rejected despite the fact that p lies within
## the confidence interval):
acc <- binom.blaker.acc(59,355,p)
rbind(p,acc)[,207:211]
##           [,1]       [,2]       [,3]       [,4]       [,5]
## p   0.20600000 0.20700000 0.20800000 0.20900000 0.21000000
## acc 0.06606867 0.05759836 0.05014189 0.04999082 0.04330283
##
## Modified acceptability is above 0.05 at p = 0.05 (so that
## hypothesis p = 0.05 is not rejected by the modified test):
acc1 <- binom.blaker.acc(59,355,p,type="unimod")
rbind(p,acc1)[,207:211]
##            [,1]       [,2]       [,3]       [,4]       [,5]
## p    0.20600000 0.20700000 0.20800000 0.20900000 0.21000000
## acc1 0.06608755 0.05759836 0.05014189 0.05000009 0.04331409
##
## Applying cummax to unmodified acceptabilities guarantees unimodality 
## but lacks accuracy, leaving the value at p = 0.209 below 0.05:
m <- max(which(p <= 59/355))
acc2 <- acc
acc2[1:m] <- cummax(acc2[1:m])
acc2[1001:(m+1)] <- cummax(acc2[1001:(m+1)])
rbind(p,acc2)[,207:211]
##            [,1]       [,2]       [,3]       [,4]       [,5]
## p    0.20600000 0.20700000 0.20800000 0.20900000 0.21000000
## acc2 0.06606867 0.05759836 0.05014189 0.04999082 0.04330283