Anova para diseños factoriales con 2 factores (efectos fijos).

Editado el 10/5/2013 para añadir, al final,  código del ANOVA con 3 factores

Otra larga pausa sin publicar. Como en la canción de Lennon, “Life is what happens to you while you’re busy making other plans“.  Tengo varias entradas del blog y dos capítulos más de los apuntes, preparados pero que no acaban de ver la luz, porque varias cosas se cruzan por el camino. Y la entrada de hoy es, de hecho, una de esas cosas que se cruzan. Varios alumnos me han pedido ayuda en relación con un análisis de datos que involucra una tabla Anova con dos factores (y efectos fijos). En mi – ingente – lista de tareas para el blog está previsto tratar ese tema con más profundidad, a lo largo de varias entradas, y con algún artista invitado, si puede ser.  Porque las generalizaciones del ANOVA que vemos en el curso son, de hecho, una de las técnicas estadísticas más utilizadas en la práctica. Pero la realidad impone sus ritmos, así que vamos a hacer una introducción, “quick and dirty” al tema; el lector queda avisado: esta no es, ni mucho menos, la historia completa.

Lo que sigue se inspira en el Ejemplo 12A del Capítulo 12 (pág. 334) de “Biostatistical Design and Analysis Using R“, de M. Logan. A su vez, ese ejemplo se basa en los datos del Ejemplo 9.4, del Capítulo 9 (pág. 224) de “Experimental Design and Data Analysis for Biologists“, de G.P. Quinn y M.J. Keough. En este ejemplo, se usan los datos de un estudio (Quinn, G. P. (1988). “Ecology of the intertidal pulmonate limpet Siphonaria Diemenensis”, Journal of experimental marine biology and ecology, vol. 117, no2, pp. 115-136.) sobre los efectos que tienen

  1. la estación del año y
  2. la densidad de ejemplares adultos

sobre la producción de huevos en la especie Siphonaria Diemenensis.

El fichero de datos que vamos a usar es quinn.csv. Lo leemos con R (esta vez, como novedad, en lugar de descargarlo primero, lo vamos a descargar directamente desde R) y le echamos un vistazo, no es excesivamente largo. 

> fichero="http://www.fernandosansegundo.es/Estadistica/quinn.csv"
> (quinn=read.table(fichero,header=T,sep=","))
 DENSITY SEASON EGGS
1 8 spring 2.875
2 8 spring 2.625
3 8 spring 1.750
4 8 summer 2.125
5 8 summer 1.500
6 8 summer 1.875
7 15 spring 2.600
8 15 spring 1.866
9 15 spring 2.066
10 15 summer 0.867
11 15 summer 0.933
12 15 summer 1.733
13 30 spring 2.230
14 30 spring 1.466
15 30 spring 1.000
16 30 summer 1.267
17 30 summer 0.467
18 30 summer 0.700
19 45 spring 1.400
20 45 spring 1.022
21 45 spring 1.177
22 45 summer 0.711
23 45 summer 0.356
24 45 summer 0.711

Este conjunto de datos es típico de un ANOVA con dos factores (y efectos  fijos). Tenemos una variable (respuesta) cuantitativa, que es la producción de huevos (la tercera columna de los datos). Y queremos estudiar la relación de esa variable con dos variables (explicativas), que son la la densidad de adultos y la estación (primera y segunda columnas, respectivamente).  

Se trata de dos variables cualitativas o factores. Esto es especialmente evidente en el caso de la variable SEASON, que tiene dos niveles (spring y summer). Podemos confirmarlo preguntándole a R de qué tipo de objeto se trata, así:

> class(quinn$SEASON)
[1] "factor"

Por otro lado, si hacemos lo mismo con DENSITY:

> class(quinn$DENSITY)
[1] "integer"

veremos que R ha entendido esta variable como un vector de enteros. Pero nosotros queremos tratar a DENSITY como un factor con cuatro niveles, que son 8, 15, 30 y 45. Esta es una situación típica, cuando tenemos un factor cuyos niveles se representan numéricamente. Para asegurarnos de que las cosas son como queremos que sean, hacemos:

> quinn$DENSITY = factor(quinn$DENSITY)
> class(quinn$DENSITY)
[1] "factor"

