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.

Anuncios

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.

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.

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.

14ª sesión con Rcmdr: Contraste de hipótesis y muestras precocinadas.

Muestras precocinadas, al gusto del cliente

No voy a usar esta entrada para cantar las alabanzas de R como plataforma de análisis de datos. Eso, a estas alturas, es un clamor que ha llegado hasta lugares tan inesperados como las páginas del New York Times. Pero sí quiero hablar del uso de R en la enseñanza de la Estadística, porque eso ayudará a entender la primera parte de esta entrada:. R no se diseñó pensando en la enseñanza, y eso se le nota (“pero mucho”, me parece oír decir a algunos de mis alumnos…).

Por ejemplo, cuando los estudiantes aprenden el concepto de contraste de hipótesis, los ejercicios típicos de los libros de texto se parecen a este:

Se ha recibido un envío de latas de conserva, de las que se afirma que el peso medio son 1000gramos. Al examinar una muestra aleatoria de 5 latas se obtuvo un peso medio muestral de 995 gramos, con una cuasivarianza muestral s²= 19.6. Al nivel de confianza 95%, ¿se puede aceptar que el peso medio son 1000 gramos? Obtener también el p-valor de este contraste.

Naturalmente, al tratarse de un contraste sobre la media, con una muestra de tamaño pequeño, tendremos que usar la distribución t de Student. Y, como no podía ser de otra manera, R incluye una función, llamada t.test, diseñada para realizar este tipo de contrastes. Pero si se intenta utilizar la función t.test para este ejercicio, descubrimos que no es demasiado fácil. La ayuda de R para esta función dice que la función se usa escribiendo

t.test(x, ...)

donde x es un vector de datos numéricos, que contiene los datos de la muestra, claro. Pero no hay nada previsto para introducir directamente los estadísticos muestrales (la media \bar X y la cuasivarianza s^2, y el tamaño n de la muestra), y pedirle a R que haga el contraste usando esos valores sin disponer de la muestra.

Hay que suponer que los programadores de R piensan seguramente que el punto de partida siempre son los datos. Y que nadie hace (en el mundo real, fuera de las aulas y los libros de texto) cálculos estadísticos sin partir de los datos.  Es una suposición muy razonable para el mundo real, pero choca con la forma habitual de enseñar la Estadística.

Otros programas ofrecen la posibilidad de hacer directamente este tipo de ejercicios. ¿Es eso bueno? Podríamos enzarzarnos en una discusión pedagógica, pero lo cierto es que no me interesa demasiado. Una de las principales razones (aparte, naturalmente, del hecho de que es gratuito) por las que usamos R en las clases es porque creemos que es una inversión rentable para nuestros alumnos, porque van a poder seguir usándolo en el futuro, en aplicaciones profesionales. Y, a su vez, una de las razones que explican su éxito en entornos profesionales es su enorme flexibilidad. Si hay algo que echas en falta en R, siempre puedes programarlo tu mismo.  Y eso es lo que ha hecho alguien por nosotros en este caso. Para poder hacer los problemas de los libros de texto directamente, vamos a utilizar una función de R que prepara muestras “al gusto del consumidor”. Es decir, nosotros le decimos a R cuáles son  los valores de

\bar X, s^2 \mbox{ y } n

y R nos “cocinará” una muestra

x_1,x_2,\ldots,x_n

con exactamente esos parámetros, para que podamos dárselos como alimento a la función t.test. La función que hace esto por nosotros es mvrnorm. Por ejemplo, estos comandos de R fabrican una muestra de una población aproximadamente normal, cuya longitud, media y cuasivarianza se han elegido de antemano. Y a continuación comprobamos que esos valores son los que queríamos.

# Fijamos los parámetros deseados de la muestra,
n=10
mediaMuestral=8.2
cuasiVarianzaMuestral=1.3
# formamos la muestra,
(muestra=mvrnorm(n, mediaMuestral,cuasiVarianzaMuestral, empirical = TRUE))
# y comprobamos que responde a nuestros deseos.
mean(muestra)
var(muestra)

La parte interesante del resultado es la comprobación, claro:

> mean(muestra)
[1] 8.2
> var(muestra)
[,1]
[1,]  1.3

Con esto estamos preparados para ver como funciona la función t.test de R, y el resto de funciones que realizan contrastes de hipótesis.

Contrastes de hipótesis en R (t.test y sus parientes)

