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.Date("2020-02-14"),]

ggplot(tmp, aes(x = fecha, y = log10(casos), col = país)) + geom_line()

Los datos están extraídos de aquí, por si alguien quiere reemplazar casos por defunciones o recuperados.

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. O, en plata, si dizque hay 613 confirmados oficiales, ¿cuántos podría llegar a haber realmente?

Aquí se propone un método para el análisis de datos que resume

Consta de dos procesos divergentes,

• la exploración de los datos y
• la modelización

y dos convergentes,

• la síntesis y
• la narración, que concluye el análisis.

En el enlace anterior se describe el proceso con más detalle. Eso sí, mis comentarios. El primero es que cada vez veo menos diferencia entre explorar y modelar. No entiendo ninguna exploración que no esté motivada por un modelo implícito; p.e., representar las medias por grupo no es otra cosa que una ANOVA para pobres. Crear árboles de decisión sobre los datos brutos es muy indicativo de por dónde van los tiros en los datos, qué variables son más importantes, cuáles son irrelevantes, etc. Obviamente, el modelo final no va a ser ninguno de estos protomodelos, pero sí que contienen su germen.

Aquí se recomienda, con muy buen criterio, no realizar clasificación pura, i.e., asignando etiquetas 0-1 (en casos binarios), sino proporcionar en la medida de lo posible probabilidades. Y llegado el caso, distribuciones de probabilidades, claro.

La clave es, por supuesto:

The classification rule must be reformulated if costs/utilities or sampling criteria change.

Contexto:

modelo <- lm(dist ~ speed, data = cars)

Intervalos de confianza:

head(predict(modelo, interval = "confidence"))
#        fit        lwr       upr
#1 -1.849460 -12.329543  8.630624
#2 -1.849460 -12.329543  8.630624
#3  9.947766   1.678977 18.216556
#4  9.947766   1.678977 18.216556
#5 13.880175   6.307527 21.452823
#6 17.812584  10.905120 24.720047

Intervalos de predicción:

head(predict(modelo, interval = "prediction"))
#        fit       lwr      upr
#1 -1.849460 -34.49984 30.80092
#2 -1.849460 -34.49984 30.80092
#3  9.947766 -22.06142 41.95696
#4  9.947766 -22.06142 41.95696
#5 13.880175 -17.95629 45.71664
#6 17.812584 -13.87225 49.49741

Creo que la diferencia (y el significado) es claro. Para todos los demás, esto.

Así arranca este artículo, que presenta una extensión de XGBoost para predicciones probabilísticas. Es decir, un paquete que promete no solo una estimación del valor central de la predicción sino de su distribución.

La versión equivalente de lo anterior en el mundo de los random forests está descrito aquí, disponible aquí y mucho me temo que muy pronto voy a poder contar por aquí si está a la altura de las expectativas.

Sigo con ajustes robustos. Y cosas que como matemático, me ponen muy nervioso.

Una de las maneras de hacer ajustes robustos es la de sustituir la función cuadrática por la biweight. Es decir, utilizar la función que aparece la derecha en

en lugar de la de la izquierda. O, dicho de otra manera, en lugar de tratar de minimizar

$$ \sum_i \rho(y_i - f_\alpha(x_i))$$

usando $latex \rho(x) = x^2$, que es la función que se representa a la izquierda y a la que estamos acostumbrados, usar la de la derecha. Que es la función biweight de Tukey.

Los (o ciertos) usuarios de R de Galicia están organizando una conferencia alrededor del mundo R de la mano de satRdays. Serán el sábado 12 de septiembre (de 2020) y los interesados en saber más al respecto, harán bien en visitar esta página.

De todos modos, si quieres presentar una charla o taller, el plazo límite parece ser el día 15 de abril.

… que la de “Algoritmos” y acatarrantes definiciones de “justicia”. Que es casi una versión de la anterior reduciendo la varianza de las betas.

Las dos poblaciones de interés tienen una tasa de probabilidad (o de riesgo, en la terminología del artículo original) de .4 y .6 respectivamente. Aproximadamente el 40% de los primeros y el 60% de los segundos tienen y = 1.

El modelo (el algoritmo) es perfecto y asigna a los integrantes del primer grupo un scoring de .4 y a los del segundo, de .6.