y, como puede verse, DENSITY pasa a ser considerado como un factor. Ya he avisado de que no vamos a tratar el tema en profundidad, pero merece la pena subrayar que, como se explica en el texto de Quinn y Keough, estamos considerando ambos como efectos fijos (y no como efectos aleatorios). En particular, en el diseño de este experimento, esos valores se han fijado por el experimentador  (si alguien quiere repetir este experimento, volverá a utilizar los valores 8, 15, 30 y 45 de DENSITY, y no otras densidades al azar; con densidades al azar sería otro diseño experimental). Para entender con más detalle esta diferencia entre efectos fijos y aleatorios recomiendo:

(a) leer los libros citados (o, para mis alumnos de Alcalá, el capítulo 10 de “Estadística para Biología y Ciencias de la Salud”, de S.Milton, que está disponible en nuestra biblioteca)

(b) esperar a que terminemos de escribir el capítulo correspondiente de los apuntes :).

¡No es una cuestión menor! El análisis de los datos es distinto dependiendo de esta diferencia entre efectos fijos y aleatorios. Vamos a ocuparnos del caso de efectos fijos porque es el más sencillo para empezar.

El modelo que vamos a estudiar, la relación de la variable (respuesta) EGGS con las variables (explicativas) DENSITY  y SEASON, se describe así en R:

modelQuinn=lm(EGGS ~ DENSITY * SEASON,data=quinn)

Esta notación es una generalización de la notación ~ que ya vimos al tratar el análisis unifactorial.  El asterisco * que aparece entre DENSITY y SEASON le indica a R que estamos interesados en la posible interacción entre esas dos variables. Volveremos sobre esto más abajo.

Otra pregunta importante que tenemos que hacernos, al principio del análisis, es si el diseño es equilibrado (balanced). Es decir, si todas las muestras son del mismo tamaño. Debemos tener en cuenta, entre otras consideraciones, que las hipótesis que justifican la validez del ANOVA son especialmente  sensibles a la existencia de muestras de tamaños distintos (como ya mencionamos al hablar  del análisis unifactorial). En nuestro ejemplo podemos comprobar si el diseño es equilibrado así:

> replications(modelQuinn,data=quinn)
DENSITY SEASON DENSITY:SEASON
6 12 3
> (isBalanced=!is.list(replications(modelQuinn, data=quinn)))
[1] TRUE

La primera instrucción usa la función replications para obtener un vector, o una lista, de las veces que cada nivel de cada factor aparece repetido en el experimento.  De hecho,  si el experimento es no-equilibrado, la función replications devuelve como resultado una lista, no un vector. La segunda instrucción utiliza esto para preguntarle a R si el diseño es equilibrado o no. Para ver esto con algo más de detalle, en este enlace hay un fichero de instrucciones que, al ejecutarlo, descargará una versión modificada de los datos, no equilibrada. Es una buena idea ejecutar este código y ver la diferencia en los resultados de replications.

Ahora que ya sabemos que el diseño es equilibrado, construimos la tabla ANOVA:

> quinn.aov = aov(EGGS ~ SEASON * DENSITY, data = quinn)
> summary(quinn.aov)
Df Sum Sq Mean Sq F value Pr(>F)
SEASON 1 3.250 3.250 17.842 0.000645 ***
DENSITY 3 5.284 1.761 9.669 0.000704 ***
SEASON:DENSITY 3 0.165 0.055 0.301 0.823955
Residuals 16 2.915 0.182
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

El resultado (y los propios comandos que hemos utilizado) son similares a los del caso unifactorial. Hay una línea para cada factor, y una para el término de interacción (SEASON:DENSITY). Los p-valores, y los códigos de asteriscos que los acompañan, nos indican que:

  • se obtienen resultados significativos para las dos variables SEASON y DENSITY, cada una por separado. Es decir, que la variable EGGS depende, en efecto, de estas dos variables.
  • el término de interacción no es significativo: no hay interacción entre las dos variables.

Podemos ilustrar esto último gráficamente mediante la función interaction.plot de R:

interaction.plot(quinn$DENSITY, quinn$SEASON, quinn$EGGS)

que produce este resultado:ANOVA2-InteractionPlot

