¿Cómo configurar e interpretar los contrastes ANOVA con el paquete del automóvil en R?

15

Digamos que tengo un simple experimento factorial 2x2 en el que quiero hacer ANOVA. Así, por ejemplo:

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. En ausencia de una interacción significativa, por defecto (es decir, contr.treatmentel resultado de Anova()es la importancia general de asobre todos los niveles de by bsobre todos los niveles de a, ¿es correcto?

  2. ¿Cómo debería especificar un contraste que me permitiera poner a prueba la importancia del efecto acon el bque se mantiene constante en el nivel B1, de efecto acon el bque se mantiene constante en el nivel B2, y de la interacción a:b?

f1r3br4nd
fuente

Respuestas:

18

Su ejemplo conduce a tamaños de celdas desiguales, lo que significa que los diferentes "tipos de suma de cuadrados" importan, y la prueba de los efectos principales no es tan simple como lo dice. Anova()usa la suma de cuadrados tipo II. Ver esta pregunta para empezar.

Hay diferentes formas de probar los contrastes. Tenga en cuenta que los tipos de SS no importan ya que finalmente estamos probando en el diseño factorial asociado. Sugiero usar los siguientes pasos:

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

Ahora compruebe manualmente el resultado para el primer contraste.

> 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
lince
fuente
33
¡¡¡GRACIAS!!! Acaba de responder una pregunta que dos semestres de estadísticas de nivel de posgrado no tienen. Incluso consideré usar anova unidireccional antes, pero no pude encontrar ninguna confirmación de que este fuera un enfoque legítimo.
f1r3br4nd
@ f1r3br4nd Es legítimo ya que el error MS es igual en el diseño unidireccional asociado y el diseño bidireccional original.
caracal
Una última pregunta de seguimiento, si puedo: ¿cómo se generaliza la interacción bidireccional a las interacciones de un mayor número de variables? Si tuviera un término A B C, ¿lo construiría a partir de A: B = (A | B = 1 - A | B = 2), C: B = (C | B = 1 - C | B = 2 ), A: B: C = A: B - C: B, y así sucesivamente.
f1r3br4nd
2
@ f1r3br4nd En un diseño de 2x2x2, solo hay un contraste de interacción A B C único (como si hubiera solo uno en un caso de 2x2). Los coeficientes en un contraste de interacción A B C tienen que sumar a cero sobre filas (A), columnas (B) y planos (C) en el "cubo de diseño". Si el orden de las celdas en el diseño unidireccional asociado es a1.b1.c1, a2.b1.c1, a1.b2.c1, a2.b2.c1, a1.b1.c2, a2.b1.c2, a1.b2.c2, a2.b2.c2, entonces los coeficientes son c(1, -1, -1, 1, -1, 1, 1, -1). Si tiene más de dos grupos en sus factores, todos los contrastes que siguen la regla de suma a cero son contrastes de interacción de 3 vías.
caracal