Analyse de la variance, et participation en TD

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 https://latex.codecogs.com/gif.latex?X, 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’).

> download.file("http://freakonometrics.free.fr/notes.RData",destfile= tempfile())
> load("notes.RData")
> tail(base)
PART GENRE ORIGINE  NOTE
112 volontaire     F      R1 16.50
113   non_vol.     M      R1 11.50
114   non_vol.     F      R1 10.25
115   non_vol.     F      R1 10.75
116   non_vol.     F   autre 10.50
117 volontaire     M      R1 15.75

Avant 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)

> boxplot(base$NOTE~base$PART,ylim=c(6,20))
> abline(h=mean(base$NOTE),lty=2,col="red")
> 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

> boxplot(base$NOTE~base$GENRE,ylim=c(6,20))
> abline(h=mean(base$NOTE),lty=2,col="red")
> lines(1:2,with(base,aggregate(NOTE, by=list(GENRE),mean)$x),)

ou en fonction de l’origine des étudiants

> boxplot(base$NOTE~base$ORIGINE,ylim=c(6,20))
> abline(h=mean(base$NOTE),lty=2,col="red")
> lines(1:2,with(base,aggregate(NOTE, by=list(ORIGINE),mean)$x))

Dans l’anova, on suppose que https://latex.codecogs.com/gif.latex?X_{i,j}=\mu_j+\varepsilon_{i,j}, ou mieux encore

https://latex.codecogs.com/gif.latex?X_{i,j}=\mu+\alpha_j+\varepsilon_{i,j}où les https://latex.codecogs.com/gif.latex?j désigne les traitements possibles (ici, il y en a 3).

On va considérer le test de https://latex.codecogs.com/gif.latex?H_0:\forall%20j,\alpha_j=0, avec comme hypothèse alternative le complémentaire https://latex.codecogs.com/gif.latex?H_1:\exists%20j,\alpha_j\neq0. 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)

> sum(residuals(lm(NOTE~1,data=base))^2)
[1] 947.4979

correspondant à

> (SCR0=sum((base$NOTE-mean(base$NOTE))^2))
[1] 947.4979

Le second, on régresse sur les traitements,

> sum(residuals(lm(NOTE~PART,data=base))^2)
[1] 112.5032

ce qui revient à regarder les erreurs quand on compare à la moyenne des sous-groupes,

> Ybar=with(base,aggregate(NOTE,by=list(PART), mean)$x)
> base$moyNOTE=Ybar[as.numeric(base$PART)]
> (SCR1=sum((base$NOTE-base$moyNOTE)^2))
[1] 112.5032

La statistique de Fisher est ici

> (F=(SCR0-SCR1)*(nrow(base)-3)/SCR1/(3-1))
[1] 423.0518

Pour savoir si on est dans la région d’acceptation, ou de rejet de l’hypothèse https://latex.codecogs.com/gif.latex?H_0, on peut regarder la valeur critique, correspondant au quantile à 95% de la loi de Fisher,

> qf(.95,3-1,nrow(base)-3)
[1] 3.075853

Comme on est très largement dessus de ce seuil, on va rejeter https://latex.codecogs.com/gif.latex?H_0. On peut aussi calculer la p-value,

> 1-pf(F,3-1,nrow(base)-3)
[1] 0

qui est ici (numériquement) nul. Ca correspond à ce qu’on obtient avec la fonction

> anova(lm(NOTE~PART,data=base))
Analysis of Variance Table

Response: NOTE
Df Sum Sq Mean Sq F value    Pr(>F)
PART        2 834.99  417.50  423.05 < 2.2e-16 ***
Residuals 114 112.50    0.99
---

ou encore

> aov(NOTE~PART,data=base)
Call:
aov(formula = NOTE ~ PART, data = base)

Terms:
PART Residuals
Sum of Squares  834.9946  112.5032
Deg. of Freedom        2       114

Residual standard error: 0.9934135
Estimated effects may be unbalanced

que l’ont peut résumer avec

> anova(aov(NOTE~PART,data=base))
Analysis of Variance Table

Response: NOTE
Df    Sum Sq Mean Sq F value Pr(>F)
PART        2 834.99  417.50  423.05 < 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,

> TukeyHSD(aov(NOTE~PART,data=base))
Tukey multiple comparisons of means
95% family-wise confidence level

Fit: aov(formula = NOTE ~ PART, data = base)

$PART
diff       lwr      upr    p adj
non_vol.-non_part.   0.60416 -0.04784 1.2561 0.07539
volontaire-non_part. 6.66379  5.92912 7.3984 0.00000
volontaire-non_vol.  6.05962  5.54078 6.5784 0.00000

On 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,

> with(base,pairwise.t.test(NOTE,PART))

Pairwise comparisons using t tests with pooled SD

data:  NOTE and PART

non_part. non_vol.
non_vol.   0.03      -
volontaire <2e-16    <2e-16

Si 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,

> with(base,t.test(NOTE[PART=='volontaire'], NOTE[PART!='volontaire']))

Welch Two Sample t-test

