El blog de los erreros

Una animación con R

2012-02-12T19:15:00.008+01:00

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

2012-02-08T19:13:00.008+01:00

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

2011-12-31T07:18:00.004+01:00

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

2011-11-05T19:16:00.002+01:00

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

2011-04-12T09:46:00.007+02:00

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

2011-03-10T11:45:00.000+01:00

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)

2011-02-12T00:20:00.004+01:00

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.

La velocidad de la luz (2)

2010-11-19T20:05:00.001+01:00

Figura 1. Diferentes histogramas con los 66 datos de Newcomb.

La primera observación evidente sobre los datos de Newcomb es que éstos varían. ¿Porqué? Se supone que, sobre todo, cada observación se hizo con el mismo patrón, medida por un observador cualificado y con el mismo aparato cada vez. Newcomb sabía, y nosotros también, que las medidas más precisas casi siempre varían. El entorno de cada medida es ligeramente diferente. El aparato cambia un poco con la temperatura, la densidad de la atmósfera también cambia de un día a otro, y así con otros detalles. Newcomb hizo lo que pudo para eliminar las fuentes de variación que podía anticipar y su experimento 0btuvo unas observaciones con menos variación que las de los científicos anteriores. Pero incluso el mejor experimento produce resultados variables. Por ello Newcomb tomó muchas medidas y no sólo una. La media de los 66 valores es menos variable que el resultado de una única medida, la media no depende tanto de la temperatura y la densidad atmosférica como una sola observación.

Si nos ponemos ahora en el lugar de Newcomb, estaremos tentados por calcular la media de los datos, convertir esta velocidad en una nueva y mejor estimación de la velocidad de la luz y salir corriendo para incrementar nuestra reputación al publicar el resultado. La tentación de hacer un cálculo rutinario y anunciar la respuesta se presenta siempre que los datos han sido recogidos para contestar una cuestión específica. Las computadoras son muy buenas para los cálculos rutinarios, de manera que sucumbir a la tentación es un camino fácil. El primer paso hacia la sofisticación estadística es resistir la tentación de calcular sin pensar. Como la variabilidad está presente en cualquier conjunto de datos, debemos primero examinar la naturaleza de la variación y, en algunos casos, quedaremos sorprendidos por lo que veremos.

La distribución

El patrón de variación de una variable es su distribución. La distribución presenta los valores numéricos que toma una variable y la frecuencia con la que cada valor ocurre.

La mejor manera de observar una distribución es representarla gráficamente. Conocemos varias técnicas gráficas para representar los datos de Newcomb. En particular el histograma.

Para hacer un histograma, se procede en tres pasos:

Se divide el rango de los datos en clases, casi siempre de la misma longitud.

Se cuenta el número de observaciones que caen en cada clase, es decir, las frecuencias absolutas y se construye la tabla de frecuencias.

Se dibuja el histograma.

Como el propósito de un histograma es mostrar la forma de una distribución, debemos estar muy atentos a los aspectos más visuales del gráfico. Nuestros ojos responden al área de las barras, de forma que ésta ha de ser proporcional a la frecuencia de la clase.

Ahora bien, antes de dibujar un histograma hay que decidir el número de clases, y eso es cuestión de opiniones.

La figura 1 se ha creado con el siguiente código:

tiempo <- Newcomb$Time*10^9 - 24800

png("histogramas.png") 
oldpar <- par(mfrow=c(2,2))    

hist.a <- hist(tiempo, breaks=-60+5*1:25, main=NULL, xlab="Tiempo (a)", ylab="Frecuencia", lim=c(0,30), col="grey") 
text(hist.a$mids, hist.a$counts, labels=ifelse(hist.a$counts==0,"",hist.a$counts), pos=3)     

hist.b <- hist(tiempo, breaks=-52.5+5*1:25, main=NULL, xlab="Tiempo (b)", ylab="Frecuencia", ylim=c(0,30), col="grey") 
text(hist.b$mids, hist.b$counts, labels=ifelse(hist.b$counts==0,"",hist.b$counts), pos=3)    

hist.c <- hist(tiempo, breaks=c(-44,-33.5,-23,-12.5,-2,8.5,19,29.5,40,50.5), main=NULL, xlab="Tiempo (c)", ylab="Frecuencia", ylim=c(0,45), col="grey") 
text(hist.c$mids, hist.c$counts, labels=ifelse(hist.c$counts==0,"",hist.c$counts), pos=3)    

hist.d <- hist(tiempo, breaks=-60+10*1:12, main=NULL, xlab="Tiempo (d)", ylab="Frecuencia", ylim=c(0,45), col="grey") 
text(hist.d$mids, hist.d$counts, labels=ifelse(hist.d$counts==0,"",hist.d$counts), pos=3)    

par(oldpar) 
dev.off()

Los histogramas de la figura 1 son diferentes en función de les clases que formemos.

En el caso (a) el histograma está formado con los intervalos
]-60,-55], ]-55,-50],... de anchura 5 y marca de clase en los puntos –57.5,
-52.5,..., -2.5, 2.5,..., 52.5, 57.5.

En (b) el histograma está formado con los intervalos
]-52.5,-47.5], ]-47.5,-42.5],..., también con una longitud de 5 y con las marcas de clase en los valores -50, -45,..., -5, 0, 5,...,60.

Si nosotros no decidimos los límites de los intervalos, la macro de Excel “Análisis de datos:Histograma” proporciona el gráfico (c). La anchura de los intervalos es 10.5 y el límite inferior del primer intervalo es el valor mínimo –44. El programa de Estadística SPSS construye un histograma con las marcas de clase -40, -30,..., -10, 0, 10,..., 40 y, por tanto, intervalos de longitud 10.

Cuando consideramos intervalos de anchura 10 y centrados con marcas de clase -55, -45,..., -5, 5,..., 55, tenemos el histograma (d).

Los histogramas anteriores tienen la característica común de que todos los intervalos son de la misma longitud. Con estos datos podríamos considerar un primer intervalo, como por ejemplo ]-45,-15] que tiene una frecuencia absoluta de 2, y dibujar una barra con la altura de manera que el área sea proporcional a 2.


hist.e <- hist(tiempo, breaks=c(-45,15,20,25,30,35,40), main=NULL, xlab="Tiempo (e)", ylab="Frecuencia", ylim=c(0,0.09), col="grey")
text(hist.e$mids, hist.e$intensities, labels=ifelse(hist.e$counts==0,"",hist.e$counts), pos=3)
lines(density(tiempo))

Otro tipo de gráfico que podemos utilizar con estos datos es el gráfico de tallo y hojas (stem and leaf):


stem(tiempo)

El resultado es, más o menos, así:

La construcción de este gráfico no es tan arbitraria como el histograma y permite la reproducción de los datos originales.

Cualquier científico de laboratorio sabe que incluso las medidas más precisas varían, pero espera una distribución simétrica y con un único pico. Entonces el mejor estimador del verdadero valor de la cantidad medida es el centro de la distribución. Los histogramas y el gráfico de tallo y hojas muestran que Newcomb no fue afortunado, ya que los datos contienen valores atípicos (outliers). ¿Qué hay que hacer con ellos?

(continuará)

La velocidad de la luz (1)

2010-11-16T10:04:00.009+01:00

Simon Newcomb

Los viajes de la luz son rápidos, pero no instantáneos. La luz tarda alrededor de un segundo en llegar desde la luna y sobre 10 billones de años desde el objeto más distante que se ha observado en nuestro expansivo universo. Como la radio y el radar también viajan a la velocidad de la luz, un valor ajustado de esta velocidad es muy importante para la comunicación con los astronautas y los satélites en órbita. Un valor ajustado de la velocidad de la luz también es muy importante para los diseñadores de computadoras, ya que las señales eléctricas viajan a esta velocidad.
La primera medida razonablemente ajustada de la velocidad de la luz se hizo hace más de 100 años gracias a los experimentos de A. A. Michelson y Simon Newcomb. La siguiente tabla contiene las 66 medidas hechas por Newcomb entre julio y septiembre de 1882.

28 22 36 26 28 28
26 24 32 30 27 24
33 21 36 32 31 25
24 25 28 36 27 32
34 30 25 26 26 25
-44 23 21 30 33 29
27 29 28 22 26 27
16 31 29 36 32 28
40 19 37 23 32 29
-2 24 25 27 24 16
29 20 28 27 39 23
Tabla. Medidas de Newcomb del lapso de la luz

Un conjunto de datos como estos no tiene sentido sin la información sobre su contexto. Debemos contestar algunas preguntas iniciales sobre cualquier conjunto de datos. Primero, ¿qué variable se está midiendo? Newcomb midió cuanto tiempo tardaba la luz en ir y volver desde su laboratorio en el río Potomac hasta un espejo en la base del monumento a Washington, una distancia de aproximadamente 7400 metros. Tal como nosotros calculamos la velocidad de un coche como el tiempo necesario para recorrer un kilómetro, Newcomb calculó la velocidad de la luz a partir del tiempo de ese viaje.
Contestar la pregunta “¿Qué variable se está midiendo?” requiere una descripción del instrumento utilizado para hacer la medida. Entonces podremos juzgar si la variable medida es apropiada para nuestro propósito. Este juicio frecuentemente pide un conocimiento de experto en el campo particular de estudio. Por ejemplo, Newcomb inventó un nuevo y complicado aparato para medir el lapso de la luz. Nosotros como estadísticos aceptamos el juicio de los físicos en el sentido de que este instrumento es apropiado para la tarea encomendada y más preciso que los instrumentos anteriores.
El estudio de Newcomb de la velocidad de la luz mide una variable claramente definida y fácilmente comprensible. Preguntas sobre las medidas son mucho más difíciles de contestar en los campos de las Ciencias Sociales y Económicas que en las Físicas. Nosotros nos podemos poner de acuerdo fácilmente en el tipo de medida apropiado para calcular la altura de una persona, pero ¿como medimos la inteligencia?
Los usuarios de datos deberían ser conscientes de que considerar los números con su valor nominal, sin pensar en la variable medida y el proceso utilizado para medirla, puede producir serias malinterpretaciones.

Las dos preguntas que faltan sobre cualquier conjunto de datos deberían contestarse con mayor senzillez: ¿Cuales son las unidades de medida? y ¿Como están registrados los datos?

El paquete MASS contiene la base de datos newcomb con los datos de la tabla anterior, pero desordenados. Es mejor obtener los datos originales en la siguiente dirección:

http://people.reed.edu/~jones/141/Newcomb.html

La primera medida de Newcomb del lapso de la luz fué 0,000024828 segundos. De manera que su unidad de medida fueron los segundos. Pero los valores de la tabla inicial no se parecen a 0,000024828. Estos números son incómodos de escribir y difíciles de tratar aritméticamente. En consecuencia nosotros hemos movido la coma decimal nueve posiciones a la derecha, esto es 24828, y hemos registrado únicamente la desviación respecto a 24800. Así pues, 28 es el resumen de 0,000024828 y –2 significa 0,000024798. Este procedimiento se conoce como codificación de los datos. Se debe codificar cuando los datos originales contienen muchas cifras de las cuales sólo algunas varían de observación en observación. Los datos codificados son más fáciles de leer. Además, si utilizamos una calculadora o una computadora, reducir el número de cifras mejora los cálculos aritméticos y su precisión.
Ahora que hemos entendido lo que significan los datos de la tabla y su procedencia, podemos empezar a mirarlos más de cerca y estudiarlos más a fondo. En cualquier caso, y aunque parezca que los datos ya se entienden, conviene siempre contestar las tres preguntas preliminares.

(continuará...)

Bibliografía

D.S. Moore & G.P. McCabe, Introduction to the Practice of Statistics, W.H. Freeman & Company.

Instalar R y Bioconductor en Ubuntu

2010-11-01T17:10:00.006+01:00

Estos días he tenido que instalar Ubuntu en algunos ordenadores portátiles del departamento. Como es habitual los portátiles vienen con Windows 7 instalado, pero a la mayoría de compañeros y compañeras nos gusta trabajar con Linux y su distribución más popular es Ubuntu o Kubuntu (versión de Ubuntu con el escritorio KDE). La verdad es que ya quedan muy lejos los días en los que era muy difícil instalar una distribución de Linux. Ahora con Ubuntu es realmente muy sencillo.
Después de instalar el sistema operativo y sus aplicaciones por defecto, llega la hora de instalar las aplicaciones para trabajar en serio y, entre ellas, nuestro amado R. La cosa es tan sencilla como ir al Centro de software de Ubuntu del menú Aplicaciones y buscar en el apartado de Ciencia e ingeniería -> Matemáticas el famoso programa R Commander o, como alternativa RKWard (mejor en Kubuntu). Se instalan con un sólo click. Fabuloso.
Sin embargo, si deseamos tener siempre la última versión, es mejor añadir el repositorio oficial del R-project a nuestro conjunto de fuentes de programas. Para ello seguiremos las instrucciones del siguiente enlace

http://cran.es.r-project.org/bin/linux/ubuntu/

Si no hay ninguna dificultad, así tendremos la última versión, tanto para 32 como para 64 bits, y cuando salga la siguiente el propio sistema de Ubuntu nos avisará para actualizarnos.
El único detalle que puede sorprender al principio es la utilización de un APT seguro con las claves de Vincent Goulet. Basta con ejecutar las dos instrucciones que se muestran en una Terminal y ya está.

Una vez instalados todos los paquetes de R, ya que siempre se necesita alguno de ellos, queremos instalar Bioconductor.
Bioconductor es un conjunto de paquetes de R, más de 400, para el análisis de datos de genómica que ha tenido un brillante desarrollo en los últimos años. De hecho, el proyecto de Bioconductor es uno de los más activos y ha impulsado notablemente el desarrollo del propio R.
Instalar el Bioconductor básico es muy sencillo. En una sesión de R (mejor si la hemos iniciado como administrador, con "sudo R"), debemos introducir las instrucciones

> source("http://bioconductor.org/biocLite.R")
> biocLite()

