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)
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[-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.