R

Durante la charla de Carlos Ortega del pasado jueves sobre el paquete caret y sus concomitancias, se planteó el asunto de la optimización de los parámetros de un modelo usando rejillas (grids) de búsqueda. Cuando un determinado algoritmo depende de, p.e., cuatro parámetros, se puede definir una rejilla como en gbmGrid <- expand.grid(interaction.depth = c(1, 5, 9), n.trees = (1:30)*50, shrinkage = 0.1, n.minobsinnode = 20) y caret se encarga de ajustar el modelo bajo todas esas combinaciones de parámetros (90 en el ejemplo) para ver cuál de ellas es, con las debidas salvedades, óptima. ...

Sí, de drogas de las que mantienen despierto al lumpenazgo. Porque he encontrado (aquí) un conjunto de datos muy interesante sobre la valoración que una serie de personas, unas 900, da a una serie de drogas más o menos legales que se llaman —me acabo de enterar— nootrópicos. El gráfico extraído de la página enlazada más arriba resume parte de los resultados. No obstante, es sabido entre los que se dedican a los sistemas de recomendación que hay usuarios que tienden a valorar sistemáticamente por encima de la media y otros, por debajo. En los manuales de la cosa suelen recogerse mecanismos más o menos sofisticados para mitigar ese efecto y normalizar las valoraciones entre usuarios. Generalmente, solo exigen matemáticas de bachillerato. Y son meras aproximaciones que no tienen en cuenta circunstancias tales como que puede que un usuario dé valoraciones bajas solo porque evalúa productos malos, etc. ...

Imaginemos un banco que construye modelos para determinar si se concede o no un crédito. Este banco tiene varias opciones para crear el modelo. Sin embargo, en algunos países el regulador exige que el banco pueda explicar el motivo de la denegación de un crédito cuando un cliente lo solicite. Esa restricción impediría potencialmente usar modelos de caja negra como el que construyo a continuación: library(randomForest) raw <- read.table("http://archive.ics.uci.edu/ml/machine-learning-databases/credit-screening/crx.data", sep = ",", na.strings = "?") dat <- raw dat$V14 <- dat$V6 <- NULL # me da igual dat <- na.omit(dat) # ídem modelo <- randomForest(V16 ~ ., data = dat) Fijémonos en el sujeto 100, a quien se le deniega el crédito (suponiendo, ¡mal hecho!, que el punto de corte de la probabilidad para concederlo es el 50%), y la variable $V8$. Podemos ver cuál sería el score del cliente modificando esa variable entre su valor real y el máximo del rango dejando las demás tal cual: ...

Leyendo por ahí, he encontrado un comentario sobre el paquete RJSONIO de R en el que se recomendaba no usarlo por no ser libre. El paquete, aparentemente, está liberado bajo una licencia BSD. Pero su pecado es que dentro de uno de los ficheros que contiene, src/JSON_parser.c, dice The Software shall be used for Good, not Evil. Más información, aquí. No sé qué pensar sobre toda esta historia.

En varias de las ediciones de las Jornadas de Usuarios de R he formado parte del comité organizador, que se encarga, fundamentalmente, de la logística de la cosa. Este año, para variar, estoy en el comité científico. Como integrante del cual, es labor mía tratar de animaros a que enviéis alguna propuesta de participación, que puede tener alguno de los siguientes formatos: Una presentación de unos 20 minutos, mostrando alguna aplicación de R. Y no necesariamente en el mundo académico. Son también bienvenidas y apreciadas las aplicaciones en empresas e instituciones de todo tipo. De hecho, una de las presentaciones más recordadas del año pasado, la de Antonio Sánchez Chinchón, trató de aplicaciones ludicomatemáticas de R. Un taller (típicamente de 2 horas) para enseñar a utilizar alguna herramienta particular de R. En el pasado las ha habido de mapas, de gráficos, de Spark… ¡y recuerdo uno sobre plyr y reshape2 que impartí en 2010 cuando esos paquetes eran una novedad y una rareza! Hay tiempo hasta el primero de mayo para realizar las propuestas. Los detalles pueden consultarse aquí. ...

