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.

Actualizar R en Windows

Y una entrada mínima para complementar la información de la anterior, esta vez sobre la actualización en Windows. Las instrucciones, en inglés, están aquí.

Una advertencia: si, por ejemplo, has instalado R y R-Commander usando el paquete R-UCA de la universidad de Cadiz (como expliqué en esta entrada), entonces es posible que te hayas acostumbrado a comunicarte con R a través de R-Commander, y entonces apenás usarás la ventana RGui. En tal caso, y para evitar problemas, no ejecutes la instalación desde R-Commander.  Es mejor que cierres la ventana de R-Commander y ejecutes la instalación desde esa ventana RGui.

Actualizar R a la última versión en Ubuntu 12.04

Una entrada breve, y de nuevo al margen de la serie principal de entradas del blog. Si no usas Ubuntu, huye de aquí.

Trabajo la mitad del tiempo en máquinas Linux, de hecho, desde hace años, básicamente en Ubuntu. Buena parte del software científico que uso se instala y funciona con más naturalidad en Linux. A raíz del asunto de los ficheros csv grandes en R de la anterior entrada. y tratando de reproducir esas operaciones en Ubuntu 12.04 (Precise), he descubierto que la versión de R que se incluye en los repositorios de Ubuntu es la 2.14.1-1, mientras que la última versión de R disponible es la 2.15.3. No parece una gran diferencia, hasta que tratas de instalar una librería que no está disponible en la versión más antigua (por ejemplo la librería sqldf que necesitaba para el asunto de los csv grandes). Así que he tenido que buscar una solución para actualizar R en Ubuntu, y esta entrada sirve tan sólo para apuntar la dirección de la solución, que es esta (de nuevo StackOverflow al rescate).

Ficheros csv grandes en R

Una advertencia preliminar: esta entrada se ha escrito para describir un problema concreto con R (el uso de ficheros de datos csv especialmente grandes), y por lo tanto está al margen de la serie principal de entradas del blog, que pretenden ir describiendo progresivamente el uso de R.  Si tu interés principal es ese, aprender poco a poco a usar R, puedes saltarte esta entrada y, en todo caso, regresar a ella cuando se te plantee el problema.

Size matters?

Aunque R es el protagonista de este blog, eso no impide que mantengamos un ojo atento a otras plataformas de software estadístico. Por esa razón, para mejorar mi conocimiento   de una de las alternativas más populares, estoy inscrito en el curso Passion Driven Statistics de la plataforma Coursera. Se trata de un curso básico de Estadística, que utiliza -primordialmente- el programa SAS como herramienta (más concretamente, SAS OnDemand for Academics).

Como me  siento mucho más cómodo trabajando con R, antes de abordar una tarea del curso con SAS, a menudo hago primero las cuentas necesarias con R. De esa manera puedo saber fácilmente si los resultados que obtengo son correctos. Esta semana,  trabajando en una de esas tareas del curso, me he topado con el problema de tener que utilizar  en R un fichero de datos (en formato csv) relativamente grande (15.5 Mb aprox.). De hecho, una de las ventajas  que a menudo se citan, cuando se compara SAS con otros programas, es su capacidad para trabajar con grandes ficheros de datos. Y, de hecho, el primer intento de cargar el fichero en R no funcionó correctamente. Afortunadamente, como sucede a menudo con R, la solución estaba a la vuelta de…  unas pocas esquinas. Pero no adelantemos acontecimientos, y vamos paso a paso. 

El fichero csv de marras es una base de datos sobre cráteres (¡más de trescientos mil!) en la superficie del planeta Marte, procedente de la tesis doctoral (University of Colorado) de  Stuart  Robbins (aunque no es relevante para lo que aquí nos ocupa, el fichero se describe con más detalle en este fichero).

Manos a la obra con el primer intento. Descargo el fichero al Escritorio de la máquina en la que voy a hacer las cuentas (s. operativo Windows 7), arranco R y ejecuto este código, que crea un data frame (llamado mars) de R a partir del fichero csv, y a continuación comprueba el número de filas de ese data frame para verificar que hemos leído tantos datos como esperábamos. Se muestra la salida de los comandos:

> setwd("C:/Users/Fernando/Desktop")
> mars=read.table(file="marscrater_pds.csv",sep=",",header=TRUE)
> nrow(mars)
[1] 169587

Ahora es evidente que algo ha ido mal, porque como hemos dicho el fichero contiene datos de más de trescientos mil cráteres. Mi primera impresión, al ver el resultado, fue que el fichero se había truncado en el proceso de lectura. Y lo peor es que no había aparecido ningún mensaje de advertencia…

La librería sqldf

A estas alturas del problema,  si eres, como yo, un poco maniático (¿un poco, en serio?), te habrás lanzado ya sin frenos por la pendiente de los ficheros de gran tamaño en R. La mala noticia es que me lo hubiera podido ahorrar, porque en realidad el problema se debe a que estaba usando la función read.table de forma incorrecta. La buena noticia es que, por el camino, doblando varias de esas esquinas de las que hablaba, he descubierto una librería de R que permite trabajar con ficheros verdaderamente grandes. Y en este caso no me refiero a unos pocos megas, sino a un fichero csv de ¡cerca de un gigabyte de datos!

La fuente de inspiración de lo que sigue es esta discusión de StackOverflow, que se ha consolidado para mi como uno de los oráculos más fiables de Internet. Para empezar bien las cosas vamos a usar el propio R para crear un fichero de datos verdaderamente enorme. He tenido  que modificar un poco el código de StackOverflow, porque estoy usando una modesta versión de 32 bits de R en una modesta máquina Windows 7, en la que R se quedaba sin memoria para ejecutar directamente el código original. Antes de ejecutar este código asegúrate de que tienes espacio en tu máquina para alojar un fichero de aprox. un gigabyte:

bigdf = data.frame(dim=sample(letters, replace=T, 2e7), fact1=rnorm(2e7), fact2=rnorm(2e7, 20, 50))
write.csv(bigdf, 'bigdf.csv', quote = F)

En mi máquina R tarda unos minutos (<10) en completar estos comandos, y el resultado,  en mi escritorio,  es un enorme fichero csv de 20 millones de líneas. Ahora ha llegado el momento de leer este fichero desde R. Porque una cosa es escribirlo, y otra bien distinta es cargarlo en memoria y poder hacer operaciones con él.  La clave, como explicaban en StackOverflow, es instalar la librería sqldf de R (en la séptima sesión con Rcmdr explicamos como instalar librerías adicionales).

Para cargar el  fichero, reiniciamos  R,  instalamos esa librería, y ejecutamos este código:

