// pour compiler ce document : // $ asciidoc -a encoding=ISO-8859-1 contrasts.txt Contrastes sous R ================= Christophe Lalanne Mai 2008 ****************************************************************************** Cette note a été révisée en septembre 2008. Les principales mises à jour sont des compléments sur la formulation des constrastes à l'aide des packages +Design+ et +multcomp+ dans le cadre des comparaisons multiples, ainsi que les exemples en fin de document. ****************************************************************************** À 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 2^3^, 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) ------------------------------------------------------------------------------ image::1.png[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 ------------------------------------------------------------------------------ // http://pidgin.ucsd.edu/pipermail/r-lang/2007-May/000025.html 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 : - {2} vs. {1} - {3} vs. {1,2} - {4} vs. {1,2,3} - {5} vs. {1,2,3,4} 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)) ------------------------------------------------------------------------------ [caption="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] ------------------------------------------------------------------------------ Si on veut comparer les traitements deux à deux ({1} vs. {4}, {2} vs. {4} et {3} vs. {4}), on spécifie l'option lors de l'appel à +lm()+ : ------------------------------------------------------------------------------ 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 ---------- + Chambers, J. M. and Hastie, T. J. (1992). _Statistical models in S_. Wiley. + Dunnett, C. W. (1955). A multiple comparisons procedure for comparing several treatments with a control. _Journal of the American Statistical Association_, _50_, 1096-1121. + Harrell, F. E. (2001). _Regression Modeling Strategies. With Applications to Linear Models, Logisitic Regression, and Survival Analysis_. Springer. + Alzola, C. and Harrell, F. E. (2006). http://biostat.mc.vanderbilt.edu/twiki/pub/Main/RS/sintro.pdf[_An Introduction to S and the Hmisc and Design Libraries_]. [pdf] + Howell, D. C. (2004). _Méthodes Statistiques en Sciences Humaines_. De Boeck.