report-template.Rmd

---
author: "LSHTM CMMID COVID-19 Working Group"
date: "`r format(Sys.Date(),'%d %b %Y')`"
output:
  bookdown::pdf_document2:
    fig_caption: yes
    toc: false
    citation_package: natbib
    latex_engine: xelatex
    fig_crop: no
mainfont: Arial
geometry:
 - top=10mm
 - left=12mm
 - right=12mm
 - bottom=15mm
params:
  location: Default
  unmitigated: "`r data.frame()`"
  peaks: "`r data.frame()`"
  accs: "`r data.frame()`"
  alls: "`r data.frame()`"
  contactmatrices: "`r matrix()`"
  poppyra: "`r integer()`"
  plothelpers: "`r character()`"
  is_analogy: "`r NA_character_`"
title: "`r sprintf('Modelling projections for COVID-19 epidemic in %s', params$location)`"
bibliography: COVID.bib
biblio-style: unsrtnat
header-includes:
  - \usepackage{booktabs}
  - \usepackage{longtable}
  - \usepackage{array}
  - \usepackage{multirow}
  - \usepackage{wrapfig}
  - \usepackage{float}
  - \usepackage{colortbl}
  - \usepackage{pdflscape}
  - \usepackage{tabu}
  - \usepackage{threeparttable}
  - \usepackage[normalem]{ulem}
  - \usepackage{makecell}
  - \usepackage[square,numbers,sort&compress]{natbib}
  - \newcommand{\blandscape}{\begin{landscape}}
  - \newcommand{\elandscape}{\end{landscape}}
  - \usepackage{titling}
  - \pretitle{\vspace{\bigskipamount}\begin{center}\LARGE}
  - \posttitle{\end{center}}
editor_options: 
  chunk_output_type: console
---
```{r echo=FALSE, include=FALSE}
# TODO for images, in pretitle: \LARGE\includegraphics[width=12cm]{my_graphic.png}
knitr::opts_chunk$set(echo=FALSE, warning=FALSE, message=FALSE)
```

```{r echo=FALSE, include=FALSE}
suppressPackageStartupMessages({
  require(data.table)
  require(ggplot2)
  require(ggthemr)
  require(cowplot)
  require(knitr)
  require(tidyverse)
  require(RColorBrewer)
  require(kableExtra)
})

for (ph in params$plothelpers) {
  if (grepl("rda$", ph)) {
    load(ph)
  } else if (grepl("R$", ph)) {
    source(ph)
  }
}

alls.dt <- params$alls[scen_id %in% c(1,scens)]

milestones <- c(90,180,270,360)
sizerefs <- c(M=6, K=3, 0)
mk <- function(v) {
  i <- which.max(log10(v) >= sizerefs)
  sprintf("%g%s",signif(v/(10^sizerefs[i]),2),names(sizerefs)[i])
}
pasteAnd <- function(v, oxford = TRUE) {
  if (length(v) <= 2) {
    paste(v, collapse = " and ")
  } else {
    opening <- paste(head(v,-1), collapse = ", ")
    if (oxford) opening <- sprintf("%s,",opening)
    pasteAnd(c(opening, tail(v,1)))
  }
}
boldL <- function(kbl) column_spec(row_spec(kbl, 0, bold = TRUE), 1, bold = TRUE)
typicalfmt <- "%s~~(%s)"
```

# Overview

This report provides simulation-based estimates for COVID-19 epidemic scenarios in `r params$location`. In this model, we use a population of `r mk(sum(params$poppyra$pop))` in `r params$location`, with an average age of `r signif(params$poppyra[, sum(per_by_group[-.N]*seq(2,by=5,length.out = .N-1))+per_by_group[.N]*87.5], 3)`, based on [@wpp;@wpprepo].

- We use a transmission model to project the COVID-19 epidemic in `r params$location`; it accounts for the local age structure, as well as uncertainty in transmission due to uncertainty in $R_0$ and randomness in the process.
- We consider an unmitigated epidemic (*i.e.*, no intervention), and then interventions combining  social distancing and shielding interventions to decrease transmission.
- We estimate predicted peak counts and timing for `r params$location` in new symptomatic cases and deaths (incidence), as well as peaks in daily hospital demand (prevalence). We also calculated totals for these values, and compute the impact of the interventions relative to an unmitigated epidemic.
- The totals for hospital demand are calculated in terms of person-days, which can be translated into the total required expenditure on hospitalised care at different levels.
- Populations are based on demographics of `r params$location` and estimated contact patterns, and assume a single national population
- The transmission model is based on global outbreak data, which have primarily been observed in Asia, Europe, and North America.
- Severity rates are rescaled from the global setting. We assume that due to comorbidities, lower-middle income populations will have less benefit from age and will in general have higher rates of symptomatic disease.
- We start the epidemic simulation with 50 initial infections, proportionally distributed according to the age demography of `r params$location`. We refer to this day 0 as the *initial introduction*. This corresponds roughly to when there were roughly 50 active infections in `r params$location`; due to under-reporting and subclinical disease, this cannot be directly determined from reported cases.
- These reports will be updated as information changes, and if `r params$location`-specific information on hospitalisation rates become available.
- All numerical results are reported to two significant figures only, and are reported as interquartile (IQR) and 95% intervals.
- Detailed modelling methods are given at the end of the report.