setwd("C:/Users/Fernando/Desktop")
library(sqldf)
f =file("bigdf.csv")
bigdf = sqldf("select * from f", dbname = tempfile(), file.format =list(header = T, row.names = F))
nrow(bigdf)

En mi máquina ha tardado menos de diez minutos, lo suficiente para estirar las piernas. Y al final la función nrow(bigdf) ha contestado, como se esperaba, 20000000. Y ahí tenemos un hermoso data frame con veinte millones de filas. Es bastante impresionante ejecutar comandos sobre esos datos y ver como R apenas acusa el esfuerzo. Por ejemplo, he probado a hacer  estas cosas:

> table(bigdf$dim)
a b c d e f g h i j k l m n o p q r s t u v w x y z
768523 770070 770127 768994 769183 770936 768517 769162 769874 768586 769646 768160 767658 768757 768536 770234 769834 769125 768700 769659 768125 769524 769147 768973 770505 769445
> mean(bigdf$fact1)
[1] 0.0007351861
> sd(bigdf$fact1)
[1] 0.9999446
> hist(bigdf$fact1)

Y todas se han ejecutado casi instantáneamente. El resultado del histograma es el que se esperaba, claro:HistogramaBigData

pero lo importante es verlo aparecer casi en el acto, teniendo en cuenta que representa veinte millones de observaciones…

De vuelta a Marte

Naturalmente, cuando se usa la librería sqldf, el fichero de datos de los cráteres de Marte se carga sin problemas. Así que “problema resuelto”. ¿O no? Me seguía rondando la cabeza la idea de que, en la presunta lectura truncada con read.table, no había habido ningún mensaje de advertencia (¿he dicho ya lo de que soy un poco maniático?).  Así que un poco más de búsqueda por la red (ver, por ejemplo, esta discusión) me ha servido para descubrir que, en realidad, lo que sucedía es que estaba usando mal la función read.table, que es muy capaz por si misma, de leer ficheros bastante grandes, de varios cientos de miles de filas. La clave, cuando como en mi caso se usa un separador como

sep=","

es incluir la opción quote=”” en la función read.table. De modo que el segundo intento de usar read.table es este (se muestra la salida):

</pre>
> rm(list=ls())
<span style="line-height: 1.5;">> setwd("C:/Users/Fernando/Desktop")</span>
<span style="line-height: 1.5;">> mars=read.table(file="marscrater_pds.csv",sep=",",header=TRUE,quote="")</span>
<span style="line-height: 1.5;">> nrow(mars)</span>
[1] 384343
<span style="line-height: 1.5;">

Y, como puede verse, ahora las cosas funcionan correctamente, La razón por la que es necesario usar quote=”” es bastante cabalística, y tiene que ver con la presencia de un número muy reducido de comillas simples ( es decir, el carácter ‘ ) en fichero csv de cráteres de Marte. Esas comillas afectan a la forma en que R lee el fichero y le confunden sobre donde empiezan y terminan las líneas. Así que, desde el punto de vista de R, la primera lectura del fichero había sido un éxito. Y eso explica la  ausencia de mensajes de error…

Es un comportamiento muy extraño, que nos recuerda la dificultad intrínseca que tienen siempre estas operaciones de lectura y/o conversión de datos entre distintos formatos. De hecho, la motivación para escribir esta entrada es, básicamente, para que me sirva (a mi, o a otros) de referencia en el futuro, porque estoy seguro de que tardaré en tropezarme otra vez con este dichoso problema de las comillas simples. Al menos, por el camino, hemos aprendido a cargar ficheros enormes en R.

Llevo muy poco tiempo trabajando con SAS, como para formarme una opinión definitiva. Pero me alegra, habiendo apostado por R desde hace tiempo, comprobar que los ficheros grandes tampoco suponen un problema. Y en cuanto a otro de los supuestos inconvenientes de R, su pronunciada “curva de aprendizaje inicial”… en fin, ya digo que tengo muy poca experiencia con SAS. Pero estoy empezando a entender porque alguien lo resumió diciendo “SAS=semicolon, always semicolon”. Me da la sensación de que La sintaxis de SAS es, al menos, tan poco intuitiva como pueda serlo la de R.

Muchas gracias por la atención.

Varias cosas colaterales: tasa de paro personal y libro sobre Numb3rs

Hasta la entrada de ayer, más de un mes sin publicar. El final de la temporada de exámenes, como siempre, en lugar de ser un alivio nos enfrenta a todas las tareas urgentes que hemos retrasado para ocuparnos de las inaplazables. Afortunadamente las cosas se van encauzando, y espero que eso se traduzca en una actividad más regular del blog. Tengo muchas ideas en cartera, así que temas no nos faltarán.

Tasa de paro personalizada

De momento, una entrada sobre uno de los temas colaterales, la visualización de datos. He descubierto hace poco esta herramienta TTP: Tu Tasa de Paro, que permite situarse a uno mismo en relación con los datos de la EPA. Llegué a ella buscando información, después de algunas conversaciones sobre la reforma de la Universidad,  motivadas, entre otras cosas, por la reciente publicación del informe de la Comisión de Expertos para la Reforma del Sistema Universitario (nota de prensa del ministerio), y la evidencia – anecdótica, no tengo todavía los datos  – de que las titulaciones tradicionales de Ciencias y muchas Ingenierías se están “despoblando”.

Libro sobre Numb3rs

Mis (sufridos) alumnos me han visto emplear en más de una ocasión algún vídeo de la serie Numb3rs para ilustrar algún tema de la clase. En particular, este fragmento me gusta mucho como ilustración del Problema de Monty Hall. Precisamente ahora, estoy terminando de leer el libro The numbers behind Numb3rs: Solving Crime with Mathematics.  El libro es, en general, un complemento interesante. a la serie. Y me llama la atención la presencia prominente de la Estadística y en general el análisis de datos, como uno de los temas centrales.

Contrastes de hipótesis y primeros pasos en la programación: estructuras If-else.

Contraste de hipótesis sobre igualdad de varianzas

Aunque esta entrada se presenta como continuación de la 14ª sesión con Rcmdr, en la que hablábamos de contrastes de hipótesis, en realidad es más bien una excusa para empezar a hablar en el blog de las estructuras de  programación de R. Si ya tienes alguna experiencia escribiendo programas (en cualquier lenguaje de programación), la lectura de esta entrada te resultará muy, muy fácil. Pero si nunca lo has hecho, entonces esta entrada es para ti, he tratado de hacerla lo más fácil de digerir posible.

