Contrastes sous R ================= 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) ------------------------------------------------------------------------------ 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 pas 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"))) ------------------------------------------------------------------------------ // Note: // Revoir comparaison par rapport à grande moyenne... 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) Références ^^^^^^^^^^ + Chambers, J. M. and Hastie, T. J. (1992). _Statistical models in S_. Wiley.