*Please direct correspondence regarding these reports to [carl.pearson@lshtm.ac.uk](mailto:carl.pearson@lshtm.ac.uk). Source available at [https://github.com/cmmid/covidm_reports](https://github.com/cmmid/covidm_reports).*

# Unmitigated COVID-19 Epidemic Trajectory

The panels in Fig. \@ref(fig:tsunmit) show aggregate and by-age group outcomes for daily incidence of COVID-19 cases and deaths in `r params$location`, and daily demand for hospital care. The ranges for peak timing and values for those outcomes are summarised in Table \@ref(tab:unpeakstable) (by 50% and 95% forecast intervals); Table \@ref(tab:uncumulative1) covers total counts by `r pasteAnd(milestones/30)` months from initial introduction in `r params$location`.

```{r}
unmit.q <- meltquantiles(params$alls[scen_id == 1])
unmit.q$age <- factor(unmit.q$age, levels(unmit.q$age)[c(6,1:5)])
reage <- function(a, lvls = levels(unmit.q$age)) factor(a, levels = lvls, ordered = TRUE)
unmit.q$compartment <- factor(unmit.q$compartment, levels = levels(unmit.q$compartment)[c(3,1,6,4:5,2)], ordered = TRUE)
recomp <- function(cmp, lvls = levels(unmit.q$compartment)) factor(cmp, levels = lvls, ordered = TRUE)

reord <- function(dt) dt[, age := reage(age) ][, compartment := recomp(compartment) ]

unmit.pks <- params$alls[scen_id == 1 & compartment != "E"][,{
  res <- apply(.SD, 2, max)
  names(res) <- colnames(.SD)
  as.list(res)
}, by=.(compartment, age, scenario = int.factorize(scen_id))]
reord(unmit.pks)

quotepeaks <- unmit.pks[age == "all"]

up <- function(v, sf = 2) {
  if (length(v) && v!=0) {
  ref <- max(floor(log10(v))-(sf-1),0)
  ceiling(v/(10^ref))*(10^ref) } else v
}

down <- function(v, sf = 2) {
  if (length(v)) {
    if (v>0) {
      ref <- max(floor(log10(v))-(sf-1),0)
      floor(v/(10^ref))*(10^ref)
    } else if (v<0) {
      -up(-v, sf)
    } else v
  } else v
}

numfmt <- function(n) if(n>=1e6) {
  sprintf("%sM",format(n/1e6, big.mark = ",", scientific = FALSE))
} else if(n>=1e3) {
  sprintf("%sK",format(n/1e3, big.mark = ",", scientific = FALSE))
} else {
  format(n, big.mark = ",", scientific = FALSE)
}

rngprint <- function(r, prefix="") {
  same <- down(r[1])==up(r[2])
  ifelse(
    same,
    numfmt(signif(r[1],2)),
    sprintf(
      "%s%s - %s", prefix,
      numfmt(down(r[1])),
      numfmt(up(r[2]))
    )
  )
}
bullettext <- function(filterdt) {
  paste0(rngprint(filterdt[,c(lo,hi)],"between ")," (95%: ",rngprint(filterdt[,c(lo.lo,hi.hi)]),")")
}
tabltext <- function(filterdt) {
  sprintf(typicalfmt, rngprint(filterdt[,c(lo,hi)]), rngprint(filterdt[,c(lo.lo,hi.hi)]))
}
outcomesc <- params$alls[
  scen_id == 1 & age == "all" & compartment != "E"
][,{
  vs <- colSums(.SD)
  names(vs) <- colnames(.SD)
  as.list(vs)
  },by=compartment,.SDcols = -c("age", "t", "scen_id")
]
outcomest <- params$alls[
  scen_id == 1 & age == "all" & compartment != "E"
][,{
  vs <- t[apply(.SD, 2, which.max)]
  names(vs) <- rev(colnames(.SD))
  as.list(vs)
  },by=compartment,.SDcols = -c("age", "t", "scen_id")
]
```

In an unmitigated epidemic, we anticipate that in `r params$location`, one year after the initial introduction, the total number of symptomatic cases will have been `r bullettext(outcomesc[compartment == "cases"])`, the total deaths `r bullettext(outcomesc[compartment == "death_o"])`, and total hospitalised person-days `r bullettext(outcomesc[compartment == "hosp_p"])` with `r bullettext(outcomesc[compartment == "icu_p"])` of those person-days requiring critical care facilities. Table \@ref(tab:sumtable) shows the size and timing of peaks for these main outcomes.

```{r sumtable}
tabl <- data.table(
  `Outcome` = c(
    "Incidence of Symptomatic Cases",
    "All Hospital Demand",
    "General Hospital Demand",
    "Critical Care Demand",
    "Incidence of Deaths"
  ),
  `Peak Day IQR (95\\% interval)` = c(
    tabltext(outcomest[compartment == "cases"]),
    tabltext(outcomest[compartment == "hosp_p"]),
    tabltext(outcomest[compartment == "nonicu_p"]),
    tabltext(outcomest[compartment == "icu_p"]),
    tabltext(outcomest[compartment == "death_o"])
  ),
  `Peak Number IQR (95\\% interval)` = c(
    tabltext(quotepeaks[compartment == "cases"]),
    tabltext(quotepeaks[compartment == "hosp_p"]),
    tabltext(quotepeaks[compartment == "nonicu_p"]),
    tabltext(quotepeaks[compartment == "icu_p"]),
    tabltext(quotepeaks[compartment == "death_o"])
  )
)
boldL(kable_styling(
  kable(
    tabl,
    caption=sprintf("Overview of estimated peak impacts based on the population in %s. Rates of symptoms and severe outcomes are based on global outbreak data, which have primarily been observed in Asia, Europe, and North America. Timing is in days from initial introduction to %s.", params$location, params$location),
    booktabs = TRUE, escape = FALSE
  ),
  latex_options = c("striped","hold_position")
))
```

\blandscape

```{r tsunmit, fig.height=6, out.width="10in", fig.width=10, fig.cap=sprintf("Unmitigated epidemics in %s. Time series and peak sizes for new symptomatic cases and deaths, general ward hospital demand, critical care hospital demand, and total hospital demand (general ward plus critical care). Both summary and by age group results are provided; population and underlying age groups from demography of %s. The different line transparency denotes different simulation quantiles; the lightest lines correspond to the 95\\%% interval, then darker the IQR, then the solid line the median. Cases and deaths reflect new daily events (incidence), while hospital demand reflects the number in hospital on that day (prevalence).", params$location, params$location)}

age.annotate <- function(l) sprintf("%s\n%s",l,ifelse(l=="all","ages","years"))
stagger <- function(l) {
  l[seq(1,length(l),by=2)] <- sprintf("%s\n", l[seq(1,length(l),by=2)])
  l[seq(2,length(l),by=2)] <- sprintf("\n%s", l[seq(2,length(l),by=2)])
  l
}

xscale <- function(
  name=sprintf("Months since initial introduction in %s", params$location),
  breaks = milestones,
  labels = function(b) sprintf("%i\n", b/30)
) scale_x_t(name=name, labels = labels, breaks = breaks)

ftsize = 8

sharedelements <- list(
  scale_color_discrete(guide = "none"),
  scale_alpha_quantile(),
  theme_minimal(base_size = ftsize),
  scale_y_continuous(labels = scales::label_number_si()),
  theme(
    panel.spacing.y = unit(12, "pt")
  )
)

plot_grid(ggplotqs(unmit.q[compartment != "E"], aes(
  alpha = variable,
  color = "Unmitigated"
)) +
  facet_grid(
    compartment ~ age,
    scale = "free_y", labeller = fct_labels(
      age = age.annotate
    ),
    switch = "y"
  ) +
  sharedelements +
  xscale() +
  coord_cartesian(ylim = c(0, NA), expand = F) +
  theme(
    strip.placement = "outside",
    axis.title.y = element_blank(),
    plot.margin = margin(l = 12)
  ),
int.peaks.all(unmit.pks, aes(x=age), range.scale = c(1, 1.5, 2)) + facet_grid(
  compartment ~ "\nPeak Values", scales = "free_y"
) +
  scale_x_discrete("Age Group", labels = age.annotate) +
  sharedelements +
  theme(
      axis.title.x = element_text(),
      axis.text.x = element_text(),
      axis.title.y = element_blank(),
      strip.text.y = element_blank(),
      plot.margin = margin(r=12)
    ),
  axis = "l", align = "v", rel_widths = c(3, 2)
)
```

\elandscape

```{r unpeakstable}
tabpeak <- alls.dt[scen_id == 1 & compartment != "E",{
  v.ll = max(lo.lo); t.ll = t[which.max(lo.lo)]
  v.lo = max(lo); t.lo = t[which.max(lo)]
  v.md = max(med); t.md = t[which.max(med)]
  v.hi = max(hi); t.hi = t[which.max(hi)]
  v.hh = max(hi.hi); t.hh = t[which.max(hi.hi)]
  .(
    t=sprintf(typicalfmt, rngprint(c(t.hi,t.lo)), rngprint(c(t.hh, t.ll))),
    v=sprintf(typicalfmt, rngprint(c(v.lo,v.hi)), rngprint(c(v.ll,v.hh)))
  )},
  by=.(age, compartment)
]
reord(tabpeak)

tabpeak <- tabpeak[order(compartment, age),.(
  `Peaks`=c(
    cases="Incident Cases",
    death_o="Incident Deaths",
    hosp_p="Hospital Demand",
    nonicu_p="Non-critical Demand",
    icu_p="Critical Demand"
  )[as.character(compartment)],
  `Age Group`=sprintf("%s %s",as.character(age),ifelse(age=="all","ages","years")),
  `Peak Day, IQR (95\\% interval)`=t,
  `Peak Number, IQR (95\\% interval)`=v
)]

mkkbl <- function(tbl, b = TRUE, ...) boldL(row_spec(
  collapse_rows(
  kable(
    tbl,
    ...
    ), columns = 1, latex_hline = "major"
  ),
  seq(1, by=length(unique(tbl[[2]])), length.out = length(unique(tbl[[1]]))), bold = b
))

mkkbl(
  tabpeak,
  caption=sprintf("Peak timing and values for main outcomes in %s for the unmitigated scenario, by age group. Timing is days from the date of initial introduction.", params$location),
  escape = FALSE
)
```

```{r uncumulative1}
outcome_cumulative <- c(
    cases="Cases",
    death_o="Deaths",
    hosp_p="Hospital\nPerson-Days",
    nonicu_p="Non-critical\nPerson-Days",
    icu_p="Critical\nPerson-Days",
    E="Infections"
  )

unmitacc <- params$unmitigated[compartment != "E"][
  order(t),.(t, value=cumsum(value)),by=.(run, age, compartment)
][t %in% milestones,{
  qs <- quantile(value, probs = refprobs)
  iqr <- rngprint(qs[c(2,4)])
  i95 <- rngprint(qs[c(1,5)])
  sprintf(typicalfmt, iqr, i95)
},by=.(age, compartment, t)]

reord(unmitacc)

unmitacctable1 <- dcast(unmitacc, age + compartment ~ t, value.var="V1")[
order(compartment, age),.(
  `Totals`=linebreak(outcome_cumulative[as.character(compartment)]),
  `Age Group`=sprintf("%s %s",as.character(age),ifelse(age=="all","ages","years")),
  `3 months`=`90`,
  `6 months`=`180`,
  `9 months`=`270`,
  `12 months`=`360`
)
]

landscape(mkkbl(
    unmitacctable1,
    caption=sprintf("Total counts for main outcomes in %s for the unmitigated scenario, by age group; counts are evaluated at %s months from initial introduction.", params$location, pasteAnd(milestones/30)),
    escape = FALSE
))
```

# Intervention Scenarios

The panels in Fig. \@ref(fig:tsintervention) compare unmitigated epidemics with 5 different potential interventions in `r params$location`. The IQR and 95% intervals for peak timing and values for those outcomes are summarised in Table \@ref(tab:intpeakstable). Fig. \@ref(fig:cumtsint) shows the relative trajectories in total outcomes, Table \@ref(tab:intcumulative1) shows cumulative reductions due to the interventions at `r pasteAnd(milestones/30)` months from initial introduction. Note that these reduction may rise and decline with time: as Fig. \@ref(fig:cumtsint) shows, some interventions have large initial impact by delaying an epidemic, but ultimately many of the cases still happen.

- We assumed these interventions are applied at the country level.
- General physical distancing was implemented as a reduction in contacts, and thus transmission, for all interactions outside the household. We assumed no change in transmission within the household.
- Shielding was implemented by stratifying the population in one shielded and one unshielded compartment. Shielding applies to those aged 60+ years. We assume shielding of this population has a coverage fraction, and reduction fraction; we refer to shielding interventions by coverage / reduction, so "80/80 shielding" is 80% coverage of 80% contact reduction. We show results for 40/80 and 80/80 here. Other results available upon request. In `r params$location`, 80% coverage corresponds to `r mk(params$poppyra[age == "60+", sum(pop)*.8])` individuals  or `r mk(params$poppyra[age == "60+", round(unique(per_by_age)*.8*100)])`% of the population, with half of those values for 40% coverage.
- To be effective, shielding must be maintained until the epidemic is over and must ensure that there is not substantially increased mixing amongst the shielded population.

The % net reduction to cases, deaths, and hospitalised person-days in these scenarios, one year after initial introduction in `r params$location` is summarised in Table \@ref(tab:intsumtabl):

```{r intsumtabl}
effref <- params$accs[compartment != "E"][metric == "effectiveness" & age == "all" & t == 360 & scen_id %in% scens]
effref[, scenario := int.factorize(scen_id)]
reord(effref)
#redref <- params$accs[compartment != "E"][metric == "reduction" & age == "all" & t == 360 & scen_id %in% scens]
#reord(redref)
#redref[, longname := int.factorize(scen_id) ]
efftext <- function(r) {
  v <- round(r*100)
  sprintf("%i-%i\\%% (%i-%i\\%%)",v[2],v[4],v[1],v[5])
}
efftab <- dcast(
  effref[,.(etext=efftext(c(lo.lo,lo,med,hi,hi.hi))),keyby=.(compartment, scenario)],
  linebreak(outcome_cumulative[as.character(compartment)]) ~ linebreak(gsub("%","\\%",gsub(" &\n",",\n", scenario, fixed = T), fixed = T)), value.var = "etext"
)
colnames(efftab)[1] <- "Total"
setcolorder(efftab, c(1,2,4,3,5,6))
boldL(kable_styling(
  kable(
    efftab, booktabs = TRUE, caption = "Scenario Effectiveness Summary",
    escape = FALSE, 
  )
  , latex_options = c("striped", "hold_position")
))
```

\blandscape

```{r tsintervention, fig.height=6, out.width="10in", fig.width=10, fig.cap=sprintf("Epidemics in %s with interventions for all ages. Time series and peak counts for new symptomatic cases and deaths, general ward hospital demand, critical care hospital demand, and total hospital demand (general ward plus critical care). The different line transparency denotes different simulation quantiles; the lightest lines correspond to the 95\\%% interval, then darker the IQR, then the solid line the median. Cases and deaths reflect new daily events (incidence), while hospital demand reflects the number in hospital on that day (prevalence).", params$location)}
pks <- alls.dt[age=="all" & compartment != "E",{
  res <- apply(.SD, 2, max)
  names(res) <- colnames(.SD)
  as.list(res)
}, by=.(compartment, age, scenario = int.factorize(scen_id))]
reord(pks)

sharedelements <- list(
  scale_alpha_quantile(),
  coord_cartesian(ylim = c(0, NA), expand = F),
  scale_y_continuous(labels = scales::label_number_si()),
  theme_minimal(base_size = ftsize),
  theme(
    panel.spacing.y = unit(12, "pt")
  )
)

int.qs <- meltquantiles(alls.dt[age=="all" & compartment != "E"])
int.qs[, scenario := int.factorize(scen_id) ]
reord(int.qs)

plot_grid(
ggplotqs(int.qs, aes(
   color=scenario, group=variable, alpha=variable
)) +
  facet_grid(compartment ~ scenario, scale = "free_y", switch = "y", labeller = fct_labels(
    #scen_id = { res <- levels(int.factorize(c(1, scens))); names(res) <- c(1, scens); res }
  )) +
  xscale() +
  scale_color_discrete(guide = "none") +
  sharedelements +
  theme(
    strip.placement = "outside",
    axis.title.y = element_blank(),
    plot.margin = margin(l = 12)
  ),
int.peaks.all(pks, range.scale = c(1, 1.5, 2)) + facet_grid(
  compartment ~ "\nPeak Values", scales = "free_y", labeller = labeller(compartment = function(l) rep("",length(l)))
) + sharedelements +
  scale_x_discrete("\n", labels = function(bs) rep("", length(bs))) +
  theme(
      axis.title.x = element_text(),
      axis.text.x = element_text(),
      axis.title.y = element_blank(),
      legend.key.size = unit(24, "pt"),
      strip.text.y = element_blank(),
      plot.margin = margin()
    ),
  axis = "l", align = "v", rel_widths = c(3, 2)
)
```

\elandscape

\blandscape

```{r cumtsint, fig.height=6, out.width="10in", fig.width=10, fig.cap=sprintf("Trajectories for cumulative outcomes in %s. The ribbons correspond to the 95\\%% intervals (lightest) and IQR (darker), with solid lines for the median. Note that some interventions have large initial impact by delaying an epidemic, but ultimately many of the cases still happen; when comparing interventions to an unmitigated epidemic, this leads to large initial reduction in total cases, which later declines.", params$location)}

accs.dt <- params$alls[compartment != "E" & age == "all" & scen_id %in% c(1, scens)]

cum.qs <- accs.dt[order(t),{
  .(t, lo.lo=cumsum(lo.lo), lo=cumsum(lo), med=cumsum(med), hi=cumsum(hi), hi.hi = cumsum(hi.hi))
}, keyby=.(scenario = int.factorize(scen_id), compartment, age)]

reord(cum.qs)

pks <- cum.qs[t %in% milestones, {
  res <- apply(.SD, 2, max)
  names(res) <- colnames(.SD)
  as.list(res)
}, by=.(compartment, age, scenario, t)]
reord(pks)

sharedelements <- list(
  scale_y_continuous(labels = scales::label_number_si()),
  theme_minimal(base_size = ftsize),
  theme(
    panel.spacing.y = unit(12, "pt")
  )
)

range.scale = c(1, 1.5, 2)

plot_grid(
ggplot(cum.qs) + aes(
  t, med,
  fill = scenario
) +
  geom_ribbon(aes(ymin=lo.lo, ymax=hi.hi, alpha = "r95")) +
  geom_ribbon(aes(ymin=lo, ymax=hi, alpha = "r50")) +
  geom_line(aes(alpha = "med", color = scenario)) +
  facet_grid(compartment ~ "Cumulative Total Count", scale = "free_y", switch = "y", labeller = fct_labels(compartment = c(
    cases="Total Cases",
    hosp_p="All Hospital\nPerson-Days",
    nonicu_p="General Hospital\nPerson-Days",
    icu_p="Critical Care\nPerson-Days",
    death_o="Total Deaths"
  ))) +
  coord_cartesian(ylim = c(0, NA), xlim = c(0, 360), expand = F) +
  xscale() +
  scale_color_discrete(guide = "none", aesthetics = c("fill","color")) +
  scale_alpha_manual(values=c(r95=0.2,r50=0.2,med=1), guide = "none") +
  sharedelements +
  theme(
    strip.placement = "outside",
    axis.title.y = element_blank(),
    plot.margin = margin(l = 12)
  ),
ggplot(pks) +
    aes(as.numeric(scenario)+(t/30-7.5)/18, med, color=scenario) +
    geom_linerange(aes(ymin=lo.lo, ymax=hi.hi, alpha="hi.hi"),  size = range.scale[1]) +
    geom_linerange(aes(ymin=lo, ymax=hi, alpha="hi"), size = range.scale[2]) +
    geom_point(aes(y=med, alpha="med"), size = range.scale[3]) +
    scale_color_discrete(NULL, guide = "none") +
    coord_cartesian(ylim = c(0, NA), expand = F) + facet_grid(
  compartment ~ "Total Values at 3, 6, 9, and 12 Months", scales = "free_y", labeller = labeller(compartment = function(l) rep("",length(l)))
) + sharedelements +
  coord_cartesian(ylim = c(0, NA), expand = F, clip = "off") +
  scale_x_continuous("", breaks = 1:c(length(scens)+1), labels = int.factorize(c(1,scens)) ) +
  scale_alpha_manual(values = c(0.2, 0.2, 1), guide = "none") +
  theme(
      axis.title.y = element_blank(),
      legend.key.size = unit(24, "pt"),
      strip.text.y = element_blank(),
      plot.margin = margin(r=unit(24,"pt"))
    ),
  axis = "l", align = "v", rel_widths = c(2, 3)
)
```

\elandscape

```{r intpeakstable}
tabpeak <- alls.dt[age=="all" & compartment != "E",{
  v.ll = max(lo.lo); t.ll = t[which.max(lo.lo)]
  v.lo = max(lo); t.lo = t[which.max(lo)]
  v.md = max(med); t.md = t[which.max(med)]
  v.hi = max(hi); t.hi = t[which.max(hi)]
  v.hh = max(hi.hi); t.hh = t[which.max(hi.hi)]
  .(
    t=sprintf(typicalfmt, rngprint(c(t.hi,t.lo)), rngprint(c(t.hh, t.ll))),
    v=sprintf(typicalfmt, rngprint(c(v.lo,v.hi)), rngprint(c(v.ll,v.hh)))
  )},
  by=.(scenario=int.factorize(scen_id), compartment, age)
]
reord(tabpeak)

inttab <- tabpeak[order(scenario, compartment)][,.(
  Scenario=linebreak(gsub("%","\\%",gsub(" &\n",",\n", scenario, fixed = T), fixed = T)),
  `All ages`=c(
    cases="Incident Cases",
    death_o="Incident Deaths",
    hosp_p="Hospital Demand",
    nonicu_p="Non-critical Demand",
    icu_p="Critical Demand"
  )[as.character(compartment)],
  `Peak day IQR (95\\% interval)`=t,
  `Peak number IQR (95\\% interval)`=v
)]

mkkbl(inttab, caption=sprintf("Peak timing and values for main outcomes in %s, by intervention scenario. Timing is from the date of ongoing community spread with 50 infections.", params$location), b = FALSE, escape = FALSE)
```

```{r intcumulative1}
intacc <- params$accs[compartment != "E"][metric == "reduction" & age == "all" & scen_id %in% scens][,{
  iqr <- rngprint(c(lo, hi))
  i95 <- rngprint(c(lo.lo,hi.hi))
  sprintf(typicalfmt, iqr, i95)
},by=.(age, compartment, t, scenario = int.factorize(scen_id))]

reord(intacc)

intacctable1 <- dcast(intacc, scenario + compartment ~ t, value.var="V1")[order(scenario, compartment)][
,.(
  Scenario=linebreak(gsub("%","\\%",gsub(" &\n",",\n",scenario, fixed = T), fixed = T)),
  Outcome=gsub("\\n"," ",outcome_cumulative)[as.character(compartment)],
  `3 months`=`90`,
  `6 months`=`180`,
  `9 months`=`270`,
  `12 months`=`360`
)
]

landscape(mkkbl(
  intacctable1,
  caption=sprintf("Net reductions for main outcomes in %s by intervention scenario; reductions are evaluated at %s months. Note that net reductions may decline with time; as shown in Fig. 3, some interventions have a large initial impact, but ultimately still allow many cases.", params$location, pasteAnd(milestones/30)),
  escape = FALSE,
  b = FALSE
))
```

# Methods, data and assumptions

We used the stochastic age-structured dynamic transmission model reported in [@davies2020].

## Dynamic transmission model 

We used a stochastic compartmental model stratified into 5-year age bands, with individuals classified according to current disease status (Fig. \@ref(fig:modelpars))) and transmission between groups based on social mixing patterns [@daviesage2020;@mossong2008social]. After infection with SARS-CoV-2 in the model, susceptible individuals pass through a latent period before becoming infectious, either with a preclinical and then clinical infection, or with a subclinical infection, before recovery or isolation. We refer to those infections causing few or no symptoms as subclinical. We assume older individuals are more likely to show clinical symptoms [@daviesage2020].

