Les séries temporelles

add.png

\bullet Présentation:

Les séries temporelles ou chronologiques sont des processus stochastiques définies par une variable aléatoire continue X = (X_1, \cdots, X_T). Les séries temporelles sont accompagnées également d’une batterie d’outils statistique visant à modéliser et fournir des indicateurs sur leur nature.

Les séries temporelles ont connu leur premier fait d’arme au début du vingtième siècle. C’est en 1970, grâce à l’oeuvre de George Box et Gwilym Jenkins qu’elles deviennent très populaires et régulièrement utilisées.

Etudier les séries temporelles revient à étudier les trois principaux éléments suivants,

– la stationnarité soit l’évolution de la structure de la série temporelle en fonction du temps,

– la saisonnalité soit le cycle de reproduction de la série temporelle au travers de l’autocorrélation,

– la tendance globale de la série.

Dans l’objectif de modéliser une série temporelle aussi bien pour étudier sa reproductibilité qu’à des fins prévisionnelles et prédictives, une large variété de modèles ont vus le jour et dont les principaux sont:

– le processus AR (Auto Regressive) basé sur les valeurs antérieurs,

– le processus MA (Moving-Average) qui est souvent confondu avec les moyennes mobiles (ou Moving Average, le trait d’union faisant toute la différence) classiques, qui n’a rien à voir et qui est basé sur les résidus,

– le processus ARMA (Auto Regressive Moving Average) limité aux séries stationnaires et linéaires et mixant les caractéristiques du modèle AR et du modèle MA,

– le processus ARIMA (Auto Regressive Integrated Moving Average) qui est une généralisation du modèle ARMA aux séries non stationnaires et linéaires,

– le processus SARIMA (Seasonnal AutoRegressive Integrated Moving Average)  qui est une extension du modèle ARIMA pour les séries avec saisonnalité fixe,

– le processus ARFIMA (Autoregressive Fractionally Integrated Moving Average) qui généralise le modèle ARIMA aux processus à mémoire longue,

– le processus VAR (Vector Auto Regressive) qui permet la gestion de plusieurs séries endogènes linéaires,

– les processus ARCH (Auto Regressive Conditionnal Heteroscedasticity), GARCH (Generalized Regressive Conditionnal Heteroscedasticity) et leur multiple variante (NGARCH, IGARCH, EGARCH, GARCH-M, QGARCH, GJR-GARCH, TGARCH, fGARCH et COGARCH) pour la gestion de phénomènes non linéaires avec volatilité.

Enfin, les limites des séries chronologiques ont donné naissance aux modèles combinant séries chronologiques et chaînes de Markov, représentant une nouvelle branche en plein essor.

\bullet Les modèles:

De manière générale, nous poserons B, l’opérateur retard, de formule,

B \cdot X_t = X_{t - 1}

Et,

\Phi(B) = 1 - \phi_1 B - \cdots - \phi_l B ^l

Le modèle AR(p)

Le modèle:

Un processus autorégressif d’ordre p est un processus stationnaire tel que,

X_t = \sum_{i = 1} ^p \phi_i X_{t - i} + \epsilon_t, \forall t \in Z

En passant par l’opérateur retard, nous pouvons réécrire le processus sous la forme,

\Phi(B) X_t = \epsilon_t

Estimation des coefficients (\phi_1, \cdots, \phi_p):

La résolution du modèle se fait au travers des équations de Yule-Walker,

\begin{pmatrix} \phi_1 \\ \vdots \\ \phi_p \end{pmatrix} = \begin{pmatrix} 1 & \rho(1) & \rho(2) & \rho(3) & \ldots & \rho(p-1) \\ \rho(1) & 1 & \rho(1) & \rho(2) & \ldots & \rho(p - 2) \\ \rho(2) & \rho(1) & 1 & \rho(1) & \ldots & \rho(p - 3) \\ \vdots & \vdots & \vdots & \vdots & \cdots & \vdots \\ \rho(p - 1) & \rho(p - 2) & \rho(p - 3) & \rho(p - 4) & \ldots & 1 \end{pmatrix} ^{-1} \begin{pmatrix} \rho(1) \\ \vdots \\ \rho(p) \end{pmatrix}

Avec \rho(i) le coefficient d’autocorrélation.

Le modèle moyenne mobile MA(q)

Le modèle:

Un processus moyenne mobile d’ordre q est tel que,

X_t = \epsilon_t - \sum_{i = 1} ^q \theta_i \epsilon_{t - i}, \forall t \in Z

En passant par l’opérateur retard, nous pouvons réécrire le processus sous la forme,

X_t = \Theta ^{-1} (B) \epsilon_t

Estimation des coefficients (\theta_1, \cdots, \theta_q):

L’estimation des paramètres d’un modèle MA(q) nécessite la transformation du processus en un modèle AR(\infty) en inversant la seconde forme ci-dessus,

X_t = \Theta (B) \epsilon_t = - \sum_{i = 1} ^{\infty} \pi_i X_{t - i} + \epsilon_t, \sum_{i = 1} ^{\infty} | \pi_i | < + \infty

Ce système se résout alors par ajustement itératif non linéaire.

Le modèle ARMA(p,q)

Le modèle:

Un processus moyenne mobile autorégressive d’ordre (p, q) est un processus stationnaire tel que,

X_t = \theta_0 + \epsilon_t + \sum_{i = 1} ^p \phi_i X_{t - i} - \sum_{i = 1} ^q \theta_i \epsilon_{t - i}

En passant par l’opérateur retard, nous pouvons réécrire le processus sous la forme,

\Phi(B) X_t = \theta_0 + \Theta(B) \epsilon_t

Estimation des coefficients (\phi_1, \cdots, \phi_p, \theta_1, \cdots, \theta_q):

Plusieurs familles de méthodes existent pour l’estimation des coefficients:

– celles basées sur la fonction d’autocovariance,

– celles basées sur le pseudo-maximum de la vraisemblance selon l’un des trois critères MCC (moindres carrés conditionnel), MCN (moindres carrés non conditionnel) ou MV (maximum de vraisemblance),

– celles basées sur la résolution directe des équations de vraisemblance approchées,

– celles basées sur des approches itératives d’estimation.

Le modèle ARIMA(p,d,q)

Le modèle:

Un processus moyenne mobile intégrée autoregressive d’ordre (p,d,q) est un processus non stationnaire tel que,

X_t = X_{-S + t} + \sum_{i = 1} ^{t + S - 2} \epsilon_{t - i}, t \leq -S + 2, S \in N

En passant par l’opérateur retard, et en posant \Delta ^d X_t = (1 - B) ^d la différence d’ordre d de X_t, nous pouvons réécrire le processus sous la forme,

\Phi(B) \Delta ^d X_t = \Theta(B) \epsilon_t

Estimation des coefficients (\phi_1, \cdots, \phi_p, \theta_1, \cdots, \theta_q):

L’Estimation des coefficients se fait de la même manière que pour un processus ARMA(p + d,q).

\bullet La saisonnalité:

L’Autocorrélation d’une série temporelle est en générale le premier point à étudier.

Elle se base sur le coefficient de corrélation de Pearson et doit son nom au fait qu’elle consiste à étudier le lien entre la série elle-même et sa version décalée de h périodes. Nous parlons alors d’autocorrélation d’ordre h. Elle permet donc de visualiser les liens entre les différents temps de X et de détecter la présence de saisonnalité.

L’Autocorrélation partielle représente la corrélation entre X_t et X_{t+h} une fois que nous avons retiré l’influence de X entre ces deux temps. Elle permet ainsi de détecter les fausses corrélations et d’améliorer l’information associée à la détection de vrai ou fausse saisonnalité.

– L’Autocovariance d’ordre h:

cov(X_t, X_{t+h}) = \frac{1}{T} \sum_{t = 1} ^{T-h} (X_t - \overline{X}) (X_{t + h} - \overline{X})

– L’Autocorrélation d’ordre h:

\rho(h) = \frac{cov(X_t, X_{t+h})}{V(X_t)}

– L’Auto-corrélation partielle d’ordre h:

p_h = \frac{det(P_h ^*)}{det(P_h)}

Avec,

P_h ^* = \begin{pmatrix} 1 & \rho(1) & \rho(1) & \rho(1) & \cdots & \rho(1) & \rho(1) \\ \rho(1) & 1 & \rho(2) & \rho(2) & \cdots & \rho(2) & \rho(2) \\ \rho(2) & \rho(1) & 1 & \rho(3) & \cdots  & \rho(3) & \rho(3) \\ \rho(3) & \rho(2) & \rho(1) & 1 & \cdots & \rho(4) & \rho(4) \\ \vdots & \vdots & \vdots & \vdots & \cdots & \vdots & \vdots \\ \rho(h-2) & \rho(h-3) & \rho(h-4) & \rho(h-5) & \cdots & 1 & \rho(h-1) \\ \rho(h-1) & \rho(h-2) & \rho(h-3) & \rho(h-4) & \cdots & \rho(1) & \rho(h) \end{pmatrix}

NB: P_1 ^* = \rho(1)

P_h = \begin{pmatrix} 1 & \rho(1) & \rho(2) & \rho(3) & \cdots & \rho(h-1) \\ \rho(1) & 1 & \rho(3) & \rho(2) & \cdots & \rho(h-2) \\ \rho(2) & \rho(3) & 1 & \rho(1) & \cdots & \rho(h-3) \\ \vdots  & \vdots & \vdots & \vdots & \ldots & \vdots \\ \rho(h-2) & \rho(h-3) & \rho(h-4) & \rho(h-5) & \cdots & \rho(1) \\ \rho(h-1) & \rho(h-2) & \rho(h-3) & \rho(h-4) & \cdots & 1 \end{pmatrix}

NB: P_1 = 1

– Les tests pour la détection de la saisonnalité:

Tout une gamme de tests existe pour déterminer si les autocorrélations ou autocorrélations partielles sont significatives. Parmi eux nous avons les tests du Portemanteau (test de Box-Pierce et de Ljung-Box), le test de Durbin-Watson et celui de Breusch-Godfrey.

\bullet La tendance:

La tendance (que nous retrouvons souvent dans la littérature sous le terme anglosaxon: trend) représente l’orientation globale de la série X de son premier temps mesuré au dernier.

Elle s’étudie généralement via une régression linéaire pour sa forme générale et via une analyse de variance (ANOVA) pour sa forme regroupée (par mois, jour, etc).

Les informations obtenues suite à l’analyse de la tendance de la série temporelle permettent entre autres de déterminer si oui ou non elle se délocalise dans le temps et ainsi d’adapter le modèle.

\bullet La stationnarité:

La notion de stationnarité d’un processus joue un rôle primordiale dans l’étude des séries chronologiques. Elle consiste à observer l’évolution du processus X au fur et à mesure que t \in [1, T] croît. Si elle ne change pas alors nous parlons de processus stationnaire, dans le cas contraire nous parlons de processus non stationnaire.

Enfin, étudier la stationnarité d’un processus revient également à étudier si sa racine est unitaire.

– Le test de Dickey-Fuller augmenté

Présentation:

Le test de Dickey-fuller augmenté , développé par David A. Dickey et Wayne A. Fuller en 1979, permet de détecter la non-stationnarité d’une série temporelle X = (X_1, \cdots, X_T). Il s’agit d’une version étendue du test de Dickey-Fuller et qui était restreint aux processus AR(1).

Le test:

L’idée du test est d’étudier le modèle linéaire,

Z ^1 = X|_{t \in [K, T]} + 1 + t_{\in [K, T]} + (Z ^2, \cdots, Z ^K)

Avec,

K = E[(T - 1) ^{\frac{1}{3}}] + 1, le nombre de décalage à considérer et E[.] désigne ici la partie entière et non l’espérance,

\mathbf{Z} = (X|_{t \in [K,T]}, X|_{t \in [K-1,T-1]}, \cdots, X|_{t \in [1,T-K+1]}), la matrice des décalages successifs de X,

1, le vecteur unitaire de longueur (T - K).

Nous appliquons alors une régression linéaire sur ce modèle et retenons le coefficient associé à X|_{t \in [K, T]}, que nous noterons \beta, ainsi que son écart-type sd_{\beta}.

La statistique de test de Dickey-Fuller augmenté est alors,

ADF = \frac{\beta}{sd_{\beta}}

Elle suit une loi de Dickey-Fuller et l’hypothèse H_0 est: « \beta = 0 / le processus est non-stationnaire ».

– Le test de Kwiatowski-Phillips-Schmidt-Shin

Présentation:

Le test de Kwiatowski-Phillips-Schmidt-Shin, publié en 1992 par D. Kwiatkowski, P. C. B. Phillips, P. Schmidt et Y. Shin et noté test KPSS, est un test de stationnarité.

Le test:

Le calcul de la statistique de test se fait selon les trois étapes suivantes:

– Étape 1: Appliquer une régression linéaire de X sur t, \in [1,T] et en sortir le vecteur des résidus \epsilon_t, t \in [1, T].

– Étape 2: Calculer le vecteur S des sommes cumulées des résidus et de forme,

S = (\epsilon_1, \sum_{t = 1} ^2 \epsilon_t, \cdots, \sum_{t = 1} ^T \epsilon_t)

– Étape 3: Calculer la statistique de test KPSS,

KPSS = \frac{\frac{1}{T ^2} \sum_{t = 1} ^T S_t}{\frac{1}{T} \sum_{t = 1} ^T \epsilon_t}

La statistique de test suit une loi de Kwiatowski-Phillips-Schmidt-Shin et l’hypothèse H_0 est: « \frac{1}{T ^2} \sum_{t = 1} ^T S_t = 0 / le processus est stationnaire ».

– Les autres tests pour la détection de la stationnarité:

Toute une série de tests existent dans l’objectif de détecter la stationnarité d’une série temporelle.

Parmi eux, nous relevons le test de Phillips-Perron, le test de Schmidt-Philipps, le test de Elliott-Rothenberg-Stock, le test de Leybourne-McCabe, le test de Hasza-Fuller  (pour le cas d’un SARIMA), le test de Osborn-Chui-Smith-Birchenhall (pour le cas d’un SARIMA), le test de Hylleberg-Engle-Granger-Yoo (pour le cas d’un SARIMA) et le test de Franses (pour le cas d’un SARIMA).

\bullet Exemple:

Soit l’échantillon suivant,

add1

Dans un premier temps, dessinons la série de données,

add1

Visuellement, notre série semble présenter une tendance décroissante ainsi qu’une saisonnalité avec ses deux pics à T = 6 et T = 16. En terme d’autocorrélation, nous pouvons constater qu’après un pic fort, la série chute pour finalement remonter petit à petit.

Étude de saisonnalité:

Nous fixons un pas de h = 1 et avons \overline{X} = 3.991.

Dans un premier temps, calculons l’autocovariance,

cov(X_t, X_{t+1}) = \frac{1}{20} \sum_{t = 1} ^{19} (X_t - 3.991) \times (X_{t+1} - 3.991)

= \frac{1}{20} \times (8.805951 + 17.526621 + \cdots + 85.773141 + 51.154461)

= \frac{1}{20} \times 275.5481

= 13.7774

cov(X_t, X_{t + 2}) = \frac{1}{20} \sum_{t = 1} ^{18} (X_t - 3.991) \times (X_{t+2} - 3.991)

= \frac{1}{20} \times (20.393571 + 4.404351 + \cdots + 55.979631 + 74.798121)

= \frac{1}{20} \times (-76.07523)

= -3.803762

\vdots

