Rapid spread of SARS-CoV-2 virus across the planet has precipitated a global 2 pandemic, infecting millions and killing hundreds of thousands. Governments around 3 the world have undertaken unprecedented interventions aimed at curtailing the spread 4 and lethality of the virus. These interventions and more recently proposals for relaxing 5 them have relied heavily on predictions of COVID-19 case growth and mortality. 6 COVID-19 prediction models can be roughly classified as mathematical models that 7 simulate disease within a population or statistical models that fit a function to observed 8 data and extrapolate forward into the future. We will discuss the features of both types 9 of models. Most COVID-19 models are compartmental models , a type of 10 mathematical model used by epidemiologists to simulate infectious disease epidemics for 11 over a century. Compartmental models divide a population into mutually exclusive 12 compartments that denote disease status and supply a set of differential equations that 13 define the flow of the population between compartments [62] . Traditionally they are 14 named after their compartments with the SIR (susceptible-infectious-recovered) [63] and 15 SEIR (susceptible-exposed-infectious-recovered) models classic examples. 16 In an infectious disease compartmental model, S(t) is the number of susceptible 17 individuals at time t, and new infections are represented by the flow of individuals out 18 of the S compartment. This is governed by the first derivative of S(t) with respect to 19 time, dS(t)/dt. Classic SIR and SEIR models express this as proportional to the 
where I(t) is the number of infectious individuals at time t. The rate β is often 22 interpreted as disease transmissibility and may be expressed as a function of the 23 reproductive number R 0 -the expected number of individuals infected by an infectious 24 person-and contact rates between individuals. It may also be normalized in Eq 1 by 25 division with the total population size. 26 The simplest approach for simulating infections is to assume a value for β or its 27 constituent parts from the literature or other prior information [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] . While this is 28 convenient, the predictive accuracy can suffer. Another approach that has been used by 29 other studies is to estimate β (or a related quantity) by fitting a statistical model or 30 other optimization procedure to observed data . This empirical approach can 31 make these models more realistic, but they still may be limited in their ability to 32 accurately model the COVID-19 pandemic. Disease transmission rates in COVID-19 33 have changed substantially over time depending upon the political and societal responses 34 and possibly other factors [54] . As a result, modelers operating within this framework network that models time-varying transmission by estimating intervention efficiency 48 from reported data as a function of time [42] . 49 Recovery, death, and other states (e.g., hospitalization) may be incorporated into the 50 model as separate compartments. Solutions to the differential system provide values for 51 each compartment at each time, allowing for easy joint modelling of disease states once 52 their derivative is specified. This is an advantage of compartmental models over many 53 other approaches, which may require separate models for each quantity. 54 A number of agent-based COVID-19 models have been developed or adapted from 55 influenza pandemic models to simulate the individuals of a population and their 56 interactions [64] [65] [66] [67] [68] . This provides a mechanism for modelling interventions that target 57 contacts between individuals and does not assume the population exists in homogeneous 58 compartments as compartmental models generally do, but also requires a number of 59 assumptions to be made on the behavior and interactions within a population as well as 60 the infectivity of COVID-19. 61 Serial growth models for COVID-19 simulate an epidemic by expressing the number 62 of new infections at a given time as a weighted sum of new infections on previous days 63 usually scaled by the reproductive number, which may be time-varying [69] [70] [71] [72] [73] . The 64 weights are sampled from a probability distribution defining the amount of time between 65 an individual being infected and infecting another person. Deaths or other outcomes 66 may be modelled as a second step, for example using a negative binomial model that 67 predicts daily deaths conditional on the number of infections in recent days [70] . 68 Statistical models often eschew deterministic population dynamics and fit the 69 observed data as a function of time and possibly other covariates in a regression (or 70 equivalent) framework. Log-linear [74] , generalized Richards [75] , ARIMA [76, 77] , 71 exponential [78] , Gaussian CDF [79] , and logistic [80] [81] [82] models, which all 72 accommodate the generally sigmoidal shape of the cumulative infection count that is 73 often observed in epidemics, as well as various other models [83] [84] [85] [86] including machine 74 learning algorithms [87] [88] [89] have been proposed for COVID-19. Murray et al. and 75 Woody et al. take similar approaches for modeling COVID-19 deaths using the error 76 function (ERF) [90, 91] . Modeling deaths is appealing, because COVID-19 deaths have 77 been more reliably reported than infections. However, because deaths lag infections by 78 some amount of time, it may not enable projections to incorporate the latest 79 information on disease spread. 80 Within the framework of a statistical (or other regression-like) model, it is easier to 81 fit observed data, assuming an appropriate functional form is selected, but it may be 82 challenging to accurately project the future trajectory of an epidemic. Time-varying 83 covariates like mobile phone tracking data [91] , Google trends [89, 92] , and social 84 media [93] are easily incorporated into such a model and may be quite predictive, but 85 they are often unknown in the future, requiring assumptions to be made regarding their 86 future values when forecasting. Such approaches can only model one outcome (e.g., infections or deaths) and additional steps must be taken to predict other quantities. 88 Here we present an approach to projecting COVID-19 cases and deaths that employs providing location-specific trajectories that extrapolate well within a full probability 92 model that includes uncertainty quantification. We use a random forest algorithm for 93 the death transition function that learns the relationship between COVID-19 cases and 94 population characteristics to predict deaths. The remaining sections of this paper lay 95 out the SIRD compartmental model, the Bayesian mixed model for case velocity, the 96 random forest death model, and close with results and a discussion. 97 Materials and methods 98 The SIRD Compartmental Model 99 We model the spread and progression of COVID-19 through a population using a SIRD compartmental model named after the four compartments into which it partitions the population: S for susceptible, I for infectious, R for recovered, and D for dead. The number of population members in each compartment is a function of time, t, and these functions are linked by a system of ordinary differential equations (ODEs) that govern the flow of the population through the different disease states: will not substantially alter inference in the short term.
The transition rates between compartments are determined by the functional forms 105 and parameter values in Eq 2. Given these and initial conditions for the system, S(t 0 ), 106 May 15, 2020 4/19 I(t 0 ), R(t 0 ), and D(t 0 ), the system of ODEs in Eq 2 is deterministic, but in general 107 does not accommodate analytical solutions. Consequently, we compute numerical 108 solutions using the lsoda solver in the deSolve package [94] of R [95] .
Due to the novelty of the SARS-CoV-2 virus and a desire to empirically ground the 110 compartmental model, we fit transition functions that can vary in time and incorporate 111 covariates and other information. The transition between S and I is determined by ξ(t), 112 which describes the number of individuals becoming confirmed COVID-19 cases and is 113 presented in the next section. In many locations there is no reliable data on the 114 transition of individuals out of I into R, except for hospitalized patients in some places. 115 We follow the standard SIR approach and model transition out of I with a rate 116 parameter ρ that is the inverse of the number of days a person is expected to be 117 infected. The destination of individuals transitioning out of I is determined by the 118 death model θ(t).
Initial conditions for the model were constructed by stepping the system through the 120 observed case count data and then projecting forward. This approach should be more 121 accurate than simply beginning the simulation on the last observed day, because the 122 distribution of cases into compartments I, R, and D is not observed. were enacted, which have gradually curbed case growth. We investigate the dynamics of 128 COVID-19 case growth by modeling the velocity, i.e., the first derivative with respect to 129 time, of the log cumulative case count. This is the instantaneous rate of new cases to 130 cumulative cases at a given time, and is very similar to the reproductive number.
Calculating the reproductive number at a particular time, however, requires knowing 132 the number of active infections. There is currently no reliable data on this, as most 133 infections resolve on their own outside of a clinical or otherwise supervised setting in 134 which their transition from active case to recovered might be recorded. The velocity of 135 log cumulative cases on the other hand is readily estimated from the data and presents 136 itself as an appealing target of analysis.
A very crude estimate of the derivative can be obtained using first differences, but 138 smoothing allows for more precise estimates, as calculating the derivative requires some 139 notion of function smoothness [96] . We estimate the velocity by fitting a cubic spline to 140 the observed log cumulative case count and then evaluating its derivative at the 141 observed time points. Since there is relatively little noise in the cumulative accounts, we 142 assume any uncertainty introduced by this procedure is negligible. The velocity of a cumulative function cannot be negative, since cumulative functions 152 are monotonically increasing. Consequently, we employed a log link to map velocity to 153 the entire real line and modeled its mean with a linear predictor. We selected a 154 Bayesian mixed model to obtain location-specific estimates of the trajectory. . CC-BY-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.
The copyright holder for this preprint this version posted May 18, 2020. . https://doi.org/10.1101/2020.05.15.20102608 doi: medRxiv preprint Let u i (t) denote the cumulative case count for location i at time t, and y i (t) the first 160 derivative with respect to time of its log transformation, i.e., y i (t) = d log u i (t)/dt.
Log-transformed velocity is modeled as a linear combination of random effects and 162 Gaussian noise,
where y i (t i ) is a vector of length n i denoting the velocity evaluated at times
is a vector indicating whether each time occurred 165 after intervention, i.e., z ij = I(t ij ≥ π i ); α i is the random intercept, β i the random 166 slope, γ i the random effect for intervention, and δ i is the random slope for intervention. 167 The elements in the error vector i are assumed to be independently Gaussian 168 distributed with mean zero and precision (inverse variance) τ ij , where j = 1, ..., n i . In 169 vector notation,
where N ni is an n i dimensional Gaussian distribution, its mean 0 is a vector of zeros, 171 and its covariance matrix is diagonal with elements τ −1
Obvious heteroskedasticity is apparent in the observed case counts with variation 173 decreasing with time. To account for this and to allow for the variance to differ between 174 locations, the precision vector for location i, τ i , was modeled as linear in time with 175 location-specific random effects,
May 15, 2020 6/19
. CC-BY-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.
The copyright holder for this preprint this version posted May 18, 2020. . Each of the velocity and precision random effects were assigned Gaussian priors,
with means µ α , µ β , µ γ , µ δ , µ η , µ ζ , and variances σ 2 α , σ 2 β , σ 2 γ , σ 2 δ , σ 2 η , σ 2 ζ . These means 177 were themselves given Gaussian hyperpriors. The prior mean and variance values used 178 for the predictions presented here are listed in S2 Table of the supplemental material. 
where t ini is the final observed time point for location i. The lognormal model for
The 
where the addition of a superscript (m) to a parameter denotes its m-th posterior . CC-BY-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.
The copyright holder for this preprint this version posted May 18, 2020. classification trees into an ensemble [98] . It is a very commonly used algorithm known 219 to perform well at a variety of predictive tasks [99] . Covariate importance scores were computed using permutation variable importance. 238 Briefly, permutation importance can be thought of as the decrease in predictive 239 accuracy (in terms of mean squared error (MSE)) between the original model and when 240 each variable is randomly permuted thus obscuring any signal it has with the outcome 241 variable. The idea is that if a covariate is important in terms of prediction, obscuring its 242 signal should result in a decrease in predictive accuracy. The most important of the 243 lagged cases was that at t − 10, indicating that the model was able to discern a lag time 244 between positive tests and COVID-19 fatalities. Additional lagged cases beyond 14 days 245 were not included in the model, even though they may have been informative, because 246 each additional lagged day reduces the number of available training observations at each 247 location.
Fitting the model to data collected through April 29, 2020, resulted in a very high 249 out-of-bag R 2 of 0.90. This is an overly optimistic estimate of prediction error, due to 250 the within-location and temporal dependence of the data [101] , but more significantly 251 due to the lagged case counts being very informative covariates. Lagged cases were far 252 more important than the demographic characteristics, which is not surprising . CC-BY-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.
The copyright holder for this preprint this version posted May 18, 2020. testing was quite limited. The random forest predictions were capped at a percentage of 256 the new cases to avoid unrealistically high death predictions, which can occur when 257 there are relatively few new cases. This upper bound was set to be equivalent to a 15% 258 case fatality rate for the first 30 days of the epidemic and reduced to 7% subsequently, 259 with the higher initial death rate motivated by the relative severity of early confirmed 260 cases due to limited testing.
Predictive Accuracy
We assessed the predictive accuracy of our model by training it on case and death 263 counts collected through April 15, 2020, and predicting through April 29. We quantified 264 prediction error for each state on each day using the mean absolute scaled error (MASE) 265 of the posterior median number of new cases and deaths. MASE is computed by 266 dividing the mean absolute prediction error by the in-sample mean absolute error of a 267 naive random walk forecast,
where Y = (Y 1 , ..., Y n ) is the training data outcome, Y * = (Y * 1 , ..., Y * m ) is the observed 269 outcome in the evaluation set andŶ = (Ŷ * 1 , ...,Ŷ * m ) is the prediction for Y * to be 270 evaluated [102] . A MASE of one indicates that the predictions were on average equally 271 accurate to the mean accuracy of a random walk forecast in the training data. A useful 272 feature of MASE for our purposes is its scale invariance, which makes comparisons of 273 predictive accuracy between states with epidemics on different scales more meaningful. 274 Figure 4 depicts the distribution of MASE across states for cases and deaths over the . CC-BY-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.
The copyright holder for this preprint this version posted May 18, 2020. 
The incorporation of covariates into the case growth and mortality models allow for 295 alternate trajectories to be explored, which is the subject of ongoing research. The 296 cumulative case counts eventually plateau for each state as its case velocity decreases 297 toward zero. The predicted daily death counts are not smooth because the random 298 forest algorithm averages many discontinuous segments into a prediction.
New York, especially New York City with its large, dense population, has been the 300 epicenter of the largest COVID-19 outbreak in the United States with over 300,000 301 confirmed cases by late April. Initial exponential case growth was slowly curbed by 302 public interventions, leading to a consistent decrease in case velocity and peaks in active 303 cases and deaths in mid April. Case growth being well past its peak translates into a 304 plateauing cumulative case curve and a relatively narrow interval estimate compared to 305 states that peaked more recently or have yet to peak.
Like many other states, Ohio has had many fewer cases than New York with 307 approximately 18,000 cases and appears to be peaking near the end of April. Its interval 308 estimates are relatively wider than New York, because there is less uncertainty in the 309 estimated trajectory farther past the peak. Ohio also exhibits more relative variation in 310 its daily death counts than New York because of the smaller number.
Mississippi with fewer than 7,000 cases through the end of April illustrates the 312 estimated trajectories of a relatively rural, Southern state that has not yet peaked.
With cases farther from their plateau, there is correspondingly more relative uncertainty 314 in its trajectory and a much wider range of dates over which its peak may occur. . CC-BY-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.
The copyright holder for this preprint this version posted May 18, 2020. Cumulative Confirmed Cases q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0 5,000 Cumulative Confirmed Cases q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q The figures include 95% credible intervals around the median indicating that 95% of 316 simulation results fell within this region. These intervals are not true credible intervals 317 in the Bayesian sense, because random forest is not a probability model. Nevertheless, 318 they represent a reasonable account of model uncertainty, as they incorporate credible 319 intervals from the Bayesian case model and uncertainty around the duration of illness. 320 Despite the many strengths of the current approach, it is not without limitations.
The projections produced here assume states continue upon their current trajectories.
Changes in policy interventions, for example, could result in substantial deviation from 323 this. Projecting outcomes under different or changing intervention scenarios is the 324 subject of ongoing work.
Considering COVID-19 cases and death over large areas can obscure variation on a 326 smaller scale. It is possible for a generally positive trajectory at the state-level to mask 327 a burgeoning outbreak in some locale within the state until that outbreak contributes 328 sufficiently many cases to influence the state-wide trajectory. A more granular approach 329 that models COVID-19 at a finer resolution may be able to identify such an outbreak 330 earlier.
There is substantial interest in estimating the proportion of the population that has 332 or will have recovered from COVID-19 in the hopes that these individuals have acquired 333 at least temporary immunity to the virus and can be the vanguard to economic recovery. 334 Since we focus on modeling confirmed cases and deaths, our model does not predict the 335 true number of recovered individuals. It is well known that, especially in the U.S., 336 confirmed cases are a substantial undercount for the true number of COVID-19 337 infections. As a result, estimating the number of recovered individuals requires 338 additional information beyond predictions of confirmed cases and deaths. Attempts to 339 quantify recovery using serology testing are underway in the U.S. and elsewhere.
There is residual temporal autocorrelation in the case velocity not captured by the 341 random effects for intercept and slope in the mixed model. We expect this has minimal 342 impact on our posterior trajectories, as we are primarily concerned with the trend in 343 mean over time, but incorporating a more sophisticated approach for temporal 344 dependence could be used to explicitly model this autocorrelation.
Finally, one could consider more elegant methods for incorporating lagged case 346 counts into a death model than simply inserting them as covariates into random forest. 347 However, many approaches to lag estimation are only good retrospectively and thus are 348 insufficient for the current task. 