Le test de Shapiro-Wilk

add

Samuel Sanford Shapiro (à gauche) et Martin Bradbury Wilk (à droite)

(attention MAJ V2 en cours!!!)

\bullet Historique:

\begin{tabular}{|l|c|c|} \hline Bloc & 11/05/2013-V1 & 11/04/2019-V2 \\ \hline Historique &  & Cr\'eation \\ \hline Sommaire &  & Cr\'eation \\ \hline Pr\'esentation & Cr\'eation & MAJ \\ \hline Le test & Cr\'eation, appelation: le test / indice-corr\'elation & Calcul p-valeur, Shapiro-Francia, Ryan-Joiner \\ \hline Tendance lorsque... & Cr\'eation & MAJ \\ \hline Annexe th\'eo... & Cr\'eation &  \\ \hline Exemple & Cr\'eation & Shapiro-Francia, Ryan-Joiner \\ \hline Appli... info... & Cr\'eation & ??? \\ \hline Bibliographie & Cr\'eation & ??? \\ \hline \end{tabular}

\bullet Sommaire:

  • Présentation
  • Le test
    • La table de Shapiro-Wilk
    • Calcul de la p-valeur exacte
    • Condition pour le rejet de H_0
    • Le test de Shapiro-Francia
    • Le test de Ryan-Joiner
  • Tendance lorsque n \rightarrow + \infty
  • Annexe théorique
  • Exemple
    • Test de Shapiro-Wilk
    • Test de Shapiro-Francia
    • Test de Ryan-Joiner
  • Application informatique
    • Procédure SAS
    • Fonction et package R
  • Bibliographie

\bullet Présentation:

Publié en 1965 par Samuel Sanford Shapiro et Martin Bradbury Wilk, le test de Shapiro-Wilk est une approche non paramétrique permettant de tester si une variable continue X suit une loi normale.

A l’origine, le test ne pouvait s’appliquer que sur des échantillons de taille comprise entre 3 et 50. En 1972, les travaux de Shapiro et Francia ont permis d’étendre cette contrainte à des échantillons inférieur à 100 observations. C’est ensuite en 1982, grâce à J. P. Royston que le test devient utilisable sur des tailles comprises entre 3 et 2 000. Enfin, en 1997, M. Mahibbur Rahman et Z. Govindarajulu rendent son utilisation jusqu’à des échantillons de taille 5000.

Actuellement, le test de Shapiro-Wilk et ses variantes, qui seront présentées dans cet article, demeurent les tests de détection de la normalité les plus efficaces pour des échantillons incluant moins de 5000 observations.

\bullet Le test:

Hypothèse préliminaire: Variable X continue.

Posons a_i le coefficient de pondération associé à l’individu i, 1 \leq i \leq n, n la taille de notre variable X. En notant E[c] la partie entière de c et en partant du principe que X est ordonné, la statistique du test W de Shapiro-Wilk est alors:

W = \frac{(\sum_{i = 1} ^{E[\frac{n}{2}]} a_i \times (X_{n - i + 1} - X_i))^2}{\sum_{i = 1} ^n (X_i - \overline{X}) ^2}

Elle suit une loi de Shapiro-Wilk et l’hypothèse H_0 est:

X suit une loi normale / F_X = F_{L(\mu,\sigma)}

Avec W_{1 - \alpha} la valeur seuil de la distribution de la statistique de test  W pour une confiance \alpha, l’hypothèse alternative est alors,

H_1: F_X \neq F_{L(\mu,\sigma)}, soit W < W_{1 - \alpha}

La table de Shapiro-Wilk:

Dans un premier temps, on trouvera la table des coefficients a_i en fonction de n,

FIGURE23Puis celle de la loi de Shapiro-Wilk:

add

Calcul de la p-valeur exacte:

Le calcul de la p-valeur exacte associée au test de Shapiro-Wilk se fait au travers d’un algorithme proposé par J. P. Royston et connu sous le nom de AS 181.2 et qui s’articule autour de la résolution du polynôme algébrique d’ordre 3 et 4. En notant W_{obs} la statistique de test de Shapiro-Wilk et en fonction de la taille d’échantillon:

