[Rd] fisher.test() gives wrong confidence interval (PR#4019)

[Kurt Hornik]

>>>>> jerome  writes:

> The problem occurs when the sample odds ratio is Inf, such as in the 
> following example. Given the fact that both upper bounds of the two 95% 
> confidence intervals are Inf, I would have expected that the two lower 
> bounds be equal, but they aren't.

> x <- matrix(c(9,4,0,2),2,2)
> x
> #     [,1] [,2]
> #[1,]    9    0
> #[2,]    4    2
> rbind("two.sided.95CI"=fisher.test(x)$conf.int,
> "greater.95CI"=fisher.test(x,alt="greater")$conf.int)
> #                    [,1] [,2]
> #two.sided.95CI 0.2985103  Inf
> #greater.95CI   0.4625314  Inf

> Using the noncentral hypergeometric distribution, we can calculate the
> probability mass of each possible table with same marginals as x.
> Ref.: Alan Agresti (1990). Categorical data analysis. New York: Wiley.
> Page 67.

> Hence, the result below suggests that the two-sided confidence interval 
> has a confidence level of 97.5% as opposed to 95%.

> n11 <- 7:9
> theta <- 0.2985103
> choose(9,n11)*choose(15-9,13-n11)*theta^n11/
>   sum(choose(9,n11)*choose(15-9,13-n11)*theta^n11)
> #[1] 0.67344877 0.30154709 0.02500414

> The 95% confidence interval with one-sided (greater) alternative appears 
> to be correct.

> theta <- 0.4625314
> choose(9,n11)*choose(15-9,13-n11)*theta^n11/
>   sum(choose(9,n11)*choose(15-9,13-n11)*theta^n11)
> #[1] 0.5608724 0.3891316 0.0499960

Right.  When going for a confidence interval that is to include the
parameter estimate and this is already at the boundary of the parameter
range, then there is no point finding mass even further away as there
obviously cannot be any.

I changed the code (r-devel) to

        CINT <- switch(alternative,
                       less = c(0, ncp.U(x, 1 - conf.level)),
                       greater = c(ncp.L(x, 1 - conf.level), Inf),
                       two.sided = {
                           if(ESTIMATE == 0)
                               c(0, ncp.U(x, 1 - conf.level))
                           else if(ESTIMATE == Inf)
                               c(ncp.L(x, 1 - conf.level), Inf)
                           else {
                               alpha <- (1 - conf.level) / 2
                               c(ncp.L(x, alpha), ncp.U(x, alpha))
                           }
                       })

which seems to fix the problem.

Thanks,

-k