R

El otro día escribí sobre visitas a la Wikipedia. El otro día (posiblemente otro) oí hablar de prophet.

Hoy con

library(wikipediatrend)
library(prophet)
library(ggplot2)

visitas <- wp_trend(
    "R_(lenguaje_de_programaci%C3%B3n)",
    from = "2010-01-01", to = Sys.Date(),
    lang = "es")

mis.visitas <- visitas[, c("date", "count")]
colnames(mis.visitas) <- c("ds", "y")

pasado <- mis.visitas[1:1500,]
m <- prophet(pasado)

futuro <- make_future_dataframe(m,
    periods = nrow(mis.visitas) - 1500)
prediccion <- predict(m, futuro)

pred.plot <- plot(m, prediccion)
pred.plot +
    geom_line(data = mis.visitas[1501:nrow(mis.visitas),],
        aes(x = ds, y = y), col = "red", alpha = 0.2) +
    xlab("fecha") + ylab("visitas") +
    ggtitle("Predicción de visitas a la página de R\nen la Wikipedia con prophet")

construyo

Hace un tiempo probé el paquete wikipediatrend de R ya no recuerdo para qué. Desafortunadamente, el servicio que consulta debía de estar caído y no funcionó. Ahí quedó la cosa.

Una reciente entrada de Antonio Chinchón en su blog me ha invitado a revisitar la cuestión y ahora, al parecer, stats.grok.se vuelve a estar levantado. Por lo que se pueden hacer cosas como:

visitas <- wp_trend("R_(lenguaje_de_programaci%C3%B3n)",
    from = "2010-01-01", to = Sys.Date(),
    lang = "es")

[Aquí ahorro al lector unos párrafos de pésima literatura.]

El paquete hunspell de R permite procesar texto utilizando como soporte la infraestructura proporcionada por Hunspell, el corrector ortográfico que subyace a muchas aplicaciones en R.

Existe una viñeta que ilustra el uso del paquete pero, como siempre, en inglés. En español las cosas son parecidas pero, como siempre, nunca exactamente iguales. En esta entrada, por lo tanto, voy a repasar partes de la viñeta aplicándolas a nuestra tan frecuentemente maltratada mas por ello no menos querida por algunos como yo (pausa) lengua.

Es esta:

156.67 * 100
# 15667
as.integer(156.67 * 100)
#15666

Claro, hay que leer ?as.integer para enterarte de que, en realidad, la función que quieres usar es round.

Una mala manera de perder un par de horas.

El número efectivo de partidos es el nombre de una página de la Wikipedia, que contiene la fórmula

$$ N = \frac{1}{\sum_i p_i^2}$$

y excipiente alrededor.

Aplicada a España (usando datos del CIS como proxy),

Como casi siempre, el código:

library(rvest)
library(rvest)
library(reshape2)
library(plyr)
library(zoo)

url <- "http://www.cis.es/cis/export/sites/default/-Archivos/Indicadores/documentos_html/sB606050010.html"

raw <- read_html(url)
tmp <- html_nodes(raw, "table")
tmp <- html_table(tmp[[2]], fill = TRUE)

colnames(tmp)[1] <- "partido"

tmp <- tmp[!is.na(tmp$partido),]
tmp <- tmp[1:30,]

tmp <- melt(tmp, id.vars = "partido")
tmp <- tmp[tmp$value != ".",]
tmp$value <- as.numeric(tmp$value)

tmp$variable <- gsub("ene", "01-", tmp$variable)
tmp$variable <- gsub("abr", "04-", tmp$variable)
tmp$variable <- gsub("jul", "07-", tmp$variable)
tmp$variable <- gsub("oct", "10-", tmp$variable)

tmp$variable <- gsub("-0", "-200", tmp$variable)
tmp$variable <- gsub("-1", "-201", tmp$variable)
tmp$variable <- gsub("-9", "-199", tmp$variable)

tmp$variable <- paste0("01-", tmp$variable)

tmp$variable <- as.Date(tmp$variable, format = "%d-%m-%Y")

dat <- tmp

dat <- ddply(dat, .(variable), transform, total = value / sum(value))
res <- ddply(dat, .(variable), summarize, enp = 1 / (sum(total^2)))

res <- zoo(res$enp, order.by = res$variable)

plot(res, main = "Número efectivo de partidos\nen España(1996-2016)",
        xlab = "", ylab = "número efectivo de partidos")