Vamos a ver como usar la función t.test de R para realizar un contraste de hipótesis sobre la media. Usaremos este ejemplo (que apreciará cualquiera que haya estado parado en un paso a nivel canadiense, esperando a que termine de pasar el mercancías de turno):

Una compañía ferroviaria canadiense afirma que sus trenes de mercancías no bloquean los pasos a nivel durante más de 8 minutos, en promedio. Una muestra aleatoria de 10 tiempos de bloqueo dio como resultado estos valores (en minutos):

10.1, 9.5, 6.5, 8.0, 8.8, 12,  7.2, 10.5, 8.2, 9.3

Empezamos por observar que en este caso tenemos todos los valores de la muestra. Llamando μ al tiempo medio de bloqueo, queremos usar estos valores para contrastar la hipótesis nula:

H_0=\{ \mu \leq \mu_0=8 \}

Y naturalmente la hipótesis alternativa es:

H_a=\{ \mu > \mu_0=8 \}

El estadístico del contraste es

\dfrac{\bar X-\mu_0}{\dfrac{s}{\sqrt{n}}}

Hagamos primero los cálculos del contraste “a mano”, sin recurrir a t.test. Vamos a fijar un nivel de significación del 5%, es decir, 0.05. Puesto que se trata de una muestra pequeña (n=10), usaremos la distribución t de Student para el cálculo del p-valor.

datos=c(10.1, 9.5, 6.5, 8.0, 8.8, 12,  7.2, 10.5, 8.2, 9.3)
mu0=8
(xBarra=mean(datos))
(s=sd(datos))
(n=length(datos))
EstadisticoContraste = ( xBarra - mu0 ) / ( s / sqrt(n) )
(pValor=1-pt(EstadisticoContraste,df=n-1))

El resultado es este:

> datos=c(10.1, 9.5, 6.5, 8.0, 8.8, 12,  7.2, 10.5, 8.2, 9.3)
> mu0=8
> (xBarra=mean(datos))
[1] 9.01
> (s=sd(datos))
[1] 1.631938
> (n=length(datos))
[1] 10
> (EstadisticoContraste = ( xBarra - mu0 ) / ( s / sqrt(n) ))
[1] 1.957121
> pValor=1-pt(EstadisticoContraste,df=n-1)
> (pValor=1-pt(EstadisticoContraste,df=n-1))
[1] 0.04101152

De manera que, con un nivel de significación 0.05 (mayor que el p-valor), tenemos evidencia empírica para rechazar la hipótesis nula y concluir que los trenes bloquean el paso a nivel más tiempo del que dice la empresa. Veamos ahora como hacer el contraste usando t.test:

datos=c(10.1, 9.5, 6.5, 8.0, 8.8, 12,  7.2, 10.5, 8.2, 9.3)
mu0=8
t.test(datos,mu=mu0, alternative = c("greater"), conf.level = 0.95)

Como se ve, aparte del vector de datos, le hemos indicado a R el valor de \mu_0, y mediante las opciones  alternative = c(“greater”)conf.level = 0.95 hemos seleccionado un contraste de cola derecha (greater, o mayor que) y el nivel de significación deseado (conf.level, R usa la misma terminología que para los intervalos de confianza) .

El resultado es este:

> t.test(datos,mu=mu0, alternative = c("greater"), conf.level = 0.95)
One Sample t-test
data:  datos
t = 1.9571, df = 9, p-value = 0.04101
alternative hypothesis: true mean is greater than 8
95 percent confidence interval:
8.063996      Inf
sample estimates:
mean of x
9.01

Y como se ve, la respuesta contiene tanto el valor del estadístico de contraste ( t = 1.9571), como el p-valor = 0.04101. Además, para que la interpretación del resultado sea más fácil, y para que podamos comprobar que estamos haciendo lo que deseamos, R describe la hipótesis alternativa del contraste. Como subproducto se obtiene lo que R llama un intervalo de confianza, pero que en el caso de un contraste unilateral como este, no tiene la forma que seguramente espera un estudiante de un curso de introducción a la Estadística. Para obtener la forma del intervalo de confianza que aparece en los libros introductorios, basta con ejecutar t.test para hacer un contraste bilateral (que es la opción por defecto):

> t.test(datos,mu=mu0, alternative = c("two.sided"), conf.level = 0.95)
One Sample t-test
data:  datos
t = 1.9571, df = 9, p-value = 0.08202
alternative hypothesis: true mean is not equal to 8
95 percent confidence interval:
7.842582 10.177418
sample estimates:
mean of x
9.01

