Programación literaria en R con knitr y markdown.

Knitr y Rmarkdown.

Publicado el 2016-05-05

ültima actualización: 2016-05-06

Esta entrada es la primera de una serie que quiero dedicar a describir algunas herramientas del ecosistema R que en los últimos tiempos he ido incorporando a mi forma de trabajar. Quiero escribir sobre los siguientes temas, con la intención confesa de estimular la curiosidad de quien pueda leer estas entradas:

• Knitr y Markdown.
• Programación literaria.
• Shiny.
• Reproducible Research.
• El paquete exams de R.

Tiempo habrá, en otras entradas de esa serie, para adentrarse en los detalles y la historia. Lo que quiero hacer hoy es simplemente empezar. Y no hay mejor manera de empezar que con un ejemplo muy sencillo.

Las cosas serán más fáciles si ya eres usuario de RStudio. Si no lo eres: (a) te lo recomendamos una vez más y (b) puedes encontrar algo más de información e instrucciones de instalación en nuestro Tutorial00 (pdf). Primero explicaremos como proceder en RStudio, y luego daremos algunas indicaciones para quienes no lo usen.

Usando RStudio para crear un documento Rmarkdown.

En RStudio empieza por ir al menú Tools y dentro de ese menú selecciona Global Options. Se abrirá una ventana de opciones como esta:

En la que, como hemos indicado, debes seleccionar Sweave (más adelante, en otra entrada, hablaremos más de Sweave). Al hacerlo aparecerán las opciones de la siguiente figura:

Y como muestra la flecha debes asgurarte de que has seleccionado Knitr bajo la opción Weave Rnw files using. Puedes dejar el resto de las opciones como están por el momento. Haz click en Ok.

Ya estamos casi listos. Ahora, usando el menú File de RStudio, selecciona la entrada New File y, dentro de esta, elige R Markdown. Aparecerá esta ventana:

Elige un título para el documento, añade la información de autoría y pulsa en Ok. Yo he elegido como título Comenzando con Rmarkdown. ¡Pero cuidado! El título que has elegido es sólo eso: un título (en particular puede contener espacios, acentos, etc). Tu documento todavía no está guardado en disco. Enseguida lo guardaremos. Pero antes fíjate en que en la ventana del Editor de RStudio ha aparecido un documento con un texto inicial. Más abajo vamos a analizar la estructura de ese texto.

Pero si te fijas verás que el documento sigue siendo Untitled, inidcando que aún no lo hemos guardado. Hazlo ahora, usando el menú File/Save as y guarda el documento en una carpeta adecuada. Tendrás que elegir el nombre del fichero, que no tiene porque tener nada en común con el título que has elegido antes. Eso sí, como de costumbre, te recomiendo que no uses espacios, acentos o caracteres especiales en el nombre del fichero. No te preocupes por la extensión: RStudio añadirá de forma automática la extensión Rmd al fichero. Yo he elegido como nombre testRmarkdown01.

¿Ya lo has guardado? Estupendo. Pues ahora, antes de hacer nada más, vamos a ver hacia dónde nos encaminamos: pulsa el botón Knit HTML que aparece sobre el editor de RStudio, y que puedes ver en la siguiente figura:

Tras pulsarlo verás que, en el espacio que normalmente ocupa la consola de RStudio aparece fugazmente una nueva pestaña titulada R markdown que muestra mensajes con información sobre el proceso que RStudio está llevando a cabo con nuestro documento. Y digo fugazmente porque, con un documento tan breve, en apenas unos instantes verás aparecer una nueva ventana, el visor de RStudio, que muestra el resultado que hemos obtenido al procesar ese documento.

Por si tienes algún problema, este enlace muestra el resultado que deberías obtener.

¿Y si no uso RStudio?

En ese caso puedes descargarte el fichero de código para Knitr desde este enlace:
testRmarkdown01.Rmd
Guárdalo en tu directorio de trabajo con R y ábrelo con el editor de texto que uses habitualmente con ficheros de código R.

Para procesarlo (el equivalente del botón Knit de RStudio), ejecuta este código en la consola de R:

rmarkdown::render("testRmarkdown01.Rmd")

Al hacerlo, en tu directorio de trababjo aparecerá un fichero llamado testRmarkdown01.html. Ábrelo con tu navegador de internet para ver el resultado que ha producido Knitr. Debería ser parecido a este.

En lo que resta de esta entrada, cada vez que yo diga que hay que pulsar el botón Knit de RStudio, tú deberás ejecutar de nuevo esa instrucción para volver a generar el documento HTML.

El documento HTML producido por Knitr.

El documento que estás viendo es un documento HTML, como los que se usan en las páginas Web. Es un documento con un formato visual simple, pero que contiene elementos visuales como texto en cursiva, negrita, títulos con dierente tamaño de letra. Además, se incluyen enlaces reconocibles, como en cualquier página web. Pero, como puedes ver, lo más interesante es que en ese documento aparecen recuadros sombreados con código R y a continuación el resultado de la ejecución de ese código. El primero de esos bloques usa la función summary con el conjunto de datos cars y el resultado aparece justo debajo. Además aparece un gráfico, también realizado con R, aunque en este caso el código que se ha usado para fabricar el gráfico es invisible. De hecho, la última frase del documento hace referencia a esto y explica de forma sucinta cómo se ha conseguido.

Hay muchas cosas ocurriendo a la vez en esta ventana, ya lo sé. Pero tendremos ocasión de ir desgranando poco a poco los detalles técnicos. lo más importante en este momento del proceso es captar la idea general: tenemos un documento inicial (el fichero testRmarkdown01.Rmd) y RStudio lo ha procesado (usando Knitr) para dar como resultado esta página web que contiene código R y los resultados que produce ese mismo código. Como pronto verás, esas páginas web son un documento extremadamente flexible que nos va a permitir crear documentos muy sofisticados de forma muy productiva. No sólo eso, sino que al tratarse de documentos de texto muy sencillos, podemos escribir programas en R que se encarguen de crear y publicar de manera automática esas páginas web cuando, por ejemplo, los datos que nos interesan hayan cambiado.

Estructura del documento en Rmarkdown.

