[Rd] fisher.test() gives wrong confidence interval (PR#4019)
Kurt Hornik
hornik at ci.tuwien.ac.at
Sat Aug 30 13:00:16 MEST 2003
>>>>> 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
More information about the R-devel
mailing list