Stata : analyse de la variance

Table des matières

L’analyse de variance (ANOVA), souvent retrouvée dans l’exploitation d’un plan d’expérience, peut être dérivée dans la plupart des cas à partir des modèles de régression vus au chapitre 3 (et plus tard, aux chapitres 9 et 10). Toutefois, une présentation de la décomposition des sources de variance et du calcul des sommes de carré associées, indépendante de l’approche de régression linéaire, est proposée dans ce chapitre, d’autant que Stata dispose de commandes spécifiques pour les modèles d’ANOVA, en particulier oneway et anova.

Dupont [2] consacre plusieurs chapitres de son excellent ouvrage à l’analyse de variance, incluant le cas des mesures répétées.

Analyse de variance à un facteur

Modèle d’ANOVA

Dans l’exemple suivant, tiré de [2], on s’intéresse à la relation entre l’âge des participants à l’étude au moment du diagostique et un polymorphisme de séquence. Les données sont disponibles dans le fichier polymorphism.dta :

use data/polymorphism
list in 1/5
set more off
use data/polymorphism
list in 1/5

     +---------------------+
     | id   age   genotype |
     |---------------------|
  1. |  1    43    1.6/1.6 |
  2. |  2    47    1.6/1.6 |
  3. |  3    55    1.6/1.6 |
  4. |  4    57    1.6/1.6 |
  5. |  5    61    1.6/1.6 |
     +---------------------+

Un résumé descriptif des données de groupe (ici défini par le type de polymorphisme) peut être obtenu très rapidement à l’aide de tabstat :

tabstat age, by(genotype) stat(mean sd n)
tabstat age, by(genotype) stat(mean sd n)

Summary for variables: age
     by categories of: genotype (Genotype)

genotype |      mean        sd         N
---------+------------------------------
 1.6/1.6 |  64.64286  11.18108        14
 1.6/0.7 |  64.37931  13.25954        29
 0.7/0.7 |    50.375  10.63877        16
---------+------------------------------
   Total |  60.64407  13.49427        59
----------------------------------------

Un résumé graphique de la distribution conditionnelle de l’âge dans les trois groupes s’obtient tout aussi simplement à l’aide de graph box :

set scheme plotplain
graph box age, over(genotype)
graph export "fig-04-boxplot-age-genotype.eps", replace

fig-04-boxplot-age-genotype.png

Figure 1 : Distribution de l’âge de diagnostic en fonction du polymorphisme

Il est également possible d’afficher ces distributions sous forme d’histogrammes d’effectifs ou de fréquences, par exemple :

histogram age, by(genotype, col(3)) freq
graph export "fig-04-hist-age-genotype.eps", replace

fig-04-hist-age-genotype.png

Figure 2 : Distribution de l’âge de diagnostic en fonction du polymorphisme

Le tableau d’ANOVA peut être constuit à l’aide de oneway dans le cas où il n’y a qu’une seule avriable explicative. L’option tabulate fournit exactement le même résumé numérique que la commande tabstat ci-dessus, tandis que l’option bonferonni ajoute la comparaison de l’ensemble des paires de moyennes avec une correction de Bonferroni pour les tests multiples. Si l’on ne souhiate que ces derniers, il suffit d’ajouter l’option noanova, mais en règle générale il est préférable de réaliser le test d’ensemble avant d’examiner les tests spécifiques, sauf si l’on a une hypothèse de recherche très particulière.

oneway age genotype
oneway age genotype

                        Analysis of Variance
    Source              SS         df      MS            F     Prob > F
------------------------------------------------------------------------
Between groups      2315.73355      2   1157.86678      7.86     0.0010
 Within groups      8245.79187     56   147.246283
------------------------------------------------------------------------
    Total           10561.5254     58   182.095266

Bartlett's test for equal variances:  chi2(2) =   1.0798  Prob>chi2 = 0.583

Conditions d’application du test

Il est généralement recommendé de vérifier les hypothèses du modèle, en particulier celle portant sur l’égalité des variances parentes, de manière graphique. Notons que oneway propose d’office un test formel d’égalité des variances (test de Bartlett). Le test de Levenne est disponible via la commande robvar (statistique W50) :

robvar age, by(genotype)
robvar age, by(genotype)

            |     Summary of Age at Diagnosis
   Genotype |        Mean   Std. Dev.       Freq.
