Jak wykonać ANOVA SS typu III w R z kodami kontrastu?

26

Proszę podać kod R, który pozwala przeprowadzić ANOVA między podmiotami z kontrastami -3, -1, 1, 3. Rozumiem, że jest debata na temat odpowiedniego rodzaju sumy kwadratów (SS) do takiej analizy. Jednak jako domyślny typ SS używany w SAS i SPSS (typ III) jest uważany za standard w mojej okolicy. Dlatego chciałbym, aby wyniki tej analizy idealnie pasowały do ​​tego, co generują te programy statystyczne. Aby zostać zaakceptowanym, odpowiedź musi bezpośrednio wywołać aov (), ale inne odpowiedzi mogą zostać poddane pod głosowanie (szczególnie jeśli są łatwe do zrozumienia / użycia).

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

Edycja: Proszę zauważyć, że żądany kontrast nie jest prostym kontrastem liniowym lub wielomianowym, ale kontrastem wyprowadzonym z teoretycznej prognozy, tj. Rodzaju kontrastów omówionych przez Rosenthala i Rosnowa.

russellpierce
źródło
5
Rozumiem, że potrzebujesz sumy typu III, ale ten artykuł ( stats.ox.ac.uk/pub/MASS3/Exegeses.pdf ) jest dobrym materiałem do przeczytania. Ilustruje kilka interesujących punktów.
suncoolsu,
Jeśli chodzi o twoje pytanie, możesz zainteresować się następującą dyskusją: stats.stackexchange.com/questions/60362/… Wybór między ANOVA typu I, II i III nie jest tak łatwy, jak się wydaje.
phx
Poparłem twoje pytanie jako przydatne, ponieważ wywołało kilka wyuczonych odpowiedzi, ale zauważam również, że zgodziłeś się z respondentem, który w zasadzie stwierdził, że przesłanka pytania była nieprawidłowa. Mam nadzieję, że podsumowuję stanowisko StaGuy'a, mówiąc, że zdefiniowane kontrasty były z definicji „typem I”, a dyskusja na temat innych typów stała się istotna tylko podczas oceny statystyk częściowej regresji, prawdopodobnie najważniejszej przy pozwoleniu „maszynie na prowadzenie” za pomocą metod automatycznych.
DWin
@DWin: Nie jestem pewien, czy całkowicie cię śledzę. Można legalnie korzystać z innych typów SS, nie pozwalając „maszynie prowadzić” (przynajmniej tak rozumiem to zdanie). Mogę być trochę zardzewiały, ale jeśli pamięć służy, inne typy mogą być istotne, gdy nie stosuję częściowej regresji. Na przykład SS typu III nie eliminuje głównych efektów interakcji. Różnica między typami ma tam znaczenie właśnie dlatego, że Typ III nie jest częściowy, podczas gdy Typ I ma. Problem, jak stwierdzono, obejmował tylko jeden kontrast, a zatem rozróżnienie między rodzajami SS było / jest dyskusyjne.
russellpierce
Zrozumiałem, że uzasadnieniem podanym przez SAS przy wyborze SSS typu III (i wydaje się, że to dlatego ludzie uważają, że preferowany jest typ III) jest to, że lepiej wspiera proces selekcji do tyłu i do przodu.
DWin

Odpowiedzi:

22

Suma kwadratów typu III dla ANOVA jest łatwo dostępna poprzez Anova()funkcję z pakietu samochodowego .

Kontrast kodowanie może być wykonywane na wiele sposobów, przy użyciu C(), na contr.*rodziny (jak wskazano przez @nico) lub bezpośrednio contrasts()funkcji / argumentu. Jest to szczegółowo opisane w § 6.2 (s. 144–151) Modern Applied Statistics with S (Springer, 2002, wydanie 4). Zauważ, że aov()jest to tylko funkcja otoki dla tej lm()funkcji. Jest to interesujące, gdy chce się kontrolować termin błędu modelu (jak w projekcie wewnątrz przedmiotu), ale w przeciwnym razie oba dają te same wyniki (i niezależnie od tego, w jaki sposób dopasujesz swój model, nadal możesz wygenerować ANOVA lub LM- jak streszczenia z summary.aovlub summary.lm).

Nie mam SPSS do porównania dwóch wyjść, ale coś w tym rodzaju

> 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 

warto spróbować w pierwszej instancji.

O kodowaniu czynnikowym w R vs. SAS: R uważa poziom bazowy lub referencyjny za pierwszy poziom w porządku leksykograficznym, podczas gdy SAS uważa ostatni. Tak więc, aby uzyskać porównywalne wyniki, albo masz do wykorzystania contr.SAS()lub do relevel()swojej czynnika R.

chl
źródło
1
Nie sądzę, aby w tej odpowiedzi użyto podanego przeze mnie kontrastu -3, -1,1,3, ani nie wydaje się, aby zapewniał test kontrastu 1 df.
russellpierce
@drknexus Tak, masz rację. Napisałem zbyt szybko. Coś takiego Anova(lm(DV ~ C(IV, c(-3,-1,1,3),1), data=sample.data), type="III")powinno być lepsze. Daj mi znać, jeśli ci to odpowiada.
chl
Dzięki! Wygląda na to, że poprawię to jutro w stosunku do SPSS i skontaktuję się z Tobą.
russellpierce
1
BTW, spójrz na pakiet ez ( cran.r-project.org/web/packages/ez/index.html ) do owijania kodu Anova ...
Tal Galili
2
@drknexus: Gdyby tylko była prośba o funkcję i strona przesyłania problemów dla ez ... github.com/mike-lawrence/ez/issues :)
Mike Lawrence
11