```{r modelpars, echo=FALSE, fig.height=2.5, fig.cap=sprintf("Population Structure%s", ifelse(is.na(params$is_analogy), "", paste0("; contact matrix based on ",params$is_analogy))) }
plotter <- function(pop) {
  p <- ggplot(pop) +
    aes(x = age, y = pop / 1000) +
    geom_col(fill = "#cc3366") +
    coord_flip() +
    labs(x = NULL, y = "Population (thousands)") +
    theme_minimal()
  return(p)
}
sp_contact_matrices <- function(mat, legpos = "right", is_analogy = params$is_analogy) {
  mat <- Reduce(function(l, r) l + r, mat)
  data <- data.table(reshape2::melt(mat))
  names(data) <- c("From", "To", "Contacts")
  ggplot(data) +
    geom_raster(aes(x = To, y = From, fill = Contacts)) +
    scale_fill_viridis_c() +
    theme_bw() +
    theme(
      axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5),
      strip.background = element_rect(fill = "white", colour = "white"),
      legend.position = legpos,
    )
}
plot_grid(plotter(poppyra), sp_contact_matrices(params$contactmatrices), ncol = 2)
```

## Key model parameters 

As documented in [@davies2020], we collated multiple sources of evidence to estimate key model parameters. In a meta-analysis, we estimated that the basic reproduction number, $R_0$, was 2.7 (95% credible interval: 1.6–3.9) across settings without substantial control measures in place ($R_0$ describes the average number of secondary infections caused by a typical primary infection in a completely susceptible population). We derived age-stratified case fatality ratios (CFR) to estimate a CFR that ranged substantially across age groups, from 0.1% in the 20–29 age group to 7.7% in the over-80 age group.