En aquella entrada hicimos algún contraste de hipótesis usando, por un lado, funciones de R específicas para ese tarea (como t.test), y por otro lado, haciendo “a mano” el contraste. Es decir, calculando los valores necesarios para el contraste paso a paso. Y al final de aquella entrada dejamos sin cubrir el caso del contraste de igualdad de varianzas en dos poblaciones normales. En parte para no hacerla demasiado extensa, y en parte porque quería reservar ese caso para esta entrada. Vamos con ello, entonces. En uno de los últimos exámenes que hemos hecho, pedíamos a nuestros alumnos que hicieran este ejercicio:

Dos marcas de analgésicos a firman en su publicidad que sus productos proporcionan un
alivio muy rapido del dolor de cabeza. Para decidir cual de los dos es más rápido, hemos
seleccionado un grupo de 22 pacientes y los hemos dividido aleatoriamente en dos grupos. Cada grupo ha recibido tratamiento con uno de los dos analgésicos. Los tiempos, en minutos, que los pacientes de cada grupo tardaron en sentir aliviado su dolor fueron estos:
Grupo 1={ 9.9, 8.0, 9.5, 5.9, 12.2, 13.5, 9.6, 11.5, 12.5, 9.5, 10.3, 11.9 }
Grupo 2={ 8.2,17.3, 10.1, 10.2, 9.1, 10.5, 9.7, 9.0, 15.2, 11.6 }
Utiliza estos datos para decidir cuál de los dos medicamentos parece producir efecto más
rápidamente. Ten en cuenta que no se sabe nada sobre la posible igualdad de varianzas en las poblaciones involucradas. Utiliza un nivel de confi anza del 95% en todos los contrastes de hipótesis que realices, y calcula además el p-valor de esos contrastes.

Es un ejercicio de contraste de medias en dos poblaciones (que se asumen normales), pero, como indica la frase subrayada, y dado que se trata de muestras pequeñas, es conveniente llevar a cabo previamente un contraste de igualdad de varianzas, con hipótesis nula:

H_0=\{\sigma_1^2=\sigma_2^2\}

Que es precisamente el caso que habíamos dejado pendiente.

El estadístico que vamos a utilizar en este caso es el cociente de las cuasivarianzas muestrales:

\dfrac{s_1^2}{s_2^2}

Utilizamos R para calcular el valor del estadístico en las muestras de ambas poblaciones:

> muestra1=c(9.9, 8.0, 9.5, 5.9, 12.2, 13.5, 9.6, 11.5, 12.5, 9.5, 10.3, 11.9)
> muestra2=c(8.2,17.3, 10.1, 10.2, 9.1, 10.5, 9.7, 9.0, 15.2, 11.6)

> (s1=sd(muestra1))
[1] 2.115509
> (s2=sd(muestra2))
[1] 2.914504

> (n1= length(muestra1))
[1] 12
> (n2= length(muestra2))
[1] 10

> (Estadistico=s1^2/s2^2)
[1] 0.5268664

Si la hipótesis nula fuera cierta, y las varianzas de las dos poblaciones fuesen iguales,  la distribución muestral de este estadístico sería una F de Fisher.

Así que para calcular el p-valor del contraste basta con hacer

> (pValor=2*(pf(Estadistico,df1=n1-1,df2=n2-1)))
[1] 0.3139952

Las preguntas estándar que se hace cualquier aprendiz de Estadística (y que tienen una contestación más extensa en los apuntes de teoría) son:  ¿por qué el 2? y ¿por qué usamos pf en lugar de 1-pf? Multiplicamos por dos porque se trata de un contraste bilateral, y debemos tener en cuenta la probabilidad de las dos colas de la distribución. Y usamos pf porque estamos calculando la probabilidad correspondiente a la cola izquierda de la distribución. Y ese es el punto clave (que hace de este el ejemplo ideal para esta entrada, como veremos) y que muchas veces causa dificultades a los principiantes. ¿Cómo sabemos que hay que usar la cola izquierda y no la cola derecha? Una gráfica de la distribución de Fisher nos ayudará en la discusión:2012_13_EstadisticaBiologia-ExamenEnero-Ejercicio01-ContrasteSi la hipótesis nula fuera cierta, los valores del estadístico cercanos a uno serían los más probables (porque el estadístico vale uno cuando las cuasivarianzas son iguales). Los valores que ponen en aprietos a la hipótesis nula son los que se alejan de uno, hacia una u otra de las colas de la distribución. Y como nosotros hemos obtenido un estadístico aproximadamente igual a 0.527, que es menor que uno, debemos usar la cola izquierda. Para tratar de dejarlo aún más claro, si intercambiamos los nombres de las muestras (como si las hubiéramos ordenado al revés), tendríamos:

> muestra2=c(9.9, 8.0, 9.5, 5.9, 12.2, 13.5, 9.6, 11.5, 12.5, 9.5, 10.3, 11.9)
> muestra1=c(8.2,17.3, 10.1, 10.2, 9.1, 10.5, 9.7, 9.0, 15.2, 11.6)

> (s1=sd(muestra1))
[1] 2.914504
> (s2=sd(muestra2))
[1] 2.115509

> (n1= length(muestra1))
[1] 10
> (n2= length(muestra2))
[1] 12

> (Estadistico=s1^2/s2^2)
[1] 1.898014

Y naturalmente, \dfrac{1}{1.898014}=0.5268665. Y si usamos  esta segunda ordenación de las muestras, puesto que ahora el valor del estadístico es mayor que uno, tenemos que usar la cola derecha de la distribución, como ilustra esta figura (no es la misma F de antes; los grados de libertad también han intercambiado sus posiciones):2012_13_EstadisticaBiologia-ExamenEnero-Ejercicio01-Contraste-aY ahora usamos 1-pf  para calcular el p-valor. En cualquier caso el resultado es, naturalmente,el mismo:

> (pValor=2*(1-pf(Estadistico,df1=n1-1,df2=n2-1)))
[1] 0.3139952

Y ahora que hemos discutido esto, vamos a lo que más nos interesa en esta entrada. Hemos visto que, dependiendo del valor del estadístico (es decir, dependiendo de cómo se ordenen las muestras), tendremos que decidir si usar la cola derecha o la cola izquierda para obtener el p-valor.  ¿Podemos automatizar este proceso de decisión? En el fondo, de eso trata la programación: de automatizar las decisiones. Por eso en el próximo apartado vamos a dar nuestro primer paso en el mundo de la programación en R.

Tomando decisiones. Estructuras If-else.

Si el objetivo (confesable) de todas las sesiones precedentes era acompañar al lector que está usando R en un curso de introducción a la Estadística, a partir de esta entrada vamos a  adentrarnos en un terreno que excede ese objetivo (tan modesto, por otra parte).  Si R ha llegado a la posición que ahora ocupa, y si su radio de acción no deja de crecer, es sobre todo porque no se limita a ser un programa de Estadística más al uso, con cada vez más opciones en los menús. Ante todo, R es un lenguaje de programación, con un fuerte sesgo -desde luego- hacia el análisis de datos. Para sacarle todo el partido a R hay que aprender algo de programación. Y otro de los objetivos de este blog es allanar ese camino a quien quiera hacerlo.