cov(X_t, X_{t + 13}) = \frac{1}{20} \sum_{t = 1} ^7 (X_t - 3.991) \times (X_{t+13} - 3.991)

= \frac{1}{20} \times (-7.679199 + 4.156761 + \cdots - 133.350849 + 22.100811)

= \frac{1}{20} \times (-102.7504)

= -5.13752

Maintenant que nous avons définit la covariance pour les différentes combinaisons possibles, définissons,

cov(X_t, X_{t+0}) = var(X_t) = \frac{1}{20} \times \sum_{t = 1} ^{20} (X_t - 3.991) ^2 = \frac{955.6434}{20} = 47.78217

Nous pouvons calculer l’autocorrélation en chaque instant mesurable,

\rho(1) = \frac{cov(X_t, X_{t+1})}{var(X_t)} = \frac{13.7774}{47.78217} = 0.2883377

\rho(2) = \frac{cov(X_t, X_{t+2})}{var(X_t)} = \frac{-3.803762}{47.78217} = -0.07960631

\vdots

\rho(13) = \frac{cov(X_t, X_{t+13})}{var(X_t)} = \frac{-5.13752}{47.78217} = -0.1075196

Nous obtenons alors l’autocorrélogramme suivant,

add1.png

Ce graphe nous permettent de retrouver nos deux pics saisonniers. Le second, notamment, semble corrélé aux temps qui le précèdent et celui qui le suivent. Les temps intermédiaires sont, quand à eux, corrélés les uns aux autres.

Il nous reste alors à calculer l’autocorrélation partielle.

– Pour p_1. Nous avons,

P_1 = 1 \Rightarrow |P_1| = 1

Et,

P_1 ^* = \rho(1) = 0.2883377 \Rightarrow |P_1 ^*| = 0.2883377

Par conséquent,

p_1 = \frac{0.2883377}{1} = 0.2883377

– Pour p_2. Nous avons,

P_2 = \begin{pmatrix} 1 & \rho(1) \\ \rho(1) & 1 \\ \end{pmatrix} = \begin{pmatrix} 1 & 0.2883377 \\ 0.2883377 & 1 \\ \end{pmatrix} \Rightarrow |P_2| = 0.9168614

Et,

P_2 ^* = \begin{pmatrix} 1 & \rho(1) \\ \rho(1) & \rho(2) \\ \end{pmatrix} = \begin{pmatrix} 1 & 0.2883377 \\ 0.2883377 & -0.07960631 \\ \end{pmatrix} \Rightarrow |P_2 ^*| = -0.1627449

Par conséquent,

p_2 = \frac{-0.1627449}{0.9168614} = -0.1775022

– Pour p_3. Nous avons,

P_3 = \begin{pmatrix} 1 & \rho(1) & \rho(2) \\ \rho(1) & 1 & \rho(1) \\ \rho(2) & \rho(1) & 1 \\ \end{pmatrix} = \begin{pmatrix} 1 & 0.2883377 & -0.07960631 \\ 0.2883377 & 1 & 0.2883377 \\ -0.07960631 & 0.2883377 & 1 \\ \end{pmatrix}

\Rightarrow |P_3| = 0.814441

Et,

P_3 ^* = \begin{pmatrix} 1 & \rho(1) & \rho(1) \\ \rho(1) & 1 & \rho(2) \\ \rho(2) & \rho(1) & \rho(3) \\ \end{pmatrix} = \begin{pmatrix} 1 & 0.2883377 & 0.2883377 \\ 0.2883377 &1 & -0.07960631 \\ -0.07960631 & -0.07960631 & 0.015 \\ \end{pmatrix}

\Rightarrow |P_3 ^*| = 0.08556691

Par conséquent,

p_3 = \frac{0.08556691}{0.814441} = -0.1050621

– Pour p_4. Nous avons,

P_4 = \begin{pmatrix} 1 & \rho(1) & \rho(2) & \rho(3) \\ \rho(1) & 1 & \rho(3) & \rho(2) \\ \rho(2) & \rho(3) & 1 & \rho(1) \\ \rho(3) & \rho(2) & \rho(1) & 1 \end{pmatrix} = \begin{pmatrix} 1 & 0.2883377 & -0.07960631 & 0.015 \\ 0.2883377 & 1 & 0.015 & -0.07960631 \\ -0.07960631 & 0.015 & 1 & 0.2883377 \\ 0.015 & -0.07960631 & 0.2883377 & 1 \\ \end{pmatrix}

\Rightarrow |P_4| = 0.823916

Et,

P_4 ^* = \begin{pmatrix} 1 & \rho(1) & \rho(1) & \rho(1) \\ \rho(1) & 1 & \rho(2) & \rho(2) \\ \rho(2) & \rho(1) & 1 & \rho(3) \\ \rho(3) & \rho(2) & \rho(1) & \rho(4) \\ \end{pmatrix} = \begin{pmatrix} 1 & 0.2883377 & 0.2883377 & 0.2883377 \\ 0.2883377 & 1 & -0.07960631 & -0.07960631 \\ -0.07960631 & 0.2883377 & 1 & 0.015 \\ 0.015 & -0.07960631 & 0.2883377 & -0.140 \\ \end{pmatrix}

\Rightarrow |P_4 ^*| = -0.1805203

Par conséquent,

p_4 = \frac{-0.1805203}{0.823916} = -0.2191004

\vdots

En continuant ainsi pour les différents temps, nous pouvons tracer l’autocorrélogramme partiel,

add1.png

Le graphe des autocorrélations partielles présentent une nouvelle information et met en évidence que nos corrélations entre les temps intermédiaires aux deux saisonnalités sont en réalité peut avérées et qu’il s’agisse plus d’une corrélation en forme de saut que d’une corrélation entre temps adjacents.

Étude de la stationnarité:

Étudions maintenant la stationnarité de X. Pour se faire, nous allons donc utiliser le test de Dickey-Fuller et le test KPSS.

– Le test de Dickey-Fuller augmenté:

Dans un premier temps, construisons la matrice \mathbf{Z} de taille,

(\sharp \lbrace X \rbrace - 1) \times (E[(\sharp \lbrace X \rbrace - 1) ^{\frac{1}{3}}] + 1) = 19 \times  (E[2.668402] + 1) = 19 \times (2 + 1) = 19 \times 3, (K = 3)

Nous obtenons alors,

\mathbf{Z} = \begin{pmatrix} X_3 & X_2 & X_1 \\ \ldots & \ldots & \ldots \\ X_{20} & X_{19} & X_{18} \end{pmatrix} = \begin{pmatrix} -4.77 & 3.62 & -0.45 \\ \ldots & \ldots & \ldots \\ 0.98 & 3.54 & -3.89 \end{pmatrix}

La régression linéaire de Z ^1 sur le modèle,

X|_{t \in [3,20]} + 1 + t_{\in [3,20]} + (Z ^2, Z ^3)

, nous donne le coefficient \beta associé à X|_{t \in [3,20]} et son écart-type suivants,

\beta = -1.104306, sd_{\beta} = 0.5171055

Nous pouvons alors calculer la statistique du test de Dickey-Fuller augmenté,

Z = \frac{-1.104306}{0.5171055} = -2.135552

Reportée à la table de Dickey-Fuller, nous obtenons une p-valeur de 0.3196>5\%. Nous ne pouvons pas rejeter H_0 soit que X est non stationnaire.

– Le test KPSS:

Dans un premier temps, déterminons les résidus issues de la régression linéaire de X sur t \in [1,T]. Nous obtenons alors le vecteur,

e = (-2.0628571, -1.9587669, \cdots, -1.4151429)

Nous avons alors,

S = (e_1, \sum_{t = 1} ^2 e_t, \sum_{t = 1} ^3 e_t, \dots, \sum_{t = 1} ^n e_t)

= (-2.0628571, -2.0628571 - 1.9587669, -2.0628571 - 1.9587669 + 2.2153233, \dots)

= (-2.628571, -4.0216241, -1.8063008, \cdots, 0)

Reste à calculer le numérateur et le dénominateur de la statistique du test KPSS.

Nous avons pour le numérateur,

\frac{\sum_{t = 1} ^{20} S_t ^2}{20 ^2} = \frac{951.8915}{400} = 2.379729

, et au dénominateur,

\frac{\sum_{t = 1} ^{20} e_t ^2}{20} = \frac{751.4778}{20} = 37.57389

Nous pouvons enfin calculer la statistique du test KPSS,

Z = \frac{2.379729}{37.57389} = 0.06333465

Si nous reportons cette valeur à la table de la loi de KPSS, nous obtenons une p-valeur de 0.1>5\% mais pas à 10\% soit que nous pouvons rejeter H0 au second seuil et que X n’est pas stationnaire.

Étude de la tendance:

Intéressons-nous ensuite à la tendance de notre série.

Pour cela nous allons recourir à une régression linéaire simple de X sur le temps t \in [1,20] sur X. Nous obtenons alors la modélisation suivante,

add1.png

Nous en déduisons que nous sommes bien en présence d’une tendance décroissante.

Modélisation par un processus AR(2):

Tentons de modéliser X par un processus auto-régressif d’ordre p = 2. L’écriture de ce modèle est,

X_t = \sum_{i = 1} ^2 \phi_i X_{t - i} + \epsilon_t = \phi_1 X_{t - 1} + \phi_2 X_{t - 2} + \epsilon_t

Afin de déterminer les coefficients (\phi_1, \phi_2) nous passons par les équations de Yule-Walker,

\begin{pmatrix} \phi_1 \\ \phi_2 \end{pmatrix} = \begin{pmatrix} 1 & \rho(1) \\ \rho(1) & 1 \end{pmatrix} ^{-1} \times \begin{pmatrix} \rho(1) \\ \rho(2) \end{pmatrix}

= \begin{pmatrix} 1 & 0.2883377 \\ 0.2883377 & 1 \end{pmatrix} ^{-1} \times \begin{pmatrix} 0.2883377 \\ -0.07960631 \end{pmatrix}

= \begin{pmatrix} 1.0906774 & -0.3144834 \\ -0.3144834 & 1.0906774 \end{pmatrix} \times \begin{pmatrix} 0.2883377 \\ -0.07960631 \end{pmatrix}

= \begin{pmatrix} 0.3395183 \\ -0.1775022 \end{pmatrix}

Nous obtenons ainsi le modèle,

X_t = 0.3395183 \times X_{t - 1} - 0.1775022 \times X_{t - 2} + \epsilon_t

Nous pouvons alors comparer les valeurs réels de X (en noir) versus celles prédites par le modèle (en rouge),

add1

\bullet Application informatique:

– Procédure SAS: https://support.sas.com/documentation/onlinedoc/ets/132/arima.pdf

– Fonction R: http://stat.ethz.ch/R-manual/R-devel/library/stats/html/arima.html

\bullet Bibliographie:

– Séries chronologiques de Jean-Jacques Droesbeke, Bernard Fichet, Philippe Tassi

– Econométrie Appliquée. Séries Temporelles de Christophe HURLIN

– Introduction to Statistical Times Series de Wayne A. Fuller

– Stationarity Issues in Times Series de David A. Dickey

– Testing the nulle hypothesis of stationarity against the alternative of a unit root de D. Kwiatkowski, P. C. B. Phillips, P. Schmidt et Y. Shin

– Modèles de prévision Séries temporelles de Arthur Charpentier

– Le site web: http://www.minitab.com/fr-fr/Published-Articles/ARIMA%C2%A0–Comment-analyser-des-donn%C3%A9es-de-s%C3%A9rie-chronologique-diff%C3%A9remment%C2%A0-/

Les tests pour la détection d’outliers

add.png

\bullet Introduction:

Les outliers ou observations aberrantes sont des valeurs d’une variable continue X douteuses au sens d’un écart fort avec le reste de la distribution de X. L’intérêt au problème des outliers débute en 1750 avec les travaux de Roger Joseph Boxcovitch, Pierre Simon de Laplace et Adrien-Marie Legendre. 

Deux approches sont à suggérer pour le traitement des outliers: les conserver en adaptant le choix de l’outil statistique ou alors leur suppression. La principale motivation à leur détection est qu’ils peuvent étirer la moyenne ou l’écart-type à tort apportant ainsi un biais au calcul des estimateurs d’intérêt.

La plupart du temps les outliers ont pour origine une erreur de saisie ou de mesure ce qui privilégierait leur suppression, néanmoins le manque de certitude sur leur origine fait qu’il n’y a pas vraiment de méthode universelle car un outlier peut également être une réelle valeur et simplement inattendue. C’est d’ailleurs sur ce point précis que la communauté scientifique s’oppose lors de leur traitement.

Toute une batterie de tests et de critères statistiques a été conçue afin de déterminer au sens statistique du terme si une observation est un outlier ou non, le plus populaire de tous restant le critère décisionnel de Tukey.

add1

\bullet Les différents tests:

– Le critère de Tukey

Qui et quand: John Wilder Tukey, 1977

Le critère:

L’Outil possède l’avantage de pouvoir déterminer un seuil à partir duquel les observations, en fonction de leur valeur, doivent être considérées comme aberrante relativement à la distribution de X. Permettant ainsi la détection d’un ou plusieurs outliers simultanément.

L’algorithme décisionnel du critère de Tukey est, \forall i \in [1,n],

si X_i < Q_1 - 1.5 IQR ou X_i > Q_3 - 1.5 IQR

, alors l’observation i est un outlier.

Nous noterons Q_1 le premier quartile (soit le quantile à 25\%), Q_3 le troisième quartile (soit celui à 75\%) et IQR l’interquartile. Leurs formules respectives sont,

Q_1 = X_{E[\frac{n+1}{4}]} + 0.75 (X_{E[\frac{n+1}{4}] + 1} - X_{E[\frac{n+1}{4}]}),

Q_3 = X_{3 E[\frac{n+1}{4}]} + 0.75 (X_{3 E[\frac{n+1}{4}] + 1} - X_{3 E[\frac{n+1}{4}]}),

IQR = Q_3 - Q_1.

E[.] désigne la partie entière de ..

– Le test Q de Dixon

Qui et quand: Wilfrid Joseph Dixon, 1951.

Le test:

Le test Q de Dixon se base principalement sur la distance entre X_i, l’observation suspectée d’être aberrante, et celle qui la précède afin de déterminer s’il s’agit d’un outlier ou non.

Supposons X trié, la formule du test Q de Dixon est,

– si X_i \geq mediane(X) alors Q = \frac{X_{i + 1} - X_i}{X_n - X_1},

– si X_i \leq mediane(X) alors Q = \frac{X_i - X_{i - 1}}{X_n - X_1}.

Nous reportons ensuite Q à la table de Dixon ci-dessous selon un niveau de confiance \alpha fixé,

add

Si Q > Q_{\alpha \%}, alors l’observation i peut être considérée comme un outlier.

– Le critère de Chauvenet

Qui et quand: William Chauvenet, 1863.

Le critère:

Le critère de Chauvenet se base sur une approche stochastique pour la détermination d’outliers et considère les observations les unes après les autres. La procédure pour départager tous les outliers consiste à supprimer les observations en partant de la valeur maximale ou minimale et de relancer le test à chaque nouvelle observation à vérifier.

Pour une observation X_i suspectée d’être un outlier, le critère de Chauvenet est,

si n \cdot erfc (\frac{|X_i - \overline{X}|}{\sigma_X}) < 0.5 alors X_i est un outlier

, où que erfc est la fonction d’erreur complémentaire.

– Le test de Grubbs

Qui et quand: Frank E. Grubbs, 1969.

Le test:

Le test de Grubbs, également connu sous le nom de test du maximum des résidus normalisés, recherche à déterminer si la valeur maximale ou minimale de X est un outlier. Lorsque nous suspectons un lot d’observations, il faudra se tourner vers sa généralisation: le test de Tietjen-Moore.

La formule du test de Grubbs est,

G_{min} = \frac{\overline{X} - min_{i \in 1, \cdots, n} (X)}{s} et G_{max} = \frac{max_{i \in 1, \cdots, n} (X) - \overline{X}}{s}

Le critère de décision, selon un seuil de confiance \alpha fixé et permettant de déterminer si la valeur maximale et/ou la valeur minimale sont des outliers est,

Si G_{min}, G_{max} > \frac{n - 1}{\sqrt{n}} \sqrt{\frac{t_{\frac{\alpha}{n},n-2} ^2}{n - 2 + t_{\frac{\alpha}{n},n - 2} ^2}} alors ce sont des outliers

, t_{\delta,df} représentant la fonction quantile associée à la distribution de Student pour le seuil de confiance \delta et pour df degrés de liberté.

– Le test de Tietjen-Moore

Qui et quand: Gary Tietjen et Roger Moore, 1972

Le test:

Le test de Tietjen-Moore est une généralisation de celui de Grubbs permettant de détecter un lot d’outliers et plus uniquement si la valeur maximale et minimale en sont. A noter que les valeurs testées doivent se succéder au niveau de la distribution. De plus, le test de Tietjen-Moore reste très dépendant du paramètre I à fixer et représentant le nombre d’observations à tester, ce qui peut nuire à la qualité de la détection s’il est mal paramétré.

Trois version du test de Tietjen-Moore existe en fonction de ce que nous cherchons à déterminer (les données seront supposées triées),

– Pour la ou les valeurs minimale(s) de la première à la I-ème:

L_{I,min} = \frac{\sum_{i = 1} ^{n-I} (X_i - \overline{X_{-[1,I]}}) ^2}{\sum_{i = 1} ^n (X_i - \overline{X}) ^2}

– Pour la ou les valeurs maximale(s) de la I-ème à la n-ème:

L_{I,max} = \frac{\sum_{i = I + 1} ^n (X_i - \overline{X_{-[I,n]}}) ^2}{\sum_{i = 1} ^n (X_i - \overline{X}) ^2}

, où X_{-[.,.]} représente le vecteur X privé des éléments inclus dans l’intervalle [.,.],

– Pour la ou les valeurs minimale(s) et maximale(s) (approche symétrique):

E_I = \frac{\sum_{i = 1} ^{n - I} (z_i - \overline{z_{-[I,n]}}) ^2}{\sum_{i = 1} ^n (z_i - \overline{z_i}) ^2}

, où z est le vecteur X ordonné en fonction du score résiduel \forall i \in [1,n], |X_i - \overline{X}|.

Afin de déterminer la valeur seuil à partir de laquelle nous pouvons tester notre hypothèse il faut simuler selon une loi normale centrée-réduite un nombre L de vecteurs de taille identique n à celle de X et appliquer le test de Tietjen-Moore pour chacun. La valeur seuil correspond au quantile à 5\% de la distribution des tests réalisés.

Si la valeur L_{I,min}, L_{I,max}, E_I est inférieur à la valeur seuil obtenu alors nous pouvons conclure que les observations considérées sont des outliers.

– Le test généralisé de la déviation extrême de Student

Qui et quand: Bernard Rosner, 1983

Le test:

Le test généralisé ESD (Extrem Studentized Deviation) permet la détection simultanée d’un ou plusieurs outliers et ne nécessite aucun paramétrage.

La première étape pour pouvoir calculer le test généralisé ESD est de ranger par ordre croissant X et déterminer le vecteur Z tel que,

Z_1 = \frac{max_{i \in [1, n]} (| X_i - \overline{X} |)}{\sigma_X}

\forall i \in [2,n],

— Suppression de l’observation qui maximise max_I(|X_i - \overline{X_I}|) et mise à jour de l’ensemble des observations à considérer I et dépourvu d’elle.

— Calcul de Z_i = \frac{max_{i \in I} (| X_i - \overline{X_I} |)}{\sigma_{X_I}}

, où X_I représente le vecteur X dépourvu des observations dans l’intervalle I indiqué.

En seconde étape, il faut pour chaque observation du vecteur Z, calculer la statistique de test \lambda,

\forall i \in [1,n], \lambda_i = \frac{(n - i) t_{1 - \frac{\alpha}{2 (n - i + 1)},n - i - 1}}{\sqrt{(n - i - 1 + t_{1 - \frac{\alpha}{2 (n - i + 1)},n - i - 1} ^2) (n - i + 1)}}

, où t_{\delta,df} représente la fonction quantile associée à la distribution de Student pour le seuil de confiance \delta et pour df degrés de liberté pour un seuil \alpha fixé.

Si Z_i > \lambda_i alors l’observation i peut être considérée comme un outlier.

– Le test \tau modifié de Thompson

Qui et quand: David J. Thompson, 1950.

Le test:

Le test \tau modifié de Thompson repose sur un algorithme permettant de vérifier directement quelles observations au sein de X peuvent être considérées comme des outliers.

Il faut dans un premier temps calculer,

\forall i \in [1,n], \delta_i = | X_i - \overline{X}|

Et pour l’observation maximisant \delta_i, pour un seuil de confiance \alpha fixé, la comparer à la valeur de référence,

\tau = \frac{t_{\frac{\alpha}{2},n-2} (n - 1)}{\sqrt{n (n - 2 + t_{\frac{\alpha}{2},n-2} ^2)}}

, où t_{\beta,df} représente la fonction quantile de Student pour le seuil de confiance \beta et pour df degrés de liberté.

Si \delta_i > \tau \sigma_X alors l’observation peut être considérée comme un outlier. 

La procédure doit être réitérée jusqu’à ce que \delta_i \leq \tau \sigma_X en supprimant à chaque nouvelle itération l’observation détectée comme un outlier et en mettant à jour la moyenne et l’écart-type.

– Le critère de Peirce

Présentation: Benjamin Peirce, 1877

Le test:

Le critère de Peirce présente de nombreux avantages: plus rigoureux que la majorité des autres méthodes de détection, il permet de déterminer directement l’ensemble des outliers au sein de X et ne repose pas sur un critère arbitraire pour le rejet d’une observation. De plus, il reste très simple à mettre en oeuvre.

L’algorithme de Peirce est le suivant,

1- Déterminer \overline{X}, \sigma_X,

2- A l’aide de la table de Peirce, calculer R pour un outlier suspecté sur n observations de X,

3- Calculer | X_i - \overline{X} |_{max} = \sigma_X R,

4- Rechercher les observations tel que \forall i \in [1,n], | X_i - \overline{X} | > |X_i - \overline{X} |_{max}, celles répondant à ce critère pourront être classées comme outlier,

5- Si au moins une observation répond à ce critère alors nous la supprimons de X et recalculons | X_i - \overline{X} |_{max} pour les même valeurs de n, \overline{X}, \sigma_X, seul R change et doit être relevé sur la table de Peirce pour deux outliers suspectés sur n observations. L’Opération est à renouveler à chaque fois que nous détectons un nouvel outlier.

6- Si plus aucune observation ne répond au critère, alors il faut mettre à jour n, \overline{X}, \sigma_X et reprendre depuis l’étape une à 5 tant qu’à l’étape 4 au moins une observation est détectée comme outlier.

Ci-dessous, la table de Peirce pour la détermination de R,

add

\bullet Exemple:

Soit l’échantillon suivant,

add2

Le boxplot associé à X est,

add3.png

– Le critère de Tukey

Nous cherchons donc à calculer le seuil, selon le critère de Tukey, à partir duquel une observation peut être considérée comme un outlier.

Nous avons donc,

Q_1 = 3.1547 + 0.75 \times (0.35896 - 3.1457)

= 3.1547 + 0.75 \times 0.4349

= 3.1547 + 0.326175

= 3.480875

Q_3 = 9.3500 + 0.25 \times (9.4578 - 9.3500)

= 9.3500 + 0.25 \times 0.1078

= 9.3500 + 0.02695

= 9.37695

IQR = 9.37695 - 3.480875 = 5.896075

Par conséquent,

S_1 = Q_1 - 1.5 \times IQR

= 3.480875 - 1.5 \times 5.896075

= 3.480875 - 8.44112

= - 5.363237

S_2 = Q_3 + 1.5 \times IQR

= 9.37695 + 1.5 \times 5.896075

= 9.37695 + 8.844112

= 18.22106

Nous n’avons pas d’observation pour laquelle sa valeur est inférieure à S_1, cependant nous avons une observation (X_{12}) dont la valeur est supérieure à S_2. Nous en déduisons que X_{12} est un outlier.

– Le test Q de Dixon

Nous cherchons à voir si la plus petite et la plus grande valeur de X sont des outliers. Si nous regardons la table de Dixon, la valeur seuil pour n = 20 et au risque de 5\% est: 0.450. Nous avons donc,

Q = \frac{2.1475 - 1.9593}{19.1245 - 1.9593} = \frac{0.1882}{17.1652} = 0.1882 < 0.450,

Q = \frac{19.1245 - 9.8308}{19.1245 - 1.9593} = \frac{9.2937}{17.1652} = 0.5414268 > 0.450.

Nous en déduisons, qu’au risque de 5\%, la valeur minimale (1.9593) n’est pas un outlier contrairement à la valeur maximale (19.1245) qui en est un.

– Le critère de Chauvenet

Nous cherchons à voir si la plus petite et la plus grande valeur de X sont des outliers. Dans un premier temps calculons la moyenne et l’écart-type,

\mu = 6.605515,

\sigma = 4.095388.

Nous pouvons maintenant calculer les éléments à soumettre au critère de Chauvenet. Nous avons alors,

\frac{1.9593 - 6.605515}{4.095388} = \frac{-4.646215}{4.095388} = -1.134499

\Rightarrow 20 \times erfc(|-1.134499|) = 20 \times 0.1086206525 = 2.1724130491

\Rightarrow 2.1724130491 > 0.5, nous en concluons que la valeur minimale n’est pas un outlier.

\frac{19.1245 - 6.605515}{4.095388} = \frac{12.51899}{4.095388} = 3.05685

\Rightarrow 20 \times erfc(|3.05685|) = 20 \times 0.0000153895 = 0.0003077894

\Rightarrow 0.0003077894 < 0.5, nous en concluons que la valeur maximale est un outlier.

– Le test de Grubbs

Nous cherchons à voir si la plus petite et la plus grande valeur de X sont des outliers. Dans un premier temps calculons la moyenne et l’écart-type,

\mu = 6.605515,

\sigma = 4.095388.

Nous avons donc,

G_{min} = \frac{6.605515 - 1.9593}{4.095388} = \frac{4.646215}{4.095388} = 1.134499

G_{max} = \frac{19.1245 - 6.605515}{4.095388} = \frac{12.51899}{4.095388} = 3.056851

Nous devons comparer ces deux valeurs au seuil,

\frac{20 - 1}{\sqrt{20}} \times \sqrt{\frac{t_{\frac{0.05}{20},20-2} ^2}{20 - 2 + t_{\frac{0.05}{20},20-2} ^2}} = \frac{19}{4.472136} \times \sqrt{\frac{t_{0.00125,18} ^2}{18 + t_{0.00125,18} ^2}}

= 4.248529 \times \sqrt{\frac{(-3.196574) ^2}{18 + (-3.196574) ^2}}

= 4.248529 \times \sqrt{\frac{10.21809}{28.21809}}

= 4.248529 \times 0.6017568

= 2.556381

Nous constatons que G_{min} < 2.556381 et G_{max} > 2.556381 et en concluons que seul la valeur maximale de X peut être considérée comme un outlier.

– Le test de Tietjen-Moore

Nous cherchons à déterminer si la valeur maximale et minimale sont des outliers, nous fixons donc I = 2. Nous appliquons donc la troisième version du test de Tietjen-Moore.

Dans un premier temps, nous avons les valeurs résiduelles,

\forall i \in[1,20], |X_i - \overline{X}| = (4.646215, 4.058315, \cdots, 1.482015, 1.340115)

Les valeurs maximales sont alors 4.646215, 12.518985, qui correspondante au minimum de X et à son maximum. Le vecteur z ordonné selon les valeurs résiduelles et contenant les valeurs de X est alors,

z = (7.5456, 5.2654, 5.2575, 5.1235, 8.1457, 8.9854, 4.1493, 4.1254

, 9.3500, 9.4578, 9.5965, 9.6160, 3.5896, 9.8308, 3.1547, 3.1386, 2.5472, 2.1475, 1.9593, 19.1245)

Nous avons donc,

E_2 = \frac{\sum_{i = 1} ^{20 - 2} (z_i - \overline{z_I}) ^2}{\sum_{i = 1} ^n (z_i - \overline{z}) ^2} = \frac{\sum_{i = 1} ^{18} (z_i - 6.168139) ^2}{\sum_{i = 1} ^n (z_i - 6.605515) ^2} = \frac{136.9162}{318.6719} = 0.4296463

Maintenant, simulons le seuil. Pour L = 100000 simulations, nous obtenons un quantile à 5\% égal à 0.4150498. Nous constatons que E_2 > 0.4150498, par conséquent nous ne pouvons conclure que la valeur minimale et maximale sont des outliers.

– Le test généralisé de la déviation extrême de Student

Nous cherchons à déterminer au sein de X si les trois observations aux valeurs les plus grandes et les plus petites peuvent être considérées comme des outliers.

Dans un premier temps déterminons les différents éléments du vecteur Z,

\overline{X} = 6.605515, sigma_X = 4.095388,

Z_1 = \frac{max_{i \in [1,20]} (| X_i -  6.605515|)}{4.095388} = \frac{12.51899}{4.095388} = 3.05685

L’observation pour laquelle max_{i \in [1,20]} (|X_i - \overline{X}|) est atteint est la valeur maximale de X, 19.1245. Elle sera donc supprimée pour le calcul de Z_2.

\overline{X_{-20}} = 5.946621, \sigma_{X_{-20}} = 5.946621,

Z_2 = \frac{max_{i \in [1,19]} (| X_i -  5.946621|)}{2.92212} = \frac{3.987321}{2.92212} = 1.36453

L’observation pour laquelle max_{i \in [1,19]} (|X_i - \overline{X_{-20}}|) est atteint est la valeur minimal de X_{-20}, 1.9593. Elle sera donc supprimé pour le calcul de Z_3.

\overline{X_{-\lbrace 1,20 \rbrace }} = 6.168139, \sigma_{X_{-\lbrace 1,20 \rbrace }} = 2.837938,