Pero, tras esperar un buen rato que se bajen, descompriman y compilen un motón de paquetes, nuestra frustración será grande cuando veamos que algunos paquetes no se han instalado con un mensaje final del tipo

installation of package 'XML' had non-zero exit status

El problema es que la instalación por defecto de Ubuntu no permite compilar algunos paquetes de Bioconductor. La solución es instalar previamente todo el software necesario.
En el blog de James Reid está muy bien explicado. Los paquetes de Ubuntu necesarios se instalan así:

sudo apt-get install build-essential curl graphviz-dev  libcurl4-gnutls-dev libboost-dev libgd2-xpm-dev libglu1-mesa-dev  libgtk2.0-dev libmysqlclient-dev libxml2-dev mesa-common-dev  sun-java6-jdk unixodbc-dev

De paso, dejo estas instrucciones como recordatorio para nuestras propias instalaciones.

El paquete HistData

2010-06-08T18:16:00.005+02:00

Gracias a dos alumnos del postgrado de Bioestadística de la UOC he descubierto el paquete HistData que contiene, entre otros, los datos de Florence Nightingale y que utilicé en el artículo anterior.
El paquete HistData recoge algunos de los más famosos conjuntos de datos de la historia de la Estadística, como los datos de Sir Francis Galton que sirvieron para entrever la normal bivariante y los conceptos de correlación y regresión. Pero además, el paquete contiene la reproducción con R de gráficos famosos asociados a esos datos. De esta forma los docentes podemos hacer memoria histórica estadística, que siempre va bien.
Así, debo rectificar parcialmente mi afirmación en el sentido de que no hay una función de R que reproduzca el coxcomb o rosa de Nightingale. Si bien es estrictamente cierta, el código que acompaña los datos de Nightingale permite generar el gráfico que veis al principio de este artículo y que compara las frecuencias de soldados muertos por diversas causas, antes y después de aplicar las mejoras sanitarias.
Este código utiliza el lenguaje gráfico del paquete ggplot2 que es una auténtica maravilla. Aprender este lenguaje no es trivial, pero si se domina se pueden hacer gráficos que combinan los mejores aspectos de los gráficos base de R y los lattice. Ya tengo otro buen tema de estudio para este verano.

data(Nightingale)
# For some graphs, it is more convenient to reshape death rates to long format
# keep only Date and death rates
require(reshape)
Night<- Nightingale[,c(1,8:10)]
melted <- melt(Night, "Date")
names(melted) <- c("Date", "Cause", "Deaths")
melted$Cause <- sub("\\.rate", "", melted$Cause)
melted$Regime <- ordered( rep(c(rep('Before', 12), rep('After', 12)), 3), levels=c('Before','After'))
Night <- melted

require(ggplot2)
cxc <- ggplot(Night, aes(x = factor(Date), y=Deaths, fill = Cause)) +
# do it as a stacked bar chart first
geom_bar(width = 1, position="identity", color="black") +
# set scale so area ~ Deaths
scale_y_sqrt() +
facet_grid(. ~ Regime, scales="free", labeller=label_both)
# A coxcomb plot = bar chart + polar coordinates
cxc + coord_polar(start=3*pi/2) +
opts(title="Causes of Mortality in the Army in the East") +
xlab("")

La idea fundamental es que un coxcomb es un diagrama de barras con áreas proporcionales a los datos y en coordenadas polares. Justamente su definición. ¡Fantástico!

Los datos de Florence Nightingale

2010-06-03T10:03:00.015+02:00

Florence Nightingale (Florencia, 12 de mayo de 1820 - Londres, 13 de septiembre de 1910), fue una enfermera británica considerada una de les pioneras en la práctica de la enfermería moderna y creadora del primer modelo conceptual de enfermería. Pero además destacó desde muy joven en matemáticas, aplicando después la estadística a la epidemiología y explotando la estadística sanitaria. Fue la primera mujer admitida en la Royal Statistical Society británica, y miembro honorario de la American Statistical Association.

Sin embargo, esta mujer, estudiada y reconocida en enfermería, es poco conocida entre los estudiantes de estadística y de matemáticas. Con motivo del pasado Día Internacional de la Mujer (8 de marzo) compartí con Carmina Olivé un acto de CCOO en la Universitat Politècnica de Catalunya (conjunto con la UB) en la que ella explicó su importancia en la enfermería moderna y yo glosé sus dotes científico-matemáticas-estadísticas. Voy a explicar brevemente mis argumentos.

A través de su trabajo como enfermera en guerra de Crimea, Florence Nightingale fue pionera al establecer la importancia de la higiene en los hospitales, en particular los hospitales de campaña. Ella reunió muchos datos en relación al coste en vidas por la falta de limpieza y, con sus gráficos y razonamientos matemáticos, también fue una avanzada en estadística aplicada. Con todo ello consiguió convencer a las autoridades de la época y, en especial, a los entonces muy obtusos militares para mejorar las condiciones sanitarias de los hospitales militares y civiles.

Entre los resultados estadísticos de Nightingale el más famoso es este gráfico, conocido como coxcomb o rosa de Nightingale, en el que se comparan los datos de soldados muertos por diversas causas, antes y después de aplicar las mejoras sanitarias en los hospitales.

Una versión más moderna y dinámica de este gráfico se puede hallar en el siguiente enlace:

http://understandinguncertainty.org/node/213

y las matemáticas y los datos para generar dicho gráfico se pueden consultar aquí:

http://understandinguncertainty.org/node/214

Por desgracia y que yo sepa, R no dispone de este tipo de gráfico. Llegados a este punto y, aunque hay algunas técnicas univariantes descriptivas para ayudar a Florence Nightingale, vamos a utilizar R y un par de técnicas multivariantes para describir sus datos.

La enfermera Nightingale recogió el número de soldados muertos por diversas causas: Zymotic diseases (ZD), Wounds & injuries (WI) y All other causes (OC), con ellos calculó unos índices de mortalidad anual relativa por cada 1000. Son precisamente esos índices los que vamos a utilizar como variables para un análisis de componentes principales.

dades <- read.table(file="dades.csv",header=T,sep=";")
dades$treatment <- factor(c(rep(0,12),rep(1,12)))
levels(dades$treatment) <- c("No", "Yes")
dades$active <- dades$size - (dades$zymotic_diseases + dades$wounds.injuries + dades$all_other)
dades$zymotic_diseases_ar <- dades$zymotic_diseases*12000/dades$size
dades$wounds.injuries_ar <- dades$wounds.injuries*12000/dades$size
dades$all_other_ar <- dades$all_other*12000/dades$size
attach(dades)

library(FactoMineR)
dades.res <- dades[8:11]
rownames(dades.res) <- idx
res.pca <- PCA(dades.res, quali.sup = 1)
plot(res.pca,habillage=1)

Con este código se obtiene el siguiente gráfico que se explica con mucha facilidad:

En él observamos que, si bien el primer eje admite la típica explicación de tamaño, en este caso depende básicamente de la causa evitable ZD y de OC (menos importante). El otro eje se explica con la variable WI.

En el gráfico con los 24 meses observados sobre las dos primeras componentes principales, donde los 12 primeros son antes de aplicar sus nuevos métodos de cuidado en los hospitales militares, observamos que los meses con el tratamiento se sitúan claramente a la izquierda, indicando la mejora en las muertes por causas evitables.

Por otra parte, y como los datos iniciales son frecuencias, podemos realizar un análisis de correspondencias tomando los meses como perfiles.

tabla <- data.frame(zymotic_diseases, wounds.injuries, all_other, active)
library(ca)
my.ca <- ca(tabla)
plot(my.ca, map="rowprincipal", xlim=c(-0.2,0.8), ylim=c(-0.2,0.2))

El resultado es similar. La explicación de los ejes que proporcionan los vértices (causas de mortalidad) es la misma que con el PCA e indica que la mejora en el tratamiento de las enfermedades en los hospitales hace que dichos meses estén a la izquierda. El único detalle es que en el gráfico se ha tenido que hacer un zoom sobre los puntos, ya que todos ellos quedaban fuertemente agrupados junto al vértice de "activos" que es el más pesado (el triángulo rojo del gráfico).

Discriminador lineal de Fisher

2010-05-16T22:02:00.008+02:00

Supongamos que queremos discriminar los datos entre dos poblaciones normales con diferente vector de medias y la misma matriz de covarianzas. Vamos a ver como hacerlo en el caso bivariante ya que así podemos representar el problema gráficamente.
En vez de realizar el problema con datos simulados o reales, lo vamos a tratar con dos distribuciones teóricas concretas.
Así pues, consideremos la distribución de una población normal bivariante con media mu1=(2.5,4) y matriz de covarianzas Sigma=(2,1,1,2) y otra población con la misma matriz de covarianzas, pero de media mu2=(6,3). El gráfico de arriba representa las densidades bivariantes de estas dos poblaciones. Un segmento de puntos une los dos puntos medios de las poblaciones.

La idea de Fisher consistió en hallar una combinación lineal de las variables originales (X1,X2) de la forma w1X1+w2X2 y tal que discrimine "el máximo posible" las dos poblaciones. Él mismo definió el criterio de máxima discriminación como maximizar la razón entre la suma de cuadrados entre grupos y la suma de cuadrados dentro de los grupos sobre la combinación lineal (discriminant scores). La solución viene dada por el vector w (o es proporcional al vector w)

w = Sigma^{-1} * (mu2 - mu1)

Es decir, entre todas las combinaciones lineales de la forma w1X1+w2X2, la que mejor discrimina es justamente esa. En el ejemplo propuesto, el vector solución w se ha dibujado en el gráfico como una flecha desde el punto mu1. Observemos que w tiene en cuenta la diferencia entre las medias pero también la relación de dependencia entre las variables.

Para visualizar este resultado he dibujado las densidades de las dos poblaciones si consideramos el vector w=(1,0), es decir, que la combinación lineal considerada sea simplemente X1.

Esta no es la solución óptima.

La solución exacta es w=(8/3,-11/6)

En el gráfico inicial, la recta que discrimina las dos poblaciones tiene un vector director ortogonal al vector w. Esa recta es la recta de puntos que equidistan de los dos puntos medios según la distancia de Mahalanobis.

El código para hacer todos los cálculos y los gráficos es el siguiente:


mu1 <- c(2.5,4)
mu2 <- c(6,3)
Sigma <- matrix(c(2,1,1,2),ncol=2)

dnormbv <- function(x1,x2,mu,Sigma) {
           sigma1 <- sqrt(Sigma[1,1])
           sigma2 <- sqrt(Sigma[2,2])
           rho <- Sigma[1,2]/(sigma1*sigma2)
           cte <- 1/(2*pi*sigma1*sigma2*sqrt(1-rho^2))
           cte * exp((-1/2) * (1/(1-rho^2)) * (
           (x1-mu[1])^2/sigma1^2 - 2*rho*(x1-mu[1])*(x2-mu[2])/(sigma1*sigma2) + (x2-mu[2])^2/sigma2^2)
           )}
          
x1<-seq(0,6,0.1)
x2<-seq(1,7,0.1)
z <- matrix(rep(0,length(x1)*length(x2)),ncol=length(x2))
for (i in 1:length(x1)) {
for (j in 1:length(x2)) {
z[i,j] <- dnormbv(x1[i],x2[j],mu1,Sigma)
}
}
plot(outer(x1,x2),xlim=c(0,9),ylim=c(0,9),xlab="",ylab="",type="n")
contour(x1,x2,z,nlev=4,add=T)

x1<-seq(2,9,0.1)
x2<-seq(0,7,0.1)
z <- matrix(rep(0,length(x1)*length(x2)),ncol=length(x2))
for (i in 1:length(x1)) {
for (j in 1:length(x2)) {
z[i,j] <- dnormbv(x1[i],x2[j],mu2,Sigma)
}
}
contour(x1,x2,z,nlev=4,add=T)

segments(mu1[1],mu1[2],mu2[1],mu2[2],lty="dotted")

SigmaInv <- matrix(c(2/3,-1/3,-1/3,2/3),ncol=2)

w <- SigmaInv %*% (mu2-mu1)

arrows(mu1[1],mu1[2],mu1[1]+w[1],mu1[2]+w[2],lwd=2)

pm <- (mu1+mu2)/2
a <- w[1] * pm[1] / w[2] + pm[2]
b <- -w[1]/w[2]
abline(a,b,lty="dashed")


ee <- 3.5
x1 <- seq(mu1[1]-ee,mu1[1]+ee,by=0.1)
x2 <- dnorm(x1,mean=mu1[1],sd=sqrt(Sigma[1,1]))
plot.new()
plot.window(xlim=c(-1,9),ylim=c(0,0.5))
axis(1)
axis(2)
lines(x1,x2)
x1 <- seq(mu2[1]-ee,mu2[1]+ee,by=0.1)
x2 <- dnorm(x1,mean=mu2[1],sd=sqrt(Sigma[1,1]))
lines(x1,x2)



E1 <- t(w) %*% mu1
E2 <- t(w) %*% mu2
var.w <- as.numeric(t(w) %*% Sigma %*% w)

ee <- 7
x1 <- seq(E1-ee,E1+ee,by=0.1)
x2 <- dnorm(x1,mean=E1,sd=sqrt(var.w))
plot.new()
plot.window(xlim=c(-6,16),ylim=c(0,0.5))
axis(1)
axis(2)
lines(x1,x2)
x1 <- seq(E2-ee,E2+ee,by=0.1)
x2 <- dnorm(x1,mean=E2,sd=sqrt(var.w))
lines(x1,x2)

Con datos reales o simulados podemos utilizar la función lda() del paquete MASS.

Bibliografía

Lattin, J. et al, Analyzing Multivariate Data, Ed. Brooks/Cole, Belmont (2003).

Análisis canónico de poblaciones

2010-03-27T19:26:00.027+01:00