Sin más preámbulos, ¿cómo se hace esto? La decisión más básica consiste siempre en elegir entre dos caminos posibles (en general, dos opciones; pero la imagen del camino ayuda a pensar).  Y el esquema de una decisión básica siempre es este:

Si se cumple cierta condición, entonces vamos por el camino A. Y si no se cumple, vamos por el camino B.

Las frases que hemos marcado en azul son las que cambiarán en cada caso concreto, porque en cada ejemplo la condición será distinta, y distintos serán también los caminos A y B. Pero la forma de la frase es la misma siempre. Para traducir una frase de esta forma  al lenguaje de R disponemos de un estructura especial que, esquemáticamente, dejando todavía sin traducir las partes en azul de la frase, funciona así:

if(se sumple cierta condicion){
  vamos por el camino A
}else{
  vamos por el camino B
}

Veamos un ejemplo de decisión más concreto. Vamos a lanzar un dado (usando la función sample de R).  Si el resultado es menor o igual que  tres, entonces tú ganas. Y si eso no se cumple (es decir, si es mayor que tres),  yo gano. He coloreado en azul las partes correspondientes para resaltar la estructura de la frase. En R esto se convierte en:

(dado=sample(1:6,1))

if(dado<=3){
decision='tu ganas'
}else{
decision='yo gano'
}

decision

La parte interesante es el grupo central de instrucciones; el fragmento de código empieza con el lanzamiento del dado (y se muestra el resultado para ayudar a entender lo que sucede a continuación). Y el código termina mostrando el resultado, la decisión que se ha tomado.  Si es la primera vez que ves una de estas estructuras if-else, te conviene ejecutar el código varias veces para hacerte una idea de como funciona.

Volvamos al contraste sobre igualdad de varianzas. ¿Cuál es la estructura de la decisión que teníamos que tomar? Es decir, ¿cuál es la condición, y cuáles son los dos caminos? Vamos a intentar escribirlo de forma que la estructura sea evidente:

Si el estadístico es menor que uno, entonces calculamos el p-valor usando la cola izquierda de la F. Y si eso no se cumple (si el estadístico es >=1), entonces calculamos el p-valor usando la cola derecha.

Traducido a R el código que continúa con el ejemplo del examen es así (repetimos la sentencia en la que se calcula el Estadísitico para que se vea cómo encaja  esto con lo anterior):

(Estadistico=s1^2/s2^2)

if(Estadistico<1){
  pValor=2*pf(Estadistico,df1=n1-1,df2=n2-1)
}else{
  pValor=2*(pf(Estadistico,df1=n1-1,df2=n2-1))
)

pValor

Al ejecutar este código se obtiene el resultado correcto ¡sea cual sea el orden en que se coloquen las muestras! Esto nos permite escribir un fichero de comandos R que calcula el p-valor correcto, sin tener que preocuparnos del orden en que se colocan las muestras. Es muy razonable decir que es un programa en lenguaje R. Muy sencillo, pero un programa al fin y al cabo.

En próximas entradas vamos a aprender más sobre estos grupos if-else, y vamos a presentar más estructuras de programación,  que nos van a permitir crear programas en R cada vez más sofisticados, que por un lado nos permitirán sacar el mayor partido posible del sistema, y por otro lado nos harán capaces de entender los programas que otros han escrito y tal vez adaptarlos a nuestras necesidades concretas.

En la próxima entrada de esta serie, y antes de avanzar más, vamos a asegurarnos de que el terreno que pisamos es firme. Para ello, nos entretener en la forma de escribir la condición que usamos en una estructura de tipo if-else.  ¿Qué significa exactamente una condición? ¿Y como se escriben en R? Pronto lo veremos.

Muchas gracias por la atención.

Correlación vs. causalidad…

Por aquí estamos todos, alumnos y profesores, sumergidos en plena temporada de exámenes. Así que me ha parecido que podría venir bien interrumpir un par de minutos el trabajo (en mi caso, una montaña de exámenes por corregir), y relajarnos con esta viñeta

:

Me atrevo a traducir, por si alguien se lía con el inglés:

1. “Yo antes creía que la correlación implicaba causalidad.”

2. “Pero entonces hice un curso de Estadística. Y ahora ya no lo creo.”

3. “Parece que el curso sirvió para algo”. “Bueno, quizás.”

Si no conocéis xkcd.com, os lo recomiendo. Humor inteligente, casi siempre con trasfondo científico (el tipo de viñetas que lee Sheldon Cooper, no se las vayáis a recomendar a cualquiera…)  A mis estudiantes de Biología tal vez les haga gracia esta que apareció hace poco.

15ª sesión con Rcmdr: Factores y ANOVA unifactorial en R.

(Modificado el 20/01/2013 para corregir un error en la primera aparición de tapply)

Factores

En la octava sesión (se recomienda una relectura rápida de esa entrada antes de leer esta) hablamos de los datos de tipo character en R, que podían usarse, por ejemplo, para representar valores de variables cualitativas. Una situación frecuente en Estadística es esta: hemos medido una serie de valores de una variable cuantitativa X (es decir, los valores de X son números), pero esos valores aparecen agrupados de manera natural. Por ejemplo, cuando estamos midiendo la respuesta X (un número) de una serie de pacientes frente a varios tratamientos, es evidente que lo natural es agrupar los resultados según el tratamiento empleado. Una manera de hacer esto, en R, es usar un data.frame (ver la octava sesión) que contenga los valores de X, junto con los valores de una variable que indique el tipo de tratamiento empleado. Vamos a llamar

T_1, T_2, \ldots, T_k

a los distintos tratamientos. La terminología estándar en Estadística es decir que el tratamiento T es un factor, y que sus valores son los niveles del factor. Podríamos entonces pensar en recoger los resultados en un fichero de datos como este, llamado experimento-01.csv. y trabajar así, usando variables de tipo character para los factores. Pero eso tendría un impacto negativo en el rendimiento de R. Y como los factores son omnipresentes en Estadística, R ha optado por incluir un tipo especial de vector para representar los factores y sus niveles. Y, de hecho, cuando leemos un fichero como el anterior usando read.table (de nuevo, remitimos al lector a la octava sesión). el comportamiento por defecto de R es convertir las variables cualitativas en factores.

