I believe that the resulting interval will sometimes fail to satisfy the
confidence statement for values of alpha other than 0 or Inf.
alpha?
('sometimes' rather than always, because of the discreteness of the
distribution)
I think you need
if(ESTIMATE == 0)
c(0, ncp.U(x, 1 - conf.level / 2))
else if(ESTIMATE == Inf)
c(ncp.L(x, 1 - conf.level / 2 ), Inf)
else {
alpha <- (1 - conf.level) / 2
c(ncp.L(x, alpha), ncp.U(x, alpha))
}
I don't think this is right either. The only "correct" fix is to invert
the two-sided *test*. This may or may not have to include mass in the
opposite tails, so any attempt to base the logic on the value of
ESTIMATE is going to fail. Unfortunately, this inversion is tricky
because adjusting the parameter causes points to enter and leave the
critical region in a discontinuous manner, or, put otherwise, the
p-value is an awkward function of the parameters.
For your edification consider the following
binom.test(3,15)
binom.test(3,15,p=0.04331201) # p = 0.025
binom.test(3,15,p=0.48089113) # p > 0.025, but < 0.05
# Define B(x) = two-sided p value for param == x
# Needs to be vectorized so that curve works
B <- function(x) sapply(x,function(x)binom.test(3,15,p=x)$p.value)
curve(B,from=0,to=1,n=5000)
uniroot(function(x)B(x)-.05,c(.2,1))
uniroot(function(x)B(x)-.05,c(0,.2))
In this case, it actually works and gives an improved CI from 0.05685
to 0.4657785, but if you look closely at the curve of B, you'll see
that it is discontinuous, and in fact not even monotone. There is, by
the way, a particularly nasty case with binom.test(0,5) where a
discontinuity hits exactly at p=0.5, which is also the endpoint of the
CI.
I think that the real issue here is whether we should take our chances
and implement the uniroot technique, knowing that it might fail in
some (probably few) cases.
O__ ---- Peter Dalgaard Blegdamsvej 3
c/ /'_ --- Dept. of Biostatistics 2200 Cph. N
(*) \(*) -- University of Copenhagen Denmark Ph: (+45) 35327918
~~~~~~~~~~ - (p.dalgaard@biostat.ku.dk) FAX: (+45) 35327907