El análisis canónico de poblaciones es una técnica de análisis multivariante que tiene el objetivo de representar varios grupos de individuos de forma óptima mediante unos ejes canónicos ortogonales. Eso se consigue de manera que la dispersión entre esos grupos sea máxima con relación a la dispersión dentro de cada grupo. Además, en esa representación, la distancia euclídea entre dos individuos coincide con la distancia de Mahalanobis entre ellos en las variables originales.
Se trata pues de una técnica de representación de datos en dimensión reducida, normalmente los dos primeros ejes, que puede acompañar gráficamente un MANOVA de un factor (la población).

Vamos a poner un ejemplo muy conocido. Se trata de los datos sobre cráneos de varones egipcios de cinco épocas históricas que se pueden obtener en el siguiente enlace:

http://lib.stat.cmu.edu/DASL/Datafiles/EgyptianSkulls.html

También podemos bajarlos desde la página del libro de Everitt (2005) y así los podremos cargar directamente en R con la siguientes instrucciones:

skulls <- source("/(path)/chap5skulls.dat")$value

str(skulls)

attach(skulls)

Donde el path debe ser la dirección a la carpeta donde hemos dejado el archivo una vez descomprimido. Así tendremos la base de datos skulls con cinco variables. La primera variable es el factor EPOCH y las otras cuatro son las medidas biométricas del cráneo estudiadas.

En primer lugar podemos realizar un MANOVA para contrastar la diferencia de medias entre los niveles del factor o poblaciones. No entraremos aquí en la comprobación de las hipótesis de normalidad y de igualdad de las matrices de covarianzas.

skulls.manova <- manova(cbind(MB,BH,BL,NH) ~ EPOCH)

summary(skulls.manova, test="Wilks") # test="Pillai" or "Hotelling" or "Roy"

El test rechaza la igualdad de medias y, por lo tanto, justifica el análisis canónico de poblaciones.

El siguiente paso es obtener el paquete candisc para un Canonical discriminant analysis ya que los ejes canónicos también sirven ese tipo de análisis discriminante.

library(candisc)

La función candisc realiza el análisis canónico discriminante, pero necesita como objeto principal un modelo lineal:

skulls.mod <- lm(cbind(MB,BH,BL,NH) ~ EPOCH)

Anova(skulls.mod, test="Wilks") # Manova

skulls.can1 <- candisc(skulls.mod, term="EPOCH")

plot(skulls.can1, conf=0.90, type="n")

El gráfico que se obtiene es el que vemos al principio de este artículo. Como corresponde a un ejemplo de manual, los dos primeros ejes canónicos representan muy bien a los datos y las poblaciones se separan de forma cronológica. También se representan las variables en una especie de biplot que explica mejor las diferencias entre las poblaciones. Todo muy bonito, pero no es un auténtico gráfico del análisis canónico de poblaciones, donde los círculos deben ser regiones de confianza con los datos sin escalar. Son círculos y no elipses los que representan a las poblaciones, ya que los ejes son independientes.
Por suerte, entre los resultados del objeto skulls.can1 disponemos de los coeficientes sin escalar que proporcionan las variables canónicas. Podemos pues aprovecharlos para calcular los scores de todos los datos sin escalar. Observemos también la utilización de la función aggregate para calcular las medias de cada población por separado como expliqué en el artículo Aplicar una función según el grupo.

# raw scores

scores <- as.matrix(skulls[,-1]) %*% skulls.can1$coeffs.raw
plot(scores[,1],scores[,2],xlim=c(-1,6),ylim=c(20,25), xlab="1er. eje canónico",ylab="2o. eje canónico",pch=16)

medias<-aggregate(skulls[,-1],skulls["EPOCH"],mean)

scores.medias <- medias %*% skulls.can1$coeffs.raw
text(scores.medias[,1],scores.medias[,2],1:5,pch=15,col="red")

Ahora ya sólo queda añadir los círculos de confianza para un nivel de confianza especificado, por ejemplo el 90%. Para ello vamos a calcular los radios de cada círculo según el tamaño de la muestra de cada población y la función symbols. Observemos que el objeto r es un vector.

resumen <- table(EPOCH)

g <- length(resumen) # número de poblaciones

p <- dim(skulls[,-1])[2] # número de variables

n <- as.vector(resumen) # tamaño de las muestras en cada población

radios <- function(g,p,n,conf.level=0.95) {

N <- sum(n)

F <- qf(conf.level,p,N-g-p+1)

sqrt(F*(N-g)*p/((N-g-p+1)*n))

}

r <- radios(g,p,n,0.90)

plot.new()

plot.window(xlim=c(0.5,4),ylim=c(21,24.5))

axis(1)

axis(2)

box()

title(main="Análisis canónico", xlab="1er. eje", ylab="2o. eje")

text(scores.medias[,1],scores.medias[,2],labels=levels(EPOCH),col="red")

symbols(scores.medias[,1],scores.medias[,2],circles=r,inches=FALSE,add=T,lwd=2,fg="red")

Bibliografía

Cuadras, C.M. (2010). Nuevos métodos de Análisis Multivariante. CMC Editions, Barcelona.

Everitt, B. (2005). An R and S-Plus Companion to Multivariate Analysis, Springer, London.

Arthur Thomson & R. Randall-Maciver (1905). The Ancian Races of the Thebaid. Oxford University Press.

Aplicar una función según el grupo

2010-02-14T12:01:00.000+01:00

Para aplicar una función a unas columnas de datos, separadamente según las clases o grupos que queramos, en R disponemos de dos funciones: tapply() y aggregate().
La función tapply() se aplica sobre un único vector de datos, cuyos valores se agrupan en función de una o unas variables o factores. Con los datos mtcars, por ejemplo, para saber el máximo valor de mpg según el número de cilindros hacemos

> data(mtcars)
> attach(mtcars)
> tapply(mpg, cyl, max)
4 6 8
33.9 21.4 19.2

En este caso el resultado es un vector, aunque también puede ser una lista si el resultado de la función aplicada no es un escalar:

> tapply(mpg,cyl,range)
$'4'
[1] 21.4 33.9

$'6'
[1] 17.8 21.4

$'8'
[1] 10.4 19.2

Acceder a estos resultados significa utilizar las llamadas a elementos de una lista:

> rangos <- tapply(mpg, cyl, range)
> rangos[[1]]
[1] 21.4 33.9

> rangos$'4'
[1] 21.4 33.9

> rangos[['4']]
[1] 21.4 33.9

> rangos[["4"]]
[1] 21.4 33.9

Cuando se utiliza más de una variable de agrupación y el resultado de la función a aplicar no es un escalar, el valor que retorna tapply() es más difícil de gestionar.

> rangos2 <- tapply(mpg, mtcars[c("cyl","am")], range)
> rangos2
am
cyl 0 1
4 Numeric,2 Numeric,2
6 Numeric,2 Numeric,2
8 Numeric,2 Numeric,2

> rangos2["4","0"]
[[1]]
[1] 21.5 24.4

Al llamar la función tapply() sin el argumento función resulta un índice de agrupación, según los grupos determinados, que puede ser útil en algunos casos:

> idx <- tapply(mpg,mtcars[c("cyl","am")]) > idx
[1] 5 5 4 2 3 2 3 1 1 2 2 3 3 3 3 3 3 4 4 4 1 3 3 3 3 4 4 4 6 5 6 4

Por otra parte, si lo que se pretende es resumir una o más columnas de un data.frame o matriz mediante un estadístico o función escalar, entonces utilizaremos la función aggregate(), donde el segundo argumento debe ser una lista.

