Une implémentation en javascript de la version non centrée du test exact de Fisher

Une implémentation en javascript de la version non centrée du test exact de Fisher

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 étonné 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
}\]

\[\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!

By | 2017-01-16T09:25:09+00:00 9 janvier 2017|Categories: Javascript, Statistiques, Test|0 Commentaires

About the Author:

Informaticien de formation, j'occupe mon temps à faire parler les données biologiques… tous les moyens sont bons!

Laisser un commentaire