La trace de Pillai-Bartlett

add\bullet Présentation:

Créé suite aux travaux de K. C. Sreedharan Pillai et de Maurice Stevenson Bartlett, la trace de Pillai-Bartlett est une approche paramétrique permettant de tester si les sous-échantillons \mathbf{X|_{Y = 1}}, \cdots, \mathbf{X|_{Y = K}}, d’une matrice \mathbf{X} de p variables continues qui suivent une loi gaussienne et restreinte aux K groupe d’une variable qualitative Y, ont même barycentre.

La trace de Pillai-Bartlett est également utilisée afin de valider les hypothèses de l’Analyse de Variance multivariée (MANOVA) en travaillant directement sur les matrices résiduelles.

La trace de Pillai-Bartlett a été conçue dans le même objectif que ceux du \lambda de Wilks, de la trace de Hotelling-Lawley, du T ^2 de Hotelling et de la plus forte valeur propre de Roy. De ces cinq outils, il demeure le plus robuste aux violations des hypothèses du modèle.

\bullet Le test:

Hypothèse préliminaire: Matrice \mathbf{X} de p variables continues qui suivent une loi gaussienne et une variable qualitative à K modalités.

La formule de la trace de Pillai est:

V = tr(\mathbf{SCE} \cdot ((\mathbf{SCE} + \mathbf{SCR}) ^{-1}))

\mathbf{SCE} est la matrice de dispersion inter-groupes et \mathbf{SCR} celle de dispersion intra-groupes, de formules respectives:

\mathbf{SCE} = \sum_{k = 1} ^K n_k (\mathbf{\overline{X|_{Y = k}}} - \mathbf{\overline{X}})' \cdot (\mathbf{\overline{X|_{Y = k}}} - \mathbf{\overline{X}})

\mathbf{SCR} = \sum_{k = 1} ^K \mathbf{\hat{X|_{Y = k}}}' \cdot \mathbf{\hat{X|_{Y = k}}}

, où \mathbf{\hat{X|_{Y = k}}} est la matrice centrée non réduite de \mathbf{X|_{y = k}}.

Notons que par la formule de décomposition de la variance issue du théorème de Huygens:

\mathbf{SCT} = \mathbf{\hat{X}} ' \cdot \mathbf{\hat{X}} = \mathbf{SCE} + \mathbf{SCR}

, permettant de relier \mathbf{SCE} et \mathbf{SCR}.

La statistique de test issue de la trace de Pillai-Bartlett est alors:

F = \frac{2 \cdot u + s + 1}{2 \cdot t + s + 1} \cdot \frac{V}{s - V}

, où u, s, t valent:

u = \frac{n - K - p - 1}{2}

s = min(p, K - 1)

t = \frac{|p - K + 1| - 1}{2}

Cette statistique de test suit une loi de Fisher-Snedecor de paramètres (s \cdot (2 \cdot t + s + 1),s \cdot (2 \cdot u + s + 1)) et l’hypothèse H_0 est: « Egalité des barycentres / \mathbf{\overline{X}}|_{Y = 1} = \cdots = \mathbf{\overline{X}}|_{Y = K}« .

Ci-dessous la table de Fisher-Snedecor:

addTendance pour rejeter  H_0:

Plus la statistique F est forte et plus nous avons de chance de rejeter H_0. Ce qui revient étudier le rapport,

F = \frac{2 \cdot u + s + 1}{2 \cdot t + s + 1} \cdot \frac{V}{s - V} \rightarrow \infty

\Rightarrow \frac{V}{s - V} \rightarrow \infty

\Rightarrow \frac{tr(\mathbf{SCE} \cdot \mathbf{SCT} ^{-1})}{s - tr(\mathbf{SCE} \cdot \mathbf{SCT} ^{-1})} \rightarrow \infty

\Rightarrow \frac{s - tr(\mathbf{SCE} \cdot \mathbf{SCT} ^{-1})}{tr(\mathbf{SCE} ^{-1} \cdot \mathbf{SCT} ^{-1})} \rightarrow 0