Vamos a ver esto en funcionamiento. Guardaremos el contenido del fichero en un data.frame llamado experimento, con dos campos:

  1. experimento$Respuesta, de tipo numérico que almacena los valores de la variable X (la respuesta individual de cada uno de los pacientes).
  2. experimento$Tratamiento, de tipo factor que almacena los valores de la variable T (el tipo de tratamiento empleado).

Usaremos la función read.table para leer (y mostrar) el fichero, y luego la función class para ver el tipo de objeto que corresponde a cada uno de los campos del data.frame:

(experimento=read.table(file="experimento-01.csv",header =T, sep=" "))
class(experimento$Respuesta)
class(experimento$Tratamiento)

La ventana de resultados de Rcmdr, aparte de mostrar el contenido del fichero, indica que el campo experimento$Tratamiento es, en efecto, un factor.

¿Cómo se trabaja con factores en R? La contestación no puede ser completa, porque los factores son un ingrediente clave del lenguaje de R, e intervienen en muchísimas construcciones del sistema. Pero, en cualquier caso, podemos empezar por lo más sencillo.

¿Qué aspecto tiene un vector de tipo factor? En el caso del data.frame experimento que acabamos de crear, al mostrar el campo experimento$Tratamiento se obtiene esta salida:

> experimento$Tratamiento
[1] T1 T1 T4 T1 T2 T5 T5 T2 T5 T2 T3 T1 T5 T3 T3 T4 T1 T4 T1 T5 T4 T3 T4 T4 T1 T4 T1 T2 T3 T5 T5 T4
[33] T2 T4 T4 T5 T1 T3 T3 T1 T3 T2 T2
Levels: T1 T2 T3 T4 T5

Como se ve, R muestra el contenido del vector, junto con sus niveles. Si sólo queremos mostrar los niveles del factor, podemos usar la función levels, cuya salida es esta:

> levels(experimento$Tratamiento)
[1] "T1" "T2" "T3" "T4" "T5"

¿Qué significan esas comillas? Para entenderlo, es preciso conocer la diferencia que hay entre el vector experimento$Tratamiento y un vector de tipo character “con el mismo contenido“. Podemos apreciar la diferencia pidiendo a R que convierta el factor en character mediante la función as.character:

 > as.character(experimento$Tratamiento)
[1] "T1" "T1" "T4" "T1" "T2" "T5" "T5" "T2" "T5" "T2" "T3" "T1" "T5" "T3" "T3" "T4" "T1" "T4" "T1"
[20] "T5" "T4" "T3" "T4" "T4" "T1" "T4" "T1" "T2" "T3" "T5" "T5" "T4" "T2" "T4" "T4" "T5" "T1" "T3"
[39] "T3" "T1" "T3" "T2" "T2"

Las comillas que aparecen en este caso indican precisamente que se trata de un vector de tipo character. Aunque los factores representan variables cualitativas, R los gestiona internamente mediante códigos numéricos. Al mostrarlos, a petición nuestra, el uso de las comillas permite distinguir entre character y factor.
Por cierto, cuando R carga un data.frame mediante read.table, los campos alfanuméricos (que contienen cadenas de texto, como el campo tratamientos de nuestro ejemplo) se convierten, por defecto en factores. Pero puede que, a veces, queramos preservar esos datos como variables de tipo character. Si es eso lo que queremos, basta con usar la opción stringsAsFactors=FALSE en la función read.table.

Los factores permiten modular el comportamiento de algunas funciones de R, de manera que la información que se obtiene tiene en cuenta esa organización en niveles de los valores que estamos analizando. Por ejemplo, en la quinta sesión vimos como usar la función summary para obtener un resumen descriptivo de un vector de datos (media y precentiles, básicamente). Ahora, que tenemos un data.frame con una clasificación por tratamientos, nos gustaría seguramente obtener esa misma información, pero para cada uno de los distintos tratamientos por separado. Lo primero que queremos hacer notar es que no sirve aplicar summary directamente al data.frame:

> summary(experimento)
Respuesta    Tratamiento
Min.   :2.600   T1:10
1st Qu.:4.540   T2: 7
Median :5.040   T3: 8
Mean   :4.791   T4:10
3rd Qu.:5.250   T5: 8
Max.   :6.060

Como se ve, se obtiene un resumen por separado para cada uno de los campos del data.frame: con Respuesta por un lado, y Tratamiento por otro, sin reconocer el vínculo entre ambas. Naturalmente, podemos ir separando cada uno de los grupos de tratamiento, y aplicarles la función summary a los vectores resultantes. Sería algo como esto (mostramos, por ejemplo el segundo grupo de tratamiento):ever

> (tratamiento1=experimento[experimento$Tratamiento=="T3",1])
[1] 5.41 4.87 5.68 5.34 5.71 5.24 5.06 4.52
> summary(tratamiento1)
Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
4.520   5.012   5.290   5.229   5.478   5.710

Pero como se puede apreciar, trabajar así es engorroso. De hecho, la razón para incluir aquí esto es para poder mencionar el símbolo == que aparece ahí. En una próxima entrada, examinaremos con mucho más detalle ese símbolo, y todo lo que tiene que ver con los valores y operadores lógicos, como primer paso hacia la programación en R. Pero, para obtener lo que queremos, hay una manera mucho más natural de proceder, usando la función tapply. Veamos como funciona:

> attach(experimento)
> tapply(Respuesta,Tratamiento, summary)
$T1
Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
4.700   4.955   5.080   5.079   5.228   5.370
$T2
Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
2.600   3.110   3.180   3.214   3.335   3.830
$T3
Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
4.520   5.012   5.290   5.229   5.478   5.710
$T4
Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
4.400   4.605   4.975   5.058   5.400   6.060
$T5
Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
4.560   4.935   5.015   5.039   5.182   5.390

El comando attach(experimento) le dice a R que vamos a usar el dataframe experimento, y nos permite usar las variables Respuesta y Tratamiento  directamente por sus nombres, en lugar de tener que andar escribiendo experimento$Respuesta y experimento$Tratamiento todo el rato.

El resultado, como se ve, es una tabla (de ahí la t inicial de tapply) que contiene, ordenados por filas, los resultados de aplicar la función summary a los valores de Respuesta, agrupadas según los distintos niveles del factor Tratamiento.

Otro ejemplo de las ventajas del lenguaje de factores se obtiene al pensar en la función boxplot. Si aplicamos esa función directamente al data.frame, ejecutando

>
boxplot(experimento)

se obtiene un gráfico que no es, desde luego, el más apropiado, porque de nuevo se tratan respuestas y tratamientos por separado:

BoxPlotDataFrameCompleto

Para obtener lo que queremos basta con escribir:

>
boxplot(Respuesta~Tratamiento)

