Hawkes process modeling of COVID-19 with mobility leading indicators and spatial covariates =========================================================================================== * Wen-Hao Chiang * Xueying Liu * George Mohler ## Abstract Hawkes processes are used in machine learning for event clustering and causal inference, while they also can be viewed as stochastic versions of popular compartmental models used in epidemiology. Here we show how to develop accurate models of COVID-19 transmission using Hawkes processes with spatial-temporal covariates. We model the conditional intensity of new COVID-19 cases and deaths in the U.S. at the county level, estimating the dynamic reproduction number of the virus within an EM algorithm through a regression on Google mobility indices and demographic covariates in the maximization step. We validate the approach on both short-term and long-term forecasting tasks, showing that the Hawkes process outperforms several benchmark models currently used to track the pandemic, including an ensemble approach and an SEIR-variant. We also investigate which covariates and mobility indices are most important for building forecasts of COVID-19 in the U.S. Keywords * COVID-19 forecasting * Hawkes processes * Mobility indices * Spatial covariate ## 1. Introduction Mathematical modeling and forecasting are playing a pivotal role in the ongoing SARS-CoV-2 (COVID-19) pandemic. In mid-March 2020, a report out of Imperial College London [1] forecasted severe consequences in the U.S. and U.K. without significant public health interventions. In both nations, governments responded by closing schools, non-essential businesses and releasing general stay-at-home (shelter-in-place) orders. In the U.S., state and local policymakers are using mathematical models and projections to inform decisions about when and how to relax public health measures that have been put in place. By and large, compartmental models that explicitly incorporate transmission characteristics of infectious diseases have been favored over machine learning approaches. High profile Susceptible-Exposed-Infected-Removed (SEIR) models include those out of Columbia University [2], MIT [3], The Johns Hopkins University [4], and UCLA [5] (in the case of the UCLA model, an SEIR-variant with an unreported compartment is fit using least-squares to reported infection and recovery data). A major exception is the well-known IHME model [6], which employs Gaussian curve fitting to COVID-19 case and death count time series in locations further along (e.g., China, Europe) to estimate curves in locations where the outbreak is more recent (e.g., the United States). The IHME model has been called into question by epidemiologists because it lacks explicit transmission dynamics in the model [7]. Our goal in this paper is to show that Hawkes processes, widely used in the machine learning community to model contagion patterns in event data, are well suited for modeling and forecasting COVID-19 case and mortality data. They have several advantages that we will highlight, including being highly flexible in accommodating auxiliary spatio-temporal features such as county-level demographics and temporal mobility patterns, yet mathematically they are connected to compartmental models [8] and allow for explicit incorporation of transmission dynamics (which we briefly review in the following section). Furthermore, extensive research has been conducted in the past several years on a couple of machine learning techniques with the point process framework. Non-parametric Hawkes processes can be constructed where the triggering kernel is learned [9] and more recently, fully neural network based point processes have been developed [10, 11]. Sparse linear combinations of Hawkes processes were a winning solution in the 2017 NIJ Crime Forecasting Challenge [12]. In certain cases a mixture of Hawkes processes may be needed to model more complex event contagion with high dimensional marks through Dirichlet processes [13, 14]. Hawkes processes can also be used for causal inference on networks [15] and recent efforts have also focused on training point processes through reinforcement learning [16, 17]. We believe all of these methods have potential applications to modeling infectious diseases that are highly complex due to heterogeneity in the population, environment, and disparate public policies across regional and local jurisdictions. Despite these advantages, to our knowledge, the only U.S. state where a Hawkes process is being used to inform COVID-19 policy is in New Jersey (a collaboration with Facebook AI Research) 1. The outline of the paper is as follows. In Section 2, we introduce our Hawkes process model whose productivity (reproduction number) is dynamic and depends on spatio-temporal covariates. Unlike recently introduced models that incorporate covariates into the background rate of a Hawkes process [18, 19], our Hawkes process model may be viewed as a convolution of lagged mobility with an inter-infection time distribution to estimate the intensity of secondary infections in the future. This is important as phased reopening in the U.S. leads to mobility changes, the effects of which are not realized in the case and mortality data until days or weeks later. Hence the model can be used to forecast changes in transmission and new cases in real-time as mobility changes (see Figure 1). We estimate the intensity along with the dynamic reproduction number of the virus within an EM algorithm through a regression on Google mobility indices and demographic covariates in the maximization step. In Section 3, we validate the approach on both short-term and long-term forecasting tasks, showing that the Hawkes process outperforms several benchmark models from “The COVID-19 Forecast Hub Network 2.” The benchmark models include SEIR models from the Columbia University [2], Johns Hopkins University Applied Physics Lab [4], and an ensemble model from Berkeley that uses combined linear and exponential predictors with spatial covariates [20]. We also investigate which covariates and mobility indices are most important for building forecasts of COVID-19 in the U.S. In Section 4, we discuss directions for future research and how the machine learning community may be able to help improve Hawkes process models of COVID-19 as the pandemic continues to unfold. ![Figure 1:](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2020/12/20/2020.06.06.20124149/F1.medium.gif) [Figure 1:](http://medrxiv.org/content/early/2020/12/20/2020.06.06.20124149/F1) Figure 1: The framework of the Hawkes process model for COVID-19 transmission. Demo-graphic features at the county level impact the reproduction number of the Hawkes process. Lagged changes in mobility impact future secondary infections through a convolution with the inter-infection distribution *w*(*t*). The output of the model includes: (1) forecasts of future cases and mortality through simulation of the Hawkes process intensity, (2) an estimate of the dynamic reproduction number of the virus, and (3) regression results that allow for interpretation of the covariates that influence transmission differences across counties. ## 2. Hawkes process model of COVID-19 transmission In this section, we introduce a Hawkes process with spatio-temporal covariates for modeling COVID-19 case and death data. We then discuss the connection of the model to compartment models used in epidemiology and develop an expectation-maximization algorithm for inference. ### 2.1. Incorporating covariates into the Hawkes process We propose a novel Hawkes process model that simultaneously estimates the intensity of events and tracks the dynamic reproduction number of the virus. Given the timestamps (or dates), *𝒯* = {*t*1, *t*2, … *t**n*}, of daily reported positive test cases or deaths, we model the rate of new cases (or deaths) in each country *c* as follows: ![Formula][1] where *µ**c* is the background rate modeling imported infections, *w*(*t*) is the inter-infection time distribution, ![Graphic][2] are mobility indices on day *t*, and ***d****c* = [*d*1, *d*2, *…*]T are static demographic features. The time-varying reproduction number ![Graphic][3] can be interpreted as the average number of secondary infections caused by a primary infection. Because we are modeling reported infections rather than time of exposure, we introduce the parameter Δ to capture a potential lag between a mobility change and the time *t**j* of a reported primary infection. We then model the dynamic reproduction number ![Graphic][4]through a Poissonregression, ![Formula][5] Here we have combined the spatial and temporal covariates to simplify notation in the rest of the paper. Our approach is related to those in recent preprints that incorporate mobility into compartment models [21, 4], however those approaches typically involve large-scale Monte Carlo simulations when performing inference. As we will show, the Hawkes process likelihood can be maximized without simulation via an efficient expectation-maximization algorithm. ### 2.2. Mathematical connection between Hawkes processes and compartmental models Here we briefly review several variations of the Hawkes process in Equation 1 that can be connected to SEIR-type compartment models. The first variant is the SIR-Hawkes process. This model captures the long-term evolution of a pandemic by incorporating a pre-factor that accounts for the dynamic decrease in the number of susceptible individuals [8]: ![Formula][6] Here *I**c*(*t*) is the cumulative number of infections that have occurred up to time *t* and *N* is the total population size. The point process governed by Equation 3 is a continuous time analog of a discrete stochastic SIR model when *w*(*t*) is specified to be exponential [8]. When *w*(*t*) is chosen to be gamma distributed, the Hawkes process also can approximate staged compartment models, like SEIR, if the average waiting time in each compartment is equal [22]. More complex parametric (or non-parametric) inter-infection time distributions *w*(*t*) may be employed within the Hawkes process framework in situations where disease dynamics cannot be captured by a SIR or SEIR model. In the early exponential growth stage of an epidemic, before finite population effects play a role (which is the case with current U.S. data), the Hawkes process in Equation 1 without the prefactor can be used to model new infections arising from SIR and SEIR models, as ![Graphic][7] will be small. While a pre-factor in the Hawkes process involving the cumulative number of infections, i.e., *I**c*(*t*) in Equation 3, it is necessary to model long-term disease dynamics [8]. In the early stages of transmission, a linear Hawkes process can be used (as the prefactor will be close to 1), ![Formula][8] To illustrate this, we simulate a SEIR differential equation as the following: ![Formula][9] The parameters are chosen similar to those of COVID-19: *γ* = .1, *R* = 2, *µ* = 1, and *N* = 5·108. We then fit the linear Hawkes process model in Equation 4 to new infections, *µE*, generated by the SEIR model. We use a non-parametric histogram estimator for *w*(*t*) and find a close fit between the Hawkes process and the SEIR model in Figure 2. ![Figure 2:](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2020/12/20/2020.06.06.20124149/F2.medium.gif) [Figure 2:](http://medrxiv.org/content/early/2020/12/20/2020.06.06.20124149/F2) Figure 2: (Main figure) The red plot shows SEIR differential equation ![Graphic][10], where *β* = *γR*, *γ* = .1, *R* = 2, *µ* = 1, and *N* = 5 · 108. The blue squares show that linear Hawkes process ![Graphic][11] fit to the SEIR curve of new infections. Inset: Non-parametric histogram estimate for *w*(*t*). ### 2.3. EM algorithm for parameter inference We use an expectation–maximization (EM) algorithm to estimate the model in Equation 1. First, we introduce latent random variables, *p**c*(*i, j*), that represent the event that secondary infection *i* is caused by primary infection *j* in county *c*. We let *p**c*(*i, i*) represent the event that case *i* is imported. The complete data log-likelihood is then given by, ![Formula][12] Here we use a Weibull distribution with shape *α* and scale *β* to model inter-infection times, which is used in other studies of epidemics [23, 24, 25] and we find accurately captures transmission in the present data [26]. As the branching structure of the process is unobservable, we optimize the complete data log-likelihood in Equation 6 by iteratively alternating between an expectation step where the branching probabilities *p**c* are estimated and a maximization step where model parameters are updated by maximizing Equation 6. The EM-algorithm is equivalent to a projected gradient ascent on the likelihood of the Hawkes process [27]. #### 2.3.1. Expectation step During the expectation step, we estimate the latent variables *p**c*(*i, j*) for each county. Given the parameters *θ, α, β*, and *µ**c* estimated from the last iteration, the probabilities that case *i* was caused by case *j* (Equation 7a) or was imported (Equation 7b) are given by: ![Formula][13] ![Formula][14] Note that the rate *λ**c*(*t**i*) in Equation 1 is considered to be an aggregation of triggering kernels from all previous historical events (i.e., all *t < t**i*) and the background rate *µ**c*. Therefore, we can consider the probability of case *i* caused by case *j, p**c*(*i, j*), as the contribution of primary infection *j* in the event rate at time *t**i*, i.e., *λ**c*(*t**i*), and *p**c*(*i, i*) can be seen as the contribution of the background rate. #### 2.3.2. Maximization step We then maximize the complete data log-likelihood with respect to the model parameters, conditioned on the estimated branching structure *p**c*(*i, j*). During estimation we do not include event pairs (*i, j*) when *j* is within Ψ = 14 days 3 of the last day of the dataset, as the offspring events *i* have not yet been realized and the inclusion of these incomplete data biases parameter estimates. Given the latent variable *p**c*(*i, j*), the maximization of Equation 6 can be decoupled into three independent optimization problems. Starting with the coefficient *θ* from Poisson regression, the maximization of likelihood function can be rewritten as the following: ![Formula][15] Because the last Ψ days are removed from the dataset and we assume that all possible offspring pairs (i, j) have been observed, we can therefore approximate the integrals for the inter-infection time *w*(*t*) in Equation 6 as is done in [28] by noting that ![Graphic][16]. The optimization problem is therefore a Poisson regression, where we regress the observations ![Graphic][17] against the covariates ![Graphic][18]. ![Formula][19] The same optimization strategy can be applied on the shape and scale parameters, *α* and *β*. The optimization problem can then be solved as a weighted maximum likelihood estimation for the Weibull shape and scale parameters: ![Formula][20] where *p**c*(*i, j*) is the weight of each inter-infection time observation *t**i* *− t**j*. Third, the background rate *µ**c* is determined analytically: ![Formula][21] Pseudo code for the EM algorithm is presented in the Algorithm 1. Algorithm 1 ### EM algorithm optimization ![Figure3](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2020/12/20/2020.06.06.20124149/F3.medium.gif) [Figure3](http://medrxiv.org/content/early/2020/12/20/2020.06.06.20124149/F3) We note that the EM algorithm of the Hawkes process is also connected to the dynamic reproduction number estimator of Wallinga and Teunis [29], as the latter can be viewed as a 1-iteration EM algorithm where a histogram estimator is used for ![Graphic][22] with initial guess ![Graphic][23]. More details are discussed in the following section. ### 2.4. Connection of EM algorithm for Hawkes Process and dynamic R estimator of Wallinga and Teunis Here we make the connection between the EM algorithm for the Hawkes process and the popular dynamic reproduction number estimator of Wallinga and Teunis [30, 29, 23]. The dynamic *R* estimator of Wallinga and Teunis is constructed as follows. The probability that individual *i* at time *t**i* was infected by individual *j* at time *t**j* is defined to be, ![Formula][24] where the distribution of inter-infection times *w*(*t**i* *− t**j*) is typically modeled as Weibull, Gamma, or log-normal [23]. The expected total number of individuals that *j* infects is then given by: ![Formula][25] Wallinga and Teunis then obtain an estimate of the dynamic reproduction number *R*(*t*) by averaging *R**j* over all observed cases *j* where the time of infection *t**j* occurred on day *t*: ![Formula][26] (here *N**t* is the number of observed infections on day *t*). On the other hand, for the Hawkes process the intensity (rate) of infections is modeled as ![Formula][27] where *w*(*t*) and *R*(*t*) are the inter-infection time distribution and dynamic re-production number respectively. Rather than modeling *R*(*t*) as dependent on mobility, we can instead model *R*(*t*) as a piece-wise constant function: ![Formula][28] Here the *I**k* are intervals discretizing time, *B* is the number of such intervals, and *r**k* is the estimated reproduction rate in interval *k*. Given initial guesses for the model parameters and *r**k*, the EM algorithm for the Hawkes process iteratively updates the parameters and branching probabilities by alternating between the **E-step update**: ![Formula][29] ![Formula][30] and **M-step update:** ![Formula][31] ![Formula][32] ![Formula][33] where *T* is the total length of the observation period, *N**k* is the total number of events in interval *k*, and the *w*(*t*) is estimated via weighted MLE (for either a Gamma, Weibull or log-normal) using the inter-event times as observations and branching probabilities as weights. Finally, we can compare Equation 17 to Equation 12. The dynamic *R*(*t*) estimator in Equation 12 is what you obtain with 1 step of the EM algorithm in Equation 17 with initial guess *R*(*t*) *≡* 1, *µ* = 0 and 1 day chosen as the bin width for the histogram estimator. ### 2.5. Hawkes process forecasting Finally, we forecast future events using the branching process representation of the Hawkes process. In particular, for each event in the history of the process we simulate a Poisson random variable with mean ![Graphic][34] representing the number of secondary infections caused by event *j*. For each of these infections we simulate the time of infection by drawing inter-event times from the estimated Weibull distribution. Events falling in the future (past the forecasting date) are then used to update the forecasted intensity through Equation 1. We simulate multiple realizations of this process to estimate a mean intensity forecast along with confidence intervals. ## 3. Experiments and Results In this section we first provide details on the datasets and baseline models used in our experiments. We then discuss the experimental results of several COVID-19 retrospective forecasting tasks at the U.S. county level. The source code and dataset are included in the supplemental material and are available online in a anonymous repository, respectively 4. ### 3.1. Datasets #### 3.1.1. Covid-19 daily cases and deaths reported by The New York Times The New York Times (NYT) [31] 5 releases a daily report of the cumulative numbers of COVID-19 cases in the United States at the county level and over time. While NYT data closely tracks data aggregated by a project at Johns Hopkins University [32], NYT county level reporting started earlier and is therefore used in this study. In total, there are 3, 217 counties with cases and/or deaths in the dataset. The time series data are compiled from state and local government health departments. In order to have sufficient data for statistical inference, we select the counties with confirmed cases greater than and equal to 10 (denoted by *𝒟*conf) and the counties with at least 1 death (denoted by *𝒟*death) by 11/10/2020 when the dataset is curated. In total, there are 2, 824 and 2, 545 counties in these two datasets. Parameter sharing may improve models in counties with less data through variance reduction, but can potentially bias estimates in more populated counties with more cases. We therefore assess model performance over different subsets of counties grouped by case volume. We first rank counties by the number of confirmed cases and deaths by the cut-off date, 11/30/2020, and we then evaluate forecasting accuracy on the top-10% of counties (denoted by ![Graphic][35]), the top-25% counties (denoted by ![Graphic][36]), and counties between the top-25% and top-50% quantiles (denoted by ![Graphic][37]). In Figure 3a and 3b, we present the distribution of the cumulative confirmed cases and deaths at three different quantiles up to the cute-off date 11/30/2020. As the counties at the top-50% have more than 1,000 confirmed cases and 10 deaths, some urban counties, mostly at the top-10%, had already surpassed 10,000 confirmed cases and accumulated more than 300 deaths. In Figure 4a and 4b, we show the daily reported confirmed cases and deaths of top-3 counties in ![Graphic][38] and ![Graphic][39] from *𝒟*conf and *𝒟*death, respectively. Given different demo-graphics and different COVID-19 regulations, each state went through different phases. For example, while Cook, IL seemed to contain the first spike after May, the confirmed cases in Los Angeles, CA seem steadily increase and only slow down after July. The daily death toll of Maricopa, AZ only hit its record high only after August unlike Los Angeles, CA, which had already had their first wave in terms of deaths in April. Overall, the deaths are increasing as the U.S. heads into the winter months. Such differences in infection rates suggest that different public health and social measures may need to be tailored county by county. Therefore, the proposed county-level forecasting model may aid local government policymakers in understanding the demographic and mobility factors that play a role in local reproduction of the virus. ![Figure 3:](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2020/12/20/2020.06.06.20124149/F4.medium.gif) [Figure 3:](http://medrxiv.org/content/early/2020/12/20/2020.06.06.20124149/F4) Figure 3: Distribution of cumulative cases reported at 11/30/2020 at different quantiles ![Figure 4:](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2020/12/20/2020.06.06.20124149/F5.medium.gif) [Figure 4:](http://medrxiv.org/content/early/2020/12/20/2020.06.06.20124149/F5) Figure 4: Example of the daily # of confirmed cases/deaths #### 3.1.2. Google mobility index reports We use Google daily mobility index reports at the county level [33] to estimate a dynamic reproduction number that tracks changes in movement patterns due to stay at home orders (and their staged removal). In total, there are 6 mobility types, including grocery & pharmacy, parks, transit stations, retail & recreation, residential and workplaces. Mobility indices for each category and county are calculated with respect to a baseline value for that day of the week 6. We drop “workplace” mobility from our analysis due to high collinearity with “residential” mobility. Some mobility data are missing when data is sparse for a given date. To deal with missing values, we adopt multivariate feature imputation 7, which estimates each missing mobility entry as a function of other mobility types on the same day in the same county. We show some heatmaps of mobility patterns across counties and time in the Figure 5, where a major change can be observed coinciding with stay at home orders (the first state-wide stay-at-home order was issued at 03/21/2020). Also, the reopening phase in most of the counties can be seen after May. For counties hit by COVID-19 the most (i.e., those in the top-10 %), we can also observe some strict regulations in the “Retail and recreation” areas and better compliance with stay-at-home orders based on high mobility in “Residential” area. ![Figure 5:](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2020/12/20/2020.06.06.20124149/F6.medium.gif) [Figure 5:](http://medrxiv.org/content/early/2020/12/20/2020.06.06.20124149/F6) Figure 5: Heat map of mobility indices across counties in *𝒟*conf and over time. #### 3.1.3. County-level demographic covariates We incorporate spatial demographic features that may be predictive of symptomatic cases of COVID-19 (which are more likely to result in testing and mortality). The dataset is available in a curated form [20] and is derived from CDC and census datasets. The data is at the county level and includes population, median age, number of hospitals and ICU beds, percentage of smokers and diabetes, and heart disease mortality. In Figure 6, we present two examples of spatial demographic features at the county-level used to model variations in the reproduction number. In Figure 6a we observe that both the east and west coasts of the United States are more densely populated compared to midwestern and western regions. Diabetes percentage (shown in Figure 6b), on the other hand, is mostly higher in southern regions of the U.S. ![Figure 6:](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2020/12/20/2020.06.06.20124149/F7.medium.gif) [Figure 6:](http://medrxiv.org/content/early/2020/12/20/2020.06.06.20124149/F7) Figure 6: Examples of spatial demographic and health features at the county-level. ### 3.2. Baseline models and experimental setup for retrospective forecasting comparison We compare the Hawkes process model in Equation 1 with several bench-marks including an SEIR model used in a pandemic tracking dashboard 8 out of Columbia University [2] (denoted by **PROJ**), an geospatial SEIR Model from the Johns Hopkins University Applied Physics Lab [4] (denoted by **BUCKY**), and an ensemble model with linear and exponential predictors from University of California, Berkeley [20] (denoted by **CLEP**). Note that all three bench-marks are tested directly from the released source code and we follow the same experimental protocol as for our proposed model. A simplified Hawkes process, denoted by **Hawkes**, where the reproduction number is held constant is used for comparison to demonstrate the effectiveness of tracking the reproduction number dynamically. We also compare our full Hawkes process model, denoted by ![Graphic][40], to a Hawkes process, **HkPR**m, with only mobility features to determine the marginal improvement of adding demographics. We backtest the six competing models on the 𝒟conf and 𝒟death datasets using the “walk-forward” validation approach. In particular, for 7-day forecasts we first train the models based on cases and deaths before the first cut-off date, 04/15/2020, and then forecast through 04/21/2020. We then slide the forecasting window, training on data before 4/22/2020 and forecasting from 04/22/2020 to 04/28/2020. We repeat this process until the final date of 05/19/2020 (a similar approach is used for 14 and 28 day forecasts). The multivariate imputation models are also trained in the same walk forward fashion to avoid possible data leakage. The hyper-parameter of the lag parameter Δ ranges from 7, 14, 21, and 28 days in our experiments. We evaluate the models according to mean absolute error, MAE, averaged across counties and forecasting windows of the same length, along with percentage error, PE. Mean absolute error (MAE) and the percentage error (PE) are calculated as follows: ![Formula][41] where ![Graphic][42], and ![Graphic][43] are the number of reported events and predicted events, respectively. We also compare the ranking quality of the competing models using Normalized Discounted Cumulative Gain (NDCG) [34], which can be used to evaluate the power of recommendations for counties with potential COVID-19 spikes in the near future. ### 3.3. Experimental results In Table 1 and Table 2, we present the experimental results for 7, 14, and 28 days window forecasts of MAE for all models applied to both confirmed cases (*𝒟*conf) and deaths (*𝒟*death), and in Table 3 and Table 4, we report the results for PE. In terms of MAE and PE, both of our proposed models, **HkPR**m and ![Graphic][44], outperform the benchmarks, **PROJ** and **CLEP**, by a large margin in all three forecasting periods and across quantile subsets of the data.The improvements of MAE and PE can also be seen in the simplistic baseline Hawkes process, **Hawkes**. This suggests that the Hawkes process approach has a good potential on modeling infectious disease due to the self-exciting properties that lie in the COVID-19 cases. View this table: [Table 1:](http://medrxiv.org/content/early/2020/12/20/2020.06.06.20124149/T1) Table 1: MAE on predicted confirmed cases *𝒟*conf View this table: [Table 2:](http://medrxiv.org/content/early/2020/12/20/2020.06.06.20124149/T2) Table 2: MAE on predicted confirmed cases *𝒟*death View this table: [Table 3:](http://medrxiv.org/content/early/2020/12/20/2020.06.06.20124149/T3) Table 3: PE on predicted confirmed cases *𝒟*conf View this table: [Table 4:](http://medrxiv.org/content/early/2020/12/20/2020.06.06.20124149/T4) Table 4: PE on predicted confirmed cases *𝒟*death We found that adding mobility indices improves **Hawkes**, where forecasting accuracy of **HkPR**m also increases across the subsets and all forecasting window. For example, the improvements on MAE over **Hawkes** can go up to 13%, 11%, and 18% for 28 days forecast when **HkPR**m is applied to ![Graphic][45], and ![Graphic][46] in *𝒟*conf, respectively. Similar decrease on MAE can be observed when **HkPR**m is applied to three quantile subsets in death, where **HkPR**m outperforms **Hawkes** by 29%, 25%, and 13% in MAE, respectively. In terms of PE, **HkPR**m stays ahead of **Hawkes** with only one exception at ![Graphic][47] of 𝒟death in 28 days forecasting. This shows that by modeling the reproduction number through daily mobility indices we can enhance the forecasting accuracy and obtain more precise estimation on the spikes in the future. By adding demographic features, we can marginally boost the MAE and PE of ![Graphic][48] over **HkPR**m in some cases. In general, the variation, ![Graphic][49], also shows similar improvements over the benchmarks. In particular, ![Graphic][50] has the best PE enhancement over **HkPR**m at ![Graphic][51] in *𝒟*death for 28 days forecast, which is 20%. This demonstrates that the major forecasting power comes from the joint modeling of mobility indices in the reproduction number while the choices of the background rate and inter-infection distribution may only play a minor part. Moreover, we notice that model **BUCKY** is a competitive baseline in *𝒟*death where it has better accuracy in a few cases, such as MAE and PE ![Graphic][52] and ![Graphic][53] for 28 days forecast. Possible explanation for its advantage could be the CDC-recommended parameters that has been introduced to aid the model training especially for recovery and deaths compartments in its SEIR model. Those parameters include case fatality ratio, case hospitalization ratio, time between death and reporting, etc. However, introducing such pre-trained parameters from CDC may not be practical in real-time forecasting and may potentially bring in the data leakage issue. In Table 5 and Table 6, we present the NDCG results for the ranking evaluation. Generally, the proposed models **HawkPR** have a better NDCG performance when applied to confirmed cases for most of the quantile subsets. In terms of NDCG on the 𝒟death dataset, the baseline Hawkes process, **Hawkes**, performs better in some cases but proposed method consistently comes in second for most of the forecasting window. By generating rankings with good qualities, **HawkPR** can serve as a recommender system for the hotspot counties and the public health policymakers can tailor strategies specifically for each region to contain the virus. We also note that in our model, we are estimating inter-event distributions of observed cases (ignoring asymptomatic cases) and therefore these are observed or “effective” inter-event distributions, rather than true inter-infection distributions based on longitudinal data (we will clarify the language in the revision). We believe this approach is justified by the performance of the model in forecasting observed cases (and this approach is taken in other applications, like seismology where some earthquakes are not observed). View this table: [Table 5:](http://medrxiv.org/content/early/2020/12/20/2020.06.06.20124149/T5) Table 5: NDCG on predicted confirmed cases *𝒟*conf View this table: [Table 6:](http://medrxiv.org/content/early/2020/12/20/2020.06.06.20124149/T6) Table 6: NDCG on predicted confirmed cases *𝒟*death #### 3.3.1. Importance of covariates In Table 7 and Table 8, we show the dynamic reproduction number coefficients of ![Graphic][54] estimated from the Poisson regression component (Equation 2) when applied to 𝒟conf and 𝒟death, respectively. The absolute value of the coefficients indicates the magnitude of the correlation between the reproduction number and the features. With the exception of population estimation in 𝒟death, the coefficients of all variables are statistically significant at the 10*−*7 level or below. The dynamic reproduction number is positively correlated with “Retail and recreation” while negatively correlated with “Residential”, meaning that as mobility shifted away from commercial areas towards residences, the reproduction number decreased. In terms of spatial covariates, the reproduction number is positively correlated with “Population density” and “# of ICU beds.” This suggests that the regions hit the hardest by COVID-19 are mostly urban areas, where most of intensive treatment units are situated. The reproduction number is also negatively correlated with percent of the population with” Diabetes” and “Heart disease mortality rate.” Several possible explanations for this observation include high-risk individuals are being more cautious or that they tend to live in areas with less cases, potentially with less population. View this table: [Table 7:](http://medrxiv.org/content/early/2020/12/20/2020.06.06.20124149/T7) Table 7: Model coefficients (*𝒟*conf) View this table: [Table 8:](http://medrxiv.org/content/early/2020/12/20/2020.06.06.20124149/T8) Table 8: Model coefficients (*𝒟*death) #### 3.3.2. COVID-19 forecasting and reproduce number analysis In Figure 7, we present an example of 28 days projection made through ![Graphic][55] from 10/28/2020 - 11/25/2020 for both 𝒟conf and 𝒟death. We can observe that ![Graphic][56] has very promising results in making projections, especially for the short term future, When the number of forecasting windows increases, the forecasting error increase as the task also being more difficult. Moreover, the narrow confidence interval calculated through 100 Hawkes processes simulations suggests that the the proposed model can make relatively stable forecasting. Lastly, based on the projections, as the number of confirmed cases soon would hit over 500,000 in the top counties including Los Angeles, CA, Cook, IL, and etc. It is imperative to have a robust framework to help governments to design strategies to combat COVID-19 or even more, prioritize vaccine distribution. ![Figure 7:](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2020/12/20/2020.06.06.20124149/F8.medium.gif) [Figure 7:](http://medrxiv.org/content/early/2020/12/20/2020.06.06.20124149/F8) Figure 7: Forecasting for 28 days from 10/28/2020 - 11/25/2020 In Figure 8 and Figure 9, we find that the estimated dynamic reproduction number closely tracks lagged mobility, where the optimal lag parameter is determined as Δ = 14 days for *𝒟*conf and Δ = 21 days for *𝒟*death. The top-2 counties in ![Graphic][57] have estimated reproduction number initially above 2.5. After stay-at-home orders (around 04/11/2020), mobility in residential areas increased. On the other hand, mobility in retail and recreation decreased and the reproduction number fell to around 1, which explains why curves were relatively “flat” in many areas in the U.S. after the lockdown. However, as most of states reopened and lifted up the restrictions, the reproduction number increased after a large population resumed their daily routine, which can be also be observed by the increased mobility in retail and recreation after July. Lastly, to validate the reproduction numbers, we also compare our results to the ones estimated by Stanford University 9 and our estimation match to their findings, which are around 1.5-2.5 initially and 0.5-1.5 up to the beginning of December in 2020. ![Figure 8:](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2020/12/20/2020.06.06.20124149/F9.medium.gif) [Figure 8:](http://medrxiv.org/content/early/2020/12/20/2020.06.06.20124149/F9) Figure 8: Estimated R of confirmed cases *𝒟*conf and lagged mobility changes (Δ = 14 days) ![Figure 9:](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2020/12/20/2020.06.06.20124149/F10.medium.gif) [Figure 9:](http://medrxiv.org/content/early/2020/12/20/2020.06.06.20124149/F10) Figure 9: Estimated R of deaths *𝒟*death and lagged mobility changes (Δ = 21 days) ## 4. Conclusion We showed how Hawkes processes can be combined with spatio-temporal covariates to accurately model COVID-19 transmission and forecast future cases and deaths. The model is competitive with several benchmark models used to forecast the pandemic, achieving improved MAE and NDCG scores on a majority of the experiments we conducted. Our hope is that this work will encourage more research into Hawkes process models of disease spreading that incorporate more advanced features and machine learning principles. As vaccinations are rolled out across the U.S. (given recent FDA approval), local impacts on dynamic reproduction can be flexibly accommodated by our model and used to obtain more accurate and timely forecasts. One potential direction for future research is extending the work here to neural network based point process models [10, 11]. These models may be able to capture more complicated relationships between mobility patterns, demographics, and transmission. The challenges of such an approach include the potential for over-fitting with added parameters and determining how best to realistically model transmission in a neural point process (analagous to the SIR-Hawkes process), which will be important if neural point processes are to be used in long-term forecasting. Many of the preprints and models currently being released on academic archives and websites present a single model without model evaluation, good-ness of fit analysis, or comparison to baselines. Here we believe the “common task framework” [35] could be beneficial in model selection and validation. The machine learning community can contribute to pandemic modeling efforts by performing careful benchmarking of methodologies, creating standardized datasets and tasks, and comparing competing models that come from different fields such as epidemiology, statistics, and machine learning. ## Data Availability The source code and pruned dataset is available online. [https://github.com/chiangwe/HawkPR](https://github.com/chiangwe/HawkPR) ## Footnotes * Email addresses: chiangwe{at}iupui.edu (Wen-Hao Chiang), xl17{at}iupui.edu (Xueying Liu), gmohler{at}iupui.edu (George Mohler) * 1 [https://ai.facebook.com/blog/using-ai-to-help-health-experts-address-the-covid-19-pandemic](https://ai.facebook.com/blog/using-ai-to-help-health-experts-address-the-covid-19-pandemic) * 2 [https://covid19forecasthub.org/community/](https://covid19forecasthub.org/community/) * 3 We choose Ψ = 14 as the incubation period for COVID-19 is thought to extend to 14 days given by the Clinical Care Guidance from the CDC: [https://www.cdc.gov/coronavirus/2019-ncov/hcp/clinical-guidance-management-patients.html](https://www.cdc.gov/coronavirus/2019-ncov/hcp/clinical-guidance-management-patients.html). * 4 [https://anonymous.4open.science/r/d425dcf9-3cfb-4f82-a08c-ee583ab36291/](https://anonymous.4open.science/r/d425dcf9-3cfb-4f82-a08c-ee583ab36291/) * 5 [https://github.com/nytimes/covid-19-data](https://github.com/nytimes/covid-19-data) * 6 The baseline is the median value, for the corresponding day of the week calculated during the 5-week period, 01/03/2020 to 02/06/2020. * 7 [https://scikit-learn.org/stable/modules/impute.html#multivariate-feature-imputation](https://scikit-learn.org/stable/modules/impute.html#multivariate-feature-imputation) * 8 [https://covid19forecasthub.org/community/](https://covid19forecasthub.org/community/) * 9 [https://web.stanford.edu/chadj/Covid/Dashboard.html](https://web.stanford.edu/chadj/Covid/Dashboard.html) * Received June 6, 2020. * Revision received December 19, 2020. * Accepted December 20, 2020. * © 2020, Posted by Cold Spring Harbor Laboratory This pre-print is available under a Creative Commons License (Attribution-NonCommercial-NoDerivs 4.0 International), CC BY-NC-ND 4.0, as described at [http://creativecommons.org/licenses/by-nc-nd/4.0/](http://creativecommons.org/licenses/by-nc-nd/4.0/) ## References 1. [1]. N. M. Ferguson, D. Laydon, G. Nedjati-Gilani, N. Imai, K. Ainslie, M. Baguelin, S. Bhatia, A. Boonyasiri, Z. Cucunubá, G. Cuomo-Dannenburg, et al., Impact of non-pharmaceutical interventions (npis) to reduce covid-19 mortality and healthcare demand. 2020, DOI 10 (2020) 77482. 2. [2]. S. Pei, J. Shaman, Initial simulation of sars-cov2 spread and intervention effects in the continental us, medRxiv. 3. [3]. M. Smith, K. Yourish, S. Almukhtar, K. Collins, D. Ivory, A. McCann, et al., Coronavirus in the us: Latest map and case count, The New York Times. 4. [4]. M. Kinsey, K. Tallaksen, O. R.F., L. Asher, C. Costello, M. Kelbaugh, S. Wilson, Bucky model: a spatial seir model for simulating covid-19 at the county level. URL [https://buckymodel.com/](https://buckymodel.com/) 5. [5]. D. Zou, L. Wang, P. Xu, J. Chen, W. Zhang, Q. Gu, Epidemic model guided machine learning for covid-19 forecasts in the united states, medRxiv. 6. [6]. I. Covid, C. J. Murray, et al., Forecasting covid-19 impact on hospital bed-days, icu-days, ventilator-days and deaths by us state in the next 4 months, medRxiv. 7. [7]. N. P. Jewell, J. A. Lewnard, B. L. Jewell, Caution warranted: using the institute for health metrics and evaluation model for predicting the course of the covid-19 pandemic, Annals of Internal Medicine. 8. [8]. M.-A. Rizoiu, S. Mishra, Q. Kong, M. Carman, L. Xie, Sir-hawkes: linking epidemic models and hawkes processes to model diffusions in finite populations, in: Proceedings of the 2018 World Wide Web Conference, 2018, pp. 419–428. 9. [9]. K. Zhou, H. Zha, L. Song, Learning triggering kernels for multi-dimensional hawkes processes, in: International Conference on Machine Learning, 2013, pp. 1301–1309. 10. [10]. H. Mei, J. M. Eisner, The neural hawkes process: A neurally self- modulating multivariate point process, in: Advances in Neural Information Processing Systems, 2017, pp. 6754–6764. 11. [11]. T. Omi, K. Aihara, et al., Fully neural network based model for general temporal point processes, in: Advances in Neural Information Processing Systems, 2019, pp. 2120–2129. 12. [12]. G. Mohler, M. D. Porter, Rotational grid, pai-maximizing crime forecasts, Statistical Analysis and Data Mining: The ASA Data Science Journal 11 (5) (2018) 227–236. 13. [13]. H. Xu, H. Zha, A dirichlet mixture model of hawkes processes for event sequence clustering, in: Advances in Neural Information Processing Systems, 2017, pp. 1354–1363. 14. [14]. N. Du, M. Farajtabar, A. Ahmed, A. J. Smola, L. Song, Dirichlet-hawkes processes with applications to clustering continuous-time document streams, in: Proceedings of the 21th ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, 2015, pp. 219–228. 15. [15]. H. Xu, M. Farajtabar, H. Zha, Learning granger causality for hawkes processes, in: International Conference on Machine Learning, 2016, pp. 1717– 1726. 16. [16]. U. Upadhyay, A. De, M. G. Rodriguez, Deep reinforcement learning of marked temporal point processes, in: Advances in Neural Information Processing Systems, 2018, pp. 3168–3178. 17. [17]. S. Li, S. Xiao, S. Zhu, N. Du, Y. Xie, L. Song, Learning temporal point processes via reinforcement learning, in: Advances in neural information processing systems, 2018, pp. 10781–10791. 18. [18]. G. Mohler, J. Carter, R. Raje, Improving social harm indices with a modulated hawkes process, International Journal of Forecasting 34 (3) (2018) 431–439. 19. [19]. A. Reinhart, J. Greenhouse, Self-exciting point processes with spatial covariates: modelling the dynamics of crime, Journal of the Royal Statistical Society: Series C (Applied Statistics) 67 (5) (2018) 1305–1329. 20. [20]. N. Altieri, R. L. Barter, J. Duncan, R. Dwivedi, K. Kumbier, X. Li, R. Netzorg, B. Park, C. Singh, Y. S. Tan, et al., Curating a covid-19 data repository and forecasting county-level death counts in the united states, arXiv preprint arxiv:2005.07882. 21. [21]. A. C. Miller, N. J. Foti, J. A. Lewnard, N. P. Jewell, C. Guestrin, E. B. Fox, Mobility trends provide a leading indicator of changes in sars-cov-2 transmission, medRxiv. 22. [22]. A. L. Lloyd, Destabilization of epidemic models with the inclusion of realistic distributions of infectious periods, Proceedings of the Royal Society of London. Series B: Biological Sciences 268 (1470) (2001) 985–993. 23. [23]. T. Obadia, R. Haneef, P.-Y. Boëlle, The r0 package: a toolbox to estimate reproduction numbers for epidemic outbreaks, BMC medical informatics and decision making 12 (1) (2012) 147. 24. [24]. B. J. Cowling, M. S. Lau, L.-M. Ho, S.-K. Chuang, T. Tsang, S.-H. Liu, P.-Y. Leung, S.-V. Lo, E. H. Lau, The effective reproduction number of pandemic influenza: prospective estimation, Epidemiology (Cambridge, Mass.) 21 (6) (2010) 842. [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1097/EDE.0b013e3181f20977&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=20805752&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2020%2F12%2F20%2F2020.06.06.20124149.atom) [Web of Science](http://medrxiv.org/lookup/external-ref?access_num=000282600600015&link_type=ISI) 25. [25]. J. Hellewell, S. Abbott, A. Gimma, N. I. Bosse, C. I. Jarvis, T. W. Russell, J. D. Munday, A. J. Kucharski, W. J. Edmunds, F. Sun, et al., Feasibility of controlling COVID-19 outbreaks by isolation of cases and contacts, The Lancet Global Health. 26. [26]. A. L. Bertozzi, E. Franco, G. Mohler, M. B. Short, D. Sledge, The challenges of modeling and forecasting the spread of covid-19, Proceedings of the National Academy of Sciences 117 (29) (2020) 16732–16738. arXiv:[https://www.pnas.org/content/117/29/16732.full.pdf](https://www.pnas.org/content/117/29/16732.full.pdf), xdoi: 10.1073/pnas.2006520117. URL [https://www.pnas.org/content/117/29/16732](https://www.pnas.org/content/117/29/16732) [Abstract/FREE Full Text](http://medrxiv.org/lookup/ijlink/YTozOntzOjQ6InBhdGgiO3M6MTQ6Ii9sb29rdXAvaWpsaW5rIjtzOjU6InF1ZXJ5IjthOjQ6e3M6ODoibGlua1R5cGUiO3M6NDoiQUJTVCI7czoxMToiam91cm5hbENvZGUiO3M6NDoicG5hcyI7czo1OiJyZXNpZCI7czoxMjoiMTE3LzI5LzE2NzMyIjtzOjQ6ImF0b20iO3M6NTA6Ii9tZWRyeGl2L2Vhcmx5LzIwMjAvMTIvMjAvMjAyMC4wNi4wNi4yMDEyNDE0OS5hdG9tIjt9czo4OiJmcmFnbWVudCI7czowOiIiO30=) 27. [27]. E. Lewis, G. Mohler, A nonparametric em algorithm for multiscale hawkes processes, Journal of Nonparametric Statistics 1 (1) (2011) 1–20. 28. [28]. F. P. Schoenberg, Facilitated estimation of etas, Bulletin of the Seismological Society of America 103 (1) (2013) 601–605. 29. [29]. J. Wallinga, P. Teunis, Different epidemic curves for severe acute respiratory syndrome reveal similar impacts of control measures, American Journal of epidemiology 160 (6) (2004) 509–516. [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1093/aje/kwh255&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=15353409&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2020%2F12%2F20%2F2020.06.06.20124149.atom) [Web of Science](http://medrxiv.org/lookup/external-ref?access_num=000223938000001&link_type=ISI) 30. [30]. S. Cauchemez, P.- Y. Boëlle, C. A. Donnelly, N. M. Ferguson, G. Thomas, G. M. Leung, A. J. Hedley, R. M. Anderson, A.-J. Valleron, Real-time estimates in early detection of sars, Emerging infectious diseases 12 (1) (2006) 110. 31. [31].The new york times, coronavirus in the u.s.: Latest map and case count (Mar 2020). URL [https://www.nytimes.com/interactive/2020/us/coronavirus-us-cases.html](https://www.nytimes.com/interactive/2020/us/coronavirus-us-cases.html) 32. [32]. E. Dong, H. Du, L. Gardner, An interactive web-based dashboard to track covid-19 in real time, The Lancet Infectious Diseases. 33. [33].Google, Covid-19 community mobility report (2020). URL [https://www.google.com/covid19/mobility/](https://www.google.com/covid19/mobility/) 34. [34]. K. Järvelin, J. Kekäläinen, Cumulated gain-based evaluation of ir techniques, ACM Transactions on Information Systems (TOIS) 20 (4) (2002) 422–446. 35. [35]. D. Donoho, 50 years of data science, Journal of Computational and Graphical Statistics 26 (4) (2017) 745–766. [1]: /embed/graphic-2.gif [2]: /embed/inline-graphic-1.gif [3]: /embed/inline-graphic-2.gif [4]: /embed/inline-graphic-3.gif [5]: /embed/graphic-3.gif [6]: /embed/graphic-4.gif [7]: /embed/inline-graphic-4.gif [8]: /embed/graphic-5.gif [9]: /embed/graphic-6.gif [10]: F2/embed/inline-graphic-5.gif [11]: F2/embed/inline-graphic-6.gif [12]: /embed/graphic-8.gif [13]: /embed/graphic-9.gif [14]: /embed/graphic-10.gif [15]: /embed/graphic-11.gif [16]: /embed/inline-graphic-7.gif [17]: /embed/inline-graphic-8.gif [18]: /embed/inline-graphic-9.gif [19]: /embed/graphic-12.gif [20]: /embed/graphic-13.gif [21]: /embed/graphic-14.gif [22]: /embed/inline-graphic-10.gif [23]: /embed/inline-graphic-11.gif [24]: /embed/graphic-16.gif [25]: /embed/graphic-17.gif [26]: /embed/graphic-18.gif [27]: /embed/graphic-19.gif [28]: /embed/graphic-20.gif [29]: /embed/graphic-21.gif [30]: /embed/graphic-22.gif [31]: /embed/graphic-23.gif [32]: /embed/graphic-24.gif [33]: /embed/graphic-25.gif [34]: /embed/inline-graphic-12.gif [35]: /embed/inline-graphic-13.gif [36]: /embed/inline-graphic-14.gif [37]: /embed/inline-graphic-15.gif [38]: /embed/inline-graphic-16.gif [39]: /embed/inline-graphic-17.gif [40]: /embed/inline-graphic-18.gif [41]: /embed/graphic-30.gif [42]: /embed/inline-graphic-19.gif [43]: /embed/inline-graphic-20.gif [44]: /embed/inline-graphic-21.gif [45]: /embed/inline-graphic-22.gif [46]: /embed/inline-graphic-23.gif [47]: /embed/inline-graphic-24.gif [48]: /embed/inline-graphic-25.gif [49]: /embed/inline-graphic-26.gif [50]: /embed/inline-graphic-27.gif [51]: /embed/inline-graphic-28.gif [52]: /embed/inline-graphic-29.gif [53]: /embed/inline-graphic-30.gif [54]: /embed/inline-graphic-31.gif [55]: /embed/inline-graphic-32.gif [56]: /embed/inline-graphic-33.gif [57]: /embed/inline-graphic-34.gif