The ANMs are used for identifying causal effect of a factor or an intervention measure on the number of new cases or an intervention measure [29, 34]. Assume no confounding, no selection bias, and no feedback. Let Y be the average number of new cases or new deaths from COVID-19 and X be a scalar factor or an intervention measure such as gender, population density, ethnic group, among others. Consider a bivariate additive noise model X → Y where Y is a nonlinear function of X and independent additive noise E Y : Y = f Y ( X ) + E Y ,       X ∐ E Y , (1) where X and E Y are independent. Then, the density P X , Y is said to be induced by the additive noise models (ANM) from X to Y [35]. In some cases, we may have the following alternative direction of the ANMs: Y → X : X = f X ( Y ) + E X ,       Y ∐ E X , (2) where Y and E X are independent. If the density P X , Y is induced by the ANM X → Y , but not by the ANM Y → X , then the ANM X → Y is identifiable.

Assume that n + m state data were sampled. Divide the dataset into a training data set by specifying D 1 = { Y n ,   X n } ,   Y n = [ y 1 ,   … ,     y n ] T   ,   X n = [ x 1 , … ,   x n ] T   for fitting the model and a test data set D 2 = { Y ˜ m ,   X ˜ m } ,   Y ˜ m = [ y ˜ 1 , … ,   y ˜ m ] T ,     X ˜ m = [ x ˜ 1 , … , x ˜ m   ] T for testing the independence, where n is not necessarily equal to m.

Procedures for using the ANM to assess causal relationships between two variables are summarized below [34].

Step 1. Regress Y on X using the training dataset D 1 and non-parametric regression methods:

There are many non-parametric methods that can be used to regress Y on X or regress X on Y. For example, we can use neural networks [36], smoothing spline regression methods [37], B-spline [38] and local polynomial regression [39]. In this paper, the smoothing spline regression method was used to fit the regression models.

