Calculer un chi 2 sur un tableau de contingence

Links: notebook, html ., PDF, python, slides ., GitHub

\chi_2 et tableau de contingence, avec numpy, avec scipy ou sans.

from jyquickhelper import add_notebook_menu
add_notebook_menu()
run previous cell, wait for 2 seconds

formule

Le test du \chi_2 (wikipedia) sert à comparer deux distributions. Il peut être appliqué sur un tableau de contingence pour comparer la distributions observée avec la distribution qu’on observerait si les deux facteurs du tableau étaient indépendants. On note M=(m_{ij}) une matrice de dimension I \times J. Le test du \chi_2 se calcule comme suit :

  • M = \sum_{ij} m_{ij}
  • \forall i, \; m_{i \bullet} = \sum_j m_{ij}
  • \forall j, \; m_{\bullet j} = \sum_i m_{ij}
  • \forall i,j \; n_{ij} = \frac{m_{i \bullet} m_{\bullet j}}{N}

Avec ces notations :

T = \sum_{ij} \frac{ (m_{ij} - n_{ij})^2}{n_{ij}}

La variable aléatoire T suit asymptotiquement une loi du \chi_2 à (I-1)(J-1) degrés de liberté (table). Comment le calculer avec numpy ?

tableau au hasard

On prend un petit tableau qu’on choisit au hasard, de préférence non carré pour détecter des erreurs de calculs.

import numpy
M = numpy.array([[4, 5, 2, 1],
                 [6, 3, 1, 7],
                 [10, 14, 6, 9]])
M
array([[ 4,  5,  2,  1],
       [ 6,  3,  1,  7],
       [10, 14,  6,  9]])

calcul avec scipy

Evidemment, il existe une fonction en python qui permet de calculer la statistique T : chi2_contingency.

from scipy.stats import chi2_contingency
chi2, pvalue, degrees, expected = chi2_contingency(M)
chi2, degrees, pvalue
(6.1685985038926212, 6, 0.40457120905808314)

calcul avec numpy

N = M.sum()
ni = numpy.array( [M[i,:].sum() for i in range(M.shape[0])] )
nj = numpy.array( [M[:,j].sum() for j in range(M.shape[1])] )
ni, nj, N
(array([12, 17, 39]), array([20, 22,  9, 17]), 68)

Et comme c’est un usage courant, numpy propose une façon de faire sans écrire une boucle avec la fonction sum :

ni = M.sum(axis=1)
nj = M.sum(axis=0)
ni, nj, N
(array([12, 17, 39]), array([20, 22,  9, 17]), 68)
nij = ni.reshape(M.shape[0], 1) * nj / N
nij
array([[  3.52941176,   3.88235294,   1.58823529,   3.        ],
       [  5.        ,   5.5       ,   2.25      ,   4.25      ],
       [ 11.47058824,  12.61764706,   5.16176471,   9.75      ]])
d = (M - nij) ** 2 / nij
d.sum()
6.1685985038926212