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