Le test M de Box

add3

\bullet Présentation:

Publié en 1949 suite aux travaux de George Edward Pelham Box, le test M de Box est une approche paramétrique permettant de tester si les matrices de covariance \mathbf{\Sigma}_{\mathbf{X}|_{Y = 1}}, \cdots, \mathbf{\Sigma}_{\mathbf{X}|_{Y = K}}, associées à \mathbf{X} = (X ^1, \cdots, X ^P) restreint aux K \geq 2 groupes d’une variable qualitative Y, sont égales.

Le test M de Box requiert que les différents sous-échantillons \mathbf{X}|_{Y = 1}, \cdots, \mathbf{X}|_{Y = K} suivent un loi normale.

A noter que malgré son utilité majeure pour l’utilisation de nombreux autres tests nécessitant que l’hypothèse d’homoscédasticité (égalité des matrices de covariance) soit réalisée, il reste un test dont la fiabilité est sujette à caution.

Enfin, le test M de Box est la généralisation du test de Bartlett pour le cas multidimensionnel.

\bullet Le test:

Hypothèse préliminaire: Matrice de P variables continues et une variable qualitative à K \geq 2 modalités, multinormalité.

La statistique du test M de Box est:

M = (n - K) log | \mathbf{W} | - \sum_{k = 1} ^K (n_k - 1) log | \mathbf{\Sigma}_{\mathbf{X}|_{y = k}} |

Où,

\mathbf{W} = \frac{1}{n - K} \sum_{k = 1} ^K (n_k - 1) \mathbf{\Sigma}_{\mathbf{X}|_{Y = k}}

Afin de pouvoir reporter la statistique de test M de Box à la table de la loi du \chi ^2, il convient de lui appliquer la transformation suivante,

\chi ^2 = M (1 - \delta)

Où,

\delta = \frac{2 p ^2 + 3 p + 1}{6 (p + 1) (K - 1)} [(\sum_{k = 1} ^K \frac{1}{n_k - 1}) - \frac{1}{n - K}]

La statistique du test suit alors une loi du \chi ^2 à \frac{p (p + 1) (K - 1)}{2} degrés de libertés et l’hypothèse H_0 est: « Egalité des matrices de covariance / \mathbf{\Sigma}_{\mathbf{X}|_{Y = 1}} = \cdots = \mathbf{\Sigma}_{\mathbf{X}|_{Y = K}}« .

Ci-dessous le tableau de la loi du \chi ^2.

add

Tendance pour le rejet de H_0:

Plus la statistique de test \chi ^2 est grande et plus nous avons de chance de rejeter H_0, ce qui revient à dire que,

\chi ^2 \rightarrow \infty \Rightarrow M (1 - \delta) \rightarrow \infty \Rightarrow M \rightarrow \infty

Le test M de Box se résume à comparer le logarithme de la somme des dispersions intra-groupes à la somme des logarithmes des dispersions intra-groupes. Par simplicité, intéressons-nous au cas où M \rightarrow 0.

Supposons que \mathbf{\Sigma}_{\mathbf{X}|_{Y = 1}} = \cdots = \mathbf{\Sigma}_{\mathbf{X}|_{Y = K}}, si nous reprenons la formule nous avons:

(n - K) log | \mathbf{W} | - \sum_{k = 1} ^K (n_k - 1) log | \mathbf{\Sigma}_{\mathbf{X}|_{Y = k}} | = (n - K) log | \frac{1}{n - K} \sum_{k = 1} ^K (n_k - 1) \mathbf{\Sigma}_{\mathbf{X}|_{Y = k}} | - \sum_{k = 1} ^K (n_k - 1) log | \mathbf{\Sigma}_{\mathbf{X}|_{Y = K}} |

= (n - K) log(\frac{1}{n - K} | \sum_{k = 1} ^K (n_k - 1) \mathbf{\Sigma}_{\mathbf{X}|_{Y = k}} | ) - \sum_{k = 1} ^K (n_k - 1) log | \mathbf{\Sigma}_{\mathbf{X}|_{Y = k}} |

