Logiciel R, séances de perfectionnement
en 4 demi-journées
Séance 2 : Statistiques inférentielles (tests)
gilles.hunault "at" univ-angers.fr
Table des matières cliquable
1. Fonctions usuelles de R pour les tests
2. Comparaisons de pourcentages
3. Comparaison d'une moyenne à une référence
4. Comparaisons de moyennes et ANOVA, paramétrique ou non paramétrique
Il est possible d'afficher toutes les solutions via ?solutions=1 et de toutes les masquer avec via ?solutions=0.
1. Fonctions usuelles de R pour les tests
Rappeler quels sont les "grands tests classiques", paramétriques ou non. Quelles sont les fonctions usuelles de R associées ?
Comment récupère-t-on la p-value ? Et la valeur de la statistique de test ? Et son intervalle de confiance ?
Application : Donner l'intervalle de confiance de la variable Length pour les données lea.dar.
Solution : afficher la solution
Pour comparer deux QT, on utilise un test d'hypothèse paramétrique ou non paramétrique suivant la normalité de la distribution des échantillons : pour des données normales non appariées, on utilise le test t de Student qui est équivalent à une ANOVA à un facteur ; pour des données normales appariées, on utilise le test t de Student adapté aux données appariées ; pour des données non normales non appariées, on utilise le test de Mann-Whitney ; pour des données non normales appariées, on utilise le test de Wilcox.
Voici le tableau résumé de ces divers cas pour comparer deux QT :
Données non appariées Données appariées Données normales test t de Student (ou anova) test t de Student adapté Données non normales test de Mann-Whitney test de Wilcox Pour comparer plus de deux QT, on utilise comme tests d'hypothèse paramétriques des anova, ancova, manova et autre mancova (donc pour des données normales) ; dans le cas non paramétrique, on utilise le test de Kruskal-Wallis pour des données non normales non appariées et le test de Friedman pour des données non normales appariées.
Voici le tableau résumé de ces divers cas pour comparer plus de deux QT :
Données non appariées Données appariées Données normales anova, ancova etc. anova, ancova etc. Données non normales test de Kruskal-Wallis test de Friedman En termes de graphiques, il est d'usage de tracer des boites à moustaches simultanées des QT, avec éventuellements des encoches (qui illustrent le test d'égalité des médianes).
Il faut noter que, pour plus de deux QT, lorsque les tests d'hypothèses montrent des différences significatives entre les QT, on doit effectuer des tests post hoc pour déterminer quelles QT sont significativement différentes. De plus, ces comparaisons s'effectuent plutôt sur des groupes à l'intérieur d'une même QT que sur des QT différentes.
Pour étudier conjointement deux QL, on effectue un tri croisé puis un test d'indépendance du khi-deux sur ce tri croisé si les effectifs présents le permettent. En termes de graphiques, il est d'usage de tracer les histogrammes des effectifs des deux tris à plats et du tri croisé dans les deux sens afin de détecter d'éventuelles dépendances (changements de profils).
Voici un tableau résumé avec les principaux tests et les fonctions R associées :
t.test() test de Student wiki correspondant wilcox.test() test de Mann-Whitney wiki correspondant wilcox.test() test de Wilcox wiki correspondant kruskal.test() test de Kruskal-Wallis wiki correspondant friedman.test() test de Friedman wiki correspondant chisq.test() test du chi-deux wiki correspondant Les tests non paramétriques s'utilisent lorsque les conditions de normalité, d'égalité des variances homoscédasticité etc. ne sont pas respectées. Leur utilisation est semblable à celles des test paramétriques. Ainsi wilcox.test(x,y) se substitue à t.test(x,y) et friedman.test(y ~ x) est semblable à anova(lm(y ~ x)).
Lorsqu'on connait la fonction de R qui réalise le test demandé, il faut regarder le nom des composantes de la liste obtenue par application de cette fonction. Les noms sont obtenus via la fonction names(). On peut alors extraire avec le symbole dollar la composante désirée :
# lecture des données url <- "http://www.info.univ-angers.fr/~gh/Datasets/lea.dar" lea <- read.table(url,header=TRUE,row.names=1) lng <- lea$length # intervalle de confiance via la fonction t.test() print(t.test(lng)) # récupération de cet intervalle rdt <- t.test(lng) print(names(rdt)) cat("\n l'intervalle de confiance est donc ",rdt$conf.int,"\n") # récupération de la p-value pval <- rdt$p.value cat("\nla p-value est ",pval,"\n") # avec nos fonctions : # chargement des fonctions gh source("http://forge.info.univ-angers.fr/~gh/wstat/statgh.r",encoding="latin1") icmQT(lng)One Sample t-test data: lng t = 38.4601, df = 772, p-value < 2.2e-16 alternative hypothesis: true mean is not equal to 0 95 percent confidence interval: 195.1897 216.1868 sample estimates: mean of x 205.6882 [1] "statistic" "parameter" "p.value" "conf.int" "estimate" "null.value" "alternative" "method" "data.name" l'intervalle de confiance est donc 195.1897 216.1868 la p-value est 1.376861e-181 (gH) version 4.87 fonctions d'aides : lit() fqt() fql() ic() fapprox() chi2() fcomp() datagh() taper aide() pour revoir cette liste [1] 195.1897 216.1868
2. Comparaisons de pourcentages
La parité Femme/Homme est-elle respectée dans les données elf.dar ? On trouvera des informations sur ces données à l'URL elf.htm.
Peut-on considérer que les différents règnes sont représentés équitablement dans les données lea.dar ? Oon trouvera des informations sur ces données à l'URL lea.htm.
Solution : afficher la solution
L'équirépartition peut se faire via un test de chi-deux sous hypothèse d'indépendance des deux modalités Homme/Femme de la variable qualitative SEXE. Intuitivement, la répartition en à peu près deux tiers de femmes, un tiers d'hommes laisse à penser qu'il n'y a pas égalité des proportions. Attention : cette intuition peut être trompeuse puisque le test utilise le nombre total d'individus. Voir par exemple la fin de la page pourquoi.
Voici comment avoir rapidement le détail des calculs pour le vérifier avec notre fonction chi2Adeq() ou avec les fonctions classiques de R nommées prop.test() et binom.test().
# chargement de nos fonctions source("http://www.info.univ-angers.fr/~gh/statgh.r",encoding="latin1") # lecture des données elf <- lit.dar("http://www.info.univ-angers.fr/~gh/Datasets/elf.dar") # tri à plat de la variable SEXE decritQL("SEXE",elf$SEXE,c("Homme","Femme"),TRUE) sexe <- factor(elf$SEXE,levels=0:1,labels=c("Homme","Femme")) decritQLf("SEXE",sexe,TRUE,testequi=TRUE) # test de l'égalité des proportions (3 moyens) chi2Adeq(vth=c(49,50),vobs=c(35,64)) print( binom.test(x = 64,n = 99) ) print( prop.test(x = 64,n = 99) ) # graphiques pour comparer les distributions # théoriques et pratiques countRes <- t(cbind(table(elf$SEXE),c(49,50))) echantillons <- c("ELF","Théorique") row.names(countRes) <- echantillons modalites <- c("Homme","Femme") colnames(countRes) <- modalites couleurs <- c("blue","red") barplot(countRes,beside = TRUE,col=couleurs,main="Répartition SEXE",ylim=c(0,100)) legend(x=1,y=100,fill=couleurs,legend=echantillons)(gH) version 4.87 fonctions d'aides : lit() fqt() fql() ic() fapprox() chi2() fcomp() datagh() taper aide() pour revoir cette liste TRI A PLAT DE : SEXE (ordre des modalités) Homme Femme Total Effectif 35 64 99 Cumul Effectif 35 99 99 Frequence (en %) 35 65 100 Cumul fréquences 35 100 100 VARIABLE : SEXE (par fréquence décroissante) Femme Homme Total Effectif 64 35 99 Cumul Effectif 35 99 99 Frequence (en %) 65 35 100 Cumul fréquences 65 100 100 [1] "Homme" "Femme" TRI A PLAT DE : SEXE (ordre des modalités) Homme Femme Total Effectif 35 64 99 Cumul Effectif 35 99 99 Frequence (en %) 35 65 100 Cumul fréquences 35 100 100 VARIABLE : SEXE (par fréquence décroissante) Femme Homme Total Effectif 64 35 99 Cumul Effectif 35 99 99 Frequence (en %) 65 35 100 Cumul fréquences 65 100 100 [1] "Homme" "Femme" CALCUL DU CHI-DEUX D'ADEQUATION Valeurs théoriques 49.5 49.5 Valeurs observées 35 64 Valeur du chi-deux 8.494949 Chi-deux max (table) à 5 % 3.841459 pour 1 degrés de liberté ; p-value 0.003561337 décision : au seuil de 5 % on peut rejeter l'hypothèse que les données observées correspondent aux valeurs théoriques. Contributions au chi-deux Ind. The Obs Dif Cntr Pct Cumul 1 49.5 35 14.5 4.247475 50 4.247475 2 49.5 64 -14.5 4.247475 50 8.494949 CALCUL DU CHI-DEUX D'ADEQUATION Valeurs théoriques 49 50 Valeurs observées 35 64 Valeur du chi-deux 7.92 Chi-deux max (table) à 5 % 3.841459 pour 1 degrés de liberté ; p-value 0.004889127 décision : au seuil de 5 % on peut rejeter l'hypothèse que les données observées correspondent aux valeurs théoriques. Exact binomial test data: 64 and 99 number of successes = 64, number of trials = 99, p-value = 0.004641 alternative hypothesis: true probability of success is not equal to 0.5 95 percent confidence interval: 0.5439647 0.7399466 sample estimates: probability of success 0.6464646 1-sample proportions test with continuity correction data: 64 out of 99, null probability 0.5 X-squared = 7.9192, df = 1, p-value = 0.004891 alternative hypothesis: true p is not equal to 0.5 95 percent confidence interval: 0.5432921 0.7381205 sample estimates: p 0.6464646Pour les règnes dans lea.dar le plus simple est de passer par notre fonction decritQLf() avec son paramètre de test d'équirépartition mis à vrai (TRUE) :
# lecture des données url <- "http://www.info.univ-angers.fr/~gh/Datasets/lea.dar" lea <- read.table(url,header=TRUE,row.names=1) reign <- lea$reign # description du facteur et test d'équirépartition decritQLf(titreQL = "Règne",nomFact = reign,graphique = TRUE,testequi = TRUE)TRI A PLAT DE : Règne (ordre des modalités) Alveolata Bacteria Euryarchaeota Fungi Metazoa Parabasalidea Viridiplantae Total Effectif 1 38 1 11 23 1 698 773 Cumul Effectif 1 39 40 51 74 75 773 773 Frequence (en %) 0 5 0 1 3 0 90 99 Cumul fréquences 0 5 5 7 10 10 100 99 VARIABLE : Règne (par fréquence décroissante) Viridiplantae Bacteria Metazoa Fungi Alveolata Euryarchaeota Parabasalidea Total Effectif 698 38 23 11 1 1 1 773 Cumul Effectif 1 39 40 51 74 75 773 773 Frequence (en %) 90 5 3 1 0 0 0 99 Cumul fréquences 90 95 98 100 100 100 100 99 [1] "Alveolata" "Bacteria" "Euryarchaeota" "Fungi" "Metazoa" "Parabasalidea" "Viridiplantae" CALCUL DU CHI-DEUX D'ADEQUATION Valeurs théoriques 110.4286 110.4286 110.4286 110.4286 110.4286 110.4286 110.4286 Valeurs observées 1 38 1 11 23 1 698 Valeur du chi-deux 3657.928 Chi-deux max (table) à 5 % 12.59159 pour 6 degrés de liberté ; p-value 0 décision : au seuil de 5 % on peut rejeter l'hypothèse que les données observées correspondent aux valeurs théoriques. Contributions au chi-deux Ind. The Obs Dif Cntr Pct Cumul 1 110.4286 1 109.42857 108.4376 2.964455 108.4376 2 110.4286 38 72.42857 47.5049 1.298683 155.9425 3 110.4286 1 109.42857 108.4376 2.964455 264.3802 4 110.4286 11 99.42857 89.5243 2.447405 353.9045 5 110.4286 23 87.42857 69.2190 1.892301 423.1235 6 110.4286 1 109.42857 108.4376 2.964455 531.5611 7 110.4286 698 -587.57143 3126.3665 85.468245 3657.9276On en déduit (évidemment !) qu'il n'y a pas équirépartition, le règle très majoritaire étant celui des plantes (Viridiplantae).
3. Comparaison d'une moyenne à une référence
Après avoir analysé la variable Length des données lea.dar (on trouvera des informations sur ces données à l'URL lea.htm), peut-on affirmer qu'en moyenne une protéine LEA comporte 361 résidus d'acides aminés ? Au passage, d'où provient cette valeur 361 ?
Solution : afficher la solution
Avant de réaliser ce test de conformité de moyenne à une valeur théorique, cela ne fait pas de mal d'analyser la variable Length :
# lecture des données source("http://www.info.univ-angers.fr/~gh/statgh.r",encoding="latin1") lea <- lit.dar("http://www.info.univ-angers.fr/~gh/Datasets/lea.dar") # description des données decritQT("Longueur des protéines de LEAPDB",lea$length ,"aa",TRUE) # test de conformité print(t.test(lea$length,mu=361))(gH) version 4.87 fonctions d'aides : lit() fqt() fql() ic() fapprox() chi2() fcomp() datagh() taper aide() pour revoir cette liste DESCRIPTION STATISTIQUE DE LA VARIABLE Longueur des protéines de LEAPDB Taille 773 individus Moyenne 205.6882 aa Ecart-type 148.6924 aa Coef. de variation 72 % 1er Quartile 130.0000 aa Mediane 168.0000 aa 3eme Quartile 236.0000 aa iqr absolu 106.0000 aa iqr relatif 63.0000 % Minimum 68 aa Maximum 1864 aa Tracé tige et feuilles The decimal point is 2 digit(s) to the right of the | 0 | 77777888888899999999999999999999999999999999999999999999999999999999 1 | 00000000000000000000000000000000000000000000000111111111111111111111+323 2 | 00000000000000000000000000001111111111111111222222222222223333333333+107 3 | 000000011111222223333344444445556666788889999 4 | 0011111112223444556666777778889 5 | 0122455667888 6 | 22344588 7 | 234 8 | 5 9 | 10 | 11 | 12 | 4 13 | 14 | 3 15 | 1 16 | 17 | 18 | 6 One Sample t-test data: lea$length t = -29.0406, df = 772, p-value < 2.2e-16 alternative hypothesis: true mean is not equal to 361 95 percent confidence interval: 195.1897 216.1868 sample estimates: mean of x 205.6882La p-value du test de conformité très inférieure à 0,05 indique qu'on doit rejeter l'hypothèse selon laquelle en moyenne une protéine LEA comporte 361 résidus d'acides aminés.
La valeur 361 (qui est en fait une valeur de médiane et non pas de moyenne, ce qui est sans doute plus "adapté") provient de l'article très intéressant de Brocchieri et Karlin (Nucleic Acids Research, 2005) -- dont une copie locale est ici -- dans lequel on peut lire : The median length of the proteins annotated among Eukaryotes (361 amino acids) is much higher than in Bacteria (267 amino acids).
4. Comparaisons de moyennes et ANOVA, paramétrique ou non paramétrique
On voudrait comparer les distributions d'ages pour les femmes et les hommes dans les données elf.dar (on trouvera des informations sur ces données à l'URL elf.htm). Quels tests envisager ? Comment choisir ici entre tests paramétriques et non paramétriques ?
Solution : afficher la solution
Il y a bien sûr plusieurs façons de constituer le vecteur des ages des hommes et celui des femmes ou les sous-ensembles correspondants :
# lecture des données source("http://www.info.univ-angers.fr/~gh/statgh.r",encoding="latin1") elf <- lit.dar("http://www.info.univ-angers.fr/~gh/Datasets/elf.dar") # décomposition en sous-groupes d'hommes et de femmes # 1. par filtre explicite dans le vecteur AGE agesFemmes <- elf$AGE[ elf$SEXE==1 ] agesHommes <- elf$AGE[ elf$SEXE==0 ] cat(" ",length(agesFemmes)," femmes et ",length(agesHommes)," hommes\n\n") # 2. à l'aide de split on crée une liste avec une composante pour les hommes # et une composante pour les femmes decoupe <- with(elf, split(x=AGE,f=SEXE) ) names(decoupe) <- c("Hommes","Femmes") print(sapply(X=decoupe,F=length)) # 3. avec subset # sous-ensemble 1 cat("\nFemmes\n") elfFemmes <- subset(elf,SEXE==1) print(dim(elfFemmes)) # sous-ensemble 2 cat("\nHommes\n") elfHommes <- subset(elf,SEXE==0) print(dim(elfHommes))(gH) version 4.87 fonctions d'aides : lit() fqt() fql() ic() fapprox() chi2() fcomp() datagh() taper aide() pour revoir cette liste 64 femmes et 35 hommes Hommes Femmes 35 64 Femmes [1] 64 6 Hommes [1] 35 6Mais en fait, grâce à la notation ~ (tilde) qui se lit «en fonction de», il n'y a pas besoin de décomposer quoi que ce soit :
# conversion explicite du sexe en variable qualitative ("factor") elf$SEXE <- factor(elf$SEXE,levels=0:1,labels=c("Hommes","Femmes")) # test paramétrique print( t.test(AGE~SEXE,data=elf,var.equal=FALSE) ) print( t.test(AGE~SEXE,data=elf,var.equal=TRUE) ) print( anova(lm(AGE~SEXE,data=elf) )) # test non-paramétrique print( wilcox.test(AGE~SEXE,data=elf) )Welch Two Sample t-test data: AGE by SEXE t = 0.2431, df = 74.04, p-value = 0.8086 alternative hypothesis: true difference in means is not equal to 0 95 percent confidence interval: -6.363368 8.132118 sample estimates: mean in group Hommes mean in group Femmes 36.40000 35.51562 Two Sample t-test data: AGE by SEXE t = 0.2385, df = 97, p-value = 0.812 alternative hypothesis: true difference in means is not equal to 0 95 percent confidence interval: -6.475006 8.243756 sample estimates: mean in group Hommes mean in group Femmes 36.40000 35.51562 Analysis of Variance Table Response: AGE Df Sum Sq Mean Sq F value Pr(>F) SEXE 1 17.7 17.696 0.0569 0.812 Residuals 97 30176.4 311.097 Wilcoxon rank sum test with continuity correction data: AGE by SEXE W = 1195, p-value = 0.5854 alternative hypothesis: true location shift is not equal to 0Il serait certainement plus raisonnable de commencer par analyser la variable AGE globalement et par sous-groupe avant de décider quoi que ce soit :
# lecture des données source("http://www.info.univ-angers.fr/~gh/statgh.r",encoding="latin1") elf <- lit.dar("http://www.info.univ-angers.fr/~gh/Datasets/elf.dar") # 1. analyse globale de la variable AGE decritQT("AGE (global)",elf$AGE,"ans",TRUE) # 2. analyse par sous-groupe decritQT("AGE (femmes)",elf$AGE[elf$SEXE==1],"ans",TRUE) decritQT("AGE (hommes)",elf$AGE[elf$SEXE==0],"ans",TRUE)(gH) version 4.87 fonctions d'aides : lit() fqt() fql() ic() fapprox() chi2() fcomp() datagh() taper aide() pour revoir cette liste DESCRIPTION STATISTIQUE DE LA VARIABLE AGE (global) Taille 99 individus Moyenne 35.8283 ans Ecart-type 17.5529 ans Coef. de variation 49 % 1er Quartile 22.0000 ans Mediane 29.0000 ans 3eme Quartile 48.5000 ans iqr absolu 26.5000 ans iqr relatif 91.0000 % Minimum 11 ans Maximum 78 ans Tracé tige et feuilles The decimal point is 1 digit(s) to the right of the | 1 | 1223455567778889999 2 | 0011222344555566667778888888999 3 | 0011123556799 4 | 0123446778899 5 | 0002239 6 | 0011222345 7 | 033368 DESCRIPTION STATISTIQUE DE LA VARIABLE AGE (femmes) Taille 64 individus Moyenne 35.5156 ans Ecart-type 18.0273 ans Coef. de variation 50 % 1er Quartile 21.0000 ans Mediane 29.5000 ans 3eme Quartile 48.2500 ans iqr absolu 27.2500 ans iqr relatif 92.0000 % Minimum 11 ans Maximum 76 ans Tracé tige et feuilles The decimal point is 1 digit(s) to the right of the | 1 | 123455577889999 2 | 11234455566778889 3 | 011156799 4 | 01446889 5 | 00039 6 | 01123 7 | 03336 DESCRIPTION STATISTIQUE DE LA VARIABLE AGE (hommes) Taille 35 individus Moyenne 36.4000 ans Ecart-type 16.8927 ans Coef. de variation 46 % 1er Quartile 25.5000 ans Mediane 29.0000 ans 3eme Quartile 48.0000 ans iqr absolu 22.5000 ans iqr relatif 78.0000 % Minimum 12 ans Maximum 78 ans Tracé tige et feuilles The decimal point is 1 digit(s) to the right of the | 1 | 2678 2 | 00225667888899 3 | 0235 4 | 23779 5 | 22 6 | 02245 7 | 8Le peu de données dans le groupe des hommes et la non normalité apparente des données mène à la conclusion qu'il n'est pas conseillé d'utiliser une analyse paramétrique. On utilise donc un test de Wilcox et on trace des boites à moustaches.
# test non-paramétrique print( wilcox.test(AGE~SEXE,data=elf) ) # boites à moustaches associées with(elf, boxplot(AGE~factor(SEXE,levels=0:1,labels=c("Hommes","Femmes")) ,col="yellow",main="AGES (dossier ELF)",ylim=c(10,80)) ) # pour information with(elf, boxplot(AGE~factor(SEXE,levels=0:1,labels=c("Hommes","Femmes")), col="yellow",main="AGES (dossier ELF)",ylim=c(10,80), plot=FALSE ) # fin de boxplot ) # fin de withWilcoxon rank sum test with continuity correction data: AGE by SEXE W = 1195, p-value = 0.5854 alternative hypothesis: true location shift is not equal to 0 $stats [,1] [,2] [1,] 12.0 11.0 [2,] 25.5 21.0 [3,] 29.0 29.5 [4,] 48.0 48.5 [5,] 78.0 76.0 attr(,"class") Hommes "integer" $n [1] 35 64 $conf [,1] [,2] [1,] 22.99095 24.06875 [2,] 35.00905 34.93125 $out numeric(0) $group numeric(0) $names [1] "Hommes" "Femmes"La conclusion est qu'il n'y a pas de différence significative détectée au seuil de 5 % entre les hommes et les femmes pour la variable AGE dans ce petit échantillon d'une centaine de personnes dont deux tiers de femmes et un tiers d'hommes.
Code-source php de cette page. Retour à la page principale du cours.
Retour à la page principale de (gH)