El documento testRmarkdown01.Rmd es nuestro primer ejemplo de documento Rmarkdown. Vamos a analizar la forma en la que está construido. A grandes rasgos el documento contiene tres tipos de elementos:

La cabecera.

Que aparece justo al principio del documento, comprendida entre las dos líneas que contienen ---. Ese bloque de cabecera es, en este documento de muestra:

title: "Comenzando con RMarkdown."
author: "www.postdata-statistics.com"
date: "4 de mayo de 2016"
output: "html_document"

Este bloque está escrito en su propio lenguaje, denominado YAML. Pero no te preocupes demasiado. Por el momento todo lo que necesitamos saber es esto:

• La función de esas líneas es evidente. Aquí se determina el título, autor, fecha y, en la última línea, cuál es el tipo de documento que queremos producir. En este caso HTML para producir un documento para la web. Pero más adelante veremos cómo producir, por ejemplo, un pdf o un documento docx para Microsoft Word.
• Además, este bloque es opcional. Si borras todo el bloque entre las dos líneas con --- (incluidas esas líneas), Knitr seguirá siendo capaz de procesar el documento resultante.
• Si sabes algo sobre creación de documentos HTML es posible que el ejemplo anterior te haya parecido muy básico. Es cierto. Pero cuando aprendamos un poco más, podremos usar este bloque YAML de la cabecera para incluir, por ejemplo, nuestras propias hojas de estilo CSS. Se pueden llegar a hacer cosas bastante sofisticadas.

Texto en markdown.

Además, el documento contiene párrafos de texto escritos en markdown. El lenguaje markdown sirve para representar el formato de un documento de una forma muy sencilla. A menudo se dice que la sintaxis básica de markdown se aprende en un minuto. Por ejemplo, para escribir un texto en negrita basta con rodearlo con dos asteriscos a cada lado. En el documento testRmarkdown01.Rmd tienes un ejemplo. El segundo párrafo contiene esta frase:

When you click the **Knit** button a document will be generated ... 

Como ves, la palabra Knit aparece con dos asteriscos a cada lado. Al procesar el documento con Knitr, el resultado es que la palabra Knitr aparece en negrita, como puedes ver en esta figura:

También puedes ver, en el primer párrafo de texto, lo fácil que resulta insertar un enlace a una página web. Basta con escribir

  <http://rmarkdown.rstudio.com>

para crear ese enlace. Si quisieras hacer lo mismo usando HTML tendrías que escribir:

  <a href="http://rmarkdown.rstudio.com">http://rmarkdown.rstudio.com</a>

Como ves, bastante más complicado. Esa simplificación es la clave que convierte a markdown en una herramienta extremadamente productiva. Porque puede sonar a exageración, pero te aseguro que el tiempo que dediques a aprender a usar markdown será una de las inversiones más rentables que puedes hacer en términos de productividad. Un buen ejemplo puede ser esta página:

http://www.markitdown.net/markdown

en la que puedes ver a la izquierda ejemplos de código markdown y a la derecha su traducción cuando se convierten en HTML para producir una página web. Además ese documento tiene la ventaja de que proporciona una demostración interactiva del funcionamiento de markdown. Haz un cambio a la izquierda y verás su efecto a la derecha. Es muy conveniente que dediques algún tiempo a hacer experimentos allí para familiarizarte con los ingredientes básicos del markdown. Tal vez no sea un minuto, pero tampoco vas a necesitar una hora, desde luego. A cambio, después podrás usar markdown en la creciente colección de sitios web y programas compatibles con su uso: WordPress, GitHub, los foros de StackOverflow y sus asociados, etc. Usar markdown no es la panacea universal, pero para muchos de nosotros se ha convertido en una de esas herramientas que te llevan a preguntarte cómo podías pasarte sin ellas antes.

Bloques de código R.

Pero hemos dejado para el final lo mejor. El documento contiene además bloques de codigo R (en inglés se habla de R code chunks). Por ejemplo, el primero de esos bloques es muy breve, y sólo contiene una llamada a una función de R:

summary(cars)

La siguiente figura muestra la manera de indicarle a Knitr los límites de uno de esos bloques de código:

---

---

Cuando Knitr procesa el fichero testRmarkdown01.Rmd va buscando estos bloques de código y los ejecuta secuencialmente, uno tras otro. Y a medida que los ejecuta va insertando las respuestas de R en el documento html que produce, dentro de unos bloques de resultados que aparecen recuadrados (ahora con fondo blanco) bajo el correspondiente bloque de código.

---

---

Prueba a ejecutar summary(cars) en tu consola de R y verás que obtienes el mismo resultado web que aparece en el documento elaborado por Knitr.

Un inciso sobre Programación Literaria.

Este ejemplo es tan sencillo que puede que no resulte obvio lo que acabamos de explicar. Por eso voy a insistir en el punto clave. Tradicionalmente, ejecutaríamos ese código en R y entonces copiaríamos y pegaríamos el resultado en algún documento, informe, página web, etc. Esa es la forma tradicional de incorporar los resultados de la ejecución del código a las publicaciones científicas o los informes técnicos. En particular, el código y nuestras explicaciones o argumentos basados en ese código vivirían en ficheros separados. Por supuesto, si el código cambia, es muy posible que esas explicaciones y argumentos ya no sean válidos, con lo que las ocasiones para el error aumentan. Pero con Knitr eso ya no es necesario: el código y el texto van unidos en un mismo documento. Esa es una de las ideas centrales de la Programación Literaria, de la que hablaremos más extensamente en alguna otra entrada.

Figuras y control del código que se muestra.

El segundo bloque de código R que aparece en testRmarkdown01.Rmd es también interesante, porque ilustra:

• la forma de obtener un gráfico.
• cómo podemos controlar qué partes del código R se muestran y cuáles se ocultan.

Ese bloque es:

---

---

