À propos des facteurs sous R

Les fonction utiles pour créer et gérer des variables de type facteur (variables de classification, ordonnées ou non) sont essentiellement les fonctions factor et gl. À la différence de factor qui peut être utilisé pour spécifier un vecteur de conditions, gl ne s'applique qu'aux plans équilibrés (même nombre de répétition des niveaux du facteur). Il existe également l'opérateur C que l'on peut appliquer directement dans la formule définissant un modèle de test (voir section suivante), mais on l'utilise plutôt pour manipuler des contrastes. Par la suite, nous utiliserons indifféremment “groupe” ou “condition” pour désigner les niveaux d'un facteur.

Par défaut, des contrastes de type “traitement” sont utilisés. Dans cette approche proposée par Dunnett pour les comparaisons de différents groupes à un groupe de référence (baseline condition), chaque moyenne de groupe est comparée à la moyenne du groupe de référence. Par défaut également, R considère comme groupe de référence le premier niveau, défini par ordre lexicographique, du facteur considéré. Attention : Certains logiciels, comme SAS, ont une approche différente. Bien entendu, on peut changer le niveau de référence avec la commande relevel. Par exemple, pour un plan factoriel 23, si on définit un facteur grp comme

grp <- gl(3,2,labels=letters[1:3])

alors le niveau de référence est "a" (cf. l'affichage des niveaux en tapant simplement grp en ligne de commande). C'est normal car nous l'avons défini explicitement avec l'option labels=.

Maintenant, si l'on considère des niveaux de type "a","A","b","B" (en supposant que cela aît un sens), on constate que l'ordre des niveaux considéré par R respecte l'ordre lexicographique et ce sont les caractères en capitales qui apparaissent parmi les premiers niveaux. Dans l'exemple suivant, nous aurons donc "A" comme niveau de référence, quelque soit l'ordre de la séquence d'apparition des niveaux du facteur.

grp <- factor(sample(rep(c("a","A","b","B"),2),replace=F))

On peut redéfinir le niveau de référence d'un facteur comme suit :

grp <- relevel(grp,"a")

Le niveau de référence sera désormais "a" et non plus "A".

Gestion des contrastes dans les modèles linéaires

On peut définir les contrastes à utiliser dans les modèles linéaires (aov(), lm(), lme(), etc.) en les précisant dans la fonction options(), e.g.

op <- options(contrasts=c("contr.helmert", "contr.poly"))
summary(aov1 <- aov(y~x))
coef(aov1)
options(op)

Si l'on applique un modèle de type aov(y~x), les coefficients de régression sont indiqués grâce à la commande coef précédente ou simplement avec la fonction summary.lm. Celle-ci peut s'utiliser lorsque l'on utilise la fonction aov plutôt que lm. On vérifiera que ces coefficients correspondent bien aux différences de moyennes estimées entre chaque condition et la condition de référence. Il y a donc toujours k-1 coefficients, avec k modalités ou niveaux pour le facteur.

Pour retrouver l'approche ANOVA dans un modèle de régression, c'est-à-dire l'écart des moyennes de groupes à la moyenne générale, et non la comparaison par rapport à un groupe de référence, on peut spécifier directement le contraste lors de l'appel à la fonction lm(), e.g.

summary(aov1b <- lm(y ~ x,data=xy,contrasts=list(x="contr.sum")))

Autre solution, beaucoup plus rapide : il suffit de supprimer le terme d'ordonnée à l'origine (intercept), et les coefficients de régression traduisent alors les écarts à la moyenne générale, et non à la moyenne du groupe de référence.

summary(aov1c <- lm(y ~ -1 + x, data=xy))

Nous voyons donc l'importance de comprendre le codage et l'utilisation des contrastes sous R, sous peine de fournir une interprétation erronée des résultats. Simulons quelques données artificielles pour illustrer ces différences.

x <- gl(4,10,labels=paste("a",1:4,sep=""))
y <- as.numeric(x)*rnorm(40,10,2/(as.numeric(x)/4)) + rnorm(40)
dat1 <- data.frame(y=y,x=x)
boxplot(y~x, dat1)
1.png

Maintenant, appliquons le modèle d'ANOVA vu dans les paragraphes précédents. En premier lieu, summary(lm(y~x,dat1)) donne

Call:
lm(formula = y ~ x, data = dat1)

Residuals:
     Min       1Q   Median       3Q      Max
-13.1688  -4.0018   0.2298   3.9516  13.5267

Coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept)    9.143      2.073   4.410 8.97e-05 ***
xa2           12.007      2.932   4.095 0.000228 ***
xa3           26.285      2.932   8.965 1.06e-10 ***
xa4           34.594      2.932  11.798 6.28e-14 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 6.556 on 36 degrees of freedom
Multiple R-squared: 0.8197,     Adjusted R-squared: 0.8047
F-statistic: 54.57 on 3 and 36 DF,  p-value: 1.798e-13

La valeur 12.007 correspond bien à la différence entre la moyenne du 2ème groupe et celle du 1er groupe, comme on peut le vérifier très facilement :

grp.means <- with(dat1, tapply(y,x,mean))
grp.means[2:4]-grp.means[1]

Lorsqu'on supprime l'intercept (summary(lm(y~x-1,dat1))), on obtient :

Call:
lm(formula = y ~ x - 1, data = dat1)

Residuals:
     Min       1Q   Median       3Q      Max
-13.1688  -4.0018   0.2298   3.9516  13.5267

Coefficients:
    Estimate Std. Error t value Pr(>|t|)