------------+------------------------------------
    1.6/1.6 |   64.642857   11.181077          14
    1.6/0.7 |    64.37931   13.259535          29
    0.7/0.7 |      50.375   10.638766          16
------------+------------------------------------
      Total |   60.644068   13.494268          59

W0  =  0.83032671   df(2, 56)     Pr > F = 0.44120161

W50 =  0.60460508   df(2, 56)     Pr > F = 0.54981692

W10 =  0.79381598   df(2, 56)     Pr > F = 0.45713722

Concernant la normalité des distributions (parentes), à partir de l’échantillon il est toujours possible de tracer la fonction de répartition des observations, un histogramme ou un diagramme de type boîte à moustaches (cf. section suivante). La commande swilk fournit le test de Shapiro-Wilks, mais comme les méthodes graphiques restent à privilégier, voici une façon de vérifier la normalité des résidus du modèle (dans le cas de l’ANOVA à un facteur, cela ne change pas grand-chose par rapport au fait de travailler sur les observations directement, mais dans le cas à plus d’un facteur en interaction c’est la méthode recommendée) :

quietly anova age genotype
predict r, resid
qnorm r
graph export "fig-04-qnorm-age-genotype-resid.eps", replace

fig-04-qnorm-age-genotype-resid.png

Figure 3 : Distribution des résidus du modèle d’ANOVA

Commandes graphiques

Stata n’offre pas vraiment de commandes graphiques de post-estimation à l’exception de celles applicables au cas du modèle linéaire, telles que les graphiques de résidus en fonction des valeurs prédites ou observées. Il existe toutefois de nombreux packages disponibles sur le site SSC, en particulier la commande anovaplot (une approche manuelle pour construire un graphique d’interaction est présentée dans la section suivante).

Une référence intéressante concernant les approches graphiques pour le diagnostique des modèles de régression, au sens large, est l’article de Nick Cox [1].

Comparaisons multiples

Lorsque le résultat du test d’ANOVA est significatif, il est souvent intéressant, voire nécessaire d’aller inspecter le résultat des comparaisons de chaque paire de moyennes. Comme indiqué plus haut, il est déconseillé de procéder à de telles comparaisons multiples lorsque l’ANOVA n’est en elle-même pas significative puisque dans ce cas il s’agit un peu de ce que l’on appelle du data fishing. Cela dit, dans certains cas de figure, il est possible de pré-spécifier les comparaisons d’intérêt et de n’analyser que ces contrastes d’intérêt. Dans ce cas, le modèle d’ANOVA est même facultatif.

Il y a deux manières de procéder pour réaliser des comparaisons de paires de moyennes sous Stata (outre la comparaison manuelle à l’aide de ttest): (1) la commande oneway offre une option permettant de tester les paires de moyennes et d’appliquer une correction pour borner le risque d’erreur d’ensemble (FWER) ; (2) les commandes pwcompare, contrast ou lincome (ou test). La commande margins fournit plus de souplesse dans l’estimation conditionnelle ou marginale des effets et permet en outre la représentation graphique des résultats grâce à marginsplot. En dehors de pwcompare, les autres commandes restent applicables dans le cas du modèle linéaire dans son ensemble.

Voici le résultat des tests multiples à partir de oneway :

oneway age genotype, bonferroni noanova
oneway age genotype, bonferroni noanova

                  Comparison of Age at Diagnosis by Genotype
                                (Bonferroni)
Row Mean-|
Col Mean |    1.6/1.6    1.6/0.7
---------+----------------------
 1.6/0.7 |   -.263547
         |      1.000
         |
 0.7/0.7 |   -14.2679   -14.0043
         |      0.007      0.001

On arriverait naturellement aux mêmes conclusions en construisant manuellement les tests de Student correspondants, moyennant la prise en compte de la bonne somme de carrés (variance « poolée ») :

quietly ttest age if genotype != 1, by(genotype)
display r(p)*3
quietly ttest age if genotype != 1, by(genotype)
00228566

Analyse de variance à plusieurs facteurs

La commande oneway est limité au cas à un facteur explicatif. La commande anova est plus générale et couvre : les plans factoriels et emboîtés, les plans équilibrés ou non (cf. calcul des sommes de carrés), les mesures répétées, l’analyse de covariance. Dans le cas à un facteur à effet fixe, on retrouvera évidemment les mêmes résultats que plus haut :