data:  NOTE[PART == "volontaire"] and NOTE[PART != "volontaire"]
t = 29.511, df = 50.73, p-value < 2.2e-16
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
5.749719 6.589231
sample estimates:
mean of x mean of y
16.66379  10.49432

On 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èle

https://latex.codecogs.com/gif.latex?Y_i%20=%20\mu+\alpha_j+\varepsilon_{i,j}avec https://latex.codecogs.com/gif.latex?\varepsilon_{i,j}\sim\mathcal{N}(0,\sigma^2), correspond à un modèle homoscédastique, au sens où https://latex.codecogs.com/gif.latex?\sigma^2 ne dépend pas du traitement https://latex.codecogs.com/gif.latex?j. 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%

> bartlett.test(base$NOTE,base$PART)

Bartlett test of homogeneity of variances

data:  base$NOTE and base$PART
Bartlett's K-squared = 0.5524, df = 2, p-value = 0.7587

On 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 https://latex.codecogs.com/gif.latex?A+B, et on peut inclure des interactions on considérant https://latex.codecogs.com/gif.latex?A+B+A\times%20B

De manière générale, on a ici un modèle de la forme

https://latex.codecogs.com/gif.latex?X_{i,j,k}=\mu+\alpha_j+\beta_k+\gamma_{j,k}+\varepsilon_{i,j,k}

https://latex.codecogs.com/gif.latex?A=(\alpha_j) désigne les écarts à la moyenne des moyennes des grandeurs observées, par traitement, et où https://latex.codecogs.com/gif.latex?B=(\beta_k) 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 ici

https://latex.codecogs.com/gif.latex?%20\widehat{\mu}=\overline{x}

correspondant à la moyenne globale

https://latex.codecogs.com/gif.latex?%20\widehat{\alpha}_j=\overline{x}_j-\overline{x}

correspondand à la moyenne pa r traitement (ou plutôt son écart par rapport à la moyenne globale),

https://latex.codecogs.com/gif.latex?\widehat{\beta}_k=\overline{x}_k-\overline{x}

et enfin

https://latex.codecogs.com/gif.latex?\widehat{\gamma}_{j,k}=\overline{x}_{j,k}-\widehat{\mu}-\widehat{\alpha}_j-\widehat{\beta}_k

soit

https://latex.codecogs.com/gif.latex?\widehat{\gamma}_{j,k}=\overline{x}_{j,k}-\overline{x}_j-\overline{x}_k+\overline{x}

où on utilise abusivement https://latex.codecogs.com/gif.latex?\overline{x}_{\cdot} 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 https://latex.codecogs.com/gif.latex?m individus par traitement et par groupe)

https://latex.codecogs.com/gif.latex?SCR_0=\sum_{i,j,k}[x_{i,j,k}-\overline{x}]^2

pour le modèle contraint,

https://latex.codecogs.com/gif.latex?SCR_{\alpha}=[m%20J]\sum_{j}[\overline{x}_{j}-\overline{x}]^2https://latex.codecogs.com/gif.latex?SCR_{\beta}=[m%20K]\sum_{k}[\overline{x}_{k}-\overline{x}]^2https://latex.codecogs.com/gif.latex?SCR_{\gamma}=m\sum_{j,k}[\overline{x}_{j,k}-\overline{x}_j-\overline{x}_k+\overline{x}]^2

https://latex.codecogs.com/gif.latex?J et https://latex.codecogs.com/gif.latex?K désignent le nombre de traitement et le nombre de groupes, et où

https://latex.codecogs.com/gif.latex?SCR_{\varepsilon}=\sum_{i,j,k}\widehat{\varepsilon}_{i,j,k}^2

La formule de décomposition de la variance donne ici

https://latex.codecogs.com/gif.latex?SCR_0=SCR_{\alpha}+SCR_{\beta}+SCR_{\gamma}+SCR_{\varepsilon}

On va passer sur les calculs manuels,

> aov(NOTE~PART*GENRE,data=base)
Call:
aov(formula = NOTE ~ PART * GENRE, data = base)

Terms:
PART    GENRE PART:GENRE Residuals
Sum of Squares  834.9946  20.9618     3.4398   88.1017
Deg. of Freedom        2        1          2       111

Residual standard error: 0.8909034
Estimated effects may be unbalanced

que l’on peut résumer avec la fonction suivante

> anova(aov(NOTE~PART*GENRE,data=base))
Analysis of Variance Table

Response: NOTE
Df Sum Sq Mean Sq  F value    Pr(>F)
PART         2 834.99  417.50 526.0081 < 2.2e-16 ***
GENRE        1  20.96   20.96  26.4099 1.194e-06 ***
PART:GENRE   2   3.44    1.72   2.1669    0.1194
Residuals  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

> anova(aov(NOTE~PART+GENRE,data=base))
Analysis of Variance Table

Response: NOTE
Df Sum Sq Mean Sq F value    Pr(>F)
PART        2 834.99  417.50 515.364 < 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…)


One thought on “Analyse de la variance, et participation en TD”

Leave a Reply

Your email address will not be published. Required fields are marked *