Mon article précédent expliquait pourquoi la version non-centrée du test exact de Fisher est plus appropriée dans la plupart des cas rencontrés en bio-informatique. Je poursuis en présentant maintenant une implémentation de ce test en Javascript qui pourrait facilement être intégrée à une interface web.

Même si le Javascript est un langage très mal adapté à l’implémentation de méthodes statistiques, j’espère que cet article présentera tous les détails nécessaires pour simplifier l’implémentation de ce test dans d’autres langages, selon les besoins. À tout le moins, cet exemple devrait démystifier le fonctionnement interne de ce test très utile. À noter que je me limiterai seulement aux tables de contingence $2\times 2$. Avant de plonger, je dois admettre avoir été très agréablement surpris de la performance de ce code!

Les premières fonctions ci-dessous calculent différentes expressions requises par la distribution hypergéométrique non-centrée. Plus précisément, elles utilisent des log-probabilités pour éviter la manipulation de petits nombres ce qui risquerait de produire un underflow.

  function lngamm (z) {
     // Reference: "Lanczos, C. 'A precision approximation
     // of the gamma function', J. SIAM Numer. Anal., B, 1, 86-96, 1964."
     // Translation of  Alan Miller's FORTRAN-implementation
     // See http://lib.stat.cmu.edu/apstat/245
     
     var x = 0;
     x += 0.1659470187408462e-06 / (z + 7);
     x += 0.9934937113930748e-05 / (z + 6);
     x -= 0.1385710331296526     / (z + 5);
     x += 12.50734324009056      / (z + 4);
     x -= 176.6150291498386      / (z + 3);
     x += 771.3234287757674      / (z + 2);
     x -= 1259.139216722289      / (z + 1);
     x += 676.5203681218835      / (z);
     x += 0.9999999999995183;
     return (Math.log (x) - 5.58106146679532777 - z + (z - 0.5) * Math.log (z + 6.5));
  }
  function lnfact (n) {
     if (n <= 1) return (0);
     return (lngamm (n + 1));
  }
   
  function lnbico (n, k) {
     return (lnfact (n) - lnfact (k) - lnfact (n - k));
  }

La fonction « lnbico » (ln of the binomial coefficient) retourne le log du nombre de combinaisons possibles pour choisir
$k$ dans $n$:

\begin{eqnarray*}
\binom{n}{k} & = & \frac{n!}{k!(n-k)!}\\
\mbox{lnfact}(n) & = & \ln(n!)\\
\mbox{lnbico}(n,k) & = & \ln\binom{n}{k}\\
& = & \mbox{lnfact}(n) – (\mbox{lnfact}(k) + \mbox{lnfact}(n – k))
\end{eqnarray*}

La distribution de probabilités de la loi hypergéométrique non-centrée est la suivante (où $\omega$ est le seuil du odds-ratio):

\[
\frac{\binom{m_1}{x} \binom{m_2}{n-x} \omega^x}{
\sum_{y=x_{\min}}^{x_{\max}} \binom{m_1}{y} \binom{m_2}{n – y} \omega^y
}\]

Où

\[\begin{eqnarray*}
x_{\min} & = & \max(0, n – m_2)\\
x_{\max} & = & \min(n, m_1)
\end{eqnarray*}\]

Notation utilisée dans le code pour la matrice de contingence et les totaux:

\[\begin{array}{|cc||c}\hline
\bf n_{11} & \bf n_{12} & n \\
\bf n_{21} & \bf n_{22} &  \\\hline\hline
m_1 & m_2 &  \\
\end{array}\]

Pour calculer une p-value pour le test, les comptes doivent être réarrangés et j’ai utilisé la notation suivante:

\[\begin{array}{|cc||c}\hline
x & ? & n \\
? & ? &  \\\hline\hline
m_1 & m_2 &  \\
\end{array}\]

Ici, décider d’un $x$ permet de remplir la table en utilisant les totaux. La p-value du test exact de Fisher est ensuite calculée en additionnant les probabilités pour le scénario où $x$ est plus grand que la valeur observée $n_{11}$. Ceci est gardé dans « den_sum » et est calculé par la dernière boucle « for » de la fonction. Notez que toutes les sommes proviennent de probabilités ajustées ( – max_l en espace log). Cela évite de faire la somme de très petits nombres, ce qui pourrait conduire à un underflow. Comme la p-value finale est le résultat d’un ratio, le facteur d’ajustement s’annule de lui-même dans Math.exp (den_sum – sum_l). Finalement, « w » est utilisé pour passer le seuil requis pour le odds-ratio ($\omega$).

 function exact_nc (n11, n12, n21, n22, w) {
   var x = n11;
   var m1 = n11 + n21;
   var m2 = n12 + n22;
   var n = n11 + n12;
   var x_min = Math.max (0, n - m2);
   var x_max = Math.min (n, m1);
   var l = [];

   for (var y = x_min; y <= x_max; y++) {
     l[y - x_min] = (lnbico (m1, y) + lnbico (m2, n - y) + y * Math.log (w));
   }
   var max_l = Math.max.apply (Math, l);
 
   var sum_l = l.map (function (x) { return Math.exp (x - max_l); }).reduce (function (a, b) {
     return a + b; }, 0);
   sum_l = Math.log (sum_l);
 
   var den_sum = 0;
   for (var y = x; y <= x_max; y++) {
     den_sum += Math.exp (l[y - x_min] - max_l);
   }
   den_sum = Math.log (den_sum);
   return Math.exp (den_sum - sum_l);
 };

Veuillez noter que ce code a été testé, mais pas de façon extensive. J’ai confirmé qu’il retournait des résultats identiques à l’implémentation trouvée dans R (fisher.exact) sur un nombre limité de cas non triviaux. Si vous tombez sur des problèmes ou sur des erreurs, je serai heureux de les investiguer!