Una animación con R

La idea es sencilla. Se trata de componer un vídeo mediante la unión de imágenes como en una película de dibujos animados.
En este ejemplo, vamos a utilizar el parámetro de amplitud h de una función coseno

f(t) = h*cos(omega(t-rho))

con omega=1 y rho=0 fijos.

En primer lugar fijamos los parámetros y hacemos un gráfico preliminar:

h <- 1
rho <- 0
omega <- 1

x <- seq(-1.5*pi+0.5, 2*pi+0.5, length=250)
y <- h*cos(omega*(x-rho))
plot(x, y, type="l", ylim=c(-1.5,1.5))
abline(h=0, v=0, col="grey75")



Ahora vamos a crear una secuencia de imágenes variando el parámetro h entre los valores 1 y 1.5 primero, después bajándolo hasta 0.5 y luego subiéndolo otra vez hasta 1.
Para un vídeo de 8 segundos necesitamos al menos 8*25=200 imágenes.

hnew <- c(seq(1, 1.5, by=0.01),
seq(1.48, 0.5, by=-0.01),
seq(0.52, 1, by=0.01))

También debemos crear una carpeta donde guardaremos todas las imágenes separadas del resto de documentos. La carpeta debe existir realmente para que la podamos utilizar.

# imagesDir <- "C:/Mis documentos/R/video"
imagesDir <- "/home/francesc/R/video"

A continuación abrimos el dispositivo para crear múltiples imágenes en formato PNG, dentro de la carpeta escogida y con una numeración de tres dígitos correlativa.

png(file.path(imagesDir, "Rplot%03d.png"))

Ahora ya podemos generar las imágenes con un for.

for (i in 1:length(hnew)){
plot(x, y, type="l", ylim=c(-1.5,1.5),
main=expression(f(t)==h*cos(omega*t-omega*rho)), lwd=2)
abline(h=0, v=0, col="grey75")
h <- hnew[i]
lines(x, h*cos(omega*(x-rho)), col=2)
lines(c(0,0),c(0,h), col=2, lwd=2)
lines(c(-0.05,0.05),c(h,h), col=2, lwd=2)
lines(c(-0.05,0.05),c(0,0), col=2, lwd=2)
lines(c(0,1.8), c(h/2,0.85), col=2, lwd=2)
text(1.9,0.87, paste("h=", round(h,2), sep=""), adj=0, col=2)
text(1.9,1.1, "Changing the amplitude:", font=2, adj=0, col=2)
}
dev.off()

La última instrucción cierra el dispositivo png.
El último paso es crear el vídeo con algún programa específico. Una opción es ImageMagick que se puede ejecutar desde la consola con la instrucción:

convert /home/francesc/R/video/*.png video.mpeg

En este caso no ha funcionado (otras veces sí) y he utilizado el programa OpenShot que permite crear vídeos con la importación de una secuencia de imágenes con el mismo nombre que nosotros le hemos dado: Rplot%03d.png.

¿Os animáis a crear vuestros vídeos?

Sweave desde un editor de texto

En algún otro artículo ya he resaltado la importancia de utilizar Sweave para mezclar código de R con texto y producir informes automáticos con resultados estadísticos y gráficos dinámicos. Sweave se puede utilizar desde LibreOffice, no sin ciertas dificultades, con el paquete odfWeave, pero sin duda los mejores resultados se obtienen con el texto en formato LaTeX.

Supongamos que tenemos una buena instalación de R y de algún programa de LaTeX para nuestro sistema operativo. Entonces se puede ejecutar Sweave desde la consola de R o incluso desde una GUI como RStudio, pero es mucho mejor ejecutarlo desde el propio editor de LaTeX que seguramente es donde hemos escrito el documento .Rnw que queremos compilar. Este es el objetivo: ejecutar Sweave sobre el documento que estamos editando sin salir del editor.

Sweave en un sistema Windows

Supongamos que utilizamos Texmaker con la distribución MikTeX de LaTeX, aunque el procedimiento es similar si utilizamos otro editor y otra distribución de LaTeX.
Sweave viene con la distribución estándar de R y se puede ejecutar desde su consola, pero yo prefiero compilar el documento directamente desde el editor. Para ello necesitamos un archivo llamado Sweave.bat que se puede obtener en el CRAN dentro del archivo batchfiles_x.x-x.zip. Para instalar el archivo Sweave.bat debemos dejarlo en alguna carpeta que sea accesible en nuestro PATH o añadir dicha carpeta al PATH.
Ahora ya deberíamos ser capaces de ejecutar este .bat desde la consola, sólo necesitamos indicarle al editor donde encontrarlo.
Abrimos Texmaker y seleccionamos Opciones -> Configurar Texmaker e incluir la orden:

Sweave.bat --pdf %.Rnw

en alguna de las categorías, por ejemplo como Montaje rápido.

Un detalle importante: en Windows hay una dificultad añadida ya que necesitamos que LaTeX sea capaz de hallar el archivo Sweave.sty. Para ello copiaremos este archivo desde

C:\Archivos de programa\R\R-x.x.x\share\texmf

en alguna carpeta de la distribución de MikTeX y luego refrescaremos la base de datos.

Sweave en un sistema Mac OS X

Los pasos para poner a punto un ordenador Mac OS X son incluso más sencillos que los explicados para un ordenador con Windows. Necesitamos el script Sweave.sh que hallaremos en el CRAN.
Para instalarlo, debemos copiarlo en

/usr/local/bin

y, dentro del programa Terminal, escribimos

sudo chmod +x /usr/local/bin/Sweave.sh

para hacerlo ejecutable.

Ahora ya podemos ejecutar este script desde la consola. Sin embargo, para que podamos ejecutarlo desde TEXshop, sólo hay que indicarle a TEXShop donde hallar dicho Sweave.sh.
Abrimos TEXShop y buscamos TEXShop→Preferences. Hacemos click en la pestaña Misc y escribimos

/usr/local/bin/Sweave.sh −ld

en el campo LaTeX Progam Personal Script.

Sweave en un sistema Linux

Para poner a punto un ordenador con Linux, el procedimiento es también muy sencillo, aunque los detalles pueden variar ligeramente dependiendo de la distribución y el editor utilizados. El ejemplo que explicaré se basa en Kubuntu Linux y el editor es Kile.
Para ejecutar Sweave desde el editor necesitamos el script Sweave.sh que podemos hallar en el CRAN. Luego instalaremos ese script en

/usr/local/bin

entonces abrimos la Konsole y escribimos

sudo chmod +x /usr/local/bin/Sweave.sh

para hacerlo ejecutable.
El siguiente paso es indicarle a Kile donde hallar este script. Abrimos Kile y hacemos click en Settings→Configure Kile. Nos vamos a Tools→Build y añadiremos una nueva herramienta que podemos llamar Sweave. En la ventana superior escribiremos la orden Sweave.sh y

−ld '%source'

en la de abajo.
Recordemos que el archivo Sweave.sty debe ser accesible.

Más información en el artículo:
http://www.scribd.com/doc/6451985/Learning-to-Sweave-in-APA-Style

Sweave, FactoMineR y PDFs

El asunto es que estoy escribiendo un informe en LaTeX con Sweave. Lo hago así porque de esta forma me aseguro que añadir una variable o cambiar un dato no será un problema. El informe se volverá a generar en un minuto con todos los cambios.
El análisis que hago en concreto es un Análisis Factorial Múltiple o MFA (Multiple Factor Analysis) con el paquete FactoMineR. Se trata de un análisis de reducción de la dimensión cuando tenemos varios grupos de variables, tanto cuantitativas como cualitativas. Además este análisis permite estudiar los ejes principales mediante algunos gráficos y, como era de esperar, la representación de los individuos con las componentes principales.
En concreto utilizo PDFLaTeX directamente con el archivo .tex que resulta de procesar el archivo .Rnw con el Sweave. Así, los gráficos deben estar en formato PDF o PNG y no EPS. Ahora con Sweave esta es la opción por defecto.
Bueno, el caso es que los gráficos generados por el Sweave estaban vacíos. Eran archivos PDF con el nombre habitual por defecto, pero vacíos. ¿Sería un bug de Sweave? ¿Un error de mi instalación de R? ¿Un problema al cambiar de carpeta los archivos PDF?

Pues me costó un buen rato descubrir la razón del problema. Como siempre, lo primero fue aislar el error. No era un problema de Sweave, ni de mi instalación de R, ni de las opciones del Sweave (cambio de carpeta para las imágenes). Tuve que hacer unas cuantas pruebas, para descubrir que el problema estaba en los gráficos que generaba el paquete FactoMineR. En mi caso la función plot.MFA(). Tal vez, si hubiera leído mejor la ayuda de la función...

Bien, pues aquí va la explicación:
Resulta que para generar un gráfico en PDF, debemos desviar el dispositivo de salida gráfica hacia un archivo PDF. Eso se hace manualmente con la instrucción pdf() y seguramente lo debemos cerrar con la instrucción dev.off(). Eso es lo que hace automáticamente el Sweave.
Resulta que leyendo el código de la función plot.MFA() observé que, por defecto, hace una llamada a un nuevo dispositivo con la instrucción dev.new(). Ese era el problema.
Si abrimos el dispositivo pdf() y luego otro dispositivo con dev.new(), entonces el gráfico va al segundo dispositivo y al cerrarlo con dev.off(), cerramos el segundo y no el primero. La gente de FactoMineR tenía el problema previsto con un parámetro new.plot=TRUE que para mi gusto debería estar en FALSE por defecto. Ellos lo han puesto en TRUE para que podamos ver unos cuantos gráficos sin machacar el anterior, pero eso en mi caso era mortal.
En resumen, con la opción new.plot=FALSE, he podido generar el informe con PDFLaTeX sin problemas.
Moraleja: si pensáis en una función que genere varios gráficos seguidos, hay que pensar en un parámetro que regule la presentación uno a uno y que, mejor por defecto, no abra un nuevo dispositivo para cada uno.

RKWard 0.5.7 para el nuevo R 2.14.0

Los desarrolladores de RKWard han lanzado una nueva versión (23 de octubre de 2011) de este magnífico software. Lo que me ha pasado es que al actualizar R a su versión 2.14.0, el RKWard dejó de funcionar. Pero como era de esperar, los amigos de RKWard ya estaban preparados.
Como yo utilizo Kubuntu (natty) he buscado el repositorio adecuado para actualizar RKWard que es

deb http://ppa.launchpad.net/rkward-devel/rkward-stable-cran/ubuntu natty main

y automáticamente se ha actualizado sin problema ninguno. Como veis todavía sigo con la versión natty o 11.04 de Kubuntu, ya que me gusta esperar algunos días antes de pasar a la nueva 11.10.

Además de nuevos detalles y corrección de errores, esta versión de RKWard viene con una importante mejora en la forma de añadir complementos, ya que ahora estos complementos se instalan como paquetes de R.

Aunque hay una versión para Windows, yo no la recomiendo ya que requiere un escritorio KDE para Windows y todavía no está demasiado desarrollada. Una opción alternativa y que comentaré próximamente es el programa multiplataforma RStudio que viene pegando fuerte.

El código de una función

Una de las ventajas de utilizar un programa como R de código libre es la posibilidad de aprender a través de la lectura del código de las funciones incorporadas, ya sea de los paquetes oficiales del CRAN o de los que proporcionan todas las personas que contribuyen de forma particular.

La forma más sencilla de acceder al código de una función es con la llamada a su nombre (sin paréntesis). Por ejemplo, para ver el código de la función matrix():

> matrix
function (data = NA, nrow = 1, ncol = 1, byrow = FALSE, dimnames = NULL)
{
data <- as.vector(data)
if (missing(nrow))
nrow <- ceiling(length(data)/ncol)
else if (missing(ncol))
ncol <- ceiling(length(data)/nrow)
.Internal(matrix(data, nrow, ncol, byrow, dimnames))
}
< environment: namespace:base >


Desgraciadamente los comentarios de estas funciones en los paquetes instalados se han eliminado para ahorrar memoria. Si queremos acceder al código fuente original, podemos descargar directamente el código de los paquetes del CRAN y descomprimirlos. Como los archivos están escritos con texto plano, podemos utilizar cualquier editor de texto para leer su contenido. Para los paquetes básicos de R, la carpeta es $R_HOME/src/library/NombrePaquete/R/.

En algunos casos, una función llama a otra función parcialmente oculta. Esa función oculta está en un namespace y se puede acceder a ella mediante la función getAnywhere("NombreFunción"). Esta función proporciona el namespace del que procede la función requerida y entonces podemos acceder a las fuentes de paquete correspondiente. Esto es así para los métodos S3 tales como plot.factor:

> plot.factor
Error: object 'plot.factor' not found
> getAnywhere("plot.factor")
A single object matching ‘plot.factor’ was found
It was found in the following places
registered S3 method for plot from namespace graphics
namespace:graphics
with value

function (x, y, legend.text = NULL, ...)
{
...
}
< environment: namespace:graphics >


El archivo que contiene el código de plot.factor es $R_HOME/src/graphics/R/plot.R.

Como ejemplo consideremos la función t.test(). Si queremos acceder a su código tenemos:

> t.test
function (x, ...)
UseMethod("t.test")
< environment: namespace:stats >

La llamada a UseMethod() oculta el código, pero nos indica que se trata de una función S3 genérica que llama a un método específico que será el apropiado para la clase del objeto. De modo que podemos preguntar cuales son los métodos apropiados así:

> methods(t.test)
[1] t.test.default* t.test.formula*

Non-visible functions are asterisked

De modo que ahora ya conocemos los métodos y podemos ver el código:

> getAnywhere(t.test.default)
A single object matching ‘t.test.default’ was found
It was found in the following places
registered S3 method for t.test from namespace stats
namespace:stats
with value

function (x, y = NULL, alternative = c("two.sided", "less", "greater"),
mu = 0, paired = FALSE, var.equal = FALSE, conf.level = 0.95,
...)
{
AQUÍ VEMOS EL CÓDIGO
}
< environment: namespace:stats >


En el caso de fuentes S4 es conveniente trabajar con los archivos fuente del paquete. Las funciones getClass(), getGeneric() y getMethod() nos pueden ayudar.

Cuando en el código de una función de R aparecen funciones del tipo .C(), .Call(), .Fortran(), .External() o .Internal() y .Primitive() es que llaman a código compilado. En estos casos deberemos mirar las fuentes originales (en C, C++ o Fortran) si queremos comprender todo el código.

Para més información consultar el artículo de Uwe Ligges:

R Help Desk en la Rnews_2006-4

El gráfico beanplot

Ciertamente, cuando se desea representar gráficamente unos datos cuantitativos univariantes existen varias posibilidades con diversas propiedades. Podemos hacer un histograma, un diagrama de tallo y hojas (stem-and-leaf), un diagrama de caja (boxplot), una densidad estimada y muchos más. Pero cuando se trata de comparar los valores de una variable cuantitativa para varias poblaciones o tratamientos, la mayoría no sirve. Por ejemplo, comparar un conjunto de histogramas o de gráficos de tallo y hojas es muy difícil, básicamente por el espacio que requieren. Así, para comparar datos univariantes en diferentes poblaciones casi siempre se utiliza el boxplot de Tukey o un gráfico de puntos unidimensional tipo stripchart cuando hay pocos datos.
Aunque hay algunas variaciones del boxplot clásico, como incrementar la anchura de la caja en función del número de observaciones con el parámetro varwidth=T, este gráfico se basa esencialmente en los estadísticos de posición como la mediana y las bisagras.
Otra posibilidad es el gráfico de violín descrito por Hintze and Nelson (1998) en el que una densidad se combina con los cuartiles del boxplot. En este gráfico no se muestran los valores atípicos (outliers).
Algunas críticas que se pueden hacer al boxplot se deben precisamente a la utilización de los cuartiles como elementos de difícil explicación (?) para los no matemáticos y la definición arbitraria del concepto de outlier. En un gráfico de violín la distribución subyacente es más visible, pero los datos concretos no se ven y no se conoce el número de observaciones del grupo.
El beanplot es una combinación de un gráfico de densidad (doble) con las marcas de todos los datos. Dichas marcas cambian de color si se salen del interior de la doble densidad y se alargan si coinciden algunos datos con el mismo valor.
Para poder comparar los grupos, se señalan las medias de cada grupo y la media general.
Por otra parte, si en la población general hay un factor con dos niveles, como el sexo, de puede considerar un beanplot asimétrico con dos densidades distintas en función del factor.
Sin embargo, aunque la vistosidad del beanplot parece aventajar al boxplot, no debemos olvidar que el boxplot se creó con la intención de tener un gráfico sencillo y sobre todo robusto, es decir, basado en estadísticos de orden.

Como es gratis, con cualquier conjunto de datos podemos hacer los dos gráficos, incluso no es difícil añadir al boxplot los datos concretos en otro color, y representaremos el que mejor explique nuestro experimento o los dos, si lo creemos conveniente.

La velocidad de la luz (3)

En primer lugar vamos a calcular los estadísticos más comunes con la ayuda del paquete fBasics y su función basicStats:


> library(fBasics)
Loading required package: MASS
Loading required package: timeDate
Loading required package: timeSeries

Attaching package: 'fBasics'

The following object(s) are masked from 'package:base':

norm

> basicStats(tiempo)
tiempo
nobs 66.000000
NAs 0.000000
Minimum -44.000000
Maximum 40.000000
1. Quartile 24.000000
3. Quartile 30.750000
Mean 26.212121
Median 27.000000
Sum 1730.000000
SE Mean 1.322658
LCL Mean 23.570591
UCL Mean 28.853652
Variance 115.462005
Stdev 10.745325
Skewness -4.391574
Kurtosis 25.518829

A primera vista podemos ver que hay mucha diferencia entre la media y la mediana, el recorrido es muy ancho y la asimetría elevada.

Podemos calcular los límites para que un valor sea considerado como atípico (outlier) con las bisagras (según Tukey)

LowerHinge - 1.5*(UH-LH), UpperHinge + 1.5*(UH-LH)


> lh <- fivenum(tiempo)[2]

> uh <- fivenum(tiempo)[4]

> iqr <- uh - lh

> c(lh,uh) + c(-1,1)*1.5*iqr

[1] 13.5 41.5


o con los cuartiles:

Q1 - 1.5*IQR, Q3 + 1.5*IQR


> q1 <- as.numeric(quantile(tiempo,0.25))

> q3 <- as.numeric(quantile(tiempo,0.75))

> c(q1,q3) + c(-1,1)*1.5*IQR(tiempo)

[1] 13.875 40.875


En cualquier caso resulta evidente que los valores –44 y –2 son atípicos. El tratamiento de los outliers es otro de los asuntos que requieren un poco de sentido común. Algunas veces los valores atípicos tienen un especial interés como evidencia de un suceso extraordinario. Un valor atípico en la distribución del brillo observado por un satélite de vigilancia puede representar el lanzamiento de un misil. Un valor atípico en la distribución de las alturas puede mostrar a un jugador de baloncesto. En estos casos la distribución general muestra la rutina o la normalidad, mientras que los sucesos extraordinarios caen fuera. Pero Newcomb esperaba una distribución bien formada, con un claro centro y, en cambio, dos valores atípicos molestaban.

Cuando los valores atípicos son sorprendentes o inesperados, en primer lugar, hemos de buscar la causa o explicación, tal como un error del equipo de medida o un error de escritura del dato. Casi todos los grandes conjuntos de datos contienen errores, con frecuencia por erratas al entrar los datos en el archivo informático. Los valores atípicos sirven para detectar estos errores y corregirlos al repasar los datos originales. Si el equipo de medida falla o alguna otra condición anormal ha provocado el valor atípico, entonces debemos borrar el dato sin ningún problema. También el valor atípico puede evidenciar una extraordinaria incidencia o una inesperada variabilidad de los datos.

Newcomb, finalmente, despreció el más pequeño de los valores (–44) y retuvo los otros. El basó su estimación de la velocidad de la luz en el término medio (la media) de sus observaciones. La media de las 66 observaciones es 26.21; la media de las 65 retenidas es 27.29. El fuerte efecto de un único valor –44 sobre la media es un motivo para descartarlo, cuando nuestro interés es el centro de la distribución como un todo.

Para estos datos, la media recortada al 5% es 27.4.

> mean(tiempo)
[1] 26.21212

> mean(tiempo[-6])
[1] 27.29231

> mean(tiempo, trim=0.05)
[1] 27.4


La mediana es directamente 27 o puede subir a 27,5 si suprimimos los dos valores negativos.


> median(tiempo)
[1] 27

> median(tiempo[-c(6,10)])
[1] 27.5


Todavía podemos contemplar los datos de Newcomb con otro gráfico. Cuando los datos representan observaciones similares tomadas a lo largo del tiempo, es bastante sensato dibujarlas con el tiempo o el orden temporal en que las observaciones fueron hechas.



La figura dibuja los lapsos de tiempo de Newcomb con su orden de recogida. El gráfico contiene una cierta sugerencia en el sentido de que la variabilidad (la anchura vertical del gráfico) disminuye con el tiempo. En particular, las dos observaciones atípicas fueron hechas muy pronto. Puede que él ganara experiencia y se hizo más hábil en la utilización de su aparato. Los efectos del aprendizaje, como en este caso, son bastante comunes. Si dejamos a Newcomb 20 observaciones para aprender, la media de los 46 restantes es 28.15.


> mean(tiempo[-(1:20)])
[1] 28.15217


Las medidas más modernas sugieren que el "verdadero valor" del lapso de tiempo medido por el experimento de Newcomb es 33.02. Eliminar los valores atípicos o permitir un período de aprendizaje no mueven el término medio más cerca del valor cierto. En todo caso, los ajustes basados sólo en criterios subjetivos son sospechosos. Si es posible, siempre es necesario hallar la razón de un valor atípico.