xa1    9.143      2.073    4.41 8.97e-05 ***
xa2   21.150      2.073   10.20 3.64e-12 ***
xa3   35.428      2.073   17.09  < 2e-16 ***
xa4   43.736      2.073   21.09  < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 6.556 on 36 degrees of freedom
Multiple R-squared: 0.9598,     Adjusted R-squared: 0.9554
F-statistic: 215.1 on 4 and 36 DF,  p-value: < 2.2e-16

Il existe une option contrasts= pour la fonction aov qui, à la différence de la modificatio globale par la commande options, permet de spécifier des contrastes pour certains termes du modèle uniquement. On en verra un exemple d'utilisation plus tard.

La fonction lm() repose sur la matrice de design (?model.matrix), que l'on peut construire soi-même en spécifiant l'argument contrasts (sous forme de list). On peut vérifier les contrastes utilisés dans un modèle GLM comme ceci :

contrasts(model.matrix(fit))

Par défaut, R utilise des contrastes de type Helmert, e.g.

Ch <- contr.helmert(5)  ## les vecteurs sont orthogonaux 2 à 2
crossprod(Ch)           ## les éléments hors-diagonale sont nuls

Les contrastes ainsi spécifiés figurent toujours en colonne, pas en ligne ! Si on utilise lm() (ou summary.lm() avec un modèle spécifié grâce à aov()), les coefficients (hors intercept) correspondent dans l'ordre aux colonnes de la matrice de contraste.

Les comparaisons pour k=5 sont ainsi :

On a donc 4 comparaisons indépendantes, qui permettent de faire des comparaisons de moyennes de manière séquentielle en sélectionnant des niveaux adjacents.

D'autre part, ce qui est appelé contraste type traitement sous R ne constitue pas à proprement parler un vrai contraste (au sens de la théorie des modèles linéaires) car la contrainte d'orthogonalité par rapport à l'intercept n'est pas respectée. Toutefois, les vecteurs lignes sont bien orthogonaux 2 à 2.

      Helmert        Traitement
1   -1 -1 -1 -1      0  0  0  0
2    1 -1 -1 -1      1  0  0  0
3    0  2 -1 -1      0  1  0  0
4    0  0  3 -1      0  0  1  0
5    0  0  0  4      0  0  0  1

Pour avoir des résultats comparables avec ceux produits par SAS (PROC GLM), il faut utiliser contr.SAS à la place de contr.treatment pour spécifier que le niveau de base à considérer est le dernier et pas le premier. Autre solution : redéfinir les niveaux du facteur, e.g.

# facteur à 4 niveaux randomisé sur 20 sujets
x1 <- factor(sample(rep(letters[1:4],5)))
contr.treatment(levels(x1))  ## a est le niveau de base
x2 <- relevel(x1,"d")
contr.treatment(levels(x2))  ## d est le niveau de base

ou définir une matrice de contrastes appropriée, e.g.

Cm <- matrix(rbind(diag(3),rep(0,3)),ncol=3)

Dernière solution (sans doute la plus simple), spécifier l'argument base lors de l'appel à contr.treatment()

contr.treatment(levels(x1),base=nlevels(x1))
Exemple 1 : Analyse d'un plan équilibré
ni <- c(8,10,9,10)
x <- factor(rep(letters[1:4],ni))
y <- rnorm(sum(ni)) + as.numeric(x)  ## relation linéaire de pente ~ 1
summary(ols1 <- lm(y~x))             ## ou summary.lm(aov1 <- aov(y~x))

Ici l'intercept considéré est la moyenne dans le groupe "a" et les coefficients de régression représentent les différences de moyenne avec "a" :

mean(y[x=="a"])
(tapply(y,x,mean)-mean(y[x=="a"]))[-1]

summary(ols2 <- lm(y~x,contrasts=list(x="contr.sum")))
model.matrix(ols2)

On notera que la distribution des résidus reste la même.

# les données proviennent du livre de Howell
# D.C. Howell, Méthodes Statistiques en Sciences Humaines, De Boeck, 1998
# [1re ed., 6e tirage, 2004]
# http://howell.psp.ucl.ac.be/
a <- read.table("tab13-2.dat")
a$V1 <- factor(a$V1, labels = c("ages", "jeunes"))
a$V2 <- factor(a$V2, labels = c("addition", "rimes", "adjectifs", "images", "controle"))
summary(a)
# tableau d'ANOVA classique avec aov()
a.aov <- aov(V3 ~ V1 * V2, data = a)
summary(a.aov)
summary.lm(a.aov)
# Note: les tests t associés aux coefficients estimés des facteurs
# croisés n'ont aucune pertinence...
# avec le modèle linéaire et la fonction lm()
a.lm <- lm(V3~V1:V2-1, a, contrasts=list(V1=contr.sum,V2=contr.sum))
# l'option 'contrasts' est nécessaire pour obtenir une analyse de type III
# correcte
summary(a.lm)
anova(a.lm)
library(car)
Anova(a.lm, type="II")
Anova(a.lm, type="III")
library(Rcmdr)
plotMeans(V3,V1,V2)
levene.test(V3,V1:V2)

Utilisation du package Design

Références

  1. Chambers, J. M. and Hastie, T. J. (1992). Statistical models in S. Wiley.

  2. Dunnett, C. W. (1955). A multiple comparisons procedure for comparing several treatments with a control. Journal of the American Statistical Association, 50, 1096-1121.

  3. Harrell, F. E. (2001). Regression Modeling Strategies. With Applications to Linear Models, Logisitic Regression, and Survival Analysis. Springer.

  4. Alzola, C. and Harrell, F. E. (2006). An Introduction to S and the Hmisc and Design Libraries.

  5. Howell, D. C. (2004). Méthodes Statistiques en Sciences Humaines. De Boeck.