y <- c(rnorm(1000), rnorm(2000, 1, 0.5)) es una mezcla de dos normales (N(0, 1) y N(1, 0.5)) con pesos 1/3 y 2/3 respectivamente. Pero, ¿cómo podríamos estimar los parámetros a partir de esos datos? Se puede usar, p.e., flexmix, que implementa eso del EM. Pero en el librillo de este maestrillo dice library(rstan) y <- c(rnorm(1000), rnorm(2000, 1, 0.5)) codigo <- " data { int<lower=1> K; // number of mixture components int<lower=1> N; // number of data points real y[N]; // observations } parameters { simplex[K] theta; // mixing proportions real mu[K]; // locations of mixture components real<lower=0> sigma[K]; // scales of mixture components } model { real ps[K]; // temp for log component densities sigma ~ cauchy(0,2.5); mu ~ normal(0,10); for (n in 1:N) { for (k in 1:K) { ps[k] <- log(theta[k]) + normal_log(y[n],mu[k], sigma[k]); } increment_log_prob(log_sum_exp(ps)); } }" fit <- stan(model_code = codigo, data = list(K = 2, N = length(y), y = y), iter=48000, warmup=2000, chains=1, thin=10) En el código anterior no sé si queda claro cómo cada punto $y_i$ sigue una distribución (condicionada a los parámetros) con densidad $\theta_1 \phi(y_i, \mu_1, \sigma_1) + \theta_2 \phi(y_i, \mu_2, \sigma_2)$. ...

Supongo que no debería escribirlo aquí sino comunicárselo a quien mantiene ggmap. Pero ya tuve una experiencia mejorable con él y dos no serán. Así que lo cuento por acá. La mayor parte del mérito en el descubrimiento, en cualquier caso, es de una alumna de la clase de R que he dado hoy (en el momento en el que escribo, no en el que lees) en el Banco de Santander. No tengo su nombre ni tengo claro que quisiese que lo mencionase. ...

Fácil: library(maptools) library(ggmap) # un fichero bajado del Ayto. de Madrid # (catálogo de datos abiertos) rutas <- getKMLcoordinates("dat/130111_vias_ciclistas.kml") # procesando el fichero kml rutas <- lapply(1:length(rutas), function(x) data.frame(rutas[[x]], id = x)) rutas <- do.call(rbind, rutas) # mapa de Madrid mapa <- get_map("Madrid", source = "stamen", maptype = "toner", zoom = 12) # pintando los tramos sobre el mapa ggmap(mapa) + geom_path(aes(x = X1, y = X2, group = id), data = rutas, colour = "red") produce Nota: KML es esto.

Por si alguien aún no lo sabe: estamos todos citados en Albacete los días 17 y 18 de noviembre de 2016 en las VIII Jornadas de Usuarios de R. Los detalles, aquí.

Está de moda usar caret para estas cosas, pero yo estoy todavía acostumbrado a hacerlas a mano. Creo, además, que es poco instructivo ocultar estas cuestiones detrás de funciones de tipo caja-negra-maravillosa a quienes se inician en el mundo de la construcción y comparación de modelos. Muestro, por tanto, código bastante simple para la validación cruzada de un modelo con R: # genero ids ids <- rep(1:10, length.out = nrow(cars)) # Nota: da igual si nrow(df) no es múltiplo de 10 # los aleatorizo ids <- sample(ids) # esto devuelve una lista de dfs: preds.cv <- lapply(unique(ids), function(i){ preds <- predict(lm(dist ~ speed, data = cars[ids != i,]), cars[ids == i,]) data.frame( preds = preds, real = cars[ids == i,]$dist) }) # "apilo" los dfs: preds.cv <- do.call(rbind, preds.cv) # calculo el rmse rmse <- sqrt(mean((preds.cv$preds - preds.cv$real)^2)) Sí, estoy usando el RMSE aunque sea un detractor del mismo. ...