> aggregate(mtcars[c("mpg","hp","wt")], mtcars$cyl, mean)
Error in aggregate.data.frame(mtcars[c("mpg", "hp", "wt")], mtcars$cyl, :
'by' must be a list
Calls: aggregate -> aggregate.data.frame

> aggregate(mtcars[c("mpg","hp","wt")], mtcars["cyl"], mean)
cyl mpg hp wt
1 4 26.66364 82.63636 2.285727
2 6 19.74286 122.28571 3.117143
3 8 15.10000 209.21429 3.999214

> aggregate(mtcars[c("mpg","hp","wt")], mtcars[c("cyl","am")], mean)
cyl am mpg hp wt
1 4 0 22.90000 84.66667 2.935000
2 6 0 19.12500 115.25000 3.388750
3 8 0 15.05000 194.16667 4.104083
4 4 1 28.07500 81.87500 2.042250
5 6 1 20.56667 131.66667 2.755000
6 8 1 15.40000 299.50000 3.370000

Finalmente, cuando el problema es mas complejo y se trata de aplicar una función no escalar a más de un vector según una determinada agrupación, deberemos combinar dos procedimientos como split() para agrupar y sapply() o lapply() para aplicar.

> xx <- split(mtcars[c("mpg","hp","wt")], mtcars["cyl"])

> aggregate(xx$'4', list(cyl=rep(4, dim(xx$'4')[1])), min)

cyl mpg hp wt

1 4 21.4 52 1.513

> sapply(xx, cor)

4 6 8

[1,] 1.0000000 1.0000000 1.00000000

[2,] -0.5235034 -0.1270678 -0.28363567

[3,] -0.7131848 -0.6815498 -0.65035801

[4,] -0.5235034 -0.1270678 -0.28363567

[5,] 1.0000000 1.0000000 1.00000000

[6,] 0.1598761 -0.3062284 0.01761795

[7,] -0.7131848 -0.6815498 -0.65035801

[8,] 0.1598761 -0.3062284 0.01761795

[9,] 1.0000000 1.0000000 1.00000000

Otra opción es utilizar la función by() que generaliza la función tapply().

> by(mtcars[c("mpg","hp","wt")],mtcars["cyl"],cor)
cyl: 4
mpg hp wt
mpg 1.0000000 -0.5235034 -0.7131848
hp -0.5235034 1.0000000 0.1598761
wt -0.7131848 0.1598761 1.0000000
------------------------------------------------------------
cyl: 6
mpg hp wt
mpg 1.0000000 -0.1270678 -0.6815498
hp -0.1270678 1.0000000 -0.3062284
wt -0.6815498 -0.3062284 1.0000000
------------------------------------------------------------
cyl: 8
mpg hp wt
mpg 1.0000000 -0.28363567 -0.65035801
hp -0.2836357 1.00000000 0.01761795
wt -0.6503580 0.01761795 1.00000000

Aplicar una función a una matriz o array

2010-01-31T14:00:00.000+01:00

La función apply

Para aplicar una función a una matriz o array, R dispone de la función apply que se utiliza con tres parámetros: el objeto array (una matriz es un array de dos dimensiones), la dimensión sobre la que actuaremos y la función que se aplicará. Como con la función sapply, es posible añadir al final todos los argumentos que precise la función que aplicaremos. En el caso de las matrices, si el segundo argumento es un 1, significa que la función se aplicará a las filas, mientras que si es un 2, se aplicará a las columnas.
En muchos casos se puede utilizar una instrucción apply en vez de un bucle. Además, el resultado es un vector o una matriz con nombres (si los había) y dicen que a menudo es más eficiente que un bucle.
Si queremos aplicar un procedimiento, primero definiremos la función a aplicar. Por ejemplo, si queremos saber el número de datos en cada una de las variables (columnas) que forman una base de datos (matriz), primero definiremos la función

n.datos <- function(x) n <- sum(!is.na(x))

y la aplicaremos sobre la base de datos

apply(base.de.datos, 2, n.datos)

La función apply se puede aplicar sobre elementos que no son arrays, pero debemos saber que entonces trata de convertirlos mediante un as.matrix o un as.array, dependiendo de las dimensiones.
Supongamos ahora que queremos estandarizar cada una de las variables de una base de datos, pero utilizando la mediana como estadístico de centralidad y la MAD como medida de dispersión.
El truco consiste en utilizar la función scale pero con los parámetros adecuados:

datos.estand <- scale(datos, center=apply(datos, 2, median), scale=apply(datos, 2, mad))

Sin embargo, si lo que queremos es calcular la suma o la media, es mejor utilizar funciones como rowSums, colSums, rowMeans o colMeans que son más eficientes. Además tienen un argumento na.rm= que permite prescindir de los valores faltantes.

La función sweep

Otra situación se presenta cuando queremos procesar una matriz por filas o por columnas, pero cada fila o columna de forma distinta, dependiendo de los valores de un vector dado. En estos casos utilizaremos la función sweep que tiene básicamente 4 argumentos. Como en la función apply, los dos primeros son la matriz y la dimensión sobre la que aplicaremos la función. El tercero es el vector de valores para procesar cada fila o columna de forma distinta y finalmente, el cuarto argumento es la función a aplicar.
Como funciones a aplicar podemos utilizar los operadores binarios como la suma "+", la resta "-" (por defecto), el producto "*" y el cociente "/", siempre con comillas. Por ejemplo, para dividir cada columna de unos datos por su máximo haremos

maximos <- apply(datos, 2, max)

sweep(datos, 2, maximos, "/")

Ahora bien, para otro tipo de funciones deberemos asegurarnos que la función trabaje correctamente con sweep. No todas lo hacen.
Si la función a aplicar calcula correctamente lo que queremos para cada vector de la matriz, pero no lo hace con sweep, podemos pensar en la función mapply.
La función mapply es la versión multivariante de sapply.
mapply aplica la función (primer argumento) sobre los primeros elementos de cada argumento, los segundos elementos, los terceros y así sucesivamente. El resultado se trata de simplificar, normalmente es una lista o vector.

mapply(rep, 1:4, 4:1)

Si la función necesita más argumentos, se los podemos pasar así:

mapply(rep, 1:4, MoreArgs=list(x=5))

Con un data.frame datos podemos hacer

mapply("/", datos, maximos)

y obtendremos un resultado similar al de sweep.

Aplicar una función a un vector o una lista

2009-12-21T17:44:00.010+01:00

Ciertamente R está preparado para operar todos los elementos de un vector de forma casi "natural". La mayoría de funciones son vectoriales. No ocurre lo mismo con una lista. Aunque muchas funciones devuelven una lista, que permite mucha flexibilidad, no se puede aplicar una función cualquiera directamente a una lista. Para ello disponemos de las funciones lapply y sapply.
Las dos funciones admiten, básicamente, dos argumentos: el primero es el objeto vector o lista sobre el que se aplica la función expresada en el segundo argumento. La diferencia entre ambas es que mientras lapply devuelve una (l)ista, sapply siempre tratará de (s)implificar el resultado en un vector, si es posible.

> mi.lista <- list(a=1:10, b=letters[1:3], cc=c(TRUE,FALSE,TRUE,FALSE))

> length(mi.lista)
[1] 3

> lapply(mi.lista, length)
$a
[1] 10

$b
[1] 3

$cc
[1] 4

> sapply(mi.lista, length)
a b cc
10 3 4

Observamos que mientras la función length aplicada a la lista nos devuelve el número de elementos de la lista (3), aplicada mediante un lapply o sapply nos da la longitud de sus elementos.

Cuando el objeto del primer argumento no es un vector o una lista, se forzará automáticamente que sea una lista mediante la función as.list. De este modo podemos aplicar estas funciones a un data.frame que no es estrictamente una lista, pero su conversión es sencilla.

> mi.df <- data.frame(a=1:10, b=letters[1:10], ca=runif(10), cb=rnorm(10))
> class(mi.df)
[1] "data.frame"

> sapply(mi.df,class)
a b ca cb
"integer" "factor" "numeric" "numeric"

> mi.df.num <- mi.df[ ,sapply(mi.df,class)=="numeric"]

> sapply(mi.df.num, mean)
ca cb
0.59009349 0.02423536

Cuando la función a aplicar (segundo argumento) necesita fijar sus propios argumentos, éstos se pueden incluir en las funciones lapply o sapply.

> sapply(mi.df.num, mean, trim = 0.05)

Por último, como estas funciones sirven para repetir el mismo cálculo sobre los elementos de un vector o lista, podemos pensar en utilizarlas siempre que podamos en lugar de un bucle (loop).
Además, tenemos una función replicate que permite simulaciones como la siguiente:

hist(replicate(100,mean(rexp(10))))

Fórmulas matemáticas en un gráfico

2009-12-11T20:07:00.002+01:00

Algunas de las funciones que permiten introducir texto en un gráfico son text, mtext, axis, title,...
Si en el argumento que explicita el texto de alguna de estas funciones se escribe una expresión, R la interpretará como una expresión matemática y le dará formato al estilo TeX. Dichas expresiones son el resultado de una función expression().
En una expresión de R ciertos nombres se interpretan como símbolos matemáticos, por ejemplo, alpha será la letra griega, sum será el símbolo de sumatorio, etc.
Para ver las diversas posibilidades, se puede consultar la ayuda help(plotmath) o la demostración demo(plotmath).

El gráfico inicial se consigue con el siguiente código:


curve(dnorm(x),-3,3,axes=F)
box()
axis(2)
axis(1,at=0,labels=c(expression(mu)))
title("Densidad normal")
text(2, 0.35, expression(paste(f(x)==frac(1, sigma*sqrt(2*pi)), " ",
                          plain(e)^{frac(-(x-mu)^2, 2*sigma^2)})),
   cex = 1.2)

En ciertas situaciones, como por ejemplo para definir una función, se necesita combinar texto con valores y variables. En ese caso no es posible utilizar la función expression() ya que las variables se tratarían literalmente como texto. Para solucionar este caso se utiliza la función substitute().
Ejemplo:


mifunc <- function(media) {           
          text(2, 3, substitute(paste("El valor de ", bar(x), " es ", media)))
          }

Este tipo de fórmulas matemáticas se puede reproducir en cualquiera de los dispositivos gráficos de pantalla como X11, Windows y Quartz y, gracias a la información de las fuentes Adabe Tipo 1 estándar que tiene R, también en PostScript y PDF.

Clases S3 y S4

2009-12-08T22:17:00.003+01:00

El sistema de clases de los objetos en R proporciona alguno de los mecanismos de la programación orientada a objetos como el despacho del método (method dispatch) y la herencia.
El método es la implementación de un algoritmo que representa una operación o función que un objeto realiza. El conjunto de los métodos de un objeto determinan el comportamiento del objeto. En R, el despacho del método consiste en examinar la clase de los argumentos de una función para decidir (despachar) la versión adecuada para los objetos de esa clase. No todas las funciones de R tienen despacho del método. Las que sí lo tienen se llaman funciones genéricas.
La herencia permite a los programadores crear nuevas clases, similares a otras ya existentes. Únicamente deberán proporcionar métodos adecuados para las nuevas clases o mantener los heredados. Un objeto de R que hereda las propiedades de un objeto ya definido, tiene como atributo de clase un vector que contiene la clase de ese objeto (en primer lugar), junto con las clases del objeto del que hereda.

El primer mecanismo de la programación orientada a objetos en R es el conjunto de clases S3 o del viejo estilo (old-style), donde el despacho del método se produce a través de las funciones genéricas del siguiente modo:

Supongamos que tenemos un objeto con el nombre de cdr de la clase Cuadrado, el cual es una subclase de Rectangulo que a su vez es una subclase de Forma. El mecanismo de despacho del método consiste en S3 y la función UseMethod. Utilizando S3 definimos un método Area para la clase Rectangulo como

  Area.Rectangulo <- function(objeto) { attr(objeto, "ladoA") * attr(objeto, "ladoB"); }

dado que un objeto Rectangulo tiene dos atributos ladoA y ladoB (ahora se puede acceder directamente a los argumentos con el operador @, es decir, objeto@ladoA y objeto@ladoB, respectivamente). Entonces definimos la función genérica Area así:

  Area <- function(objeto, ...) UseMethod("Area");

Cuando esta función se aplica en el objeto con Area(cdr), UseMethod despachará el método basado en la clase del primer argumento, es decir, objeto. Como cdr es de la clase Cuadrado, primero buscará un método llamado Area.Cuadrado. Si no existe, como en este caso, probará con Area.Rectangulo y así sucesivamente. Si no existiera ningún método específico para cualquiera de las clases del objeto, entonces se recurrirá al método Area.default que siempre debemos tener.

Observemos que las funciones genéricas S3 se pueden reconocer por la utilización de UseMethod en su código. Esto es importante, ya que las páginas de ayuda para una combinación de método y objeto depende de su nombre completo del tipo "function.class". Por ejemplo, la página de ayuda de la función summary no explica absolutamente nada sobre su actuación cuando le pasamos un objeto de la clase factor. Será mejor buscar la ayuda de la función summary.factor, aunque para que la función actúe sólo hay que escribir summary(objeto).

Cuando se crea una clase, asociadas a ella, deberemos crear un conjunto de funciones para extraer datos o información sobre los objetos de esa clase. Dada la convención explicada sobre los nombres de las funciones en las clases S3, podemos utilizar la función apropos para hallar los métodos disponibles para una clase:

> apropos('.*\\.factor$')
 [1] "all.equal.factor"     "as.character.factor"  "as.data.frame.factor"
 [4] "as.Date.factor"       "as.factor"            "as.list.factor"  
 [7] "as.POSIXlt.factor"    "as.vector.factor"     "codes.factor"    
[10] "[<-.factor"           "[.factor"             "[[.factor"             
[13] "format.factor"        "is.factor"            "is.na<-.factor"        
[16] "length<-.factor"      "levels<-.factor"      "Math.factor"           
[19] "Ops.factor"           "print.factor"         "rep.factor"            
[22] "summary.factor"       "Summary.factor"       "xtfrm.factor"

Así descubriremos que existen algunas funciones específicas para la clase factor.

Ahora bien, como el despachado del método que proporcionan las clases S3 está limitado al primer argumento de la función, y como las convenciones sobre los nombres que hemos explicado pueden crear alguna confusión, se decidió crear un nuevo sistema de clases S4 o de nuevo estilo (new-style). Éste es ahora el sistema preferido en el desarrollo de R.
En las clases S4, las funciones genéricas se reconocen por la llamada a una función standardGeneric en su definición.
Las funciones necesarias para trabajar con las clases S4 se hallan en el paquete methods. Por ejemplo, para saber si un objeto utiliza el nuevo estilo podemos hacer:


> temp <- c(20,15,15,20,20,30,25,25,30,15,25,30,30)  
> temp <- factor(temp)  
> isS4(temp)
[1] FALSE

Para saber los métodos asociados a una clase S4 como los objetos mle, resultado de la función de estimación de máxima verosimilitud, hacemos:


> library(methods)
> showMethods(class='mle')

Descubriremos que se puede calcular la matriz de varianzas-covarianzas de un objeto mle con la función vcov.

Aunque no hay una función genérica print para las clases S4, la función show permite ver el contenido de un objeto.
Por otra parte, los elementos que componen un objeto S4 se guardan en los llamados slots. Para ver los tipos de slots en un objeto, podemos utilizar la función showClass.


> library(stats4)
> showClass("mle")
Class “mle” [package "stats4"]

Slots:
                                                                        
Name:       call      coef  fullcoef      vcov       min   details minuslogl
Class:  language   numeric   numeric    matrix   numeric      list  function
            
Name:     method
Class: character

Para acceder directamente a un slot de un objeto, utilizaremos el operador @ del mismo modo que el operador $ para una lista. La función slot también obtiene el mismo resultado.


> slot(objeto.lme,"vcov")

Primeras jornadas de R

2009-11-30T12:29:00.002+01:00

Cuando me disponía a redactar un resumen de las Primeras jornadas de R celebradas en la Universidad de Murcia los pasados días 26 y 27 de noviembre, he descubierto con agrado que se me han adelantado.
En el blog Análisis y decisión podéis leer un amplio resumen de Carlos Gil Bellosta.
En fin, un ejemplo de trabajo colaborativo, él lo escribe y yo lo cito. ¡Ejem! Ya sé que no es eso.
Simplemente debo añadir que corroboro todas las buenas impresiones que nos llevamos las personas que asistimos del Departamento de Estadística de la Universidad de Barcelona: Miquel Calvo, Esteban Vegas y yo mismo. Además de conocernos y hablar de R y del software libre en general, también han quedado algunos compromisos de colaboraciones que poco a poco seguro que iremos concretando.
Un abrazo a todos y todas los que estuvisteis allí.

Factores numéricos

2009-11-28T19:57:00.005+01:00

Población africana con SIDA

Frecuentemente es necesario convertir una variable numérica en factor ya que algunas funciones de R así lo esperan. Sin embargo, si disponemos únicamente de un factor no podremos calcular algunos estadísticos u otras operaciones numéricas, aunque los niveles (o sus etiquetas) sean aparentemente numéricos.


> temp <- c(20,15,15,20,20,30,25,25,30,15,25,30,30)
> temp <- factor(temp)
> temp
 [1] 20 15 15 20 20 30 25 25 30 15 25 30 30
Levels: 15 20 25 30
> mean(temp)
[1] NA
Warning message:
In mean.default(temp) : argument is not numeric or logical: returning NA

Si no disponemos de los datos numéricos originales y deseamos convertir el factor a datos numéricos podemos probar así:


> temp.n <- as.numeric(temp)
> temp.n
 [1] 2 1 1 2 2 4 3 3 4 1 3 4 4

Pero el resultado son los valores enteros en los que se codifica internamente el factor.
Mejor si primero convertimos el factor a caracteres con las etiquetas de los niveles y luego esos mismos a valores numéricos:


> temp.n <- as.numeric(as.character(temp))
> temp.n
 [1] 20 15 15 20 20 30 25 25 30 15 25 30 30

Por otra parte, para crear un factor a partir de una variable continua se utiliza la función cut.
Los siguientes datos corresponden a la tasa de mortalidad del SIDA por cada mil habitantes en los países africanos en el año 2007:


> tasa <-  c(7.21,1.15,10.49,2.86,2.37,3.79,2.49,4.88,0.81,4.70,2.10,1.97,0.49,0.65,
            0.89,8.96,1.30,5.84,2.53,1.29,0.65,8.76,0.62,1.38,0.80,1.70,0.47,2.46) 

> tasa.f <- cut(tasa,breaks=seq(0,12,2))   

> table(tasa.f)
tasa.f
(0,2]   (2,4]   (4,6]   (6,8]  (8,10] (10,12]
   14       7       3       1       2       1

> class(tasa.f)
[1] "factor"

Valores perdidos

2009-11-22T12:20:00.001+01:00

En R, para representar un valor perdido (missing value) o que falta se utiliza el valor NA (not available) sin comillas. Se puede asignar este valor a una variable de la forma

> x <- c(3,2,NA,6)

pero no se puede consultar así

> x[3] == NA
[1] NA

para ello tenemos la función is.na():

> is.na(x[3])
[1] TRUE

Por otra parte, al realizar un determinado cálculo, como dividir por cero o la raíz de un negativo, el resultado puede ser Inf o NaN (not a number). La función is.nan() permite dirimir si se trata de uno de estos valores.

Cuando se importan datos de algún tipo de archivo, los valores perdidos pueden ser un problema si no se tratan de forma adecuada. Una posibilidad es utilizar el argumento na.strings= de la función read.table, al que le podemos pasar un vector con todos los valores del tipo carácter que R debe considerar perdidos y dar valor NA. En todo caso vale la pena revisar el resultado de la importación. Si los datos provienen de una hoja de cálculo, siempre podemos substituir los valores perdidos por NA, antes de la exportación del archivo al formato CSV.

Trabajar con valores perdidos

Algunas funciones de R disponen de argumentos para trabajar con datos que contienen valores perdidos. Por ejemplo, la mayoría de funciones estadísticas del tipo mean, var, sum, min, max, etc.
tienen un argumento na.rm= que se puede pasar a TRUE para que suprima los valores perdidos del cálculo. Para otras funciones que no tengan este parámetro podemos eliminar dichos valores creando un nuevo vector:

> x[!is.na(x)]
[1] 3 2 6

También algunas funciones de modelización estadística, como lm, glm, etc., disponen de un argumento na.action= al que le podemos pasar la acción a realizar para los valores perdidos. La más común es na.action=na.omit con lo que prescindirá de dichos valores en su procedimiento.
Otra opción es seleccionar del vector, matriz o data.frame(s) los individuos o filas que no tienen valores perdidos mediante la función complete.cases que devuelve un vector lógico:

> complete.cases(x)
[1] TRUE TRUE FALSE TRUE

Finalmente, señalar que en la conversión de un vector a factor, los valores perdidos no se consideran.

> y <- factor(x)
> y
[1] 3 2 6
Levels: 2 3 6

Si deseamos formar un nivel con los valores perdidos haremos

> y <- factor(x, exclude=NULL)
> y
[1] 3 2 NA 6
Levels: 2 3 6 NA

De R a ODF: el paquete odfWeave

2009-11-08T19:23:00.002+01:00

En el artículo LaTeX y R ya comenté las ventajas de combinar texto y fórmulas junto con instrucciones de R para obtener un informe dinámico que contendrá el texto, las fórmulas y los resultados numéricos y gráficos generados por R. Lamentablemente esto supone saber LaTeX (y R), y utilizar LaTeX no es sencillo. Sin embargo, la potencia de la función Sweave en R se ha extendido a otros formatos más amigables. Primero fué el paquete R2HTML para el formato html y ahora el paquete odfWeave para el formato ODF u OpenDocument.
El paquete odfWeave proporciona las funciones y el entorno para escribir informes automáticos en un formato libre y gratuito (como R), que se puede modificar a posteriori de forma muy sencilla y que permite su exportación a otros formatos como html o PDF.
La idea consiste en escribir un archivo ODF con un programa editor como OpenOffice con comentarios o explicaciones y código R para generar resultados y gráficos. La función odfWeave procesa ese archivo y nos entrega el documento ODF con el mismo texto y los resultados numéricos o gráficos insertados. Ese documento se puede modificar añadiendo otros elementos: vínculos, fotografías, encabezados y piés de página, etc. y finalmente imprimir o exportar.
Por ahora sólo funciona con documentos de texto, pero más adelante se puede ampliar a presentaciones y hojas de cálculo.

Veamos un ejemplo.

El documento de nombre ejemplo.odt contiene el siguiente texto:

Ejemplo

En este documento vamos a probar el paquete odfWeave.

Primero cargamos los datos de forma oculta para el documento final.

<<>>=

# Los datos son

coches <- mtcars[ ,1:6]

Los datos a estudiar contienen \Sexpr{dim(coches)[2]} variables medidas sobre \Sexpr{dim(coches)[1]} coches.

Podemos insertar una tabla resumen de los datos:

<<>>=

medias <- apply(coches,2,mean)

odfTable(medias, horizontal = TRUE)

También podemos añadir algún bonito gráfico:

<<>>=

pairs(coches)

Y así seguiría nuestro trabajo.

Si aplicamos el código:

library(odfWeave)
odfWeave(ejemplo.odt,resultado.odt)

obtenemos el documento resultado.odt

Ejemplo

En este documento vamos a probar el paquete odfWeave.

Primero cargamos los datos de forma oculta para el documento final.

Los datos a estudiar contienen 6 variables medidas sobre 32 coches.

Podemos insertar una tabla resumen de los datos:

mpg	cyl	disp	hp	drat	wt
20.091	6.188	230.722	146.688	3.597	3.217

También podemos añadir algún bonito gráfico:

Y así seguiría nuestro trabajo.

El formato de las imágenes y su tamaño se controlan con las funciones getImageDefs y setImageDefs.
Del mismo modo, el formato del texto o de las tablas se controla con las funciones de definición de estilo getStyleDefs y setStyleDefs, y las de asignación getStyle y setStyle para determinar algún elemento.

La clase odfTable se utiliza para convertir vectores, matrices y data frames a tablas ODF. La función odfCat sirve para escribir texto directamente en el formato ODF. También hay funciones para crear listas y para insertar imágenes externas.

Una vez tenemos el documento resultado, podemos utilizar OpenOffice para convertir manualmente este documento a otros formatos. Para ello se usa el menú Guardar como... o Exportar... del menú Archivo. En algunos sistemas operativos también disponemos de programas conversores desde la consola.

Para utilizar el paquete odfWeave es necesario que nuestro ordenador disponga de un compresor/descompresor de archivos ZIP. También deberemos tener en cuenta que los gráficos se graban en formato PNG, por si nuestro sistema no tiene el dispositivo adecuado.

Test de Levene para la igualdad de varianzas

2009-11-07T12:01:00.002+01:00

El test de Levene (1960) se usa para contrastar si k muestras tienen la misma varianza, es decir, la homogeneidad de varianzas. Otros contrastes, como por ejemplo el análisis de la varianza, suponen que las varianzas son iguales para todos los grupos. De ahí la importancia de verificar con el test de Levene esa hipótesis.

Este test es una alternativa al test de Bartlett. El test de Levene es menos sensible a la falta de normalidad que el de Bartlett. Sin embargo, si estamos seguros de que los datos provienen de una distribución normal, entonces el test de Bartlett es el mejor.

El test de Levene se resuelve con un ANOVA de los valores absolutos de las desviaciones de los valores muestrales respecto a un estadístico de centralidad (media, mediana o media truncada) para cada grupo. La elección del estadístico de centralidad de los grupos detemina la robustez y la potencia del test. Por robustez se entiende la habilidad del test para no detectar erróneamente varianzas distintas, cuando la distribución no es normal y las varianzas son realmente iguales. La potencia significa la habilidad del test para señalar varianzas distintas, cuando efectivamente lo son.

El artículo original de Levene proponía la media como estadístico de centralidad. Brown y Forsythe (1974) extendieron este test al utilizar la mediana e incluso la media truncada al 10%. Sus estudios de Monte Carlo mostraron que la utilización de la media truncada mejoraba el test cuando los datos seguían una distribución de Cauchy (colas grandes) y la mediana conseguía mejorarlo cuando los datos seguían una (distribución asimétrica). Con la media se consigue el mejor test para distribuciones simétricas y con colas moderadas.

Así pues, aunque la elección óptima depende de la distribución de los datos, la definición del test basada en la mediana es la recomendación general ya que proporciona una buena robustez para la mayoría de distribuciones no normales y, al mismo tiempo, una aceptable potencia. Si conocemos la distribución de los datos, podemos optar por alguna otra de las opciones.

El paquete car Companion to Applied Regression de J. Fox trae una función levene.test con la mediana como medida de centralidad de los grupos. El mismo test está incorporado al paquete Rcmdr.

Aplicar el test de Levene a unos datos con la opción de la media truncada como medida de centralidad es un ejercicio sencillo:

datos <- c(rcauchy(50,0,10),rcauchy(50,0,20),rcauchy(50,0,30))

grupo <- gl(3,50)

boxplot(datos~grupo)

func <- function(x) mean(x,trim=0.1)
medias_trunc <- tapply(datos,grupo,func)
desv <- abs(datos-medias_trunc[grupo])
summary(aov(desv~grupo))

Referencias

Brown, M. B. and Forsythe, A. B. (1974), Journal of the American Statistical Association, 69, 364-367.

Fox, J. (2002), An R and S-PLUS Companion to Applied Regression, Sage.

Levene, H. (1960). In Contributions to Probability and Statistics: Essays in Honor of Harold Hotelling, I. Olkin et al. eds., Stanford University Press, pp. 278-292.

Contrastes de normalidad

2009-10-31T17:21:00.007+01:00

Cuando se dispone de muy pocos datos no es posible utilizar el omnipresente test ji-cuadrado para contrastar la normalidad de la muestra. Por ello quedan descartados la función chisq.test y el histograma como gráfico más apropiado.

Algunos test no paramétricos sobre la distribución de la población se basan en la función de distribución empírica EDF de la muestra. La idea es comparar sus valores (en azul en el gráfico anterior) con la función de distribución teórica (en rojo).
En el caso de la distribución normal como distribución teórica disponemos de todo un paquete llamado nortest con varios test.

El más famoso es el test de Lilliefors que es una variante del test de Kolmogorov-Smirnov. Aunque el estadístico que se obtiene con lillie.test(x) es el mismo que el que se obtiene con ks.test(x, "pnorm", mean(x), sd(x)), no es correcto utilizar el p-valor de éste último con la hipótesis de normalidad (media y varianza desconocidas), ya que la distribución del estadístico es diferente cuando estimamos los parámetros. Dicho estadístico es el valor absoluto de la máxima diferencia entre los valores de la distribución empírica y la teórica.
Sin embargo, el test de Lilliefors ha quedado superado por el test de Anderson-Darling o el de Cramer-von Mises.
El test de Anderson-Darling ad.test es el test EDF recomendado por Stephens (1986). Comparado con el test de Cramer-von Mises cvm.test (como segunda elección) da mayor peso a las colas de la distribución.

Por otra parte, el test de Shapiro-Wilk se puede calcular con la función shapiro.test. Este test se basa en el estadístico W proporcional al cuadrado de una combinación lineal de los estadísticos de orden.

El estadístico del test de Shapiro-Francia sf.test es simplemente la correlación al cuadrado entre los valores muestrales ordenados y los cuantiles (aproximados) esperados para la distribución normal estandar. El p-valor se calcula con la fórmula dada por Royston (1993).

A pesar de lo dicho al principio, el paquete nortest dispone de la función pearson.test para resolver el test ji-cuadrado. En todo caso no se recomienda por su inferior potencia comparado con los test anteriores.

Recientemente, RKward ha añadido al menú Distributions el test de normalidad de Jarque y Bera que se obtiene con el paquete tseries y la función jarque.bera.test. El estadístico de este test se basa en los valores muestrales de asimetría y curtosis. Judge et al. (1988) y Gujarati (2003) recomiendan este test.

Finalmente, como se puede observar en el gráfico inicial de este artículo, resulta muy difícil para el ojo humano apreciar si la distribución empírica se ajusta a la teórica, al menos con ese tipo de gráfico.

y <- c(-0.1, -1.8, -0.1, -0.8, -1.0, 0.5, 1.4, -0.8, -0.2, -0.3, -0.4, 0.5)
Fn12 <- ecdf(y)
plot(Fn12, col.p="blue", col.h="blue", lwd=2, main="Empirical Cumulative Distribution Function")
abline(v=knots(Fn12),lty=2,col='gray70')
curve(pnorm(x), col="red", add=T)

Por ello es mejor utilizar el gráfico qqnorm(y) que dibuja los valores muestrales en el eje Y y los cuantiles teóricos en el eje X. Observemos que en este caso los cuantiles teóricos son

sort(qqnorm(y)$x)==qnorm((1:12-0.5)/12) La recta dibujada une los puntos del primer y tercer cuartil.

Referencias

Stephens, M.A. (1986): Tests based on EDF statistics. In: D'Agostino, R.B. and Stephens, M.A., eds.: Goodness-of-Fit Techniques. Marcel Dekker, New York.

Royston, P. (1993): A pocket-calculator algorithm for the Shapiro-Francia test for non-normality: an application to medicine. Statistics in Medicine, 12, 181–184.

Thode Jr., H.C. (2002): Testing for Normality. Marcel Dekker, New York.

R 2.10.0

2009-10-29T08:26:00.003+01:00

Ya está disponible la versión 2.10.0 para las principales distribuciones de Linux, para Windows y para Mac.
Los (K)ubunteros sólo debemos esperar que nuestro actualizador nos avise del cambio.
Esta nueva versión trae importantes cambios, especialmente en el sistema de ayuda, y mejoras en multitud de funciones, además de las correcciones de algunos fallos.
Más detalles en what's new.

Interpolación de colores

2009-10-25T19:58:00.004+01:00

En algunos gráficos es preciso disponer de una rampa de colores, es decir, ir de un color a otro(s) de forma gradual. Para ello, R dispone de las funciones colorRamp() y colorRampPalette() que podemos utilizar con dos, tres o más colores.
Estas funciones devuelven funciones que interpolan un conjunto de colores para crear una rampa de nuevos colores o una paleta. La función resultado de colorRamp() va del intervalo [0,1] a la rampa de colores. En cambio colorRampPalette() proporciona una función que toma como argumento un número entero y devuelve exactamente ese número de colores interpolados.

Por ejemplo,

> colorRampPalette(c("red", "orange", "blue"), space = "rgb")(10)
[1] "#FF0000" "#FF2400" "#FF4900" "#FF6E00" "#FF9200" "#E2921C" "#AA6E54"
[8] "#71498D" "#3824C6" "#0000FF"

nos da una paleta de 10 colores interpolados del rojo al azul, pasando por el naranja.

El gráfico de este artículo se consigue con la paleta de colores jet.colors al estilo de Matlab:

> jet.colors <-
colorRampPalette(c("#00007F", "blue", "#007FFF", "cyan", "#7FFF7F", "yellow", "#FF7F00", "red", "#7F0000"))
> filled.contour(volcano, color = jet.colors, asp = 1)

Georg Cantor

2009-10-25T08:29:00.008+01:00

He leído en arieder's Blog que KDE Edu está preparando una nueva aplicación para la versión 4.4: Cantor.

Se trata de utilizar algunas de las mejores aplicaciones matemáticas libres totalmente integradas en el escritorio KDE de linux (aunque también está en Windows).
Entre estas aplicaciones de cálculo numérico, simbólico y estadístico están, por ahora, Sage, Maxima y R.

El nombre de este interfaz se ha tomado del matemático padre de la Teoría de Conjuntos: Georg Cantor.

Aquí tenéis una imagen de la nueva aplicación:

Colores

2009-10-24T18:28:00.009+02:00

En R hay varias formas de especificar un color. La más sencilla es llamar al color por su nombre, eso sí, en inglés. Por ejemplo, el color azul con su intensidad máxima se consigue con "blue".
El listado completo de nombres que R reconoce se puede obtener con la función colors() o colours():

> colors()[1:10]
[1] "white" "aliceblue" "antiquewhite" "antiquewhite1"
[5] "antiquewhite2" "antiquewhite3" "antiquewhite4" "aquamarine"
[9] "aquamarine1" "aquamarine2"

También se puede especificar un color con el código de intensidades RGB. Para ello disponemos de la función rgb() en la que debemos introducir la intensidad de los tres colores básicos: R rojo, G verde y B azul.
La intensidad se mide por defecto de 0 a 1, pero también se puede fijar entre 0 y 255.
Por ejemplo, el azul corporativo de la Universidad de Barcelona es un Azul Pantone 285, que según la equivalencia del artículo Convert Pantone Colours to RGB se puede utilizar como

rgb(58, 117, 196, max=255)

Una alternativa consiste en especificar las intensidades en forma hexadecimal. El azul UB es

"#3A75C4"

Los dos dígitos hexadecimales para cada intensidad se mueven entre el cero 00 y el máximo FF.

También disponemos de una función hsv() para especificar los colores con la tripleta HSV Hue-Saturation-Value. Sin demasiada precisión se puede decir que el matiz (hue) indica la posición del color en el arco iris, desde el rojo (0) , a través del naranja, amarillo, verde, azul, índigo, hasta el violeta (1).
La saturación fija la intensidad del color entre la palidez y la viveza.

Por último, el valor o brillo se usa para describir que tan claro u oscuro parece un color, y se refiere a la cantidad de luz percibida.

La función rgb2hsv() nos puede ayudar a convertir un color de RGB a HSV. El azul UB es

> rgb2hsv(58,117,196)
[,1]
h 0.5954106
s 0.7040816
v 0.7686275

Otro sistema para especificar un color es fijar un conjunto de colores o paleta y seleccionar uno de ellos con un número entero, es decir, su posición en la paleta. Por defecto tenemos ya definida una paleta de colores que podemos ver con la función palette()

> palette()
[1] "black" "red" "green3" "blue" "cyan" "magenta" "yellow"
[8] "gray"

En esta paleta el color azul es el número 4.
Así pues podemos definir nuestra propia paleta y utilizarla con los números de posición de los colores.

Por último decir que R proporciona algunas funciones que definen grupos de colores que tienen algún sentido de conjunto. Ejemplos: rainbow(), heat.colors(), terrain.colors(), topo.colors(), cm.colors() y grey()o gray(). En el gráfico que encabeza este artículo podemos ver los colores de estos grupos.

Símbolos para los datos

2009-10-12T21:30:00.006+02:00

R proporciona un conjunto de 26 símbolos que podemos seleccionar con el parámetro pch para un gráfico. Este parámetro puede tomar como valor un número entero entre 0 y 25 o un único carácter de texto (ver figura). Algunos de los símbolos predefinidos (entre el 21 y el 25) se pueden rellenar con un color diferente al del borde. El color se concreta con el parámetro bg.
Si pch es un carácter entonces la letra se usa como símbolo. El carácter "." es un caso especial para el que el dispositivo trata de dibujar un punto lo más pequeño posible.
El tamaño de los símbolos está relacionado con el tamaño del texto y el parámetro cex.

El gráfico de arriba está en el libro R Graphics de Paul Murrell y se consigue con el siguiente código:

library(RGraphics)
ncol <- 6
nrow <- 5
grid.rect(gp=gpar(col="grey"))
for (i in 1:nrow) {
for (j in 1:ncol) {
x <- unit(j/(ncol+1), "npc")
y <- unit(i/(nrow + 1), "npc")
pch <- (i - 1)*ncol + j - 1
if (pch > 25)
pch <- c("A", "b", ".", "#")[pch - 25]
grid.points(x + unit(3, "mm"), y,
pch=pch, gp=gpar(fill="grey"))
grid.text(pch, x - unit(3, "mm"), y, gp=gpar(col="grey"))
}
}

Gráfico de densidad condicional

2009-10-12T11:51:00.004+02:00

Si se pretende estudiar la regresión cuando la variable respuesta es dicotómica, un primer paso para comprender el posible modelo es examinar algún tipo de gráfico como el de la densidad condicional.

Con los datos de plasma del paquete HSAUR, el gráfico de arriba se consigue con las instrucciones

> library(HSAUR)
> data(plasma, package=HSAUR)
> cdplot(ESR ~ fibrinogen, data = plasma)

En el gráfico se observa que los niveles altos de la proteína se asocian con valores de ESR por encima de 20 mm/hr.

A continuación podemos ajustar un modelo de regresión logística con la función glm. Si tomamos como distribución de la variable respuesta la binomial, el código puede ser

> plasma.glm <- glm(ESR ~ fibrinogen, data = plasma, family = binomial())
> coef(plasma.glm)

Por defecto, la función de enlace (link function) para la familia binomial es la función logística.

De R a HTML: el paquete hwriter

2009-10-10T20:03:00.005+02:00

Three iris flowers

Para completar un informe como resultado de un análisis estadístico o bioinformático es muy importante la combinación de datos e hipervínculos en un documento HTML. Existen ya varias herramientas para exportar datos de R a documentos HTML. Por ejemplo, el paquete R2HTML es capaz de traducir una buena colección de objetos de R a HTML, pero con él no resulta nada fácil combinarlos en una presentación estructurada. Por otra parte, el paquete xtable puede traducir matrices de R de forma sencilla, pero tampoco puede combinar varios elementos de HTML y tiene un conjunto muy limitado de instrucciones de formato.
Sin embargo, el paquete hwriter de Gregoire Pau y Wolfgang Huber permite traducir diversos objetos de R a HTML y combinarlos en un documento con un cierto diseño y formato.

Este paquete utiliza la función hwriter que permite escribir un objeto en un documento HTML, al mismo tiempo que le añadimos algún formato o un hipervínculo. Pero la clave es que si existe una conexión abierta hacia un documento HTML ya creado, entonces las instrucciones HTML se añaden a dicho documento.

El siguiente código escribe una sencilla página con un título y una tabla:

> library(hwriter)
> hwriter("Three iris flowers", "iris.html", heading=1)
> p <- openPage("iris.html")
> himg <- hwriteImage(c("iris1.jpg","iris2.jpg","iris3.jpg"),
+ link=c("http://en.wikipedia.org/wiki/Iris_virginica",
+ "http://en.wikipedia.org/wiki/Iris_versicolor",
+ "http://en.wikipedia.org/wiki/Iris_virginica"), table=FALSE)
> mat <- rbind(himg, c("Setosa","Versicolor","Virginica"))
> rownames(mat) <- c("Image", "Species")
> hwrite(mat, p, br=TRUE, center=TRUE,
+ row.bgcolor=list(Species=c("#ffaacc", "#ff88aa", "#ff6688")),
+ col.bgcolor="#ffffaa", row.style=list(Species="text-align:center"))
> closePage("iris.html")

Observemos que la función es capaz de escribir una mezcla de imágenes y datos en una tabla, junto con hipervínculos, en la conexión abierta p.

Podemos saber más de este paquete en el artículo de RJournal_2009-1.

También podemos ver en acción todas las funciones y sus parámetros en una página de ejemplo que se muestra con la instrucción

> example(hwriter)

I Conferencia Hispana R-Project

2009-10-10T13:30:00.008+02:00

Con el subtítulo de I Jornadas de usuarios de R en castellano, se celebrará los días 26 y 27 de noviembre en Murcia una reunión para impulsar el desarrollo y la documentación de R en castellano.

El objetivo principal es conectar a los usuarios de habla hispana y organizar algún tipo de comunidad que optimize los esfuerzos y canalize la colaboración.

Estas jornadas se plantean como una reunión de distintos profesionales que tienen en común:

El uso de R como instrumento científico o docente
La necesidad de conectarse con redes de trabajo en R
Interés por el desarrollo de materiales (textos o librerías)
Deseo de colaborar en proyectos ya iniciados
Participar en la traducción de materiales desarrollados
Crear grupos de trabajo de distintas disciplinas

Se pretende la constitución de grupos de trabajo y la puesta en común de materiales desarrollados en la lengua de Cervantes y dispersos en la red.
Podéis hallar toda la información y los detalles para la inscripción en la página web: http://www.ereros.org.

¡Allí nos veremos!

Integración de datos

2009-10-01T10:08:00.005+02:00

Ferran Reverter.
Cada vez hay más demanda de integración de datos provenientes de diferentes fuentes. Especialmente en estudios médicos donde la información de la genómica (datos de microarrys) se complementa con datos clínicos. En R disponemos, entre otros, de tres packages que implementan técnicas de análisis multivariante muy adecuados para la integración. El package FactoMineR desarrolla el Análisis Factorial Múltiple (Escofier-Pagès, 1998), el CCA implementa el Análsis de Correlaciones Canónicas Regularizado (Vinod, 1976) y el VEGAN desarrolla el Análisis de Correspondencias Canónico (Ter Braak, 1986).

Cádiz y R-Commander

2009-09-27T21:03:00.008+02:00

Entre los días 23 y 25 de septiembre se ha celebrado en Cádiz la XII Conferencia Española de Biometría a la que hemos asistido un buen grupo de investigadores e investigadoras, tanto españoles como portugueses, incluso alguno/a venido/a de iberoamérica. Entre las comunicaciones se ha notado una gran presencia de cálculos y gráficos realizados con R o alguno de sus paquetes. Pero lo que yo quiero destacar aquí es la realización de un curso de R con R-Commander paralelo a las sesiones que ofrecieron los organizadores de forma gratuita a todos los asistentes. Además nos regalaron un manual de Estadística Básica con R y R-Commander en español.
La Universidad de Cádiz ha hecho una apuesta decidida por el software libre que es realmente envidiable. Y como no, la utilización de R en todos los ámbitos de la estadística. A destacar la puesta en marcha, en el curso 2006/2007, del proyecto R-UCA.
R-Commander o su paquete Rcmdr es es una interfaz gráfica para R desarrollada por John Fox. La página web de R-Commander es http://socserv.socsci.mcmaster.ca/jfox/Misc/Rcmdr/. Su objetivo es cubrir toda la estadística de un curso elemental. Lo consigue de sobras.

Por cierto, os recomiendo el restaurante "Los naranjos" en Zahara de la Sierra. Fantástico.

RKWard 0.5.1

2009-09-13T20:29:00.002+02:00

Disponible desde el 4 de agosto.
RKWard es una interfaz transparente easy to use para el lenguaje R en un escritorio KDE para linux o windows.
Aunque todavía está en desarrollo, RKWard es un potente entorno de trabajo para el lenguaje R que facilita la utilización de las instrucciones, las funciones estadísticas, la creación de gráficos y su integración en ofimática. La intención es que RKWard sea un substituto libre para algunos paquetes comerciales de estadística.

Además de la facilidad de uso, hay tres aspectos que son especialmente importantes:

Es una interfaz transparente para el lenguaje R. Es decir, no esconde la sintaxis del lenguaje, sino que la proporciona de una forma más conveniente. Expertos y novatos pueden realizar la mayoría de los trabajos habituales. Una interfaz únicamente gráfica no podría servir a todo el poder del lenguaje R. En algunos casos, los usuarios querrán ajustar algunas funciones a sus necesidades particulares, por ejemplo para automatizar algunos trabajos. RKWard muestra la sintaxis, es visible para el usuario, y lo hace posible.
Para la salida, RKWard se esfuerza en separar el contenido y el diseño. No trata de diseñar sus propias tablas o gráficos, que han de convertirse "a mano" al estilo utilizado en una determinada publicación. Actualmente RKWard utiliza HTML para su producción. Sin embargo, utilizando las definiciones de estilo, es posible modificar el formato de salida para que coincida con la demanda de una publicación. En futuras versiones RKWard incluso buscará una mayor integración con alguna de las suites de ofimática existentes.
Se basa en un lenguaje, que no sólo es muy potente, sino que también es extensible y para el que ya existen cientos de extensiones.

Las novedades de la última versión son diversas y van desde la mejora de algunas funcionalidades, a la corrección de errores, pasando por el ajuste de algunos iconos.
Ya está disponible una compilación para Debian Sid i386 y para otras distribuciones podéis ir a la página de Binaries and Build Scripts.

Estimación de una función de densidad bivariante

2009-09-06T13:34:00.007+02:00

Del mismo modo que la función density se usa para calcular un conjunto de estimadores de densidad kernel, podemos estimar una densidad bivariante con la función bkde2D del paquete KernSmooth (Wand and Ripley, 2005).

El gráfico de este artículo se ha obtenido con los datos de energía y temperatura de unas estrellas (en logaritmos) del paquete HSAUR y el siguiente código:

> install.packages(c("HSAUR","KernSmooth"))
> library(HSAUR)
> library(KernSmooth)
> data(CYGOB1, package="HSAUR")
> CYGOB1d <- bkde2D(CYGOB1, bandwidth = sapply(CYGOB1, dpik))
> persp(x = CYGOB1d$x1, y = CYGOB1d$x2, z = CYGOB1d$fhat,
+ xlab = "log surface temperature",
+ ylab = "log light intensity",
+ zlab = "estimated density",
+ theta = -35, axes = TRUE, box = TRUE)

También podemos hacer un gráfico de contornos con la función contour
> contour(x = CYGOB1d$x1, y = CYGOB1d$x2, z = CYGOB1d$fhat,
+ xlab = "log surface temperature",
+ ylab = "log light intensity")

Font: B.S. Everitt & T. Hothorn, A Handbook of Statistical Analysis Using R, Chapman & Hall/CRC, 2006.

SPSS + R = PASW Statistics Developer

2009-09-05T10:48:00.003+02:00

Según SPSS Inc., su nuevo producto PASW Statistics Developer permite que los programadores de R y Python integren procedimientos en la sintaxis de PASW® Statistics para que los usuarios puedan acceder a ellos, dicen que de forma más "amigable".

Tres son los retos que pretenden resolver con su producto: la implementación de funciones y paquetes de R, su integración en un entorno de menús y procedimientos tipo SPSS y el manejo de bases de datos de gran tamaño. Eso sí, afirman que no es una implementación comercial del lenguaje R, sino un programa, de "precio modesto", para integrar funciones de R y paquetes en un formato que permite su ejecución "fácil y eficaz".

Yo creo, confieso que sin haberlo probado, que se trata de aprovechar los nuevos procedimientos desarrollados en R y utilizarlos en el nuevo producto. Los directivos de SPSS Inc. reconocen implícita y explícitamente que en el campo de la investigación en estadística, pero también en su aplicación a las ciencias básicas y experimentales, R es el lenguaje de mayor éxito.
Si no puedes vencer a tu competidor, ...

SystemRescueCD

2009-09-01T15:37:00.002+02:00

Crea una copia de seguridad de una partición del disco duro.
No os ha pasado, que en ciertas situaciones como actualizaciones o fallos técnicos, habéis pensado: ¿Porqué no tengo una copia de seguridad?
En el caso de los archivos asociados al sistema operativo, la cosa es especialmente delicada. No es suficiente con tener alguna copia de los archivos en sí. Debemos tener una copia de todo el sistema, es decir, de la partición donde se encuentra el sistema. Normalmente, en sistemas como Linux se dispone de una partición específica "/" para los archivos del sistema y otra para los documentos "/home". En Windows también es conveniente esta separación pero no todo el mundo la hace.

SystemRescueCD es un LiveCD de Linux (también en lápiz USB) con algunos programas imprescindibles para administrar y reparar particiones del disco duro como Gparted o Partimage. No requiere ninguna instalación, simplemente arrancar el ordenador desde el CD-ROM y puede utilizarse para administrar tanto servidores, como PC's domésticos (con Linux, Windows, etc). El kernel soporta la mayoría de sistemas de archivos como ext2/ext3/ext4, reiserfs, reiser4, btrfs, xfs, jfs, vfat, ntfs (windows xp), iso9660 e incluso sistemas de red (samba y nfs).

En el blog de forat.info tenéis un tutorial con imágenes muy bueno para hacer una copia de seguridad de, por ejemplo, la partición del sistema "/" (o /dev/sda1 en el ejemplo):

Como crear y restaurar copias de seguridad con SystemRescueCD y partimage

El único detalle que puedo mejorar de ese tutorial es que creo que es necesario montar la partición donde grabaremos el archivo antes de lanzar el programa Partimage:

mount /dev/sda5 /mnt/backup

y luego escribir todo el path del archivo donde creamos o rescatamos la copia

* Image file to create/use
/mnt/backup/luiscarimage

Aquí tenéis la web oficial de donde podéis descargaros la ISO del CD. Grabarla y hacer una copia de seguridad. La tranquilidad posterior es infinita.

Emacs-snapshot

2009-08-31T17:39:00.003+02:00

He estado a punto de titular este artículo Pretty emacs, como el artículo de Alexandre Vassalotti. Creo que es muy acertado, pero he preferido no copiarlo.

Emacs-snapshot es un conjunto de paquetes de la versión CVS de GNU Emacs, el conocido editor (y otras cosas), especialmente concebidos para la distribución de linux Debian y, por extensión, para Ubuntu. Han preparado dos versiones. Una para el sistema gráfico GTK y la otra para la consola. Yo utilizo la primera.
A poco que se navegue un poco por la red, todo el mundo coincide que Emacs es un potentísimo editor (dicen que incluso una forma de vida), pero también hay coincidencia en la dificultad de su aprendizaje. Así, muy pocos serán los que se atreban a probarlo. En Kubuntu hay una alternativa que es RKWard, mucho más sencillo i intuitivo, pero también ligeramente inestable. De modo que yo me animé a probar Emacs-snapshot.

Para trabajar con R, la instalación de Emacs-snapshot es muy sencilla:

sudo apt-get install emacs-snapshot ess

El segundo paquete es Emacs speaks statistics, imprescindible para que Emacs y R se entiendan (se supone que R está instalado).
Con esta instrucción se instalaran las versiones de los repositorios oficiales de Ubuntu. Si alguien desea lo último de lo último, hay que instalar los repositorios PPA for Ubuntu Emacs Lisp.
Una vez instalados, cuando abrimos un archivo con extensión .R, Emacs-snapshot lo reconoce como tal y la primera vez que utilicemos los botones de ejecución, abrirá una ventana para solicitar la carpeta de trabajo. La fijamos y aparecerá un buffer (ventana en el lenguaje de Emacs) con R en la parte inferior (ver la figura de arriba). En esta situación ya podemos ir ejecutando las instrucciones del archivo y viendo los resultados en el buffer inferior. El sistema es sencillo y muy práctico. Y a partir de aquí, buscar un manual de Emacs o seguir su tutorial. El esfuerzo vale la pena.

Por cierto, ¿a qué viene lo de pretty emacs?

Pues, resulta que es muy fácil cambiar la fuente y poner la que más nos guste. La que viene por defecto suele ser horrorosa. Ver el post de Alexandre Vassalotti.

¿Y los que trabajamos con Windows? Pues a la izquierda tenéis un enlace a la, a mi parecer, más sencilla instalación para Windows.

Disposición de los gráficos (2)

2009-08-29T11:11:00.015+02:00

Animado por los comentarios al anterior artículo sobre la disposición de los gráficos con ayuda de la función layout, vamos a seguir explicando otras características de esta función.

Ya indiqué que es posible dejar una casilla en blanco, basta poner un cero en la matriz, pero también es posible controlar las alturas o anchuras en centímetros. Para ello se utiliza la función lcm(). Veamos un ejemplo.

Con este código

> layout(matrix(c(1,0,2), ncol=1), heights=c(2, lcm(0.5), 1))

se consigue esta figura

en el que tenemos 0.5 cm entre los dos gráficos.

Otro aspecto a tener en cuenta es la relación de la altura con la anchura. Si queremos que ésta sea de 1:1, lo podemos indicar con el parámetro respect=T.

> layout(matrix(c(1,0,2), ncol=1), heights=c(2, lcm(0.5), 1), respect = T)

El resultado es

Si en la matriz se repite algún número, entonces el gráfico ocupará esas casillas (consecutivas).
Ejemplo:

> layout(rbind(c(1,2),
              c(0,0),
              c(3,3)),
        heights=c(2,lcm(0.5),1),
        respect=T)

Observemos que en esta figura la relación de aspecto se fija en 1 para la altura de la tercera fila (gráfico 3), de modo que el parámetro respect=T hace que la anchura de una celda de la tercera fila también sea 1 (de las dos que hay).

Finalmente podemos decir que el parámetro respect también admite como valor una matriz del mismo tamaño que la que indica la layout, pero con ceros y un uno. El uno marca la posición de la celda para la que hay que respetar la relación de aspecto.

Más ejemplos e información en el libro de Paul Murrell R Graphics Ed. Chapman & Hall/CRC.

Por cierto, los gráficos que muestran la layout actual se consiguen con la instrucción

> layout.show(n)

donde n es el número máximo de gráficos.

Iconos para carpetas de R

2009-08-26T23:20:00.013+02:00

Estos días de vacaciones he aprovechado para estudiar un curso on-line de Gimp (programa GNU para crear gráficos). Es un curso básico pero muy bien hecho que podéis encontrar en el siguiente enlace:

Curso de El Gimp

Como ejercicio, para probar mis conocimientos de capas y manipulación de las mismas, he preparado un icono en formato PNG para el tema Oxigen de KDE 4. Creo que ha quedado bastante digno. Lo podéis utilizar en vuestro sistema operativo Linux con KDE 4 (yo utilizo Kubuntu).

Para aquellos que resisten con Windows XP también he preparado un conjunto de iconos de diferentes tamaños (32x32, 64x64 y 128x128) que podéis descargar del siguiente enlace

folder-R.zip Los iconos de Windows Vista y 7 son más sofisticados. ¿Alguien se anima?

Por cierto, faltaba una versión para la gente que utiliza Gnome en Ubuntu:

Disposición de los gráficos

2009-08-26T10:29:00.018+02:00

En el paquete UsingR tenemos la función gráfica simple.hist.and.boxplot que produce la figura que mostramos con el siguiente código:

> simple.hist.and.boxplot function (x, ...)
{
op <- par(no.readonly = TRUE)
on.exit(par(op))
layout(matrix(c(1, 2), 2, 1), heights = c(3, 1))
par(mai = c(1, 1, 1, 1)/2)
hist(x, xlab = FALSE, col = gray(0.95), yaxt = "n", ...)
rug(x)
boxplot(x, horizontal = TRUE)
}

Observemos especialmente que se trata de combinar dos gráficos, el histograma y el boxplot. Pero el detalle que queremos destacar en este artículo es la disposición de los gráficos.

Hay dos parámetros que sirven para disponer varios gráficos en una misma figura. Por ejemplo, en la figura del artículo Chiplot se han dibujado dos gráficos de dispersión en horizontal. Esto se consigue con la instrucción

oldpar <- par(mfrow = c(1,2))

El vector c(nr,nc) especifica el número de filas y de columnas y el parámetro mfrow su escritura (por filas). También hay un parámetro mfcol para escribir por columnas.

Recordar que la modificación de un parámetro es para toda la sesión, de manera que si queremos revertir el cambio debemos ejecutar

par(oldpar)

Sin embargo, estos parámetros sólo dibujan gráficos del mismo tamaño. ¿Cómo se consique una disposición como la del gráfico que encabeza este artículo?

La respuesta es la instrucción layout.

Básicamente, en la instrucción layout debemos especificar una matriz de números correlativos 1,2,... (puede incluir ceros) que indicarán el orden de escritura. Un cero dejará la posición en blanco. La ventaja es que podemos definir las alturas de las filas y las anchuras de las columnas.

Para lograr la disposición de la figura anterior hacemos

> mm <- matrix(c(1,2), nrow = 2)
> layout(mm, heights = c(3,1))

Observemos que la definición de las alturas es proporcional, es decir, el primer gráfico ocupará 3/(3+1) del total, mientras que el segundo 1/(3+1).

Otros detalles de la función layout los dejaremos para otro artículo.

Colores transparentes

2009-08-23T17:32:00.004+02:00

En algunas ocasiones puede ser de gran utilidad disponer de colores transparentes. En un gráfico como el de arriba, la acumulación de puntos se hace más visible utilizando un color transparente para los puntos.

La función rgb de definición de un color dispone de un parámetro alpha (entre 0 y 1, por defecto 1) que controla el nivel de transparencia. Por ejemplo, podemos definir un rojo transparente así

> rojo20 <- rgb(1,0,0,0.2)

El gráfico de la figura de arriba es

> plot(rnorm(2000), rnorm(2000), pch=16, col=rojo20, xlab="", ylab="")

También se puede definir un color mediante un código de 8 dígitos hexadecimales que empieza por # y cuyos dos últimos corresponden al valor de alpha. Por ejemplo

> rojo21 <- "#FF000050"

Chiplot

2009-08-22T21:12:00.005+02:00

Aunque un gráfico de dispersión (scatterplot) es la primera opción para observar un conjunto bivariante de datos, a veces resulta complicado decidir a primera vista si las variables son independientes o no. Como es lógico nos debemos apoyar en algún contraste estadístico, como el test Ji-cuadrado (discretizando los datos) o el test rho de Spearman, pero también sería de gran ayuda disponer de un gráfico adicional. Éste es el chi-plot de Fisher y Switzer (1985, 2001).

El chi-plot es un gráfico que representa dos variables (lambda y chi) calculadas a partir de los datos muestrales. Tiene la propiedad de que en el caso de que las variables originales sean independientes, los puntos se agrupan en una región central. En caso contrario, si las variables estan relacionadas, los puntos aparecerán fuera de esa región central.

En el gráfico de arriba observamos un conjunto de datos sobre polución que relacionan la cantidad de SO2 con la mortalidad. En el chi-plot apreciamos su dependencia con claridad (rho=0.4915, p-value<0.01).

Fuente: Los datos y la función chiplot se pueden hallar en el libro de B.S. Everitt An R and S-PLUS Companion to Multivariate Analysis y en su web.

Código:
> attach(airpoll)
> chiplot(SO2,Mortality,vlabs=c("SO2","Mortality"))
> cor.test(SO2,Mortality,method="spearman")

Importar datos de una hoja de cálculo

2009-08-21T10:06:00.004+02:00

Una de las primeras preguntas de todas las personas que empiezan a trabajar con R es ¿cómo puedo importar mis datos de Excel? La misma pregunta se hacen los que utilizan otras hojas de cálculo como Calc de OpenOffice.org, Gnumeric de Gnome, KSpread de KOffice, etc.

La solución más sencilla consiste en utilizar la función read.table (o sus variantes read.delim y read.csv) para importar una tabla de datos tipo texto. Es decir, previamente debemos guardar la tabla de datos de la hoja de cálculo como un archivo de texto delimitado con tabulaciones o con comas.
Los pasos son los siguientes:

Abrir la hoja de cálculo con nuestro programa (Excel, Calc,...) y guardar la hoja activa que contenga la tabla como un archivo de texto delimitado. Podremos elegir el tipo de delimitador (tabulación, coma, punto y coma,...) y el nombre y la extensión del archivo. Por ejemplo, podemos guardar un archivo de texto con el separador "," y el nombre "MiTabla.txt" o "MiTabla.csv" (mejor cuando el separador es "," o ";").
Recuperar la tabla desde la consola de R con la instrucción

> mitabla.df <- read.table("MiTabla.csv", header=T, sep=",")

o también con la instrucción equivalente

> mitabla.df <- read.csv("MiTabla.csv")

Cuando el separador es una tabulación, podemos utilizar la instrucción read.delim.

Para todas las personas que tienen la coma "," como separador decimal en sus programas de hojas de cálculo, al guardar el archivo de texto se debe utilizar como separador el punto y coma ";" o la tabulación. Entonces la importación se realiza con las instrucciones:

> mitabla.df <- read.table("MiTabla.csv", header=T, dec=",", sep=";")
o
> mitabla.df <- read.csv2("MiTabla.csv")

Para las tabulaciones es

> mitabla.df <- read.table("MiTabla.txt", header=T, dec=",", sep="\t")
o
> mitabla.df <- read.delim2("MiTabla.txt")

Truco: También es posible recuperar una tabla desde el portapapeles de nuestro sistema operativo con el parámetro file="clipboard" en lugar del nombre del archivo. ¡Pruebalo!

En resumen, hay que guardar la tabla como archivo de texto con separadores y tener en cuenta qué tipo de punto decimal tenemos y qué separador utilizamos.

Para los usuarios de Windows y Excel, también es posible abrir una conexión con el archivo de las hojas de cálculo y seleccionar columnas o tablas de alguna de las hojas que contiene. Para ello debemos utilizar odbcConnectExcel del paquete RODBC.

Más información en R Data Import/Export

Múltiples archivos de gráficos

2009-08-07T18:40:00.002+02:00

El otro día comentaba los múltiples formatos para guardar un gráfico generado por R (ver Devices).
Por ejemplo, en el caso de pdf(), una vez abierto, todos los gráficos producidos se guardan en el mismo archivo, hasta que lo cerremos. Pero ¿como podemos guardar cada gráfico en un archivo diferente? Pues con

pdf(onefile=F)

Esto llevará a escribir diferentes archivos del tipo Rplot001.pdf, Rplot002.pdf,...

Otros parámetros controlan el nombre del archivo y los dígitos. Para ver todas las opciones leer atentamente el help(pdf). No olvidemos cerrar el dispositivo con dev.off().

Por otra parte, en Windows, una vez que tenemos un gráfico en su ventana, lo podemos guardar a través del menú de la ventana en cualquier formato o incluso en el portapapeles. Sin embargo, en Linux es posible que no tengamos ese menú (en RKWard sí). La función savePlot() soluciona este inconveniente:

savePlot("MiFigura", type="png")

Otros formatos y posibilidades se pueden consultar en la ayuda de esta función.

polygon()

2009-08-06T17:33:00.005+02:00

Un truco para rellenar regiones a destacar en un gráfico es utilizar la función polygon().
El gráfico anterior se consigue con el siguiente código:

curve(exp,0.5,3.5,lwd=2)

t<-seq(1,3,by=0.01)
x <- c(1,t,3)
y <- c(0,exp(t),0)
polygon(x,y,density=20)

title("Exponencial entre 1 y 3")

Observemos la utilización de los vértices adecuados para reseguir la curva exponencial.
El parámetro density=20 controla la densidad de las linias (por pulgada). Tambíen es posible rellenar el polígono con un color sólido como en el siguiente ejemplo con la densidad normal:

El código para dibujarlo es:

curve(dnorm,-3,3,axes=F) # dibujamos la curva pero no los ejes
axis(1) # ahora los ejes
axis(2)

t <- seq(0,1,by=0.01)
x <- c(0,t,1)
y <- c(0,dnorm(t),0)
polygon(x,y,col="grey") # la región en gris

curve(dnorm,-3,3,lwd=2,add=T) # la curva encima
abline(h=0)
title("Probabilidad entre 0 y 1")

Devices

2009-08-05T18:34:00.005+02:00

Normalmente en R, un gráfico nuevo sustituye al anterior en la misma ventana, si está abierta. Pero a veces nos puede interesar dibujar el gráfico en una ventana nueva, por ejemplo para comparar los dos gráficos. Para abrir una nueva ventana gráfica podemos invocar la instrucción

x11() # en Linux

windows() # en Windows

quartz() # en Mac OS

Aunque creo que x11() funciona en los tres sistemas.

Una ventana no es más que uno de los posibles dispositivos (devices) de salida para los gráficos.
Para ver los dispositivos adicionales de salida gráfica de que disponemos en un momento dado, podemos hacer

dev.list()

y para elegir uno entre ellos (activamos el tercero, si existe):

dev.set(3)

Sabremos cual es la ventana activa ya que aparecerá la palabra ACTIVE en la barra de títulos.

Además, también es posible redirigir la salida gráfica a un archivo PDF. La siguiente instrucción abre un archivo PDF en la carpeta actual de trabajo con el nombre especificado:

pdf("figura1.pdf")

Mientras esté activo, todos los gráficos que generemos se enviarán sistemáticamente a dicho archivo. Para cerrar el dispositivo activo (dev.cur()) hacemos

dev.off()

Las funciones postscript(), png(), jpeg() y bmp() son similares a pdf(), pero crean gráficos en los formatos EPS, PNG, JPEG y BMP respectivamente. Para insertar un gráfico en LaTeX se recomienda el formato EPS, mientras que en un procesador de textos como Write de OpenOffice o Word de Microsoft es mejor PNG. Recordar que hay que cerrar estos dispositivos para poder acceder al archivo correspondiente.

Box plots

2009-08-04T19:13:00.007+02:00

Un boxplot es un gráfico robusto que permite describir la distribución de una muestra de valores numéricos y compararla con otras (en las mismas unidades).
La robustez se debe a la utilización de estadísticos de orden. En este sentido y como medida de centralidad se utiliza la mediana, en lugar de la media.

En cuanto a la caja, conviene aclarar que no son los cuartiles los que la definen sino las bisagras (hinges), inferior y superior. La bisagra inferior es la mediana de los valores situados por debajo de la mediana de la muestra. La bisagra superior es la mediana de los valores situados por encima de la mediana de la muestra. Para muestras grandes, la bisagra inferior prácticamente coincidirá con el primer cuartil y la bisagra superior con el tercer cuartil, pero no son iguales.
En el boxplot se utilizan las bisagras por definición. Los cuartiles (y otros percentiles) se calculan de forma muy diversa y no son adecuados para fijar una definición única.
Como es natural la función boxplot() de R tiene en cuenta estos detalles y es absolutamente rigurosa con la definición de Tukey.

Otro detalle es la longitud de los bigotes (whiskers). Éstos se extienden fuera de la caja hasta el valor muestral más lejano que no supere una proporción del rango intercuartílico (calculado con las bisagras). Esa proporción se controla con el parámetro range. Su valor por defecto es 1.5, fijado así para que el número de outliers en una población normal sea relativamente pequeño. En otros programas de software estadístico como SPSS, incomprensiblemente, este valor es 2.

Por último, decir que hay otro parámetro varwidth que permite controlar el ancho de las cajas. Si lo pasamos a TRUE, el ancho de las cajas en un boxplot múltiple será proporcional a la raíz cuadrada de los tamaños muestrales de cada población. Interesante.

R Graph Gallery

2009-08-04T18:30:00.003+02:00

En la R Graph Gallery encontraréis un buen conjunto de gráficos, todos ellos con su correspondiente código y listos para ser utilizados.
En la imagen podéis ver el resultado de la función stars() aplicado a un conjunto de variables numéricas sobre la base de datos mtcars.

LaTeX y R

2009-08-03T16:53:00.005+02:00

Si queremos elaborar un documento dinámico o programable que combine texto y resultados estadísticos, tanto numèricos como gráficos, la mejor opción es utilizar LaTeX y R.

Para quien no conozca LaTeX diré que es un sofisticado procesador de textos que se utiliza especialmente con documentos científicos. Aunque su aprendizaje no es sencillo, vale la pena dedicarle un tiempo. Cuando se domina, los otros procesadores sólo se utilizan con textos elementales.
Para aprender LaTeX se puede empezar por la web de CervanTeX y leer el manual La introducción no-tan-corta a LaTeX2e.

Pero si ya sabes un poco de LaTeX, combinarlo con R es muy fácil. La idea es escribir un archivo con extensión .Rnw con código LaTeX normal y texto y añadirle trozos de código R (chunks).

El texto normal empieza con el signo "@", mientras que el código R empieza con "<<>>="

<<Argumento1=ValorDelArgumento1, Argumento2=ValorDelArgumento2, ...>>=
Código R

@
Texto normal y código LaTeX

Una vez redactado el documento combinado, se utiliza la función Sweave() de R para procesarlo y obtener un documento LaTeX, con extensión .tex, que contiene el texto original y los resultados, cálculos y gráficos demandados por el código R insertado. Posteriormente, el archivo .tex se procesará como es habitual en LaTeX para obtener el documento definitivo: .ps, .pdf,...

La función Stangle() permite extraer todo el código R del documento combinado (los chunks) en un archivo con extensión .R.

Un sencillo manual de Sweave.
Leer más en http://www.traba.org/wikitraba/index.php/Sweave

¿En qué SO estamos?

2009-08-02T18:02:00.006+02:00

¡Vaya tontería! Todo el mundo sabe con qué sistema operativo trabaja.
Pero resulta que si queremos escribir un programa en R que se pueda utilizar en varios sistemas operativos, entonces la pregunta es pertinente. Esto no significa que las instrucciones de R dependan del SO en el que trabajemos. Pero hay situaciones en las que es conveniente cambiar alguna cosa en función del SO, por ejemplo, para fijar la carpeta de trabajo.

El objeto version nos proporciona una lista con los datos de nuestra versión de R. En particular, el SO con el que trabajamos.

El siguiente código, fija la carpeta de trabajo en función del sistema operativo:

linuxPath <- "/home/paco/estudios/"
windowsPath <- "c:/users/paco/documents/estudios/"
estudio <- "Blog"

SO <- version[["os"]]

if (SO=="linux-gnu"){
workingDir <- file.path(linuxPath, estudio)
}else{
workingDir <- file.path(windowsPath, estudio)
}

setwd(workingDir)

Cálculo matricial

2009-08-02T17:04:00.004+02:00

Trabajar con matrices en R es muy sencillo.
En primer lugar, para crear una matriz podemos transformar un vector:

> x <- 1:9
> X <- matrix(x, ncol=3)

Por defecto el relleno de la matriz se hace por columnas. Si queremos que sea por filas, hay que indicarlo:

> Y <- matrix(x, ncol=3, byrow=T)

En cuanto a las operaciones con matrices, la suma de dos matrices se consigue con el operador "+"

> X + Y

El producto ordinario se puede hacer con el operador "*", pero éste sirve para calcular el producto término a término. No es el producto matricial.
El producto matricial se obtiene con el operador "%*%":

> X %*% Y

Recordar que es necesario que el número de columnas de la primera matriz coincida con el número de filas de la segunda.

La función t() calcula la matriz traspuesta:

> A <- t(X) %*% X

La función diag() se puede utilizar en dos sentidos, para construir una matriz diagonal o para extraer la diagonal de una matriz cuadrada:

> diag(1:3)
> d <- diag(X)

El determinante de una matriz cuadrada se consigue con la función det(). En cuanto a la inversa de una matriz, si tiene, se calcula con la función solve():

> b <- c(1,0,2,-1,1,0,0,0,1)
> B <- matrix(b, ncol=3)
> det(B) # no es cero, luego es inversible
> solve(B)

Si lo que queremos son los valores y los vectores propios de una matriz, debemos utilizar la función eigen() que nos da una lista:

> r.list <- eigen(X)
> valor.propio1 <- r.list$values[1]
> vector.propio1 <- r.list$vectors[ ,1]

Finalmente, si lo que queremos es la g-inversa de una matriz, podemos utilizar la función ginv() del paquete MASS que calcula la g-inversa de Moore-Penrose:

> library(MASS)
> ginv(X)

Instalación del programa R

2009-08-01T12:26:00.010+02:00

Vamos a empezar por el principio.
Se trata de descargar e instalar el programa R en nuestro ordenador. Para ello nos debemos dirigir al servidor más cercano del CRAN. El CRAN o The Comprehensive R Archive Network es un sistema mundial de servidores que permiten acceder al mismo contenido (mirrors) en todos ellos. En el caso de España es lógico acceder a http://cran.es.r-project.org, pero podemos descargar el programa R de cualquier otro.
El siguiente paso depende del sistema operativo de nuestro ordenador.

R para Windows

Si nuestro SO es Windows (XP, Vista) de 32 bits debemos descargarnos el archivo de nombre R-X.X.X-win32.exe (donde X.X.X es el número de versión, 2.9.1 por ejemplo) que encontraremos en el subdirectorio base y que podemos guardar en alguna carpeta de nuestro disco o directamente ejecutar, para una instalación inmediata.
Durante la instalación, en general, debemos aceptar todas las opciones por defecto, pero hay una pregunta clave que debemos tener en cuenta. Cuando ofrece utilizar las opciones de configuración y por defecto nos dice que no, es mejor cambiar a si. Entonces podemos elegir la opción SDI (ventanas separadas) que está por defecto. Esto es absolutamente necesario si vamos a utilizar un interfaz gráfico como Tinn-R. En otro caso, dejaremos la opción a MDI (una ventana). El resto de las preguntas las podemos pasar con sus valores por defecto.

R para Linux

En el caso que tengamos una distribución de Linux, deberemos considerar si queremos descargarnos el código fuente y compilarlo o, mucho más fácil, descargar e instalar los paquetes ya compilados de nuestra distribución de Linux. La primera opción no es difícil, pero se considera para "expertos". La explicaremos en otro escrito más adelante. Ahora nos concentraremos en la segunda opción para una distribución Ubuntu que es una de las distribuciones de Linux más populares. Otras distribuciones para las que se dispone de paquetes específicos son Debian, Redhat i Suse.

En el caso de tener una distribución Ubuntu (en mi caso Kubuntu) debemos dirigirnos a la página correspondiente
http://cran.es.r-project.org/bin/linux/ubuntu/
donde encontraremos todos los detalles. Pero en realidad nos bastará con seguir las siguientes instrucciones:

Añadir a la lista de repositorios la entrada:
deb http://cran.es.r-project.org/bin/linux/ubuntu jaunty/
Ejecutar en una consola las instrucciones:

sudo apt-get update
sudo apt-get install r-base

Está claro que debemos apuntar al repositorio más cercano y que hay que poner la versión de Ubuntu que tengamos en nuestro ordenador (en mi caso Jaunty).
Si no añadimos este repositorio específico, también se puede instalar R con la instrucción de consola, pero instalaríamos una versión que no sería la más reciente.
Otros detalles de esta instalación los comentaremos más adelante en otro escrito.

R para Mac OS X

Para el sistema operativo de los Mac encontraremos un archivo autoinstalable del tipo R-X.X.X.dmg. La instalación es trivial.