Fíjate en que para obtener la figura simplemente hemos incluido la función plot de R que produce esa figura. Insisto otra vez en lo que no hemos hecho: no hemos fabricado la figura con R para después hacer un corta/pega. Si fuera así y más adelante cambiamos algo en nuestro código, la figura ya no se correspondería con el código y habría que volver a pasar por el ciclo completo de generar la figura correcta, copiarla y pegarla. Pero si la figura se genera e inserta automáticamente a partir del código… ¿te das cuenta del avance que supone esta forma de trabajar?
Enseguida vamos a ver otro ejemplo para abundar en esas ideas. Pero antes vamos con el segundo punto que ilustra este segundo bloque de código del ejemplo. Si vuelves a mirar el documento html que ha producido Knitr, verás que el código plot(cars) no aparece por ningún lado (a diferencia de lo que ocurría con summary(cars), que sí era visible en el documento final). Para hacer invisible ese código, la cabecera del segudo bloque de código incluye la instrucción echo = FALSE de Knitr. Hagamos un experimento. Cambia echo = FALSEpor echo = TRUE y vuelve a pulsar el botón Knit de RStudio.

---

---

Ahora, además de la figura, puedes ver el código que la produce en el documento final. Además de echo, hay muchas otras opciones que podemos introducir en la cabecera de un bloque de código para controlar el comportamiento de Knitr. Por ejemplo, con eval=TRUE o eval=FALSE puedes decidir si un bloque de código se ejecuta o no. ¿Y si no quiero que un bloque se ejecute, no es más fácil no incluirlo o comentarlo? Puede parecer que sí, a primera vista. Pero ¿y si tedigo que puedes usar algo como eval=x, donde x es una variable booleana de R y que por lo tanto puedes decidir programáticamente qué partes del código se ejecutan y cuáles no? Es decir, que si algo cambia en tus datos y ya no quieres mostrar un cierto tipo de gráfico, puedes usar eval y echo para encender/apagar partes de tu documento. Y no hemos hecho más que rozar la superficie. Las posibilidades con Knitr son inmensas.

Un ejemplo más.

Vamos a añadir algunos ingredientes más al fichero testRmarkdown01.Rmd para ilustrar las ideas anteriores. Para empezar, añade estas líneas al final del documento (exactamente como aparecen aquí, puedes copiarlas y pegarlas):

## Un ejemplo adicional de gráfica con Knitr.

Vamos a mostrar una gráfica que podremos cambiar **simplemente** 
cambiando el valor de la variable `n`.

De momento se trata simplemente de un poco más de markdown, sin código R. Cuando hayas añadido esas líneas pulsa de nuevo el botón Knit para ver el resultado. La siguiente figura muestra los cambios que deberías ver en la parte final del documento html procesado.

---

---

Vamos con un bloque de código. El código es este:

```{r }
n = 1
x = 1:100
plot(x, x^n, type="l", lwd=4, col="red")
```

De nuevo, puedes copiarlo y pegarlo al final del documento. Tras hacer esto pulsa una vez más Knit. EL resultado debería ser cómo el de la siguiente figura:

---

---

Como ves, Knitr ha procesado el código y ha incorporado al documento la figura que se obtiene, que con n = 1 es una recta. Para insistir en la idea central de Knitr, cambia ahora ese valor por n = 2 y vuelve a pulsar Knit. Sin necesidad de ningún cambio adicional por tu parte, sólo con cambiar ese número, el resultado será el de la siguiente figura.

---

---

Como puedes ver, ahora con n=2 el gráfico es una parábola. Y déjame que insista, a riesgo de ser pesado: no hemos tenido que hacer nada más que cambiar ese número para que todo se actualice en el documento final.

Comentarios finales y enlaces recomendados.

Esta entrada sólo pretende presentar Knitr (y markdown) a quienes no los conozcan todavía y, con suerte, animaros a continuar indagando. En próximas entradas profundizaremos un poco más en el funcionamiento de Knitr. Pero para quien quiera empezar ya, hay dos enlaces que le pueden servir de puerta de entrada.

• La gente de RStudio ha publicado una página de introducción a la combinación de markdown y knitr que denominan Rmarkdown:
  http://rmarkdown.rstudio.com/
  Recomiendo encarecidamente el Quick Tour y la R Markdown Cheat Sheet.
• La página oficial de Knitr, que sirve de carta de presentación a su creador Yihui Xie, es http://yihui.name/knitr/. En ella encontraréis además un enlace al libro de Knitr, que os recomiendo si queréis profundizar en esta herramienta.

Ambos en inglés, me temo. No me convencen demasiado las (escasas) presentaciones de Knitr que he encontrado en español. Aunque, naturalmente, me encantaría que hubiera alguna mucho mejor (me ahorraría el trabajo de escribir estas entradas) y estoy abierto a cualquier sugerencia que los lectores quieran hacerme llegar.

Muchas gracias por la atención.

Compendio de comandos R (y otras cosas)

Esta entrada será breve. Vamos primero  con la parte más útil. A medida que se va aprendiendo a usar R, resulta conveniente disponer de un compendio o resumen de comandos, funciones, etc., como ayuda a la memoria. Hay varios de estos resúmenes (lo que en inglés de denomina a veces cheat sheet) publicados en la red, y yo suelo tener siempre a mano este en partiicular, elaborado por Tom Short, que me ahorra trabajo a la hora de localizar la sintaxis de un comando que se me resiste. Pero hay otros, que se pueden localizar fácilmente buscando en Google con R cheat sheet, o R reference card, etc.  Conviene mirar unos cuantos, localizar el que se adapte mejor a nuestras necesidades y conocimientos, y tenerlo a mano. No conozco ninguno en español que consiga empaquetar, en un espacio reducido, la cantidad de información de los que existen en inglés. Aprovecho para pedir que, si alguno de los lectores lo conoce, me lo haga saber, y lo reflejaremos aquí.

Las otras cosas: algo sobre planes futuros

La baja actividad del blog de estas semanas no se debe al verano. Los apuntes de las asignaturas de Estadística que imparto van camino de convertirse en un libro, escrito junto con el profesor M.Marvá. Este verano el trabajo se concentra en esa dirección. Queremos, además, escribir un libro que de alguna manera pueda seguir vivo una vez publicado, lo cual significa, entre otras cosas, que este blog se convierta en una extensión natural del libro, que irá acompañado además de un sitio Moodle, y de otros recursos en la red. Por eso, y para definir mejor la relación entre todos estos elementos, mientras el libro avanza el blog frena un poco el paso, pero no se para. En concreto, esta misma semana espero poder publicar la siguiente entrada sobre programación en R, que ya se ha hecho esperar bastante.