Z_3 = \frac{(max_{i \in [2,19]} (| X_i -  6.168139|)}{2.837938} = \frac{4.020639}{2.837938} = 1.416747

Maintenant, calculons la valeur seuil qui nous permettra de déterminer si nos observations peuvent être considérées comme des outliers,

\lambda_1 = \frac{(20 - 1) \times t_{(1 - \frac{0.05}{2 (20 - 1 + 1)},20-1-1)}}{\sqrt{(20 - 1 - 1 + t_{(1 - \frac{0.05}{2 (20 - 1 + 1)},20-1-1)} ^2) \times (20 - 1 + 1)}}

= \frac{19 \times t_{(0.99875,18)}}{\sqrt{(18 + t_{(0.99875,18)}) \times 20}}

= \frac{19 \times 3.510104}{\sqrt{(18 + 3.510104) \times 20}}

= \frac{66.69198}{24.62553}

= 2.708245

\lambda_2 = \frac{(20 - 2) \times t_{(1 - \frac{0.05}{2 (20 - 2 + 1)},20-2-1)}}{\sqrt{(20 - 2 - 1 + t_{(1 - \frac{0.05}{2 (20 - 2 + 1)},20-2-1)} ^2) \times (20 - 2 + 1)}}

= \frac{18 \times t_{(0.9986842,17)}}{\sqrt{(17 + t_{(0.9986842,17)}) \times 19}}

= \frac{18 \times 3.5193}{\sqrt{(17 + 3.5193) \times 19}}

= \frac{63.3474}{23.62888}

= 2.680931

\lambda_3 = \frac{(20 - 3) \times t_{(1 - \frac{0.05}{2 (20 - 3 + 1)},20-3-1)}}{\sqrt{(20 - 3 - 1 + t_{(1 - \frac{0.05}{2 (20 - 3 + 1)},20-3-1)} ^2) \times (20 - 3 + 1)}}

= \frac{17 \times t_{(0.9986111,17)}}{\sqrt{(17 + t_{(0.9986111,17)}) \times 18}}

= \frac{17 \times 3.530648}{\sqrt{(17 + 3.530648) \times 18}}

= \frac{60.02102}{22.63578}

= 2.651599

Nous constations que Z_1 > \lambda_1, Z_2 < \lambda_2, Z_3 < \lambda_3, par conséquent nous ne pouvons conclure que seul la valeur maximale de X peut être considérée comme un outlier.

– Le test \tau modifié de Thompson

Nous cherchons à déterminer les outliers au sein de de X. Nous avons \overline{X} = 6.605515. Calculons les \delta_i, \forall i \in [1,20],

\delta = (|1.9593 - 6.605515|, \cdots, |5.2654 - 6.605515|) = (4.6426215, \cdots, 1.340115)

L’Observation pour laquelle la valeur de \delta est maximale est celle qui correspond au maximum de X, soit \delta_{12} = 12.51899.

Maintenant, déterminons la valeur de référence pour un seuil de confiance de 5\%,

\tau = \frac{t_{\frac{0.05}{2},20-2} \times (20 - 1)}{\sqrt{20 \times (20 - 2 + t_{\frac{0.05}{2},20-2} ^2)}}

= \frac{t_{0.025,18} \times 19}{\sqrt{20 \times (18 + t_{0.025,18} ^2)}}

= \frac{-2.100922 \times 19}{\sqrt{20 \times (18 +4.413873)}}

= \frac{-39.91752}{\sqrt{448.2775}}

= \frac{-39.91752}{21.17256}

= -1.885342

Nous avons \sigma_X ^2 = 16.7722 et donc que,

\delta_{12} = 12.51899 > |-1.885342| \times 16.7722 = 31.61233

, nous en concluons que la valeur maximale de X est un outlier.

Lançons la seconde itération, nous supprimons donc la valeur maximale de X et obtenons,

\overline{X_{-12}} = 5.946621 et \sigma_{X_{-12}} ^2 = 8.538788

Après calcul de \delta, il en ressort que la valeur maximale correspond à la valeur minimale de X. Nous obtenons en valeur de référence et pour n = 19 observations,

\tau = -1.881106

Nous avons alors, \delta_1 = 3.9873211 < |-1.881106| \times 8.538788 = 16.06237. Nous en déduisons que la valeur minimale ne peut être considérée comme un outlier et stoppons l’algorithme à cette itération.

– Le critère de Peirce

Déterminons l’ensemble des outliers de X par le critère de Peirce. Nous débutons l’algorithme en calculons \overline{X} = 6.6055145, \sigma_X = 4.095388.

En nous reportant à la table de Peirce pour n = 20 et un outlier suspecté, nous obtenons R = 2.209. Par conséquent,

|X_i - \overline{X} |_{max} = 4.095388 \times 2.209 = 9.046712

Nous constatons que la seule observation remplissant le critère,

|X_i - \overline{X}| > 9.046712

, est la numéro douze, soit la valeur maximale de X.

Nous recherchons la présence d’un second outlier en mettant à jour la valeur de R = 1.914 depuis la table de Peirce pour n = 20 et deux outliers suspectés, nous obtenons,

|X_i - \overline{X} |_{max} = 4.095388 \times 1.914 = 7.838573

Aucune observation ne remplit le critère,

|X_i - \overline{X}| > 7.838573

Pour cette première itération nous en déduisons que X_{12} est un outlier et le supprimons du jeu de données.

Lançons une seconde itération en mettant à jour \overline{X_{-12}} = 5.946621, \sigma_{X_{-12}} = 2.92212. Nous obtenons R = 2.185 pour un outlier suspecté et n = 19, ce qui implique,

|X_i - \overline{X} |_{max} = 2.92212 \times 2.185 = 6.384832

Il se trouve qu’aucune observation ne remplit le critère,

|X_i - \overline{X_{-12}}| > 6.384832

Par conséquent, nous pouvons arrêter l’algorithme et conclure que seule la valeur maximale de X peut être considérée comme un outlier.

\bullet Application informatique:

Procédure SAS:

Package et fonction R:

\bullet Bibliographie: 

– Statistique. Dictionnaire encyclopédique de Yadolah Dodge

– Exploratory Data Analysis de John Wilder Tukey

– Simplified Statistics for Small Numbers of Observations de R. B. Dean et W. J. Dixon

– A Manual fo Spherical and Practical Astronomy de William A. Chauvenet

– Procedures for Detecting Outlying de F. E. Grubbs

– Rejecting Outliers in Factorial Designs de W. Stefansky

– Some Grubbs-Type Statistics for the Detection of Outliers de Gary Tietjen et Roger Moore

– Percentage Points for a Generalized ESD Many-Outlier Procedure de Barnard Rosner

– Outliers de John M. Cimbala

– Criterion for the Rejection of Doubtful Observations de Benjamin Peirce

L’Analyse de survie

add.png

\bullet Introduction:

L’Analyse de survie regroupe un ensemble d’outils statistique permettant de construire des indicateurs sur la durée de vie représentée par une variable aléatoire T dans le but de quantifier le délai de survenu d’un évènement d’un instant à l’autre. Une autre façon homogène de positionner le problème est de considérer l’analyse de survie comme une approche visant à étudier une variable réponse Y binaire par rapport à une matrice de P + 1 variables explicatives (T, \mathbf{X})T déjà définie auparavant et \mathbf{X} continues et/ou qualitatives.

Universellement, nous notons S la fonction de survie à modéliser et de formule,

S(t) = P(T > t)

La notion de données censurées aura également un impact important dans une analyse de survie. Une donnée censurée (ou censure) est une unité statistique qui aura été perdue de vue lors de la période d’observation mais dont l’information récoltée avant sa disparition sera prise en compte dans le modèle en étant considérée comme n’ayant ni connu ni pas connu l’évènement à mesurer.

Le principe de l’analyse de survie est principalement de modéliser une courbe de survie. Trois types d’approche existent: non paramétrique, semi-paramétrique et paramétrique.

Enfin, rappelons que le terme de survie n’est pas à prendre au sens stricte du terme puisque que l’analyse de survie s’étend également aux phénomènes d’apparition d’un évènement. La raison étant que l’analyse de survie avait été initialement développée pour mesurer le temps survenu avant un décès et a ensuite été extrapolée à d’autres thématiques de recherche.

\bullet Estimateur de Kaplan-Meier:

Présentation:

Publié en 1958 par Edward L. Kaplan et Paul Meier, l’estimateur de Kaplan-Meier, également connu sous le nom d’estimateur produit-limite, est un estimateur non paramétrique du maximum de vraisemblance visant à estimer la fonction de survie.

L’estimateur:

Pour un temps t fixé, la formule de l’estimateur de Kaplan-Meier est,

\hat{S} (t) = \prod_{t_i \leq t} \frac{n_i - Survie_i}{n_i}

Où,

n_i représente l’effectif restant pour t \geq t_i,

Survie_i = \sharp \lbrace \mbox{Evenement(s) produit(s) a } t_i = t \rbrace.

Gestion des ex-aequos:

Pour un même temps, il peut y avoir plusieurs situations différentes (apparition de l’évènement ou non). Lorsque cette situation se produit, la démarche à suivre est de considérer les cas où l’évènement a eu lieu avant ceux où il n’a pas eu lieu.

Autres estimateurs non paramétriques:

Plusieurs estimateurs existent dans cet objectif tel que celui de Nelson-Arlen, de Breslow, de Harrington-Fleming et de la méthode actuarielle. Il n’y pas vraiment de meilleur estimateur non paramétrique que les autres, disons simplement que l’estimateur de Kaplan-Meier est le plus populaire d’entre eux et donc le plus utilisé.

\bullet La régression de Cox:

Présentation:

La régression de Cox a été publiée en 1975 par David Cox dans l’objectif de modéliser des données de survies. Elle porte également le nom de modèle de Cox ou de modèle à risque proportionnel et est considérée comme une approche semi-paramétrique pour la modélisation de la courbe de survie.

Le modèle:

La régression de Cox consiste à modéliser la probabilité que l’évènement se produise au temps t alors qu’il ne s’était pas produit au temps t-1. Le modèle s’écrit,

\lambda (t, \mathbf{X}_i) = \lambda_0 (t) e ^{\sum_{p = 1} ^P \beta_p X_i ^p}

Avec \lambda_0 (t) le risque instantané de survenu de l’évènement lorsque (X_1, \cdots, X_n) sont nulles. L’Aspect semi-paramétrique de la régression de Cox provient du fait qu’il ne requiert pas l’estimation du paramètre \lambda_0(t), finalement indépendant du temps t.

L’estimation des coefficients \beta:

Posons n_E = \sharp \lbrace \mbox{Evenements produits} \rbrace. Les paramètres du modèle de Cox se déterminent de manière classique par maximisation de la vraisemblance et de formule,

L(\beta) = \prod_{i = 1} ^{n_E} \frac{e ^{\sum_{p = 1} ^P \beta_p X_i ^p}}{\sum_{i' \in [i,n_E]} e ^{\sum_{p = 1} ^P \beta_p X_i ^p}}

La résolution se fait, par exemple, via l’algorithme de Newton-Raphson,

– Etape intiale \beta = \beta_0 = (0, \cdots, 0)

– Etape 1: Calcul de la dérivée partiel de la vraisemblance,

U(\beta) = \frac{\partial}{\partial \beta} L(\beta) = (\frac{\partial}{\partial \beta_1} L(\beta), \cdots, \frac{\partial}{\partial \beta_P} L(\beta))

\Rightarrow U(\beta) = (\sum_{i = 1} ^{n_E} [X_i ^1 - \frac{\sum_{i' \in [i,n]} X_i ^1 e ^{\sum_{p = 1} ^P \beta_p X_{i'} ^p}}{\sum_{i' \in [i,n]} e ^{\sum_{p = 1} ^P \beta_p X_{i'} ^p}}], \cdots, \sum_{i = 1} ^{n_E} [X_i ^P - \frac{\sum_{i' \in [i,n]} X_i ^P e ^{\sum_{p = 1} ^P \beta_p X_{i'} ^p}}{\sum_{i' \in [i,n]} e ^{\sum_{p = 1} ^P \beta_p X_{i'} ^p}}])

– Etape 2: Calcul de la matrice jacobienne des dérivées secondes de la vraisemblance: I(\beta). Les éléments de la diagonale sont de formules,

\frac{\partial}{\partial \beta_p ^2} L(\beta) = - \frac{\sum_{i' \in [i,n]} e ^{\sum_{p = 1} ^P \beta_p X_{i'} ^p} \sum_{i' \in [i,n]} X_{i'} ^p e ^{\sum_{p = 1} ^P \beta_p X_{i'} ^p} - (\sum_{i' \in [i,n]} X_{i'} ^p e ^{\sum_{p = 1} ^P \beta_p X_{i'} ^p}) ^2}{(\sum_{i' \in [i,n]} e ^{\sum_{p = 1} ^P \beta_p X_{i'} ^p}) ^2}

, et ceux non diagonaux de formules,

\frac{\partial}{\partial \beta_{p_1} \partial \beta_{p_2}} L(\beta) = - \frac{\sum_{i' \in [i,n]} e ^{\sum_{p = 1} ^P \beta_p X_{i'} ^p} \sum_{i' \in [i,n]} X_{i'} ^{p_1} X_{i'} ^{p_2} e ^{\sum_{p = 1} ^P \beta_p X_{i'} ^p} - \sum_{i' \in [i,n]} X_{i'} ^{p_1} e ^{\sum_{p = 1} ^P \beta_p X_{i'} ^p} \sum_{i' \in [i,n]} X_{i'} ^{p_2} e ^{\sum_{p = 1} ^P \beta_p X_{i'} ^p}}{(\sum_{i' \in [i,n]} e ^{\sum_{p = 1} ^P \beta_p X_{i'} ^p}) ^2}

– Etape 3: a l’itération k+1,

\beta ^{k+1} = \beta - U(\beta) \cdot I(\beta) ^{-1}

– Etape finale: Si || \beta ^{k+1} - \beta ^k || > \epsilon alors nous reprenons les étapes 1, 2, 3 avec \beta = \beta ^{k+1}. Sinon, l’algorithme s’arrête et nous avons notre vecteur de coefficients finaux.

Les hypothèses de validité:

Le modèle de Cox repose sur l’hypothèse forte des risques proportionnels et nécessite qu’elle soit vérifiée avant application. Pour cela il s’agit de vérifier que le temps n’a aucune influence bénéfique, nocif ou nul sur les différentes variables de \mathbf{X}.

Le risque relatif ajusté ou hazard ratio (RR):

Du modèle de Cox, il est possible de déduire le risque relatif associé à chaque variable. Deux formats sont à considérer,

– si Y est binaire ou continue alors,

RR_p = e ^{\beta_p}, p \in [1,P]

– si Y est ordinale, à valeurs dans [a, b], alors,

RR_p = e ^{\beta_p (b - a)}, p \in [1,P]

Tests associés au modèle de Cox:

Trois principaux tests existent pour juger de la fiabilité des coefficients,

– Le test du rapport de vraisemblance,

\chi_V ^2 = 2 [log L(\beta) - log(\beta_0)]

– Le test de Wald,

\chi_W ^2 = (\beta - \beta_0) I(\beta) (\beta - \beta_0)

– Le test du Score,

\chi_S ^2 = (U(\beta)) ^T (I(\beta_0)) ^{-1} U(\beta_0)

Ces trois tests suivent une loi du \chi ^2 à P degrés de liberté et l’hypothèse H_0 est \beta = \beta_0 soit que le modèle considéré avec l’information auxiliaire \mathbf{X} est aussi performant que sans.

Gestion des ex-aequos:

Dans le cas de survenue de plusieurs évènements pour un même temps, le calcul des dérivées partielles de la vraisemblance doit être modifié. Plusieurs méthodes existent, parmi les plus populaires nous avons la méthode de Breslow et celle d’Efron.

\bullet Approches paramétriques:

Une approche purement paramétrique est également possible. Il s’agit de partir du constat que la courbe de survie est identifiable à l’une des nombreuses loi de probabilité dont les principales sont,

– La loi Exponentielle,

S(t | \theta) = e ^{-\frac{t}{\theta}}

– La loi de Weibull,

S(t | \theta, \alpha) = e ^{-(\frac{t}{\theta}) ^{\alpha}}

