covid_madrid/src.R.deprecated at master · cjgb/covid_madrid · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
library(rstan)
library(reshape2)
library(plyr)
library(ggplot2)

pop_madrid <- 6.6e6

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

fecha_ini <- min(defs$fecha)

madrid <- defs[defs$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)
defunciones_diarias <- diff(madrid$defs)

madrid <- data.frame(dia = madrid$dia[-1], defs = defunciones_diarias)


fit <- stan(file = "stan.stan",
            data = list(N = nrow(madrid),
                        dia0 = which(madrid$dia == 0),
                        pop = pop_madrid,
                        dia = madrid$dia,
                        defs = madrid$defs),
            iter = 20000, warmup = 5000,
            chains = 4, thin = 10, cores = 4,
            control = list(adapt_delta = 0.99))

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.01) +
    xlab("fecha") + ylab("casos") +
    ggtitle("Casos de coronavirus en Madrid\n(¡Resultado de un modelo muy crudo y\n casi seguro con errores!)")

ggsave("Rplot.png")

tmp <- casos[casos$fecha == max(casos$fecha),]
tmp$casos <- tmp$casos / 1000
ggplot(tmp, aes(x = casos)) +
    geom_histogram(fill = "steelblue") +
    xlab("casos (miles)") +
    ylab("número de simulaciones") +
    ggtitle(paste0("casos estimados en ", Sys.Date()))

ggsave("Rplot01.png")


# just a check

# ggplot(tmp, aes(x = log10(casos))) +
#     geom_histogram(fill = "steelblue") +
#     xlab("casos (miles)") +
#     ylab("número de simulaciones") +
#     ggtitle(paste0("casos estimados en ", Sys.Date()))

# hist(tmp$casos, breaks = 30, main = "Casos 'hoy'", col = "steelblue", xlab = "casos", freq = FALSE)


tmp <- res[, c("casos_0", "r0", "letalidad")]
plot(tmp)


traceplot(fit, pars = c("casos_0", "r0", "letalidad"))

# evaluación del ajuste
pairs(fit, pars = "contagios", include = FALSE)