## Key parameters of the transmission model 

We used a serial interval of 6.5 days based on published studies [@li2020early;@bi2020epidemiology;@nishiura2020serial], and assumed that the length of the preclinical period was 30% of the total period of clinical infectiousness [@liu2020contribution]. From this, we fixed the mean of the latent period to 4 days, the mean duration of preclinical infectiousness to 1.5 days, and the mean duration of clinical infectiousness to 3.5 days. The basic reproduction number $R_0$ was estimated by synthesizing the results of a literature review (Fig. S1 in [@davies2020]). For each reported value of the basic reproduction number, we matched a flexible PERT distribution (a shifted beta distribution parameterised by minimum, maximum, and mode) to the median and confidence interval reported in each study. We sampled from the resulting distributions, weighting each study equally, to obtain estimates of $R_0$ for our simulations. The age-specific clinical fraction, $y_i$, was adopted from an estimate based on case data from 6 countries [@daviesage2020], and the relative infectiousness of subclinical cases, $f$, was assumed to be 50% relative to clinical cases, as assumed in a previous study [@daviesage2020].

## Hospital burden estimation

To calculate ICU and non-ICU beds in use through time, we scaled age-stratified symptomatic cases by age-specific hospitalisation and critical outcome probability, then summed to get the total number of hospitalised and critical cases. We then distributed hospitalised cases over time based on expected time of hospitalisation and duration admitted. We assumed gamma-distributed delays, with the shape parameter set equal to the mean, for: delay from symptom onset to hospitalisation of mean 7 days (standard deviation 2.65) [@linton2020incubation;@cao2020trial]; delay from hospitalisation to discharge / death for non-ICU patients of mean 8 days (s.d. 2.83) [@NHSDigital]; delay from hospitalisation to discharge / death for ICU patients of mean 10 days (s.d. 3.16) [@cao2020trial]; and delay from onset to death of mean 22 days (s.d. 4.69) [@linton2020incubation;@cao2020trial]. We calculated the age-specific case fatality ratio based on data from the COVID-19 outbreak in China and on the Diamond Princess cruise ship. We first calculated the naive case fatality ratio, nCFR, (*i.e.* deaths/cases) for each age group, then scaled down the naive CFR based on a correction factor estimated from data from the Diamond Princess [@russell2020] to give an adjusted CFR. We then calculated risk of hospitalisation based on the ratio of severe and critical cases to cases (18.5%) and deaths to cases (2.3%) in the early China data, which we took to imply 8.04 times more hospitalisations than deaths in each age group. We assumed all age groups had a 30% risk of requiring critical care if hospitalised [@cao2020trial].

