{"word_count": 1184, "character_count": 9270}
Nous avons vu lors du dernier cours l'analyse de la variance. Comme je l'avais expliqué, comme son nom ne l'indique pas, le but de l'anova est de comparer des moyennes. Histoire d'illustrer ce test, considerons la base suivante, qui donne la note d'étudiants à un examen final de statistique (note comprise entre 0 et 20), dans un cours de Master. Ca sera notre variable , comme nous l'avions noté en cours. La variable de "traitement" sera la participation des étudiants en TD, avec comme modalité possible 'participation volontaire' (pour aller faire des exercices au tableau), 'participation non-volontaire' (les étudiants vont au tableau quand ils y sont appelés) et enfin 'non-participation' (les étudiants qui ne participent pas, ou qui refusent de passer au tableau s'ils y sont appelés). Afin de constituer des groupes, nous disposons de deux variables. La première est le genre de l'étudiant ('F' et 'M') et la seconde est sur l'origine des étudiants (suivant qu'ils ont fait leur licence ou pas, à Rennes, 'R1' ou 'autre').&#62; download.file("http://freakonometrics.free.fr/notes.RData",destfile= tempfile())&#62; load("notes.RData")&#62; tail(base)PART GENRE ORIGINE  NOTE112 volontaire     F      R1 16.50113   non_vol.     M      R1 11.50114   non_vol.     F      R1 10.25115   non_vol.     F      R1 10.75116   non_vol.     F   autre 10.50117 volontaire     M      R1 15.75Avant de commencer avec une analyse à plusieurs facteurs, commençons par une analyse à un facteur. On peut regarder l'évolution des notes, en fonction de la variable de traitement (la participation en TD)&#62; boxplot(base$NOTE~base$PART,ylim=c(6,20))&#62; abline(h=mean(base$NOTE),lty=2,col="red")&#62; lines(1:3,with(base,aggregate(NOTE,by=list(PART), mean)$x),pch=19,lwd=3,type="b")On peut aussi regarder la note, en fonction du genre&#62; boxplot(base$NOTE~base$GENRE,ylim=c(6,20))&#62; abline(h=mean(base$NOTE),lty=2,col="red")&#62; lines(1:2,with(base,aggregate(NOTE, by=list(GENRE),mean)$x),)ou en fonction de l'origine des étudiants&#62; boxplot(base$NOTE~base$ORIGINE,ylim=c(6,20))&#62; abline(h=mean(base$NOTE),lty=2,col="red")&#62; lines(1:2,with(base,aggregate(NOTE, by=list(ORIGINE),mean)$x))Dans l'anova, on suppose que , ou mieux encoreoù les  désigne les traitements possibles (ici, il y en a 3).On va considérer le test de , avec comme hypothèse alternative le complémentaire . On va alors estimer deux modèles. Le premier est un modèle contraint, au sens où on ne tient pas compte du traitement (et on régresse juste sur la constante)&#62; sum(residuals(lm(NOTE~1,data=base))^2)[1] 947.4979correspondant à&#62; (SCR0=sum((base$NOTE-mean(base$NOTE))^2))[1] 947.4979Le second, on régresse sur les traitements,&#62; sum(residuals(lm(NOTE~PART,data=base))^2)[1] 112.5032ce qui revient à regarder les erreurs quand on compare à la moyenne des sous-groupes,&#62; Ybar=with(base,aggregate(NOTE,by=list(PART), mean)$x)&#62; base$moyNOTE=Ybar[as.numeric(base$PART)]&#62; (SCR1=sum((base$NOTE-base$moyNOTE)^2))[1] 112.5032La statistique de Fisher est ici&#62; (F=(SCR0-SCR1)*(nrow(base)-3)/SCR1/(3-1))[1] 423.0518Pour savoir si on est dans la région d'acceptation, ou de rejet de l'hypothèse , on peut regarder la valeur critique, correspondant au quantile à 95% de la loi de Fisher,&#62; qf(.95,3-1,nrow(base)-3)[1] 3.075853Comme on est très largement dessus de ce seuil, on va rejeter . On peut aussi calculer la p-value,&#62; 1-pf(F,3-1,nrow(base)-3)[1] 0qui est ici (numériquement) nul. Ca correspond à ce qu'on obtient avec la fonction&#62; anova(lm(NOTE~PART,data=base))Analysis of Variance TableResponse: NOTEDf Sum Sq Mean Sq F value    Pr(&#62;F)PART        2 834.99  417.50  423.05 &#60; 2.2e-16 ***Residuals 114 112.50    0.99---ou encore&#62; aov(NOTE~PART,data=base)Call:aov(formula = NOTE ~ PART, data = base)Terms:PART ResidualsSum of Squares  834.9946  112.5032Deg. of Freedom        2       114Residual standard error: 0.9934135Estimated effects may be unbalancedque l'ont peut résumer avec&#62; anova(aov(NOTE~PART,data=base))Analysis of Variance TableResponse: NOTEDf    Sum Sq Mean Sq F value Pr(&#62;F)PART        2 834.99  417.50  423.05 &#60; 2.2e-16 ***Residuals 114 112.50    0.99---On a ici toutes les grandeurs calculées auparavant.On voit ici que la note n'est pas indépendante du traitement. Mais on n'en sait pas beaucoup plus...On peut creuser davantage. Le test de Tukey propose un "test multiple", qui va regarder les différences de moyennes, par paires,&#62; TukeyHSD(aov(NOTE~PART,data=base))Tukey multiple comparisons of means95% family-wise confidence levelFit: aov(formula = NOTE ~ PART, data = base)$PARTdiff       lwr      upr    p adjnon_vol.-non_part.   0.60416 -0.04784 1.2561 0.07539volontaire-non_part. 6.66379  5.92912 7.3984 0.00000volontaire-non_vol.  6.05962  5.54078 6.5784 0.00000On voit ici que la différence entre entre 'non-volontaire' et 'non-participatif' est non-significativement non-nulle. Ou dit plus simplement (mais moins juste), disons on va accepter hypothèse qu'elle est nulle. En revanche, la note quand on participe 'volontairement' est significativement plus grande que quand on participe 'non-volontairement' ou qu'on ne participe pas.On peut aussi regarder les tests de Student, par paires,&#62; with(base,pairwise.t.test(NOTE,PART))Pairwise comparisons using t tests with pooled SDdata:  NOTE and PARTnon_part. non_vol.non_vol.   0.03      -volontaire &#60;2e-16    &#60;2e-16Si on regroupe les deux modalités du traitement 'non-volontaire' et 'non-participatif', et qu'on compare cette modalité, avec 'volontaire', on se retrouve à faire un test de moyennes,&#62; with(base,t.test(NOTE[PART=='volontaire'], NOTE[PART!='volontaire']))Welch Two Sample t-testdata:  NOTE[PART == "volontaire"] and NOTE[PART != "volontaire"]t = 29.511, df = 50.73, p-value &#60; 2.2e-16alternative hypothesis: true difference in means is not equal to 095 percent confidence interval:5.749719 6.589231sample estimates:mean of x mean of y16.66379  10.49432On voit que l'on accepte ici l'hypothèse que les étudiants 'volontaires' ont une note différente des autres.Avant d'aller plus loin, rappelons que dans le modèleavec , correspond à un modèle homoscédastique, au sens où  ne dépend pas du traitement . On peut tester cette hypothèse, avec un test de Bartlett (qui va tester l'homogénéité des variances), en rappelant que l'hypothèse 'variance homogène' est validée si la p-value dépasse 5%&#62; bartlett.test(base$NOTE,base$PART)Bartlett test of homogeneity of variancesdata:  base$NOTE and base$PARTBartlett's K-squared = 0.5524, df = 2, p-value = 0.7587On est bien sous les hypothèses du modèle de Fisher.Pour aller plus loin, on peut tenter une analyse de variance à deux facteurs, avec toujours notre traitement, mais aussi une variable de groupe, basé sur le genre. On parlera d'un modèle sans interraction quand la moyenne est de la forme , et on peut inclure des interactions on considérant De manière générale, on a ici un modèle de la formeoù  désigne les écarts à la moyenne des moyennes des grandeurs observées, par traitement, et où  désigne les écarts à la moyenne des moyennes des grandeurs observées, par groupe. Pour que le modèle est identifiable en rajoutant quelques contraintes (les valeurs sont centrées). Les estimateurs du maximum de vraisemblance sont icicorrespondant à la moyenne globalecorrespondand à la moyenne pa r traitement (ou plutôt son écart par rapport à la moyenne globale),et enfinsoitoù on utilise abusivement  pour parler des moyennes par traitement, ou par groupe, avec des notations (probablement) naturelles, pour faire notre analyse de la variance, avec un traitement, et un groupe, on utilise lentes grandeurs suivante (on suppose ici qu'il y a  individus par traitement et par groupe)pour le modèle contraint,où  et  désignent le nombre de traitement et le nombre de groupes, et oùLa formule de décomposition de la variance donne iciOn va passer sur les calculs manuels,&#62; aov(NOTE~PART*GENRE,data=base)Call:aov(formula = NOTE ~ PART * GENRE, data = base)Terms:PART    GENRE PART:GENRE ResidualsSum of Squares  834.9946  20.9618     3.4398   88.1017Deg. of Freedom        2        1          2       111Residual standard error: 0.8909034Estimated effects may be unbalancedque l'on peut résumer avec la fonction suivante&#62; anova(aov(NOTE~PART*GENRE,data=base))Analysis of Variance TableResponse: NOTEDf Sum Sq Mean Sq  F value    Pr(&#62;F)PART         2 834.99  417.50 526.0081 &#60; 2.2e-16 ***GENRE        1  20.96   20.96  26.4099 1.194e-06 ***PART:GENRE   2   3.44    1.72   2.1669    0.1194Residuals  111  88.10    0.79---Comme il ne semble pas y avoir d'interractions entre le traitement et le groupe, on peut le supprimer de l'analyse de la variance&#62; anova(aov(NOTE~PART+GENRE,data=base))Analysis of Variance TableResponse: NOTEDf Sum Sq Mean Sq F value    Pr(&#62;F)PART        2 834.99  417.50 515.364 &#60; 2.2e-16 ***GENRE       1  20.96   20.96  25.875 1.461e-06 ***Residuals 113  91.54    0.81---Il y a un effet du genre ici. Mais moins fort que celui de la participartion. Moralité ? Les statistiques, ça sert vraiment à quelque chose, et mieux vaut passer au tableau (volontairement), ça aide ! (même si je n'ai montré aucun lien de causalité ici...)