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