[En esta entrada deambulo peligrosamente por los límites de un NDA; sin embargo, me siento obligado a exponerme a las posibles consecuencias debido a la gravedad de las circunstancias actuales.]

En España existe un mecanismo de monitorización de la mortalidad diaria por todas las causas. Su existencia no es explícitamente pública, pero sí que existen indicios implícitos de su existencia en informes de salud pública: véanse, p.e., referencias a MoMo y EuroMOMO aquí. [Nota: MoMo es el acrónimo de mortality monitoring].

[Retomando un tema que dejé inconcluso y que tampoco remataré hoy aquí.]

Imagina que quieres saber cuánto le dura a la gente el portátil. Para eso preguntas por ahí cuándo se compraron el último.

Lo que obtienes es un conjunto de datos donde todas las observaciones están censuradas. Y no, el análisis de la supervivencia clásico no funciona.

Buscando en la literatura he encontrado, sin embargo, Survival Analysis of Backward Recurrence Times, donde se discute el problema y al que, bueno, otro día con menos penas volveré.

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.frame(fecha = x$fecha[-1],
        casos = diff(x$casos))
}

res <- ddply(tmp, .(país), foo)

res$país <- reorder(res$país, res$casos, function(x) -max(x))

res <- res[res$fecha > as.Date("2020-02-15"),]

ggplot(res, aes(x = fecha, y = casos)) +
    geom_point(size = 0.5) + geom_line(alpha = 0.3) +
    facet_wrap(~país, scales = "free_y") +
    ggtitle("Coronavirus: new daily cases") +
    theme_bw()

ggsave("/tmp/new_daily_cases.png", width = 12,
    height = 8, units = "cm")

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.