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.

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.

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:

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.

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.