Comment fait-on une ANOVA SS de type III en R avec des codes de contraste?

26

Veuillez fournir le code R qui permet d'effectuer une ANOVA entre sujets avec -3, -1, 1, 3 contrastes. Je comprends qu'il y a un débat concernant le type de somme de carrés (SS) approprié pour une telle analyse. Cependant, comme le type par défaut de SS utilisé dans SAS et SPSS (Type III) est considéré comme la norme dans ma région. J'aimerais donc que les résultats de cette analyse correspondent parfaitement à ce qui est généré par ces programmes statistiques. Pour être acceptée, une réponse doit appeler directement aov (), mais d'autres réponses peuvent être votées (surtout si elles sont faciles à comprendre / à utiliser).

sample.data <- data.frame(IV=rep(1:4,each=20),DV=rep(c(-3,-3,1,3),each=20)+rnorm(80))

Edit: Veuillez noter que le contraste que je demande n'est pas un simple contraste linéaire ou polynomial mais est un contraste dérivé d'une prédiction théorique, c'est-à-dire le type de contrastes discuté par Rosenthal et Rosnow.

r anova contrasts sums-of-squares russellpierce
la source
5
Je comprends que vous avez besoin d'une somme de type III, mais cet article ( stats.ox.ac.uk/pub/MASS3/Exegeses.pdf ) est une bonne lecture. Il illustre quelques points intéressants.
suncoolsu
Concernant votre question, vous pourriez être intéressé par la discussion suivante: stats.stackexchange.com/questions/60362/… Le choix entre ANOVA type I, II et III n'est pas aussi simple qu'il y paraît.
phx
Votre question a été jugée utile en ce qu'elle a provoqué plusieurs réponses apprises, mais je note également que vous étiez d'accord avec le répondant qui a essentiellement dit que la prémisse de la question était incorrecte. J'espère que je résume la position de StaGuy en disant que les contrastes définis étaient par définition de "type I" et que la discussion sur d'autres types n'était pertinente que lors de l'évaluation des statistiques de régression partielle, probablement plus importante lorsque l'on laissait "la machine faire la conduite" à l'aide de méthodes automatisées.
DWin
@DWin: Je ne suis pas sûr de vous suivre entièrement. On peut légitimement utiliser d'autres types de SS sans laisser la «machine faire la conduite» (du moins si je comprends bien cette phrase). Je pourrais être un peu rouillé ici, mais si la mémoire est bonne, d'autres types peuvent être pertinents lorsqu'ils n'utilisent pas de régression partielle. Par exemple, Type III SS ne supprime pas les principaux effets de l'interaction. La distinction entre les types y importe précisément parce que le type III n'est pas partiel alors que le type I le fait. Le problème, comme indiqué, ne comprenait qu'un seul contraste et, par conséquent, la distinction entre les types de SS était / est sans objet.
russellpierce
Ma compréhension était que la justification donnée par SAS pour choisir le SSS de type III (et cela semble être la raison pour laquelle les gens pensent que le type III est préféré) est qu'il soutient mieux le processus de sélection en amont et en aval.
DWin

Réponses:

22

La somme des carrés de type III pour l'ANOVA est facilement disponible via la Anova()fonction du package de voiture .

Le codage de contraste peut être effectué de plusieurs manières, en utilisant C()la contr.*famille (comme indiqué par @nico) ou directement la contrasts()fonction / l'argument. Ceci est détaillé au §6.2 (pp. 144-151) de Modern Applied Statistics with S (Springer, 2002, 4e éd.). Notez que aov()c'est juste une fonction wrapper pour la lm()fonction. C'est intéressant quand on veut contrôler le terme d'erreur du modèle (comme dans une conception intra-sujet), mais sinon ils donnent tous les deux les mêmes résultats (et quelle que soit la façon dont vous ajustez votre modèle, vous pouvez toujours générer ANOVA ou LM- comme des résumés avec summary.aovou summary.lm).

Je n'ai pas SPSS pour comparer les deux sorties, mais quelque chose comme

> library(car)
> sample.data <- data.frame(IV=factor(rep(1:4,each=20)),
                            DV=rep(c(-3,-3,1,3),each=20)+rnorm(80))
