Why do odds ratios from formula and R's fisher.test differ? Which one should one choose?

Question

In the following example

> m = matrix(c(3, 6, 5, 6), nrow=2)
> m
     [,1] [,2]
[1,]    3    5
[2,]    6    6
> (OR = (3/6)/(5/6))    #1
[1] 0.6
> fisher.test(m)        #2

    Fisher's Exact Test for Count Data

data:  m 
p-value = 0.6699
alternative hypothesis: true odds ratio is not equal to 1 
95 percent confidence interval:
 0.06390055 5.07793271 
sample estimates:
odds ratio 
 0.6155891 

I calculated the odds ratio (#1) "manually", 0.600; then (#2) as one of the outputs of the Fisher's exact test, 0.616.

Why didn't I get the same value?

Why do several ways of computing the odds-ratio exist, and how to choose the most appropriate one?

Answer 1

From the help page for fisher.test():

Note that the conditional Maximum Likelihood Estimate (MLE) rather than the unconditional MLE (the sample odds ratio) is used.

Answer 2

To add to the discussion here, it is useful to ask what exactly is conditioned on in this "conditional" likelihood. The Fisher test differs from other categorical analyses in that it considers all margins of the table to be fixed whereas the logistic regression model (and corresponding Pearson chi-square test which is the score test of the logistic model) only consider one margin to be fixed.

The Fisher test then considers the hypergeometric distribution as a probability model for the counts observed in each of the 4 cells. The hypergeometric distribution has the peculiarity that, since the distribution of the originating odds ratio is not continuous, you often obtain a different OR as a maximum likelihood estimate.

Answer 3

To answer your second question, biostats isn't my forte but I believe the reason for multiple odds ratio statistics is to account for sampling design and design of experiments.

I've found three references here that will give you a bit of understanding as to why there is a difference between conditional MLE vs unconditional for odds ratio, as well as other types.

1. Point and interval estimation of the common odds ratio in the combination of 2 × 2 tables with fixed marginals

2. The Effect of Bias on Estimators of Relative Risk for Pair-Matched and Stratified Samples

3. A Comparative Study of Conditional Maximum Likelihood Estimation of a Common Odds Ratio

Answer 4

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

1. Define x = matrix(c(3, 6, 5, 6), nrow=2)
2. Run every line of the code extracted from fisher.test in order
3. 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