y el resultado es un gráfico mucho adecuado para comparar los resultados de los tratamientos:BoxPlotDataFramePorFactorLa sintaxis Respuesta~Tratamiento que hemos usado en esta llamada es la que usa R para decir algo así como “la forma en que los valores de Respuesta dependen de los tratamientos“. En general, en R, el símbolo ~ se utiliza para indicar una dependencia o relación entre variables. En este caso, lo que estamos estudiando es si hay alguna relación entre la variable cualitativa (factor) Tratamiento, y la vaiable cuantitativa Respuesta. En particular, nos interesa saber si la respuesta media al tratamiento es distinta según el grupo de tratamiento que se considere. Se trata de un contraste de hipótesis, en la situación prototípica del método ANOVA unifactorial (¡porque sólo hay un factor, claro!) en el que, si llamamos

\mu_1,\mu_2,\ldots,\mu_k

a las respuestas medias de la población tratada con cada uno de los grupos de tratamiento, la hipótesis nula del contraste,sería

H_0=\{\mu_1=\mu_2=\cdots=\mu_k\}

Vamos a ver con más detalle algunas de las herramientas que R pone a nuestra disposición para llevar a cabo este tipo de contrastes ANOVA unifactoriales.

ANOVA unifactorial en R: hipótesis del modelo y tabla ANOVA

La teoría correspondiente se puede consultar en las notas de teoría, que están en esta entrada. Para empezar el trabajo de análisis con R, usaremos la función aov (de analysis of variance, claro) para calcular el modelo ANOVA y poder acceder a varias representaciones gráficas (por ejemplo, de los correspondientes residuos). La sintaxis sería algo como

aov(Respuesta~Tratamiento,datos)

siendo datos el dataframe que recibe como argumento, creado como hemos visto en el apartado anterior.  La ventana de resultados de Rcmdr muestra una respuesta bastante concisa, que no incluye lo que solemos llamar la tabla ANOVA, ni el p-valor del contraste de igualdad de medias:

> (datos.aov=aov(Respuesta~Tratamiento,datos))
Call:
aov(formula = Respuesta ~ Tratamiento, data = datos)
Terms:
Tratamiento Residuals
Sum of Squares    20.968566  5.693996
Deg. of Freedom           4        38

Residual standard error: 0.3870943
Estimated effects may be unbalanced

Pero si se almacenan el resultado de la función aov en una variable, como hemos hecho nosotros con datos.aov, entonces se puede llamar a la función plot así:

plot(datos.aov)

para obtener un análisis gráfico de los residuos del modelo, que proporcionan información adicional, útil para saber si se cumplen los requisitos de aplicación del ANOVA. En particular,  se obtiene una representación de los residuos frente a los valores predichos por el modelo (que son las medias), en la que debemos estar atentos a la aparición de posibles patrones en forma de cuña. En nuestro ejemplo, ese gráfico es ANOVA-ResiduosVsPredichosy no se aprecia en el ningún patrón que indique un fallo evidente en las hipótesis del modelo (lo que se aprecia, de hecho, es que la media de uno de los grupos -niveles- del factor aparece claramente separada de las cuatro restantes; el resto del análisis confirmará esto).

También se obtiene el llamado QQ-plot (de quantile vs quantile), en el que se comparan las distribuciones de los residuos concretos de nuestros datos, con los que se obtendría si la muestra  se ajustara exactamente a una distribución normal. Si las hipótesis del modelo se verifican, en este gráfico esperamos ver aparecer una nube de puntos situados muy aproximadamente a lo largo de una recta. En el ejemplo que nos ocupa, el gráfico QQ-plot es: ANOVA-QQplot-Residuosy como puede verse,  no es exactamente una recta, lo que sugiere que podemos tener algún problema con la hipótesis de normalidad en este ejemplo. ¿Qué hacemos, entonces? Bueno, por un lado, debemos tener presente que las muestras son pequeñas, así que, al igual que sucedía con los diagramas de cajas, el QQ-plot no es demasiado fiable como test de normalidad. En particular no sirve para garantizar que las poblaciones no son normales. Por otro lado, el contraste ANOVA es relativamente robusto frente a las desviaciones pequeñas de la normalidad, siempre y cuando las demás hipótesis del modelo se mantengan, y en particular cuando las varianzas de las poblaciones son homogéneas y todas las muestras son del mismo tamaño.  En este ejemplo, esa segunda condición no se cumple, lo cual nos obliga a ser especialmente cautos. Vamos a empezar por dedicar nuestra atención a la homogeneidad de las varianzas.

Un primer paso puede ser un diagrama de varianzas frente a medias. En este ejemplo se obtiene: ANOVA-VarianzasVsMediasLas varianzas de los grupos parecen ser bastante distintas, lo que podría preocuparnos, pero en este caso lo que más nos interesa es destacar que no hay relación aparente entre el tamaño de la medias y el tamaño de las varianzas. Esta información se une a la del diagrama de residuos frente a valores predichos. En este punto, y para obtener una representación alternativa de la dispersión, podemos volver a los boxplots paralelos por niveles del factor, que no contradicen a los otros gráficos. Como se ve, usamos toda la información disponible para ver si se detecta algún problema evidente, aunque en casos como este, con muestras tan pequeñas, es difícil ser concluyentes.

Antes de seguir adelante, una observación: puede pensarse en utilizar los diagramas de cajas para comprobar si la hipótesis de normalidad  de las poblaciones del modelo se verifican. Pero debe tenerse en cuenta que en el caso de muestras pequeñas, no hay tal cosa como un diagrama de cajas representativo de una distribución normal. El siguiente fragmento de código R ilustra lo que decimos. Se calcula una muestra aleatoria (la media y la desviación típica también se eligen de forma aleatoria) de tamaño 10 de una distribución normal, y se dibuja su diagrama de cajas. Si se ejecuta este código, se verá que la forma de ese diagrama de cajas varía tanto de unas muestras a otras, que realmente no tiene sentido tratar de usar la forma de ese diagrama como comprobación de normalidad.

(media=runif(1,min=-10,max=10))
#(desvTip=sample(seq(from=0, to=abs(media)/2,length.out=100),1))
(desvTip=runif(1,min=0,max=abs(media)/2))
k=10
(muestra=rnorm(k,mean=media,sd=desvTip))
boxplot(muestra)

Y si quieres verlo de otra manera, aquí tienes un applet Java (hecho con GeoGebra, del que ya hablaremos en otra entrada) que ilustra el uso de los boxplots para comprobar la hipótesis de normalidad de una población a partir de sus muestras. Un poco más abajo volveremos sobre el tema del chequeo de la hipótesis de normalidad. El diagrama de boxplots paralelos tampoco parece ser de mucha utilidad para comprobar la homogeneidad de las varianzas. Es fácil usar R para comprobar que aproximadamente el 20% de las muestras de tamaño 10 de la normal estándar tienen desviaciones típicas que se alejan de 1 (la desviación típica de la población) en más de 0.3 unidades.