QQ-plots y funciones de distribución empírica en R.

Una advertencia antes de empezar: esta entrada es, probablemente, bastante más técnica que la media del blog. Si el lector entiende bien cómo se usa un QQ-plot, seguramente pueda ahorrarse la última parte, donde se intenta explicar cómo se construyen ese tipo de gráficos. Pero en cualquier caso, recomiendo leer la parte sobre funciones de distribución empírica de las muestras, porque esa es una herramienta interesante en si misma.

Y por si luego, metidos en faena, se nos olvida, una buena banda sonora para un discusión elevada como esta puede ser Skytalk, de Beats Antique.

Introducción

Ya hemos hablado de QQ-plots (quantile vs quantile plots) en el blog,en esta entrada sobre ANOVA unifactorial. Un QQ-plot es (aparte de algo que garantiza unas risas el primer día que lo mencionas en clase) una representación gráfica, que sirve para comparar dos distribuciones y ver si coinciden.  Allí vimos un ejemplo, que repito aquí:La idea de este tipo de gráficos es que, si las dos distribuciones coinciden, veremos un gráfico muy parecido a una recta. Tanto más parecido cuanto mayor sea la coincidencia.

Cuando hacemos Inferencia, a menudo suponemos que los datos de la muestra proceden de una distribución de cierto tipo, por ejemplo de una distribución normal.  Y en otras ocasiones, al contrario, hemos recogido una muestra y sospechamos que no responde a una distribución normal. ¿Cómo podemos verificar esto? Hay, por supuesto, contrastes diseñados para esta tarea. Pero  la exploración gráfica, como siempre, es una gran ayuda. Y un QQ-plot es una de las formas más sencillas de obtener información en estas situaciones.

En R, obtener un QQ-plot es muy sencillo, y lo veremos en detalle más abajo. Pero en lugar de conformarnos con esto, quiero aprovecharlo como excusa para introducir  la idea de función de distribución empírica de una muestra, y su tratamiento en R. La excusa, desarrollada por completo, va así: «¿si R no nos ofreciera el QQ-plot, cómo podríamos hacerlo a mano?». Para entender esto necesitamos esa idea de la función de distribución, que por un lado es una idea sencilla, y por otra parte es muy interesante por si misma. Veamos de qué se trata.

Función de distribución empírica de una muestra

En el fichero empiricalCDF-01.csv hay una muestra de datos, de tamaño n=10. Es una muestra pequeña, pero la elegimos así para que las gráficas que vamos a dibujar resulten más claras. Más abajo trabajaremos con muestras de mayor tamaño. 

Uno de los valores de esa muestra es 4.6550. Supongamos que queremos responder a esta pregunta: si elegimos un valor de la muestra al azar (de manera que todos los valores tienen la misma probabilidad de ser elegidos) ¿cual es probabilidad de que el valor elegido sea menor o igual que 4.6550?  El procedimiento – bastante sencillo – para responder, es este:

• Ordenamos los valores de la muestra. Se obtiene:
  

  x1 = 0.2849
  x2 = 0.5262
  x3 = 1.5600
  x4 = 2.7790
  x5 = 3.8230
  x6 = 4.6550
  x7 = 5.0150
  x8 = 5.3990
  x9 = 5.6570
  x10 = 9.9890
  

  Obsérvese que el valor 4.6550 ocupa la posición número seis en esta lista ordenada.

• Puesto que son equiprobables, la probabilidad de elegir cada uno de ellos es 1/n, es decir, 1/10. La probabilidad de elegir un valor menor o igual que x6=4.6550 es, por lo tanto, 6/10. En general, si el valor xk ocupa la posición número k en la lista, la probabilidad de elegir un valor menor o igual que xk es k/n. Diremos que

Una vez entendido esto, la función de distribución empírica de la muestra F(x) se define de forma muy parecida, extendiendo esa idea a un número cualquiera, que no tiene por que formar parte de la muestra. Concretamente,

Y, con menos notación pero tal vez más claridad, si x está, por ejemplo, entre el quinto y el sexto valor de una muestra (ordenada) de 17 valores, entonces F(x)=5/17.

En R disponemos de la función ecdf (empirical cumulative distribution function, ver aquí el fichero de ayuda de esa función), para obtener la función de distribución empírica de una muestra cualquiera. En nuestro caso, si el fichero  empiricalCDF-01.csv  está en el Escritorio de la máquina, para leer los datos, definir y la función y dibujarla (con algunos adornos), ejecutaríamos este código:

# En Windows, las siguientes lineas
# establecen el Escritorio del usuario actual como directorio de trabajo

(dirTrabajo=paste("C:/Users/",Sys.info()["user"],"/Desktop",sep=""))
setwd(dirTrabajo)

# En sistemas Mac o Linux alguna de estas lineas (descomentar) debería tener el mismo efecto.

# setwd("~/Desktop")
# setwd("~/Escritorio")

(datos=scan(file="empiricalCDF-01.csv",sep=" "))

# Aunque en general no es necesario, en este ejemplo vamos a reordenarlos

(datos=sort(datos))

# Ahora definimos y dibujamos la gráfica de la función de distribución empírica:

ECDF=ecdf(datos)
plot(ECDF,col="red",lwd=3,xlab="datos",ylab="")
points(datos,rep(0,n),col="blue",pch=20,cex=2)
segments(datos,rep(0,n),datos,sapply(datos,ECDF),col="blue",lwd=3)

Como puede verse, hemos llamado ECDF a la función. No es una buena idea llamarla F, porque R interpreta ese símbolo como una abreviatura del valor lógico FALSE, y se pueden generar errores inesperados. La gráfica que se obtiene es esta: Como se ve, es una gráfica en forma de escalera, ascendente de izquierda a derecha, con un escalón (todos de altura 1/n) sobre cada uno de los valores de la muestra, que sube desde 0 hasta 1 a medida que recorremos los valores de la muestra. Estas son características generales de cualquier función de distribución empírica (entrada en la Wikipedia, en inglés.).

