lunes, 21 de diciembre de 2009

Aplicar una función a un vector o una lista


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

viernes, 11 de diciembre de 2009

Fórmulas matemáticas en un gráfico


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.

martes, 8 de diciembre de 2009

Clases S3 y S4


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")

lunes, 30 de noviembre de 2009

Primeras jornadas de R

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í.

sábado, 28 de noviembre de 2009

Factores numéricos

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"

domingo, 22 de noviembre de 2009

Valores perdidos

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

domingo, 8 de noviembre de 2009

De R a ODF: el paquete odfWeave

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.

sábado, 7 de noviembre de 2009

Test de Levene para la igualdad de varianzas

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 Chi-Square(4) (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.

sábado, 31 de octubre de 2009

Contrastes de normalidad


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.

jueves, 29 de octubre de 2009

R 2.10.0

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.

domingo, 25 de octubre de 2009

Interpolación de colores

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

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:

sábado, 24 de octubre de 2009

Colores

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.

lunes, 12 de octubre de 2009

Símbolos para los datos

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



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.

sábado, 10 de octubre de 2009

De R a HTML: el paquete hwriter

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


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!

jueves, 1 de octubre de 2009

Integración de datos


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


domingo, 27 de septiembre de 2009

Cádiz y R-Commander


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.

domingo, 13 de septiembre de 2009

RKWard 0.5.1

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:
  1. 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.
  2. 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.
  3. 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.

domingo, 6 de septiembre de 2009

Estimación de una función de densidad bivariante

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.

sábado, 5 de septiembre de 2009

SPSS + R = PASW Statistics Developer

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, ...

martes, 1 de septiembre de 2009

SystemRescueCD

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.

lunes, 31 de agosto de 2009

Emacs-snapshot

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.

sábado, 29 de agosto de 2009

Disposición de los gráficos (2)

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.