Casos de coronavirus en Madrid provincia: un modelo muy crudo basado en la mortalidad

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

Escribo casi al vuelo e inspirado por la idea de que:

  • Los casos verificados de coronavirus no son fiables. Ni están todos los que son, ni son todos los que están.
  • Los casos de defunciones son mucho más fiables: son menos y se les presta mucha atención.

La idea la he expresado públicamente en:

Y aquí van los resultados crudísimos de un modelo crudísimo y seguramente erróneo cuya convergencia ni he constrastado ni nada:

Y el código es:

data {
  int<lower=0> N;
  int<lower=0> dia0;
  vector[N] dia;
  vector[N] defs;
}

parameters {
  real<lower = 0, upper = 1000> casos_0; 
  vector[N] contagios;
  real<lower = 0, upper = 3> r0;  //priori cutre sobre r0, infectados diarios que genera un caso
  real<lower = 0, upper = .02> letalidad;
}

model {
  vector[N] casos;
  
  contagios[1] ~ normal(casos_0 * r0, 1);
  casos[1] = casos_0 + contagios[1];
  
  
  for (i in 2:N){
    contagios[i] ~ normal(casos[i-1] * r0, 1);
    casos[i] = casos[i-1] + contagios[i];
  }
  
  for (i in dia0:N) {
    defs[i] ~ normal(letalidad * casos[i - 22], 1);
  }
}

y

library(rstan)
library(reshape2)
library(plyr)

library(ggplot2)



dat <- read.csv("https://raw.githubusercontent.com/datadista/datasets/master/COVID%2019/ccaa_covid19_fallecidos.csv")
dat <- melt(dat, id.vars = c("cod_ine", "CCAA"))
dat$fecha <- as.Date(dat$variable, format = "X%d.%m.%Y")

fecha_ini <- min(dat$fecha)

madrid <- dat[dat$CCAA == "Madrid",]
madrid$dia <- as.numeric(madrid$fecha - fecha_ini)

madrid <- madrid[, c("dia", "value")]
colnames(madrid) <- c("dia", "defs")

tmp <- data.frame(dia = -30:-1, defs = 0)

madrid <- rbind(tmp, madrid)


fit <- stan(file = "stan.stan",
            data = list(N = nrow(madrid), dia0 = which(madrid$dia == 0), dia = madrid$dia, defs = madrid$defs), 
            iter = 10000, warmup = 2000, 
            chains = 1, thin = 10)

res <- as.data.frame(fit)

contagios <- extract(fit, pars = "contagios")$contagios

contagios <- data.frame(t(contagios))
contagios$fecha <- as.Date(madrid$dia, origin = fecha_ini)

contagios <- melt(contagios, id.vars = "fecha")

casos <- ddply(contagios, .(variable), transform, casos = cumsum(value))


ggplot(casos, aes(x = fecha, y = casos, group = variable)) +
    geom_line(alpha = 0.3) + 
    xlab("fecha") + ylab("casos") +
    ggtitle("Casos de coronavirus en Madrid\n(¡Resultado de un modelo muy crudo y\n casi seguro con errores!)")

5 comentarios sobre “Casos de coronavirus en Madrid provincia: un modelo muy crudo basado en la mortalidad

  1. Emilio 19 marzo, 2020 23:31

    Haciendo la cuenta de la vieja con los resultados de tu modelo: si en Madrid hay 6.6 millones de personas (con un 28% de la población mayor de 55 años), sale un tasa de mortalidad del 9%. No creo que ese modelo sea válido en nuestra civilización actual. Puede que valga para otros lugares con sistemas de salud más limitados o que no tomen medidas activas.

    Desde mi ignorancia -no tengo ni idea de este tema, y la ignorancia es atrevida. -, creo que hay dos parámetros del modelo que habría que repensar:
    – la r0 real debería disminuir a medida que más personas estén infectadas (aunque parece que sí se pueden volver a infectar) y aisladas.
    – la letalidad debería ir disminuyendo. Los más vulnerables fallecen antes.

    Gracias por el esfuerzo divulgativo que estás desarrollando. Mejor pecar de exceso que por defecto.

  2. Emilio 19 marzo, 2020 23:42

    Acabo de ver en tu twitter que señalas que este modelo solo es valido para la fase de inicio (exponencial). Creo que mis comentarios no se aplican a esta fase.

  3. Isidro Hidalgo 20 marzo, 2020 9:34

    Efectivamente, Emilio, es un modelo para la fase inicialmente exponencial de la epidemia. Posteriormente el virus se va encontrando con gente infectada, por lo que la curva exponencial no crece tan rápido. Para prolongar un poco más tiempo la validez del modelo puede utilizarse un exponente para calcular los casos del día anterior: casos[i] = r0 * casos[i-1] ^ r1, con r1 entre 0 y 1. De todas formas, cuando avanza la enfermedad, este modelo tampoco sirve, porque hay que tener en cuenta de alguna forma la población que queda sana…

Comenta

Your email address will not be published.

Puedes usar estas etiquetas y atributos de HTML: <a href="" title=""> <abbr title=""> <acronym title=""> <b> <blockquote cite=""> <cite> <code> <del datetime=""> <em> <i> <q cite=""> <s> <strike> <strong>

Este sitio usa Akismet para reducir el spam. Aprende cómo se procesan los datos de tus comentarios.