R

Nunca pude ser un erizo. Lo intenté durante años y acabé en el sicólogo. Si el cuerpo me hubiese dado, ahora, tal vez, como algunos compañeros de promoción, sería un experto en un área diminuta del conocimiento y corregiría exámenes los fines de semana. Descubrí con tiempo y muchas sesiones de a 60 euros la hora que había nacido para ser un zorro, un merodeador que olisquea aquí y allá.

Una lectura superficial de estas páginas que escribo puede dar la impresión de que contienen conocimientos profundos (¡ojalá!). Una lectura profunda, que contiene superficialidades (¡convengo!). Son obra de un zorro.

En ocasiones uno quiere añadir un total calculado en ciertos bloques a una tabla. Por ejemplo, en la tabla

set.seed(1234)
ventas.orig <- data.frame(
    cliente = rep(1:10, each = 5),
    producto = rep(letters[1:5], times = 10),
    importe = rlnorm(50))

tenemos clientes, productos e importes. Y nos preguntamos por el porcentaje en términos de importe que cada producto supone para cada cliente.

Una manera natural pero torpe de realizar este cálculo consiste en usar un objeto intermedio y merge:

library(plyr)
tmp <- ddply(ventas.orig, .(cliente),
    summarize, total = sum(importe))
ventas <- merge(ventas.orig, tmp)
ventas$pct.producto <- 100 * ventas$importe /
    ventas$total

No os asustéis, se puede hacer aún peor (p.e., usando sqldf). Pero existen dos maneras, cuando menos, de hacerlo mejor. La primera es usando data.table.

Hoy traigo a mis páginas cuatro enlaces que apuntan a recetarios y tutoriales sobre la solución a cuatro problemas que pueden encontrar los usuarios de R:

• Conectar R y Excel
• Importar grandes ficheros CSV (y falta LaF)
• Integrar R con C/C++
• Paralelizar código con snow

¡Espero que os resulten útiles!

Llamo sospechosos habituales a esos programas y lenguajes para el análisis de datos distintos de R cuya decreciente popularidad nos parece tan natural a los partidarios de este último. Abundan los análisis de cuotas de mercado tales como What Analytic Software are People Discussing?

¿Cuáles son estos sospechosos habituales? Pues SAS, SPSS y algún otro: Stata, Statistica, Minitab,…

Sin embargo, R tiene competidores más serios a medio plazo. Uno de ellos, el más importante, es Python. Lo veo a mi alrededor: son muchos los físicos, los ingenieros, los informáticos que tienen experiencia en ese lenguaje y, sintiéndose cómodos en él —y les alabo el gusto— quieren utilizarlo para analizar datos cuando les toca.

Acaba de publicarse Displaying Time Series, Spatial, and Space-Time Data with R, un libro de Óscar Perpiñán que, conociéndolo como lo conozco, me atrevo a recomendar sin haberlo hojeado siquiera.

Además, Óscar nos ha regalado una guía sobre cómo escribir un libro técnico con Emacs y otras herramientas libres.

El primer enlace de la selección de esta semana es The evolution of the modern census. Todos sabemos que lo que llevó a José y María a Belén hace más de 2000 años fue dizque tenían que censarse. Hay noticias de censos anteriores. Desde entonces hasta ahora ha habido muchos, muchísimos censos, pero su mismo concepto y finalidad ha ido cambiando a lo largo de la historia: ya no se trata solamente de contar, medir la riqueza o el poderío militar. Ahora nos interesan otros aspectos relacionados ya no tanto con el cuántos sino con el cómo somos.

El escritor exemplar es un experimento de escritura automática realizado por Molino de Ideas sobre una idea de Mario Tascón y con la colaboración de Carlos J. Gil Bellosta en conmemoración por los 400 años de la publicación de Las Novelas Ejemplares.

Eso reza el pie de página de El escritor exemplar un artilugio que a veces crea frases tales como

que debieran ser aleatorias, no muy distintas en estilo de las Novelas Ejemplares y, con muchísima suerte, inspiradoras.

Me llegó recientemente un artículo con una lista de veinte paquetes de R para data scientists. Y no la encuentro afortunada. Voy a agrupar esos veinte paquetes en algunas categorías y añadiré comentarios. La primera de ellas es la de manipulación de datos, tal vez la más amplia, que recoge los siguientes: sqldf, plyr, stringr (para procesar texto), lubridate (para procesar fechas),reshape2 y los paquetes de acceso a bases de datos.

Inspirado en esto aunque con la intención de mejorar el horrible código adjunto, escribí el otro día esto:

library("biOps")
library("cluster")

# leo una foto usando readJpeg de biOps
# el objeto devuelto es un array mxnx3 dimensional
# la última dimensión es el rgb de cada pixel

tmp <- tempfile()
download.file("http://blog.guiasenior.com/images/Retrato_Garber.jpg", tmp)
x <- readJpeg(tmp)

# si quieres mostrar la foto como un gráfico...
#plot(x)

# convertimos el array 3D nxmx3 en uno 2D (nm)x3
# luego buscamos 5 clústers
# esencialmente, buscamos 7 "píxels representativos"
d <- dim(x)
clarax <- clara(array(x, dim = c(d[1] * d[2], d[3])), 7)

# reemplazamos cada rgb de cada cluster por su
# "píxel representativo" (medioide) correspondiente
rgb.clusters <- clarax$medoids[clarax$cluster,]

# convertimos la matriz resultante en un array 3D
# (invirtiendo la transformación anterior)
# y representamos gráficamente
plot(imagedata(array(rgb.clusters, dim = d)))

Obviamente, podéis cambiar la foto y hacer variar el número de clústers. Pero conviene recordar que:

Después de hablar con tirios y troyanos sobre mi entrada sobre los efectos de binarizar una variable objetivo continua, he decidido tomarme la justicia por mi mano y llamar a la caballería. Es decir, utilizar random forests.

Aquí va el código:

library(randomForest)

set.seed(1234)

my.coefs <- -2:2
n <- 200
train.n <- floor(2*n/3)

test.error <- function(){
  X <- matrix(rnorm(n*5), n, 5)
  Y <- 0.2 + X %*% my.coefs + rnorm(n)
  Y.bin <- factor(Y>0)

  train <- sample(1:n, train.n)

  X <- as.data.frame(X)
  X$Y <- Y

  modelo <- randomForest(Y ~ .,
    data = X[train,])
  pred <- predict(modelo, X[-train,])
  error.cont <- length(pred) -
    sum(diag(table(pred >0, Y[-train]>0)))

  X$Y <- Y.bin
  modelo <- randomForest(Y ~ .,
    data = X[train,])
  pred <- predict(modelo, X[-train,])
  error.bin <- length(pred) -
    sum(diag(table(pred, Y.bin[-train])))

  data.frame(error.cont = error.cont,
    error.bin = error.bin)
}

errores <- do.call(rbind,
  replicate(1000, test.error(), simplify = F))

sapply(errores, fivenum)

El resultado, si te interesa, en tu pantalla.