– La loi Gamma,

S(t | \theta, \alpha) = \frac{\frac{\alpha^{\theta}}{\Gamma(\theta)} t ^{\theta - 1} e ^{\alpha t}}{1 - \frac{1}{\Gamma(\theta)} \int_0 ^{\alpha t} u ^{\theta - 1} e ^{-u} du}

\bullet Le test du log-rank:

Présentation:

Publié par Nathan Mantel et David Cox en 1966, le test du log-rank (connu aussi sous le nom de test de Mantel-Cox) est une approche non paramétrique permettant de comparer les distributions de deux séries de survie.

Le test:

Le calcul de la statistique de test peut se résumer en trois étapes itératives. Soit T l’ensemble des temps considéré pour les deux séries considérées.

– Sortie des temps pour lesquelles l’évènement que nous cherchons à mesurer n’est pas survenu au sein des deux séries d’intérêt,

T_S = T_S (serie_1) \cup T_S (serie_2)

– Pour chacun des temps retenus dans T_S = (T_S ^1, \cdots, T_S ^t, \cdots), construire le tableau,

add2

– Calcul des effectifs observés (O_t), théoriques (E_t) et de la variance (V_t) associé au tableau ci-dessus,

O_t = C

E_t = (A + B) \frac{C + D}{A + B + C + D}

V_t = \frac{(A + B) (C + D) (A + C) [(A + B + C + D) - (A + C)]}{(A + B + C + D) ^2 [(A + B + C + D) - 1]}

A noter que si V_t \rightarrow \infty alors .</span></p> <p style="text-align:justify;"><span style="color:#808080;">– Une fois les différents temps parcourus, calculer la statistique de test du log-rank,</span></p> <p style="text-align:center;"><span style="color:#808080;">[latex]Z = |\frac{(\sum_t O_t - \sum_t E_t) ^2}{\sum_t V_t}|

Elle suit une loi du \chi ^2 à un degré de liberté et l’hypothèse H_0 est: « les deux courbes de survie sont différentes ».

Autres tests pour la comparaison de deux courbes de survie:

D’autres tests existent dans cet objectif, parmi les plus populaires nous avons le test du log-rank pondéré, le test de Gehan et le test de Peto et Prentice.

Tous partagent la même condition que le test du log-rank, les courbes de survies ne doivent pas se croiser sous peine d’occasionner une forte chute de la puissance statistique.

\bullet Exemple:

Soit l’échantillon suivant,

add

, où la variable « Temps » est le moment où nous observons la production de l’évènement mesuré (« Survie = 0« ) ou non (« Survie = 1« ).

Estimateur de Kaplan-Meier:

Pour le groupe A,

– Nous avons pour « Temps = 59« ,

\hat{S} (t = 59) = \frac{6 - 1}{6} = 0.83333333333

– Puis pour « Temps \in [59, 115]« ,

\hat{S} (t \in [59, 115]) = \frac{5- 1}{5} \times \hat{S} (t = 59) = 0.8 \times 0.83333333333 = 0.666666667

– Puis pour « Temps \in [59, 115, 156]« ,

\hat{S} (t \in [59, 115, 156]) = \frac{4 - 1}{4} \times \hat{S} (t \in [59,115]) = 0.75 \times 0.666666667 = 0.5

– Puis pour « Temps \in [59, 115, 156, 431]« ,

\hat{S} (t \in [59, 115, 156, 431]) = \frac{3 - 1}{3} \times \hat{S} (t \in [59,115,156]) = 0.333333333

– Puis pour « Temps \in [59, 115, 156, 431, 448]« ,

\hat{S} (t \in [59, 115, 156, 431, 448]) = \frac{2 - 0}{2} \times \hat{S} (t \in [59, 115, 156, 431]) = 0.333333333

– Puis pour « Temps \in [59, 115, 156, 431, 448, 477]« ,

\hat{S} (t \in [59, 115, 156, 431, 448, 477]) = \frac{1 - 0}{1} \times \hat{S} (t \in [59, 115, 156, 431, 448])

= 0.333333333

En procédant aux mêmes calculs pour le groupe B, nous pouvons alors les deux courbes de Kaplan-Meier suivantes,

add

Test du Log-Rank:

La courbe ci-dessus semble laisser penser qu’il y a de grandes différences entre le groupe A et le groupe B. Vérifions le au sens statistique du terme.

Dans un premier temps, relevons les temps pour lesquelles l’évènement n’a pas lieu (« survie = 1« ). Nous avons pour les groupes A et B:

(T_S)_A = (59,115,156,431) et (T_S)_B = (464,475,563)

Quatre cas pour le groupe A et trois cas pour le groupe B. Nous obtenons alors les temps fusionnés suivants,

T_S = (59,115,156,431,464,475,563)

Procédons itérativement pour l’ensemble des temps considérés (évènement produit ou non),

– pour T = 59,

add2

\Rightarrow O_{59} = 0

\Rightarrow E_{59} = (1 + 0) \times \frac{0 + 4}{1 + 5 + 0 + 4} = 1 \times \frac{4}{10} = 0.4