> Anova(lm1 <- lm(DV ~ IV, data=sample.data, 
                  contrasts=list(IV=contr.poly)), type="III")
Anova Table (Type III tests)

Response: DV
            Sum Sq Df F value    Pr(>F)    
(Intercept)  18.08  1  21.815  1.27e-05 ***
IV          567.05  3 228.046 < 2.2e-16 ***
Residuals    62.99 76                      
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

vaut la peine d'essayer en première instance.

À propos du codage factoriel dans R par rapport à SAS: R considère le niveau de référence ou de référence comme le premier niveau dans l'ordre lexicographique, tandis que SAS considère le dernier. Donc, pour obtenir des résultats comparables, vous devez utiliser contr.SAS()ou utiliser relevel()votre facteur R.

chl
la source
1
Je ne pense pas que cette réponse utilise le contraste -3, -1,1,3 que j'ai spécifié ni ne semble fournir un test de 1 df du contraste.
russellpierce
@drknexus Oui, vous avez raison. A écrit trop vite. Quelque chose comme ça Anova(lm(DV ~ C(IV, c(-3,-1,1,3),1), data=sample.data), type="III")devrait être mieux. Veuillez me faire savoir si cela vous convient.
chl
Merci! Cela semble correct, je le validerai contre SPSS demain et je vous répondrai.
russellpierce
1
BTW, jetez un œil au paquet ez ( cran.r-project.org/web/packages/ez/index.html ) pour envelopper le code Anova ...
Tal Galili
2
@drknexus: Si seulement il y avait une page de demande de fonctionnalité et de soumission de problèmes pour ez ... github.com/mike-lawrence/ez/issues :)
Mike Lawrence
11

Cela peut ressembler à un peu d'auto-promotion (et je suppose que oui). Mais j'ai développé un package lsmeans pour R (disponible sur CRAN) qui est conçu pour gérer exactement ce genre de situation. Voici comment cela fonctionne pour votre exemple:

> sample.data <- data.frame(IV=rep(1:4,each=20),DV=rep(c(-3,-3,1,3),each=20)+rnorm(80))
> sample.aov <- aov(DV ~ factor(IV), data = sample.data)

> library("lsmeans")
> (sample.lsm <- lsmeans(sample.aov, "IV"))
 IV    lsmean        SE df   lower.CL  upper.CL
  1 -3.009669 0.2237448 76 -3.4552957 -2.564043
  2 -3.046072 0.2237448 76 -3.4916980 -2.600445
  3  1.147080 0.2237448 76  0.7014539  1.592707
  4  3.049153 0.2237448 76  2.6035264  3.494779

> contrast(sample.lsm, list(mycon = c(-3,-1,1,3)))
 contrast estimate       SE df t.ratio p.value
 mycon    22.36962 1.000617 76  22.356  <.0001

Vous pouvez spécifier des contrastes supplémentaires dans la liste si vous le souhaitez. Pour cet exemple, vous obtiendrez les mêmes résultats avec le contraste polynomial linéaire intégré:

> con <- contrast(sample.lsm, "poly")
> con
 contrast   estimate        SE df t.ratio p.value
 linear    22.369618 1.0006172 76  22.356  <.0001
 quadratic  1.938475 0.4474896 76   4.332  <.0001
 cubic     -6.520633 1.0006172 76  -6.517  <.0001

Pour confirmer cela, notez que la "poly"spécification lui ordonne d'appeler poly.lsmc, ce qui produit ces résultats:

> poly.lsmc(1:4)
  linear quadratic cubic
1     -3         1    -1
2     -1        -1     3
3      1        -1    -3
4      3         1     1

Si vous souhaitez effectuer un test conjoint de plusieurs contrastes, utilisez la testfonction avec joint = TRUE. Par exemple,

> test(con, joint = TRUE)

Cela produira un essai de "type III". Contrairement à car::Anova()cela, il le fera correctement quel que soit le codage de contraste utilisé dans l'étape d'ajustement du modèle. En effet, les fonctions linéaires testées sont spécifiées directement plutôt qu'implicitement via la réduction de modèle. Une caractéristique supplémentaire est qu'un cas où les contrastes testés sont linéairement dépendants est détecté et que les statistiques de test et les degrés de liberté corrects sont produits.

