R

Esta entrada debería ser un comentario más en esta otra, pero voy a abusar del privilegio de ser dueño de la plataforma para promocionarla.

Voy a decir cosas que son aproximadamente ciertas. Los detalles de la verdad de todo están en la ayuda y el código de sort y sus métodos.

En R hay dos métodos de ordenación: shell y radix. El primero es genérico y el segundo es mejor cuando en el vector hay muchos elementos repetidos (p.e., ordenar el censo por provincias).

Con datos bajados de aquí:

library(MicroDatosEs)
library(dplyr)
library(microbenchmark)
library(ggplot2)

censo <- censo2010("MicrodatosCP_NV_per_nacional_3VAR.txt")

censo_char <- as.data.frame(censo[,
    c("CPRO", "SEXO", "ECIVIL", "FACTOR")])
censo_factor <- censo_char
censo_factor$CPRO <- factor(censo_factor$CPRO)


foo <- function(x)
    x %>% group_by(CPRO) %>%
    summarise(res = sum((SEXO == "Mujer") *
        (ECIVIL == "Divorciado") * FACTOR) /
        sum(FACTOR) * 100)

res <- microbenchmark(
    char = foo(censo_char),
    factor = foo(censo_factor),
    times = 10
)

autoplot(res)

Da:

¿No es sorprendente? De hecho, plyr es más rápido que dplyr en este caso si no se usan factores.

Notas:

• El hilo de por qué es así en lugar de otra manera se pierde en código escrito en C++. Para otra vida (mía o de otro).
• Debo agradecer a Diego Castro el intercambio de ideas, código y perplejidades que dieron pie a todo lo de arriba.

Esta entrada es un (otro, que sumar a este o este) recordatorio de que las XI Jornadas de Usuarios de R están en marcha.

Y que serán en Madrid, del 14 al 16 de noviembre, etc. Información toda ella que los enlaces anteriores extienden debidamente.

(Además hay una tarifa reducida cuyo plazo termina, aviso, muy, muy pronto.)

He construido

que, obviamente no es la gran maravilla, basándome en Rectangular Statistical Cartograms in R: The recmap Package y usando

library(rgdal)
library(pxR)
library(recmap)

provs <- readOGR(dsn = "provincias/",
    layer = "Provincias")

pobl <- as.data.frame(read.px("2852.px",
    encoding = "latin1"), use.codes = T)
pobl2 <- as.data.frame(read.px("2852.px",
    encoding = "latin1"))

pobl$nombre <- pobl2$Provincias

pobl <- pobl[, c("Provincias", "nombre", "value")]
colnames(pobl) <- c("COD_PROV", "nombre", "poblacion")
pobl <- pobl[pobl$COD_PROV != "null",]

pobl <- pobl[!pobl$COD_PROV %in%
    c("51", "52", "38", "07", "35"),]


dat <- merge(provs, pobl,
    by = "COD_PROV", all.x = FALSE)
dat@data$NOM_PROV <- NULL
dat$z <- dat$poblacion

tmp <- as.recmap(dat)

tmp$name <- dat@data$nombre
tmp$ccaa <- dat@data$COD_CCAA

res <- recmapGA(tmp, popSize = 300,
    maxiter = 30, run = 10)

cartogram <- res$Cartogram

ccaa <- tmp[, c("name", "ccaa")]
ccaa$ccaa <- as.numeric(factor(ccaa$ccaa))
cartogram <- merge(cartogram, ccaa)

plot.recmap(cartogram, col.text = "black",
    main = "cartograma -- población\n  españa peninsular",
    col = cartogram$ccaa)

Como los datos los he bajado de por ahí y no recuerdo dónde, dejo como referencia el objeto arriba llamado tmp aquí.