El factor DENSITY se muestra en el eje horizontal, la respuesta EGGS en el eje vertical, y cada uno de los dos factores de SEASON se muestra como una línea que conecta los correspondientes valores. En un caso ideal, la ausencia completa de interacción correspondería a dos líneas perfectamente paralelas en este gráfico (y la existencia de interacción es evidente cuando las líneas se cruzan). En la práctica, el paralelismo en general está lejos de ser perfecto. No obstante, la gráfica de este ejemplo sí parece indicar que no existe interacción entre ambos factores. Y en cualquier caso, tenemos los resultados de la tabla ANOVA para corroborarlo.

El ejemplo 12B del libro de Logan ilustra otro experimento de Quinn en el que, para las mismas variables, de hecho, existe interacción. Para no repetir toda la discusión, vamos a enlazar un fichero de código R que corresponde al mismo análisis, pero para ese segundo conjunto de datos.  En este caso la tabla ANOVA es:

Df Sum Sq Mean Sq F value Pr(>F)
SEASON 1 17.148 17.148 119.845 1.34e-07 ***
DENSITY 2 4.002 2.001 13.984 0.000732 ***
SEASON:DENSITY 2 1.691 0.845 5.908 0.016363 *
Residuals 12 1.717 0.143

que muestra que el término de interacción es, en este caso, significativo. La gráfica obtenida con interaction.plot es esta: ANOVA2-InteractionPlot-1

en la que la falta de paralelismo entre ambas líneas resulta bastante más evidente.

Volviendo al ejemplo 12A, otra representación gráfica interesante se consigue utilizando el paquete effects. Sin entrar en detalles, el código para este ejemplo sería:

library(effects)
Quinn.eff=effect("DENSITY*SEASON",modelQuinn,confidence.level=0.90)
plot(Quinn.eff)

y el resultado es:

ANOVA2-effectPlotSe trata esencialmente de otra manera de presentar los resultados que ya obtuvimos en el interaction.plot, junto con intervalos de confianza para la media en cada uno de los niveles del factor DENSITY.

Añadido el 10/5/2013: varios de mis alumnos me han explicado que en sus problemas aparecen tres factores. No puedo entrar aquí en muchos detalles, que quedan para una  futura entrada del blog, pero voy a dejar enlazado un fichero de código R (que carga los datos necesarios) correspondiente a un ejemple de este tipo de análisis (ANOVA, 3 factores, efectos fijos, diseño equilibrado). El problema consiste en estudiar el rendimiento de tres variedades de centeno (I,II,III), con dos fertilizantes (A y B), y dos niveles de agua

La tabla de datos original del problema era esta

I II III
A B A B A B
Bajo 18.6 10.2 19.4 14.7 15.9 20.9
18.8 10.5 18.9 17.8 16.5 21.0
15.8 8.5 18.4 15.6 17.2 21.3
17.9 10.6 18.9 16.5 16.5 20.4
Alto 12.3 12.3 13.2 10.2 19.4 13.2
10.9 10.7 15.5 8.5 21.2 13.5
10.7 12.0 11.0 7.1 20.4 12.8
11.8 11.8 11.8 6.5 20.3 15.2

Y una de las primeras tareas del código R consiste en convertir estos datos en un data.frame manejable desde R. El problema procede de Statgraphics Plus 4, de J.Llovet, D. Delgado, J. Martínez. Pág. 237. Ed. Anaya Multimedia, ISBN 9788441509504. Es muy posible que haya que instalar en R algunas librerías adicionales para poder ejecutar el código completo. Hay instrucciones de cómo hacer eso en esta entrada.

Gracias por la atención.

Anuncios

Responder

Introduce tus datos o haz clic en un icono para iniciar sesión:

Logo de WordPress.com

Estás comentando usando tu cuenta de WordPress.com. Cerrar sesión /  Cambiar )

Google+ photo

Estás comentando usando tu cuenta de Google+. Cerrar sesión /  Cambiar )

Imagen de Twitter

Estás comentando usando tu cuenta de Twitter. Cerrar sesión /  Cambiar )

Foto de Facebook

Estás comentando usando tu cuenta de Facebook. Cerrar sesión /  Cambiar )

Conectando a %s