Quantitative Analysis of the Effectiveness of Public Health Measures on COVID-19 Transmission

2.1. Relevant background: early-time dynamic of SIR models in networks

In this section, we present relevant background on the Susceptible-Infectious-Recovered (SIR) model in networks. We refer the reader to [11] for a comprehensive analysis on epidemiological models and to [12] for the seminal paper on the original SIR model. Since we focus on the early-time dynamics of the SIR models, we can assume that the number of births and deaths are much smaller than the population, in a way that the closed population hypothesis holds.

Define as s_i(t), x_i(t), and r_i(t) the share of susceptible, infectious, and recovery persons of city i relative to the local population at time t. In a closed population, the SIR model in networks is government by the following differential equations:

∀i ∈ V and t ≥ 0. We can substitute (4) into (2), yielding:

In early time, i.e., we can assume that x_i(t) ≪ 1 and r_i(t) ≈ 0, ∀i ∈ 𝓋. Therefore, we can ignore second-order x_i(t) terms and effectively set r_i(t) to 0. With these modifications, Equation (5) becomes: in which I is the identity matrix, is the adjacency matrix A with a homogeneous perturbation of in the main diagonal, and δ_ij = 1 if i = j, and δ_ij = 0 otherwise. Equation (6) is a standard differential linear system whose solution can be written in terms of the eigenvector basis of the adjacency matrix A: in which A · v_k = λ_kv_k holds k 1, …,V. The term λ_k is the k-th eigenvalue of A, v_i,k is the i-th entry of the eigenvector associated with the k-th eigenvalue. The parameter a_i,k(0) in (7) is a scaling constant that depends on the initial condition of city i.

In early time, the growth rate of equation (7) is government by the exponent term with the largest eigen-value λ₁ = λ_max of matrix A, which is a well-known measure from spectral graph theory denominated graph spectrum [13]. Therefore:

− i.e., the growth rate is λ_maxβ − γ and the probability of contagion is proportional to the eigenvector associate with the largest eigenvalue λ_max, v₁, which corresponds to the eigenvector centrality measure of the graph, according to the spectral graph theory [13].

2.2. Relative literature

Basically, there are two strategies to prevent epidemic spreading in networks [14]. One is the efficient immunization protocols and the other is to find out relevant spreaders and activation mechanisms.

Immunization strategies are methods for identification of nodes that shall be immunized, taking into account the network structure. Immunized nodes and all the incident links can be removed from the epidemic network. Immunization can not only protect immunized individuals, but can also reduce the epidemic threshold, precluding the outbreak of the disease. Among various immunization strategies, random immunization protocol is the simplest one, where a fraction of randomly selected nodes are made immune. However, in this case, the immunization threshold tends to be 1 in heterogeneous networks, indicating that almost the whole network must be immunized to suppress the disease [15]. Target immunization protocol considers special nodes to be immunized. In [16, 15], the authors show that the immunization threshold can be exponentially small over a large range of the spreading rate if considers the immunization of a fraction of nodes with the largest degree. Other approaches consider not only the critical nodes, but also the entire prevalence curve (the so-called viral conductance) [17, 18].

Although immunization is a fundamental strategy in the epidemic study, the research community pays also much attention to find out which nodes, links and local structures are most influential or effective in the spreading process [19, 20, 21, 22, 23, 24, 25, 26]. These findings aimed at understanding network measures on nodes and links, such as degree, betweenness, K-core index, closeness, link property on spreading dynamics. Besides of finding superspreaders, some researchers also worked on the identification of how topological features influence global epidemics [27, 28, 29].

However, the above mentioned strategies require the discovery of vaccine or at least partial knowledge on the epidemic network under consideration. With the mass and quick spreading of COVID-19, neither of them is a practical method to prevent the outbreak. Therefore, global intervention methods, like social isolation, even lockdown, have already been proven to be efficient. For this reason, we study the effectiveness of public intervention methods. Our results provide strong evidence on the effectiveness of public health measures, such as quarantine and use of masks, to reduce the increasing rate of infection even without detailed information of the highly dynamical population network.

3.1. Network construction using panel data

