Me sorprende ver todavía a gente utilizar técnicas stepwise para la selección de variables en modelos. Sobre todo, existiendo herramientas como elastic net o lasso.
Otra de las técnicas disponibles es la del spike and slab (de la que oí hablar, recuerdo, por primera vez en el artículo de Varian Big Data: New Tricks for Econometrics). Es una técnica de inspiración bayesiana en cuya versión más cruda se imponen sobre las variables del modelo de regresión prioris que son una mezcla de dos distribuciones:
[Coge aire: aquí arranca una frase muy larga] Simplemente, que he creado un repositorio en GitHub para extraer información de los ficheros excel y sus muchas pestañas que componen el sistema de difusión de datos estadísticos sobre consumo de alimentos y bebidas de las familias que realiza el ministerio de como se llame ahora.
La página de ministerio es esta; el repositorio, este.
Nota: hay mucha información muy buena que merece ser más conocida y mejor explotada.
Estoy aquí analizando datos para un cliente interesado en estudiar si como consecuencia de uno de esos impuestos modennos con los que las administraciones nos quieren hacer más sanos y robustos. En concreto, le he echado un vistazo a si el impuesto ha encarecido el precio de los productos gravados (sí) y si ha disminuido su demanda (no) usando CausalImpact y me ha complacido mucho que la salida de summary(model, "report") sea, literalmente, esta:
Es un asunto tangencial que, además, se soluciona las más de las veces con density. Pero parece que tiene mucha más ciencia detrás.
Por algún motivo, acabé un día en la página del paquete logspline, que ajusta densidades usando splines. Su promesa es que puede realizar ajustes de densidades tan finos como
que está extraído de Polynomial Splines and their Tensor Products in Extended Linear Modeling, el artículo que le sirve de base teórica.
[Nota: el código relevante sigue estando en GitHub. No es EL código sino UN código que sugiere todos los cambios que se te puedan ocurrir. Entre otras cosas, ilustra cómo de dependientes son los resultados de la formulación del modelo, cosa muchas veces obviada.]
Continúo con la entrada de ayer, que contenía más errores que información útil respecto a objetivos y métodos.
Los objetivos del análisis son los de obtener una estimación del número de casos activos de coronavirus en la provincia de Madrid.
[Nota: si no sabes interpretar las hipótesis embebidas en el código que publico, que operan como enormes caveats, no hagas caso en absoluto a los resultados. He publicado esto para ver si otros que saben más que yo lo pulen y consiguen un modelo más razonable usándolo tal vez, ojalá, como núcleo.]
[Edición: He subido el código a GitHub.]
[El código de esta sección y los resultados contienen errores de bulto; consúltese el código de GitHub.
El de los k-vecinos es uno de mis métodos favoritos de modelización. Al menos, teóricamente: luego, en la práctica, es complicado construir una función de distancias decente. Pero tiene la ventaja indiscutible de ser tremendamente local: las predicciones para una observación concreta dependen únicamente de su entorno.
lme4::lmer (y sus derivados) es ya casi la lente a través de la que imagino cómo operan las variables dentro de un modelo. Desafortunadamente, es un modelo global y no gestiona particularmente bien las interacciones, cuando son muchas y complejas.
Esta entrada es casi una referencia para mí. Cada vez tiro más de lme4 en mis modelos y en uno en concreto que tengo entre manos toca simular escenarios. Para lo cual, simulate.merMod.
Véamoslo en funcionamiento. Primero, datos (ANOVA-style) y el modelo que piden a gritos:
library(plyr) library(lme4) a <- c(0,0,0, -1, -1, 1, 1, -2, 2) factors <- letters[1:length(a)] datos <- ldply(1:100, function(i){ data.frame(x = factors, y = a + rnorm(length(a))) }) modelo <- lmer(y ~ (1 | x), data = datos) El resumen del modelo está niquelado:
Desafortunadamente, el concepto de interacción, muy habitual en modelización estadística, no ha penetrado la literatura del llamado ML. Esencialmente, el concepto de interacción recoge el hecho de que un fenómeno puede tener un efecto distinto en subpoblaciones distintas que se identifican por un nivel en una variable categórica.
El modelo lineal clásico,
$$ y \sim x_1 + x_2 + \dots$$
no tiene en cuenta las interacciones (aunque extensiones suyas, sí, por supuesto).
El código usado en
Coronavirus: los nuevos casos diarios se estabilizan en muchos países menos en... pic.twitter.com/XOwxyccsZG
— Carlos Gil Bellosta (@gilbellosta) March 10, 2020 es
library(reshape2) library(ggplot2) library(plyr) url <- "https://raw.githubusercontent.com/CSSEGISandData/COVID-19/master/csse_covid_19_data/csse_covid_19_time_series/time_series_19-covid-Confirmed.csv" cvirus <- read.table(url, sep = ",", header = T) cvirus$Lat <- cvirus$Long <- NULL cvirus$Province.State <- NULL cvirus <- melt(cvirus, id.vars = "Country.Region") colnames(cvirus) <- c("país", "fecha", "casos") cvirus$fecha <- as.Date(as.character(cvirus$fecha), format = "X%m.%d.%y") tmp <- cvirus[cvirus$país %in% c("Italy", "Spain", "France", "Germany", "South Korea", "UK"),] foo <- function(x){ x <- x[order(x$fecha),] data.
Por referencia y afán de completar dos entradas que hice hace un tiempo sobre el método delta, esta y esta, dejo constar mención al paquete propagate, que contiene métodos para la propagación de la incertidumbre.
Para desavisados: si $latex x \sim N(5,1)$ e $latex y \sim N(10,1)$, ¿cómo sería la distribución de $latex x/y$? Etc.
Mi código (guarrongo) para seguir la evolución del coronavirus por país en cuasi-tiempo real:
library(reshape2) library(ggplot2) url <- "https://raw.githubusercontent.com/CSSEGISandData/COVID-19/master/csse_covid_19_data/csse_covid_19_time_series/time_series_19-covid-Confirmed.csv" cvirus <- read.table(url, sep = ",", header = T) cvirus$Lat <- cvirus$Long <- NULL cvirus$Province.State <- NULL cvirus <- melt(cvirus, id.vars = "Country.Region") colnames(cvirus) <- c("país", "fecha", "casos") cvirus <- cvirus[cvirus$país %in% c("Italy", "Spain"),] cvirus$fecha <- as.Date(as.character(cvirus$fecha), format = "X%m.%d.%y") ggplot(cvirus, aes(x = fecha, y = casos, col = país)) + geom_line() tmp <- cvirus tmp$fecha[tmp$país == "Spain"] <- tmp$fecha[tmp$país == "Spain"] - 9 ggplot(tmp, aes(x = fecha, y = casos, col = país)) + geom_line() tmp <- tmp[tmp$fecha > as.
Lo del coronavirus nos ha convertido a todos en epidemiólogos circunstanciales. Casi ninguno de vosotros tenéis acceso a los datos necesarios para hacer cosas por vuestra cuenta, pero sí, tal vez gracias a esta entrada, las herramientas necesarias para ello.
Podéis empezar por el paquete survellance de R, que implementa muchos de los métodos más modernos para la monitorización de brotes epidémicos.
En particular, puede que os interese la función bodaDelay, intitulada Bayesian Outbreak Detection in the Presence of Reporting Delays, y que implementa una serie de métodos para estimar el número real de casos cuando las notificaciones de los positivos llegan tarde.