Może to wyglądać na autopromocję (i przypuszczam, że tak). Ale opracowałem pakiet lsmeans dla R (dostępny w CRAN), który został zaprojektowany do obsługi dokładnie tego rodzaju sytuacji. Oto jak to działa na twoim przykładzie:

> 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

Jeśli chcesz, możesz określić dodatkowe kontrasty na liście. W tym przykładzie uzyskasz takie same wyniki z wbudowanym liniowym kontrastem wielomianowym:

> 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

Aby to potwierdzić, należy zauważyć, że "poly"specyfikacja kieruje go do wywołania poly.lsmc, co daje następujące wyniki:

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

Jeśli chcesz wykonać wspólny test kilku kontrastów, użyj testfunkcji za pomocą joint = TRUE. Na przykład,

> test(con, joint = TRUE)

Spowoduje to wykonanie testu „typu III”. W przeciwieństwie do car::Anova()tego zrobi to poprawnie, niezależnie od kodowania kontrastu stosowanego na etapie dopasowania modelu. Wynika to z faktu, że testowane funkcje liniowe są określone bezpośrednio, a nie domyślnie, poprzez redukcję modelu. Dodatkową cechą jest wykrycie przypadku, w którym testowane kontrasty są liniowo zależne, oraz wygenerowanie prawidłowej statystyki testowej i stopni swobody.

rvl
źródło
7

Kiedy robisz kontrasty, robisz określoną, określoną liniową kombinację średnich komórek w kontekście odpowiedniego terminu błędu. W związku z tym koncepcja „typu SS” nie ma znaczenia w przypadku kontrastów. Każdy kontrast jest zasadniczo pierwszym efektem przy użyciu SS typu I. „Rodzaj SS” ma związek z tym, co jest podzielone lub uwzględnione w innych warunkach. W przypadku kontrastów nic nie jest podzielone ani rozliczane. Kontrast jest sam w sobie.

StatGuy
źródło
Masz absolutną rację.
russellpierce
3

To, że testy typu III są używane w miejscu pracy, jest najsłabszym powodem do ich dalszego używania. SAS wyrządził poważne szkody w statystykach w tym zakresie. Egzegeza Billa Venablesa, o której mowa powyżej, jest świetnym źródłem informacji na ten temat. Po prostu powiedz „nie” typowi III; opiera się na błędnym pojęciu równowagi i ma niższą moc z powodu głupiego ważenia komórek w niezrównoważonym przypadku.

Bardziej naturalny i mniej podatny na błędy sposób uzyskiwania ogólnych kontrastów i opisywania tego, co zrobiłeś, zapewnia funkcja rmspakietu R. contrast.rmsKontrasty mogą być bardzo złożone, ale dla użytkownika są bardzo proste, ponieważ są wyrażone w kategoriach różnic wartości predykcyjnych. Obsługiwane są testy i jednoczesne kontrasty. Obsługuje to efekty regresji nieliniowej, efekty interakcji nieliniowej, częściowe kontrasty, wszelkiego rodzaju rzeczy.

Frank Harrell
źródło
To wszystko w porządku i dobre dla ciebie jako osoby o ugruntowanej reputacji. Inni nie mają siły, by nie zgodzić się z recenzentami. Ponieważ interpretacje statystyk różnią się, będziesz prosił nowych ludzi, aby zasadniczo postępowali i ponosili nieuzasadnione koszty. Mówię to jako osoba, która umarła przez mój czas na tych (i podobnych) wzgórzach. Zmiany IMO na tym froncie są obowiązkiem strażników, tj. Redaktorów i recenzentów.
russellpierce
Ludzie, którzy naprawdę dobrze posługują się danymi, mają szeroki wybór miejsc pracy i mogą mieć możliwość pracy w obszarach, w których ich umiejętności i opinie są szanowane.
Frank Harrell
1
... i to właśnie robię teraz. Ale ludzie, którzy przychodzą na to pytanie, często nie należą do tej klasy. Tak jak nie byłem 7 lat temu. Opowiadam się tylko za odrobiną empatii wobec nowicjusza i jego okoliczności.
russellpierce
2

Wypróbuj polecenie Anova w bibliotece samochodów. Użyj argumentu type = "III", ponieważ domyślnie jest to typ II. Na przykład:

library(car)
mod <- lm(conformity ~ fcategory*partner.status, data=Moore, contrasts=list(fcategory=contr.sum, partner.status=contr.sum))
Anova(mod, type="III")
jebyrnes
źródło
3
Wiem, że Moore znajduje się w bibliotece samochodowej, ale po dostarczeniu przykładowych danych łatwiej jest pytającemu zrozumieć twoją odpowiedź, jeśli użyjesz przykładowych danych.
russellpierce
0

Również autopromocja, napisałem funkcję do tego dokładnie: https://github.com/samuelfranssens/type3anova

Zainstaluj w następujący sposób:

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))

Będziesz także musiał carzainstalować pakiet.

sam_f
źródło
Jak zastosowałbyś to do części pytania dotyczącej kontrastów?
russellpierce
1
Przepraszam, nie przeczytałem pytania poprawnie. Moja funkcja uprości tylko przeprowadzanie Anova typu III. Podobnie jak powyżej StatGuy, nie widzę, gdzie SS wchodzi w grę podczas testowania określonych kontrastów.
sam_f,