= (n - K) [log | \sum_{k = 1} ^K (n_k - 1) \mathbf{\Sigma}_{\mathbf{X}|_{Y = k}} | - log(n - K)] - \sum_{k = 1} ^K (n_k - 1) log | \mathbf{\Sigma}_{\mathbf{X}|_{Y = k}} |

= (n - K) [log | \mathbf{\Sigma}_{\mathbf{X}|_{Y = 1}} | \sum_{k = 1} ^K (n_k - 1) - log(n - K)] - log | \mathbf{\Sigma}_{\mathbf{X}|_{Y = 1}} | \sum_{k = 1} ^K (n_k - 1)

= (n - K) log | \mathbf{\Sigma}_{\mathbf{X}|_{Y = 1}} | + (n - K) log (\sum_{k = 1} ^K (n_k - 1))  - (n - K) log(n - K) - log | \mathbf{\Sigma}_{\mathbf{X}|_{Y = 1}} | \sum_{k = 1} ^K

Or \sum_{k = 1} ^K (n_k - 1) = \sum_{k = 1} ^K n_k - K = n - K, par conséquent,

= (n - K) log | \mathbf{\Sigma}_{\mathbf{X}|_{Y = 1}} | + (n - K) log(n - K) - (n - K) log(n - K) - log | \mathbf{\Sigma}_{\mathbf{X}|_{Y = 1}} | (n - K)

= 0

\bullet Tendance lorsque n \longrightarrow \infty:

Nous cherchons désormais à étudier l’influence des grands échantillons sur la robustesse du test M de Box. Nous générons des matrices de 5 vecteurs gaussiens pour 3 groupes. Nous ne faisons pas varier le nombre de groupes ni le nombre de variables aléatoires mais seulement la taille d’échantillon.

Le tableau ci-dessous présente l’évolution des p-valeurs associées aux statistiques de test calculées sur plusieurs simulations dans le cas où les matrices de covariance sont différentes d’un groupe à l’autre.

add.png

Nous constatons que l’hypothèse est bien rejetée quelque soit la taille de n, ce qui correspond au résultat attendu d’après nos hypothèses.

Procédons à la même expérience mais cette fois-ci dans le cas où les matrices de covariance ne devraient pas être statistiquement différentes. Le tableau ci-dessous présente ces résultats.

add.png

Nous constatons que jusqu’à n = 100 le résultat du test est en adéquation avec nos hypothèses mais qu’entre 100 et 1000 observations ce dernier devient significatif, rejetant  H_0 à tort.

Nous en concluons que le test associé au M de Box est bien influencé par les grands effectifs.

\bullet Annexe théorique:

Cette partie présente une esquisse de la démonstration que le M de Box suit une loi du \chi ^2.

Soit s_{j_1,j_2,k} l’estimateur non biaisé de la variance ou de la covariance de A ^{j_1,j_2,k} entre la j_1-ème et la j_2-ème variable pour le k-ème groupe, basé sur la somme des carrés et produits ayant v_k degrés de liberté. Et posons,

s_{j_1,j_2} = \frac{\sum_{k = 1} ^K v_k s_{j_1,j_2,k}}{N}

, avec N = \sum_{k = 1} ^K v_k. Nous avons donc,