La función ECDF que hemos definido en R es, por supuesto, un objeto de tipo function (usar class(ECDF) para comprobarlo, y ver esta entrada si se necesita más información sobre funciones). Así que podemos calcular ECDF(x) para cualquier valor x. ¿Cómo interpretar estos valores ECDF(x)? Sabemos como hacerlo para los valores de la muestra, pero ¿y para un valor x cualquiera? Si calculamos ECDF(x) y obtenemos 0.36, podemos interpretar que el 36% de los valores de la muestra son menores o iguales que x. Visto así, esta noción está estrechamente emparentada con la idea de percentil (o cuantil).

¿Qué tienen que ver estas funciones de distribución empíricas con las funciones de distribución teóricas de una variable aleatoria (de las que ya hemos hablado en la duodécima sesión)? Recordemos que si X es una variable aleatoria (discreta o continua), su función de distribución se define como:

Es decir, la probabilidad (según la distribuye la variable X) de la cola izquierda definida por el valor x. Por ejemplo, la siguiente gráfica muestra la función de distribución de una variable Z de tipo normal estándar N(0,1) que, en R, no es otra que nuestra vieja amiga pnorm. 

Como se ve, una función que sube, de izquierda a derecha, desde el 0 hasta el 1. ¿Qué sucede si tomamos una muestra de, pongamos 200 valores de una distribución normal , los tipificamos, después los usamos para definir una función de distribución empírica, y finalmente la superponemos sobre esta teórica? Para hacerlo, basta con ejecutar este código:

set.seed(2013)

#Generamos unos datos para que sirvan de ejemplo de cdf
n=200
datos=rnorm(n)
(xbar=mean(datos))
(desvt=sd(datos))
datosNormalizados=(datos-xbar)/desvt

ft=function(x){pt(x,df=2)}

plot(pnorm,from=-5,to=5,col="red",lwd=3)

plot(ft,from=-5,to=5,col="magenta",lwd=3,add=TRUE,lty=2)
ECDF=ecdf(datosNormalizados)
plot(ECDF,col="blue",lwd=1,xlab="datos",ylab="",add=TRUE)

Hemos incluido la orden set.seed(2013) para que el lector obtenga exactamente la misma muestra, y por tanto los mismos resultados, que nosotros. Si se desea obtener muestras normales distintas, comentar esa línea. Además, hemos incluido (a trazos) la gráfica de la función de distribución teórica de una t de Student con dos grados de libertad, para que quede claro que, en este caso, hay una clara diferencia con la muestra. El resultado es esta figura: 

Y como se ve, la distribución empírica, de una muestra de este tamaño, sigue muy fielmente a la distribución teórica. El mensaje que hay que extraer de este ejemplo es que, en general, para muestras grandes, la función de distribución empírica caracteriza bastante bien a la distribución teórica subyacente. Y como se trataba de tener una forma de comparar distribuciones, esto nos deja a las puertas del QQ-plot.

Cómo construir un QQ-plots en R «a mano»

Empecemos con una muestra (ordenada)

suficientemente grande, que procede de una distribución normal . La tipificación consiste en calcular:

Es decir, que desde el punto de vista matemático, es una transformación de la forma:

cuya gráfica de y frente a x es una recta. Es importante entender esto, porque esa es la razón última por la que en el QQ-plot esperamos ver algo parecido a una recta, si hemos partido de una muestra normal. Hemos empezado con los valores x, pero ahora, después de tipificar, tenemos unos valores

que podemos considerar como una muestra de Z (la normal N(0,1)). Eso significa que su función de distribución empírica es la escalera determinada por los puntos:

Y hemos visto (esto es muy importante, revisar la última gráfica) que la escalera que forman estos puntos se parece mucho a la función de distribución teórica de Z, y que esa distribución teórica es lo que en R llamamos pnorm.

La función qnorm es la inversa de pnorm. Ahora, si empezamos con las probabilidades:

y calculamos los valores:

, ,…,

entonces es de esperar que los valores

sean muy parecidos a los .

Esperamos que el lector no haya sufrido demasiado, y no se haya perdido en la discusión. Es un poco complicada, pero ya hemos llegado. Si los se parecen mucho a los , y la gráfica de los frente a los   era una recta, entonces la gráfica de los  frente a los  se debe parecer mucho a una recta. Esta última gráfica es nuestro QQ-plot «artesanal».

A continuación incluyo el código en R para este proceso que hemos descrito. Pero antes, es necesario hacer algunos comentarios de advertencia: el código usa, además,  la función qqnorm de R, para comparar el resultado con el proceso artesanal que hemos descrito. Y hay buenas razones para que ambos no coincidan exactamente. En particular, nosotros hemos usado 1/n,…,n/n como probabilidades, mientras que R usa un conjunto distinto (ver el Capítulo 2 del libro Testing for Normality, de H.C.Thode. ) Además, R coloca la muestra original en el eje vertical (es la práctica habitual en un QQ-plot), y nosotros lo hemos colocado en el horizontal porque, de lo contrario, la discusión anterior habría sido aún más complicada de entender… eso nos obliga a hacer algunas manipulaciones a partir del resultado de qqnorm, para colocar las cosas como queremos.

set.seed(2013)
#Generamos una muestra con datos de una distribucion normal,
n=200
x=rnorm(n,mean=7,sd=1.5)

# la ordenamos,
x=sort(x)

# calculamos su media y cuasidesviacion tipicas,
(xbar=mean(x))
(desvt=sd(x))

# y hacemos una especie de tipificacion.
y=(x-xbar)/desvt

# Ahora empezamos con las probabilidades,
probs=(1:n)/n

# y las invertimos usando qnorm.

y2=qnorm(probs)

# El QQ-plot artesanal es la grafica x frente a y2.
# Para comparar con lo que hace R, usamos qqnorm

qqplot=qqnorm(x,plot.it=FALSE)

# Ademas R invierte los ejes (cambia x por y)
# asi que tenemos que colocarlos al reves para poder comparar
qqx=qqplot$x
qqy=qqplot$y

# En azul el QQ-plot artesanal, y en rojo el de R.
plot(x,y2,col="blue",pch=3)
points(qqy,qqx,col="red") # Observar el cambio de orden en los ejes