De aquí a cuando lo tenga que usar realmente, seguro que me olvido. Así que retomo el uso original de este blog, que era el de dejarme notas a mí mismo y apunto: usa [flextable`](https://cran.r-project.org/package=flextable).

¿Y por qué?, me preguntaré a mí mismo dentro de unos días. Pues por cosas como esta:

(Claro está, salvo que alguien tenga a bien proponer una alternativa mejor).

Los modelos GARCH son otra de esas cosas de las que oyes hablar y como nunca se convierten en problemas de los de carne en el asador, preocupan poco y ocupan menos (más allá de que sabes que se trata de modelos similares a los de series temporales de toda la vida donde la varianza varía de cierta forma a lo largo del tiempo). Pero comienzas a leer cosas como esta y no te enteras de nada: solo hay letras y llamadas a funciones oscuras y oscurantistas.

Son los enteros del 1 al 35 ordenados de forma que dos consecutivos en la serie suman un cuadrado perfecto. Los he obtenido así:

library(adagio)

foo <- function(n){
    desde <- 1:n
    hasta <- 1:n
    todos <- expand.grid(desde, hasta)
    todos <- todos[todos$Var1 < todos$Var2,]
    todos$sqrt <- sqrt(todos$Var1 + todos$Var2)
    todos <- todos[todos$sqrt == round(todos$sqrt),]
    todos$sqrt <- NULL
    vertices <- as.vector(t(todos))
    hamiltonian(vertices)
}

foo(35)

Notas:

• Esta entrada está inspirada en algo que he visto en Twitter (pero cuya referencia he olvidado guardar).
• Puedes probar con otros números, aunque no siempre existe un ciclo hamiltoniano.

Predicciones puntuales:

O (sub)modelos:

Y parece que ahora también distribuciones:

Notas:

• Obviamente, la clasificación anterior no es mutuamente excluyente.
• La tercera gráfica está extraída de Transformation Forests, un artículo donde se describe el paquete trtf de R.
• Los autores dicen que [r]egression models for supervised learning problems with a continuous target are commonly understood as models for the conditional mean of the target given predictors. ¿Vosotros lo hacéis así? Yo no, pero ¡hay tanta gente rara en el mundo!
• Y añaden que [a] more general understanding of regression models as models for conditional distributions allows much broader inference from such models. Que era lo que creía que todos hacíamos. Menos, tal vez, algún rarito.

Aprovechando que el paquete GameTheoryAllocation ha emergido de mi FIFO de pendientes a los pocos días de conocerse los resultados de las [adjetivo superlativizado omitidísimo] elecciones generales, voy a calcular de la manera más naíf que se me ocurre el índice de Shapley de los distintos partidos. Que es:

Al menos, de acuerdo con el siguiente código:

library(GameTheoryAllocation)

partidos <- c(123, 66, 57, 35, 24, 15, 7, 7,
              6, 4, 2, 2, 1, 1)
names(partidos) <- c("psoe", "pp", "cs", "iu",
                      "vox", "erc", "epc", "ciu",
                      "pnv", "hb", "cc", "na",
                      "compr", "prc")

coaliciones <- coalitions(length(partidos))
tmp <- coaliciones$Binary

profit <- tmp %*% partidos
profit <- 1 * (profit > 175)

res <- Shapley_value(profit, game = "profit")

res <- as.vector(res)
names(res) <- names(partidos)
res <- rev(res)

dotchart(res, labels = names(res),
          main = "naive shapley index \n elecciones 2019")

Lo del índice de Shapley, de ignorarlo, lo tendréis que consultar por vuestra cuenta. Al menos, para saber por qué no debería usarse tan frecuentemente (en problemas de atribución, entre otros).

Fueron mis modelos favoritos un tiempo, cuando modelaba visitas y revisitas de usuarios a cierto malhadado portal.

Si las visitas fuesen aleatorias (en cierto sentido), tendrían un aspecto no muy distinto del que se obtiene haciendo

library(IHSEP)

suppressWarnings(set.seed(exp(pi * complex(imaginary = 1))))

tms <- simPois(int = function(x) .1, cens = 1000)
hist(tms, breaks = 100, main = "Proceso homogéneo de Poisson",
      xlab = "", ylab = "frecuencia")

Es decir,

o bien una distribución uniforme en el tiempo. Pero bien puede ocurrir que una visita incremente la probabilidad de otra inmediatamente después, por lo que las visitas tenderían a arracimarse en determinados momentos. Con el paquete [IHSEP](https://cran.r-project.org/package=IHSEP) de R pueden simularse (y ajustarse) este tipo de modelos. Por ejemplo,

suppressWarnings(set.seed(exp(pi * complex(imaginary = 1))))
runif(1)
#[1] 0.4866672
set.seed(-1)
runif(1)
#[1] 0.4866672

Coda: ¿De qué si no creéis que iba esto?

Uno de mis proyectos permanentemente pospuestos es el del análisis de textos muy cortos. Se citarán Twitter y similares, aunque el € está en otros sitios, como los mensajes asociados a transferencias bancarias, reseñas o keywords.

Pero parece que no soy el único interesado en el tema. Otros con más tiempo y talento han desarrollado BTM, que parece ser una versión modificada de LDA para el análisis de textos cortos.

El artículo en el que está basado el paquete también es una buena referencia de técnicas y trucos cuando toca analizar este tipo de conjuntos de datos.