## Severity Shift

To account for increased prevalence of comorbidities in lower-middle income countries, we shifted the age-specific outcome probabilities by ten years. *E.g.*, the 0-4 year old category in our model has the same risk of disease as we estimated elsewhere for 10-14 year olds. We also applied a 1.5 relative risk modifier to the case fatality ratio, this has the net effect of increasing the death rates relative to hospitalisation rated and infections.

# Acknowledgements

The following authors are part of the Centre for Mathematical Modelling of Infectious Disease (CMMID) COVID-19 working group; each contributed in processing, cleaning and interpretation of data, interpreted findings, contributed to the manuscript, and approved the work for publication: Emily S Nightingale, James D Munday, Graham Medley, Hamish P Gibbs, Sam Abbott, Rein M G J Houben, Kathleen O'Reilly, Kiesha Prem, Akira Endo, Samuel Clifford, Mark Jit, Simon R Procter, Nikos I Bosse, Kevin van Zandvoort, Anna M Foss, Alicia Rosello, Quentin J Leclerc, Sebastian Funk,
Stéphane Hué, Eleanor M Rees, David Simons, Christopher I Jarvis, Carl A B Pearson, Adam J Kucharski, Petra Klepac, Joel Hellewell, Arminder K Deol, Rachel Lowe, Nicholas G. Davies, Charlie Diamond, Damien C Tully, Gwenan M Knight, Jon C Emery, Billy J Quilty, Yang Liu, W John Edmunds, Megan Auzenbergs, C Julian Villabona-Arenas, Katherine E. Atkins, Timothy W Russell, Fiona Yueqian Sun, Stefan Flasche, Rosalind M Eggo, Thibaut Jombart, Amy Gimma, and Sophie R Meakin.