Covariance can be used to measure association but cannot be used to test independence between two variables with a non-Gaussian distribution (https://en.wikipedia.org/wiki/Correlation_and_dependence). A covariance operator that is a generalization of the finite dimensional covariance matrix to infinite dimensional feature space can be used to test for independence between two variables with arbitrary distributions. Specifically, we will use the Hilbert-Schmidt norm of the cross-covariance operator or its approximation, the Hilbert-Schmidt independence criterion (HSIC) to measure the degree of dependence between the residuals and potential causal variable and test for their independence [35, 40].

The covariance operator can be defined as C o v   ( f ( X ) , g ( Y ) ) = = 〈 f , C ( X , Y ) g 〉   ℱ , where f , g are any nonlinear functions and C ( X , Y ) is the covariance operator and 〈 . 〉 ℱ is an inner product in the ℱ Hilbert space. The covariance operator is positive symmetric, which implies linearity and continuity. In addition, the covariance operator maps spaces in their dual spaces [41]. The Hilbert-Schmidt norm of the covariance operator can be used as criterion for assessing independence between two random variables and is called the Hilbert-Schmidt independence criterion (HSIC). The Hilbert-Schmidt norm of the centered covariance operator is defined as H S I C 2 ( X , Y ) = ‖ C X Y ‖ H S 2 , where ‖ . ‖ H S is the Hilbert-Schmidt norm.

We know that (Wang et al., 2018). H S I C 2 ( X , Y ) = 0 if and only if X and Y are independent. H S I C 2 ( X , Y ) can be approximated by. H S I C 2 ( X , Y ) = 1 n 2 tr ( K H L H ) , where n is a sample size, K and L = are n × n dimensional kernel matrices and H = I − 1 n 1 n T n . We used the Gaussian kernel: k ( x , y ) = e − x − y 2 2 σ 2 ,   σ > 0 and polynomial kernel. The order- d   polynomial kernel is defined as k ( x , y ) = ( x T y + c ) d . The fourth and sixth order polynomial kernels were used in this analysis. To test independence between the potential cause X and residual E Y , we calculated H S I C 2 ( X , E Y ) as follows. H S I C 2 ( X , E Y ) = 1 n 2 tr ( K X H K E Y ) , K E Y = [ k E Y ( ε 1 ,   ε 1 ) ⋯ k E Y ( ε 1 ,   ε n ) ⋮ ⋮ ⋮ k E Y ( ε n ,   ε 1 ) ⋯ k E Y ( ε n ,   ε n ) ] K X = [ k X ( x 1 ,   x 1 ) ⋯ k X ( x 1 , x n ) ⋮ ⋮ ⋮ k X ( x n , x 1 ) ⋯ k X ( x n , x n ) ] , ε i = E ^ Y ( i ) = Y i − f ^ Y ( X i ) ,   i = 1 , … , n ,

In summary, the general procedure for testing independence between the average number of new cases or new deaths and the scalar factor or intervention measure is given as follows [34, 35]:

Step 1: Divide a data set into a training data set D n = { Y n ,   X n } for fitting the model and a test data set D m = { Y ˜ n ,   X ˜ n } for testing the independence.

We do not have closed analytical forms for the asymptotic null distribution of the HSIC and hence it is difficult to calculate the p-values of the independence tests. To solve this problem, the permutation/bootstrap approach can be used to calculate the p-values of the causal test statistics. The null hypothesis is H 0 : no causations X → Y and Y → X (Both X and E Y are dependent, and Y and E X are dependent).

Calculate the test statistic T C = | H S I C 2 ( E Y , X ) − H S I C 2 ( E X , Y ) | . (6)

Assume that the total number of permutations is n p . For each permutation, we fix x i ,   i = 1 , … ,   m and randomly permutate y i ,   i = 1 ,   … ,   m . Then, fit the ANMs and calculate the residuals E X i ,   E Y i ,   i = 1 , … ,   n and test statistic T c . Repeat above procedures n p times. The p-values are defined as the proportions of the statistic T ˜ C (computed on the permuted data) greater than or equal to T ^ C (computed on the original test data D m ).

Each state was a sample. Since the sample sizes were small (only 50), the p-value for declaring significance was 0.05 without Bonferroni correction for multiple comparison.

Before performing multivariate linear Grander causality test, we first need to transform nonstationary time series to stationary time series.

Consider an m -variable VAR with p lags: Y t = μ + ∑ i = 1 p A i Y t − i + ε t , (7) where Y t is a m dimensional vector, μ is a mean, the A i   ( i = 1 , … , p ) are m × m coefficient matrices and m dimensional residual vector ε t is assumed to have mean zero ( E [ ε t ] = 0 ) , with no autocorrelation ( E [ ε t ε t − s T ] = 0 ) , but can be correlated across equations ( E [ ε t ε t T ] = Σ ) .

Vector error correction model (VECM) consists of first differences of cointegrated I ( 1 ) variables, their lags, and error correction terms: Δ Y t = μ + Π Y t − 1 + ∑ i = 1 p − 1 Φ i Δ Y t − i + ε t , (8) where matrixes Π and Φ i   ( i = 1 , … , p − 1 ) are functions of matrices A i   ( i = 1 , … , p ) .

When two non-stationary variables are cointegrated, the VAR model should be augmented with an error correction term for testing the Granger causality (Engle and Granger, 1987).

The VECM can be reduced to Δ Y t = A 0 + A ( L ) Δ Y t − 1 + α ( E C M t − 1 ) + ε t , (9) where Δ Y t = [ Δ Y 1 t ⋮ Δ Y m t ] ,   A 0 = [ A 10 ⋮ A m 0 ] ,   A ( L ) = [ A 11 ( L ) ⋯ A 1 m ( L ) ⋮ ⋮ ⋮ A m 1 ( L ) ⋯ A m m ( L ) ] ,   α = [ α 11 ⋯ α 1 k ⋮ ⋮ ⋮ α m 1 ⋯ α m k ]   , E C M t − 1 = β T Y t − 1 + C t ,   Y t − 1 = [ Y 1 , t − 1 ⋮ Y m , t − 1 ] ,   β = [ β 11 ⋯ β 1 k ⋮ ⋮ ⋮ β m 1 ⋯ β m k ] ,   C = [ c 1 ⋮ c m ] .

Consider two non-stationary time series, X t = [ X 1 t , … ,   X m , t ] T and Y t = [ Y 1 t , … ,   Y m t ] T . Let m = m 1 + m 2 .

Suppose that X t and Y t are cointegrated with the residuals E C M t . The VECM model for testing the Granger causality is given by [ Δ X t Δ Y t ] = [ A x A y ] + [ A x x ( L ) A x y ( L ) A y x ( L ) A y y ( L ) ] [ Δ X t − 1 Δ Y t − 1 ] + [ α x α y ] E C M t − 1 + [ ε x ε y ] (10) where A x and A y are m 1 and m 2 dimensional vectors of intercept terms, respectively, A x x ( L ) ,   A x y ( L ) ,   A y x ( L ) and A y y ( L ) are n 1 × n 1 ,     n 1 × n 2 ,     n 2 × n 1 and n 2 × n 2 dimensional matrices of lag polynomials, respectively, α x and α y are n 1 and n 2 dimensional coefficient vectors for the error correction term E C M t − 1 , respectively. The lag length was selected using the two-stage procedure [42].

There are four different cases of causal relationships between two vectors of time series X t and Y t [33].

(1) If A x y ( L ) is significantly different from the zero, while A y x ( L ) shows no significantly different from zero, then there exists a unidirectional Ganger causality from time series Y t to X t ;

The four statements imply that Ganger causal relationships between X t and Y t depend on the coefficients A x y ( L ) and A y x ( L ) . Therefore, the null hypotheses for testing the Ganger causality between X t and Y t are.

(1) H 0 1 :   A   x y ( L ) = 0 ,

Likelihood ratio tests for multivariate Granger causality are given by.

(1) The likelihood ration statistics for testing the null hypothesis: H 0 1 :   A   x y ( L ) = 0 is

The total number of variables to be tested was 18. The p-value for declaring significance after Bonferroni correction was 0.0028.

Data on the number of new cases and new deaths of COVID-19 across the 50 states in the US were obtained from John Hopkins Coronavirus Resource Center (https://coronavirus.jhu.edu/MAP.HTML). Google mobility indexes were downloaded from Google COVID-19 Community Mobility Reports (https://www.google.com/covid19/mobility/). Comprehensive data and insights on COVID-19’s impact on mobility, economy, and society were downloaded from the University of Maryland COVID-19 Impact Analysis Platform (https://data.covid.umd.edu) [43, 44]. All data were collected from March 5, 2020 to August 25, 2020.

The scalar variables tested for causation of the new cases and deaths from COVID-19 in the US included the number of contact tracing workers per 100,000 people, percent of population above 60 years of age, median income, population density, percentage of African Americans, percentage of Hispanic Americans, percentage of males, employment density, number of points of interests for crowd gathering per 1,000 people, number of staffed hospital beds per 1,000 people, and number of ICU beds per 1,000 people. The number of new cases and deaths were averaged over time. Each state was a sample. Since the sample sizes were small, the p-value for declaring significance was 0.05 without Bonferroni correction for multiple comparison. The p-values for testing 11 scalar potential causes of the number of new cases and deaths from COVID-19 in the US using both Gaussian kernel and polynomial kernels were summarized in Table 1 where minimum of three p-values using Gaussian kernel and fourth order and sixth order polynomial kernels were listed and only one causal direction was observed. We observed from Table 1 that population density (minimum of p-value < 0.0002, which was due to Gaussian kernel), percentage of males (minimum of p-value < 0.03, which was due to Gaussian kernel) and Percentage of Hispanic Americans (minimum of p-value < 0.0325, which was due to sixth order polynomial kernel) showed significant evidence of causing the spread of COVID-19. Employment density (minimum of p-values < 0.0223, which was due to sixth polynomial kernel), Percentage of African American (minimum of p-values < 0.024, which is due to Gaussian kernel), population density (minimum of p-values < 0.025, which was due to Gaussian kernel) and Percentage of males (minimum of p-values <0.0377, which was due to fourth order polynomial) showed significant evidence of causing deaths due to COVID-19. percentage of Hispanic Americans (minimum of p-value < 0.052, which was due to sixth order polynomial kernel) were close to significance level 0.05 for causing death.

Population density was an important risk factor for both the spread and death from COVID-19. High density resulted in closer contact, stronger interaction among residents and lower social distancing, which facilitated the spread and increased the death rate from COVID-19 [45 – 48]. However, our results were contradictory with the conclusion of Hamidi et al. [49]. Some literature also confirmed that high proportion of African Americans caused a high rate of deaths [48, 50, 51]. Our results concluded that percentage of Hispanic Americans was a risk factor for the spread and a weak risk factor for death from COVID-19, while the literature showed stronger evidence that Hispanic communities were highly vulnerable to COVID-19 [52].

The second most significant demographic risk factor for the spread of COVID-19 was percentage of males. We found higher COVID-19 morbidity in males than females. However, we did not find higher COVID-19 mortality in males than females.

It was reported that higher COVID-19 mortality in males than females can be due to the following factors [53]. The first factor was higher expression of angiotensin-converting enzyme-2 (ACE two; receptors for coronavirus) in males than females. The second factor was sex-based immunological differences due to sex hormone and the X chromosome.

Daily mobility and social distancing data from a COVID-19 impacted the analysis platform, including four categories: category A: mobility and social distancing, category B: COVID and health, category C: economic impact, and category D: vulnerable population. A total of 12 temporal metrics in four categories and 12 metrics from the COVID-19 impact analysis platform, six daily Google Community Mobility indexes and protest attendee data that captured real-time trends in movement patterns for each state in the US were included in the analysis to test for Granger causality between these risk factors, health intervention measures and the number of new cases and deaths from COVID-19 across 50 states in the US [44, 54]. The total number of variables to be tested was 18. The p-value for declaring significance after Bonferroni correction was 0.0028.

All 18 metrics except for protest attendee showed high significance in causing a reduction of the new cases of COVID-19 in 19 less affected states: VT, WY, ME, AK, NH, WV, ND, SD, NM, RI, DE, KY, KS, CT, CO, IA, WA, WI, and MS. Most of these states were less populated. However, although CA was most affected and the most populated state, all 18 metrics except protest attendance showed a strong significance in causing rapid spread of COVID-19 (Table 2 and Supplementary Table S1). To provide complete causal testing information, we listed the values of statistics for testing 18 temporal potential causes of the number of new cases of COVID-19 across 50 states in the US in Supplementary Table S2.

All 18 metrics showed no significance in causing reduction of the new cases of COVID-19 in Florida. The majority of the 18 metrics did not demonstrate evidence that they can significantly mitigate the spread of COVID-19 in most the affected states such as TX, NY, GA, IL, AZ, NJ, NC, and TN. These 10 states were in the top largest states by population in the US. Public health intervention measures such as closing schools and businesses, avoiding public gatherings, restricting traffic, placing residents to stay-at-home and adherence to guidelines were less well implemented or difficult to implement homogeneously due to large populations and geographical areas [55]. These results also explained why the number of new cases of COVID-19 in these states was high and confirmed by several studies [56 – 58].

Table 3 summarized the ranges of p-values, Supplementary Tables S3, S4 summarized all p-values and values of statistics for testing 18 temporal potential causes of the number of new deaths from COVID-19 across 50 states in the US, respectively. All 18 metrics except for protest attendance showed high significant evidence for causing a reduction of new deaths across 50 states except for Michigan (MI) in the US. Our results suggested that a cascade of causes led to the COVID-19 tragedy in the US.

Table 4 listed the most significant risk factor for the new cases of COVID-19 in each of the 50 states in the US. Active Cases/1,000 People, workplaces, number of tests completed/1,000 people, imported COVD cases, unemployment rate and unemployment claims/1,000 people, mobility trends for places of residence (residential), retail and recreation, mobility trends for places like restaurants, cafes, shopping centers, theme parks, museums, libraries, and movie theaters (retail) and test capacity were the most significant risk factors for the new cases of COVID-19 in 23, 7, 6, 5, 4, 2, 1 and 1 states of the US, respectively.

Table 5 summarized the most significant risk factor for the deaths from COVID-19 in each of the 50 states in the US. Active Cases/1,000 people, workplaces, residential, unemployment rate, imported COVID cases, unemployment claims/1,000 people, transit, test done/1,000 people, grocery, testing capacity, retail, percentage of change in consumption, percentage of working from home were the most significant risk factor for the deaths of COVID-19 in 17, 10, 4, 4, 3, 2, 2, 2, 1, 1, 1, 1 states, respectively. We also observed that the number of protest attendees showed mild significant evidence to cause increasing the number of new cases of COVID-19 in KY (p-value < 0.00012), KS (p-value < 0.00026), NH (p-value < 0.00108), MA (p-value < 0.0016) and TN (p-value < 0.0024) or to cause more deaths from COVID-19 in OR (p-value < 5.11 E-05), TX (p-value < 0.00017), ME (p-value < 0.00028), KS (p-value < 0.00061), MI (p-value < 0.0015), OH (p-value < 0.0021) and NC (p-value < 0.0023).

In some cases, two causal directions may occur. To examined whether the COVID-19 cases and deaths caused daily mobility and social distancing metrics to change, we summarized the values of statistics and p-values for testing COVID-19 new case and death potential causes of six Google mobility indexes across 50 states in the US in Supplementary Tables S5–7 and DS8, respectively. We did observe the opposite causal direction. Table 6 presented the smallest p-value across 50 states in the US for testing two causal directions: from Google mobility indexes to the number of new cases of COVID-19 and from the number of new cases of COVID-19 to the Google mobility indexes. We observed significant causation from the number of COVID-19to the Google mobility indexes in some states. However, the p-values for testing causation from the number of new cases of COVID-19 to the Google mobility indexes were much larger than that from the Google mobility indexes to the number of new cases of COVID-19. The causal pattern for the COVID-19 deaths was similar to the number of new cases of COVID-19.

To illustrate the causal relationships between the risk factors and the number of new cases and deaths from COVID-19, we plotted Figures 1 and 2. Figure 1 plotted the social distance index curves as a function of time from March 5, 2020 to August 25, 2020 in Florida (FL) and Rhode Island (RI). Figure 1 showed that the social distance index in FL was much higher than that in RI state, which resulted in the larger number of new cases of COVID-19 in FL than that in RI. Figure 2 showed the number of imported COVID-19 cases as a function of time from March 5, 2020 to August 25, 2020 in Maryland (MD) and Wyoming (WY). We observed a huge difference in the number of imported COVID-19 cases between MD and RI. The very low number of imported cases of COVID-19 in WY resulted in the very low number of deaths from COVID-19 in WY, while the high number of imported cases in MD state led to the increased deaths from COVID-19 in MD. These results were consistent with the finding in the literature. It was reported that strong interventions would substantially decrease the number of deaths (Davies et al. 2020; Gagnon et al. 2020).

The COVID-19 case curves as a function of time from March 5, 2020 to August 25, 2020 in FL, RI, MD and WY were plotted in Figure 3, respectively. The response of COVID-19 to public health intervention which were measured by Google mobility indexes were usually delayed.