R

Modelos de conteos con sobredispersión (con Stan)

R

Esta entrada muestra cómo afrontar (con Stan) un problema que encontré el otro día en un lugar que no puedo mencionar pero en el que sé que me leen (y los destinatarios sabrán que va por ellos).

El contexto es el siguiente: se hace un test A/B donde la variable de interés son unos conteos. Hay varios grupos (aquí los reduciré a dos) y los datos siguen aproximadamente (aquí omitiré la parte de la inflación de ceros) una distribución de Poisson. Pero solo aproximadamente: existe sobredispersión, es decir, la varianza de los datos excede su media.

data.tree: porque no todos los datos son tabulares

R

De acuerdo, casi todos los datos son tabulares. Digamos que el 90% de ellos. Pero muchos de ellos, no. Y data.tree es un paquete con muy buena pinta para manejar estructuras arborescentes de datos: véanse esta y esta viñeta.

Como no podía ser de otra manera, tiene funciones para recorrer, filtrar y podar los árboles de datos.

La aplicación gracias a la cual di con él es el paquete prof.tree, que es lo mismo que el Rprof de toda la vida… solo que mola más:

Siete años después, dejo la presidencia de la Comunidad R Hispano

R

Eso, que dejo la comunidad de la Comunidad R Hispano. Ocho años después, que ya son. La noticia, en todo caso, no es tanto que abandone la presidencia sino las circunstancias que me condujeron a ella. Noticias viejas, pero noticias al fin y al cabo, que sirven para entender por qué lo fui entonces y por qué dejo de serlo ahora.

La Comunidad R Hispano (por qué se llama así y no, como habría sido natural, Asociación Española de Usuarios de R, es una larga historia que tal vez cuente algún día) se fundó en Madrid hace ocho años en el seno de las III Jornadas de Usuarios de R (a todo esto, esa página la hice yo en html puro y con vi como editor), cuando la comunidad (informal) de usuarios de R ya llevaba un tiempo desarrollando actividades.

Colinealidad y posterioris

Estadística, R

En esta entrada voy a crear un conjunto de datos donde dos variables tienen una correlación muy alta, ajustar un modelo de regresión y obtener la siguiente representación de la distribución a posteriori de los coeficientes,

donde se aprecia el efecto de la correlación entre x1 y x2.

El código,

library(mvtnorm)
library(rstan)
library(psych)

n <- 100
corr_coef <- .9

x <- rmvnorm(n, c(0, 0),
  sigma = matrix(c(1, corr_coef, corr_coef, 1), 2, 2))
plot(x)

x1 <- x[,1]
x2 <- x[,2]
x3 <- runif(n) - 0.5

y <- 1 + .4 * x1 - .2 * x2 + .1 * x3 + rnorm(n, 0, .1)

summary(lm(y ~ x1 + x2 + x3))

stan_code <- "
data {
  int N;
  vector[N] y;
  vector[N] x1;
  vector[N] x2;
  vector[N] x3;
}
parameters {
  real a;
  real a1;
  real a2;
  real a3;
  real sigma;
}

model {
  a ~ cauchy(0,10);
  a1 ~ cauchy(0,2.5);
  a2 ~ cauchy(0,2.5);
  a3 ~ cauchy(0,2.5);

  y ~ normal(a + a1 * x1 + a2 * x2 + a3 * x3, sigma);
}"


datos_stan <- list(
    N = n,
    y = y,
    x1 = x1,
    x2 = x2,
    x3 = x3
)

fit2 <- stan(model_code = stan_code,
              data = datos_stan,
              iter = 10000, warmup = 2000,
              chains = 2, thin = 4)

res <- as.data.frame(fit2)
pairs.panels(res[, c("a", "a1", "a2", "a3", "sigma")])

Cuatro paquetes interesantes de R

R

Son paquetes que marcado como potencialmente relevantes pero que aún no he revisado como debiera. Tal vez alguien tenga algo más que decir sobre ellos. Tiene los comentarios, por supuesto, abiertos.

longRPart2: Particionamiento recursivo para modelos longitudinales. Extiende ctree y, por supuesto, mob del paquete party a datos de tipo longitudinal.

radiant: Más que un paquete, es un conjunto de paquetes para business analytics usando R y Shiny. Ni idea de para qué parte de ese amplio campo del business analytics puede resultar útil, pero si resulta que es precisamente el tuyo, ¡enhorabuena!