Volviendo a las hipótesis del modelo ANOVA, ¿cuál es la conclusión? La conclusión, como ya hemos dicho, es que es difícil ser concluyentes con muestras tan pequeñas. ¿Cómo vamos a proceder, entonces? Vamos a usar la tabla ANOVA para calcular el p-valor del contraste. Si ese p-valor es concluyente (en uno u otro sentido), confiaremos en la robustez del método frente a pequeñas desviaciones de las hipótesis. Naturalmente, si el p-valor se sitúa en el umbral de la región de rechazo, deberíamos ser mucho más cuidadosos, mirar con escepticismo las conclusiones y plantearnos la posibilidad de obtener muestras de mayor tamaño.

Para obtener la tabla ANOVA tenemos dos opciones. La primera es usar summary junto con la función aov que ya hemos visto. Recuérdese que hicimos

(datos.aov=aov(Respuesta~Tratamiento,datos))

La tabla ANOVA obtenida mediante summary(datos.aov) es:

> summary(datos.aov)
Df Sum Sq Mean Sq F value   Pr(>F)
Tratamiento  4 20.969   5.242   34.98 2.91e-12 ***
Residuals   38  5.694   0.150
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

La otra opción es usar la sintaxis lm, de linear model, que en R se usa de forma genérica para representar la relación entre variables. De hecho, esta estructura lm juega un papel tan importante en R que, en breve, le dedicaremos una entrada. Para obtener la tabla ANOVA de esta manera hacemos:

anova(lm(Respuesta~Tratamiento,datos))

y el resultado es

> anova(lm(Respuesta~Tratamiento,datos))
Analysis of Variance Table
Response: Respuesta
Df Sum Sq Mean Sq F value    Pr(>F)
Tratamiento  4 20.969  5.2421  34.984 2.906e-12 ***
Residuals   38  5.694  0.1498
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Como se ve, de una u otra manera, resultados son esencialmente iguales en este ejemplo. Y, en particular, el p-valor  (2.906e-12) del  contraste es extremadamente pequeño. Así que podemos rechazar la hipótesis nula, concluir que las medias son distintas, y a la vez estar relativamente seguros de que esta conclusión no depende de pequeñas desviaciones de las hipótesis del  modelo.

En cualquier caso, el método ANOVA no agota, desde luego, el repertorio de recursos para realizar un contraste de igualdad de medios. Existen otros métodos, en especial métodos no paramétricos, que no dependen de las hipótesis de normalidad y homogeneidad de varianzas en la misma medida en que lo hace ANOVA. En R podemos usar uno de estos métodos para confirmar nuestras conclusiones, usando la función oneway.test. En concreto la forma de hacerlo es esta (se muestra la ventana de resultados):

> oneway.test(Respuesta~Tratamiento,datos)
One-way analysis of means (not assuming equal variances)
data:  Respuesta and Tratamiento
F = 36.7064, num df = 4.000, denom df = 17.373, p-value = 2.946e-08

El p-valor (2.946e-08) de este método alternativo confirma la conclusión que habíamos obtenido utilizando ANOVA.  ¿Por qué, entonces, no usamos siempre este método en lugar de ANOVA? La razón principal para hacerlo es que, cuando sus hipótesis se cumplen, el método ANOVA es, en general, más potente que este otro método.

Resultado significativo. ¿Y ahora qué?

Hemos rechazado la hipótesis nula, y por lo tanto, creemos que las medias de los distintos grupos son diferentes entre sí. El siguiente paso natural es, en muchos casos, tratar de comparar por parejas las medias. En nuestro caso eso significa comparar T1 con T2, T1 con T3, etcétera, hasta T4 con T5; un total de diez comparaciones, ya que

\binom{5}{2}=10.

Pero no siempre queremos hacer esto; a veces sólo nos interesa comparar una media concreta con las restantes, en lugar de comparar todas las parejas posibles.  En R se puede hacer esto usando summary combinado con lm de esta forma:

summary(lm(Respuesta~Tratamiento,datos))

y el resultado es:

> summary(lm(Respuesta~Tratamiento,datos))
Call:
lm(formula = Respuesta ~ Tratamiento, data = datos)
Residuals:
Min      1Q  Median      3Q     Max
-0.7087 -0.1590 -0.0180  0.1710  1.0020
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept)        5.07900    0.12241  41.492  < 2e-16 ***
Tratamiento[T.T2] -1.86471    0.19076  -9.775  6.4e-12 ***
Tratamiento[T.T3]  0.14975    0.18361   0.816    0.420
Tratamiento[T.T4] -0.02100    0.17311  -0.121    0.904
Tratamiento[T.T5] -0.04025    0.18361  -0.219    0.828
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

¿Cómo se interpreta esto? La parte más importante para nosotros es la tabla que aparece  bajo el epígrafe Coefficients. En esa tabla aparecen los contrastes de igualdad de medias de T1 (el primer nivel del factor; en el orden en que nosotros los hemos ordenado, claro) con el resto de los niveles. El valor que aparece al lado de (Intercept) es simplemente la media del primer grupo (el de T1),  mientras que en las restantes filas de la columna Estimate aparecen las diferencias entre las medias de T1 y las de los restantes niveles. Por ejemplo, la diferencia entre las medias de los tratamientos T1 y T4 es -0.02100. Y el p-valor para el contraste de diferencia de medias entre T1 y T4 está en la última columna, y vale 0.904. Los resultados de esta tabla indican que la media de T1 es significativamente distinta de T2 (tal vez sea buena idea volver a mirar los boxplots paralelos), mientras que no hay diferencias significativas entre la media de T1 y la de los tratamientos T3, T4 y T5.

Naturalmente, esta forma de trabajar privilegia a T1 frente a los demás tratamientos, y no responde a preguntas como ¿hay diferencia significativa entre T3 y T5? Para obtener la información de los diez contrastes por parejas posibles utilizamos la función pairwise.t.test de R, de esta manera (se muestra la ventana de resultados):

> pairwise.t.test(datos$Respuesta,datos$Tratamiento, p.adj="bonferroni")
Pairwise comparisons using t tests with pooled SD
data:  datos$Respuesta and datos$Tratamiento
T1      T2      T3 T4
T2 6.4e-11 -       -  -
T3 1       2.9e-11 -  -
T4 1       8.7e-11 1  -
T5 1       4.3e-10 1  1
P value adjustment method: bonferroni