Contributing authors gratefully acknowledge funding of the NTD Modelling Consortium by the Bill and Melinda Gates Foundation (OPP1184344) and via other grants (INV-003174). They also acknowledge support from Elrha’s Research for Health in Humanitarian Crises (R2HC) Programme, which aims to improve health outcomes by strengthening the evidence base for public health interventions in humanitarian crises. The R2HC programme is funded by the UK Government (DFID), the Wellcome Trust, and the UK National Institute for Health Research (NIHR). This work was also supported by the Department for International Development/Wellcome Epidemic Preparedness - Coronavirus research programme (ref. 221303/Z/20/Z). Additionally, the authors acknowledge Global Challenges Research Fund (GCRF) project ‘RECAP’ managed through RCUK and ESRC (ES/P010873/1). Finally, this work was also supported by grants from HDR UK (grant: MR/S003975/1), MRC (grant: MC\_PC 19065), NIHR (16/137/109), NIHR: Health Protection Research Unit for Modelling Methodology (HPRU-2012-10096), and the European Commission (101003688).

The following funding sources are acknowledged as providing funding for the working group authors: Alan Turing Institute (AE), BBSRC LIDP (BB/M009513/1: DS), Bill & Melinda Gates Foundation (INV-003174: KP, MJ, YL; NTD Modelling Consortium OPP1184344: GM, CABP; OPP1180644: SRP; OPP1183986: ESN; OPP1191821: KO’R, MA; ), ERC Starting Grant (#757688: JV-A, KEA; #757699: JCE, RMGJH), European Commission (101003688: KP, MJ, WJE, YL), Global Challenges Research Fund (ES/P010873/1: AG, CIJ), Nakajima Foundation (AE), NIHR (16/137/109: CD, FYS, MJ, YL, BQ; HPRU Modelling Methodology: TJ; HPRU-2012-10096: NGD; PR-OD-1017-20002: AR), RCUK/ESRC (ES/P010873/1: TJ), Royal Society (Dorothy Hodgkin Fellowship: RL; RP\\EA\\180004: PK), UK DHSC/UK Aid/NIHR (ITCRZ 03010: HPG), HDR UK Innovation Fellowship (MR/S003975/1: RME), UK MRC (LID DTP MR/N013638/1: EMR, QJL; MR/P014658/1: GMK), UK Public Health Rapid Support Team (TJ), Wellcome Trust (206250/Z/17/Z: AJK, TWR; 208812/Z/17/Z: SC, SF; 210758/Z/18/Z: JDM, JH, NIB, SA, SFunk, SRM), No funding (AKD, AMF, DCT, SH).

The views expressed in this publication are those of the author(s) and not necessarily those of any of the listed funding sources.

\begin{figure}[h]
    \centering
\includegraphics[width=.2\linewidth]{wellcome-logo-black.jpg}
\includegraphics[width=.2\linewidth]{UK-AID-Standard-RGB.jpg}
\end{figure}

# References