El intervalo de confianza que obtenemos seguramente se parece más a lo que se esperaba.

Si has estudiado los  contrastes de hipótesis en el orden habitual, seguramente te estarás preguntado ¿y el contraste para la media con la Z, la normal estándar?  Lo cierto es que esos contrastes Z son casi exclusivamente “ejemplos de libro de texto”, que no se usan en las aplicaciones reales. Y no se incluyen en R por defecto (¿hemos dicho ya que R no se diseñó pensando en la enseñanza?). Pero eso no significa que no estén disponibles. Basta con cargar una librería, cuyo revelador nombre es TeachingDemos, y listos: ya tenemos disponible una función z-test, con la que podemos hacer cosas como esta: suponemos que los datos del ejemplo anterior proceden de una población normal con desviación típica σ conocida, e igual a 1.6, y realizamos el mismo contraste para la media que hicimos antes:

datos=c(10.1, 9.5, 6.5, 8.0, 8.8, 12,  7.2, 10.5, 8.2, 9.3)
z.test(datos, mu =mu0, sd=1.6, alternative = c("greater"), conf.level = 0.95)

El resultado es:

> datos=c(10.1, 9.5, 6.5, 8.0, 8.8, 12,  7.2, 10.5, 8.2, 9.3)
> library("TeachingDemos")
> z.test(datos, mu =mu0, sd=1.6, alternative = c("greater"), conf.level = 0.95)
One Sample z-test
data:  datos
z = 1.9962, n = 10.000, Std. Dev. = 1.600, Std. Dev. of the sample mean = 0.506, p-value = 0.02296
alternative hypothesis: true mean is greater than 8
95 percent confidence interval:
8.177763      Inf
sample estimates:
mean of datos
9.01

Es un ejercicio saludable repetir estas cuentas “a mano”, usando pnorm, para comprobar que somos capaces de obtener el p-valor de este contraste de esa forma, sin recurrir a z.test.

Volviendo a t.test, esa misma función puede usarse para contrastes de medias entre dos poblaciones, simplemente introduciendo sendos vectores con las muestras correspondientes. Por ejemplo, dado este ejercicio:

Una empresa quiere introducir un nuevo método de trabajo, y quiere saber si ese método influye sobre el tiempo que se emplea en llevar a cabo una determinada tarea. Se han seleccionado aleatoriamente dos grupos, cada uno de cinco empleados, a los que se les ha pedido que lleven a cabo la tarea por uno de los dos métodos. La tabla de tiempos empleados es esta:

\begin{tabular}{|l|c|c|c|c|c|} \hline\mbox{M\'etodo antiguo}&20&18&18&19&17\\ \hline \mbox{M\'etodo nuevo}&18&15&17&16&17\\ \hline \end{tabular}

¿Se puede afirmar, con un nivel de significación del 5%, que el  tiempo que se emplea en esa tarea es el mismo para ambos métodos?

El resultado de llevar a cabo este contraste en R, utilizando la función t.test, es el siguiente:

> muestra1=c(20,18,18,19,17)
> muestra2=c(18,15,17,16,17)
> t.test(muestra1,muestra2,alternative=c("two.sided"),conf.level=0.95)
  Welch Two Sample t-test
data:  muestra1 and muestra2
t = 2.4962, df = 8, p-value = 0.03716
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
0.1371168 3.4628832
sample estimates:
mean of x mean of y
18.4      16.6

Algunos comentarios sobre este ejemplo: en este caso el contraste es bilateral, y por eso la opción que usamos es alternative=c(“two.sided”). Si fuera unilateral, usaríamos, como en el ejemplo previo, o bien alternative=c(“greater”), o bien alternative=c(“less”), y en ambos casos las muestras se usan en el orden en que se escriben. Es decir, que si usáramos muestra1,muestra2,alternative=c(“greater”), estaríamos utilizando como hipótesis alternativa que la media de la población 1 es mayor (greater) que la media de la población 2.
Además, el lector atento habrá observado que no le hemos indicado a R si estamos suponiendo que las varianzas de las poblaciones se suponen iguales o distintas (son muestras de tamaño muy pequeño, así que es necesario hacer esa distinción). La opción var.equal=FALSE (o bien, TRUE, claro) nos permite elegir explícitamente el comportamiento de R. Por defecto R utiliza, lógicamente, el valor FALSE. En particular, se usa una corrección de tipo Welch para los grados de libertad. Y digo una porque esta mañana he estado echando un vistazo en mi biblioteca, y he encontrado al menos tres expresiones distintas (todas etiquetadas como “de Welch”) para el cálculo de los grados de libertad en esta situación. R usa una de ellas, y eso significa ¡ATENCIÓN! que el intervalo de confianza que calcula R puede no coincidir exactamente con el que se obtiene usando las fórmulas de tu libro o de tus apuntes de clase. SI quieres asegurarte esa coincidencia, calcula el intervalo a mano (usando qt, y la fórmula de Welch que corresponda).

Los contrastes de proporciones también pueden hacerse en R, utilizando en este caso la función prop.test. Por ejemplo, si tenemos dos muestras independientes, de tamaños n1=100 y n2=250, en las que el número de éxitos es  respectivamente, e1=33 y e2=80, la forma de usar prop.test para calcular un contraste de igualdad de proporciones se ilustra en este ejemplo:

> prop.test(c(33,80),c(100,250),alternative = c("two.sided"),conf.level = 0.95, correct = FALSE)
2-sample test for equality of proportions without continuity
correction
data:  c(33, 80) out of c(100, 250)
X-squared = 0.0327, df = 1, p-value = 0.8566
alternative hypothesis: two.sided
95 percent confidence interval:
-0.09879831  0.11879831
sample estimates:
prop 1 prop 2
0.33   0.32

Quizá lo mas peculiar es la forma en la que se introducen los datos muestrales, de manera que el primer vector contiene los tamaños de ambas muestras, y el segundo los recuentos de éxitos. El resto de las opciones son fáciles de comprender, con la excepción de correct = FALSE.  Con esta opción desactivamos una corrección de continuidad que aplica R, para que los resultados sean (menos precisos pero) más parecidos a los que se aprende a calcular en un curso de introducción a la Estadística, usando la distribución normal estándar Z.

Obsérvese que en este caso de proporciones, la sintaxis de R hace innecesario el “cocinado” de muestras que hemos explicado en la primera parte de la entrada.

También podemos llevar a cabo contrastes sobre la varianza (y, por consiguiente, la desviación típica). Para una sola población, disponemos , en la librería TeachingDemos (de la que ya hemos hablado), de la función sigma.test. En este ejemplo vemos como usarla para un contraste unilateral sobre la desviación típica de una población (se muestra la salida en la ventana de resultados):

> muestra=c(9.6,11.4,8.5,10.8,10.5,10.7,10.1,9.8,8.9,11.6,11.4,12.9)
> library("TeachingDemos")
> sigma.test(muestra,sigma=2,alternative=c("less"),conf.level=0.95)
One sample Chi-squared test for variance
data:  muestra
X-squared = 4.1842, df = 11, p-value = 0.03578
alternative hypothesis: true variance is less than 4
95 percent confidence interval:
0.000000 3.658437
sample estimates:
var of muestra
1.521515

Para los aprendices de la Estadística es bueno leer la ayuda de la función sigma.test para entender porque está función aparece “relegada” a la librería TeachingDemos.

Si tenemos que hacer un contraste de igualdad de varianzas entre dos poblaciones independientes, podemos hacerlo empleando la función var.test (no es necesario cargar librerías adicionales para usarla. ni siquiera TeachingDemos). Para no extendernos más, remitimos al lector  a la ayuda de esa función, y ponemos aquí punto final a esta entrada.

Gracias por la atención.

13ª sesión con Rcmdr: un momento de respiro, y algo de ayuda.

Esta semana, por diversas razones,  ha sido complicada, y además la mayoría de mis alumnos se examinan la semana que viene. Así que seguramente todos, ellos y yo, agradeceremos una entrada del blog más relajada que las anteriores.  Además, estoy escribiendo una entrada más larga, con los primeros pasos en R como lenguaje de programación, y todavía no está acabada. Así que, mientras ese trabajo llega a puerto, haremos una excursión breve a algunos temas sueltos que han ido quedando por el camino. El denominador común es mejorar nuestra experiencia como usuarios de R.

Lo primero es hacer justicia con quienes nos abrieron el camino

Aunque este blog aspira a proporcionar una lectura autosuficiente sobre R, lo cierto es que hay muchas fuentes de información que yo usé, en su momento,  para aprender a usar R y a las que sigo acudiendo en busca de lo mucho, muchísimo, que siempre queda por delante. Y os animo a que las exploréis. Hay gente haciendo cosas increíbles, que nos hacen soñar con llegar a ser algún día la décima parte de buenos en lo que hacemos.

La primera de esas fuentes no hay que ir a buscarla muy lejos. En cualquier instalación de R hay disponible un documento, al que todos nos referimos como R-Intro.  En Windows, por ejemplo, para localizarlo, vamos a la ventana RGui (donde aparecen los gráficos, no la de Rcmdr), y en el menú Ayuda seleccionamos Manuales (en pdf) y luego An Introduction to R, como en esta figura.

Este documento, en inglés (en este enlace hay una traducción al castellano de una versión algo más antigua),  es relativamente fácil de leer y constituye un curso de choque a algunas de las posibilidades del programa. Es el primer paso para empezar a convertirse en un usuario avanzado de R, aunque -pese a su título- a veces da la sensación de que está más pensado como una guía de referencia que como una ayuda al aprendizaje.

Aparte de este, y otros documentos parecidos, hay varios (muchos) blogs disponibles con un espíritu similar a este.  Sin embargo, no conozco ninguno, en español, que se plantee como objetivo enseñar a usar R a un completo novato, como acompañamiento para un primer curso en Estadística. Además, incluso en los blogs en inglés, muchas veces se presentan muchos detalles técnicos de R antes de empezar con las aplicaciones estadísticas. Y por esas razones, entre otras, existe este blog.

Volviendo a los blogs en inglés, hay dos que he consultado bastantes veces, y que tienen nombres casi idénticos :

  • R tutorial. Contiene a su vez varios temas, de los cuales los más relacionados con este blog son por un lado Elementary Satistics with R y por otro R Introduction. Contiene numerosos ejemplos de problemas estadísticos resueltos con R.
  • R tutorials. Probablemente el blog más parecido a este en sus intenciones. Pero no usa Rcmdr, sino la interfaz más espartana de R (y yo no me atrevo a infligir semejante castigo a mis ya suficientemente sufridos lectores). En cualquier caso, contiene un montón de información útil y bastante bien organizada.

Pero el universo de los blogs sobre R es mucho más amplio. De hecho, existe una página (R bloggers) que se dedica a recopilar información sobre blogs que hablan de R (y tienen más de 400 registrados, mayoritariamente en inglés, pero también en otros idiomas…). Así que puedes buscar el que más te guste.

Ayuda directa sobre comandos desde el propio R

Por otra parte, el propio R, como cualquier programa mínimamente decente, incluye documentación interna sobre sus propios comandos. Por ejemplo, ya conocemos el comando sample,  que hemos usado en la novena sesión para fabricar muestras. Para acceder a la ayuda interna sobre el comando, tecleamos en R su nombre precedido de un paréntesis, así:

?sample

Al ejecutar este comando se abrirá nuestro navegador de internet (pero no hace falta estar conectados, la ayuda reside en nuestro ordenador), y nos mostrará  la documentación del comando sample (en inglés):De forma equivalente, puedes usar el comando help(sample) para obtener el mismo documento.

Toda la documentación de R está disponible también en internet, con copias en múltiples servidores. Así que muchas veces, a pesar de todos los sistemas de ayuda que hemos mencionado,  lo más rápido es usar un buscador de internet, teclear algo como “Binomial distribution with R“, y hojear los resultados.

Asignaciones con flechas

Es algo que he querido mencionar en el  blog desde hace tiempo, pero que nunca parecía encontrar su sitio. Hasta hoy, claro. SI vais a consultar el código r escrito por otras personas, descubriréis que, donde yo pondría:

datos = c(2,5,12,3,4)

mucha otra gente escribe:

datos <- c(2,5,12,3,4)

Lo que estamos haciendo, en cualquier caso, es una asignación, en este caso guardando un vector de números en la variable datos. Y la notación tradicional en R para las asignaciones utiliza esa especie de flecha <-  en lugar del símbolo igual. Es una buena idea, en principio, porque el símbolo = ya se utiliza en matemáticas para indicar que dos cosas son iguales, y no con este sentido de “guarda esto en la variable…” Pero en la práctica hay tres razones que me llevan a usar = en lugar de la flecha.

  • la primera, es que el símbolo = se utiliza en muchos lenguajes de programación para las asignaciones (no en todos, pero sí en muchos de los más usados).
  • la segunda es que el símbolo elegido, la flecha <- , no me parece demasiado afortunado. En otros lenguajes se usan alternativas como := (dos puntos seguidos de un igual), que me parecen más logradas.
  • y finalmente, ¡da igual! Puedes cambiar las flechas por iguales o viceversa, y R no protestará (al menos, ninguna versión moderna de R).

Gracias por la atención.