rvl
la source
7

Lorsque vous effectuez des contrastes, vous effectuez une combinaison linéaire spécifique et indiquée de moyennes de cellules dans le contexte du terme d'erreur approprié. En tant que tel, le concept de "Type de SS" n'a pas de sens avec les contrastes. Chaque contraste est essentiellement le premier effet utilisant un SS de type I. Le «type de SS» a à voir avec ce qui est séparé ou pris en compte par les autres termes. Pour les contrastes, rien n'est parti ni pris en compte. Le contraste tient tout seul.

StatGuy
la source
Tu as tout à fait raison.
russellpierce
3

Le fait que les tests de type III soient utilisés sur votre lieu de travail est la raison la plus faible pour continuer à les utiliser. SAS a gravement endommagé les statistiques à cet égard. L'exégèse de Bill Venables, mentionnée ci-dessus, est une excellente ressource à ce sujet. Dites simplement non au type III; il est basé sur une notion erronée d'équilibre et a une puissance inférieure en raison de la pondération stupide des cellules dans le cas déséquilibré.

Un moyen plus naturel et moins sujet aux erreurs d'obtenir des contrastes généraux et de pouvoir décrire ce que vous avez fait est fourni par la fonction de rmspackage R. contrast.rmsLes contrastes peuvent être très complexes, mais pour l'utilisateur, ils sont très simples car ils sont exprimés en termes de différences de valeurs prédictives. Les tests et les contrastes simultanés sont pris en charge. Cela gère les effets de régression non linéaire, les effets d'interaction non linéaire, les contrastes partiels, toutes sortes de choses.

Frank Harrell
la source
C'est très bien et bon pour vous en tant que personne de renom de dire. D'autres n'ont pas le poids d'être en désaccord avec les examinateurs. Étant donné que les interprétations des statistiques diffèrent, vous demanderiez à de nouvelles personnes de respecter le principe et d'engager un coût déraisonnable. Je dis cela comme quelqu'un qui est mort ma part de temps au sommet de ces collines (et similaires). Le changement de l'OMI sur ce front est la responsabilité des gardiens, c'est-à-dire des éditeurs et des examinateurs.
russellpierce
Les gens qui sont vraiment bons avec les données ont un large choix d'emplois et peuvent avoir la possibilité de travailler dans des domaines où leurs compétences et leurs opinions sont respectées.
Frank Harrell
1
... et c'est ce que je fais maintenant. Mais les gens qui arrivent à cette question ne seront pas souvent de cette classe. Tout comme je ne l'étais pas il y a 7 ans. Je plaide seulement pour un peu d'empathie pour le novice et sa situation.
russellpierce
2

Essayez la commande Anova dans la bibliothèque de voitures. Utilisez l'argument type = "III", car il est par défaut de type II. Par exemple:

library(car)
mod <- lm(conformity ~ fcategory*partner.status, data=Moore, contrasts=list(fcategory=contr.sum, partner.status=contr.sum))
Anova(mod, type="III")
jebyrnes
la source
3
Je sais que Moore se trouve dans la bibliothèque de la voiture, mais lorsque des exemples de données sont fournis, il est plus facile pour le demandeur de questions de comprendre votre réponse si vous utilisez les exemples de données.
russellpierce
0

Aussi auto-promotionnel, j'ai écrit une fonction pour exactement cela: https://github.com/samuelfranssens/type3anova

Installez comme suit:

library(devtools)
install_github(samuelfranssens/type3anova)
library(type3anova)

sample.data <- data.frame(IV=rep(1:4,each=20),DV=rep(c(-3,-3,1,3),each=20)+rnorm(80))

type3anova(lm(DV ~ IV, data = sample.data))

Vous devrez également carinstaller le package.

sam_f
la source
Comment appliqueriez-vous cela à la partie contrastée de la question?
russellpierce
1
Toutes mes excuses, je n'ai pas lu la question correctement. Ma fonction ne fera que simplifier la réalisation d'Anova de type III. Comme StatGuy ci-dessus, je ne vois pas où SS entre en jeu lors du test de contrastes spécifiques.
sam_f