– Si n = 3 alors p = \frac{6}{\pi} (arcsin(W) - arcsin(\sqrt{0.75}) ;

– Si n \in [4,11] alors on applique le changement de variable: [latex]W' = log(1 - W_{obs}) et on calcul dans un premier temps le polynôme suivant:

g = -2.273 + 0.459 n

Si y > g alors p < 0.0001

Sinon on applique le second changement de variable: W'' = -log(g - W') et on calcul les deux polynômes suivants:

\mu_W = 0.544 - 0.39978 n + 0.025054 n ^2 - 0.0006714 n ^3

\sigma_W = exp(1.3822 - 0.77857 n + 0.062767 n ^3 - 0.0020322 n ^3)

Alors, p = F_{L(\mu_W,\sigma_W^2)} (W'')

- Si n \geq 12 alors on applique le changement de variable: W' = log(1 - W_{obs}) et on calcul les deux polynômes suivants:

\mu_W = -1.5861 - 0.31082 log(n) - 0.083751 log(n) ^2 + 0.0038915 log(n) ^3

\sigma_W = exp(-0.4803 - 0.082676 log(n) + 0.0030302 log(n) ^3)

Alors, p = F_{L(\mu_W,\sigma_W^2)} (W')

Condition pour le rejet de H_0: Plus la statistique de test s'éloigne de 1 et plus grande sont les chances de rejeter H_0.

Finalement W \longrightarrow 0 soit:

- quand le numérateur, c'est-à-dire la combinaison linéaire non-biaisée optimale,  \sum_{i = 1} ^{E[\frac{n}{2}]} a_i \times (X_{n - i + 1} - X_i))^2 \longrightarrow 0, ce qui équivaut à dire que la variance pondérée par les coefficients de Shapiro-Wilk ne s'équilibre pas

- quand le dénominateur \sum_{i = 1} ^n (X_i - \overline{X}) ^2 \longrightarrow \infty, c'est-à-dire quand la variance empirique est trop forte

La distribution normale peut alors être vue comme une pente de coefficient 1 sur un diagramme P-P et la statistique de test de Shapiro-Wilk permet de mesurer de combien les données s'éloignent de cette droite.

Le test de Shapiro-Francia:

S. S. Shapiro et R. S. Francia développe en 1972 une version simplifiée du test de Shapiro-Wilk en remplaçant les coefficients a_i par les quantiles normaux. Son principal apport est de permettre d'appliquer l'idée du test de Shapiro-Wilk sur des échantillons de taille plus grande ainsi qu'un gain de puissance dans certaines configuration des données.

La formule de la statistique de test de Shapiro-Francia:

z = \frac{log(1 - rho_p) - \mu}{\sigma}

Avec:

- qq = (m_1, \cdots, m_n), vecteur des moments d'ordre 1 à n dans le cadre d'une loi normale ;

- \rho_p = cor(X,qq) ^2, coefficient de corrélation de Pearson ;

- \mu = -1.2725 + 1.0521 (log(log(n)) - log(n)) ;

- \sigma = 1.0308 - 0.26758 (log(log(n)) - \frac{2}{log(n)}).

La statistique de test de Shapiro-Francia suit une loi normale centrée-réduite et l'hypothèse H_0 est:

X suit une loi normale / F_X = F_{L(\mu,\sigma)}

Avec z_{1 - \alpha} la valeur seuil de la distribution de la statistique de test  z pour une confiance \alpha, l'hypothèse alternative est alors,

H_1: F_X \neq F_{L(\mu,\sigma)}, soit z < z_{1 - \alpha}

Le test de Ryan-Joiner:

A. Ryan, Barbara Ryan et Brian L. Joiner développe en 1976 une version alternative du test de Shapiro-Francia. Le test se base sur le coefficient de corrélation est calculé entre X et la séquence de vecteur de points de probabilité: z_i = \frac{i-0.375}{i + 1 - 0.75}.

On détermine la valeur de rejet en fonction de n:

- Si n \leq 50 alors \rho_{seuil} = 1.0063 - \frac{0.1288}{\sqrt{n}} - \frac{0.6118}{n} + \frac{1.3505}{n^2} ;

- Si n > 50 alors \rho_{seuil} = 0.9995 - \frac{0.0178}{\sqrt{n}} - \frac{1.7726}{n} + \frac{3.5582{n^2}.

On rejettera alors l'hypothèse de normalité si \rho < \rho_{seuil}.

\bullet Tendance lorsque n \longrightarrow \infty:

On s'intéresse désormais à la résistance du test de Shapiro-Wilk au fur et à mesure que la taille d'échantillon croît. Nous fixons une statistique de test W_{test} = 0.998, soit un cas fictif correspondant à un ratio variance pondérée et variance empirique équivalente. Le résultat attendu sera donc que quelque soit la taille de l'échantillon, l'on ne devrait pas rejeter  l'hypothèse de similarité des fonctions de répartition: H_0. Le graphique ci-dessous montre l'évolution de la p-valeur p associée à W_{test} fixée lorsque n croît de 10 à 5000 observations.

add

De manière hâtive, on reste en adéquation avec l'hypothèse de construction de la statistique de test de Shapiro-Wilk D_{test} jusqu'à n = 1 123 (p > 20 \%). Jusqu'à n = 1 585, on se forcera à rejeter H_0 avec un risque assez fort compris entre  5 \% et 20 \%. Enfin, à n = 2 110 la p-valeur passe en dessous des 1 \%.

Cette simulation montre que le test de Shapiro-Wilk est atteint par la malédiction des grands échantillons.

\bullet Annexe théorique: 

Nous présentons la démonstration du test de Shapiro-Wilk:

Soit m^T = (m_1, \cdots, m_n) les espérances des statistiques d'ordre d'un échantillon distribué selon une loi normale centrée-réduite. Et V = (v_{ij}) la matrice n \times n de covariance de ces statistiques d'ordre.

Notons X^T = (X_1, \cdots, X_n) échantillon aléatoire. Trivialement, si X suit une loi normale alors on peut écrire chaque X_i comme une combinaison linéaire non biaisée:

X_i = \mu + \sigma \cdot Y_i \hspace{5mm} \forall i \in \lbrace 1, \cdots, n \rbrace

Le théorème et le lemme qui suivent permettent de déterminer les estimateurs non biaisés de \mu et \sigma.

Théorème des moindres carrés généralisés d'Aitken (1935) et Lloyd (1952): il est possible de définir les meilleurs estimateurs non-biaisés de \mu et \sigma comme ceux minimisant la forme quadratique (X - \mu \cdot 1 - \sigma \cdot m) ^T V^{-1} (X - \mu \cdot 1 - \sigma \cdot m) et de formule:

\hat{\mu} =\frac{m V^{-1} (m \cdot 1^T - 1 \cdot m^T) V ^{-1} X}{1 \cdot V^{-1} \cdot 1 m^T V^{-1}m - (1^T \cdot V^{-1} m)^2}

\hat{\sigma} = \frac{1 \cdot V^{-1} (m \cdot 1^T - 1 ^cdot m^T) V ^{-1} X}{1 \cdot V^{-1} \cdot 1 \cdot m^T V^{-1}m - (1^T \cdot V^{-1} m)^2}

Lemme: Dans le cas d'une configuration symétrique de la distribution, en particulier pour la distribution normale, 1^T \cdot V^{-1} m = 0.

Nous en déduisons que:

\hat{\mu} = \frac{1 \cdot V^{-1} X}{1^T \cdot V^{-1} \cdot 1}

\hat{\sigma} = \frac{m^T V^{-1} X}{m^T V^{-1} m}

Nous avons alors la statistique de test W de Shapiro-Wilk qui quantifie la  différence entre l'estimation de l'écart-type \sigma ^2 par la méthode de construction de la pente de la régression linéaire (hypothèse de normalité) et son estimation classique (aucune hypothèse). Ce ratio peut alors s'écrire:

W = \frac{R^4 \hat{\sigma}^2}{C^2 S^2}  = \frac{b^2}{S^2} = \frac{(a^T X)^2}{S^2} = \frac{\sum_{i=1}^n a_i X_i}{\sum_{i=1} ^n (X_i - \overline{X})^2}

Avec R^2 = m^T V^{-1} m et C^2 = m^T V^{-1} V^{-1} m.

Le développement de l'équation associée à W permet de déterminer la formule  générale des a_i:

a^T = (a_1, \cdots, a_n) = \frac{m^T V^{-1}}{\sqrt{m^T V^{-1} V^{-1} m}}

\bullet Exemple:

Test de Shapiro-Wilk

Soit la variable aléatoire X suivante:

\begin{tabular}{|c|} \hline X \\ \hline 0.9754 \\ \hline 1.2699 \\ \hline 1.4189 \\ \hline 1.5761 \\ \hline 2.7850 \\ \hline 4.2176 \\ \hline 4.8538 \\ \hline 5.4688 \\ \hline 6.3236 \\ \hline 7.9221 \\ \hline 8.0028 \\ \hline 8.1472 \\ \hline 9.0579 \\ \hline 9.1338 \\ \hline 9.1574 \\ \hline 9.5717 \\ \hline 9.5751 \\ \hline 9.5949 \\ \hline 9.6489 \\ \hline 9.7059 \\ \hline \end{tabular}

Ci-dessous, le QQplot construit à partir de X permet de voir que visuellement l'échantillon ne suit pas une loi normale. Nous affichons également la pente théorique à adhérer afin de s'approcher d'une loi normale  (droite en bleu) et les points projetés observés (point noir).

add

Nous avons donc  \overline{X} = 6.42039. Afin de déterminer la statistique de test de Shapiro-Wilk associée à notre exemple, nous allons présenter les différentes phases de calcul de la formule.

Notons donc le vecteur:

D = (X_{20} - X_1, X_{19} - X_2, \ldots, X_{11} - X_{10})
= (8.7305, 8.379, 8.176, 7.999, 6.7867, 4.9398, 4.28, 3.3891, 1.8236, 0.0807)

et \sum_{i = 1} ^n (X_i ^j - \overline{X ^j}) ^2 = 206.8246.

Afin de récupérer les coefficients a_i, nous nous reportons à la table de Shapiro-Wilk pour n = 20 (nombre d'observations dans notre exemple), nous trouvons alors le vecteur:

a = (0.4734, 0.3211, 0.2565, 0.2085, 0.1686, 0.1334, 0.1013, 0.0711, 0.0422, 0.0140)

Nous avons donc la statistique de test:

W = \frac{(0.4734 \times 8.7305 + \ldots 0.0807 \times 0.0140) ^2}{206.8246} = 0.837162

ce qui nous donne une p-value valant 0.0033 impliquant le rejet de l'hypothèse  H_0. Nous en concluons que l'échantillon étudié en guise d'exemple ne suit pas une loi normale, ce qui est en adéquation avec ce que l'on a pu observer sur le QQplot ci-dessus.

Test de Shapiro-Francia

Test de Ryan-Joiner

\bullet Application informatique:

Procédure SAS: http://www.stattutorials.com/SAS/TUTORIAL-PROC-UNIVARIATE.htm

Soit l'exemple suivant:

data E;
input X;
cards;
0.9754
1.2699
1.4189
1.5761
2.7850
4.2176
4.8538
5.4688
6.3236
7.9221
8.0028
8.1472
9.0579
9.1338
9.1574
9.5717
9.5751
9.5949
9.6489
9.7059
;
run;

Package et fonction R: http://stat.ethz.ch/R-manual/R-patched/library/stats/html/shapiro.test.html

Soit l'exemple suivant:

X = c(0.9754,1.2699,1.4189,1.5761,2.7850,4.2176,
4.8538,5.4688,6.3236,7.9221,8.0028,8.1472,9.0579,
9.1338,9.1574,9.5717,9.5751,9.5949,9.6489,9.7059)

\bullet Bibliographie:

- Data mining et statistiques décisionnelles: l'intelligence des données de Stéphane Tufféry

- An Analysis of Variance Test for Normality de Samuel Shapiro et Martin Bradbury Wilk

- Tests de normalité: Techniques empiriques et tests statistiques de Ricco Rakotomalala

- Comparison of Common Tests for Normality de Johannes Hain

- Algorithm AS 181: The W Test for Normality de J. P. Royston

- A modification of the test of Shapiro and Wilk for normalité de M. Mahibbur Rahman et Z. Govindarajulu

- Normal probability plots and tests for normality de Thomas A. Ryan, Barbara Ryan et Brian L. Joiner

- Conditions d’application des méthodes statistiques paramétriques : applications sur ordinateur de R. Glele Kakaï, E. Sodjinou et N. Fonton.