M = N log_e | S_{j_1, j_2} | - \sum_{k = 1} ^K (v_k log_e | s_{j_1,j_2} |) = - N log(L_1 ')

, où L_1 ' = \prod_{k = 1} ^K \lbrace \frac{| S_{j_1,j_2,k} |}{| S_{j_1, j_2} |} \rbrace ^{\frac{v_k}{N}}.

Si les degrés de liberté sont égaux, alors M, L_1 ' peuvent être reliés par le l_1 de Bishop de forme suivante,

M = -2 N log_e l_1, L_1 ' = l_1 ^2

Nous cherchons à obtenir les moments de L_1 ' sous l’hypothèse H_0. Soit c_{j_1,j_2,k} la somme des carrés et des produits basés sur les degrés de liberté v_k correspondant à s_{j_1,j_2,k} tel que,

c_{j_1,j_2,k} = v_k s_{j_1,j_2,k}, c_{j_1,j_2} = N s_{j_1,j_2}

La densité de distribution conjointe de c_{j_1,j_2,k} pour le k-ième groupe est donné par la distribution de Wishart,

p(c_{1,1,k}, \cdots, c_{P,P,k}) = Z(v_k) | c_{j_1,j_2,k} | ^{\frac{1}{2} (v_k - P - 1)} e ^{-\frac{1}{2} \sum_{j_1,j_2} E_{j_1,j_2,k} c_{j_1,j_2,k}}

Où,

Z (v_k) ^{-1} = 2 ^{\frac{1}{2} (v_k P)} \pi ^{\frac{1}{2} P (P-1)} \prod_{j_2 = 1} ^{P - 1} \Gamma(\frac{v_k - j_2}{2}) | A_k | ^{-\frac{1}{2} v_k}

Sous l’hypothèse H_0, A_k (abréviation de A_{j_1,j_2,k}) est identique \forall k \in [1,K]. Le g-ième moment de | c_{j_1,j_2} | est alors,

\int \prod_{k = 1} ^K \lbrace Z (v_k) | c_{j_1,j_2} | ^{\frac{1}{2} (v_k - P - 1)} | c_{j_1, j_2} | ^{\frac{g}{k}} e ^{- \frac{1}{2} \sum_{j_1,j_2} A_{j_1,j_2} c_{j_1,j_2,k}} \partial c_{1,1,1} \partial c_{1,2,1} \cdots \partial c_{P,P,K}

En posant v_k (1 + 2h) degrés de liberté et g = -N h, si nous intégrons sur l’espace pour lesquels c_{j_1,j_2}, c_{j_1,j_2,k} sont définies positifs, alors,

E[ \prod_{k = 1} ^K \lbrace \frac{| c_{j_1,j_2,k} |}{| c_{j_1,j_2} |} \rbrace ] ^h \prod_{k = 1} ^K [\frac{Z(v_k (1 + 2 h))}{Z (v_k)}] = \frac{Z(N (1 + 2 h))}{Z (N)}

D’où,

E[L_1'] ^{N h} = E[e ^{- M}] ^h

= \frac{Z(N (1 + 2 h))}{Z(N)} \prod_{k = 1} ^K [\frac{Z (v_k)}{Z(N)}  (\frac{N}{v_k}) ^{p v_k h}]

= \prod_{k = 1} ^K (\frac{N}{v_k}) ^{h v_k P} \prod_{j_2 = 0} ^{P - 1} [\frac{\Gamma (\frac{1}{2} (N - j_2))}{\Gamma (\frac{1}{2} (N (1 + 2 h) - j_2))} \prod_{k = 1} ^K \frac{\Gamma (\frac{1}{2} v_k (1 + 2 h) - j_2)}{\Gamma (\frac{1}{2} (v_k - j_2))}]

, qui est une équation basée sur l’identité analytique du réel h. En posant \rho \leq 1, N = v K et h = - i t \rho complexe, nous pouvons déterminer la fonction caractéristique de \rho M. Auparavant, définissons les quantités suivantes,

\mu = \rho v_k, \mu = \rho v, \rho = \mu + \beta, v_k = \mu_k + \beta_k

Alors,

\Phi(t) = \prod_{k = 1} ^K (\frac{K \mu}{\mu_k}) ^{- i t P \mu_k} \prod_{j_2 = 0} ^(P - 1) [\frac{\Gamma(\frac{1}{2} K \mu + K \beta - j_2)}{\Gamma (\frac{1}{2}(K \mu (1 - 2 i t) + K \beta - j_2))} \prod_{k = 1} ^K \frac{\Gamma (\frac{1}{2} (\mu_k (1 - 2 i t)+ \beta_k - j_2))}{\Gamma(\frac{1}{2}(\mu_k + \beta_k - j_2))}]

En passant au logarithme, nous obtenons la fonction génératrice suivante,

\Psi (t) = g(t) - g(0)

Où,

g(t) = - \sum_{k = 1} ^K it \mu_k P log(\frac{K \mu}{\mu_k}) + \sum_{j_2 = 1} ^{P - 1} [\sum_{k = 1} ^K \lbrace log (\Gamma (\frac{1}{2} (\mu_k (1 - 2 i t) + \beta_k - j_2)))] \rbrace - log (\Gamma (\frac{1}{2} (K \mu (1 - 2 i t) + K \beta - j_2))]

, et g(0) constante indépendante de t.

Grâce à la généralisation du théorème de Stirling, nous avons,

log (\Gamma(x + h)) = log (\sqrt{2 \pi + (x + h - \frac{1}{2}) log(x) - x - \sum_{i' = 1} ^n (-1) ^{i'} \frac{B_{i' + 1} (h)}{i' (i' + 1) x ^{i'}} + R_{n + 1} (x)}

, avec B_i polynôme de degré i et d’ordre un défini par,

\frac{\gamma e ^{h i'}}{e ^{i'} - 1} = \sum_{i' = 1} ^{\infty} B_{i'} (h)

En développant chaque \Gamma-fonction grâce à cet objet, nous avons,

\Psi(t) = Q - q(0) - \frac{1}{4} (K - 1) P (P + 1) log(1 - 2 i t) + \sum_{i = 1} ^n \frac{\alpha_{i'}}{\mu ^{i'}} (1 - 2 i t) ^{-i'} + R_{n + 1} (\mu, t)

Avec,

Q = \frac{P (K - 1)}{2} log(2 \pi) + \frac{P}{2} [\sum_{k = 1} ^K (v_k - \frac{P + 1}{2}) log(\frac{\mu_k}{2}) - (K v - \frac{P + 1}{2}) log(\frac{K \mu}{2})]

\alpha_{i'} = \frac{(-1) ^{i' + 1} (2 \mu) ^{i'}}{i' (i' + 1)} \sum_{j_2 = 0} ^{P - 1} [\sum_{k = 1} ^K \frac{B_{i' + 1} (\frac{B_k - j_2}{2})}{\mu_k ^{i'}} - \frac{B_{i' + 1} (\frac{K \beta - j_2}{2})}{\mu ^{i'} K ^{i'}}]

Nous avons alors,

\Psi(t) = K(1 - 2 i t) ^{-\frac{1}{2} f} \sum_{v = 0} ^{P - 1} \frac{a_v}{\mu ^v} (1 - 2 i t) ^{-v} + R_{n + 1} ' (\mu, t)

, avec K = e ^{Q - g(0)}, f = \frac{1}{2} (K - 1) P (P + 1) et a_v coefficient de \mu ^{-v} dans le développement de e ^{\sum_{i' = 1} ^n \frac{\alpha_{i'}}{\mu ^{i'}}}

De ce fait, la fonction de densité de probabilité g de \rho M est alors,

g(\rho M) = \frac{1}{2 \pi} \int_{- \infty} ^{+ \infty} e ^{-i t \rho M} \Phi (t) = K \sum_{v = 0} ^{\infty} \frac{a_v}{\mu ^v} g(\chi_{f + 2 v} ^2) + R_{n + 1} '' (\mu, t)

Plus particulièrement, la probabilité est alors donnée par,

P(M \geq M_0) = K \sum_{v = 0} ^{\infty} \frac{a_v}{\mu ^v} G_{f + 2 v} + R_{n + 1} ''' (\mu,t)

Avec,

G_{f + 2v} = \int_{\rho M_0} g(\chi_{f + 2v} ^2) \partial \chi ^2

R_{n + 1} ''' \rightarrow 0, n \rightarrow \infty

Nous en déduisons donc que le M de Box suit une loi du \chi ^2.

\bullet Exemple:

Soit l’échantillon suivant,

add

La projection des couples X ^1, X ^2 nous donne:

add

A vue d’oeil nous serions tenter de dire que les 3 dispersions des classes A, B, C semblent équivalentes. Vérifions-le au sens statistique du terme.

Dans un premier temps, calculons les matrices de covariance par groupe et totale,

\mathbf{\Sigma}_{\mathbf{X}|_{Y = 1}} = \begin{pmatrix} 4.642112 & -4.505277 \\ -4.505277 & 13.198839 \\ \end{pmatrix}

\mathbf{\Sigma}_{\mathbf{X}|_{Y = 2}} = \begin{pmatrix} 3.476513 & 1.00403 \\ 1.004030 & 11.45007 \\ \end{pmatrix}

\mathbf{\Sigma}_{\mathbf{X}|_{Y = 3}} = \begin{pmatrix} 4.758668 & 3.655501 \\ 3.655501 & 9.107011 \\ \end{pmatrix}

\Rightarrow \mathbf{\Sigma} = \frac{1}{20 - 3} \times [(7 - 1) \times \mathbf{\Sigma}_{\mathbf{X}|_{Y = 1}} + (6 - 1) \times \mathbf{\Sigma}_{\mathbf{X}|_{Y = 2}} + (7 - 1) \times \mathbf{\Sigma}_{\mathbf{X}|_{Y = 3}}]

= \frac{1}{17} \times [6 \times \mathbf{\Sigma}_{\mathbf{X}|_{Y = 1}} + 5 \times \mathbf{\Sigma}_{\mathbf{X}|_{Y = 2}} + 6 \times \mathbf{\Sigma}_{\mathbf{X}|_{Y = 3}}]

= 0.05882353 \times \begin{pmatrix} 73.78724915 & -0.07850224 \\ -0.07850224 & 191.08543334 \\ \end{pmatrix}

= \begin{pmatrix} 4.340426421 & -0.0004617779 \\ -0.004617779 & 11.240319608 \\ \end{pmatrix}

Nous pouvons désormais calculer le M de Box,

M = (20 - 3) \times log (det( \mathbf{\Sigma})) - [(7 - 1) \times log(det(\mathbf{\Sigma}_{\mathbf{X}|_{Y = 1}})) + (6 - 1) \times log(det(\mathbf{\Sigma}_{\mathbf{X}|_{Y = 2}})) + (7 - 1) \times log(det(\mathbf{\Sigma}_{\mathbf{X}|_{Y = 3}}))]

= 17 \times log (48.78776) - [6 \times log(40.97297) + 5 \times log(38.79823) + 6 \times log(29.97455) ]

= 66.08715 - (22.27748 + 18.29187 + 20.40209)

= 66.08715 - 60.97144

= 5.11571

Maintenant vient le calcul du coefficient de transformation,

\delta = \frac{2 \times 2 ^2 + 3 \times 2 + 1}{6 \times (2 + 1) \times (3 - 1)} \times [(\sum_{k = 1} ^3 \frac{1}{n_k - 1}) - \frac{1}{20 - 3}]

= \frac{13}{36} \times (0.5333333 - 0.05882353)

= 0.3611111 \times 0.4745098

= 0.1713508

En appliquant la transformation permettant de reporter la statistique de test M à la table de la loi du \chi ^2 nous obtenons,

\chi ^2 = 5.11571 \times (1 - 0.1713508) = 5.11571 \times 0.8286492 = 4.239129

Si nous reportons cette statistique de test à la table de la loi du \chi ^2 pour \frac{2 \times (2 + 1) \times (3 - 1)}{2} = \frac{12}{2} = 6 degrés de liberté, nous obtenons une p-valeur de 0.6444 > 5 \%. Nous concluons au non-rejet de H_0 et donc à l’égalité des matrices de covariance.

\bullet Application informatique:

Procédure SAS: https://support.sas.com/documentation/cdl/en/statug/63347/HTML/default/viewer.htm#statug_discrim_sect004.htm

Package et fonction R: http://finzi.psych.upenn.edu/library/biotools/html/boxM.html

\bullet Bibliographie:

– A general distribution theory for a class of likelihood criteria par George Edward Pelham Box

– Probabilité, analyse de données et Statistique de Gilbert Saporta