anova age genotype
anova age genotype

                           Number of obs =      59     R-squared     =  0.2193
                           Root MSE      = 12.1345     Adj R-squared =  0.1914

                  Source |  Partial SS    df       MS           F     Prob > F
              -----------+----------------------------------------------------
                   Model |  2315.73355     2  1157.86678       7.86     0.0010
                         |
                genotype |  2315.73355     2  1157.86678       7.86     0.0010
                         |
                Residual |  8245.79187    56  147.246283   
              -----------+----------------------------------------------------
                   Total |  10561.5254    58  182.095266

Les comparaisons par paires de moyennes s’obtiennent à l’aide de pwcompare, commande plus générale que pwmean. Les options de correction (mcompare()) incluent en plus : tukey, snk, duncan et dunnett.

pwcompare genotype, cformat(%3.2f)
pwcompare genotype, cformat(%3.2f)

Pairwise comparisons of marginal linear predictions

Margins      : asbalanced

---------------------------------------------------------------------
                    |                                 Unadjusted
                    |   Contrast   Std. Err.     [95% Conf. Interval]
--------------------+------------------------------------------------
           genotype |
1.6/0.7 vs 1.6/1.6  |      -0.26       3.95         -8.17        7.65
0.7/0.7 vs 1.6/1.6  |     -14.27       4.44        -23.16       -5.37
0.7/0.7 vs 1.6/0.7  |     -14.00       3.78        -21.57       -6.43
---------------------------------------------------------------------

Voici un exemple de plan d’expérience dans lequel on s’intéresse à la fabrication d’une batterie capable de fonctionner dans des conditions extrêmes de température [3]. Cette étude comprend deux facteurs expérimentaux ayant trois niveaux chacun : la température (°F) et un paramètre lié au design de la batterie elle-même. Il s’agit donc d’un plan factoriel \(3^2\). Les données sont disponibles dans le fichier battery.txt.

import delimited "data/battery.txt", delimiter("", collapse) varnames(1) clear
list in 1/3
 r
(3 vars, 36 obs)
list in 1/3

     +----------------------------+
     | life   material   temper~e |
     |----------------------------|
  1. |  130          1         15 |
  2. |   74          1         15 |
  3. |  155          1         15 |
     +----------------------------+

Voici les résultats pour le modèle avec interaction :

anova life material##temperature
anova life material##temperature

                           Number of obs =      36     R-squared     =  0.7652
                           Root MSE      = 25.9849     Adj R-squared =  0.6956

                  Source |  Partial SS    df       MS           F     Prob > F
    ---------------------+----------------------------------------------------
                   Model |  59416.2222     8  7427.02778      11.00     0.0000
                         |
                material |  10683.7222     2  5341.86111       7.91     0.0020
             temperature |  39118.7222     2  19559.3611      28.97     0.0000
    material#temperature |  9613.77778     4  2403.44444       3.56     0.0186
                         |
                Residual |    18230.75    27  675.212963   
    ---------------------+----------------------------------------------------
                   Total |  77646.9722    35  2218.48492

Un graphique d’interaction peut être constuit à l’aide de scatter comme suit :

preserve
collapse (mean) mean=life (sd) sd=life, by(material temperature)
list in 1/3
drop sd
reshape wide mean, i(temperature) j(material)
twoway connected mean* temperature, legend(label(1 "#1") label(2 "#2") label(3 "#3")) ytitle(Mean life)
graph export "fig-04-scatter-life-battery.eps", replace
restore

fig-04-scatter-life-battery.png

Figure 4 : Distribution de l’âge de diagnostic en fonction du polymorphisme

La commande marginplot simplifierait grandement la tâche. Il existe également les commandes rcap et serrbar pour gérer les barres d’erreur.

Voyons ce que l’on peut réaliser en termes de tests de contrastes à l’aide de pwcompare. Dans un premier temps, la comparaison de l’ensemble des paires de moyennes peut être réalisée en décomposant le terme d’interaction, et en utilisant la méthode de Tukey comme correction pour les tests multiples :

pwcompare material#temperature, mcompare(tukey)
pwcompare material#temperature, mcompare(tukey)

Pairwise comparisons of marginal linear predictions

Margins      : asbalanced

-----------------------------------
                     |    Number of
                     |  Comparisons