Y el resultado es esta gráfica. en la que el resultado de R se muestra en rojo, y el nuestro en azul. Como puede verse, ambos son muy similares, y ambos además se parecen mucho a una recta, como esperábamos al haber partido de una muestra normal.

Muchas gracias por la atención, especialmente en una entrada larga y seguramente más técnica que la media.

RStudio

En las entradas precedentes del blog  he usado, hasta ahora, R-Commander como interfaz que simplifica nuestra relación con R, especialmente al principio. Y sigo pensando que es un buen programa para principiantes, y continuaré explicando las opciones que ofrece.  Tengo además una deuda de gratitud con ese programa, porque mis primeros pasos con R habrían sido mucho más difíciles si no hubiera descubierto R-Commander (y el paquete R-UCA de la Universidad de Cádiz).

Sin embargo, en los últimos tiempos me he dejado convencer por otra interfaz gráfica de R, que se llama R-Studio (el enlace es la página principal del proyecto). Para presentarlo, aquí tenemos una captura de pantalla de RStudio, en una sesión de trabajo con el código de la entrada sobre ANOVA para diseños factoriales (con tres factores en este caso).

RStudio ofrece muchas herramientas para el usuario avezado de R, convirtiéndose en un entorno de trabajo muy productivo. Esa es la mayor diferencia con R-Commander, que tal vez está más orientado a proporcionar al usuario una experiencia  similar a la de los programas comerciales clásicos, donde los análisis están disponibles mediante un sistema de menús y cuadros de diálogo, como este:

En RStudio no encontraremos nada de eso. No hay un menú que nos diga cómo hacer un contraste de hipótesis, un intervalo de confianza. Para eso está R. En ese sentido, RStudio da por sentado que nosotros sabemos pedirle a R que haga esos análisis. Y por eso, para un principiante, puede resultar algo más intimidatorio.

Para el usuario con algo más de experiencia, como decíamos, la situación se invierte, y RStudio se muestra como una de las mejores opciones para el trabajo diario con R. Entre otras, quiero destacar estas características:

• la posibilidad de trabajar con varios ficheros de comandos R simultáneamente, agrupados por pestañas (en la ventana superior izquierda, en la anterior captura de pantalla).
• el editor de comandos es, en si mismo, una de las grandes virtudes de RStudio. Es realmente un editor pensado por y para programar en R. No echaremos en falta un editor externo, al que en general teníamos que recurrir con R-Commander.
• autocompletado de comandos. Si no recordamos la sintaxis exacta de un comando de R, basta con escribir las primeras letras, pulsar la tecla Tab, y RStudio nos ofrece información sobre ese comando. Pero también sobre las variables, ficheros y demás objetos que hayamos creado.
• gestión de los gráficos de R mucho más eficaz que, por ejemplo, en R-Commander.
• facilidad de gestión para, por ejemplo, instalar paquetes, consultar ayuda, etcétera.
• si se dispone de varias versiones de R instaladas en la máquina, es fácil seleccionar la que utilizará en cada sesión.

Pero además, RStudio ofrece la posibilidad de instalación en un servidor (bajo Linux), para posteriormente acceder al programa a través de una interfaz web, desde cualquier máquina, aunque no esté instalado R. Esa interfaz web es esencialmente idéntica a la que hemos mostrado más arriba, como puede verse en esta figura:

Esto hace posible (aunque desde luego, no muy cómodo) acceder a RStudio desde un teléfono móvil. Y desde luego es especialmente interesante para usuarios de tablets, (tanto Ipads como Android).

Para finalizar esta primera visita, en la página del proyecto RStudio se puede ver otro producto, Shiny, del que hablaremos en una futura entrada y que permite convertir los análisis que hacemos con R en páginas web interactivas, para que cualquier pueda experimentar con esos análisis, y modificar los parámetros utilizando deslizadores, cajas de diálogo, etc.  Una muestra, en este enlace.

Shiny  es una forma de avanzar en la idea de mostrar a los demás nuestro trabajo, cada vez de forma más trasparente. Una idea que está en la base del movimiento hacia lo que se ha llamado Investigación reproducible (Reproducible Research, página en la Wikipedia). En esencia, la idea es que los trabajos de investigación, a la hora de publicar sus resultados, debería acompañar la publicación tradicional con el conjunto completo de datos experimentales y los ficheros de código del software que se haya utilizado para el análisis de esos datos.

R tiene una herramienta muy poderosa para contribuir a esta idea, llamada SWeave. Y en esta entrada no puedo dejar de mencionar que RStudio (a diferencia de R-Commander) ofrece mucha ayuda para el trabajo con SWeave. Pero ese un tema que merece ser tratado con mucho más detenimiento. En una próxima entrada volveremos sobre la investigación reproducible (y sus conexiones con la política y la economía, nada menos), sobre SWeave, y sobre la noción de Programación Literaria, que hizo posible todas estas ideas. Y (spoiler alert para profesores) sobre la forma de conectar todo esto con Moodle y usarlo en nuestras clases.

Gracias por la atención.

Coursera: vídeos del curso «Computing for Data Analysis» y otros cursos interesantes

Ya he hablado aquí alguna vez de Coursera (desde mi punto de vista, un acontecimiento esencial en la evolución de la educación superior).  Y hoy tengo dos razones para volver a hablar de ellos, aunque sea brevemente.

En primer lugar, la gente del blog Revolutions se ha tomado la molestia de recopilar enlaces a todos los vídeos del curso Computing for Data Analysis, que dirige en esa plataforma el profesor Roger Peng, de la John Hopkins University. El curso se ha impartido ya (no sé si varias veces), pero está prevista una nueva iteración en Septiembre de este año. Para quienes han tenido ya un contacto inicial con R, y desean profundizar en esta herramienta (¿se atreverá alguno de mis alumnos?), será  una buena experiencia, y con la posibilidad añadida de ganar un certificado de una universidad prestigiosa.

No es el único curso interesante que se abrirá en próximas fechas. Por mencionar sólo algunos (que han despertado mi interés, claro; hay muchos más):

