I am stuck with the same problem. I searched and searched in stackExchange and Google, and I did not find anything explicitly explaining how the odds ratio in fisher.test()
is calculated. This is not a complete answer, but it provides something useful to the discussion.
One thing is clear, fisher.test()
calculates a completely different odds ratio with a different (and more complicated) method.
If you write fisher.test
in R Console and press Enter, the output will be the complete code of the function fisher.test()
. The interesting part concerning our doubt is:
m <- sum(x[, 1L])
n <- sum(x[, 2L])
k <- sum(x[1L, ])
x <- x[1L, 1L]
lo <- max(0L, k - n)
hi <- min(k, m)
support <- lo:hi #Interval of definition of Hypergeometric Distribution
logdc <- dhyper(support, m, n, k, log = TRUE) #log of Hypergeometric Probability Function
dnhyper <- function(ncp) {
d <- logdc + log(ncp) * support
d <- exp(d - max(d))
d/sum(d)
}
mnhyper <- function(ncp) {
if (ncp == 0)
return(lo)
if (ncp == Inf)
return(hi)
sum(support * dnhyper(ncp))
}
mle <- function(x) {
if (x == lo)
return(0)
if (x == hi)
return(Inf)
mu <- mnhyper(1)
if (mu > x)
uniroot(function(t) mnhyper(t) - x, c(0, 1))$root
else if (mu < x)
1/uniroot(function(t) mnhyper(1/t) - x, c(.Machine$double.eps,
1))$root
else 1
}
Well, if you
- Define
x = matrix(c(3, 6, 5, 6), nrow=2)
- Run every line of the code extracted from
fisher.test
in order
- Run
mle(x)
You get 0.6155891
. So this is how the odds ratio in fisher.test
is calculated. However, I do not understand the algorithm. I have not found any article explaining it, and it does not seems possible to translate it to a "simple" mathematical formula.
EDIT:
The odds ratio estimate obtained by fisher.test
is the conditional maximum-likelihood estimate (CLME) of the odds ratio.
Given a table like this (where $C1$ is the sum of the first column, and so on):
\begin{array}{ll|l}
a_{11} & a_{12} & R_1 \\
a_{12} & a_{22} & R_2 \\ \hline
C_1 & C_2 & n
\end{array}
The CLME odds ratio is defined as
$$ \text{arg } \max_{\psi >0} \frac{\binom{C_1}{a_{11}}\binom{C_2}{a_{12}}\psi^{a_{11}}}{\sum_{k = \max(0, R1 - C2)}^{min(R1, C1)} \binom{C_1}{k}\binom{C_2}{R_1-k}\psi^k }$$
There is no explicit formula of the CLME odds ratio, so one has to use iterative methods to find the argument of the maximum of the function. I suppose that is what fisher.test()
does, but I don't recognize the algorithm. So I programmed one by myself:
OR <- function(x) {
C1 <- sum(x[, 1])
C2 <- sum(x[, 2])
R1 <- sum(x[1, ])
a11 <- x[1, 1]
a12 <- x[1, 2]
n <- sum(x)
fobjetivo <- function(psi) {
lo <- max(0, R1 - C2)
hi <- min(R1, C1)
sumandosDenom <- numeric(n)
for (k in lo:hi) {
sumandosDenom[k] <- choose(C1,k) * choose(C2, R1-k) * psi^k
}
salida <- ( choose(C1, a11) * choose(C2, a12) * psi^a11 ) / sum( sumandosDenom )
return(salida)
}
as.numeric(optimize(fobjetivo, interval = c(0.01, 100), maximum = T, tol = 0.00001)$maximum)
}
x = matrix(c(3, 6, 5, 6), nrow=2)
OR(x)
The output is 0.6052843
. Still not the same value obtained by the fisher.test()
, but in general it appears to be closer than the sample odds ratio. I suppose the difference is explained by the different methods, by I don't know for sure.
Reference: Kenneth J. Rothman, Sander Greenland and Timothy L. Lash (2008): Modern Epidemiology, 3rd Edition, Lippincott-Raven Publishers, p. 257