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 T , 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 T , t is collected, aggregated, and normalized by Mapbox from geography information updates of 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 T , 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 (Figure 1). In the raw data, the temporal resolution is 4 h. 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 open during the day but 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 (2020). 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.

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 location; the location information of point of interest is extracted from SafeGraph data (including restaurant, gas station, and commercial complex) (SafeGraph, 2020); (2) traffic tiles are extracted by mapping the traffic network with all tiles including roads; (3) work tiles are identified based on lack of activity during the late evening hours; and (4) other tiles are located in residential areas. We assigned these four tags to each tile from social, traffic, work to other. Once a tile is assigned with a tag, it is excluded for further categorization. We categorized tiles in this way to examine the importance of human activity in each category and its relationship with the reproduction number. The analysis in this research examines human activities for social, work, and traffic tiles separately. All the residential tiles are excluded from the analysis because the human activities in these tiles have less influence on the contact level among people. Example tile maps related to each category are shown in Figure 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 T , t in a county and captures the urban spatial structure of human activities (Louail et al., 2014). The Venables distance is a metric to describe how distant the human activities are from each other. It also provides a metric for concentration activities. In examining the risk of pandemic, understanding the concentration of human activities is an important consideration. Hence, the Venables distance provides a metric to examine fluctuations in this aspect of human activities. Compared with the well-known spatial autocorrelation metric, the meaning of D v is easier to be explained and understood. Also, the meaning of D v is different from that of spatial autocorrelation, which describes the presence of systematic in spatial space. The value of spatial autocorrelation is from −1 to 1, and it could be negative if the adjacent observations tend to have very contrasting values. The Venables distance, however, is a metric to describe how far the human activities are apart from each other, and its unit, kilometer, has reality meaning. The D v is calculated using Eq. 1: D v ( t ) = ∑ T 1 ≠ T 2 a T 1 , t ⋅ a T 2 , t ⋅ d T 1 , T 2 ∑ T 1 ≠ T 2 a T 1 , t ⋅ a T 2 , t , (1) where a T 1 , t and a T 2 , t are the metrics of human activities in Tile 1 and Tile 2 at time t , respectively, and d T 1 , T 2 is the distance from the centroids between these two tiles. In Harris County, Texas, there are more than 70K unique tiles, which makes it computationally expensive to analyze all pairs of existing tiles. To reduce the computational burden, we aggregated the 100-by-100 square-meter tiles, a T , t , to square cells 2 km in length using Eq. 2: a k , t =   ∑ for all  T  in cell  k a T , t A k , (2) where a k , t is the intensity of human activity in cell k , at time t and 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 Eq. 3: D v ( t ) = ∑ k 1 ≠ k 2 a k 1 , t ⋅ a k 2 , t ⋅ d k 1 , k 2 ∑ k 1 ≠ k 2 a k 1 , t ⋅ a k 2 , t , (3) where a k 1 , t and a k 2 , t are the intensity of human activities in cells k 1 and k 2 at time t , respectively, and d k 1 , k 2 is the distance from the centroids between these two cells. In Eq. 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 . Due to high computational cost of conducting, we aggregated the data into 2 km square. The Venables distance describes how human activities are apart from each other, and the spatial resolution of 2 km is capable to capture such feature.

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 (Dietz, 1993; Gatto et al., 2020), as well as the risk of virus spread (Newman, 2002; Liu et al., 2018; Aleta et al., 2020; Giordano et al., 2020). We adopted the R0 in the analysis, assuming a complete susceptibility of the population. This is because, at the early stage of the COVID-19 pandemic (March, 2020), all people were susceptible. The effective reproduction number ( R t ) is another metric which is widely used to describe the spread speed of disease (Nishiura and Chowell,, 2009; Althaus 2014). R t assumes only partial population is susceptible and takes into account the effect of policies such as social distancing in measuring the disease spread. Thus, R t can be used in cases with populations having immune members (Delamater et al., 2019). Hence, R t would be more appropriate in analyses at the later stage of the pandemic. In the early stage of COVID-19, it was a novel virus for all people and everyone was susceptible. Hence the basic reproduction number ( R 0 ) was used in this study to represent the spread speed of the virus. R 0 is calculated using Eq. 5: c i + t = c i ⋅ R 0 t / τ , (5) where c i and c i + t represent the total confirmed cases in day i and day i + t , respectively, and τ is a constant parameter. We estimated R 0 using CDC data ( c i , j ) in Eq. 6: R 0 i , j = e τ ln ⁡ c i , j − ln ⁡ c i − t , j t , (6) where R 0 i , 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 (Zhang et al., 2020a; Zhang et al., 2020b; Fan et al., 2020), the constant parameter τ was set to 5.1 days. Accordingly, the R 0 was calculated for all 193 counties for the analysis period. The daily basic reproduction number for each county will be used to conduct the correlation analysis with proposed two indicators of human activities.

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 Eq. 7: ρ A 1 A 2 ( τ ) =   Cov ( A 1 ( t ) , A 2 ( t + δ ) ) σ A 1 σ A 2 , (7) where ρ A 1 A 2 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; σ A 1 and σ A 2 are the standard deviation of data A 1 and A 2 , respectively. Based on the definition, ρ A 1 A 2 represents the correlation between two variables and | ρ A 1 A 2 | ≤ 1 ( | ρ A 1 A 2 | = 1 happens if and only if A 1 = m A 2 + n , where m and n are constants). Then, by iteratively calculating the ρ A 1 A 2 ( δ ) with different δ , the correlation coefficient would reach its peak when δ = T , and T was determined as the time lag between two variables.

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–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.–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 nonessential 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 the 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 both reveal the insight of massive human activity patterns, which could have a quite significant influence on the spread of the virus. The correlation analysis between these two metrics and the basic reproduction number R 0 becomes critical.

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. D v and D a show different heterogeneities facing these county features. In the left side of Figure 8, we cannot detect a clear difference in D v value changes in different group. In the right side of Figure 8, however, we can observe that D a values decreased more in the counties with larger population size, higher confirmed cases number, and earlier first confirmed case, and such pattern exists for all the social, traffic, and work tiles. Based on the definition, D a represents the intensity of human activities. Based on this result shown in Figure 8, we can conclude that the intensity of human activities decreased more in counties with larger population size, higher confirmed cases number, and earlier first confirmed cases, which could imply a greater recognition and response to the pandemic risks in such counties.