covid-19projections/run-sensitivity.R at master · phebo/covid-19projections · 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
# Copyright (C) 2020, Phebo Wibbens, Wesley Koo, and Anita McGahan

#   This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.

# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.

# You should have received a copy of the GNU General Public License
# along with this program.  If not, see <https://www.gnu.org/licenses/>.


library(stats4)
library(xtable)
library(tidyverse)
library(furrr)
library(rstan)
source("functions.R")

plan(multiprocess) # Allow parallel computing
time.now <- format(Sys.time(), format='%y%m%d-%H%M%S')
print(time.now)
suppressWarnings(dir.create(file.path("output")))

dfJh <- read_csv("input/jh-database.csv")
dfEcon <- read_csv("input/econ-database.csv")
dfPop <- read_csv("input/econ-population.csv")
dfOx <- read_csv("input/oxford-policy.csv")
dfHol <- read_csv("input/holidays.csv") %>% select(-source)

l <- list(
  clean.data(dfJh, dfEcon, dfPop, dfOx, dfHol),
  clean.data(dfJh, dfEcon, dfPop, dfOx, dfHol, minPop = 1e7),
  clean.data(dfJh, dfEcon, dfPop, dfOx, dfHol, minPop = 3e6),
  clean.data(dfJh, dfEcon, dfPop, dfOx, dfHol, pOutl = 1e-2),
  clean.data(dfJh, dfEcon, dfPop, dfOx, dfHol, pOutl = 1e-4),
  clean.data(dfJh, dfEcon, dfPop, dfOx, dfHol, idgSig = 0.01),
  clean.data(dfJh, dfEcon, dfPop, dfOx, dfHol, idgSig = 0.05),
  clean.data(dfJh, dfEcon, dfPop, dfOx, dfHol, geoExclude = c("Chile", "South Africa", "Brazil", "Mexico")),
  # clean.data(dfJh, dfEcon, dfPop, dfOx, dfHol, dgSig = c(1,1e-4)),
  clean.data(dfJh, dfEcon, dfPop, dfOx, dfHol, idgLam = c(0.85, 0.8)),
  clean.data(dfJh, dfEcon, dfPop, dfOx, dfHol, dgMin = -0.05)
)
specNames <- c("Base", "Population > 10M", "Population > 3M", "P(outlier) = 0.01", "P(outlier) = 0.0001",
               "Stdev = 0.01", "Stdev = 0.05", "Excl non-EU/US", "Lower lambda", "Negative Delta-g")

m <- stan_model("model.stan")
specs <- map(l, ~ list(m = m, data = .$lData, pars = "dg", iter = 700, warmup = 500, chains = 2, thin = 2))
chains <- make.chains(specs)
fits.chain <- future_map(chains, ~ do.chain(.))
fits <- cons.fits(fits.chain, chains)
save(list = ls(), file = paste0("output/image-sens-", time.now, ".RData"))
dgs <- map(fits, ~ extract(.)$dg)
names(dgs) <- specNames
dfRaw <- map2_dfr(dgs, l, ~ expand_grid(pol = .y$p$vPol, iter = 1:nrow(.x)) %>% mutate(value = as.vector(.x)), .id = "spec")
df <- dfRaw %>% group_by(spec, pol) %>%
  summarize(estimate = median(value), low = quantile(value, probs=0.025), high = quantile(value, probs=0.975)) %>% ungroup() %>%
  mutate(spec = factor(spec, levels = specNames))
df <- df %>% left_join(df %>% filter(spec == "Base") %>% select(pol, baseEst = estimate)) %>%
  mutate(outside = baseEst < low | baseEst > high)
fSens <- df %>% ggplot(aes(x = fct_rev(spec), y = estimate, ymin = low, ymax = high, color = outside)) + geom_pointrange() +
  geom_hline(yintercept = 0, linetype = 2) +
  facet_wrap(~ pol, ncol = 4, labeller = label_wrap_gen(25)) + coord_flip() +
  xlab(element_blank()) + ylab(element_blank()) + scale_color_manual(values = c("black", "red"), guide = FALSE) +
  theme(axis.text.y = element_text(hjust=0), strip.text = element_text(size = 8))

ggsave(paste0("figures/fig-sens.png"), height = 9, width = 6.5)
ggsave(paste0("output/chart-sens-", time.now, ".pdf"), height=12, width=8)