---------------------+-------------
material#temperature |           36
-----------------------------------

----------------------------------------------------------------------
                     |                                   Tukey
                     |   Contrast   Std. Err.     [95% Conf. Interval]
---------------------+------------------------------------------------
material#temperature |
 (1  70) vs (1  15)  |      -77.5   18.37407     -139.3232   -15.67681
 (1 125) vs (1  15)  |     -77.25   18.37407     -139.0732   -15.42681
 (2  15) vs (1  15)  |         21   18.37407     -40.82319    82.82319
 (2  70) vs (1  15)  |        -15   18.37407     -76.82319    46.82319
 (2 125) vs (1  15)  |     -85.25   18.37407     -147.0732   -23.42681
 (3  15) vs (1  15)  |       9.25   18.37407     -52.57319    71.07319
 (3  70) vs (1  15)  |         11   18.37407     -50.82319    72.82319
 (3 125) vs (1  15)  |     -49.25   18.37407     -111.0732    12.57319
 (1 125) vs (1  70)  |        .25   18.37407     -61.57319    62.07319
 (2  15) vs (1  70)  |       98.5   18.37407      36.67681    160.3232
 (2  70) vs (1  70)  |       62.5   18.37407      .6768132    124.3232
 (2 125) vs (1  70)  |      -7.75   18.37407     -69.57319    54.07319
 (3  15) vs (1  70)  |      86.75   18.37407      24.92681    148.5732
 (3  70) vs (1  70)  |       88.5   18.37407      26.67681    150.3232
 (3 125) vs (1  70)  |      28.25   18.37407     -33.57319    90.07319
 (2  15) vs (1 125)  |      98.25   18.37407      36.42681    160.0732
 (2  70) vs (1 125)  |      62.25   18.37407      .4268132    124.0732
 (2 125) vs (1 125)  |         -8   18.37407     -69.82319    53.82319
 (3  15) vs (1 125)  |       86.5   18.37407      24.67681    148.3232
 (3  70) vs (1 125)  |      88.25   18.37407      26.42681    150.0732
 (3 125) vs (1 125)  |         28   18.37407     -33.82319    89.82319
 (2  70) vs (2  15)  |        -36   18.37407     -97.82319    25.82319
 (2 125) vs (2  15)  |    -106.25   18.37407     -168.0732   -44.42681
 (3  15) vs (2  15)  |     -11.75   18.37407     -73.57319    50.07319
 (3  70) vs (2  15)  |        -10   18.37407     -71.82319    51.82319
 (3 125) vs (2  15)  |     -70.25   18.37407     -132.0732   -8.426813
 (2 125) vs (2  70)  |     -70.25   18.37407     -132.0732   -8.426813
 (3  15) vs (2  70)  |      24.25   18.37407     -37.57319    86.07319
 (3  70) vs (2  70)  |         26   18.37407     -35.82319    87.82319
 (3 125) vs (2  70)  |     -34.25   18.37407     -96.07319    27.57319
 (3  15) vs (2 125)  |       94.5   18.37407      32.67681    156.3232
 (3  70) vs (2 125)  |      96.25   18.37407      34.42681    158.0732
 (3 125) vs (2 125)  |         36   18.37407     -25.82319    97.82319
 (3  70) vs (3  15)  |       1.75   18.37407     -60.07319    63.57319
 (3 125) vs (3  15)  |      -58.5   18.37407     -120.3232    3.323187
 (3 125) vs (3  70)  |     -60.25   18.37407     -122.0732    1.573187
----------------------------------------------------------------------

Le test d’un contraste spécifique peut être réalisé tout aussi simplement, en indiquant le niveau du ou des facteurs qui nous intéressent :

Des simulations de Monte Carlo permettent évaluer la puissance statistique d’un plan factoriel, en spécifiant les valeurs attendues pour les moyennes et variances dans chacune des conditions expérimentales. Il ne s’agit donc pas d’un calcul de puissance a posteriori.

Références

[1] Nick Cox. Speaking stata: Graphing model diagnostics. The Stata Journal, 4(4):449--475, 2004.
[2] W. D. Dupont. Statistical Modeling for Biomedical Researchers. Cambridge University Press, 2nd edition, 2009.
[3] Douglas C. Montgomery. Design and Analysis of Experiments. Wiley, 5 edition, 2001.