Abstract
As of 1st June 2020, the US Centers for Disease Control and Prevention reported 104,232 confirmed or probable COVID-19-related deaths in the US. This was more than twice the number of deaths reported in the next most severely impacted country. We jointly modelled the US epidemic at the state-level, using publicly available death data within a Bayesian hierarchical semi-mechanistic framework. For each state, we estimate the number of individuals that have been infected, the number of individuals that are currently infectious and the time-varying reproduction number (the average number of secondary infections caused by an infected person). We used changes in mobility to capture the impact that non-pharmaceutical interventions and other behaviour changes have on the rate of transmission of SARS-CoV-2. On 1st June, we estimated that Rt was only below one in 23 states. We also estimated that 3.7% [3.4%-4.0%] of the total population of the US had been infected, with wide variation between states, and approximately 0.01% of the population was infectious. We demonstrate good 3 week model forecasts of deaths with low error and good coverage of our credible intervals.
1 Introduction
The first death caused by COVID-19 in the United States is currently believed to have occurred in Santa Clara County, California on the 6th February [1]. Throughout March 2020, US state governments implemented a variety of non-pharmaceutical interventions (NPIs), such as school closures and stay-at-home orders, to limit the spread of SARS-CoV-2 and ensure the number of severe COVID-19 cases did not exceed the capacity of the health system. In April 2020, the number of deaths attributed to COVID-19 in the United States (US) surpassed that of Italy [2]. Courtemanche et al. [3] used an event-study model to determine that such NPIs were successful in reducing the growth rate of COVID-19 cases across US counties. We similarly seek to estimate the impact of NPIs on COVID-19 transmission, but through a semi-mechanistic Bayesian model that reflects the underlying process of disease transmission and relies on mobility data released by companies such as Google [4].
Mobility measures revealed stark changes in behaviour following the large-scale government interventions in the first stage of the epidemic, with individuals spending more time at home and correspondingly less time at work, at leisure centres, shopping, and on public transit [4, 5]. As states continued to ease the stringency of their NPIs in the end of June, policy decisions relied on the interaction between mobility and NPIs and their subsequent impact on transmission, alongside other measures to track and curtail SARS-CoV-2 transmission.
We introduced a new Bayesian statistical framework for estimating the rate of transmission and attack rates for COVID-19 in Flaxman et al. [6]. In that paper, we inferred the time-varying reproduction number, Rt, or the average number of people an infected person will infect over time. We calculated the number of new infections through combining previous infections with the generation interval (the distribution of times between infections) and chose the number of deaths to be a function of the number of infections and the infection fatality ratio (IFR). We estimated the posterior probability of our parameters given the observed data, while incorporating prior uncertainty. This made our approach empirically driven, whilst incorporating uncertainty. This approach has also been implemented for Italy [7] and Brazil [8].
In this paper, we extend the Flaxman et al. [6] framework to model transmission in the US at the state level and include reported cases in our model. We parameterise Rt as a function of several mobility types and include an autoregressive term to capture changes in transmission that are decoupled from mobility, for example hand-washing, social distancing and changes in transmission that are decoupled from mobility. We utilise partial pooling of parameters, where information is shared across all states to leverage as much signal as possible, but individual effects are also included for state and region-specific idiosyncrasies. In this paper, we infer plausible upper and lower bounds (Bayesian credible interval summaries of our posterior distribution) of the total population that had been infected by COVID-19 on 01 June 2020 (also called the cumulative attack rate or attack rate) and estimate the effective number of individuals currently infectious given our generation distribution. We also present effect sizes of the mobility covariates and make short term forecasts, which we compare with reality throughout June. Details of the data sources and a technical description of our model are found in Sections 4 and 5 respectively. General limitations of our approach are presented in the conclusions.
2 Results
2.1 Infections
The percentage of the total population across the US infected by COVID-19 was 3.7% [3.4%-4.0%] on 01 June 2020. However, this low national average masked a stark heterogeneity across the states (Table 1). New York and New Jersey had the highest estimated cumulative attack rates, of 15.9% [12.4%-19.9%] and 14.8% [11.2%-18.2%] respectively, and Connecticut and Massachusetts both had cumulative attack rates over 10%. Conversely, other states that have drawn attention for early outbreaks, such as California, Washington, and Florida, only had cumulative attack rates of around 2% and states that were only in the early stages of their epidemics, like Maine, had estimated cumulative attack rates of less than 1%.
Figure 1 shows the effective number of infectious individuals and the number of newly infected individuals on any given day up until 01 June 2020 for each of the 8 regions in our model, which are based on US census regions (see Appendix A for further descriptions of our groupings). The effective number of infectious individuals is calculated using the generation time distribution, where individuals are weighted by how infectious they are over time, see Section 5.2 for more information. The fully infectious average includes asymptotic and symptomatic individuals. On 01 June 2020, we estimate that there were 41,100 [34,500-46,800] infectious individuals across the US, which corresponds to 0.01% of the population. Table 1 shows estimates of the number of new infections across each states on 01 June 2020. By this date, the estimated number infections were beginning to increase in the Pacific (Alaska, California, Hawaii, Oregon and Washington) and Mountain (Arizona, Colorado, Idaho, Montana, Nevada, New Mexico, Utah and Wyoming) regions.
Our model includes a state-level parameter for the infection ascertainment ratio, IAR, which we define as the number of reported cases divided by the true number of infections (including asymptomatic infections). We only estimate this parameter from 11 May 2020 when more than 375,000 tests are done each day, see Appendix B for further information. Column 3 of Table 1 shows the value of the infection ascertainment ratio in our model (see Section 5) and varies significantly between state. We would not expect the infection ascertainment to be 100% because our model includes asymptomatic individuals who may not know they have COVID-19. The mean value of this ratio varies between 43% (Missouri) to 74% (Kansas and Tennessee), which suggests that states are doing very different levels of testing.
2.2 Reproduction number
The mean estimate for Rt was below one in 23 states1 on 01 June 2020 and the 95% credible intervals did not exclude one in any state (see Appendix C for Rts by state). Figure 2 depicts the geographical variation in the posterior probability that Rt was less than 1. The closer a value is to 100%, the more certain we were that the reproduction number was below 1, indicating that new infections were not increasing. The probability was less than 40% that Rt < 1 in 20 states. There was substantial geographical clustering; most states in the Midwest and the South had reproduction numbers that suggested that the epidemic was not under control. We include figures of Rt, infections and deaths over time for each state in Appendix D.
2.3 Model effect sizes
We find that decreases in the overall average number of visits to different places had a significant effect on reducing transmission. If mobility stopped entirely (100% reduction in average mobility) then Rt would be reduced by 55.1% [26.5% - 77.0%]. The country effect size estimates are given in Figure 3, with regional and state-level effects given in Appendix E. However, in the US, the average mobility covariate never approached a 100% reduction, and only about half the states had reductions below 50% of the baseline. We define the baseline as the pre-epidemic mobility for each state [4]. As an example, consider the largest reduction observed, -62% of the baseline (Minnesota on 12 April 2020). The effect on Rt was a reduction of 37% [16% - 56%] from the country level effect.
Increased time spent in residences also reduced transmission; if time spent in residences increased to 100% of the baseline, Rt would be reduced by 15.3% [-27.5% - 54.6%]. Time spent in residences increased by 20% or more from the baseline in 36 states. As an example, consider the largest reduction observed, a 33% increase from the baseline (New Jersey on 10 April 2020). The effect on Rt from this was a reduction of 5% [-10% - 20%] in New Jersey from the country level effect.
Average mobility and residential mobility are no doubt correlated—when people spend less time in public spaces, captured by our average mobility metric, they conversely spend more time at home. Due to this collinearity, our model is unable to distinguish between the independent contributions of these covariates, with most of the effect assigned to the average mobility coefficient, due to its greater explanatory power. As a check that our overall findings were not biased by this collinearity, we verified that the posterior estimates of these coefficients was not correlated.
2.4 Short term forecasts
We used our model to produce short term death forecasts. Figure 4 compares our forecasts for the three weeks after 01 June 2020 (blue line with shaded uncertainty intervals) with the recorded daily number of deaths during this period (coral bars). As expected from our Rt values, deaths were noticeably declining in the Northeastern Corridor, where Rt > 1, with particularly low error between our forecasts and reality in New York and Connecticut. In the South, we forecast a flattening or slight increase of deaths, especially in Arkansas, Texas and Florida.
We investigated the numerical accuracy of our forecast using three metrics: mean absolute error, continuous ranked probability score (CRPS) and coverage of credible intervals. We fitted our model to three end points: 1 May, 15 May and 1 June and performed three week forecasts from each end point. We compared the metric scores with a log-linear “null” model fit to 31 days of data prior to the three specified end points (see Appendix F for further information). We find our model performs similarly to the null model (1 June) or better (15 May), however our model fit to 1 May is worse than the null model because we only include cases after 11 May in our models. This suggests that including cases improves the forecasting ability of our model and further justifies our inclusion of them. The coverage of our credible intervals is good for all models, in particular our model and the null model fit to 1 June.
2.5 Model selection and sensitivity
Mobility data provided a proxy for the behavioural changes that occur in response to non-pharmaceutical interventions. Appendix H shows the mobility trends for the 50 states and the District of Columbia up until 01 June 2020 (see Section 4 for a description of the mobility dimensions). The median correlation between the observed average mobility and the timing of the introduction of major NPIs (represented as step functions) was approximately 86% (see Appendix I). We make no explicit causal link between NPIs and mobility because this relationship is plausibly causally linked by other factors. The mobility trends data suggests that substantial early outbreak in New York state may have led to substantial changes in mobility in nearby states, like Connecticut, prior to any mandated interventions in those states, which supports including regions in our model. Including both mobility trends and the timing of imposition and lifting of “stay-at-home” orders did not affect the estimated cumulative attack rates (see Appendix J).
Mobility alone cannot fully capture how transmission evolves over time. In particular, it cannot capture the impact of case-based interventions (such as testing and tracing) or behaviour changes (such as mask wearing or hand washing). We use a second-order, weekly, autoregressive process to allow our changes in transmission to be decoupled from mobility. This autoregressive process is an additional term in our parametric equation for Rt and accounts for residual effects by capturing a correlation structure where current Rt is correlated with previous weeks Rt. This means that our forecasts were equally good whatever combination of mobility covariates were used because this term could capture the unexplained behaviour. The learnt random effects from this process are shown in Appendix K for all states. We show the contributions to Rt from the mobility and autoregressive terms in Figure 19 for three example states. The autoregressive term increases Rt before lockdown in New York, which could be explained by behaviour such as panic buying. In contrast, the autoregressive term reduces Rt in Montana and could reflect behavioural changes such as hand-washing and self isolation, which can reduce transmission with maintained mobility levels. The autoregressive term remains mostly constant in Washington and suggests that mobility is sufficient to capture the behaviour there.
3 Discussion
We developed a Bayesian semi-mechanistic modelling approach to investigate the impact of NPIs on the spread of SARS-CoV-2 in the United States through changes in mobility. Our model relies on death data from the start of the epidemic and recently reported case data to inform our predictions. This enabled us to estimate a realistic infection ascertainment ratio for the three weeks before 01 June 2020 for each state, which could help inform policy as to where testing may be lacking. The mean value of this ratio varies between 43% (Missouri) to 74% (Kansas and Tennessee). Our epidemiological grounded mechanistic model links unobserved infections to reported cases and deaths, all within a principled Bayesian statistical framework. This is a significant advancement over curve-fitting models fit directly to reported cases.
Our model suggests that although initial reductions in the daily infections had plateaued in most states by 01 June 2020, the reservoir of infectious individuals still remained large with approximately 0.01% of the population being infectious on that date. Despite this, the cumulative attack rate across the US still remained low. We found our attack rate for New York was in line with those from recent serological studies [9]. There is now evidence that mild infection is able to lead to robust immunity (via T cells) but potentially not induce antibody production, which are detected in serosurveys [10]. Therefore, serosurveys might underestimate exposure, particularly in mild cases, and our model may provide an alternative way to measure population exposure. Our cumulative attack rates are, however, sensitive to the assumed values of infection fatality rate (IFR). We account for each individual state’s age structure, and further adjust for contact mixing patterns [11], but age specific modelling may be necessary to capture potential changes in the demographics of cases in states such as Texas, Florida and South Carolina where there is evidence that younger people than were infected at the start of the epidemic are being infected [12, 13].
We estimated that 23 states had a posterior mean reproduction number Rt below one on 01 June 2020 and in no states were we more than 95% confident that Rt was below one. We compared our estimates with predictions made by rt.live [14] who use a method that fits the most likely Rt curve to the daily new daily cases (see Appendix L). Overall, our estimates were weakly correlated (ρ = 0.42) with both of us estimating Rt > 1 in 23 states (red points) including Montana and Alaska. However, the rt.live estimates are slightly more pessimistic because they predict Rt > 1 in 10 states where we predict Rt < 1 (blue points). In contrast we predict Rt > 1 in 5 states where they predict Rt < 1 (green points). Both sets of reproduction numbers strongly implied that the US epidemic was not under control in many states, and that in the presence of continued migration and the loosening of interventions seen in June, increased infections were to be expected with high probability. We found that state with high reproduction numbers on 01 June 2020 were geographically clustered in the west and south US, whilst the states that had suffered high COVID-19 mortality (such as the Northeast Corridor) in the early phase of the epidemic had lower reproduction numbers. After the period covered by this study, reported cases began to increase in the US, and 7 states (Arizona, Arkansas, California, North Carolina, South Carolina, Tennessee and Texas) had recorded higher levels of hospitalisations in early July than before [15, 16]. This suggests our estimates that Rt was not less than one were accurate. More recent estimates of Rt, the number of infections, and the number of people currently infectious are presented on our website2.
Our three week forecasts of daily deaths were highly accurate, confirming the predictive validity of our modelling approach, despite our having kept mobility constant during our forecasts. These forecasts, alongside our Rt values, show that the epidemic was not under control at the end of May. The accuracy of our forecasts varied during the epidemic and could be due to our assumption that mobility is kept constant over these three weeks. Our forecast would perform worse in weeks where mobility was significantly different to the last week of our model fit. When we include cases in our model, we are able to get similar results to a simple “null” model whilst also being about to estimate effect sizes of different mobility trends. We also compared our cumulative death forecasts with those presented by Friedman et al. [17]. Friedman et al. compared the median absolute percentage error (MAPE) for SEIR and dynamic growth rate types of models for models fit to some point in June. Unlike those models, we find the MAPE of our cumulative death forecasts did not increase significantly over time and our 3 week median cumulative death MAPE across all states (9.9%) was similar to the US estimate from Friedman et al. (4.1-8.6), see Appendix G for more information.
Our model uses mobility to predict SARS-COV2 transmission. We find that the timings that non-pharmaceutical interventions were implemented was strongly correlated to changes in mobility. This is similar to findings in Abouk and Heydari [5] who find that statewide stay-at-home orders had the strongest causal impact on reducing social interaction and that these orders significantly increase the presence of individuals at home by about six fold (our “residential mobility trend”). This supports our choice of using mobility instead of the timings of NPIs in this study instead of the times of interventions as in Flaxman et al. [6]. We find that magnitude of the reductions in average mobility, and the resulting increases in residential mobility, are important in determining the size of reduction in Rt. This agrees with Wang et al. [18] who use a stochastic age- and risk-structured susceptible-exposed-asymptomatic-symptomatic-hospitalized-recovered (SEAYHR) model to considered the effect of various levels of social distancing.
They found that social distancing measures, which reduced non-household contacts by <50%, would not prevent a healthcare crisis and that only their 75% and 90% contact reduction scenarios were projected to enable metropolitan areas to remain within health care levels.
While mobility, or social distancing measures, will explain a large amount of the trend in Rt, there is likely to be substantial residual variation from other behavioural changes such as mask wearing and hand-washing. We accounted for this residual variation through a second order, weekly, autoregressive process. This stochastic process captures changes in Rt reflected in the data, but is unable to attribute these changes to other determinants of transmission or interventions. We pool parameters in our model to leverage as much signal in our data as possible and to reflect the conjoined nature of some states, in particular in the Northeastern Corridor. While this sharing can potentially lead to over or under estimation of effect sizes, it also means that a consistent signal for all states can be estimated before that signal is presented in an individual state with little data, such as Alaska and Hawaii. Pooling also increases the robustness of our models to under reporting and time lags [6, 7, 8].
4 Data
Our model uses daily real-time state-level aggregated data published by New York Times (NYT) [19] for New York State and John Hopkins University (JHU) [2] for the remaining states. We include 105,006 deaths in our model up until 1 June and 479,422 cases from 11 May to 1 June. Age specific population counts were drawn from the U.S. Census Bureau in 2018 [20] to estimate state-specific infection fatality ratio reflective of the population age structure. The timing of NPIs were collated by the University of Washington [21]. We used Google’s COVID-19 Community Mobility Report [4], which provides data on movement in the US by states and highlights the percent change in visits to:
Grocery & pharmacy: Mobility trends for places like grocery markets, food warehouses, farmers markets, speciality food shops, drug stores, and pharmacies.
Parks: Mobility trends for places like local parks, national parks, public beaches, marinas, dog parks, plazas, and public gardens.
Transit stations: Mobility trends for places like public transport hubs such as subway, bus, and train stations.
Retail & recreation: Mobility trends for places like restaurants, cafes, shopping centres, theme parks, museums, libraries, and movie theatres.
Residential:Mobility trends for places of residence.
Workplaces: Mobility trends for places of work.
The residential data includes length of stay at different places compared to a baseline, whereas the other mobility trends are based on number of visits to a certain place. These trends are therefore relative, i.e. mobility of -20% means that, compared to normal circumstances individuals are engaging in a given activity 20% less.
5 Methods
Flaxman et al. [6] introduced a Bayesian model for estimating the transmission intensity and attack rate (percentage of the population that has been infected) from COVID-19 from the reported number of deaths. This framework used the time-varying reproduction number Rt to inform a latent function for infections, and then these infections, together with probabilistic lags, were calibrated against observed deaths. Observed deaths, while still susceptible to under reporting and delays, comprise a more consistent time series than the reported number of confirmed cases, which are susceptible to changes in the probability of ascertainment over the course of the epidemic as testing strategies changed. Our model code is available on GitHub3.
We adapted the original Bayesian semi-mechanistic model of the infection cycle to all the states in the US and the District of Columbia to infer the reproduction number over time (Rt), plausible upper and lower bounds (95% Bayesian credible intervals) of the total populations infected (attack rates) and the number of people currently infected on 01 June 2020. In this paper, we also include the reported number of cases after 11 May 2020, see Appendix B. This reflects the point in time when over 375,000 tests were being done each day across the US. We include this in our likelihood but do not use them to calculate transmission directly. We parametrise Rt as a function of Google mobility data and include an autoregressive term to capture non-mobility driven behaviour. We fit our model jointly to COVID-19 data from all states to assess the attack rates and number of people who were currently infected. Finally, we use our model to forecast for 3 weeks from 01 June 2020 and compare our estimates of deaths to those recorded in the US for each state. We assume mobility remains constant at the previous value of mobility on the same day the previous week in our forecasts and the autoregressive term remains constant.
5.1 Model specifics
The true number of infected individuals, i, is modelled using a discrete renewal process. We specify a generation distribution g with density g(τ) as: Given the generation distribution, the number of infections it,m on a given day t, and state m, is given by the following discrete convolution function: where the generation distribution is discretised by for s = 2, 3, …, and . The population of state m is denoted by Nm. We include the adjustment factor St,m to account for the number of susceptible individuals left in the population.
Both deaths and cases are observed in our model. We define daily deaths, Dt,m, for days t ∈ {1, …, n} and states m ∈ {1, …, M}. These daily deaths are modelled using a positive real-valued function dt,m = 𝔼[Dt,m] that represents the expected number of deaths attributed to COVID-19. The daily deaths Dt,m are assumed to follow a negative binomial distribution with mean dt,m and variance , where ψ1 follows a positive half normal distribution, i.e. Here, 𝒩 (µ, σ) denotes a normal distribution with mean µ and standard deviation σ. We say that X follows a positive half normal distribution 𝒩 +(0, σ) if X ∼ |Y |, where Y ∼ 𝒩 (0, σ).
We link our observed deaths mechanistically to transmission as in Flaxman et al. [6]. We use a previously estimated COVID-19 infection fatality ratio (IFR, probability of death given infection) together with a distribution of times from infection to death π. Details of this calculation can be found in [22, 23]. From the above, every region has a specific mean infection fatality ratio ifrm (see Appendix M). To incorporate the uncertainty inherent in this estimate we allow the ifrm for every state to have additional noise around the mean. Specifically we assume We believe a large-scale contact survey similar to polymod [11] has not been collated for the USA, so we assume the contact patterns are similar to those in the UK. We conducted a sensitivity analysis, shown in Appendix M, and found that the IFR calculated using the contact matrices of other European countries lay within the posterior of .
Using estimated epidemiological information from previous studies, we assume the distribution of times from infection to death π (infection-to-death) to be the convolution of an infection-to-onset distribution (π ′) [23] and an onset-to-death distribution [22]: The expected number of deaths dt,m, on a given day t, for state m is given by the following discrete sum: where iτ,m is the number of new infections on day τ in state m and where, similar to the generation distribution, π is discretized via for s = 2, 3, …, and , where π(τ) is the density of π.
For every state m, we also observe daily cases Ct,m after tc = 11 May 2020. Similar to daily deaths, daily cases are modelled using a positive real-valued function that represents the expected number of symptomatic and asymptomatic cases. Again, the daily cases Ct,m are assumed to follow a negative binomial distribution but with mean ct,m and variance , where ψ2 follows a positive half normal distribution, i.e. As before, we assume the distribution of times from infection to becoming a case π (infection-to-onset) to be We add in a new link between our observed daily cases and our estimated daily infections. We use our model to estimate an infection ascertainment ratio (iarm) for each state m, which is defined as the number of reported cases divided by the true number of infections (including both symtomatic and asymptomatic infections). This follows a Beta distribution, specifically um ∼ Beta(12, 5).
The expected number of cases ct,m, on a given day t, for state m is given by the following discrete sum: where, again, cτ,m is the number of new infections on day τ in state m and where π ′ is discretized via for s = 2, 3, …, and , where π′(τ) is the density of π ′.
We parametrise Rt,m as a linear function of the relative change in time spent and number of visits (from a baseline) where f (x) = 2 exp(x)/(1 + exp(x)) is twice the inverse logit function. Xt,m,k are covariates that have the same effect for all states, Yt,m,l is a covariate that has region-specific effects, r(m) ∈ {1, …, R} is the region a state is in (see Figure 15), Zt,m is a covariate that has a state-specific effect and is a weekly AR(2) process, centred around 0, that captures variation between states that is not explained by the covariates.
The prior distribution for R0,m[24] was chosen to be where κ is the same among all states. In the analysis of this paper we chose the following covariates: (an intercept), and , where the mobility variables are from [4] and defined as follows (all are encoded so that 0 is the baseline and 1 is a full reduction of the mobility in this dimension):
is an average of retail and recreation, groceries and pharmacies, and workplaces. An average is taken as these dimensions are strongly collinear.
are the mobility trends for places of residences.
We include regional, as well as state level parameters, in our model to encapsulate the connected nature of states. This was particularly important in the Northeasten corridor where residents in New Jersey and Connecticut regularly commuted into New York, the early epicenter of the US epidemic (see Appendix A for a map of the regions). Regions are based on US Census Divisions, modified to account for coordination between groups of state governments [25].
We assume that seeding of new infections begins 30 days before the day after a state has cumulatively observed 10 deaths. From this date, we seed our model with 6 sequential days of an equal number of infections: , where τ ∼ Exponential(0.03). These seed infections are inferred in our Bayesian posterior distribution.
The weekly, state-specific effect is modelled as a weekly AR(2) process, centred around 0 with stationary standard deviation σw that, in every state, starts on the first day of its seeding of infections, i.e. 30 days before a total of 10 cumulative deaths have been observed in this state. The AR(2) process starts with , with independent priors on ρ1 and ρ2 that are normal distributions conditioned to be in [0, 1]; the prior for ρ1 is a 𝒩 (0.8, .05) distribution conditioned to be in [0, 1] and the prior for ρ2 is a 𝒩 (0.1, .05) distribution, conditioned to be in [0, 1]. The prior for σw, the standard deviation of the stationary distribution of 𝔼w is chosen as σw ∼ 𝒩 +(0, .2). The standard deviation of the weekly updates to achieve this standard deviation of the stationary distribution is . The conversion from days to weeks is encoded in wm(t). Every 7 days, wm is incremented, i.e. we set , where is the first day of seeding. We keep the AR(2) process constant over the last 7 days of observations since this is less informed by data due to the lags and also over the forecast period.
The prior distribution for the shared coefficients were chosen to be and the prior distribution for the pooled coefficients were chosen to be We estimated parameters jointly for all states in a single hierarchical model. Fitting was done in the probabilistic programming language Stan[26] using an adaptive Hamiltonian Monte Carlo (HMC) sampler.
5.2 Generated quantities
The effective number of infectious individuals, i*, on a given day considers how infectious a previously infected individual is on a given day and includes both asymptotic and symptomatic individuals. It is calculated by first re-scaling the generation distribution by its maximum, i.e. . Based on (2), the number of infectious individuals is then calculated from the number of previously infected individuals, c, using the following: where it,m is the number of new infections on day t in state m.
Data Availability
All source code and data necessary for the replication of our results and figures is available at https://github.com/ImperialCollegeLondon/covid19model.
7 Author Contributions
HJTU, SM, VB, AG, SB and SF conceived and designed the study. HJTU, SM, TM, MACV, HZ and PM performed mobility analysis. JI-H, SLF and XX contributed to statistical analysis and MH and IH did other analysis. HJTU, SM, AG, TM, HC, MACV, VW, MM, OR, SB and SF contributed to code development. HJTU, MACV and SF did the plotting. SM and FV created the website. HJTU, VB, SB and SF wrote the first draft of the paper. All authors discussed the results and contributed to the revision of the final manuscript.
8 Data and Code availability
All source code and data necessary for the replication of our results and figures is available at https://github.com/ImperialCollegeLondon/covid19model.
9 Competing Interests statement
The authors declare no competing interests.
A Model regions
We choose 8 regions in our model, see Figure 5. These regions are based on US Census Divisions, modified to account for coordination between groups of state, see Section 5 for more information.
B Ratio of reported cases to estimated infections
We compare the weekly moving average of the percentage of reported cases to mean estimated infections for the time period of the model. We use the same model formulation without cases here to illustrate two points. First, the percentage is small initially, which is why we do not use cases in the model until 11 May 2020. Second, sufficient tests are being done towards the end of May so they are informative to our model and we can estimate an infection ascertainment ratio, see Section 5 for further details. Percentages may be greater than 100% when our model estimates fewer cases than were reported. This occurs in states with little data.
C State-level Rt values
Figure 7 shows the value of Rt at two points during the epidemic. The first time point is for the week centered on the date on which each state implemented any emergency decree order such as State of Emergency, Public Health Emergency, and Public Health Disaster declarations. The second is for the week ending on 01 June 2020.
D Model predictions for all states
State-level estimates of infections, deaths and Rt. Left: daily number of deaths, brown bars are reported deaths, blue bands are predicted deaths, dark blue 50% credible interval (CI), light blue 95% CI. Middle: daily number of infections, brown bars are reported infections, blue bands are predicted infections, CIs are same as left. The number of daily infections estimated by our model drops immediately after an intervention, as we assume that all infected people become immediately less infectious through the intervention. Afterwards, if the Rt is above 1, the number of infections will start growing again. Right: time-varying reproduction number Rt dark green 50% CI, light green 95% CI. Icons are interventions shown at the time they occurred.
E Effect sizes
We plot estimates of the posterior mean effect sizes and 95% credible intervals for each mobility category. The relative % reduction in Rt metric is interpreted as follows: the larger the percentage, the more Rt decreases, meaning the disease spreads less; a 100% relative reduction ends disease transmissibility entirely. The smaller the percentage, the less effect the covariate has on transmissibility. A 0% relative reduction has no effect on Rt and thus no effect on the transmissibility of the disease, while a negative percent reduction implies an increase in transmissibility.
F Forecast evaluation
We evaluate three-week ahead model forecasts (days t ∈ {T + 1, …, T + 21}) for our model using two metrics. The first metric is the mean absolute error (MAE), which is given by where are S posterior predictive samples of daily deaths on day t in country m and yt,m is the observed value of the corresponding quantity. The continuous ranked probability score (CRPS) is a generalisation of MAE to probabilistic forecasts and can be estimated using Gneiting and Raftery [27] Figure 11 shows the mean MAE and CRPS for three different forecast periods across all time points and states. In addition to the forecast results shown in Section 2.4, we also fit our model up to 1 May and 15 May and show the error in our 21 day forecasts. We compare these forecasts to a simple log-linear “null” model fit to the 31 days of death data prior to the end dates of our model runs (1 May, 15 May and 1 June) to evaluate our predictions. We find the MAE and CRPS error were different for models fit up to different dates. Our model fit to 1 May performs worse than the “null” model because this model only includes deaths (we only include cases in our model after 11 May, see Appendix B. Our model performs similarly (1 June) or better than the null model (15 May) in both MAE and CRPS. We include a state-level break down of MAE in Figure 12 and CRPS in Figure 13.
The forecasts were also evaluated using the mean coverage of their 95% and 50% credible intervals. If model uncertainty is well-calibrated then the observed quantities will fall outside of the 95% credible intervals 5% of the time and 50% of the time for the 50% credible intervals. Recent work has highlighted that other prominent models do not meet this criterion, which suggests that their uncertainty estimates are not well-calibrated [28]. The mean coverage of the 95% credible interval in a time period starting at time t0 with length L is given by where 1(·) is the indicator function and is the z-th percentile of the samples on day t in country m. The mean coverage of the 50% credible interval is The coverage was calculated for each state individually and then the mean across all the 50 state and the District of
Columbia was computed. Again these coverage probabilities were compared with the log-linear “null” model.
These coverage results suggest that both the 95 and 50% credible intervals are well-calibrated for all the models but best for our model and the null model fit to the 1 June. We hypothesise this is because they include more in data in the forecast.
G Median absolute percent error
We compared the magnitude of our cumulative death median absolute percent error (MAPE) estimates from our forecasts from 1 June given in Table 2 with global COVID-19 forecasts found in Friedman et al. [17]. This is not a parameter we explicitly wanted to forecast but enabled us to compare with other models in literature. Friedman et al.’s paper included both SEIR type models (IHME – MS SEIR [29], Youyang Gu [30], MIT - DELPHI [31] and Imperial-LMIC [32]) and dynamic growth rate models (LANL-Growthrate [33]). The forecasts in Friedman et al. [17] were not all fit to the same date range, but to some point in June and are only presented for the whole of the US and not to each state specifically. Unlike the models compared by Friedman et al., the MAPE of our cumulative death forecasts did not increase significantly over time. Our 3 week median cumulative death MAPE across all states (9.9%) was similar to the US estimate from Friedman et al. (4.1-8.6 excluding the Imperial-LMIC model which was not calibrated for high income settings and never meant to be used for the US), but the magnitude was slightly higher in our forecasts for some states, in particular Alaska and Wyoming, which were the last two states to reach 10 cumulative deaths and so had the least data in our models.
H Mobility trends
In Figure 15 we show the Google mobility trends [4] across the 50 states of the US and the District of Columbia. In our model, we use the time spent at one’s residence and the number of visits to grocery stores, pharmacies, recreation centres, and workplaces. Our regions are based on US Census Divisions, modified to account for coordination between groups of state governments [25]. These trends are relative to a state-dependent baseline, which was calculated shortly before the COVID-19 epidemic. For example, a value of -20% in the average trend means that individuals, on average, are visiting 20% less shops, recreation places and workplaces than before the epidemic. We overlay the timing of two major state-wide NPIs (stay at home and emergency decree) on Figure 15 (see [21] for details). We note intuitive changes in mobility such as the spike on 11th and 12th April for Easter.
I Mobility regression analysis
In Figure 16 we regressed NPIs against average mobility. We parameterised NPIs as piece-wise constant functions that are zero when the intervention has not been implemented and one when it has. We evaluated the correlation between the predictions from the linear model and the actual average mobility. We also lagged the timing of interventions and investigate its impact on predicted correlation. We observed reduced correlation when lagging (forward and backwards) the timing of NPIs, which suggested an immediate impact on mobility.
J Comparison of model attack rates with and without interventions in the model
We show in Figure 17 that there is little difference between the attack rates for two models where we include the “Stay at Home” intervention in our model as well and when we just use mobility. We, therefore, choose the mobility only model to limit the numbers of parameters we are fitting.
K State-specific weekly effects
Our model includes a state-specific weekly effect Ew,m (see equations 12, 14) for each week w in the epidemic period for a state. As described in Section 5, we assign an autoregressive process with mean 0 as prior to this effect. Figure 20 shows the posterior of this effect on the same scale as in Figure 3, that is, the percent reduction in Rt with mobility variables held constant4. Values above 0 have the interpretation that the state-specific weekly effect lowers the reproduction number Rt,m, i.e. transmission for week t and state m is lower than what is explained by the mobility covariates.
In Figure 19, we compare the different terms contributing to Rt, as explained by equation (12), for three example states. We choose New York, Montana and Washington for our vignettes because the autoregressive (purple) trend behave differently in these states. In New York, the autoregressive term increases Rt before lockdown, which we hypothesise corresponds to importations and new seeding events from Europe driving transmission. Then the autore-gressive term decreases as social distancing, mask wearing and handwashing are implemented along side behavioural changes. In contrast, the autoregressive term reduces Rt at the start in Montana and the mobility trends remain very flat, adding some weight to the hypothesis of importations. This could reflect the low density of that state and the correspondingly long distances between where people live, along with behavioural changes. The autoregressive term remains mostly constant in Washington, in line with the states early and effective epidemic response, and suggests that mobility is sufficient to capture the behaviour there.
L Comparison of Rt estimations
We compare our Rt estimates on 01 June 2020 with estimates made by rt.live [14]. Overall, our estimates were weakly correlated (ρ = 0.42)with both of us estimating Rt > 1 in 23 states (red points) including Montana and Alaska. However, the rt.live estimates are slightly more pessimistic because they predict Rt > 1 in 10 states where we predict Rt < 1 (blue points). In contrast we predict Rt > 1 in 5 states where they predict Rt < 1 (green points).
M Sensitivity analysis to infection fatality ratio
Geographic-specific contact surveys are important for calculating the weighted IFR values according to the methods in [22, 23]. There is no large-scale cross-generational contact survey, similar to the polymod survey [11], implemented in the USA. Therefore, it was important to understand if the model was robust to changes in the underlying contact survey. We calculated the IFRs using three different contact matrices: UK, France and Netherlands. We believe that the USA is culturally closest to that UK out of those countries we had contact matrices for, but also considered France where we saw the greatest mixing of the elderly and the Netherlands which showed the average behaviour of the European studies used in [23]. We found that the IFR, calculated for each state using the three contact matrices, lay within the posterior of IFR in our model (Figure 21). We also noted that our results remained approximately constant when using the IFR calculated from the three different contact matrices as the mean of the prior IFR in out model, see Section 5.
Since we are using the same contact matrix across all the states, the differences in IFR are due to the population demographics and not due to differential contacts. The low IFR in Texas and Utah reflects the younger population there whereas the higher IFR in Florida and Maine is due to the older population. This is a limitation of our methods.
6 Acknowledgements
We would like to thank Amazon AWS and Microsoft Azure for computational credits and we would like to thank the Stan development team for their ongoing assistance. We would also like to thank David Joerg and Jacob Steinhardt for their comments through Open Review. This work was supported by Centre funding from the UK Medical Research Council under a concordat with the UK Department for International Development, the NIHR Health Protection Research Unit in Modelling Methodology and Community Jameel. HJTU is funded by Imperial College London through an Imperial College Research Fellowship grant. SB would like to acknowledge the NIHR BRC Imperial College NHS Trust Infection and COVID themes the Academy of Medical sciences Springboard award and the Bill and Melinda Gates Foundation.