En la tabla que aparece hay, en efecto, diez p-valores (se evitan las duplicidades mostrando sólo los valores de la parte de la tabla bajo la diagonal). Con esos p-valores podemos establecer qué parejas de medias son relamente distintas. En este ejemplo la conclusión (que ya avanzaban los boxplots) es que la media de T2 es significativamente distinta de todas las demás, pero que no hay diferencias significativas entre T1, T3, T4 y T5 (los p-valores correspondientes son – aproximadamente – 1). Sabemos (ver las notas de teoría), que es preciso ser cuidadosos a la hora de calcular esos p-valores, porque la simple repetición de diez contrastes aumenta mucho la probabilidad de cometer errores de tipo I. Así que R aplica alguno de los métodos habituales (hay varias estrategias, de nuevo remitimos a las notas de teoría y la bibliografía que allí aparece) para ajustar esos p-valores de manera que el error de tipo I se mantenga bajo control. Aquí hemos optado por pedirle a R que utilice la más clásica y conocida de esas estrategias, la de  Bonferroni.

Estos contrastes por parejas, incluso cuando se usan correcciones como la de Bonferroni, dependen de la hipótesis de normalidad de las poblaciones correspondientes a los factores. Para reducir la dependencia de nuestras conclusiones de esa hipótesis,  podemos utilizar algún método no paramétrico, de forma similar a lo que hicimos con la tabla ANOVA.  En este caso, el método más conocido para realizar el contraste es el de Tukey. Sin entrar en detalles, la forma de aplicarlo en R es esta (se muestra la ventana de resultados):

> summary(glht(datos.aov, linfct = mcp(Tratamiento = "Tukey")))
Simultaneous Tests for General Linear Hypotheses
Multiple Comparisons of Means: Tukey Contrasts
Fit: aov(formula = Respuesta ~ Tratamiento, data = datos)
Linear Hypotheses:
Estimate Std. Error t value Pr(>|t|)
T2 - T1 == 0 -1.86471    0.19076  -9.775   <1e-05 ***
T3 - T1 == 0  0.14975    0.18361   0.816    0.924
T4 - T1 == 0 -0.02100    0.17311  -0.121    1.000
T5 - T1 == 0 -0.04025    0.18361  -0.219    0.999
T3 - T2 == 0  2.01446    0.20034  10.055   <1e-05 ***
T4 - T2 == 0  1.84371    0.19076   9.665   <1e-05 ***
T5 - T2 == 0  1.82446    0.20034   9.107   <1e-05 ***
T4 - T3 == 0 -0.17075    0.18361  -0.930    0.883
T5 - T3 == 0 -0.19000    0.19355  -0.982    0.861
T5 - T4 == 0 -0.01925    0.18361  -0.105    1.000
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Adjusted p values reported -- single-step method)

Que es, en esencia, de nuevo una tabla con los diez p-valores, que viene a confirmar lo que habíamos obtenido anteriormente (observar los astericos a la derecha de la última columna, y su significado bajo la tabla; hay cuatro contrastes significativos, todos con T2 y uno de los otros grupos).

Representaciones gráficas. Bargraphs y otras opciones.

Después de un análisis como el que hemos hecho, los gráficos son una de las formas más comunes de presentar los resultados.  Una de las representaciones gráficas más populares, pese a sus inconvenientes son los gráficos de columnas adaptados para incluir líneas de error, como este para el ejemplo que nos ocupa: ANOVA-BargraphParaMedias

Pero hay que tener en cuenta que los gráficos de columnas se concibieron para representar frecuencias, no medias. Y en estos gráficos de columnas adaptados aparece un nivel de referencia (la base de las columnas, que en este caso R ha situado de forma automática en 3), nivel que a menudo no tiene ningún significado real en el contexto del problema que estamos estudiando.  Como puede verse, las barras, cuya altura se corresponde con la media de cada uno de los grupos, se acompañan de una línea o barra de error adicional. Un problema añadido de estos diagramas es que, en muchos casos, no está claro si esas líneas de error representan intervalos de confianza (y cuál es en ese caso el nivel de confianza), o errores estándar, o alguna otra medida de dispersión.  Hay que subrayar que debe dejarse siempre claro, en el texto que acompañe a la figura, cuál es la medida de dispersión que se ha utilizado. Una alternativa es eliminar las columnas, y representar siempre en paralelo intervalos de confianza para la media, como en este gráficoque corresponde a nuestro ejemplo (la idea procede del capítulo 7 de  Introductory Statistics with R, 2nd Ed, de P.Daalgard):ANOVA-GraficoIntervalosConfianzaSe muestra la nube de puntos correspondiente a cada muestra, y barras que corresponden con intervalos de confianza para cada una de esas medias. En este caso la única información necesaria para la correcta interpretación del gráfico es el nivel de confianza de esos intervalos.

Formato del fichero de datos y conclusión

Durante todo el trabajo anterior hemos usado un fichero de datos, llamado experimento-01.csv,.en el que los datos de la respuesta por grupos se presentan en dos columnas. La primera columna contiene el valor numérico de la respuesta, y la segunda el nivel (grupo) del factor Tratamiento.  En otras ocasiones esa misma información se puede presentar de otra manera, con cada uno de los grupos/niveles ocupando su propia columna, como en este fichero experimento-01a.csv, que contiene los mismos datos del ejemplo que hemos usado. El capítulo correspondiente de las notas de teoría contiene un fichero de código (y esa será la versión que en el futuro trataré de mantener actualizada) para el análisis ANOVA unifactorial, que permite realizar todos los pasos del análisis que hemos descrito en esta entrada. Este fichero trata de reconocer automáticamente el formato de entrada del fichero de datos (en las dos situaciones más comunes, las que hemos descrito),  para realizar el análisis correcto. Es conveniente comprobar, en cualquier caso, que los datos se han reconocido correctamente y que no aparece ningún mensaje de error en la lectura. Ese reconocimiento automático utiliza estructuras de programación de R. La existencia de esas estructuras de programación es una de las razones que hacen de R una herramienta tan versatil y potente, así que volveremos sobre ellas en próximas entradas.

Gracias por la atención.

Apuntes de Estadística

Esta entrada sirve simplemente para alojar la versión (actualizada y) completa de mis apuntes de teoría de las clases de Estadística. Se pueden descargar desde este enlace:

000-CursoEstadistica

El documento es un fichero pdf (de unos 6.5 megas), con índice navegable, y que contiene como adjuntos los ficheros de código (en R , Calc, ejemplos GeoGebra en html, etc). Los capítulos que aparecen marcados con asteriscos en el índice todavía no están disponibles, se irán añadiendo en el futuro (¿lejano…?)

ULTIMA ACTUALIZACIÓN: 15 de enero de 2013.