¿Cómo se hace un ANOVA SS tipo III en R con códigos de contraste?

Proporcione un código R que le permita a uno realizar un ANOVA entre sujetos con -3, -1, 1, 3 contrastes. Entiendo que hay un debate con respecto al tipo apropiado de Suma de cuadrados (SS) para dicho análisis. Sin embargo, como el tipo predeterminado de SS utilizado en SAS y SPSS (Tipo III) se considera el estándar en mi área. Por lo tanto, me gustaría que los resultados de este análisis coincidan perfectamente con lo que generan esos programas de estadísticas. Para ser aceptada, una respuesta debe llamar directamente a aov (), pero se pueden votar otras respuestas (especialmente si son fáciles de entender / usar).

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

Editar: Tenga en cuenta que el contraste que solicito no es un simple contraste lineal o polinomial, sino que es un contraste derivado de una predicción teórica, es decir, el tipo de contraste discutido por Rosenthal y Rosnow.

La suma de cuadrados de tipo III para ANOVA está fácilmente disponible a través de la Anova()función del paquete del automóvil .

La codificación de contraste se puede hacer de varias maneras, utilizando C()la contr.*familia (como lo indica @nico) o directamente la contrasts()función / argumento. Esto se detalla en §6.2 (págs. 144-151) de Estadísticas modernas aplicadas con S (Springer, 2002, 4a ed.). Tenga en cuenta que aov()es solo una función de contenedor para la lm()función. Es interesante cuando uno quiere controlar el término de error del modelo (como en un diseño dentro del tema), pero de lo contrario ambos producen los mismos resultados (y cualquiera que sea la forma en que se ajusta a su modelo, aún puede generar ANOVA o LM- como resúmenes con summary.aovo summary.lm).

No tengo SPSS para comparar las dos salidas, pero algo así como

> 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 

Vale la pena probarlo en primera instancia.

Acerca de la codificación de factores en R vs. SAS: R considera la línea base o nivel de referencia como el primer nivel en orden lexicográfico, mientras que SAS considera el último. Entonces, para obtener resultados comparables, ya sea que tenga que usar contr.SAS()o relevel()su factor R.

Esto puede parecer un poco de autopromoción (y supongo que lo es). Pero desarrollé un paquete lsmeans para R (disponible en CRAN) que está diseñado para manejar exactamente este tipo de situación. Así es como funciona para su ejemplo:

> 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

Si lo desea, puede especificar contrastes adicionales en la lista. Para este ejemplo, obtendrá los mismos resultados con el contraste polinómico lineal incorporado:

> 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

Para confirmar esto, tenga en cuenta que la "poly"especificación lo dirige a llamar poly.lsmc, lo que produce estos resultados:

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

Si desea hacer una prueba conjunta de varios contrastes, use la testfunción con joint = TRUE. Por ejemplo,

> test(con, joint = TRUE)

Esto producirá una prueba de "tipo III". A diferencia car::Anova(), lo hará correctamente independientemente de la codificación de contraste utilizada en la etapa de ajuste del modelo. Esto se debe a que las funciones lineales que se prueban se especifican directamente en lugar de implícitamente a través de la reducción del modelo. Una característica adicional es que se detecta un caso en el que los contrastes que se prueban dependen linealmente y se produce la estadística de prueba correcta y los grados de libertad.

Es posible que desee echar un vistazo a esta publicación de blog:

La obtención de los mismos resultados ANOVA en R que en SPSS - las dificultades con las sumas de cuadrados de Tipo II y Tipo III

( Spoiler: agregue options(contrasts=c("contr.sum", "contr.poly"))al comienzo de su secuencia de comandos)

Cuando está haciendo contrastes, está haciendo una combinación lineal específica y establecida de medias de celda dentro del contexto del término de error apropiado. Como tal, el concepto de "Tipo de SS" no es significativo con los contrastes. Cada contraste es esencialmente el primer efecto que usa un SS Tipo I. El "Tipo de SS" tiene que ver con lo que está parcializado o explicado por los otros términos. Para los contrastes, nada se divide ni se explica. El contraste se destaca por sí mismo.

El hecho de que las pruebas de tipo III se usen en su lugar de trabajo es la razón más débil para seguir usándolas. SAS ha hecho un daño importante a las estadísticas a este respecto. La exégesis de Bill Venables, mencionada anteriormente, es un gran recurso para esto. Solo di no al tipo III; se basa en una noción defectuosa de equilibrio y tiene una potencia menor debido a la tonta ponderación de las células en el caso desequilibrado.

La función de rmspaquete R proporciona una forma más natural y menos propensa a errores de obtener contrastes generales y poder describir lo que hizo contrast.rms. Los contrastes pueden ser muy complejos, pero para el usuario son muy simples porque se expresan en términos de diferencias en los valores predictivos. Se admiten pruebas y contrastes simultáneos. Esto maneja efectos de regresión no lineal, efectos de interacción no lineal, contrastes parciales, todo tipo de cosas.

Pruebe el comando Anova en la biblioteca del automóvil. Utilice el argumento type = "III", ya que el valor predeterminado es type II. Por ejemplo:

library(car)
mod <- lm(conformity ~ fcategory*partner.status, data=Moore, contrasts=list(fcategory=contr.sum, partner.status=contr.sum))
Anova(mod, type="III")

También autopromocionado, escribí una función para exactamente esto: https://github.com/samuelfranssens/type3anova

Instalar de la siguiente manera:

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

También necesitará tener el carpaquete instalado.