The global spread of COVID-19, the disease caused by the novel coronavirus SARS-CoV-2, has affected nearly everyone's lives on the planet. Even the largest economies' resources have been strained due to the large infectivity and transmissibility of COVID-19. As the number of cases of COVID-19 continues increasing, understanding finer-grained spatio-temporal dynamics of this disease as well as some of the leading factors affecting disease transmissions, such as community mobility and population, is critical to helping officials make policy decisions and curb the further spread of the disease.
Most of the previous research aimed at studying the spread of COVID-19 has focused on two key measurements: the number of confirmed cases and the number of deaths. Cases going up or down over time shed light on the rate of spread of COVID-19 at a given point in time -but it is only valid if enough people get tested. The limited testing ability resulted in a serious underestimation of COVID-19 cases in the pandemic's early stages [20] . For example, when there was not enough testing capacity, as was the case in New York City in March 2020, the number of cases reported was clearly an undercount of true cases, estimated to be much larger (up by a factor of 10) [14, 24] . Some studies have circumvented the problem of underestimation by considering the case positivity rate, which measures the percentage of total COVID-19 tests conducted that are positive. However, most of the widely-used COVID-19 data sets, such as the COVID Tracking Project [29] , only collect the total number of people with a completed polymerise chain reaction (PCR) test that returns positive as reported by the state or territory, which has a much lower spatial resolution (state-level) in comparison with the cases and deaths data (county-level). Such coarse-grained testing numbers would introduce extra noise to our model and would most likely be incapable of improving the confirmed cases prediction accuracy at the county level. Deaths are also an important metric that most people care about when it comes to the virus's ultimate epidemiological impact. In contrast to the number of confirmed cases, the number of deaths is a good and accurate indicator for evaluating how serious a burden this pandemic is causing, not only on health care systems but also on the general public's mental health and well-being. Some epidemiological studies, such as [19] , also recommend tracking deaths, even though deaths lag behind new cases, typically by two weeks to a month.
A large amount of fine-grained data offers a unique opportunity to study the disease's spread dynamics from a micro-level view. For the United States, a number of teams have been working on collecting comprehensive COVID-19 tracking data, including counts of cases and deaths at the county-level on a daily basis. Such kind of data gives us a general picture of how the virus is spreading across metropolitan and micropolitan counties and how such dynamics are evolving over time. Aside from considering the cases and the deaths, we also aim to study other important local factors in the transmission of the COVID-19. Recent studies [23] on the spread of COVID-19 show that besides the distance to the epicenter, other factors, such as subway and airport, are positively connected with the virus transmission. Moreover, both urban areas and population density are positively associated with the spread of COVID-19 after the outbreak. The proportion of the elderly population has also been identified as a key factor in the death rate. Therefore, we consider the within-community mobility and two critical demographic factors by taking advantage of the COVID-19 Community Mobility Reports [12] and the American Community Survey (ACS) [3] . These two data sets are publicly available and include detailed county-level statistics that provide insights into what has changed in response to policies aimed at combating the COVID-19 and what factors may affect the disease's transmission. As illustrated in Fig. 1 , in our model, we assume the numbers of cases and deaths in each county to depend on the neighboring counties and major metropolitan areas in the U.S., which we refer as hubs in spreading the disease. Local community mobility and demographic factors, including population and elderly population, are considered local covariates in the model, which also play a crucial role in the final number of deaths.
In this paper, we use a data-driven method incorporating a large-scale data set from multiple sources to predict the deaths and the confirmed cases of COVID-19 at the county level in the United States. To improve our model's reliability and performance, we emphasize recovering county-level deaths' trajectories Figure 1 : An example of spatio-temporal covariates in our model for Coconino County, Arizona. Based on the counties in the U.S. as fundamental units, we assume the number of confirmed cases and deaths of COVID-19 reported in a given county are jointly related to the numbers reported in its adjacent counties (they are Kane, San Juan, Navajo, Gila, Yavapai, Mohave for Coconino in this example) and ten selected nationwide hubs (including San Francisco, Los Angeles, Seattle, Chicago, Atlanta, Miami, Washington, D.C., Boston, and New York). These two numbers also depend on some local covariates, such as community mobility level and some county's demographic factors. and focus less on improving confirmed cases' predictive accuracy. Our method's most notable contribution is considering the spatial structure among hubs and neighboring counties in modeling the cross-correlation between cases and deaths. We also present the effect of a wide variety of geographic community mobility and social demographic factors on the spread of COVID-19. Our approach drastically differs from previous studies [7, 2, 1] , in which the number of cases and deaths, and other covariates, including the community mobility and social demographic factors, are inter-linked through a vector autoregressive process. Our model shows that these hubs play a pivotal role in spreading the disease. We also find that both cases and deaths are significantly related to the total population at the local level and that deaths are also positively associated with the proportion of the elderly population. Additionally, we found that confirmed cases are not significantly related to the proportion of the elderly population, which may prove that the disease was mostly circulating among young people in its later stage. In particular, while we identify a spike in cases since the beginning of the summer, we do not observe a clear spike in deaths. This may be explained by the fact that a larger proportion of young people, who are generally at lower risk of death, were infected in the more recent pandemic stages. The remainder of the article is organized as follows. We first review related works in the rest of this section, followed by describing the data sets we have used in Section 2. Section 3 presents our proposed vector autoregressive model with spatial structure incorporated. We demonstrate the effectiveness of our model and discuss its interpretation in Section 4. Lastly, the article concludes with discussions and future research directions in Section 5.
Related work Compartmental models have been widely used in infectious disease epidemiological studies. In the SIR model [13] , one of the simplest compartmental models, the population is assigned to three components: S (susceptible), I (infectious), and R (recovered). These variables (S, I, and R) represent the number of people in each compartment. The transition between different compartments is modeled using a set of coupled differential equations. Based on the SIR model, many variants have been proposed in the last decades, including the SIRD model [10, 4] that considers deceased individuals, and the SEIR model [15, 33, 16 ] that considers exposed period during which individuals have been infected but are not yet infectious themselves, to name a few. The total population is usually assumed to be fixed in the compartmental models; therefore, it works well when modeling nationwide data. However, in our high-resolution modeling, each county's population is of high variability due to dynamics across the county. Therefore, we use a spatio-temporal model instead to capture the influence of major big cities and neighboring counties without fixing each county's population.
Studies evaluating the spatial spread of the COVID-19 pandemic are scarce [27] . However, understanding the spatial spread of the COVID-19 outbreak is critical to predicting local outbreaks and developing public health policies during the early stages of COVID-19. Previous studies have described the spatial spread of severe acute respiratory syndrome (SARS) in Beijing and mainland China [25, 9, 18, 17, 8] using limited or localized data. One study also considered the different types of connections between a few cities to calculate the spatial association [25] . There is also prior work using the multivariate Hawkes process to model the conditional intensity of new COVID-19 cases and deaths in the U.S. at the county level [5] , without considering the influence of the big cities and other important demographic factors.
There are also various efforts studying impacts from other aspects, such as temperature, humidity [30, 28] , age, gender [8] , and travel restrictions [6, 21] . Most of these studies are constrained on a relatively small scale because of limited data at the early stage of the pandemic.
We have used three comprehensive datasets in this study, including confirmed cases and deaths of COVID-19, community mobility data, and demographic census data. These datasets play an important role in helping us understand the spatio-temporal dynamics of COVID-19 transmission.
Confirmed cases and deaths of COVID- 19 We used the dataset from The New York Times [32] , which are based on reports from state and local health agencies. The data is the product of dozens of journalists working across several time zones to monitor news conferences, analyze data releases, and seek public officials' clarification on how they categorize cases. The data includes two parts: (i) confirmed cases are counts of individuals whose coronavirus infections were confirmed by a laboratory test and reported by a federal, state, territorial, or local government agency. Only tests that detect viral RNA in a sample are considered confirmatory. These are often called molecular or reverse transcription-polymerase chain reaction (RT-PCR) tests; (ii) confirmed deaths are individuals who have died and meet the definition for a confirmed COVID-19 case. Some states reconcile these records with death certificates to remove deaths from their count, where COVID-19 is not listed as the cause of death. This data have removed non-COVID-19 deaths among confirmed cases according to the information released by health departments, i.e., in homicide, suicide, car crashes, or drug overdose. All cases and deaths are counted on the date they are first announced. In practice, we have observed periodic weekly oscillations in daily reported cases and deaths, which could have been caused by testing bias (higher testing rates on certain days of the week). To reduce such bias, we aggregate the number of cases and deaths of each county by week.
Community mobility As global communities respond to COVID-19, we've heard from public health officials that the same type of aggregated, anonymized insights we use in products such as Google Maps could be helpful as they make critical decisions to combat COVID-19. The COVID-19 Community Mobility Reports [12] aim to provide insights into what has changed in response to policies aimed at combating COVID-19.
The reports record people's movement by county daily, across different categories such as retail and recreation, groceries and pharmacies, parks, transit stations, workplaces, and residential. The data shows how visitors to (or time spent in) categorized places change compared to the baseline days (in percentage). The negative percentage represents the level of mobility is lower than the baseline, and the positive percentage represents the opposite. A baseline day represents a normal value for that day of the week. The baseline day is the median value from the five-week period from January 3rd to February 6th, 2020. To match the temporal Demographic census Data from the American Community Survey (ACS) [3] , provided by the U.S. Census Bureau, is the comprehensive source for comprehensive information about the population, demographic, and economic status of each zip code region in the U.S. Unlike the census data, which takes place every ten years, the ACS is conducted every year. The latest ACS data are available in the year 2018. Some demographic factors are particularly useful for us in understanding how population distribution affects the spread of disease (by correlating the local socio-economic profile with its confirmed cases and deaths), and these factors contain essential information about the development and economic growth of different areas. To match the spatial resolution with the COVID-19 data, we aggregate the zipcode regions' demographic data in the same county. We selected two leading factors that affect the spread and the infection of the disease, i.e., total population, and the proportion of the elderly with an age of 65 or older [26] .
This section presents our statistical model that captures the spatio-temporal dynamics of the spread of COVID-19. We begin with a brief description of the problem setup and notations, then jointly model confirmed cases and deaths as a vector autoregressive process in Section 3.2. The essential notations defined in this section are also summarized in Table 1 .
Consider confirmed cases and deaths of the COVID-19 in N counties and T weeks (recall that we aggregated these numbers by week to reduce bias). Let I = {i = 1, . . . , N } be the set of counties and T = {t = 1, . . . , T } be the set of weeks starting from February 9th until August 9th. We assume there is a set of counties I = {i = 1, . . . , N } ⊂ I playing a significant role in spreading the disease due to their high population density and well-developed transportation network connecting to other major cities in the U.S.. We refer to these counties as hubs, and the selected hubs are marked in Fig. 1 . Denote the number of confirmed cases and deaths in county i ∈ I and week t ∈ T as c i,t ∈ Z + and d i,t ∈ Z + , respectively. In our setting, T = 27, N = 3144, and N = 10.
We also consider K mobility categories and L demographic factors as covariates of the model, where K = 6 and L = 2. Let K = {k = 1, . . . , K} be the set of community mobility categories and L = {l = 1, . . . , L} be the set of demographic factors. Denote the mobility score in category k ∈ K for county i ∈ I in week t ∈ T as m i,k,t ∈ R, and denote the data of demographic factor l ∈ L for county i ∈ I as z i,l ∈ R + . Let c t := [c 1,t , . . . , c N,t ] and d t := [d 1,t , . . . , d N,t ] denote the confirmed cases and deaths in week t ∈ T , Set of demographic factors. c i,t ∈ Z + Number of confirmed cases for county i ∈ I in week t ∈ T . d i,t ∈ Z + Number of deaths for county i ∈ I in week t ∈ T . z i,l ∈ R + Data of demographic factor l ∈ L for county i ∈ I. m i,k,t ∈ R Data of mobility category k ∈ K for county i ∈ I in week t ∈ T . 3.2 A = {(i, j) : i, j ∈ I} Set of all county pairs in the U.S. that i, j are adjacent to each other or one of i, j is a hub.
Coefficients of past confirmed cases between county i, j ∈ I for τ weeks ago.
Coefficients of past deaths between county i, j ∈ I for τ weeks ago.
Coefficients of past confirmed cases between county i, j ∈ I for τ weeks ago. µ k,τ ∈ R Coefficient for mobility category k ∈ K in the past τ -th week w.r.t. the number of cases. ν k,τ ∈ R Coefficient for mobility category k ∈ K in the past τ -th week w.r.t. the number of deaths. υ l ∈ R Coefficient for demographic factor l ∈ L w.r.t. the number of cases. ζ l ∈ R Coefficient for demographic factor l ∈ L w.r.t. the number of deaths.
respectively. Let m k,t := [m 1,k,t , . . . , m N,k,t ] denote the score of community mobility category k ∈ K for all counties i ∈ I in week t ∈ T . Let z l := [z 1,l , . . . , z N,l ] denote the data of demographic factor l ∈ L for all counties i ∈ I.
We consider a linear spatio-temporal autoregressive model where the number of confirmed cases (c t ) and deaths (d t ) is a time series regressed on their previous values and the mobility covariate m k,t and demographic covariate z l . Denote the time window's length that we consider in the past (the memory depth) as p. Based on previous studies [31] , it is known that the COVID-19 virus has an incubation period of around two weeks. Therefore, we choose p = 2 throughout this paper.
Define the augmented observation vector as (which contains both confirmed case and death counts):
Then our spatio-temporal model can be written as a vector autoregressive (VAR) process:
where ⊗ is the Kronecker product and
In our model (1), the first-term captures the dependence on past confirmed cases and deaths; the second term captures the influence of past local community mobility; the third term captures the influence of local demography, which is held constant over time. Specifically, B τ and A τ contain the autoregressive coefficients for the number of confirmed cases and deaths, respectively; H τ describes the dependence of the current number of deaths on the number of confirmed cases in τ weeks ago. These three matrices share the same sparse spatial structure, where entry at (i, j) is zero if county i and county j is not adjacent, and none of them is the hub. Formally, the set of adjacency pairs is defined by A = {(i, j) ∈ I : (i, j) is an edge of the graph G}; each node of G denotes a county, and there is an edge between two nodes whenever the corresponding counties are geographically adjacent or one of them is a hub. The µ k,τ , ν k,τ , υ l , and ζ l are four scalar coefficients.
To be specific, µ k,τ , ν k,τ represent the coefficients for the local community mobility score in category k in τ weeks ago with respect to the corresponding number of confirmed cases and deaths, respectively. Similarly, υ l , ζ l represent the coefficients for local demographic factor l with respect to the corresponding number of confirmed cases and deaths, respectively. The spatial covariance matrix between the noise at counties i and j is denoted as the (i, j)-th entry of Σ η ; it is a function of their Euclidean distance s ij and is parameterized by η. Some commonly used spatial models include: Gaussian model [22] , Exponential model [11] , and Matérn model [11] . Here we adopt the exponential spatial covariance model Σ η (i, j) = η exp{−ηs ij }, where η is a pre-specified parameter, which controls the rate of spatial decay. In this paper, we specify a reasonable value of the parameter η = 10 3 . We aim to fit the model (1) for confirmed cases and deaths jointly by minimizing the prediction error. Define the set of parameters θ = {Λ, ω, γ} ∈ Θ, where Θ is the set containing all feasible values. For a pre-specified hyper-parameter δ ∈ [0, 1], the loss function is defined as a weighted combination of quadratic loss functions for death and confirmed case residuals:
where ε t,c denotes the confirmed case prediction residual
and ε t,d denotes the death prediction residual
The hyper-parameter δ controls the proportion of death prediction loss. In practical terms, we emphasize the importance of death, and hence we choose δ = 0.9 empirically. The reason is that it is known that the confirmed cases are quite noisy and can depend on the capacity of testings. The parameters θ can be estimated by solving the following optimization with a regularization function:
where λ 1 ≥ 0 is a parameter that controls the importance of the regularization term, and R(θ) is the elastic net type regularization function (with hyper-parameter λ 2 ∈ [0, 1]) given by
where 1 A {x} is the indicator function, i.e., taking the value 1 if x ∈ A otherwise 0; λ 2 is the 1 penalty ratio in the regularization function; α i,j,τ , β i,j,τ , h i,j,τ are the entries of matrices A τ , B τ , H τ , respectively.
Our model's salient feature is that we consider the underlying spatio-temporal structure between the number of confirmed cases and deaths. If there is no specific structure in coefficient matrices, our methods looks on the surface to be a naive linear model but requires to solve a large-scale high-dimensional optimization problem, which contains 79,077,916 parameters (variables in the optimization problem) with only 84,888 data points. Instead of solving such difficult problems directly, we tackle this challenge by exploiting the sparse spatial structure and only consider the correlation between adjacent counties and hubs, which leads to a significant reduction in the number of parameters (less than 80,000). In addition, the lower triangular structure of Λ τ matrix (including B τ , H τ , and A τ ) captures the causal relationship we believe exists in the confirm case to the death count, but not the other way around. To be exact, we assume the number of confirmed cases in the past will result in the change of both the confirmed cases and deaths in the future, while the number of deaths only relates to the future's deaths. The regularization term we devised in Section 3.2 also plays a big part in achieving the ideal results. This elastic net-based method linearly combines the lasso and ridge regression penalties on B τ , H τ , and A τ to encourage sparse spatial correlation and stabilize the solution at the same time. The hyper-parameters λ 1 and λ 2 in the regularization term are chosen by 5-fold cross-validation, where the optimal choices are λ 1 ≈ 10 2 and λ 2 ≈ 10 −1 for the fitted model.
Here we solve the optimization problem by gradient descent. To fit the model, we first standardize the data of covariates and feed all the data as a single batch in one iteration, then descend the gradients of the parameters with respect to the loss defined in (2) until the model converges. To perform a one-week ahead prediction, we feed all the data before that week as a single batch in one iteration and follow the same gradient descent procedure described above. The model normally takes about 500 iterations to reach the convergence with (θ) ≈ 1.41 × 10 3 .
Now we report the results of our study. First, we describe the validation of the proposed modeling method by evaluating the in-sample estimation error. We also compare our approach against a trivial reference model, which is the persistence model. Then we demonstrate the explanatory components of our model by showing the spatio-temporal dynamics between the number of COVID-19 cases and other covariates discovered from our fitted model. 
To evaluate our method, we first compared the county-level in-sample estimation on deaths and confirmed cases using our model. The county-level in-sample estimation includes all the confirmed cases and deaths for 3,144 counties and the past 27 weeks. We fit the model using the entire data set from February 9th, 2020 to August 9th, 2020, and estimate the parameters according to (1) for each county at each week afterward. To be specific, the upper half of x t on the left-hand side of (1) is the estimated number of confirmed cases at week t, and the lower half is the estimated number of deaths at week t. For deaths count, we only present the results from April 26th, 2020, because data before May is accurate because of the low testing rate. Here, for the ease of presentation, we only focus on the mainland and do not consider Hawaii, Alaska, and other unincorporated territories of the United States in this paper.
We report the results of in-sample estimation by aggregating the county-level estimated numbers in the same state. As shown in Fig. 4 , we select eight major states with the highest total number of deaths in the U.S.. The shaded area indicates the true number of deaths reported in the COVID-19 data set, and the red line indicates the in-sample estimated deaths by our model. In addition to the in-sample estimation, we also show the one-week ahead prediction in death represented by blue dash lines. The prediction at one week is obtained by fitting an individual model using the data in the past of that week and estimating the fitted model next week. We observe that the red lines perfectly match the true death trajectories and the blue dash line present promising predictive performance in the number of deaths. We also report the state-level in-sample estimated confirmed cases in Fig. 5 . Similarly, the shaded area indicates the true number of confirmed cases reported in the COVID-19 data set, and the green line indicates the in-sample estimated confirmed cases by our model. The results further confirm that our model can also "recover" the confirmed cases' trajectories with relatively high accuracy.
More quantitative results are summarized in Fig. 6 . The heatmaps show the mean absolute error (MAE) of the county-level estimation within a particular state and at certain weeks. The average MAEs over states We can also observe that: for the confirmed cases, our model tends to achieve better performance for the states with larger populations, such as Florida, New York, Texas, etc.; for the deaths, our model has a balanced performance in each state, and the MAE is getting better (smaller) and becoming more stable after the summer surge of the COVID-19 (from June to July).
Our study focuses on exploring the in-sample explanatory content of predetermined factors in our model. We fit the model using the entire data set collected from three data mentioned above sources in Section 2, and interpret the model by examining its fitted coefficients.
The experimental results demonstrate a distinctive underlying spatio-temporal pattern between confirmed cases and deaths of the COVID-19. In Fig. 7 , we report the coefficients of five representative hubs in Λ 1 . To be specific, hubs' coefficients in B 1 , H 1 , A 1 indicate their spatial dependencies between each pair of past cases and current cases, past cases and current deaths, and past deaths and current deaths, respectively. As we can see, hubs have very strong "radiating" power on most of the United States regions and contribute a great deal to promote or curb the spread of the COVID-19. However, the rural area with lower population density in the Central United States is not significantly influenced by the hubs. Also, the hubs situated in the Northern United States (e.g., Chicago, New York) are negatively related to the spreading diseases to the other regions (in blue), which appear to have better controls on the expansion of the virus. In contrast, the hubs in the Southern United States (e.g., Dallas, Houston, Miami) usually are positively related to the increases of both cases and deaths in other regions (in red). In addition, the result presents some other interesting findings: some hubs show two opposite influences on the cases and deaths in the same region. For example, we see that on the one Figure 7 : Examples of coefficients for hubs in matrix Λ 1 . For instance, the color depth of any county i in (a) represents the value of coefficient α i,j,1 in A 1 , where county j is Chicago. Counties in blue indicate their current number of deaths is positively related to its number of deaths in the last week; counties in red are the opposite; counties in white represent no discernable correlation between the two numbers. Coefficients of different hubs show the various spatial pattern in "spreading" or "controlling" the disease.
hand, Fig. 7 (m) shows that the number of deaths in Miami is negatively related to the deaths in the New England area of the United States. On the other hand, Fig. 7 (n) shows that the number of cases in Miami is positively related to the cases in the same area. Some hubs contribute to the increase of cases or deaths in one region, reducing the cases or deaths in other regions. For example, Fig. 7 (g) and (j) show that Dallas and Houston have a positive impact on the New England area in the United States and have a negative impact on Florida and California. In Fig. 8 , we also present three typical pairs of comparisons for coefficients between one-week lag and two-week lag: coefficients of Atlanta in A 1 and A 2 , coefficients of Seattle in B 1 and B 2 , and coefficients of Los Angeles in H 1 and H 2 . All three comparisons share one thing in common: coefficients of different time lag have a similar spatial pattern, but the overall coefficients of two-week lag are relatively smaller than corresponding ones of one-week lag. This indicates that the last week has a stronger influence.
Dependence on local covariates Table 2 summarizes the fitted coefficients of local covariates in the model. As we can see, most of the covariates are statistically significant, with small p-values (< .05) except for the proportion of the elderly population age 65 and older. The positive coefficients are in bold, which indicates a positive correlation between the covariates and the cases or deaths. In particular, we observe that, for the cases, the coefficients of mobility in workplaces, retail, and recreation, transit stations have large positive values (> 9 × 10 2 ), which indicates that the increase of mobility in these areas led to the rapid spread of the COVID-19. However, things are the opposite for the deaths, where the coefficients of mobility in grocery and pharmacies, parks, residential have large positive values (> 2 × 10 2 ). Moreover, the population's coefficient for the cases is significantly larger than the other covariates, and it confirms that the population density is the dominant factor in spreading the disease. Last, we have found the proportion of the elderly population is significantly related to the deaths and has no clear connection to the cases.
While still in the development stages, the proposed spatio-temporal model has shown immense promise in modeling and predicting the deaths and confirmed cases of COVID-19 in the United States. Nevertheless, there remains numerous open questions and rooms for improvements. For example, the uncertainty in the count data commonly exists and can affect accuracy. It would be interesting to incorporate the serology data as an additional data source to calibrate our model. To avoid the issue of output negative output, we may adapt the current problem into a Poisson regression with log-linear model. It assumes the response variable x t has a Poisson distribution, and assumes the logarithm of its expected value can be modeled by the linear model defined in (1) . In particular, this adaption plays a vital role in predicting states with fewer confirmed cases and deaths, such as Hawaii and Delaware.