• Computational Molecular Evolution (empieza el 24 de Junio de 2013).
• Network Analysis in Systems Biology (empieza el 13 de Octubre de 2013).
• Dynamical Modeling Methods for Systems Biology (empieza el 13 de Octubre de 2013).
• Mathematical Biostatistics Boot Camp (empieza el 24 de Septiembre de 2013).

Este último lo puedo recomendar con conocimiento de causa, porque participé en una edición anterior. Una experiencia intensa en trabajo; que nadie se crea que estos cursos son una broma, las horas de trabajo que indican son absolutamente necesarias en la mayoría de ellos, y en ocasiones se quedan cortas. Y se aprende, desde luego, tanto si uno es profesor, como si es alumno.

En el futuro seguiré escribiendo entradas para avisar de los cursos que me parezcan interesantes y estén relacionados con la temática del blog.

Tablas de contingencia y contraste de independencia (test chi cuadrado)

Contraste de independencia

El Barómetro del CIS  (Centro de Investigaciones Sociológicas) permite obtener datos sobre las creencias religiosas de la población en España. Una pregunta que puede interesarnos es ¿hay alguna diferencia al respecto entre hombres y mujeres? Vamos a utilizar los datos del Barómetro para intentar contestar esa pregunta.

Por ejemplo, en el mes de enero de 2013 el Barómetro} recoge las respuestas de n=2452 personas sobre sus creencias religiosas. (En realidad son 2483, pero para simplificar vamos a eliminar de nuestra consideración a las 19 mujeres y a los 12 hombres que decidieron no contestar). Observa que, como de costumbre, vamos a usar n para el número total de personas encuestadas. Agrupamos a todos los creyentes de distintas religiones por un lado y a los que se declaran no creyentes o ateos por otro. Y así tenemos una tabla de doble entrada:

	Hombres	Mujeres	Total
Creyentes	??	??	1864
No creyentes	??	??	588
Total	1205	1247	2452

En la que hemos dejado pendientes los valores interiores de la tabla, y sólo se muestran los valores marginales. Si las creencias religiosas no dependieran del género, ¿cuáles son los valores que esperaríamos ver en esta tabla? Los que se deducen de los marginales, claro. Por ejemplo, el número de hombres creyentes se puede calcular multiplicando el número de hombres por la proporción de creyentes en la población. Esto es:

Con R podemos obtener estos valores esperados fácilmente. Como es la primera vez que hacemos algo así en el blog, vamos a dar los detalles. Empezamos con dos vectores, los de los valores marginales:

creencia=c(1864,588)
genero=c(1205,1247)

Y ahora los multiplicamos como matrices, trasponiendo el vector creencia, que va en vertical en la tabla, y que por tanto hace el papel de vector columna. Para trasponer un vector basta con usar la función t(). El producto matricial se representa mediante el símbolo %*%, así que la tabla (matriz producto) es (mostramos la salida):

> creencia %*% t(genero)
        [,1]    [,2]
[1,] 2246120 2324408
[2,]  708540  733236

Naturalmente, aún tenemos que dividir por el total de personas (el número n). El resultado completo es esta tabla de valores esperados:

> (tablaEsperada=creencia %*% t(genero) / sum(genero))
         [,1]     [,2]
[1,] 916.0359 947.9641
[2,] 288.9641 299.0359

Está bien, pero sería más fácil de interpretar si añadimos nombres a las filas y columnas, y además incluimos los valores marginales con los que habíamos empezado. Hacemos esto con los comandos:

colnames(tablaEsperada)=c("H","M")
rownames(tablaEsperada)=c("CREY","NO_CREY")
(addmargins(tablaEsperada))

y el resultado es mucho más agradable de ver:

                H         M  Sum
CREY     916.0359  947.9641 1864
NO_CREY  288.9641  299.0359  588
Sum     1205.0000 1247.0000 2452

De acuerdo, eso es lo que esperaríamos, si las creencias y el género fueran independientes. ¿Cuáles son los datos reales que se obtuvieron en el Barómetro? Aquí está la tabla de valores observados:

	Hombres	Mujeres	Total
Creyentes	849	1015	1864
No creyentes	356	232	588
Total	1205	1247	2452

Introducimos esta tabla en R directamente usando cbind, y le añadimos la información necesaria para hacerla más útil:

(tablaObservada=rbind(c(849,1015),c(356,232)))
colnames(tablaObservada)=c("H","M")
rownames(tablaObservada)=c("CREY","NO_CREY")
(addmargins(tablaObservada))

El resultado es:

           H    M  Sum
CREY     849 1015 1864
NO_CREY  356  232  588
Sum     1205 1247 2452

Comparando las dos tablas, resulta evidente que los valores observados no coinciden con los  esperados. De hecho, el número de hombres no creyentes es más alto de lo que habíamos estimado a partir de la población en conjunto (y,  lógicamente, el número de mujeres no creyentes es más bajo que la estimación). Pero ese número de hombres no creyentes, ¿es significativamente más alto?

Para contrastar la hipótesis de independencia entre género y creencias, usaremos la función chisq.test de R (se muestra la salida):

> (chisq.Observada=chisq.test(tablaObservada,correct=FALSE))

	Pearson's Chi-squared test

data:  tablaObservada
X-squared = 40.2253, df = 1, p-value = 2.263e-10

Como en otros casos, hemos incluido la opción correct=FALSE para que R (lo haga un poco peor de lo que sabe y) nos devuelva el mismo valor que aparecería en un curso de Introducción a la Estadística.  Se recomienda leer la ayuda de chisq.test para ver lo que implica esa opción.

En este caso el p-valor es muy pequeño, así que con esos datos deberíamos rechazar la hipótesis nula, y concluir que el género y las creencias no son independientes.

El test que estamos usando consiste, para este ejemplo, en calcular el siguiente estadístico:

donde los valores son los de la tabla observada, y los son los valores esperados. Cada término de la forma

es un residuo, y el estadístico es la suma de los residuos al cuadrado. En R,  al aplicar el test, el programa devuelve,  dentro de la variable chisq.Observada, una lista con todo lo  necesario. Accedemos a las componentes de esa lista con el operador $. Por
ejemplo, las tablas de valores esperados y observados se obtienen así (se muestra la salida):