Consider the weighted directed graph ℊ = 𝓋, ℯ in which 𝓋 is the set of vertices and ℰ is the set of links. There are V = 𝓋 vertices and E = ℰ links in the network. In our epidemiological application, vertices can represent cities, states, countries, or any well-defined entity or geographical circumscription (neighborhood, street, house etc.). For simplicity and with no loss of generality, we denominate the vertices as cities. We assume as given the set of cities/vertices 𝓋. In contrast, links between cities i and j connote potential COVID-19 transmission from i to j and are a priori unknown. In the context of cities, city-to-city contagion could happen for a series of reasons, such as when infectious persons visit or migrate or even from intercity transportation of supplies covered in surfaces that the SARS-CoV-2 is viable for long periods. Therefore, the network 𝒢 encodes all potential transmission paths between cities be through organic or non-organic media. The goal of this section is to estimate the set of links ℰ, i.e., the intercity COVID-19 transmission channels.

Let x(t) = [x₁(t), x₂(t), …, x_V (t)] denote the vector with shares of infectious persons relative to the local population of every city i ∈ 𝒱 in the network at discrete time t ≥ 0. Specifically, we denote as x_i(t) ∈ [0, 1] the share of infectious persons within city i at time t. That is, we take the ratio between the number of infectious persons to the total local population in the city. When x_i(t) = 1, then all population in the city is infectious. When x_i(t) = 0, none is infectious. In-between values represent partial shares of infectious population. Define the column vector x_i = [x_i(0), x_i(1),…, x_i(T)]^′ as the COVID-19 time series evolution in city i up to time T, in which the superscript^′ is the transpose operator. Since we perform an early-time analysis of the epidemics, T is likely to not be large. Let also the matrix X = [x₁, x₂, …, x_V], dim(X) = T × V, be all the cities’ time series with the shares of infectious persons stacked in columns over all period with available data (panel data).

To construct the network, we consider the temporal ordering of the COVID-19 spread across different cities. We attempt to describe the current share of infectious persons vector x_t with the same vector immediately at the previous time step, i.e., x_t−1 as follows:

∀t ∈ {0, 1, …, T. The term κ, dim(κ) = V 1, is an intercept column vector; A, dim(A) = V V, is the adjacency matrix encoding the set of links ℰ of the graph; and ∊_t (0, Σ_∊) is the unobservable zero mean white noise vector process (serially uncorrelated or independent) with time-invariant covariance ∈ matrix Σ_∊. Let A_ij be the (i, j)-entry of A, i, j ∈ 𝒱. When A_ij > 0, then city i can spillover COVID-19 to city j. The larger A_ij is, the stronger is such contagion. Then, the set of links is given by ℰ = {i, j ∈ 𝒱 : A_ij > 0}.

The terms κ, A in Equation (9) are unknown and are estimated using a fitting process to the observed data X. Equation (9) describes a VAR(1) model. To ensure that the system is stable, the companion matrix must have roots inside the complex unit circle. To guarantee such property, our variables x_i, i ∈ 𝒱, must be stationary. Since they are lower- and upper-bounded—i.e., x_i ∈ [0, 1]—then they are stationary by construction. Specifically, we minimize the following regularized loss function L [7] using the coordinate descent algorithm [30]: in which λ ≥ 0 is the elastic net regularization term and α ∈ [0, 1] is the tradeoff parameter between Lasso (L₁ norm) and Ridge (L₂ norm) regularizations. We notate ∥·∥_F, ∥·∥₁, ∥·∥₂ as the Frobenius, L₁, and L₂ norms, respectively. Larger values of λ encourage sparser networks. The first term represents minimization of the error term ∊_t, ∀t ∈ 0, 1, …, T, and ensures that the estimated adjacency matrix A better reflects the COVID-19 transmission dynamics over time. The second term is a regularization term over the adjacency matrix A introduced to prevent overfitting and ensure that the estimation is numerically tractable. We do not regularize the intercept vector κ because it conceptually adapts to the city-specific average values of our data.

There is an empirical challenge in fitting the adjacency matrix A to the panel data X when we are dealing with large-scale networks in which the number of cities V largely surpasses the number of available time points T, i.e., when V ≫ T. Such problem is aggravated when we only have early-time information about the disease, i.e., T is small. In this case, we would incur in overparametrization and overfitting is a concern. The regularization term in (10) mitigates such concern. We opt for an elastic net regularization scheme because it is a robust regularizator that combines positive features of Lasso and Ridge regularizations [30].

Due to the temporal dependency of the panel data, the usual k-fold cross-validation is not well-suited for our model selection procedure. Following [7], we optimize the penalty parameters λ and α in (10) using a h-step ahead mean-square forecast error (MSFE). Due to data availability, we keep h = 1 so as to minimize further data losses. We divide the data into three equally-spaced and contiguous periods: (i) initialization (t ∈ {0,…, T₁}), (ii) training (t ∈, {T₁ + 1, …, T₂}), and (iii) forecast evaluation (t ∈ {T₂ + 1, …, T}), in which and . We also use a rolling validation process as follows. We first fit the model using all data up to time T₁ and forecast , in which λ_c and α_c are fixed candidate penalty terms. We then sequentially add one observation at a time and repeat this process until T₂ 1. Then, we choose the penalty terms λ and α that minimize the one-step ahead MSFE given by:

Finally, we estimate the one-step ahead forecast accuracy using data points in t ∈ {T₂, …, T}, which have not been used in the model selection procedure. To better assess the potentiality of the network in amplifying contagion across different municipalities, we remove the self-loops in the estimated network, which correspond to the influence of the local infectious population on its own future value.

3.2. Estimating transmission rate in early-time epidemics networks

In this section, we assume the network structure 𝒢 = ‹𝒱, ℰ› as given, i.e., the set of vertices and links are already established in accordance with the network construction described in Section 3.1. We start from the results of the early-time dynamic of SIR models in networks described in Section 2.1. Therein, we show that the growth rate at early time is determined by λ_maxβ − γ (see Equation (8)). Therefore, the graph spectrum λ_max modulates the transmission rate parameter by either amplifying or dampening the contagion speed.

If λ_maxβ > γ, then Equation (8) grows exponentially, while it decays when λ_maxβ < γ. Therefore, the reproduction number (critical point) is . Recall that the reproduction number in the SIR model without network is [12]. Therefore, the reproduction numbers of both models differ by the graph spectrum of A, λ_max.

Equation (8) assumes that every city in the model has a single growth rate dynamics dictated by the term λ_maxβ γ. Changes in the epidemics spreading for each city would then be fully determined by their eigenvector centralities, because growth rates are identical across cities (see Equation (8)). However, studies show that the transmission rate parameter β is dependent on local aspects of cities [11]. In contrast, the recovery rate parameter is much less variable across different places. As mentioned earlier in this study, WHO indicates an average recovery time of 14 days for mild cases. Therefore, we consider a different transmission rate for each city in the network β_i while letting fixed the recovery rate γ for all cities. We can still apply the classical framework of SIR in networks because, even though transmission rates are city specific, they tend to be normally distributed around some mean natural value. That is, large deviations are unusual. We empirically find this fact using our application to the Brazilian case. Mathematically, we rewrite (8) as follows:

We can linearize (12) by simply taking the log(.) at both sides of the equation for each city i in the network:

∀i ∈ 𝒱. The LHS and RHS are always non-negative, because x(0) ≥ 0 and is non-decreasing (early-time assumption), e^(λ1βi−γ)t ≥ 0 (asymptotically speaking), and v_1,i ≥ 0 [13]. We can then apply the log(.) without any restrictions. We can estimate (13) for all cities i at once by adding dummies for the constant and time-dependent term for each city in the model (2 dummies per each city). We end up with a set of 2n – 1 dummy variables, because the last one is the reference dummy. Since we have a panel data with temporal dependencies (the same city appears multiple times), we use a linear panel-data estimation model [31] as follows:

∀i ∈ 𝒱, in which α_i and ρ_i are the constant and time-variant dummy terms for city i, and ε_i(t) is the residual from the least square estimation with dummies. We cluster the errors at the city level, such as to mitigate concerns with heteroskedasticity and serial correlation, which could bias our coefficient estimates. Equations (13) and (14) are linked by the following identities:

Given the recovery rate γ—which is assumed to not change over time nor across cities—we can fully identify the eigenvector centrality and the local transmission rate of every city i using (15) and (16), respectively. We only take city-specific estimations of v_1,i and β_i that are statistically significant at the 10% level. Otherwise, we set the estimated coefficients to zero.

3.3. Assessing efficiency of health policy measures in epidemics spreading

With our framework, we can analyze the speed of the epidemics spreading through the network at early time by simply inspecting the graph spectrum λ_max = λ₁ for different time horizons using the methodology described in Section 3.1. Since the reproduction number of the epidemics is proportional to the graph spectrum, then large graph spectra indicate a higher speed of contagion. Any changes of the graph spectrum can be attributed to a “net effect” of public policies of the government in the entire network. Since we use the share of infectious persons of each Brazilian city, then this “net effect” comprises not only federal policies, but also state- and even city-level policies.

Moreover, we can estimate the human impact of these policies in terms of changes in the number of infectious persons at the peak by running the SIR model described in Section 3.2 for each estimated city-specific transmission rate parameter β_i defined in (16) and for different values of the graph spectrum. We use a conservative approach and compare the largest observed graph spectrum with the most recent graph spectrum in our dataset. We assume that the largest graph spectrum occurs when public policies were still latent and were not having effects in the epidemics spreading. Most recent values of the graph spectrum are assumed to represent transmission dynamic after public policies were in, as was the case in Brazil who adopted quarantine and recommended the use of masks in the period that we have available data.

4.1. Data

We use daily data on the number of infectious persons per each city in Brazil using COVID-19 epidemiological bulletins of 27 State Health Departments from February 25, 2020, to May 8, 2020.⁸ Each Brazilian state compiles local reports from cities inside their geographical circumscription. We end up with 60,021 city-time epidemiological bulletins comprising 2,754 (out of 5,570) cities affected by COVID-19 in Brazil.

Our data is representative because local hospitals are required by law to register any COVID-19 events to the local government while cities and states must notify the federal government. However, there may be substantial sub-notifications due to persons that acquire the COVID-19 and recover unnoticed or without hospitalization.

We also collect city-level population estimates in the Brazilian Institute of Geography and Statistics (IBGE), which is the agency responsible for official collection of statistical, geographic, cartographic, geodetic and environmental information in Brazil. We evaluate the share of infectious persons by taking the ratio of COVID-19 cases reported in the local health bulletin and the local population size. The use of shares in our estimation models is important because it is a stationary variable.

We apply a three-day smoothing filter on the number of infectious persons in each municipality to alleviate concerns with late contamination reports or short-term rectifications by the local health government that could compromise our estimations. In our network construction procedure (see Section 3.1), we keep only cities that reported COVID-19 cases in at least 20% of the available time frame. Our results remain qualitatively the same if we do not apply this filtering criterion. In our estimation of the SIR parameters (see Section 3.2), we center all time points in relation to the occurrence of the first death in the city.

Figures 3a–3b portray the COVID-19 evolution in six of the most affected cities in Brazil relatively to the first reported death in terms of the number of COVID-19 cases and as a share of the local population size, respectively. São Paulo (SP) has the highest number of infectious persons. However, there is strong size effect: São Paulo (SP) has almost 12.2 million residents while the second largest city, Rio de Janeiro (RJ), has almost half of that (6.8 million). To get a sense of the local COVID-19 criticality, we can look at its evolution as a share of the local population. In this case, we note that COVID-19 transmission speed is much larger in Manaus (AM) and Fortaleza (CE). Brasília (DF) and Porto Alegre (RS) have smaller transmission rates and local COVID-19 criticality. However, mortality rates may not follow such incidence criticality, because they correlate with local health quality and demography characteristics.

• Download figure
• Open in new tab

Figure 3:

COVID-19 evolution in six of the most affected cities in Brazil (a) in absolute terms (number of infectious persons) and (a) as a share of the local population size. Horizontal axis represent the relative day in terms of the first observed death due to the COVID-19.

4.2. Results

This section presents the main empirical results of the paper. We first build the COVID-19 intercity transmission network and analyze its propensity of amplifying the COVID-19 in different cities. Then, we analyze the net effectiveness of public health measures adopted by the Brazilian government.

4.2.1. Intercity COVID-19 transmission network in Brazil

Figure 4 shows the graph spectrum of the COVID-19 intercity transmission network of Brazil over time. For each time point (horizontal axis), we run the network construction through the fitting process in Section 3.1 with data from the beginning of the sample up to that specific time point. Even though our sample starts in February 25, 2020, we start the fitting process from March 13, 2020, such as to have enough data for the fitting process. That is, we start with 18 time points for each Brazilian city. Therefore, we initially divide the panel data in three equally-sized groups with 6 time points for model training, model selection (parameters and penalty terms), and model evaluation. These group sizes increase as we add more time points. We perform the network construction estimation daily from March 13 to May 8, 2020, in an independent manner.

• Download figure
• Open in new tab

Figure 4:

Graph spectrum of the COVID-19 intercity transmission network of Brazil (see Section 3.1). The shaded area indicating the timing window in which quarantine measures were adopted by the most affected Brazilian states. The red dashed line indices the beginning of the quarantine in São Paulo, the COVID-19 epicenter in Brazil. The blue dashed line indices the beginning of the use of masks recommendation by the Federal Health Ministry. For each time point (horizontal axis), we build the network with city-specific shares of infectious persons with data up to that point.

In Figure 4, we add a shaded area indicating the timing window in which quarantine measures were adopted by the most affected Brazilian states. Since São Paulo is the COVID-19 epicenter in Brazil as it encloses 57.4% of all the COVID-19 infections in Brazil, we also add a vertical dashed red line indicating the beginning of the quarantine adopted by the São Paulo State Government. We also draw the use of masks recommendation beginning date by the Federal Health Ministry in Brazil as a dashed blue line enacted. While quarantine measures are at the state level, the use of masks recommendations goes at the federal level and encompasses all the 5,570 cities and 27 states in Brazil. São Paulo is the most central city in the transmission network. Therefore, it practically shapes the graph spectrum of the intercity transmission network.

We observe a reduction in the growth rate of the graph spectrum after the quarantine measures precisely two days after the measure. However, the growth rate still persisted at positive rates, indicating that the COVID-19 transmission speed kept increasing after such measure, but with a slower pace. After the incubation period following the use of masks recommendation, we observe a drastic change in the graph spectrum. The growth rate changed sign and started to reduce, showing that the set of health policy measures taken by the government was efficient. However, after April 23, 2020, the graph spectrum again started to increase. This can be due to several factors, such as social confusion in following health guidelines in view of the political disarray that Brazil is facing, or even non-compliance with quarantine and use of masks measures. Our model does not permit to have an isolated causal impact of the use of masks recommendation nor of the quarantine measures. However, it enables us to understand how the set of all policy measures affected the COVID-19 transmission rate across cities over time. Combining Figures 2 and 4, it seems that the reduction in the COVID-19 growth rates after the use of masks recommendation was more apparent in cities with relative low social distancing indices. This may be due to the fact that these cities have more potential close human-to-human contact and therefore the use of masks is crucial to detain the COVID-19 transmission.

To understand the topological aspects of the COVID-19 intercity transmission network, Figure 5 plots the PageRank centrality for the top 5 most central cities in each of the five regions in Brazil. We normalize the PageRank with respect to the most central city: São Paulo (SP) on May 8, 2020. As the city centrality becomes higher, the more it contributes to spreading the COVID-19 throughout the network. The top 5 most central cities in the country are the following state capitals (in decreasing order): (i) São Paulo (SP), (ii) Rio de Janeiro (RJ), (iii) Fortaleza (CE), (iv) Recife (PE), and (v) Manaus (AM). These cities all have airports and are strongly interconnected to the remainder of cities in Brazil through roadways and are likely to be the hubs for the COVID-19 spread to other nearby cities in Brazil, especially countryside municipalities. The centrality of São Paulo (SP) in the Southeast monotonically increases over the entire sample. The same roughly occurs with Manaus (AM) in the North, Fortaleza (CE) in the Northeast. Porto Alegre (RS) in the South and Brasília (DF) in the Midwest have the highest centralities in their region but with a negative growth rate in the last days of the sample. Overall, there is a very heterogeneous profile of the city centralities over time, showing the underlying non-trivial patterns in the COVID-19 transmission network.

• Download figure
• Open in new tab

Figure 5:

Evolution of the normalized PageRank centrality measure in the COVID-19 transmission network (see Section 3.1 for the network construction details). We only report the top 5 cities with highest PageRank at each Brazilian region. For each time point (horizontal axis), we build the network with city-specific shares of infectious persons with temporal data up to that point. Each label is composed of the city name followed by its state inside parentheses.

4.2.2. Measuring the human impact of health policy measures to mitigate the COVID-19 propagation

In this section, we run the SIR in networks (see Equations (1)–(4)) with different transmission rate parameters for each city in Brazil, in accordance with (16). We first estimate the city-specific ρ_i using the panel-data information on counts of the share of infectious persons in each city in Brazil via (14).

Then, we estimate the transmission rate parameter β_i of each city i ∈ 𝒱 in Brazil by fixing the recovery rate parameter as γ = 1/14. We use the remaining parameter λ_max—the graph spectrum—to evaluate the effectiveness of the set of health policy measures in detaining the COVID-19 in Brazil. We take as baseline model the graph spectrum reached in April 10, 2020, which is the maximum observed value. We assume that this graph spectrum would have not changed afterwards in case the set of health policy measures were not taken.⁹ We then run several SIR models with the observed graph spectrum values in Figure 4 after the graph spectrum maximum in April, 10, 2020.

Figure 6a shows a comparison of the infectious peaks of the baseline SIR model—i.e., the hypothetical scenario in which health policy measures were not introduced—and the ones with graph spectrum values observed daily after that maximum. The vertical axis shows the relative change in these infectious peaks of the baseline and the observed model day by day, which can be interpreted in terms of the potential share of spared infections at the infectious peak due to the introduction of the set of health policy measures. Since we have data from each city affected by the COVID-19, we plot the median, percentiles 75% (0.25 distant from the median) and 90% (0.40) of this distribution. In the Supplementary Material, we provide the effectiveness of public health policies for each affected municipality in Brazil. In April 10, 2020, the share of spared infections in the epidemics peak is zero, because the baseline model is compared with itself. Then, as we move forward in time and use smaller graph spectrum values, as shown in Figure 4, the potential share of spared infectious increases. The share of spared infectious persons in the epidemics peak reaches a median value 40% lower than that of the baseline model when we use the graph spectrum in April, 24, 2020, suggesting high effectiveness of the quarantine and use of masks health policies. After this point, the share of spared persons decreases— reflecting the increase in the graph spectrum in Figure 4—giving more room for the spread of the COVID-19. The effectiveness of the health policy measures, however, remains positive throughout the entire sample.

• Download figure
• Open in new tab

Figure 6:

Distribution of the efficiency of health policy measures along all affected cities in Brazil over time. We plot the efficiency distribution as a function of (a) time and (b) the city’s distance to the capital within the same state it resides. Since states in Brazil have substantial differences in their sizes, we normalize the city’s distance to the capital to the most distant city within the state.

The first case of the COVID-19 in Brazil was reported in São Paulo (SP) on February 25, 2020. After that, it spread to several Brazilian state capitals probably through air transportation (most of the airports in Brazil are in the state capitals and capitals are far from each other). The epidemics took some time before reaching the first case in countryside cities. Figure 6b displays the distribution of the potential share of spared infections in terms of the city distance to the state capital. Since Brazilian states are very different in size, we normalize the distance to the most distant city within the state. We observe a positive relationship between potential share of spared infectious and distance to the capital, suggesting that health public policies are most effective in cities that are distant from the capital. This may reflect not only the temporal delay of the COVID-19 in reaching the countryside, which puts the local COVID-19 at very early time in these regions, but also demography aspects, such as lower population density, and agricultural economic activities that do not require large conglomerates of persons.

Figure 7a shows the effectiveness of the set of public health measures for six of the most affected Brazilian capitals. In particular, Brasília (DF) reaches a 50% lower share of infectious persons at the peak when we compare peaks reached with the graph spectrum value on April 24, 2020 (against the baseline in April 10, 2020). Figure 7b shows the number of potential spared infectious persons due to the set of health policy measures. This figure is constructed by simply multiplying the share of spared infectious with the local population size of each of the six cities. Since São Paulo (SP) is the largest city, it would potentially spare more persons when the COVID-19 epidemics reach its peak.

• Download figure
• Open in new tab

Figure 7:

Efficiency of public health measures over time as a function of (a) the share of the spared local population and (b) the spared number of persons (in millions). We depict curves only for six capital cities that are being substantially affected by the COVID-19: Beleém (PA), Fortaleza (CE), Rio de Janeiro (RJ), Brasília (DF), Manaus (AM), and São Paulo (SP).