R

Fui un pájaro mañanero con [plyr](http://cran.r-project.org/web/packages/plyr/index.html). Probé una vez [data.table](http://cran.r-project.org/web/packages/data.table/index.html) y no me convenció. Volví a él cuando realmente lo necesitaba y ahora es la prolongación de mis dedos. Aún no me he puesto con [dplyr](http://cran.r-project.org/web/packages/dplyr/index.html) aunque he visto el suficiente código escrito con él que no creo que me cueste mucho comenzar a usarlo. Pero tengo la sensación de que tenemos un cisma como el de vi contra emacs en ciernes. Comienza a haber, parece, partidarios acérrimos de tirios y troyanos. Así que abro la sección de comentarios para que opines sobre estos paquetes. A mí y a muchos otros lectores nos gustaría conocer tu opinión al respecto. ¿Cuál utilizas? ¿Qué te gusta de cada cual? ¿Cuál recomendarías?

Voy a escribir sobre un artículo como no debe hacerse: sin haberlo leído. Los bayesianos dirían que esta opinión que aquí voy a vertir es mi prior para cuando encuentre el tiempo y bajo la cual matizaré lo que en el se diga. Lo advierto, en todo caso, para que quien me lea no renuncie al sanísimo escepticismo. Voy a hablar de Inferring causal impact using Bayesian structural time-series models y del paquete de R que lo acompaña, CausalImpact, cuyos autores trabajan en Google. ...

Tuve que saltarme una diapositiva en el DataBeers de Madrid del pasado jueves. (A propósito, aquí están las 1+20 diapositivas.) La decimonona, de la que trata la entrada, viene a hablar de lo siguiente. Tenemos una base de datos con sujetos (ids) que hacen cosas en determinados momentos. No es inhabitual calcular la frecuencia de esos sujetos así: select id, count(*) as freq from mytabla where fecha between current_date - 7 and current_date group by id ; Esa variable se utiliza frecuentemente ya sea como descriptor de los sujetos o como alimento de otros modelos. ...

Hoy he encontrado una solución decente a un problema que venía arrastrando desde hace un tiempo en R. Tengo una tabla muy grande (decenas de millones de registros) con su id. Me interesa quedarme con el subconjunto de la tabla original en que para cada id el valor de una determinada variable es mínimo. Un caso de uso: esa variable adicional mide la distancia de la observación a los centroides de unos clústers. El registro con el menor valor proporciona la asignación del sujeto a su grupo. ...

Jesús Fernández Villaverde escribió en Nada es Gratis sobre R el otro día. Publicó cuatro vídeos: Una entrevista a Hadley Wickham La conferencia plenaria de John Chambers en la useR! 2014 Una entrevista a Yihui Xie, el creador de knitr La conferencia plenaria de Martin Mächler en la useR! 2014, Good Practices in R Programming Todos muy interesantes. Pero todavía lo es más que desde unas páginas del impacto de Nada es Gratis tengan a bien contribuir al conocimiento (¡y al uso!) de R.

Una de las mayores contrariedades de estar sentado cerca de alguien que es más matemático que un servidor (de Vds., no de silicio) es que oye siempre preguntar por qué. Una letanía de preguntas me condujo a leer papelotes que ahora resumo. Primero, unos datos: set.seed(1234) n <- 100 x1 <- rnorm(n) x2 <- rnorm(n) x3 <- rnorm(n) y <- 0.3 + 0.2 * x1 + 0.5 * (x2 > 0) + 0.2 * rnorm(n) Luego, un modelo: modelo <- lm(y ~ x1) summary(modelo) # Call: # lm(formula = y ~ x1) # # Residuals: # Min 1Q Median 3Q Max # -0.9403 -0.2621 0.0420 0.2299 0.6877 # # Coefficients: # Estimate Std. Error t value Pr(>|t|) # (Intercept) 0.55632 0.03364 16.538 < 2e-16 *** # x1 0.21876 0.03325 6.579 2.34e-09 *** # --- # Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 # # Residual standard error: 0.3323 on 98 degrees of freedom # Multiple R-squared: 0.3063, Adjusted R-squared: 0.2992 # F-statistic: 43.28 on 1 and 98 DF, p-value: 2.341e-09 Pocos que no entiendan cómo se han generado los datos advertirían lo malo de su especificación: hemos omitido una variable explicativa cuyo efecto ha ido a incrementar el error de manera que los tests habituales de bondad de ajuste no advierten. ...

Dos motivos me han tenido missing estas últimas semanas. Uno es una estancia en la Universidad de Santa Catalina del Burgo de Osma. Oportunamente ubicada en las estribaciones de la muy generosa en caldos de calidad Ribera del Duero, ha sido reconvertida a la sazón en un hotel propicio para la evasión y la agrafía. El segundo es que en horas intempestivas he estado purgando de missings unas matrices enormes y de la, se conoce, mayor trascendencia. Es un asunto delicado, jamás bien resuelto, para el que el paquete [VIM](http://cran.r-project.org/web/packages/VIM/index.html) puede proporcionar ayuda. Sobre todo en los aspectos gráficos. ...

Vi esto y me dije: yo también quiero. Así que dicho y hecho: Por si acaso, cada diagrama representa la descomposición en números primos de un número del 1 al 100. El código (que no he adecentado lo que suelo) es un pequeño ejercicio con el paquete grid y unos elementos de recursividad (como en Grid, Scala y arbolitos fractales): library(grid) library(gmp) plot.factors <- function(n, new.plot = TRUE){ if(new.plot) grid.newpage() divisors <- sort(as.integer(factorize(n)), decreasing = T) foo <- function(divs){ if(length(divs) == 0){ grid.circle(x = 0.5, y = 0.5, r = 0.5, gp=gpar(fill="black")) return() } n <- divs[1] x <- (Re(exp( 2 * pi *(1:n) * 1i /n))) / 4 + 0.5 y <- (Im(exp( 2 * pi *(1:n) * 1i /n))) / 4 + 0.5 for(i in 1:n){ tmp <- viewport(x = x[i], y = y[i], w = 2/(3 + n), h = 2/(3 + n)) pushViewport(tmp) #grid.rect(gp = gpar(col = "grey")) foo(divs[-1]) popViewport() } } foo(divisors) } plot.factors(25) grid.newpage() nrow <- 10 ncol <- 10 for(y in 1:nrow){ for(x in 1:ncol){ tmp <- viewport(x = x / (1 + ncol), y = 1 - y / (1 + nrow), w = 1/(1 + ncol), h = 1/(1+ncol)) pushViewport(tmp) #grid.rect(gp = gpar(col = "grey")) plot.factors(x + y * ncol - ncol, new.plot = FALSE) popViewport() } }

Lo que escribí hace un par de días sobre procesos puntuales, ahora me doy cuenta, podía haberse resuelto con nuestro viejo amigo glm. Ejecuto el código del otro día y obtengo (para un caso nuevo) mu alfa verosimilitud delta 1 0.4493158 0.50000000 340.6141 1 2 0.2675349 0.40457418 307.3939 2 3 0.1894562 0.28917407 293.4696 3 4 0.1495654 0.22237707 287.0784 4 5 0.1243791 0.18079703 281.3900 5 6 0.1142837 0.14913172 284.9227 6 7 0.1217504 0.12150745 288.5448 7 8 0.1214365 0.10424818 289.3282 8 9 0.1204605 0.09148817 290.9081 9 10 0.1315896 0.07857330 295.3935 10</code> que significa que el parámetro óptimo es delta = 5, mu = 0.124 y alfa = 0.18. Ahora hago cuantos.previos <- function(i, muestra, delta){ indices <- Filter(function(x) x < i & x > i - delta, 1:n) cuantos <- sum(muestra[indices]) } fit.glm <- function(delta){ prev <- sapply(1:length(muestra), cuantos.previos, muestra, delta) dat <- data.frame(muestra = muestra, prev = prev) res.glm <- glm(muestra ~ prev, data = dat, family = poisson(link = "identity")) c(delta, res.glm$coefficients, summary(res.glm)$aic) } res.glm <- sapply(1:10, fit.glm) res.glm <- as.data.frame(t(res.glm)) colnames(res.glm) <- c("delta", "mu", "alfa", "aic") y obtengo ...

Tengo una serie de datos que se parecen a lo que cierta gente llama procesos puntuales y que se parecen a los que se introducen (muuuuy prolijamente) aquí. Gráficamente, tienen este aspecto: Sobre un determinado periodo de tiempo (eje horizontal) suceden eventos y los cuento por fecha. Pero no suceden independientemente (como si generados por un proceso de Poisson) sino que tienden a agruparse: el que suceda un evento tiende a incrementar la probabilidad de que suceda otro poco después. El proceso, en una mala traducción, se autoexcita. ...