Skip to content

Commit 6eb44cd

Browse files
nicholasdaviesnicholasdavies
committed
Revisions
1 parent ea01647 commit 6eb44cd

File tree

5 files changed

+218
-31
lines changed

5 files changed

+218
-31
lines changed

‎.gitignore

+6-2
Original file line numberDiff line numberDiff line change
@@ -2,14 +2,18 @@ Icon?
22

33
.DS_Store
44

5-
*.qs
6-
75
*.pdf
86

97
*.png
108

119
*.csv
1210

11+
*.o
12+
1313
output/
1414

1515
misc/
16+
17+
submission/
18+
19+
covidm/build/*

‎UK-view.R

+157-6
Original file line numberDiff line numberDiff line change
@@ -2,6 +2,7 @@
22
# UK model: plot outputs
33
# - - - - - - - - - - - - - - - - - - - - - - -
44

5+
library(ggplot2)
56
library(cowplot)
67
library(stringr)
78
library(rlang)
@@ -66,8 +67,7 @@ deaths, peak, week
6667
beds_icu, peak, t
6768
beds_nonicu, peak, t
6869
cases, peak_time, week
69-
trace_lockdown, lockdown_duration, t
70-
subclinical, total, t")
70+
trace_lockdown, lockdown_duration, t")
7171

7272
median_ci = function(x, conf = 0.95)
7373
{
@@ -77,7 +77,7 @@ median_ci = function(x, conf = 0.95)
7777
return (y)
7878
}
7979

80-
make_table = function(d)
80+
make_table = function(d, table_spec = table_spec)
8181
{
8282
d[, week := t %/% 7]
8383
results = NULL
@@ -93,6 +93,18 @@ make_table = function(d)
9393
res = res[, median_ci(x), by = .(scenario, region, statistic)];
9494
stat_nice = paste("Total", comp);
9595
} else if (stat == "peak") {
96+
# res = d[region == "United Kingdom" & compartment == comp, .(x = sum(value)), by = c("scenario", "run", "t", "region")];
97+
# if (time == "week") {
98+
# res = merge(res, res[, .(peak_day = t[which.max(x)]), by = .(scenario, run, region)]);
99+
# res[, in_peak_week := abs(t - peak_day) <= 3];
100+
# res = res[in_peak_week == T, .(x = sum(x)), by = .(scenario, run, region)];
101+
# } else {
102+
# res = res[, .(x = max(x)), by = .(scenario, run, region)];
103+
# }
104+
# res[, statistic := paste(stat, comp)];
105+
# res = res[, median_ci(x), by = .(scenario, region, statistic)];
106+
# stat_nice = ifelse(time == "t", paste("Peak", comp, "required"),
107+
# paste(comp, "in peak week"));
96108
res = d[region == "United Kingdom" & compartment == comp, .(x = sum(value)), by = c("scenario", "run", time, "region")];
97109
res = res[, .(x = max(x)), by = .(scenario, run, region)];
98110
res[, statistic := paste(stat, comp)];
@@ -241,6 +253,10 @@ plot_epi = function(d0, t, quant, ymd_start, ymd_truncate = "2050-01-01", colour
241253
qrun = t[scenario == "Base", sum(total), by = run][, which(rank(V1) %in% round(1 + (.N - 1) * quant))];
242254
d = d[run %in% qrun];
243255

256+
if (d[, max(run) == 5]) {
257+
mrun = 5;
258+
}
259+
244260
# Merge intervention traces
245261
trace_school = d[compartment == "trace_school", .(run, t, trace_school = value - 1, scenario)]
246262
d = merge(d, trace_school, by = c("run", "t", "scenario"), all.x = T);
@@ -297,12 +313,15 @@ plot_epi = function(d0, t, quant, ymd_start, ymd_truncate = "2050-01-01", colour
297313
theme_set(theme_cowplot(font_size = 7, line_size = 0.25))
298314

299315
# load data
316+
covid_uk_path = paste0(covid_uk_path, "/output/may7")
300317
d1 = reflow_dynamics(qread(paste0(covid_uk_path, "/1-dynamics.qs")));
301318
t1 = reflow_totals(qread(paste0(covid_uk_path, "/1-totals.qs")));
302319
d2.1 = reflow_dynamics(qread(paste0(covid_uk_path, "/2.1-dynamics.qs")));
303320
t2.1 = reflow_totals(qread(paste0(covid_uk_path, "/2.1-totals.qs")));
304321
d2.2 = reflow_dynamics(qread(paste0(covid_uk_path, "/2.2-dynamics.qs")));
305322
t2.2 = reflow_totals(qread(paste0(covid_uk_path, "/2.2-totals.qs")));
323+
d2.2v =reflow_dynamics(qread(paste0(covid_uk_path, "/2.2V-dynamics.qs")));
324+
t2.2v = reflow_totals(qread(paste0(covid_uk_path, "/2.2V-totals.qs")));
306325
d3 = reflow_dynamics(qread(paste0(covid_uk_path, "/3-dynamics.qs")));
307326
t3 = reflow_totals(qread(paste0(covid_uk_path, "/3-totals.qs")));
308327
d4 = reflow_dynamics(qread(paste0(covid_uk_path, "/4-dynamics.qs")));
@@ -353,6 +372,11 @@ f = plot_grid(pla1, plb, pla2, plR,
353372
nrow = 2, ncol = 2, rel_widths = c(3, 2), labels = c("a", "b", "", "c"), label_size = 9, align = "hv", axis = "bottom")
354373
ggsave(paste0(covid_uk_path, "/fig-12week.pdf"), f, width = 20, height = 12, units = "cm", useDingbats = F);
355374

375+
ggsave(paste0(covid_uk_path, "/fig-12week-a.pdf"), plot_grid(pla1, pla2, ncol = 1, align = "hv", axis = "bottom"),
376+
width = 12, height = 12, units = "cm", useDingbats = F);
377+
ggsave(paste0(covid_uk_path, "/fig-12week-b.pdf"), plb, width = 8, height = 6, units = "cm", useDingbats = F);
378+
ggsave(paste0(covid_uk_path, "/fig-12week-c.pdf"), plR, width = 8, height = 6, units = "cm", useDingbats = F);
379+
356380
# ANALYSIS 2 - TRIGGERS
357381
d2.1[scenario != "Base", scenario := paste(scenario, "national")]
358382
d2.2[scenario != "Base", scenario := paste(scenario, "local")]
@@ -377,7 +401,7 @@ save_table(tb2, paste0(covid_uk_path, "/table-triggers.csv"));
377401
pl2 = plot_attackrate(t2)
378402
pl3 = plot_epi(d2, t2, (0:10)/10, "2020-01-29")
379403
f = plot_grid(pl1, pl2, pl3, ncol = 1, rel_heights = c(6, 6, 10), labels = c("a", "b", "c"), label_size = 9);
380-
ggsave(paste0(covid_uk_path, "/COVID-UK/full-2.pdf"), f, width = 20, height = 22, units = "cm", useDingbats = F);
404+
ggsave(paste0(covid_uk_path, "/full-2.pdf"), f, width = 20, height = 22, units = "cm", useDingbats = F);
381405

382406
pla1 = plot_epi(d2[compartment != "deaths" & compartment != "beds_nonicu" & compartment != "beds_icu" &
383407
scenario %like% "Local"], t2, (0:10)/10, "2020-01-29", "2020-8-31", exclude = "Base");
@@ -407,15 +431,18 @@ plc = plot_grid(plc1, plc2, plc3, nrow = 3, align = "v", axis = "bottom", rel_he
407431
f = plot_grid(pla, plb, plc,
408432
ncol = 3, labels = c("a", "b", "c"), label_size = 9, rel_widths = c(2, 1, .5))
409433
ggsave(paste0(covid_uk_path, "/fig-triggers.pdf"), f, width = 20, height = 8, units = "cm", useDingbats = F);
434+
ggsave(paste0(covid_uk_path, "/fig-triggers-a.pdf"), pla, width = 20*4/7, height = 8, units = "cm", useDingbats = F);
435+
ggsave(paste0(covid_uk_path, "/fig-triggers-b.pdf"), plb, width = 20*2/7, height = 8, units = "cm", useDingbats = F);
436+
ggsave(paste0(covid_uk_path, "/fig-triggers-c.pdf"), plc, width = 20*1/7, height = 8, units = "cm", useDingbats = F);
410437

411438
# ANALYSIS 3 - LOCKDOWN
412439
d3[scenario == "Intensive Interventions NA lockdown", scenario := "Intensive Interventions"];
413440
d3[scenario == "Intensive Interventions 1000 lockdown", scenario := "Lockdown 1000-bed trigger"];
414441
d3[scenario == "Intensive Interventions 2000 lockdown", scenario := "Lockdown 2000-bed trigger"];
415442
d3[scenario == "Intensive Interventions 5000 lockdown", scenario := "Lockdown 5000-bed trigger"];
416-
#d3[compartment == "subclinical"]$value = d3[compartment == "subclinical"]$value + d3[compartment == "cases"]$value
443+
d3[compartment == "subclinical"]$value = d3[compartment == "subclinical"]$value + d3[compartment == "cases"]$value
417444

418-
tb3 = make_table(d3)
445+
tb3 = make_table(d3, rbind(table_spec, data.table(compartment = "subclinical", stat = "total", time = "t")))
419446
pl1 = plot_table(tb3)
420447
save_table(tb3, paste0(covid_uk_path, "/table-lockdown.csv"));
421448
pl2 = plot_attackrate(t3)
@@ -461,6 +488,11 @@ plR = ggplot(r0s2) +
461488
f = plot_grid(pla1, plb, pla2, plR,
462489
nrow = 2, ncol = 2, rel_widths = c(3, 2), labels = c("a", "b", "", "c"), label_size = 9, align = "hv", axis = "bottom")
463490
ggsave(paste0(covid_uk_path, "/fig-lockdown.pdf"), f, width = 20, height = 12, units = "cm", useDingbats = F);
491+
ggsave(paste0(covid_uk_path, "/fig-lockdown-a.pdf"), plot_grid(pla1, pla2, ncol = 1, align = "hv", axis = "bottom"),
492+
width = 12, height = 12, units = "cm", useDingbats = F);
493+
ggsave(paste0(covid_uk_path, "/fig-lockdown-b.pdf"), plb, width = 8, height = 6, units = "cm", useDingbats = F);
494+
ggsave(paste0(covid_uk_path, "/fig-lockdown-c.pdf"), plR, width = 8, height = 6, units = "cm", useDingbats = F);
495+
464496

465497
# ANALYSES 4,6 - GRANDPARENTS AND SPORTS/LEISURE
466498
# Grandparents
@@ -504,8 +536,127 @@ pla = plot_table(tb6[scenario != "Base" & statistic != "Deaths in peak week"]) +
504536
guides(colour = guide_legend(nrow = 3, byrow = TRUE)) + labs(colour = NULL)
505537
f = plot_grid(pla, plb, nrow = 2, labels = c("a", "b"), label_size = 9, align = "hv", axis = "bottom");
506538
ggsave(paste0(covid_uk_path, "/fig-misc.pdf"), f, width = 9, height = 12, units = "cm", useDingbats = F);
539+
ggsave(paste0(covid_uk_path, "/fig-misc-a.pdf"), pla, width = 9, height = 6, units = "cm", useDingbats = F);
540+
ggsave(paste0(covid_uk_path, "/fig-misc-b.pdf"), plb, width = 9, height = 6, units = "cm", useDingbats = F);
541+
542+
543+
# TOTAL INFECTIONS
544+
1 - quantile(d1[scenario == "Base" & compartment == "S" & t == 702 & region == "United Kingdom"]$value /
545+
d1[scenario == "Base" & compartment == "S" & t == 0 & region == "United Kingdom"]$value, c(0.025, 0.5, 0.975))
546+
547+
# CASES AMONG INFECTIONS
548+
cases = d1[scenario == "Base" & compartment == "cases" & region == "United Kingdom", sum(value), by = run]
549+
subclin = d1[scenario == "Base" & compartment == "subclinical" & region == "United Kingdom", sum(value), by = run]
550+
quantile(cases$V1 / (cases$V1 + subclin$V1), c(0.025, 0.5, 0.975))
507551

552+
# STOPPED EPIDEMICS FROM COMBINED INTERVENTION (DURING CONTROL PERIOD)
553+
r0s[scenario == "Combination", sum(R0 < 1)]
554+
555+
# SPORTS/LEISURE IMPACT
556+
quantile(1 - d6[scenario == "Spectator sports banned" & compartment == "cases" & region == "United Kingdom", sum(value), by = run]$V1 /
557+
d6[scenario == "Background" & compartment == "cases" & region == "United Kingdom", sum(value), by = run]$V1, c(0.025, 0.5, 0.975))
558+
559+
quantile(1 - d6[scenario == "Leisure reduced by 75%" & compartment == "cases" & region == "United Kingdom", sum(value), by = run]$V1 /
560+
d6[scenario == "Background" & compartment == "cases" & region == "United Kingdom", sum(value), by = run]$V1, c(0.025, 0.5, 0.975))
561+
562+
# INTENSIVE INTERVENTIONS IMPACT (DEATHS)
563+
# SPORTS/LEISURE IMPACT
564+
quantile(1 - d3[scenario == "Intensive Interventions" & compartment == "deaths" & region == "United Kingdom", sum(value), by = run]$V1 /
565+
d1[scenario == "Base" & compartment == "deaths" & region == "United Kingdom", sum(value), by = run]$V1, c(0.025, 0.5, 0.975))
508566

509567
# NINGBO EST.
510568
x = rbeta(100000, 6, 140)/rbeta(100000, 126, 1875)
511569
cm_mean_hdi(x)
570+
571+
# COMPARISON WITH UK DEATHS
572+
ukdeaths = fread(
573+
"Area name Area code Area type Reporting date Daily hospital deaths Cumulative hospital deaths
574+
United Kingdom K02000001 UK 2020-03-27 181 759
575+
United Kingdom K02000001 UK 2020-03-26 115 578
576+
United Kingdom K02000001 UK 2020-03-25 41 463
577+
United Kingdom K02000001 UK 2020-03-24 87 422
578+
United Kingdom K02000001 UK 2020-03-23 54 335
579+
United Kingdom K02000001 UK 2020-03-22 48 281
580+
United Kingdom K02000001 UK 2020-03-21 56 233
581+
United Kingdom K02000001 UK 2020-03-20 33 177
582+
United Kingdom K02000001 UK 2020-03-19 41 144
583+
United Kingdom K02000001 UK 2020-03-18 32 103
584+
United Kingdom K02000001 UK 2020-03-17 16 71
585+
United Kingdom K02000001 UK 2020-03-16 20 55
586+
United Kingdom K02000001 UK 2020-03-15 14 35
587+
United Kingdom K02000001 UK 2020-03-14 21
588+
United Kingdom K02000001 UK 2020-03-13 8
589+
United Kingdom K02000001 UK 2020-03-12 8
590+
United Kingdom K02000001 UK 2020-03-11 6
591+
United Kingdom K02000001 UK 2020-03-10 6")
592+
593+
deaths = d1[scenario == "Base" & compartment == "deaths" & region == "United Kingdom", cm_mean_hdi(value), by = t]
594+
ggplot(deaths[t < 60]) +
595+
geom_ribbon(aes(x = ymd("2020-01-29") + t, ymin = lower, ymax = upper), alpha = 0.25) +
596+
geom_line(aes(x = ymd("2020-01-29") + t, y = mean)) +
597+
geom_point(data = ukdeaths, aes(x = ymd(`Reporting date`), y = `Daily hospital deaths`), colour = "red") +
598+
labs(x = "Date", y = "Daily hospital deaths")
599+
600+
# VARIATION SENSITIVITY
601+
# `Combination` = list(contact = c(1.0, 0.5, 0.0, 0.5, 1.0, 0.25, 0.0, 0.25, 0), fIs = rep(0.65, 16))
602+
603+
dists = data.table(x = (0:100)/100);
604+
dists[, d50 := dbeta(x, 10, 10)];
605+
dists[, d25 := dbeta(x, 5, 15)];
606+
dists[, d65 := dbeta(x, 13, 7)];
607+
608+
pll = ggplot(dists, aes(x = x)) +
609+
geom_line(aes(y = d50/max(d50), colour = "Work, other contacts in under-70s"), linetype = "dashed") +
610+
geom_line(aes(y = d25/max(d25), colour = "Work, other contacts in over-70s"), linetype = "dashed") +
611+
geom_line(aes(y = d65/max(d65), colour = "Infectiousness of symptomatic individuals"), linetype = "dashed") +
612+
geom_pointrange(aes(x = 0.50, ymin = 0, ymax = 1, y = 1, colour = "Work, other contacts in under-70s"), size = 0.25, fatten = 0.2) +
613+
geom_pointrange(aes(x = 0.25, ymin = 0, ymax = 1, y = 1, colour = "Work, other contacts in over-70s"), size = 0.25, fatten = 0.2) +
614+
geom_pointrange(aes(x = 0.65, ymin = 0, ymax = 1, y = 1, colour = "Infectiousness of symptomatic individuals"), size = 0.25, fatten = 0.2) +
615+
labs(colour = NULL, x = NULL, y = "Normalised density") +
616+
guides(colour = guide_legend(nrow = 3, byrow = TRUE, reverse = TRUE)) +
617+
theme(legend.position = "bottom")
618+
619+
d7 = rbind(
620+
d2.2[scenario == "Base"],
621+
d2.2[scenario == "Combination 28 shift local"],
622+
d2.2v[scenario == "Combination"]
623+
)
624+
t7 = rbind(
625+
t2.2[scenario == "Base"],
626+
t2.2[scenario == "Combination 28 shift local"],
627+
t2.2v[scenario == "Combination"]
628+
)
629+
d7[scenario == "Combination 28 shift local", scenario := "No variation"]
630+
d7[scenario == "Combination", scenario := "Variation"]
631+
t7[scenario == "Combination 28 shift local", scenario := "No variation"]
632+
t7[scenario == "Combination", scenario := "Variation"]
633+
634+
tb7 = make_table(d7)
635+
pl1 = plot_table(tb7[scenario != "Base"])
636+
save_table(tb7, paste0(covid_uk_path, "/table-variation.csv"));
637+
pl2 = plot_attackrate(t7)
638+
pl3 = plot_epi(d7, t7, (0:10)/10, "2020-01-29")
639+
f = plot_grid(pl1, pl2, pl3, ncol = 1, rel_heights = c(6, 6, 10), labels = c("a", "b", "c"), label_size = 9);
640+
ggsave(paste0(covid_uk_path, "/full-7.pdf"), f, width = 20, height = 22, units = "cm", useDingbats = F);
641+
642+
pla = plot_epi(d7[compartment != "deaths" & compartment != "beds_nonicu"], t7, (0:10)/10, "2020-01-29", "2020-10-15",
643+
colours = cols4[1:3], exclude = "Base");
644+
tb7[statistic == "Cases in peak week", statistic := "Cases in\npeak week"];
645+
tb7[statistic == "Peak ICU beds required", statistic := "Peak ICU beds\nrequired"];
646+
tb7[statistic == "Peak non-ICU beds required", statistic := "Peak non-ICU beds\nrequired"];
647+
tb7[statistic == "Time to peak cases (weeks)", statistic := "Time to peak\ncases (weeks)"];
648+
plb = plot_table(tb7[statistic != "Deaths in peak week" & statistic != "Total subclinical" & scenario != "Base"]) + theme(legend.position = "bottom") +
649+
guides(colour = guide_legend(nrow = 3, byrow = TRUE)) + labs(colour = NULL)
650+
651+
googen = fread(paste0(covid_uk_path, "/../../data/Global_Mobility_Report.csv")); # Google global mobility report
652+
googen = googen[country_region_code == "GB"];
653+
googen = googen[date >= "2020-04-01" & date <= "2020-04-07", mean(workplaces_percent_change_from_baseline), by = sub_region_1]
654+
plc = ggplot(googen) +
655+
geom_histogram(aes(x = 1 + V1/100)) +
656+
xlim(0, 1) +
657+
labs(x = "Relative workplace\nvisits, 1-7 April 2020", y = "Count")
658+
659+
f = plot_grid(pll, pla, plb, plc,
660+
nrow = 1, ncol = 4, rel_widths = c(1.0, 1.3, 1.6, 0.8), labels = c("a", "b", "c", "d"), label_size = 9, align = "hv", axis = "bottom")
661+
ggsave(paste0(covid_uk_path, "/fig-variation.pdf"), f, width = 24, height = 6, units = "cm", useDingbats = F);
662+
ggsave(paste0(covid_uk_path, "/fig-variation.png"), f, width = 24, height = 6, units = "cm");

0 commit comments

Comments
 (0)