\Rightarrow V_{59} = \frac{(1 + 5) \times (0 + 4) \times (1 + 0) \times ((1 + 5 + 0 + 4) - (1 + 0))}{(1 + 5 + 0 + 4) ^2 \times ((1 + 5 + 0 + 4) - 1} = \frac{6 \times 4 \times 1 \times 9}{10 ^2 \times 9} = \frac{216}{900} = 0.24

– pour T = 115,

add2

\Rightarrow O_{115} = 0

\Rightarrow E_{115} = (1 + 0) \times \frac{0 + 4}{1 + 4 + 0 + 4} = 1 \times \frac{4}{9} = 0.444444444

\Rightarrow V_{115} = \frac{(1 + 4) \times (0 + 4) \times (1 + 0) \times ((1 + 4 + 0 + 4) - (1 + 0))}{(1 + 4 + 0 + 4) ^2 \times ((1 + 4 + 0 + 4) - 1} = \frac{5 \times 4 \times 1 \times 8}{9 ^2 \times 8} = \frac{160}{648} = 0.2469136

\cdots

– pour T = 563,

add2

\Rightarrow O_{563} = 1

\Rightarrow E_{563} = (0 + 1) \times \frac{1 + 0}{0 + 0 + 1 + 0} = 1 \times \frac{1}{1} = 1

\Rightarrow V_{563} = \frac{(0 + 0) \times (1 + 0) \times (0 + 1) \times ((0 + 0 + 1 + 0) - (0 + 1))}{(0 + 0 + 1 + 0) ^2 \times ((0 + 0 + 1 + 0) - 1} = \frac{0 \times 1 \times 1 \times 0}{1 ^2 \times 0} = 0

Nous pouvons alors calculer la statistique du test du Log-Rank,

O = \sum_{t \in [59,115,156,431,464,475,563]} O_t =0 + 0 + 0 + 0 + 1 + 1 + 1 = 3

E = \sum_{t \in [59,115,156,431,464,475,563]} E_t

= 0.4 + 0.444444444 + 0.5 + 0.5 + 0.75 + 0.6666667 + 1

= 4.261111

V = \sum_{t \in [59,115,156,431,464,475,563]} V_t

= 0.24 + 0.2469136 + 0.25 + 0.25 + 0.1875 + 0.2222222 + 0

= 1.138737

Par conséquent,

|Z| = |\frac{(3 - 4.261111) ^2}{1.396636}| = 1.138737

, que nous reportons à la loi du \chi ^2 pour un degré de liberté et nous obtenons p = 0.2926337. Nous ne pouvons rejeter H_0 et donc nous en concluons que les deux courbes ne sont pas identiques.

Régression de Cox:

Passons maintenant à la modélisation de l’évènement « Survie » pour les groupes A et B fusionnés. Nous initialisons l’algorithme de Newton-Raphson à \beta ^0 = 0.

Nous avons pour la première itération,

U(\beta ^0) = \sum_{i = 1} ^{10} Survie_i - \sum_{i = 1} ^{10} \frac{\sum_{i' \in [i,10]} Survie_{i'} e ^{0 \times Survie_{i'}}}{\sum_{i' \in [i,10]} e ^{0 \times Survie_{i'}}}

= 7 - (\frac{7}{10} + \frac{6}{9} + \cdots + \frac{1}{1})

= 7 - 6.746429

= 0.253571

I(\beta ^0) = - \sum_{i = 1} ^{10} \frac{\sum_{i' \in [i,10]} e ^{0 \times Survie_{i'}}\times \sum_{i' \in [i,10]} Survie_{i'} ^2 e ^{0 \times Survie_{i'}} - (\sum_{i' \in [i,10]} Survie_{i'} e ^{0 \times Survie_{i'}}) ^2}{(\sum_{i' \in [1,10]} e ^{0 \times Survie_{i'}}) ^2}

= - \frac{70 - 49}{100} + \frac{54 - 36}{81} + \cdots + \frac{1 - 1}{1}

= - 2.03344

Nous avons alors,

\beta ^1 = 0 - \frac{0.253571}{- 2.03344} = 0.1247005

Pour la seconde itération et en considérant \beta ^1 à la place de \beta ^0, nous obtenons,

U(\beta ^1) = 0.004441066

I(\beta ^1) = -1.960438

\beta ^2 = 0.1269661

Pour la troisième itération et en considérant \beta ^2 à la place de \beta ^1, nous obtenons,

U(\beta ^2) = 0.000001609696

I(\beta ^2) = -1.959017

\beta ^3 = 0.1269669

Pour la quatrième itération et en considérant \beta ^3 à la place de \beta ^2, nous obtenons,

U(\beta ^2) = 0.0000000000002113865

I(\beta ^2) = -1.959016

\beta ^4 = 0.1269669

Nous avons \beta ^4 = \beta ^3, nous en déduisons que l’agorithme de Newton-Raphson a convergé. Le coefficient associé à la variable « Survie » est donc \beta = 0.1269669.

\bullet Application informatique:

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

– Package R: https://stat.ethz.ch/R-manual/R-devel/library/survival/html/survfit.html

\bullet Bibliographie: 

– Comprendre et utiliser les statistiques dans les sciences de la vie de Bruno Falissard

– Nonparametric estimation from incomplete observations de Edward L. Kaplan et Paul Meier

– Evaluation of survival data and two new rank order de Nathan Mantel

– Asymptotically Efficient Rank Invariant Test Procedures de Richard Peto

– Linear Rank Tests in Survival Analysis de David Harrington

– Introduction à l’analyse de survie de Philippe Saint Pierre

– Partial Likelihood de David R. Cox

Le test de Jonckheere

add

\bullet Présentation:

Publié en 1954 par Aimable Robert Jonckheere, le test de Jonckheere est une approche non paramétrique permettant de tester si les variables continues appariées X ^1, \cdots, X ^T sont appariées.

\bullet Le test:

Hypothèse préliminaire: Variables appariées continues pour T \geq 2.

La formule de la statistique de test de Jonckheere est,

J = \sum_{i = 1} ^n U_i

Où,

U_i = \sum_{t_1 = 1} ^{T-1} \sum_{t_2 = t_1 + 1} ^T I_{(X_i ^{t_1} < X_i ^{t_2})}

, est la statistique de Mann-Whitney adaptée aux échantillons appariés avec I(.) indicatrice valant notamment \frac{1}{2} si nous sommes en présence de deux termes ex-aequos.

A noter que le lien linéaire entre le U de Mann-Whitney et le test de Wilcoxon implique que la statistique de test de Jonckheere-Terpstra peut également se réécrire en fonction du W de Wilcoxon.

Si n est suffisamment grand, alors nous travaillons sur la statistique de test,

Z = | \frac{J - E[J]}{\sqrt{V(J)}} |

Avec,

E[J] = \frac{1}{4} n T (T - 1)

V(J) = \frac{n T (T - 1) (2 T + 5)}{72}

Ci-dessous la table de la loi normale centrée-réduite.addTendance pour le rejet de H_0:

Plus la statistique de test Z est grande et plus nous avons de chance de rejeter H_0. Ce qui revient à dire que,

Z \rightarrow \infty \Rightarrow J \rightarrow \infty

J est la somme des U_i de Mann-Withney qui explose si pour chaque observation nous pouvons constater que les valeurs sont toutes inférieures ou supérieures d’un temps successif à l’autre. Soit qu’au fur et à mesure que T croît, les valeurs de X augmentent ou diminuent marquant là un lien entre les deux informations.

\bullet Tendance lorsque n \longrightarrow \infty:

– Nous avons fait tourner cinq simulations afin d’étudier la robustesse du test de Jonckheere. Nous générons des échantillons de taille 10, 10^2, 10^3, 10 ^4 puis 10 ^5 de telle manière à ce que les différents sous-échantillons soient bien distincts et nous étudions si le fait d’augmenter le paramètre n a une influence sur le fait d’accepter H_0 à tort.

add

Nous constatons que nous rejetons à chaque fois l’hypothèse H_0, ce qui est en adéquation avec les hypothèses fixées.

– Nous avons fait tourner cinq simulations afin d’étudier la robustesse du test de Jonckheere. Nous générons des échantillons de taille 10, 10^2, 10^3, 10 ^4 puis 10 ^5 de telle manière à ce que les différents sous-échantillons soient légèrement confondus et nous étudions si le fait d’augmenter le paramètre n a une influence sur le fait de rejeter H_0 à tort.

add2

Nous constatons que passé un échantillon de taille N = 10 000 nous rejetons à tort H_0.

Nos simulations montrent bien que le test de Jonckheere est influencé par la taille d’échantillon.

\bullet Annexe théorique:

Cette partie présente la démonstration liant le U de Mann-Whitney et le test de Wilcoxon (W).

La statistique de Mann-Whitney est: U = \sum_{i_1 = 1} ^{n_1} \sum_{i_2 = 1} ^{n_2} 1_{(X|_{g1})_{i_1} > (X|_{g2})_{i_2}} et plus particulièrement:

– U_{X|_{g1}} = \sum_{i_1 = 1} ^{n_1} \sum_{i_2 = 1} ^{n_2} (1_{(X|_{g1})_{i_1} > (X|_{g2})_{i_2}} + \frac{1}{2} \times 1_{((X|_{g1})_{i_1} = (X|_{g2})_{i_1})})

U_{X|_{g2}} = \sum_{i_1 = 1} ^{n_1} \sum_{i_2 = 1} ^{n_2} (1_{((X|_{g1})_{i_1} < (X|_{g2})_{i_2})} + \frac{1}{2} \times 1_{(X|_{g1})_{i_1} = (X|_{g2})_{i_1}})

Il faut voir le calcul de la somme des rangs de manière inverse, c’est-à-dire qu’il faut comprendre que \forall l, (R|_{g1})_l est égal au \sharp \lbrace X|_{g2} < X|_{g1} \mbox{ a partir de l'indice } l \rbrace , remarquons que nous appliquons un décalage de l’indice l à chaque calcul de  (R|_{g1})_l. Dés lors nous avons que:

\sum_{l = 1} ^{n_1} (R|_{g1})_l = \sum_{i_1 = 1} ^{n_1} \sum_{i_2 = 1} ^{n_2} (1_{(X|_{g1})_{i_1} < (X|_{g2})_{i_2}} + \frac{1}{2} \times 1_{(X|_{g1})_{i_1} = (X|_{g2})_{i_2}}) + \sum_{l = 1} ^n l
= U_{X|_{g2}} + \frac{n \cdot (n + 1)}{2}

Or, U_{X|_{g1}} + U_{X|_{g2}} = \sum_{i_1 = 1} ^{n_1} \sum_{i_2 = 1} ^{n_2} (1_{(X|_{g1})_{i_1} > (X|_{g2})_{i_2}} + \frac{1}{2} \times 1_{(X|_{g1})_{i_1} = (X|_{g2})_{i_1}}) + \sum_{i_1 = 1} ^{n_1} \sum_{i_2 = 1} ^{n_2} (1_{(X|_{g1})_{i_1} < (X|_{g2})_{i_2}} + \frac{1}{2} \times 1_{(X|_{g1})_{i_1} = (X|_{g2})_{i_1}})

= \sum_{i_1 = 1} ^{n_1} \sum_{i_2 = 1} ^{n_2} (1_{(X|_{g1})_{i_1} > (X|_{g2})_{i_2}} + 1_{X|_{g1})_{i_1} > (X|_{g2})_{i_2}} + 1_{(X|_{g1})_{i_1} = X|_{g2})_{i_2}}) = n_1 \times n_2

Puisque 1_{(X|_{g1})_{i_1} > (X|_{g2})_{i_2}} + 1_{(X|_{g1})_{i_1} > (X|_{g2})_{i_2}} + 1_{(X|_{g1})_{i_1} = X|_{g2})_{i_2}} = 1.

Par conséquent, \sum_{l = 1} ^{n_1} (R|_{g1})_l = U_{X|_{g2}} + \frac{n \cdot (n + 1)}{2} = n_1 \cdot n_2 - U_{X|_{g2}} + \frac{n \cdot (n + 1)}{2} et comme \sum_{l = 1} ^{n_1} (R|_{g1})_l est la somme des rangs provenant du test de Wilcoxon, le lien entre les deux statistiques est fait.

\bullet Exemple:

Soit les cinq échantillons appariés X ^1, X ^2, X ^3, X ^4, X ^5 suivants.

add

Ci-dessous, les boxplots associés à (X ^1, X ^2, X ^3, X ^4, X ^5).

add

Les boxplots montrent que X varie énormément en fonction du temps t \in \lbrace 1, 2, 3, 4, 5 \rbrace. Prouvons-le statistiquement.

Tout d’abord, commençons par déterminer la matrice du nombre de valeurs supérieurs à chaque temps t \in [1,5], nous obtenons,

U = \begin{pmatrix} 4 & 1 & 1 & 1 & 0 \\ 4 & 1 & 1 & 1 & 0 \\ 3 & 3 & 1 & 1 & 0 \\ 4 & 1 & 1 & 1 & 0 \\ 4 & 1 & 1 & 1 & 0 \\ 1 & 3 & 1 & 0 & 0 \\ 1 & 2 & 1 & 0 & 0 \\ 1 & 1 & 2 & 0 & 0 \\ 2 & 1 & 1 & 0 & 0 \\ 1 & 1 & 1 & 0 & 0 \\ 1 & 2 & 2 & 0 & 0 \\ 1 & 1 & 2 & 0 & 0 \\ 1 & 2 & 2 & 0 & 0 \\ 1 & 3 & 2 & 0 & 0 \\ 1 & 3 & 2 & 0 & 0 \\ 1 & 3 & 2 & 1 & 0 \\ 1 & 2 & 2 & 1 & 0 \\ 1 & 2 & 2 & 1 & 0 \\ 1 & 3 & 2 & 1 & 0 \\ 1 & 3 & 1 & 1 & 0 \\ \end{pmatrix}

Maintenant, nous pouvons calculer la statistique de test,

J = \sum_{i = 1} ^n \sum_{t = 1} ^T U_{i,t} = 4 + 1 + 1 + 1 + 0 + \cdots + 1 + 3 + 1 + 1 + 0 = 114

Reste à déterminer l’espérance et la variance de J. Nous avons,

E[J] = \frac{20 \times 5 \times (5 - 1)}{4} = \frac{400}{4} = 100

Et,

V(J) = \frac{20 \times 5 \times (5 - 1) \times (2 \times 5 + 5)}{72} = \frac{6000}{72} = 83.33333

Nous avons donc,

Z = | \frac{114 - 100}{\sqrt{83.33333}} | = | \frac{14}{9.128709} | = | 1.533623 | = 1.533623

Si nous reportons cette valeur de la statistique de test à la table de la loi normale centrée-réduite, nous obtenons une p-valeur de 0.06256121 > 5 \%. Nous en concluons que nous ne pouvons rejeter H_0 et qu’il n’y a d’influence du temps T sur X au sens statistique du terme.

\bullet Application informatique:

Procédure SAS:

%macro Jonckheere_test (TAB=);

/* La macro exécute le test de Jonckheere pour données appariées sur la table de données TAB contenant les variables continues */

/* Options pour nettoyer la log et les sorties */
options nonotes spool;
ods noresults;

/* Récupération du nombre d’observations n */
proc sql noprint;
select count(*) into: n from &TAB.;
quit;

proc contents data = &TAB. out = vars;
run;

/* Récupération du nombre de temps T et de la liste des variables appariées */
proc sql noprint;
select distinct(name) into: list_T separated by  »  » from vars;
select count(*) into: T from vars;
quit;

/* Initialisation de la statistique de test de Jonckheere */
%let J = 0;

%let Tmu = %eval(&T. – 1);

%do i = 1 %to &n.;

%do t1 = 1 %to &Tmu.;

%let tt1 = %eval(&t1. + 1);
%let c = 0;

/* Récupération de la valeur de référence à comparer à toutes les autres */
data _null_;
set &TAB.;
if _n_ = &i.;
call symputx(« X_ref »,%scan(&list_T.,&t1.));
run;

%do t2 = &tt1. %to &T.;

/* Récupération des valeurs à comparer à la valeur de référence */
data _null_;
set &TAB.;
if _n_ = &i.;
call symputx(« X_comp »,%scan(&list_T.,&t2.));
run;

/* Si inférieur alors on incrémente de 1 selon la définition de la fonction indicatrice */
%if &X_ref. < &X_comp. %then %do;
%let c = %sysevalf(&c. + 1);
%end;

/* Si égal alors on incrémente de 0.5 selon la définition de la fonction indicatrice */
%if &X_ref. = &X_comp. %then %do;
%let c = %sysevalf(&c. + 0.5);
%end;

%end;

/* Mise à jour de la statistique de test de Jonckheere */
%let J = %sysevalf(&J. + &c.);

%end;

%end;

/* Création de la table de résultats */
data Résultats;
Effectif = &n.;
Temps = &T.;
J = &J.;
Esperance = &n.*&T.*(&T.-1)/4;
Variance = &n.*&T.*(&T.-1)*(2*&T.+5)/72;
Z = abs((J-Esperance)/sqrt(Variance));
p_valeur = 1-probnorm(Z);
run;

/* Suppression des tables inutiles */
proc datasets lib = work nolist;
delete vars;
run;

/* Reset options */
options notes nospool;
ods results;

%mend;

Package et fonction R:

Jonckheere.test = function(TAB) {

# La fonction procède au test de Jonckheere pour données appariées sur la matrice TAB qui doit être composée de variables continues

# Suppression des données manquantes et récupération du nombre d’observations n et du nombre de temps T
TAB = na.omit(TAB)
n = dim(TAB)[1]
T = dim(TAB)[2]

# Création de la matrice synthétisant le nombre de valeur supérieur ou égale
TT = matrix(0,n,T)
for (i in 1:n) {
for(t1 in 1:(T-1)) {
c = 0
for (t2 in (t1+1):T) {
# Si la valeur est supérieur alors on incrémente de 1 selon la définition de l’indicatrice
if (TAB[i,t1]<TAB[i,t2]) {c = c + 1}
# Si la valeur est égale alors on incrémente de 0.5 selon la définition de l’indicatrice
if (TAB[i,t1]==TAB[i,t2]) {c = c + 0.5}
}
TT[i,t1] = c
}
}

# Calcul de la statistique de test de Jonckheere
J = sum(TT)
# Calcul de l’espérance
E = n*T*(T-1)/4
# Calcul de la variance
V = n*T*(T-1)*(2*T+5)/72

# Calcul de la statistique de test centrée-réduite
Z = abs((J-E)/sqrt(V))
# Calcul de la p-valeur selon la loi normale centrée-réduite
p = pnorm(Z,0,1,lower.tail=FALSE)

# Mise en forme et retour des résultats
return(data.frame(Pop = n, Temps = T, J = J, Esp = E, Var = V, Z = Z, p_valeur = p))

}

\bullet Bibliographie:

– Comparaison de populations. Tests non paramétriques de Ricco Rakotomalala

Le test de Fligner-Policello

add

\bullet Présentation:

Publié en 1981 par Michael A. Fligner et George E. Policello II, le test de Fligner-Policello, également appelé test des rangs robustes de Fligner-Policello, est une approche non paramétrique permettant de tester si X|_{Y = 1}, X|_{Y = 2}, les sous-échantillons d’une variable continue X restreinte aux deux groupes de Y, ont même médiane.

A noter que la manière dont a été pensé le test impose que passer une certaine taille d’échantillon, le test de Fligner-Policello est particulièrement consommateur en temps de calcul.

\bullet Le test:

Hypothèse préliminaire: X continue et Y binaire.

La formule de la statistique de test de Fligner-Policello est,

Z = \frac{\sum_{i_2 = 1} ^{n_2} (P_2)_{i_2} - \sum_{i_1 = 1} ^{n_1} (P_1)_{i_1}}{2 \sqrt{V_1 + V_2 + \overline{P_1} \overline{P_2}}}

Où,

\forall i_1 \in [1,n_1], (P_1)_{i_1} = \sum_{i_2 = 1} ^{n_2} [I_{(X|_{Y = 2})_{i_2} < (X|_{Y = 1})_{i_1}} + \frac{1}{2} I_{(X|_{Y = 2})_{i_2} = (X|_{Y = 1})_{i_1}}]

\forall i_2 \in [1,n_2], (P_2)_{i_2} = \sum_{i_1 = 1} ^{n_1} [I_{(X|_{Y = 1})_{i_1} < (X|_{Y = 2})_{i_2}} + \frac{1}{2} I_{(X|_{Y = 1})_{i_1} = (X|_{Y = 2})_{i_2}}]

Il s’agit en réalité des statistiques de test U_{X|_{g_1}}, U_{X|_{g_2}} de Mann-Whitney. Le lien linéaire existant entre elle et le W de Wilcoxon implique que le test de Fligner-Policello peut se réécrire en fonction de cette dernière. Enfin,

\forall k \in [1,2], \overline{P_k} = \frac{\sum_{i_k = 1} ^{n_k} (P_k)_{i_k}}{n_k}

\forall k \in [1,2], V_k = \sum_{i_k = 1} ^{n_k} [(P_k)_{i_k} - \overline{P_k}] ^2

La statistique de test suit une loi normale centrée-réduite et l’hypothèse H_0 est: « Les médianes sont égales / \theta_1 = \theta_2« .

Ci-dessous la table de la loi normale centrée-réduite.

add

Tendance pour le rejet de H_0:

Plus la statistique de test Z est grande et plus nous avons de chance de rejeter H_0. Ce qui implique,

\sum_{i_2 = 1} ^{n_2} (P_2)_{i_2} - \sum_{i_1 = 1} ^{n_1} (P_1)_{i_1} \rightarrow \infty

\Rightarrow \sum_{i_2 = 1} ^{n_2} (P_2)_{i_2} >>>> \sum_{i_1 = 1} ^{n_1} (P_1)_{i_1}

\Rightarrow \sum_{i_2 = 1} ^{n_2} \sum_{i_1 = 1} ^{n_1} [I_{(X|_{Y = 1})_{i_1} < (X|_{Y = 2})_{i_2}} + \frac{1}{2} I_{(X|_{Y = 1})_{i_1} = (X|_{Y = 2})_{i_2}}] >>>> \sum_{i_1 = 1} ^{n_1} \sum_{i_2 = 1} ^{n_2} [I_{(X|_{Y = 2})_{i_2} < (X|_{Y = 1})_{i_1}} + \frac{1}{2} I_{(X|_{Y = 1})_{i_1} = (X|_{Y = 2})_{i_2}}]

\Rightarrow \forall i_1, i_2, \sum_{i_2 = 1} ^{n_2} \sum_{i_1 = 1} ^{n_1} I_{(X|_{Y = 1})_{i_1} < (X|_{Y = 2})_{i_2}} >>>> \sum_{i_1 = 1} ^{n_1} \sum_{i_1 = 2} ^{n_2} I_{(X|_{Y = 2})_{i_2} < (X|_{Y = 1})_{i_1}}

\Rightarrow \forall i_1, i_2, (X|_{Y = 2})_{i_2} > (X|_{Y = 1})_{i_1}

Soit que les valeurs de X sont ordonnées en fonction des groupes de Y et donc que les médianes sont différentes.

\bullet Tendance lorsque n \longrightarrow \infty:

Nous proposons ici de vérifier si le test de Fligner-Policello est sensible aux grands échantillons ou non.

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 médianes sont différentes d’un groupe à l’autre.

add2

Globalement, quelque soit la taille de l’échantillon, le test statistique rejette H_0, ce qui est en accords avec nos hypothèses.

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

add2

Jusqu’à N = 100 nous restons cohérent avec nos hypothèses, cependant nous voyons que pour un échantillon de taille N = 1000, le test rejette H_0 à tort.

Nous en déduisons que le test de Fligner-Policello est influencé par la taille de l’échantillon.

\bullet Annexe théorique:

Nous présentons ici une esquisse de la démonstration du lien entre le U de Mann-Whitney et le W de Wilcoxon.

La statistique de Mann-Whitney est: U = \sum_{i_1 = 1} ^{n_1} \sum_{i_2 = 1} ^{n_2} 1_{(X|_{g1})_{i_1} > (X|_{g2})_{i_2}} et plus particulièrement:

– U_{X|_{g1}} = \sum_{i_1 = 1} ^{n_1} \sum_{i_2 = 1} ^{n_2} (1_{(X|_{g1})_{i_1} > (X|_{g2})_{i_2}} + \frac{1}{2} \times 1_{((X|_{g1})_{i_1} = (X|_{g2})_{i_1})})

U_{X|_{g2}} = \sum_{i_1 = 1} ^{n_1} \sum_{i_2 = 1} ^{n_2} (1_{((X|_{g1})_{i_1} < (X|_{g2})_{i_2})} + \frac{1}{2} \times 1_{(X|_{g1})_{i_1} = (X|_{g2})_{i_1}})

Il faut voir le calcul de la somme des rangs de manière inverse, c’est-à-dire qu’il faut comprendre que \forall l, (R|_{g1})_l est égal au \sharp \lbrace X|_{g2} < X|_{g1} \mbox{ a partir de l'indice } l \rbrace , remarquons que nous appliquons un décalage de l’indice l à chaque calcul de  (R|_{g1})_l. Dés lors nous avons que:

\sum_{l = 1} ^{n_1} (R|_{g1})_l = \sum_{i_1 = 1} ^{n_1} \sum_{i_2 = 1} ^{n_2} (1_{(X|_{g1})_{i_1} < (X|_{g2})_{i_2}} + \frac{1}{2} \times 1_{(X|_{g1})_{i_1} = (X|_{g2})_{i_2}}) + \sum_{l = 1} ^n l
= U_{X|_{g2}} + \frac{n \cdot (n + 1)}{2}

Or, U_{X|_{g1}} + U_{X|_{g2}} = \sum_{i_1 = 1} ^{n_1} \sum_{i_2 = 1} ^{n_2} (1_{(X|_{g1})_{i_1} > (X|_{g2})_{i_2}} + \frac{1}{2} \times 1_{(X|_{g1})_{i_1} = (X|_{g2})_{i_1}}) + \sum_{i_1 = 1} ^{n_1} \sum_{i_2 = 1} ^{n_2} (1_{(X|_{g1})_{i_1} < (X|_{g2})_{i_2}} + \frac{1}{2} \times 1_{(X|_{g1})_{i_1} = (X|_{g2})_{i_1}})

= \sum_{i_1 = 1} ^{n_1} \sum_{i_2 = 1} ^{n_2} (1_{(X|_{g1})_{i_1} > (X|_{g2})_{i_2}} + 1_{X|_{g1})_{i_1} > (X|_{g2})_{i_2}} + 1_{(X|_{g1})_{i_1} = X|_{g2})_{i_2}}) = n_1 \times n_2

Puisque 1_{(X|_{g1})_{i_1} > (X|_{g2})_{i_2}} + 1_{(X|_{g1})_{i_1} > (X|_{g2})_{i_2}} + 1_{(X|_{g1})_{i_1} = X|_{g2})_{i_2}} = 1.

