[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:

summary(modelo)

# Linear mixed model fit by REML ['lmerMod']
# Formula: y ~ (1 | x)
# Data: datos
#
# REML criterion at convergence: 2560.3
#
# Scaled residuals:
#     Min      1Q  Median      3Q     Max
# -3.6798 -0.6442 -0.0288  0.6446  3.3582
#
# Random effects:
#     Groups   Name        Variance Std.Dev.
# x        (Intercept) 1.5197   1.2328
# Residual             0.9582   0.9789
# Number of obs: 900, groups:  x, 9
#
# Fixed effects:
#     Estimate Std. Error t value
# (Intercept) -0.009334   0.412212  -0.023

En particular,

El primero es el chino. Es el que se aplicó a otras crisis víricas (SARS, etc.), a la viruela y a la polio. Consiste en aplicar medidas drásticas hasta que el virus desaparezca. De hecho, hay provincias en china que llegaron a tener un número importante de casos,

pero donde ya no quedan casos activos:

El otro es el escenario RU: el virus va a seguir entre nosotros y todos, en algún momento u otro vamos a pasar por él (o, más propiamente, a la inversa). En cuyo caso:

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).

[Este es un aviso para todos aquellos que depositan una excesiva fe en lo que nos cuenta el INE.]

La causa de muerte no es la causa de muerte. Al menos, necesariamente. Lo que el INE llama causa de muerte es una imagen distorsionada de la causa de muerte por culpa de un embudo administrativo.

Comiendo con unos epidemiólogos en el ISCIII hace un tiempo, me decían, con cierta envidia, cómo en otros países como Dinamarca, se registraban hasta ocho causas de muerte: la última, la concomitante, la… Y bromeaban diciendo que, al final, todos nos morimos de parada cardiorrespiratoria.

Este buscador de oro busca pepitas en su tramo de río. El río arrastra piedrecitas, muchas piedrecitas, y pepitas de oro, pocas pepitas de oro.

Tiene un artilugio que toma barro del río y que hace lo siguiente:

• Descarta casi todas las piedrecitas (y el resto las mete en una caja)
• Detecta casi todas las pepitas (y las mete en la misma caja)

Al final del día, ¿qué encontrará en la caja?

Existe una respuesta naive a la pregunta anterior que es la que tal vez en la que estéis pensando.

Pero, ¡ah!, existe también otra ciencia, la política, que no lo es menos que otras, en las que la evidencia dice cosas tales como:

Our results show that voters significantly reward disaster relief spending, holding the incumbent presidential party accountable for actions taken after a disaster. In contrast, voters show no response at all, on average, to preparedness spending, even though investing in preparedness produces a large social benefit.

[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.