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
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
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
> 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
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
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:
<<
@
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
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.
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
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
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
sábado, 24 de octubre de 2009
Colores
lunes, 12 de octubre de 2009
Símbolos para los datos
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
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:
> 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
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
- 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.
domingo, 6 de septiembre de 2009
Estimación de una función de densidad bivariante
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
martes, 1 de septiembre de 2009
SystemRescueCD
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
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)
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.