Par conséquent, \sum_{l = 1} ^{n_1} (R|_{g1})_l = U_{X|_{g2}} + \frac{n \cdot (n + 1)}{2} = n_1 \cdot n_2 - U_{X|_{g2}} + \frac{n \cdot (n + 1)}{2} et comme \sum_{l = 1} ^{n_1} (R|_{g1})_l est la somme des rangs provenant du test de Wilcoxon, le lien entre les deux statistiques est fait.

\bullet Exemple:

Soit la variable aléatoire X distribuée selon deux groupes d’une variable Y:

add

Ci-dessous, les densités associées aux distributions de X selon nos deux groupes. Nous pourrons remarquer que nos deux sous-échantillons ont la même distribution et donc que les médianes semblent égales.

addDans un premier temps, déterminons les vecteurs,

P_1 = P((X|_{Y = 1}))

= (\sum_{i_2 = 1} ^{10} [I_{(X|_{Y = 2})_{i_2} < 8.1472} + \frac{1}{2} \times I_{(X|_{Y = 2})_{i_2} = 8.1472}], \cdots, \sum_{i_2 = 1} ^{10} [I_{(X|_{Y = 2})_{i_2} < 9.6489} + \frac{1}{2} \times I_{(X|_{Y = 2})_{i_2} = 9.6489}])

= (6,6,0,6,4,0,2,4,8,9)

\Rightarrow \sum_{i_1 = 1} ^{10} (P_1)_{i_1} = 6 + \cdots + 9 = 45

P_2 = P((X|_{Y = 2}))

= (\sum_{i_1 = 1} ^{10} [I_{(X|_{Y = 1})_{i_1} < 1.5761} + \frac{1}{2} \times I_{(X|_{Y = 1})_{i_1} = 1.5761}], \cdots, \sum_{i_1 = 1} ^{10} [I_{(X|_{Y = 1})_{i_1} < 9.5949} + \frac{1}{2} \times I_{(X|_{Y = 1})_{i_1} = 9.5949}])

= (2,10,8,3,5,2,3,8,5,9)

\Rightarrow \sum_{i_2 = 1} ^{10} (P_2)_{i_2} = 2 + \cdots + 9 = 55

A partir de ces deux objets, nous pouvons déterminer les éléments manquants,

\overline{P_1} = \frac{45}{10} = 4.5

\overline{P_2} = \frac{55}{10} = 5.5

V_1 = V_{X|_{Y = 1}}

= \sum_{i_1 = 1} ^{10} [(P_1)_{i_1} - 4.5] ^2

= (6 - 4.5) ^2 + \cdots + (9 - 4.5) ^2

= 86.5

V_2 = V_{X|_{Y = 2}}

= \sum_{i_2 = 1} ^{10} [(P_2)_{i_2} - 5.5] ^2

= (2 - 5.5) ^2 + \cdots + (9 - 5.5) ^2

= 82.5

Il ne nous reste plus qu’à calculer la statistique de test,

Z = \frac{55 - 45}{2 \times \sqrt{86.5+82.5 + 4.5 \times 5.5}} = \frac{10}{2 \times \sqrt{193.75}} = \frac{10}{27.83882} = 0.3592106

Si nous reportons la valeur de la statistique de test à la table de la loi normale centrée-réduite, nous obtenons une p-valeur de 0.7465 >>> 5 \%. Nous en concluons que nous ne pouvons rejeter H_0 et donc que les médianes sont égales.

\bullet Application informatique:

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

Package et fonction R: http://finzi.psych.upenn.edu/R/library/RVAideMemoire/html/fp.test.html

\bullet Bibliographie:

– Robust Rank Procedures for the Behrens-Fisher Problem de M. A. Fligner et G. E. Policello

Le test de Fisher-Yates-Terry-Hoeffding

add

\bullet Présentation:

Publié en 1938 par Ronald Aylmer Fisher et Frank Yates, puis étendu par les travaux de Wassily Hoeffding en 1950 et de Terry en 1952, le test de Fisher-Yates-Terry-Hoeffding, également noté test de FYTH ou encore test des scores normés, est une approche non paramétrique permettant de tester si X|_{Y = 1}, X|_{Y = 2}, les sous-échantillons d’une variable continue X restreinte aux deux groupes d’une variable binaire Y, ont même distribution.

Lorsque les échantillons sont suffisamment grand, le test de Fisher-Yates-Terry-Hoeffding est réputé pour être bien plus efficace que le test de Wilcoxon-Mann-Whitney et garde l’immense avantage de ne pas être affecté par les outliers. Cette dernière propriété en fait un test particulièrement robustesse lorsque les données ne suivent pas une loi normale.

\bullet Le test:

Hypothèse préliminaire: X continue et Y binaire.

La formule de la statistique de test de Fisher-Yates-Terry-Hoeffding est,

C = \sum_{i = 1} ^{n_1} \Phi ^{-1} (F_i)

Où,

– le premier groupe sera celui ayant le plus petit effectif,

F_i = \frac{R_i - 0.375}{n + 0.25},

\Phi ^{-1} est la fonction inverse de la loi normale centrée-réduite.

Afin de pouvoir reporter la statistique de test à la table de la loi normale centrée-réduite, il faut lui appliquer la transformation,

Z = | \frac{C}{\frac{n_1 n_2}{\sqrt{n (n - 1)} \sum_{i = 1} ^n [\Phi ^{-1} (F_i)] ^2}} |

L’hypothèse H_0 est: « les deux sous-échantillons ont même distribution / F_{X|_{Y = 1}} = F_{X|_{Y = 2}} »

Ci-dessous la table de la loi normale centrée-réduite.

add

Tendance pour le rejet de H_0:

Plus la statistique de test Z est grande et plus nous avons de chance de rejeter H_0. Ce qui implique,

Z \rightarrow \infty \Rightarrow C \rightarrow \infty \Rightarrow \sum_{i = 1} ^{n_1} \Phi ^{-1} (\frac{R_i - 0.375}{n + 0.25}) \rightarrow \infty

Le test de Fisher-Yates-Terry-Hoeffding s’appuie sur la symétrie de la loi normale centrée-réduite pour étudier le lien entre Y et X. Ainsi, si les valeurs de X|_{Y = 1} sont toutes inférieures à celles de X|_{Y = 2} alors les valeurs renvoyées par la fonction inverse de la loi normale centrée-réduite seront toutes négatives et la somme sera maximale.

Dans le cas inverse, si nous sommes dans la situation d’une distribution aléatoire, alors ces valeurs vont s’annuler par symétrie et la somme tendra vers 0

\bullet Tendance lorsque n \longrightarrow \infty:

Nous proposons ici de vérifier si le test de Fisher-Yates-Terry-Hoeffding est sensible aux grands échantillons ou non.

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 distribution sont différentes d’un groupe à l’autre.

add.png

Globalement, quelque soit la taille de l’échantillon, le test statistique rejette H_0, ce qui est en accords avec nos hypothèses.

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

add1.png

Jusqu’à N = 100 nous restons cohérent avec nos hypothèses, malheureusement nous voyons qu’à partir d’un échantillon de taille 1000, le test rejette H_0 à tort.

Par conséquent, nous en déduisons que le test Fisher-Yates-Terry-Hoeffding est sensible à la taille de l’échantillon.

\bullet Annexe théorique:

Nous présentons ici quelques rappels sur les propriétés de symétrie de la loi normale centrée-réduite, principale caractéristique sur laquelle s’appuie le test de Fisher-Yates-Terry-Hoeffding.

La fonction inverse de la loi normale centrée-réduite correspond en réalité à la probabilité que l’évènement survienne,