Planes de búsqueda y rescate con R

Estadística, R

Existe un paquete muy curioso en CRAN, rSARP para diseñar, optimizar y comunicar la evolución de planes de búsqueda y/o rescate (p.e., de un niño desaparecido en un monte).

Es particularmente interesante porque este tipo de problemas lo tienen todo: desde distribuciones a priori (sobre dónde es más probable encontrar lo que se busca) hasta la decisión final (explórese tanto aquí y tanto allá) teniendo en cuenta restricciones de tiempo y recursos.

Disponible el fichero de datos abiertos más goloso de ambas castillas: las rutas de Bicimad

R

Albricias, el ayuntamiento de Madrid ha liberado el fichero más goloso de ambas castillas: el de las rutas de usuarios de Bicimad, viaje a viaje, con su estación de origen, estación de destino, tiempo de recorrido, etc. Tiempo os falta para echarle un vistazo y hacer cosas chulas con él.

Los datos están aquí.

Se puede leer con código no muy distinto de este:

library(RJSONIO)

raw <- readLines("201808_Usage_Bicimad.json")
dat <- iconv(raw, "latin1", "utf8")
dat <- sapply(dat, fromJSON)

A bote pronto, se me ocurren algunas cosas que se pueden hacer con esos datos:

Los datos están histogramizados... ¿quién los deshisotogramizará?

Números, R

Hace un tiempo quise hacer cosas malísimas con datos fiscales de España y Dinamarca. Pero los datos estaban histogramizados:

Gracias a Freakonometrics di con binequality. Adaptando su código, escribo

library(rvest)
library(plyr)

dk <- read_html("http://www.skm.dk/english/facts-and-figures/progression-in-the-income-tax-system")
tmp <- html_nodes(dk, "table")
tmp <- html_table(tmp[[2]])

header <- tmp[1,]
tmp <- tmp[-c(1, 2),]
colnames(tmp) <- header

# elimino declaraciones negativas
tmp <- tmp[-1,]

# elimino el total
tmp <- tmp[-(nrow(tmp)),]

colnames(tmp) <- c("rango", "contribuyentes",
    "X1", "income", "tax1", "tax2", "pct")

irpf_dk <- tmp[, c("rango", "contribuyentes",
    "income", "tax1", "tax2")]

irpf_dk$contribuyentes <- as.numeric(irpf_dk$contribuyentes)
irpf_dk$income <- as.numeric(irpf_dk$income)
irpf_dk$tax1 <- as.numeric(irpf_dk$tax1)
irpf_dk$tax2 <- as.numeric(irpf_dk$tax2)

irpf_dk$tax <- irpf_dk$tax1 + irpf_dk$tax2
irpf_dk$tax1 <- irpf_dk$tax2 <- NULL
irpf_dk$pct <- irpf_dk$tax / irpf_dk$income


irpf_dk$desde <- c(0, 25, 50, 75, 100, 125, 150,
    200, 250, 300, 350, 400, 500, 750, 1000)
irpf_dk$hasta <- c(irpf_dk$desde[-1], Inf)

irpf_dk$desde <- irpf_dk$desde / 7.44
irpf_dk$hasta <- irpf_dk$hasta / 7.44
irpf_dk$income <- irpf_dk$income / 7.44
irpf_dk$tax    <- irpf_dk$tax / 7.44

irpf_dk$mean_income <- irpf_dk$income /
        irpf_dk$contribuyentes * 1000

irpf_dk$rango <- NULL

para bajar y preprocesar los datos y después

Contraargumentando (materialmente) sobre la falacia del fiscal

R

Hace un par de días hablé de la falacia del fiscal y granos de arroz. La entrada iba acompañada de

y la lección era: es raro no encontrar ningún clúster cuando se tiran al azar granos de arroz sobre una superficie. De lo que se derivaban más cosas que es ocioso repetir aquí.

Pero el gráfico no es desconocido para los viejos del lugar: se parece mucho al de la página 319 de ESL. Para los que no lo tengáis a mano, la parte donde se habla de un algoritmo que se llama igual que un general de Reus con calle en Méjico DF: PRIM.

¿Por que slt-ear si puedes stR-ear?

Estadística, R

La función stl (véase aquí un ejemplo de uso). Pero tiene sus limitaciones.

El paquete stR la extiende y permite, entre otras cosas, introducir distintos tipos de estacionalidades (p.e., anuales y semanales).