14.4: Ejecución de un ANOVA en R

Estoy bastante seguro de que sé lo que estás pensando después de leer la última sección, sobre todo si seguiste mis consejos e intentaste escribir todos los comandos en ti mismo... hacer los cálculos ANOVA tú mismo apesta. Hay bastantes cálculos que necesitábamos hacer en el camino, y sería tedioso tener que hacer esto una y otra vez cada vez que quisieras hacer un ANOVA. Una posible solución al problema sería tomar todos estos cálculos y convertirlos tú mismo en algunas funciones R. Aún tendrías que escribir mucho, pero al menos solo tendrías que hacerlo una vez: una vez que hayas creado las funciones, puedes reutilizarlas una y otra vez. No obstante, escribir tus propias funciones es mucho trabajo, por lo que este es una especie de último recurso. Además, es mucho mejor si alguien más hace todo el trabajo por ti...

Usar la función aov () para especificar su ANOVA

Para hacerte la vida más fácil, R proporciona una función llamada aov (), que —obviamente— es un acrónimo de “Análisis de varianza”. ²⁰⁵ ¿Si escribes? aov y echa un vistazo a la documentación de ayuda, verás que hay varios argumentos a la función aov (), pero los únicos dos que nos interesan son fórmula y datos. Como hemos visto en algunos lugares anteriormente, el argumento de fórmula es lo que usas para especificar la variable de resultado y la variable de agrupación, y el argumento data es lo que usas para especificar el marco de datos que almacena estas variables. En otras palabras, para hacer el mismo ANOVA que calculé laboriosamente en la sección anterior, usaría un comando como este:

aov( formula = mood.gain ~ drug, data = clin.trial ) 

En realidad, esa no es toda la historia, como verás en cuanto mires la salida de este comando, que por el momento he ocultado para no confundirte. Antes de entrar en detalles, debo señalar que cualquiera de estos comandos hará lo mismo:

aov( clin.trial$mood.gain ~ clin.trial$drug ) 
aov( mood.gain ~ drug, clin.trial ) 

En el primer comando, no especificé un conjunto de datos, y en su lugar confié en el operador $ para decirle a R cómo encontrar las variables. En el segundo comando, dejé caer los nombres de los argumentos, lo cual está bien en este caso porque formula es el primer argumento de la función aov (), y data es el segundo. Independientemente de cómo especifique el ANOVA, puedo asignar la salida de la función aov () a una variable, como esta por ejemplo:

my.anova <- aov( mood.gain ~ drug, clin.trial ) 

Esto casi siempre es algo bueno que hacer, porque hay muchas cosas útiles que podemos hacer con la variable my.anova. Entonces supongamos que es este último comando el que utilicé para especificar el ANOVA que estoy tratando de ejecutar, y como consecuencia tengo esta variable my.anova sentada en mi espacio de trabajo, esperando que haga algo con ella...

Comprender lo que produce la función aov ()

Ahora que hemos visto cómo usar la función aov () para crear my.anova, será mejor que echemos un vistazo a lo que realmente es esta variable. Lo primero que hay que hacer es verificar para ver qué clase de variable hemos creado, ya que es algo interesante en este caso. Cuando hacemos eso...

class( my.anova )

## [1] "aov" "lm"

