Jak skonfigurować i zinterpretować kontrasty ANOVA z pakietem samochodowym w R?

15

Powiedzmy, że mam prosty eksperyment czynnikowy 2x2, na którym chcę wykonać ANOVA. W ten sposób na przykład:

d   <- data.frame(a=factor(sample(c('a1','a2'), 100, rep=T)),
                  b=factor(sample(c('b1','b2'), 100, rep=T)));
d$y <- as.numeric(d$a)*rnorm(100, mean=.75, sd=1) +
       as.numeric(d$b)*rnorm(100, mean=1.2, sd=1) +
       as.numeric(d$a)*as.numeric(d$b)*rnorm(100, mean=.5, sd=1) +
       rnorm(100);
  1. W przypadku braku znaczącej interakcji, domyślnie (tj. contr.treatment) Wynikiem Anova()jest ogólne znaczenie dla awszystkich poziomów bi dla bwszystkich poziomów a, czy to prawda?

  2. Jak powinienem określić kontrast, który pozwoliłby mi przetestować istotność efektu aprzy butrzymywaniu stałej na poziomie b1, efektu aprzy butrzymywaniu stałej na poziomie b2 i interakcji a:b?

f1r3br4nd
źródło

Odpowiedzi:

18

Twój przykład prowadzi do nierównych rozmiarów komórek, co oznacza, że ​​różne „rodzaje sum kwadratów” mają znaczenie, a test na główne efekty nie jest tak prosty, jak to określiłeś. Anova()wykorzystuje sumę kwadratów typu II. Zobacz to pytanie na początek.

Istnieją różne sposoby testowania kontrastów. Zauważ, że typy SS nie mają znaczenia, ponieważ ostatecznie testujemy w powiązanym projekcie jednoczynnikowym. Sugeruję wykonanie następujących kroków:

# turn your 2x2 design into the corresponding 4x1 design using interaction()
> d$ab <- interaction(d$a, d$b)       # creates new factor coding the 2*2 conditions
> levels(d$ab)                        # this is the order of the 4 conditions
[1] "a1.b1" "a2.b1" "a1.b2" "a2.b2"

> aovRes <- aov(y ~ ab, data=d)       # oneway ANOVA using aov() with new factor

# specify the contrasts you want to test as a matrix (see above for order of cells)
> cntrMat <- rbind("contr 01"=c(1, -1,  0,  0),  # coefficients for testing a within b1
+                  "contr 02"=c(0,  0,  1, -1),  # coefficients for testing a within b2
+                  "contr 03"=c(1, -1, -1,  1))  # coefficients for interaction

# test contrasts without adjusting alpha, two-sided hypotheses
> library(multcomp)                   # for glht()
> summary(glht(aovRes, linfct=mcp(ab=cntrMat), alternative="two.sided"),
+         test=adjusted("none"))
Simultaneous Tests for General Linear Hypotheses
Multiple Comparisons of Means: User-defined Contrasts
Fit: aov(formula = y ~ ab, data = d)

Linear Hypotheses:
              Estimate Std. Error t value Pr(>|t|)
contr 01 == 0  -0.7704     0.7875  -0.978    0.330
contr 02 == 0  -1.0463     0.9067  -1.154    0.251
contr 03 == 0   0.2759     1.2009   0.230    0.819
(Adjusted p values reported -- none method)    

Teraz ręcznie sprawdź wynik pierwszego kontrastu.

> P       <- 2                             # number of levels factor a
> Q       <- 2                             # number of levels factor b
> Njk     <- table(d$ab)                   # cell sizes
> Mjk     <- tapply(d$y, d$ab, mean)       # cell means
> dfSSE   <- sum(Njk) - P*Q                # degrees of freedom error SS
> SSE     <- sum((d$y - ave(d$y, d$ab, FUN=mean))^2)    # error SS
> MSE     <- SSE / dfSSE                   # mean error SS
> (psiHat <- sum(cntrMat[1, ] * Mjk))      # contrast estimate
[1] -0.7703638

> lenSq <- sum(cntrMat[1, ]^2 / Njk)       # squared length of contrast
> (SE   <- sqrt(lenSq*MSE))                # standard error
[1] 0.7874602

> (tStat <- psiHat / SE)                   # t-statistic
[1] -0.9782893

> (pVal <- 2 * (1-pt(abs(tStat), dfSSE)))  # p-value
[1] 0.3303902
karakal
źródło
3
DZIĘKUJĘ CI!!! Właśnie odpowiedziałeś na pytanie, którego nie uzyskały dwa semestry statystyk na poziomie magisterskim. Wcześniej rozważałem użycie jednokierunkowej anova, ale nie znalazłem potwierdzenia, że ​​było to uzasadnione podejście.
f1r3br4nd
@ f1r3br4nd Jest to uzasadnione, ponieważ błąd MS jest równy w powiązanym jednokierunkowym i oryginalnym dwukierunkowym projekcie.
caracal
Ostatnie pytanie, jeśli mogę: w jaki sposób interakcja dwukierunkowa uogólnia się na interakcje większej liczby zmiennych? Gdybym miał termin A B C, czy zbudowałbym go z A: B = (A | B = 1 - A | B = 2), C: B = (C | B = 1 - C | B = 2 ), A: B: C = A: B - C: B i tak dalej?
f1r3br4nd
2
@ f1r3br4nd W projekcie 2x2x2 istnieje tylko jeden unikalny kontrast interakcji A B C (tak jak jest tylko jeden w przypadku 2x2). Współczynniki w kontraście interakcji A B C muszą sumować się do zera ponad wierszami (A), kolumnami (B) i płaszczyznami (C) w „sześcianie projektowym”. Jeśli kolejność komórek w powiązanym projekcie jednostronnym jest a1.b1.c1, a2.b1.c1, a1.b2.c1, a2.b2.c1, a1.b1.c2, a2.b1.c2, a1.b2.c2, a2.b2.c2, wówczas współczynniki są c(1, -1, -1, 1, -1, 1, 1, -1). Jeśli masz więcej niż dwie grupy w swoich czynnikach, wszystkie kontrasty według zasady sumowania do zera są kontrastami interakcji 3-kierunkowej.
caracal