> (E=chisq.Observada$expected)
               H        M
CREY    916.0359 947.9641
NO_CREY 288.9641 299.0359

> (O=chisq.Observada$observed)
          H    M
CREY    849 1015
NO_CREY 356  232

Y a partir de aquí el cálculo de los residuos es inmediato:

> (residuos=(O-E)/sqrt(E))
                H         M
CREY    -2.214885  2.177266
NO_CREY  3.943532 -3.876553

Lo hemos hecho de esta manera para que quede claro cómo se definen, pero los residuos también se pueden obtener directamente con:

chisq.Observada$res

De cualquier manera, para calcular a mano el estadístico a partir de los residuos, podemos ahora elevar los residuos al cuadrado y sumarlos usando addmargins, de esta manera:

> (tablaChi2=residuos^2)
                H         M
CREY     4.905714  4.740486
NO_CREY 15.551448 15.027663

> addmargins(tablaChi2)
                H         M      Sum
CREY     4.905714  4.740486  9.64620
NO_CREY 15.551448 15.027663 30.57911
Sum     20.457163 19.768148 40.22531

El estadístico es el valor de la esquina inferior derecha de esta tabla:

> (estadistico=addmargins(tablaChi2)[3,3])
[1] 40.22531

Y podemos ver que coincide con lo que obtuvimos antes usando chisq.test. El p-valor se obtiene a partir del estadístico usando pchisq como en cualquier contraste (con un grado de libertad, porque es una tabla 2×2). Usamos la cola derecha, porque los valores del estadístico que contradicen a la hipótesis nula se sitúan ahí:

> (pValor=1-pchisq(estadistico,df=1))
[1] 2.262972e-10

Y como antes, el p-valor coincide con lo que obtuvimos en chisq.test.

Una posible representación gráfica de este contraste es mediante un gráfico de columnas apiladas (la altura indica el porcentaje), con una columna por  género. Para obtener esa representación empezamos por fabricar una tabla de proporciones (porcentajes), en la que los porcentajes son sobre el total de la muestra. Añadimos las sumas marginales para que sea más fácil entender esta tabla (recordamos que se muestra la salida):

> tablaProporciones1=prop.table(tablaObservada)
> addmargins(tablaProporciones1)
                H          M       Sum
CREY    0.3462480 0.41394780 0.7601958
NO_CREY 0.1451876 0.09461664 0.2398042
Sum     0.4914356 0.50856444 1.0000000

Para que cada columna de la tabla nos de las alturas de las columnas del gráfico, tenemos que dividir la tabla por las sumas por columnas, calculadas con colSums. Además, vamos a redondearla para mejorar la presentación.

> (tablaProporciones2=round(tablaObservada/colSums(tablaObservada),digits=1))
          H   M
CREY    0.7 0.8
NO_CREY 0.3 0.2

Ya podemos usar barplot para dibujar el grafico, con colores y etiquetas en los ejes. Aunque con esto ya se obtiene un gráfico adecuado, la interpretación de un gráfico como este es mucho más fácil, si se rotulan las columnas con el porcentaje de una de las categoría en cada columna. Primero fabricamos los rótulos, y luego los colocamos en sus posiciones con text. Estos son los comandos:

(bp_t2=barplot(tablaProporciones2,col=c("lightseagreen","darkgoldenrod1"),xlab="Género",ylab="Creyente/No creyente",main="Creencias religiosas por género",font=2))

rotulos=paste(format(100*signif(tablaProporciones2[1,],digits=2)),"% creyentes",sep="")

text(bp_t2, tablaProporciones2[1,]+0.05,rotulos , xpd = TRUE, font=2)

Y este es el resultado:

Un problema como este, con una tabla de contingencia 2×2, también se puede abordar con un contraste de igualdad de proporciones para dos muestras, como los que vimos en esta entrada. Vamos a hacerlo aquí (de nuevo, con correct=FALSE), para que se vea que R devuelve el mismo valor del estadístico, y el mismo p-valor.

> (creyentesObs=tablaObservada[1,])
   H    M
 849 1015
> prop.test(creyentesObs,genero,correct=F)

	2-sample test for equality of proportions without continuity correction

data:  creyentesObs out of genero
X-squared = 40.2253, df = 1, p-value = 2.263e-10
alternative hypothesis: two.sided
95 percent confidence interval:
 -0.14300581 -0.07577254
sample estimates:
   prop 1    prop 2
0.7045643 0.8139535

Por último, si tenemos dudas sobre las hipótesis de normalidad subyacentes al test , podemos hacer un test no paramétrico para confirmar lo anterior.

> fisher.test(tablaObservada)

	Fisher's Exact Test for Count Data

data:  tablaObservada
p-value = 2.745e-10
alternative hypothesis: true odds ratio is not equal to 1
95 percent confidence interval:
 0.4490685 0.6611693
sample estimates:
odds ratio
 0.5452314

que, en este ejemplo, confirma que el resultado del test es significativo.

Código en R y observaciones finales

En este enlace hay un fichero con el código R que hemos usado en este ejemplo. Naturalmente, el caso 2×2 no agota las posibilidades del test de independencia. Estoy terminando de escribir el correspondiente capítulo de los apuntes de teoría (avisaré en un post cuando acabe), que incluye ejemplos y código para el caso general, con tablas de contingencia con un número mayor de filas y columnas.

Aparte de esto, y como se explica en ese capítulo, existe otra variante del test, que se aplica a problemas de bondad de ajuste, es decir a la pregunta de si los datos observados se corresponden con una distribución teórica. Concretamente, ejemplos de este tipo de preguntas son: ¿cómo podemos saber si un dado está cargado? O ¿corresponden estos datos a una población normal?.  Próximamente le dedicaremos una entrada del blog a este tema. Y, puesto que aquí ha surgido, siquiera sea de pasada, el tema de las operaciones con matrices, también le dedicaremos una entrada a ese tema.

Gracias por la atención.

Y una última novedad (si P. Krugman se atreve, no voy a  ser yo quien se quede atrás): Recomendación musical de esta entrada Headlock de Imogen Heap,

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:

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: 

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:

Se 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:

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.