P(X \leq a) = p, p \in [0,1[

Par définition et donc symétrie, nous avons également,

P(X \leq -u_{\alpha}) = P(X \geq u_{\alpha}) = \frac{\alpha}{2}

Et,

P(X \leq u_{\alpha}) = 1 - P(X \geq u_{\alpha}) = 1 - \frac{\alpha}{2}

La figure ci-dessous synthétise ces éléments.

add2.png

\bullet Exemple:

Soit la variable aléatoire X distribuée selon deux groupes d’une variable Y:

add

Ci-dessous, les densités associées aux distributions de X selon nos deux groupes. Nous pourrons remarquer que nos deux sous-échantillons ont la même distribution.

add

Dans un premier temps, nous déterminons le vecteur R des rangs associés à X,

R = (12,13,2,14,9,1,5,8,17,19,4,20,16,7,11,3,6,15,10,18)

Puis le vecteur F,

F = (\frac{12 - 0.375}{20 + 0.25}, \cdots, \frac{18 - 0.375}{20 + 0.25}) = (0.57407407,\cdots,0.87037037)

Nous calculons \Phi ^{-1} la transformation de F par la fonction inverse de la loi normale centrée-réduite,

\Phi ^{-1} = (1.8675612, \cdots 1.12814365)

Nos deux groupes ont même taille, par conséquent nous prendrons celui restreint à Y = 1. Nous pouvons désormais calculer la statistique de test de Fisher-Yates-Terry-Hoeffding,

C = \sum_{i = 1} ^{n_1} \Phi ^{-1} (F_i)

= 0.1867561 + 0.3145723 - 1.4034126 + 0.4477675 - 0.1867561 - 1.8682417 - 0.7441427 - 0.03145723 + 0.9191355 + 1.4034126

= -1.24581

Il nous reste à déterminer la variance,

\frac{10 \times 10}{20 (20 - 1)} \times \sum_{i = 1} ^{20} [\Phi ^{-1} (F_i)] ^2

= \frac{100}{380} \times (0.034877849 + 0.098955727 + 1.969567026 + 0.200495751 + 0.034877849 + 3.490326881 + 0.553748421 + 0.098955727 + 0.844810108 + 1.969567026 + 0.844810108 + 3.490326881 + 0.553748421 + 0.200495751 + 0.003835526 + 1.272708084 + 0.347458138 + 0.347458138 + 0.003835526 + 1.272708084)

= 0.2631579 \times 17.63357

= 4.640413

Enfin,

Z = | \frac{-1.245481}{\sqrt{4.640413}} | = | \frac{-1.245481}{2.154162} | = | 0.5781745 | = 0.5781745

Si nous reportons cette valeur de statistique de test à la table de la loi normale centrée-réduite, nous obtenons une p-valeur de 0.5631463 >>>> 5 \%. Nous déduisons que nous ne pouvons rejeter H_0 et donc que les deux distributions sont les mêmes.

\bullet Application informatique:

Procédure SAS:

%macro fyth_test(TAB=,Y=);

/* La macro suivante procède au test de Fisher-Yates-Terry-Hoeffding sur la table de données TAB en argumentant la variable groupe Y */

/* Options pour éviter de saturer la log */
options nonotes spool;
ods noresults;

/* Récupération des effectifs par classe */
proc freq data = &TAB.;
tables &Y. / out = n;
run;

/* Récupération du nom de la classe 1 et de son effectif */
data _null_;
set n;
if _n_ = 1;
call symputx(« n1 »,count);
call symputx(« cl1 »,&Y.);
run;

/* Récupération de l’effectif de la classe 2 */
data _null_;
set n;
if _n_ = 2;
call symputx(« n2 »,count);
run;

proc contents data = &TAB. out = varX;
run;

/* Récupération du nom de la variable continue */
data _null_;
set varX;
if name ne « &Y. »;
call symputx(« varX »,name);
run;

/* Calcul des rangs associés à la variable continue */
proc rank data = &TAB. out = Phi;
var &varX.;
ranks R;
run;

/* Création de la feuille de calcul de Phi */
data Phi;
set Phi;
F = (R-0.375)/(&n1.+&n2.+0.25);
Phi = quantile(‘NORMAL’,F);
Phi2 = Phi**2;
run;

/* Calcul de la somme des phis pour les observations de la classe 1 et de la somme des phis au carré */
proc sql noprint;
select sum(Phi) into: Phi1 from Phi where &Y. = « &cl1. »;
select sum(Phi2) into: Phi2 from Phi;
quit;

/* Création de la table de résultats */
data Resultats;
Pop = &n.;
FYTH = &Phi1.;
Var = &Phi2.*(&n1.*&n2.)/((&n1.+&n2.)*(&n1.+&n2.-1));
Z = abs(FYTH/sqrt(Var));
p = 1 – probnorm(Z);
run;

/* Suppression des tables inutiles */
proc datasets lib = work nolist;
delete n varX Phi;
run;

/* Reset des options */
options notes nospool;
ods results;

%mend;

Package et fonction R:

FYTH.test = function(TAB) {

# La fonction suivante applique le test de Fisher-Yates-Terry-Hoeffding sur la base de données TAB contenant en première colonne la variable des groupes et en seconde colonne la variable continue

# Epuration des données manquantes
TAB = na.omit(TAB)
# Récupération du nombre d’osbservations n
n = dim(TAB)[1]
BiblioY = summary(factor(TAB[,1]))
# Récupération du nom de la classe de référence Cl1
Cl1 = names(BiblioY)[1]
# Récupération des effectifs par classe
n1 = BiblioY[[1]]
n2 = BiblioY[[2]]

# Application de la fonction de répartition normale inverse
Phi = matrix(0,n,1)
R = rank(TAB[,2])
for (i in 1:n) {Phi[i] = qnorm((R[i]-0.375)/(n+0.25))}

# Calcul de la statistique de test de FYTH
C = sum(Phi[which(TAB[,1]==Cl1)])
# Calcul de la variance
D = sum(Phi^2)*(n1*n2)/(n*(n-1))

# Calcul de la statistique de test
Z = abs(C/sqrt(D))
# Calcul de la p-valeur
p = 2*(1-pnorm(Z,0,1))

# Mise en forme et impression des résultats
return(data.frame(Pop = n, FYTH = C, Var = D, Z = Z, p_valeur = p))

}

\bullet Bibliographie:

– Theory of rank tests de Zybnek Sidak, Prenab K. Sen et Jaroslav Hajek

– Comparaisons de population. Tests non paramétriques de Ricco Rakotomalala

– Le document: http://www.maths-france.fr/Terminale/TerminaleS/Cours/18-loi-normale.pdf

Le test de Quade

add

\bullet Présentation:

Publié en 1979 par Dana Quade, le test de Quade, appelé également ANalyse De VAriance (ANOVA) de Quade, est une approche non paramétrique permettant de tester si plusieurs variables continues appariées X ^1, \cdots, X ^T sont liées.

Le test de Quade est considéré comme une alternative à l’ANOVA de Fisher et demeure plus puissant que l’ANOVA de Friedman lorsque T < 5.

\bullet Le test:

Hypothèse préliminaire: Variables continues appariées.

Le calcul de la statistique de test nécessite au préalable la création de la matrice des rangs en ligne \mathbf{R} associée à \mathbf{X} = (X ^1, \cdots, X ^T) ainsi que le calcul du vecteur Q des rangs associés au vecteur des différences,

(max(X_i ^t) - min(X_i ^t)), \forall i \in [1,n]

La formule de la statistique de test de Quade est alors,

F = \frac{(n - 1) SS_T}{SS_{total} - SS_T}

Où,

SS_{total} = \sum_{i = 1} ^n \sum_{t = 1} ^T (S_i ^t) ^2

SS_T = \frac{1}{n} \sum_{t = 1} ^T (\sum_{i = 1} ^n S_i ^t) ^2

S_i ^t = Q_i (R_i ^t - \frac{T + 1}{2})

La statistique de test de Quade suit une loi de Fisher-Snedecor à (T - 1, (T  - 1) (n - 1)) degrés de liberté. L’hypothèse H_0 est: « Il n’y a pas d’influence du temps T sur X / \sigma_1 = \cdots \sigma_T« .

Ci-dessous la table de Fisher-Snedecor:

addTendance pour le rejet de H_0:

Plus la statistique de test F est grande et plus nous avons de chance de rejeter H_0. Ce qui revient à dire,

F = \frac{(n - 1) SS_T}{SS_{total} - SS_T} \rightarrow \infty \Rightarrow SS_T \rightarrow \infty ou SS_{total} - SS_T \rightarrow 0

L’élément qui nous intéresse est SS_T puisque SS_{total} est indépendant de la configuration des données. En effet, ces deux objets sont fonctions des S_i ^t = Q_i (R_i ^t - \frac{T + 1}{2}) dont la réelle valeur informative provient de la considération de la somme de ces éléments en fonction de T, approche considérée exclusivement lors du calcul de SS_T.

Par conséquent, si nous étudions le numérateur, SS_T est maximal si,

\sum_{i = 1} ^n S_i ^1 = 1 \times n, \sum_{i = 1} ^n S_i ^2 = 2 \times n, \cdots, \sum_{i = 1} ^n S_i ^{n - 1} = (T - 1) \times n, \sum_{i = 1} ^n S_i ^n = T \times n

(où l’inverse), soit le cas où il y a une dépendance entre X et T.

L’ Autre cas de figure est celui où SS_T \rightarrow SS_{total}. Soit que la dispersion inter-traitements est également à la dispersion totale au travers de l’inégalité triangulaire,

SS_T = \frac{1}{n} \sum_{t = 1} ^T (\sum_{i = 1} ^n S_i ^t) ^2 \leq \frac{1}{n} \sum_{t = 1} ^T \sum_{i = 1} ^n (S_i ^t) ^2 = \frac{1}{n} SS_{total} \Rightarrow SS_T \leq \frac{1}{n} SS_{total}

Or nous avons,

SS_{total} = \sum_{i = 1} ^n \sum_{t = 1} ^T (S_i ^t) ^2 = n \sum_{t = 1} ^T t ^2

Et si nous nous retrouvons dans le cas idéal décrit ci-dessus,

SS_T = \sum_{t = 1} ^T (\sum_{i = 1} ^n S_i ^t)  ^2 = \sum_{t = 1} ^T (n \cdot t) ^2 = n ^2 \sum_{t = 1} ^T t ^2

Nous obtenons alors l’égalité SS_T = \frac{1}{n} SS_{total} lorsqu’il y a une influence totale de T sur X. En fin de compte, SS_{total} traduit cette configuration des données auquel nous cherchons à comparer la configuration observée.

\bullet Tendance lorsque n \longrightarrow \infty:

Nous proposons ici de vérifier si le test de Quade est sensible aux grands échantillons ou non.

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 distribution sont différentes d’un temps à l’autre.

add1

Globalement, quelque soit la taille de l’échantillon, le test statistique rejette H_0, ce qui est en accords avec nos hypothèses.

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

add2

Jusqu’à N = 100 nous restons cohérent avec nos hypothèses, cependant nous voyons qu’à partir d’un échantillon de taille 1 000, le test rejette H_0 à tort.

Nous en déduisons que le test de Quade est influencé par la taille de l’échantillon.

\bullet Annexe théorique:

Nous présentons ici une esquisse de la démonstration de l’inégalité triangulaire.

Trivialement nous avons,

\forall j_1 < j_2, X ^{j_1} X ^{j_2} \leq | X ^{j_1} X ^{j_2}|

\Rightarrow \forall j_1 < j_2, 2 X ^{j_1} X ^{j_2} \leq 2 | X ^{j_1} X ^{j_2} |

\Rightarrow 2 \sum_{j_1 < j_2} X ^{j_1} X ^{j_2} \leq 2 \sum_{j_1 < j_2} X ^{j_1} X ^{j_2}

\Rightarrow \prod_j (X ^j) ^2 + 2 \sum_{j_1 < j_2} X ^{j_1} X ^{j_2} \leq \prod_j (X ^j) ^2 + 2 \sum_{j_1 < j_2} | X ^{j_1} X ^{j_2} |

\Rightarrow \prod_j (X ^j) ^2 + 2 \sum_{j_1 < j_2} X ^{j_1} X ^{j_2} \leq \prod_j | X ^j | ^2 + 2 \sum_{j_1 < j_2} | X ^{j_1} X ^{j_2} |

\Rightarrow (\sum_j X ^j) ^2 \leq (\sum_j |X ^j|) ^2

\Rightarrow | \sum_j X ^j | ^2 \leq (\sum_j |X ^j|) ^2

\Rightarrow | \sum_j X ^j | \leq \sum_j |X ^j|

\bullet Exemple:

Soit les cinq échantillons appariés X ^1, X ^2, X ^3, X ^4, X ^5 suivants.

add

Ci-dessous, les boxplots associés à (X ^1, X ^2, X ^3, X ^4, X ^5).

add

Les boxplots montrent que X varie énormément en fonction du temps t \in \lbrace 1, 2, 3, 4, 5 \rbrace. Prouvons-le statistiquement.

Dans un premier temps, nous déterminons le tableau des rangs en ligne associés à X ^1, X ^2, X ^3, X ^4, X ^5,

add

Ensuite, il nous faut déterminer le vecteur des différences entre la valeur maximale et la valeur minimale pour chaque block,

(9.2858 - 0.8970, \cdots, 20.1656 - 9.5949) = (8.3888, \cdots, 10.5707)

Et qui nous permet de déterminer l’objet Q,

Q = (8,9,10,2,3,7,5,6,4,1,18,16,19,14,13,15,17,20,12,11)

Nous pouvons désormais calculer la matrice pondérée de Quade,

S = \begin{pmatrix} 8 \times (1 - \frac{5 + 1}{2}) & 8 \times (4 - \frac{5 + 1}{2}) & 8 \times (3 - \frac{5 + 1}{2}) & 8 \times (2 - \frac{5 + 1}{2}) & 8 \times (5 - \frac{5 + 1}{2}) \\ 9 \times (1 - \frac{5 + 1}{2}) & 9 \times (4 - \frac{5 + 1}{2}) & 9 \times (3 - \frac{5 + 1}{2}) & 9 \times (2 - \frac{5 + 1}{2}) & 9 \times (5 - \frac{5 + 1}{2}) \\ 10 \times (2 - \frac{5 + 1}{2}) & 10 \times (1 - \frac{5 + 1}{2}) & 10 \times (4 - \frac{5 + 1}{2}) & 10 \times (3 - \frac{5 + 1}{2}) & 10 \times (5 - \frac{5 + 1}{2}) \\ 2 \times (1 - \frac{5 + 1}{2}) & 2 \times (4 - \frac{5 + 1}{2}) & 2 \times (3 - \frac{5 + 1}{2}) & 2 \times (2 - \frac{5 + 1}{2}) & 2 \times (5 - \frac{5 + 1}{2}) \\ 3 \times (1 - \frac{5 + 1}{2}) & 3 \times (4 - \frac{5 + 1}{2}) & 3 \times (3 - \frac{5 + 1}{2}) & 3 \times (2 - \frac{5 + 1}{2}) & 3 \times (5 - \frac{5 + 1}{2}) \\ 7 \times (4 - \frac{5 + 1}{2}) & 7 \times (1 - \frac{5 + 1}{2}) & 7 \times (3 - \frac{5 + 1}{2}) & 7 \times (5 - \frac{5 + 1}{2}) & 7 \times (2 - \frac{5 + 1}{2}) \\ 5 \times (4 - \frac{5 + 1}{2}) & 5 \times (2 - \frac{5 + 1}{2}) & 5 \times (3 - \frac{5 + 1}{2}) & 5 \times (5 - \frac{5 + 1}{2}) & 5 \times (1 - \frac{5 + 1}{2}) \\ 6 \times (4 - \frac{5 + 1}{2}) & 6 \times (3 - \frac{5 + 1}{2}) & 6 \times (1 - \frac{5 + 1}{2}) & 6 \times (5 - \frac{5 + 1}{2}) & 6 \times (2 - \frac{5 + 1}{2}) \\ \ldots & \ldots & \ldots & \ldots & \ldots \\ \end{pmatrix}

= \begin{pmatrix} -16 & 8 & 0 & -8 & 16 \\ -18 & 9 & 0 & -9 & 18 \\ -10 & -20 & 10  & 0 & 20 \\ -4 & 2 & 0 & -2 & 4 \\ -6 & 3 & 0 & -3 & 6 \\ 7 & -14 & 0 & 14 & -7 \\ 5 & -5 & 0 & 10 & -10 \\ 6 & 0 & -12 & 12 & -6 \\ 0 & 4 & -4 & 8 & -8 \\ 1 & 0 & -1 & 2 & -2 \\ 18 & -18 & -36 & 36 & 0 \\ 16 & 0 & -32 & 32 & -16 \\ 19 & -19 & -38 & 38 & 0 \\ 14 & -28 & -14 & 28 & 0 \\ 13 & -26 & -13 & 26 & 0 \\ 15 & -30 & -15 & 0 & 30 \\ 17 & -17 & -34 & 0 & 34 \\ 20 & -20 & -40 & 0 & 40 \\ 12 & -24 & -12 & 0 & 24 \\ 11 & -22 & 0 & -11 & 22 \\ \end{pmatrix}

Nous avons alors,

SS_{Total} = \sum_{i = 1} ^{20} \sum_{t = 1} ^5 (S_i ^t) ^2 = (-16) ^2 + 8 ^2 + 0 ^2 + \cdots + 0 ^2 + (-11) ^2 + 22 ^2 = 28700

Et,

SS_T = \frac{1}{20} \times \sum_{t = 1} ^5 (\sum_{i = 1} ^{20} S_i ^t) ^2 = \frac{120 ^2 + (-217) ^2 + (-241) ^2 + 173 ^2 + 165 ^2}{5} = \frac{176724}{20} = 8836.2

Nous avons désormais tous les éléments pour le calcul de la statistique de test de Quade,

F = \frac{(20 - 1) \times 8836.2}{28700 - 8836.2} = 8.451948

Si nous reportons la valeur de la statistique de test à la loi de Fisher-Snedecor pour (5 - 1, (5 - 1) \times (20 - 1)) = (4,76) degrés de liberté, nous obtenons une p-valeur de 0.00001073 <<<<< 5 \%. Nous rejetons alors H_0 et pouvons en conclure que le temps T a une influence sur X.

\bullet Application informatique:

Procédure SAS:

%macro quade_test (TAB=);

/* Macro pour procéder au calcul du test de Quade sur une matrice à T colonnes et de format continue */

/* Options pour épurer les sorties */
options nonotes spool;
ods noresults;

/* Création de la matrice des rangs en ligne */
proc transpose data = &TAB. out = R;
run;

/* Aparthé pour la futur création de Q */
proc means data = R min max;
output out = min_max;
run;

proc rank data = R out = R;
run;

proc transpose data = R out = R;
run;

/* Création de Q */
data min_max;
set min_max;
if _stat_ in (« MIN » « MAX »);
drop _type_ _freq_;
run;

proc transpose data = min_max out = min_max;
id _stat_;
run;

data Q;
set min_max;
D = max – min;
keep D;
run;

proc rank data = Q out = Q;
run;

/* Récupération du paramètre n: nombre de blocks */
proc sql noprint;
select count(*) into: n from &TAB.;
quit;

proc contents data =  &TAB. out = nb_vars;
run;

/* Récupération de la liste des variables appariées et du paramètre T: nombre de temps */
proc sql noprint;
select distinct(name) into: list_vars separated by  »  » from nb_vars;
select count(*) into: T from nb_vars;
quit;

/* Initialisation de la matrice S */
data S;
run;

/* Remplissage de la matrice S avec la formule adéquate */
%let SStot = 0;

%do i = 1 %to &n.;

data ajout;
run;

data _null_;
set Q;
if _n_ = &i.;
call symputx(« Qi »,D);
run;

%do tt = 1 %to &T.;

data _null_;
set R;
if _n_ = &i.;
call symputx(« Rit »,%scan(&list_vars.,&tt.));
run;

%let Sit = %sysevalf(&Qi.*(&Rit. – ((&T. + 1)/2)));

data new;
%scan(&list_vars.,&tt.) = &Sit.;
run;

data ajout;
merge ajout new;
run;

/* Calcul du SStotal par méthode itérative */
%let SStot = %sysevalf(&SStot. + (&Sit.)**2);

%end;

data S;
set S ajout;
run;

%end;

/* Epuration de la première ligne d’initialisation de S */
data S;
set S;
if _n_ ne 1;
run;

/* Calcul du SST */
proc transpose data = S out = SST;
run;

data SST;
set SST;
drop _name_;
run;

proc contents data = SST out = sum;
run;

proc sql noprint;
select (name) into: list_sum separated by « , » from sum;
quit;

data SST;
set SST;
Sum = sum(&list_sum.);
run;

proc transpose data = SST out = SST;
run;

/* Calcul du SST par méthode itérative */
data SST;
set SST;
SST = 0;
%do tt = 1 %to &T.;
SST = SST + COL&tt.**2;
%end;
SST = SST/&n.;
run;

data _null_;
set SST;
if _n_ = &n. + 1;
call symputx(« SST »,SST);
run;

/* Mise en forme des résultats */
data Resultats;
n = &n.;
T = &T.;
F = %sysevalf(((&n. – 1) * &SST.)/(&SStot. – &SST.));
SST = &SST.;
SStotal = &SStot.;
df1 = &T. – 1;
df2 = (&T. – 1)*(&n. – 1);
p = 1 – Probf(F,df1,df2);
run;

/* Suppression des tables annexes inutiles */
proc datasets lib = work nolist;
delete R Q min_max nb_vars ajout new SST sum;
run;

/* Reset des options */
options notes nospool;
ods results;

%mend;

Package et fonction R: https://stat.ethz.ch/R-manual/R-devel/library/stats/html/quade.test.html

\bullet Bibliographie: 

– Using Weighted Rankings in the Analysis of Complete Blocks with Additive Block Effects de Dana Quade

– New Statistical Procedures for the Social Sciences: Modern Solutions To Basic Problems de Rand R. Wilcox

– Biostatistics and Microbiology: A Survival Manual de Daryl S. Paulson