Recibí un mensaje el otro día sobre polinomios monótonos. Mejor dicho, sobre el ajuste de datos usando polinomios monótonos. Frente a un modelo del tipo y ~ x (x e y reales) donde la relación entre las dos variables es

• manifiestamente no lineal y
• necesariamente monótina, p.e., creciente (por consideraciones previas),

cabe considerar ajustar un polinomio monótono, i.e., realizar una regresión polinómica con la restricción adicional de que el polinomio de ajuste resultante sea monótono.

Un año, el 2016, mueren 1160 personas en accidentes de tráfico. El anterior, 1131, i.e., 29 menos. Ruido estadístico aparte, ¿aumentan?

Comenzamos a optar. Primera elección subjetiva: son muestras de una Poisson de parámetro desconocido. La pregunta: ¿el mismo?

Una manera de estudiar lo anterior es plantear

1160 ~ poisson(lambda * (1 + incr))
1131 ~ poisson(lambda)

y estudiar la distribución de incr. Que a saber qué distribución tendrá (teóricamente). Pero, ¿importa?

Mejor que rebuscar a ver qué distribución podría tener la cosa, basta con envolverlo en un poco de seudo-C++,

Vivimos en un mundo opaco: como en los cuentecillos de Asimov, somos usuarios de tecnologías que ni conocemos ni controlamos. Parametrizamos nuestras máquinas y las echamos a correr. Poco más podemos hacer que fiarnos de quienes nos las proporcionan.

Luego pasan cosas como que, de repente, resulta que Stan, en las últimas versiones, ha estado produciendo muestras sesgadas. ¿Qué resultados condicionará eso río abajo?

Un caso mucho más famoso es el de la resonancia magnética (fMRI): un error en el software concomitante pone bajo sospecha hasta 40000 artículos sobre estudios del cerebro. Precisamente, por lo mismo.

Stan. Modelo multinivel. Variable categórica. Codificación con ceros y unos. Matriz. Coeficiente vector[n_ccaa] Cccaa. Sin priori.

Catástrofe:

(Coeficientes hasta 15000. Sin tasa, con tiempo. Los valores desorbitados, en ceros de la dummy).

Priori.

for (i in 1:n_ccaa)
    Cccaa[i] ~ cauchy(0, 20);

¿Por qué no?

Tachán:

(¿Para qué verbos?)

Hay una entrada mía, esta, que me ronda la cabeza y con la que no sé si estoy completamente de acuerdo. Trata de justificar la codificación por impacto de variables categóricas en modelos lineales (generalizados o no) y cuanto más la releo, menos me la creo. O, más bien, comienzo a cuestinarme más seriamente contextos en los que funciona y contextos en los que no.

Pero comencemos por uno simple: los árboles. Es moda pensar que, dado un predictor categórico, un árbol explora todas las permutaciones posibles de categorías y que por eso algunas implementaciones de, por ejemplo, bosques aleatorios no permiten variables categóricas de más de cierto número no particularmente generoso de niveles.

Según las últimas noticias, en España, en 2016 murieron 1160 personas en accidentes de tráfico, 29 más que en 2015. Usa el poiss.test para ver si la diferencia es o no significativa.

Me pasa un exalumno (no os perdáis su blog y su otro blog, especialmente si os interesa la versión friqui de la fotografía) de mi curso de R de KSchool, un gráfico que ha realizado con R para su empresa:

Addenda: En los comentarios hay más enlaces relevantes (proporcionados por el autor de la imagen), que conducen al código, etc.

Acabo de sugerir a mis alumnos de KSchool una lista de problemas después de sus 10 primeras horas de contacto con R.

Uno de ellos, advierto, y les he advertido (porque yo, no siendo rector de universidad pública alguna, no cuento el plagiar entre mis vicios) es una versión de otro publicado aquí.

Ejercicio

Construye una matriz que dada una entrada del tipo

a <- c('NAME:Maria /COUNTRY:uruguay /EMAIL:mariaUY@gmail.com',
'NAME:Paul/COUNTRY:UK /EMAIL:PaulUK@gmail.com',
'NAME:Jhon /COUNTRY:USA /EMAIL:JhonUSA@gmail.com',
'NAME:Carlos /COUNTRY:Spain /EMAIL:CarlosSP@gmail.com')

devuelva un tabla con columnas name, country y email (con los datos correspondientes).