\Rightarrow \frac{s}{tr(\mathbf{SCE} \cdot \mathbf{SCT} ^{-1})} - 1 \rightarrow 0

\Rightarrow tr(\mathbf{SCE} \cdot \mathbf{SCT} ^{-1}) \rightarrow s = min(p, K - 1)

Ce résultat est logique. Plus la matrice de dipersion inter-groupes (\mathbf{SCE}) s’éloigne de la matrice nulle et plus la matrice dispersion intra-groupes (\mathbf{SCR}) s’en approche, plus le rapport:

\mathbf{SCE} \cdot (\mathbf{SCE} + \mathbf{SCR}) ^{-1} \rightarrow Id

, où Id est la matrice identité de taille p \times p et donc de trace égale à p (puisque somme des éléments de la diagonale de taille p).

En d’autres termes, une dipersion maximale entre les groupes et une dispersion minimale au sein des différentes groupes impliquent que les barycentres des groupes sont éloignés. Cette idée est la base des tests liés à l’analyse multivariée comme l’analyse canonique des corrélations ou encore d’outils tels que le clustering hiérarchique.

Le dernier cas que nous n’avons pas encore vu est celui où s = K - 1, un cas un peu tordu mais qui peut se présenter et pour lesquels les barycentres sont également distincts. Celui où la dispersion totale est en réalité proportionnel à la dispersion inter-groupes. Un tel cas implique nécessairement dans son paroxisme (puisque nous raisonnons en terme de convergence) que:

\mathbf{SCT} = \frac{1}{K - 1} \mathbf{SCE}

, nous obtenons alors:

tr(\mathbf{SCE} \cdot \mathbf{SCT} ^{-1}) = tr(\frac{\mathbf{SCE}}{\frac{1}{K - 1} \mathbf{SCE}}) = K - 1

\bullet Tendance lorsque n \longrightarrow \infty:

Nous cherchons désormais à étudier l’influence des grands échantillons sur la robustesse du test associé à la trace de Pillai-Bartlett. 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 barycentres sont différents d’un groupe à l’autre.

addNous 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 un cas où les barycentres ne devraient pas être statistiquement différents. Le tableau ci-dessous présente ces résultats.

addNous 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é à la trace de Pillai-Bartlett est bien influencé par les grands effectifs.

\bullet Annexe théorique:

Nous proposons ici une justification de la loi que suit la statistique de test de la trace de Pillai-Bartlett.

Rappelons dans un premier temps la définition suivante:

Une matrice \mathbf{M} de taille p \times p a une distribution de Wishart W_p (n, \mathbf{\Sigma}) si \mathbf{M} peut s’écrire:

\mathbf{M} = \mathbf{X'} \cdot \mathbf{X}

, où \mathbf{X} matrice de taille n \times p aléatoire et composées de p vecteurs aléatoires gaussiens de taille n et de même loi N_p(0, \mathbf{\Sigma}) indépendants.

Cette définition permet ainsi de justifier le fait que les matrices de dispersion inter-groupes et intra-groupes suivent des lois de Wishart de paramètres respectifs: W_p (K - 1, \mathbf{\Sigma}) et W_p (n - K - p, \mathbf{\Sigma}), et de mettre l’accent sur l’hypothèse de multinormalité recquise pour l’application du test de Pillai-Bartlett.

Remarquons que la loi de Wishart et sa version inversée peuvent être assimilées toutes les deux à des lois du \chi ^2. Et par définition le rapport de deux lois du \chi ^2 donne une loi de Fisher, nous pouvons donc en déduire que le rapport de deux lois de Wishart de paramètres respectifs (K - 1, \mathbf{\Sigma}) et (n - K - p, \mathbf{\Sigma}) est une loi de Fisher de paramètres F(K - 1, n - K - p).

La loi de Wishart de paramètres (n, p,\mathbf{\Sigma}) a pour densité:

F(\mathbf{M}) = \frac{| \mathbf{M} | ^{n - p - 1} e ^{-\frac{1}{2} \cdot Tr(\mathbf{\Sigma ^{-1}} \mathbf{M})}}{2 ^{\frac{n \cdot p}{2}} | \mathbf{\Sigma} | ^{\frac{n}{2}} \Gamma_p (\frac{n}{2}} = | \mathbf{\Sigma} | ^{\frac{p + 1}{2}} \cdot \frac{| \mathbf{\Sigma ^{-1}} \mathbf{M} | ^{\frac{n - p - 1}{2}} e ^{\- \frac{1}{2} Tr (\mathbf{\Sigma ^{-1}} \mathbf{M})}}{2 ^{\frac{np}{2}} \Gamma_p (\frac{n}{2})}

Où, \Gamma_p(\frac{n}{2}) est la généralisation de la fonction \Gamma(.) au cas multidimensionnelle et de formule:

\Gamma_p (\frac{n}{2}) = \pi ^{\frac{p \cdot (p - 1)}{4}} \prod_{j = 1} ^p \Gamma(\frac{n}{2} + \frac{1 - j}{2})

En rappelant la formule: e ^{Tr(\mathbf{X})} = det(e ^{\mathbf{X}}), nous obtenons en enlevant le passage à la valeur absolue:

\mathbf{\Sigma} ^{\frac{p + 1}{2}} \cdot \frac{(\mathbf{\Sigma ^{-1}} \mathbf{M}) ^{\frac{n - p - 1}{2}} e ^{\- \frac{1}{2} \mathbf{\Sigma ^{-1}} \mathbf{M}}}{2 ^{\frac{np}{2}} \Gamma_p (\frac{n}{2})}

Grossièrement, en effectuant les changements de variables: x = \mathbf{\Sigma ^{-1}} \mathbf{M} et k - 2 = n - p - 1, nous obtenons:

f(X) = \frac{1}{C} \cdot \mathbf{\Sigma} ^{\frac{p + 1}{2}} \cdot \frac{x ^{\frac{k}{2} - 1} e ^{- \frac{x}{2}}}{2 ^{\frac{k}{2}} \Gamma_p (\frac{k}{2})}

, où C est une constante dépendante de 2 ^{(.)} et \prod \Gamma (.). La densité d’une loi de Wishart se rapproche alors de celle d’un \chi ^2

La dernière étape étant de montrer que la loi inverse de Wishart suit également une \chi ^2. Sa densité étant:

\frac{| \mathbf{\Sigma} | ^{\frac{n}{2}}}{2 ^{\frac{n \cdot p}{2}} \Gamma_p (\frac{n}{2})} | \mathbf{M} | ^{- \frac{n + p + 1}{2}} e ^ {- \frac{1}{2} Tr (\mathbf{\Sigma} \mathbf{X ^{-1}})}

Ce qui est automatique en procédant au changement de variable \mathbf{C} = \mathbf{X ^{-1}}, ce qui permet également d’établir que la densité d’une loi inverse de Wishart se rapproche de celle d’un \chi ^2.

\bullet Exemple:

Soit l’échantillon suivant,

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

add

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

Pour le calcul de la trace de Pillai-Bartlett il nous faut calculer la dispersion intra-groupe \mathbf{SCR} et la dispersion inter-groupe \mathbf{SCE}.

\mathbf{SCR} =

add

\mathbf{SCE} =

add

Maintenant que nous avons calculé les deux principaux termes nécessaires, nous pouvons calculer:

\mathbf{SCE} \cdot (\mathbf{SCE} + \mathbf{SCR}) ^{-1} =

addLa trace de Pillai-Bartlett s’obtient en sommant les éléments de la diagonale, nous avons alors:

V = 0.88661435 + 0.05998179 = 0.9465961

Avant de calculer la statistique de test associée à la trace de Pillai-Bartlett, calculons les degrés de liberté. Ainsi,

s = min(2,3-1) = 2

t = \frac{| 2 - 3 + 1 | - 1}{2} = - \frac{1}{2}

u = \frac{20 - 3 - 2 - 1}{2} = 7

ddl_1 = 2 \times (2 \times (- \frac{1}{2}) + 2 + 1) = 4

ddl_2 = 2 \times (2 \times 7 + 2 + 1) = 34

Nous avons donc la statistique de test qui vaut:

F = \frac{2 \times 7 + 2 + 1}{2 \times (- \frac{1}{2}) + 2 + 1} \times \frac{0.9465961}{2 - 0.9465961} = 7.63816

En appliquant la loi de Fisher de paramètres (4, 34) nous obtenons une p-valeur de 0.0001671334 < 5 \%. Nous pouvons donc rejeter H_0 et en déduire que les barycentres sont bien différents au sens statistique du terme.

\bullet Application informatique:

Procédure SAS: http://support.sas.com/documentation/cdl/en/statug/63033/HTML/default/viewer.htm#statug_glm_sect035.htm

Package et fonction R:

Trace.Pillai = function (BDD) {

    # La fonction suivante calcul la trace de Pillai-Bartlett. Il faut mettre en colonne 1 la variable réponse Y et ensuite les différentes variables explicatives.

obs.class = summary(as.factor(BDD[,1]))
nb.class = length(obs.class)
nb.var = dim(BDD)[2]

 # Calcul du SCT
X = BDD[,2] – mean(BDD[,2])
for (p in 3:nb.var) {X = cbind(X,BDD[,p] – mean(BDD[,p]))}
SCT = t(X) %*% X

# Calcul du SCR
SCR = list(nb.class)
for (k in 1:nb.class) {
K.encours = BDD[which(BDD[,1] == names(obs.class)[k]),2] – mean(BDD[which(BDD[,1] == names(obs.class)[k]),2])
for (p in 3:nb.var) {K.encours = cbind(K.encours,BDD[which(BDD[,1] == names(obs.class)[k]),p] – mean(BDD[which(BDD[,1] == names(obs.class)[k]),p]))}
SCR[[k]] = t(K.encours)%*%K.encours
}
SCR.sum = SCR[[1]]
for (k in 2:nb.class) {SCR.sum = SCR.sum + SCR[[k]]}
SCR = SCR.sum

# Calcul du SCE
SCE = list(nb.class)
for (k in 1:nb.class) {SCE[[k]] = (colMeans(BDD[which(BDD[,1] == names(obs.class)[k]),2:nb.var]) – colMeans(BDD[,2:nb.var])) %*% t(colMeans(BDD[which(BDD[,1] == names(obs.class)[k]),2:nb.var]) – colMeans(BDD[,2:nb.var])) * obs.class[k]}
SCE.sum = SCE[[1]]
for (k in 2:nb.class) {SCE.sum = SCE.sum + SCE[[k]]}
SCE = SCE.sum

 # Calcul de la trace de Pillai-Bartlett
S = SCE %*% solve(SCT)
V = sum(diag(S))

# Calcul des paramètres de la loi de Fisher-Snedecor
s = min(nb.var – 1,nb.class – 1)
t = (abs(nb.var – 1 – nb.class + 1) – 1)/2
u = (sum(obs.class) – nb.class – (nb.var – 1) – 1)/2
ddl1 = s * (2 * t + s + 1)
ddl2 = s * (2 * u + s + 1)

# Calcul de la statistique de test F de la p-valeur associée selon la loi de Fisher-snedecor de paramètres (ddl1, ddl2)
F = ((2 * u + s + 1)/(2 * t + s + 1)) * (V/(s – V))
p = (1 – pf(F,ddl1,ddl2))

# Renvoi résultats d’intérêt
return(paste(« Trace de Pillai-Bartlett = »,V, »/ Statistique de test = »,F, »/ p-valeur = »,p))

}

\bullet Bibliographie:

– Composantes de la variabilite. Plans d’experiences de Jean-Marc Azaïs

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

– Chapitre 4 – MANOVA de Catherine d’Aubigny

– Comparaison de populations, tests paramétriques de Ricco Rakotomalala