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

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.

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.

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.

Możesz rzucić okiem na ten post na blogu:

Uzyskanie tej samej ANOVA skutkuje R jak w SPSS - trudności z sumami kwadratów typu II i typu III

( Spoiler: dodaj options(contrasts=c("contr.sum", "contr.poly"))na początku skryptu)

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.

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.

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

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.