The objective of this study is to reveal and evaluate early indicators for COVID-19 spread risk in cities during the initial stages of the outbreak using measures of human activities derived from digital trace data. An arguably unprecedented global pandemic, the coronavirus disease 2019 (COVID- 19) has infected millions of people worldwide with a mortality rate of 6.6% and a high infection rate [1] [2] . Since the spread of COVID-19 is highly dependent on human activities, incidence of infection could be contained by restricting human activities and mobility [3] . Many countries and authorities have implemented various non-pharmaceutical interventions (e.g., shelter-in-place orders, regional lockdowns, and travel restrictions), which were undertaken to slow the spread of disease by disrupting transmission chains by restricting human mobility and activities. Such social distancing and activity reduction interventions have proven to be critical in slowing down the spread of pandemics both in previous epidemics [4] and during COVID-19 [5] [6] [7] .
While reduction in human activities is considered an effective measure for containing epidemics and pandemics, there are limited reliable, proven, real-time leading indicators related to human activities that could provide early insights about the risk of disease spread in a region to inform proactive policy making. One reason for this limitation has been the absence of quantitative measures and data that could be examined to proactively evaluate human activities. With advancements in location intelligence data technologies, however, information derived from cellular devices offers a large depository of digital trace data related to human activities have been have increasingly been adapted and analyzed to promote understanding of and to quantify human activity and mobility in pandemic analysis, as well as in other applications [8] [9] [10] . For example, in the context of COVID-19, the radius of gyration, which captures the mobility of individuals using human movement trajectories, was adopted to analyze the COVID-19 spread in Japan [11] . Daily step-counts (gathered from smartphones) were used to estimate and predict decreased movement of individuals within the United States during COVID-19 [3] . Two of the most important aspects of human activities during an epidemic are agglomeration of activities and intensity of activities.
Although previous research reveals insights regarding human activities in the context of COVID-19, the relationship between human activities and disease-spreading risk has not been fully explored, and leading indicators of human activities to proactively assess the risk of disease spread during the early stages of pandemics are lacking. The majority of research studies [12] [13] [14] [15] focus on quantifying and analyzing the changes in human activities as a consequence of the outbreak of the virus and in response to protective policies (such as shelter-in-place policies). The time-lag relationship between these human activity metrics and the spread of virus, which can be generally described by the basic reproduction number (R 0 ), has not yet been fully examined. The basic reproduction number, R 0 , is defined as the number of secondary cases produced by one previous case in a completely susceptible population [16] . Although research studies [17] [18] have focused on leading indicators obtained from users' online search behavior, the decrease of online search frequencies may not have direct impact on the spread of virus. Hence, the previous indicators cannot be utilized for proactive assessment of disease spread risk in a proactive manner.
In addition, the nature of human activities-such as activities in public indoor venues versus in residences-may have differing effects on the spread risk of a disease. The measures of human activities should distinguish between the nature of activities to provide useful leading insights for decision making and policy formulation.
In this study, we adopted the Venables distance (D v ) index [19] and also created the activity density (D a ) index to serve as two real-time indicators to examine the spatial and temporal patterns of human activities across 193 counties in the United States using Mapbox high-resolution temporal-spatial activity index data from January 1 to March 31, 2020. The Venables distance captures the average distance (i.e., concentration) of human activities across a city or county (less distance between persons might indicate a greater contact risk). The activity density captures the intensity level of overall activities in a county or city (higher activity levels might indicate a greater spread risk). Human activities were examined in four categories-social, traffic, work, and other-based on the location and time of activities. Accordingly, we analyzed the correlation between the two metrics (D v and D a ) and the basic reproduction number for 193 counties with the highest number of confirmed COVID-19 cases. The rest of this paper is organized into three sections. The first section discusses the description of the two datasets (Mapbox data and total confirmed cases number data), as well as the analysis methods. The second section describes the results of time-lag correlation analysis between the two metrics and the basic reproduction number. The last section presents the results and the implications of the findings for future work.
In this section, we describe the two datasets-Mapbox data and total confirmed cases number data-and the procedures for human activity categorization. Also covered in this section are definitions and equations related to the Venables distance (D v ), the activity density (D a ), and the basic reproduction number (R 0 ). The time-lag cross-correlation analysis method is presented at the end of this section.
We utilized digital trace telemetry data obtained from Mapbox from January 1 to March 30, 2020. The dataset contains a metric of telemetry-based human activity, a tile,t , which varies across spatial tiles and time t. The partition of tiles is based on Mercantile, a Python library, which is capable of creating spatial-resolution grids all worldwide. The a tile,t is collected, aggregated, and normalized by Mapbox from geography information updates from users' cell phone locations by time flows. The more users located in a tile at time t, the higher the human activity (i.e., a tile,t ). The dataset comprises the United States and the District of Columbia; however, in this study, we examined only 193 counties whose cumulative confirmed cases were greater than 100 as of March 31, 2020 (shown in 1). In the raw data, the temporal resolution is 4 hours. Each tile represents about 100 by 100 square meters for spatial resolution. Since the data is derived from cell phone activity, data may not exist for all cells at all times. For example, a park opens during the day but is closed at night would not generate any data at midnight. Also, for protecting users' privacy and the data aggregation process, tiles with a small number of users are reported without any activity data. It is also noteworthy that the data is aggregated and normalized for each month, so the absolute values of activity indices for different months cannot be directly compared. To reveal the time-lag relationship between metrics and spread of the virus, the total number of confirmed cases was used. We obtained the data from the COVID-19 Data Repository by the Center for Systems Science and Engineering (CSSE) at Johns Hopkins University [20] . The data in this repository were gathered and aggregated from various sources, such as the World Health Organization (WHO) and the Centers for Disease Control and Prevention (CDC). We extracted the total number of confirmed cases c i,j from the CSSE repository, where i represents each date and j represents each county. Figure 1 : 193 selected counties whose cumulative confirmed cases were greater than 100 as of March 31, 2020.
The nature of an activity might put its participants at a higher risk level for contracting the virus. For example, activities in public common areas, such as grocery stores or gyms, would lead to greater risk of disease spread compared to the activities in residential areas, such as working from home or walking a dog in the community. The fine granularity of the spatial resolution enables classification of each tile into one of the four categories: (1) social tiles, (2) traffic tiles, (3) work tiles, and (4) other tiles. Categorization is based on the following characteristics: (1) social tiles are the location of at least one point of interest; (2) traffic tiles include a road; (3) work tiles show no activity in the evening; and (4) other tiles are located in residential areas. The analysis in this paper examines human activities for social, work, and traffic tiles separately. Example tile maps related to each category are shown in 2 for Harris County, Texas. 
To quantify the agglomeration of human activities, we used the Venables distance (D v ) as a weighted average distance of human activities. The Venables Distance aggregates the spatial distribution of a tile,t in a county and captures the urban spatial structure of human activities (Louail et al. 2014 ). The D v is calculated using Equation 1:
where, a tile1,t and a tile2,t are the metrics of human activities in tile 1 and tile 2 at time t, respectively, and d tile1,tile2 ) is the distance from the centroids between these two tiles. In Harris County, Texas, there are more than 70K unique tiles, which makes it computational expensive to analyze all pairs of existing tiles. To reduce the computational burden, we aggregated the 100 by 100-square-meter tiles, a tile,t , to square cells 2 kilometers in length using Equation 2:
where, a k,t is the intensity of human activity in cell k, at time t, A k is the area of the cell k. By aggregating human activity into a larger spatial cell, we reduced the computational efforts, maintaining a meaningful spatial resolution without losing important characteristics in the raw data. Accordingly, the modified Venables distance is derived as shown in Equation 3 :
where, a k1,t and a k2,t are the intensity of human activities in cells k 1 and k 2 at time t, respectively, and d k1,k2 is the distance from the centroids between these two cells. In Equation 3 , the values of the activity intensity (a k,t ) are used as weights to calculate a human activity-weighted distance for the whole area. In other words, the relative values of a k,t were used to examine changes in agglomeration of activities. We calculated D v (t) for each county j, which is denoted as D v (j, t) for all cells in the county j.
Although D v captures the agglomeration of human activities, the density of activities is also critical for examining population contact. To make the raw data (a tile,t ) from Mapbox comparable among different months, we de-normalized the activity index to the contact activity metric ca tile,t for each tile and each month (where t is time). In the denormalization process, the assumption was that, in each month, the minimum human activity intensity among tiles is the same. First, the tile with minimum human activity index (a tile,t ) in each month is found. Then we set the value of contact activity ca tile,t in that low-activity tile as 5, and de-normalized the values for other tiles based on this value. For each county, the activity density at time t (D a (t)) was calculated using Equation 4 :
The basic reproduction number (R 0 ) (the number of secondary cases arising from one previous case) is a critical parameter in epidemic modeling for understanding the speed of disease spread [ 
where, c i and c i+t represent the total confirmed cases in day i and day i + t, respectively; τ is a constant parameter. We estimated R 0 using CDC data (c i,j ) in Equation 6 :
where R 0i,j is the basic reproduction number at date i in county j. Because the CDC total confirmed cases data fluctuates within the course of a week (i.e., more reported cases in weekdays and less reported cases during weekends), the time interval t was set to 7 days. Based on the existing literature and simulation models related to COVID-19 [26] [27] [28] , the constant parameter τ was set to 5.1 days. Accordingly, the R 0 was calculated for all 193 counties for the analysis period.
In the next step, we examined the correlation between the two human activity indicators and the basic reproduction number across all counties. Since these variables are a time series, we used time-lagged cross-correlation analysis to assess the synchrony of time series data sets. The cross-correlation coefficient was calculated using Equation 7 :
where ρ A1,A2 is the cross-correlation coefficient for two time series data A 1 and A 2 ; δ is the time offset of A 2 ; Cov(X, Y ) is the function calculating the covariance of two variables; σ A1 and σ A2 are the standard deviation of data A 1 and A 2 , respectively. Based on the definition, ρ A1,A2 represents the correlation between two variables and |ρ A1,A2 | ≤ 1 ( |ρ A1,A2 | = 1 happens if and only if A 1 = mA 2 + n, where m and n are constants). Then, by iteratively calculating the ρ A1,A2 (δ) with different δ, the correlation coefficient would reach its peak when δ = T , and T was determined as the time lag between two variables.
This section presents the results related to the calculation of the two human activity metrics and their time-lagged correlation with the basic reproduction number across 193 counties during the initial stage of the COVID-19 outbreak in the United States.
In this study, the Venables distance (D v ), and the activity density (D a ) were calculated to assess the human activities at the county level using data from Mapbox. The very first four weeks (January 1 to 28) were considered as the baseline, and D v and D a values were compared with the average baseline in corresponding weekdays. For example, the D v values on March 1 (Sunday) were compared with the mean value of D v values on Sundays between January 1 and 28. Three different ways of daily activity aggregation were used: peak, average, and noon. The peak (largest), average, or the noon (11 a.m. to 3 p.m.) values of human activities a tile,t were selected and set as the representative value of each tile at each day. Figure 3 and Figure 4 show the percentage of D v and D a change for social, traffic, and work activity categories and for three types of daily tile aggregation. The increasing trend of Venables distance (D v ) implies declining concentration and rising distance among people, and the decreasing trend of activity density (D a ) indicates less human activity compared with the beginning of this year. Due to the shelter-in-place policies, residents changed their daily activity patterns. More and more people reduced the non-essential outdoor activities (e.g., shopping in supermarkets, exercising in gyms, and eating at restaurants). Such changes in daily human activity patterns led to the change of D v and D a . For the three categories, significant change can be seen in both social and traffic tiles, while the change in work tiles is not obvious. This is because the activities in work tiles could be essential activities. The D v increased the most, around 15%, in social tiles, while the D a fell the most in traffic tiles, which is around 25%. Differences among the three types of daily tile aggregation were not significant. The percentage change of peak values is slightly greater than the other two values, indicating that peak values are influenced more by COVID-19, while average values are more stable. The following analysis uses peak values to calculate D v and D a . Histograms of average percentage change of D v and D a for each county are shown in Figure 5 . The average percentage change is calculated during the last week of March. The D v values in the majority of counties increased, and D a values decreased for social and traffic categories, while the work category shows more even distribution around 0% for both D v and D a . These histogram plots are consistent with the claim that human activities in work tiles are more essential than other two (social and traffic), which did not show significant change during the COVID-19 study period. While D v and D a describe the different global characteristics of human activity-the D v captures spatial distribution of human activity, and the D a focuses on the intensity of human activity-they all reveal the insight of massive human activity patterns, which could have a quite significant influence to the spread of the virus. The correlation analysis between these two metrics and the basic reproduction number R 0 becomes critical.
The spread of the coronavirus is closely related to the human activity patterns. In the previous section, we showed that the average distance between human activities (D v ) increased by 10% to 15%, and the average human activity intensity (D a ) decreased 5% to 10% for social tiles during March 2020 compared with the baseline period of January 2020. This result provides a good indication of the reduction in human activities in response to social distancing policies.
In this section, we examine the extent to which the change in human activity metrics was related with the change in reproduction rate of corona virus in the 193 counties under study. Hence, we conducted the time-lagged correlation analysis for the two human activity metrics calculated based on social, traffic, and work tile activity categories. Figure 6 shows the time offset results between the Venables distance (D v ) change and the basic reproduction number (R 0 ) for social, traffic, and work activity categories. Since the number of confirmed cases follows a skewed distribution during March 2020, the log scale is used to illustrate the results. The results show that, in majority of typical counties, the decline in the basic reproduction number (R 0 ) happens 20 to 40 days after the increase in Venables distance. This result is consistent across all three activity categories. In the right column of Figure 6 , the bar charts show the correlation between the offset D v and R 0 within different P-value intervals. The average correlation coefficients (with P-values less than 0.05) are around 0.8 for each category, indicating a significant relationship between the increased distance among human activities and the decline in the virus spread speed. For the P-values greater than 0.05 (indicating no sufficient evidence to prove the correlation between two variables), the correlation indices are correspondingly smaller. The number of counties in each P-value interval show that about 50% of the counties have P-values less than 0.001 for Venables Distance calculated based on social and traffic activity tiles. The results related to work tiles, however, show a significant correlation between the two variables in a smaller number of counties. Figure 7 shows the time offset result between activity density (D a ) change and reproduction number (R 0 ) for the three activity categories. The results show that the decline in the basic reproduction number happens 6 to 17 days after the reduction of the activity intensity (D a ); a similar result exists in all activity categories. The time lag is less than the one obtained for the Venables distance (D v ), which means that the spread of virus responds to human activity intensity reduction more quickly than to human agglomeration reduction. In the right column of Figure 7 , the bar charts show the correlation of offset D a and R 0 in different P-value intervals. The average correlation indices for P-values less 
In the next step, we examined the variation of findings across counties with different population sizes, number of confirmed cases, and date of first confirmed cases. The goal is to examine the extent to which the correlation between the two metrics of human activities and the reproduction number is sensitive to these county features. The 193 counties were divided into three uniform categories according to population size and confirmed cases (on March 18, 2020) labeled high, medium, and low. Similarly, the first case dates were labeled as early, mid-range, and late for each one third of counties. Then, the changes in D v and D a (on March 31, 2020) were examined for each label in each tile category, and the results are plotted in Figure 8 . As shown in Figure 8 , for all three tile activity categories, the D a declined more in counties with larger population, more confirmed cases, and earlier first-case date. This result indicates a greater recognition and response to the pandemic risks in more populous counties with early confirmed cases. 
This study shows the utility of two human activity metrics (the Venables distance (D v ) and the activity density (D a )) as leading indicators for the spread speed of COVID-19 in the early stages of the outbreak. These metrics were derived from digital trace data obtained from Mapbox high-resolution temporal-spatial datasets. The results provide statistical evidence regarding the time-lag correlation between these two metrics and the basic reproduction number (R 0 ) in the context of COVID-19. The results regarding the significant leader-follower relationship between human activities and the rate of spread of viral infections could provide valuable implications for authorities to monitor and control the transmission of COVID-19 and future pandemics. For example, time lag indicates that the spread of virus responds to human activity intensity reduction more quickly than to human agglomeration reduction. Hence, the proposed indicators can be calculated using digital trace telemetry data in near real time to proactively assess the risk of virus spread. This study has some limitations which need to be improved in future studies. First, the tile activity categorization-social, work, and traffic-is not precise. One tile could be labeled as both social and work. In this study and due to data availability limitations, however, we classified tiles into only one of the three categories. Second, the CDC confirmed-cases data had limitations due to testing availability. In this study, we did not adjust the confirmed case data based on the extent of testing in different counties. A lack of testing in some areas resulted in the underestimation of the total cases. 