... descubrimos que my.anova en realidad tiene dos clases! La primera clase nos dice que es un objeto aov (análisis de varianza), pero la segunda nos dice que también es un objeto lm (modelo lineal). Más adelante, veremos que esto refleja una relación estadística bastante profunda entre ANOVA y regresión (Capítulo 15 y significa que cualquier función que exista en R para lidiar con regresiones también se puede aplicar a objetos aov, lo cual es limpio; pero me estoy adelantando a mí mismo. Por ahora, quiero señalar que lo que hemos creado es un objeto aov, y también hacer el punto de que los objetos aov son en realidad bestias bastante complicadas. No voy a tratar de explicar todo sobre ellos, ya que va mucho más allá del alcance de un tema introductorio de estadística, sino para darte una pequeña pista de algunas de las cosas que R almacena dentro de un objeto aov, vamos a pedirle que imprima los nombres () de todos las cantidades almacenadas...

names( my.anova )

##  [1] "coefficients"  "residuals"     "effects"       "rank"         
##  [5] "fitted.values" "assign"        "qr"            "df.residual"  
##  [9] "contrasts"     "xlevels"       "call"          "terms"        
## [13] "model"

A medida que repasemos el resto del libro, espero que algunos de estos se vuelvan un poco más obvios para ti, pero ahora mismo eso va a parecer bastante opaco, maldito. Eso está bien. No necesitas conocer ninguno de los detalles al respecto en este momento, y la mayor parte no necesitas en absoluto... lo que sí necesitas entender es que la función aov () hace muchos cálculos por ti, no solo los básicos que delineé en las secciones anteriores. Lo que esto significa es que generalmente es una buena idea crear una variable como my.anova que almacene la salida de la función aov ()... porque más adelante, puedes usar my.anova como entrada a muchas otras funciones: esas otras funciones pueden sacar bits y piezas del aov, y calcula varias otras cosas que podrías necesitar.

En ese momento. Lo más sencillo que puedes hacer con un objeto aov es imprimirlo (). Cuando hacemos eso, nos muestra algunas de las cantidades clave de interés:

print( my.anova )

## Call:
##    aov(formula = mood.gain ~ drug, data = clin.trial)
## 
## Terms:
##                     drug Residuals
## Sum of Squares  3.453333  1.391667
## Deg. of Freedom        2        15
## 
## Residual standard error: 0.3045944
## Estimated effects may be unbalanced

De manera específica, imprime un recordatorio del comando que usaste cuando llamaste a aov () en primer lugar, te muestra las sumas de valores de cuadrados, los grados de libertad, y un par de otras cantidades que no nos interesan realmente en este momento. Observe, sin embargo, que R no usa los nombres “entre-grupo” y “dentro del grupo”. En cambio, trata de asignar nombres más significativos: en nuestro ejemplo particular, la varianza entre grupos corresponde al efecto que tiene el fármaco sobre la variable de resultado; y la varianza dentro de los grupos corresponde a la “sobrante” variabilidad, por lo que llama a eso los residuales. Si comparamos estos números con los números que calculé a mano en la Sección 14.2.5, se puede ver que son idénticos... las sumas entre grupos de cuadrados es SS _b =3.45, las sumas de cuadrados dentro de grupos es SS _w =1.39, y los grados de libertad son 2 y 15 repectivamente.

Ejecución de las pruebas de hipótesis para el ANOVA

Bien, entonces hemos verificado que my.anova parece estar almacenando un montón de los números que estamos buscando, pero la función print () no nos dio el resultado que realmente queríamos. ¿Dónde está el valor F? ¿El valor p? Estos son los números más importantes en nuestra prueba de hipótesis, pero la función print () no los proporciona. Para obtener esos números, necesitamos usar una función diferente. En lugar de pedirle a R que imprima () el objeto aov, deberíamos haber pedido un summary () del mismo. ²⁰⁶ Cuando hacemos eso...

summary( my.anova )

##             Df Sum Sq Mean Sq F value   Pr(>F)    
## drug         2  3.453  1.7267   18.61 8.65e-05 ***
## Residuals   15  1.392  0.0928                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

... obtenemos todos los números clave que calculamos antes. Obtenemos las sumas de cuadrados, los grados de libertad, los cuadrados medios, el estadístico F y el valor p en sí. Todos estos son idénticos a los números que calculamos nosotros mismos al hacerlo de la manera larga y tediosa, e incluso se organiza en el mismo tipo de tabla ANOVA que mostré en la Tabla 14.1, y luego rellenada a mano en la Sección 14.2.5. Lo único que es incluso ligeramente diferente es que algunos de los nombres de filas y columnas son un poco diferentes.