care interventions. One that caused significant controversies is the use of masks among the general public (1) (2) (3) . Initially, it was assumed that the primary mode of transmission of SARS-CoV-2 was through coughing or contact with surfaces. Extensive shortages of personal protective equipment for health workers led to an initial recommendation for the general population not to wear masks. While no extant studies show surgical masks reduce transmission of SARS-CoV-2 in humans, surgical masks do reduce viral shedding of other coronaviruses (4) , and transmission of SARS-CoV-2 in animal models (5) . Also, cloth masks might filter SARS-CoV-2 (6).
The Centers for Disease Control and Prevention (CDC) changed its guidelines on April 3, 2020 , and recommended the widespread use of masks (7) . According to the CDC, the rationale behind its policy change was the increase in evidence that asymptomatic and presymptomatic people are infectious (8) (9) (10) (11) (12) (13) and that there are many undetected cases (14) . On June 5, the World Health Organization (WHO) changed its guidelines and recommended governments to encourage the general public to wear masks in specific situations, like grocery stores (15) . However, WHO specifies: "At present, there is no direct evidence (from studies on COVID19 and in healthy people in the community) on the effectiveness of universal masking of healthy people in the community to prevent infection with respiratory viruses, including COVID-19" (15) . Meanwhile, according to The European Centre for Disease Prevention and Control (ECDC), there is no evidence that masks can prevent people from being infected (16) , and they suggest the use of masks may increase the risk of infection. The reason for the WHO and ECDC statements is that there is no work on the widespread use of masks by the general public as a non-pharmaceutical intervention (NPI). In this paper, we aim to fill this gap by testing the hypothesis that the policy change regarding masks by the CDC (and local governments) decreased the number of positive cases in the states of Connecticut (CT), Massachusetts (MA), New York (NY), Rhode Island (RI), and Virginia (VA).
To assess causality, we need to evaluate both branches of an intervention (17) : one in which the intervention did happen, and one in which it did not. The gold standard to do so is the double-blind randomized control trial (RCT) paradigm. Although possible (18) , RTCs are not the norm in public health epidemiological intervention. Even if implemented, there is no such thing as a "placebo arm" for travel restrictions or a double-blind school closure. Since a placebo or double-blind trials are not possible, there is an indirect causal path between the treatment and the outcome (19, 20) . For example, people in zip codes with open schools might be more careful with their hygiene because they know that they are at a higher risk than people in zip codes where the schools are closed. When the second branch of the intervention did not happen, it is called a counterfactual ("contrary to the facts"). One option to measure the direct effect of an intervention, and the one that we use in this work, is to estimate or simulate the counterfactual branch (21, 22) .
Here, we present a framework to analyze data from the COVID-19 epidemic that can simulate counterfactual scenarios in which specific NPI did not occur. In this framework, we use the odds of a positive test as the dependent variable, rather than the number of positive tests (23) or deaths (24) . We motivate a linear equation for the evolution of this dependent variable using the Susceptible-Infected-Recovered (SIR) model (25, 26) . Finally, we carry out a LASSO (27) regression to fit the data, to obtain a piecewise-linear fit to the logarithm of the odds with the smallest number of breaks (see the supplementary material). This regression finds the times when interventions started, allowing us to simulate alternative scenarios where these interventions did not happen and assess their net impact.
The daily number of cases and tests are highly variable (Fig. 5 ). To reduce this variability, we compute the log-odds, the logarithm of the number of positive tests over the number of negative ones. Doing this, reveals a piecewise linear pattern that we fit using the LASSO regression (Fig. 6 ). These regressions show three breakpoints in the log-odds in NY, two in CT, MA, MI, and RI, and one in VA. We should stress the fact that these breaks are not an input of the user. On the contrary, this is the result of applying the LASSO regularization. These changepoints happen after different NPIs. The first change in CT and MA, and the first and second in NY, are due to mobility restrictions (school closure, ban mass gatherings, restriction the non-essential workforce, and stay-at-home orders). In these states, the last break happens after the CDC changed its guidelines regarding masks. In MI, RI, and VA, the stay at home orders and the CDC recommendation to wear masks happened closer in time, making it hard to disentangle their effects. However, the MI government only enforced the use of masks in closed areas (like grocery stores), and the VA government never recommended the use of masks. We assume that the lack of local orders correlates negatively with local compliance with the CDC guidelines, and that explains why we do not see the masking effect on MI and VA.
We use the slopes of the regression to compute the R t , the instantaneous reproduction number (28) . This parameter expresses the average number of people infected from one positive case on day t. In Fig. 7 , we show R t as a function of time. The NY plot, in Fig. 7 , shows R t dropping down from 2.1 to 1.6 and then to 0.72 on March 30, 8 days after the closure of all nonessential business. Taking them together, this translates to a reduction by 65 % on R t due to mobility restrictions. There is a third drop from 0.72 to 0.44 on April 14, 11 days after the CDC changed their guidelines and recommended to wear masks, and two days after NY enforced the use of masks for public employees. In CT, the stay-at-home orders reduced the value of R t by 51%. Moreover, after the new CDC recommendation on masks, it dropped by 40%. Remarkably, in MA, after the stay at home order the R t value dropped from 1.9 to 1.1, still above 1. Only after the recommendation of wearing masks it fell to 0.66, below 1. As we already mentioned, we do not see the effect of masks in MI or VA, and we attribute this to the lack of local compliance. In MI, the government only enforced the use of masks in enclosed areas, and the VA government never ruled on the use of masks. The data from RI is harder to interpret because stay-at-home orders and masks guidelines happened close in time, and data from before April are unreliable (with less than 500 tests a day). Nonetheless, there is an effect from the CDC guidelines and from local governments making masks mandatory.
Finally, one advantage of the sparsifying framework is that we can simulate counterfactual scenarios by removing a breakpoint. Then, the regressed line would have continued at the previous slope. Take the case of the public wearing masks. From the fit shown in Fig. 4 , we observe that on April 14 in NY, R t changed from 0.72 to 0.44. We interpret that the counterfactual to this intervention is that if the public had not used masks, R t would have stayed at 0.72, or in causal inference jargon do(no masks). cases (95% CI 65000 to 70000) in MA between Apr. 13 and May 19, and by 29000 cases (95% CI 27000 to 31000) in CT from Apr. 14 to May 17.
In conclusion, we found that masks reduced the spread of the virus in CT, MA, and NY.
In those states, our calculations showed that the intervention reduced the R t by 40%, and we estimate that masking prevented 170000 cases (95% CI 160000 to 180000) from the moment they were adopted until the end of the stay-at-home orders. These results are consistent with recently published results by Mitze and colleagues that found the same effect in Germany (29) .
Also, we estimate that in New York City alone, masks reduced the number of cases by 17900 (95% CI 15800 to 20000) between April 17 and May 9. This number is below the one estimated by Zhang and colleagues who estimated a reduction in 66000 cases during the same period in New York City (30) . We believe the discrepancy can be accounted for by noting that Zhang and colleagues did not consider the increase in testing during the stay-at-home order in, leading to a higher difference after the masking order at which point the testing was more stable.
The framework that we presented is data-driven, and therefore it relies on only a handful of hypotheses as compared to other methods For example, the counterfactual analysis relies on one hypothesis: the log-odds are piecewise linear (see Eq. 6 in the supplementary material)without the need to assume any of the hypotheses of the SIR model. Based on the goodness of fit, we are confident that this hypothesis holds for the datasets presented here (see Fig. 2 and the supplementary material). Also, to put the framework to the test, we apply the method to synthetic data, and we found that it was able to find the corresponding breakpoints and slopes (see the supplementary text). Nevertheless, when there are no reliable figures on negative test results, our framework fails to fit the data (see the supplementary material). More detailed data will be necessary to build better models in the future. Ideally, the data would be organized in a case-by-case fashion and it would contain information on the sample criterion. Fundamentally, this information should also be available also for negative tests.
Overall, we found evidence that masks reduce the spread of the SARS-CoV-2 and prevent new infections. We hope that our findings will persuade local authorities and intergovernmental institutions to strongly recommend the use of masks to prevent the spread of SARS-CoV-2. We arrived at this conclusion by merging two different traditions: causal inference and regularized regression. We believe that the union of these techniques will be fruitful in other contexts where the causes and effects are sparse in time.
peting interests: The authors declare no competing interests. Data and materials availability: all the data and codes to reproduce our analysis are publicly available at https: //github.com/ababino/corona.
Figs. S1 to S5 Tables S1 to S4
References (31-35) 
Figs. S1 to S5 Tables S1 to S4
We collected data from States that offer raw data on the number of tests and positive cases each day. We found that 16 States offer that information. Some of them offer an Abstraction Protocol Interface (API). Others serve a file with the information. Many states offer a visual dashboard with information on testing, but if they do not offer raw data, we did not use it. There is at least one project that collects data on testing from all the states: https://covidtracking. com/. This aggregator builds its database based on snapshots of the dashboards published by the states. The information in these snapshots is averaged three or four times a day. This process makes the accumulated number of tests and positives reliably, but not their daily change. That is why we did not use information from this aggregator, and we only use direct information from official sources. We provide links to each dataset in Table 1 In the main text, we show results for six datasets with the highest R 2 . To show the robustness and also limitations of the framework here, we show an analysis for all the states where we found data on the daily number of cases and tests. In the main text, we limit the analysis to the time in which NYS, MA, and CT ended their stay-at-home order. To show the robustness of the framework, here we show the results when we apply the method to a bigger timespan. Since, due to backlog, some states have a delay in reporting of about one week, we included data until the last day we find reliable data.
As can be seen in Fig. 5 the number of daily positive tests, Positive t , oscillates in synchrony with the number of tests. To overcome this source of noise, we propose to use the odds of a positive test:
where Negative t is the number of negative tests on day t.
We show the number of positives and the number of tests for each dataset in Fig. 5 We show the evolution of the Odds in Fig. 6 . The noise due to the variation in the number of tests is reduced, and a trend emerges.
As shown previously (28) , under the SIR model hypotheses, the number of newly infected individuals in a given day, k t , can be approximated as:
where R t is the instantaneous reproduction number (28) , and γ −1 is the average infectious period (26) estimated as 7.5 days (95% CI 5.3 to 195.3) according to (31) (in agreement with (32) , but higher than reported in (33) ).
Since we do not have access to the total number of infected individuals, but only to the tested population, we have to use some statistical assumptions about this population. If we assume that the people being tested is a random sample of the population with COVID-19-like symptoms, we can state that 1 :
where P t (I|symptoms) is the probability of a patient being positive for SARS-CoV-2 given that she is symptomatic, P t (symptoms) is the probability of having COVID-19-like symptoms, N is the total population, and f t is the fraction of people with symptoms that are selected to be tested (this number can be different each day, for example, if the number of tests available changes). Similarly:
where P t (notI|symptoms) is the probability of a patient being SARS-CoV-2 negative given he has COVID-19-like symptoms.
Now, if we assume that P t (symptoms|I) = cte, we can use Bayes theorem to show that: P t (I|symptoms)P t (symptoms) ∝ P t (I) = k t N Then:
Finally, if we assume that P t (notI|symptoms)P t (symptoms) = cte:
We used four sets of hypotheses. First, we use the assumptions of the SIR model. Second, we use that the tested population is a random sample from the population with COVID-19like symptoms ( Eqs. 1 and 2) . This assumption does not hold, for example, if the basis for testing someone is that she was in contact with a confirmed case. Third, we assume that P t (notI|symptoms)P t (symptoms) = cte. This hypothesis is equivalent to say that the number of people with COVID-19-like symptoms but without the SARS-CoV-2 (for example, people with the flu) is constant, or its change rate is negligible compared with the change rate in the number of symptomatic people with SARS-CoV-2. Fourth, we use that the symptoms show up instantaneously and that the tests are performed and processed on the same day (Eq. 3). This last hypothesis is not true, and it is the reason why, in our analysis, the effects of the interventions show a delay to onset between 8 and 11 days.
We write Eq. 4 as a linear function of the rate of change of R t . Defining
We can write Eq. 4 as:
Now, instead of using b t as the parameters to estimate we decompose each b t as follows:
The a i s represent the rate of change of the variable b t in logarithmic scale. Next, we replace the 7 in 6
Now, we linearize this result, and we generalize it to the case where t = 0 using the max function:
We can write 8 as a linear problem with the following definitions:
Importantly, the SIR hypotheses are only necessary to draw the connection to R t (Eq. 5).
However, Eq. 8 might hold even if the SIR hypotheses do not. What would change is the interpretation of the parameters.
Since in Eq. 9, we have as many regressors as samples, and we assume that the changes in a are only due to top-down interventions we use a LASSO regression to fit the data (27) . This regression minimizes the loss function:
This approach finds a sparse set of β i . We add two extra steps to sparsify even further this set of parameters. If there are contiguous β i = 0, we set to zero all of them but the first in the chunk.
Then, we fit the selected regressors using ordinary least squares, and we recursively remove the Table 2 , and the value and statistics of each parameter in Table 3 .
To compute the value of R t from fitter parameters, we have to use equations 12, 7 and 5. From these equations, we arrive at the following equality:
Where most of the β i values are zero. Using this formula, we arrive at the values presented in Fig. 3 and Table 1 (main text). We show the R t values for all states with data in Fig. 7 .
To test the framework, we show that it is able to estimate the parameters from simulated data.
We use the SIR equations:
where S is the number of susceptible individuals, I the number of infected individuals, R the number of recovered (or death) individuals, β is the contact rate and γ the inverse of the average infectious period. Also, R 0 = β/γ
In our simulation we change the value of β at times t = 21, t = 41, and t = 61, from a initial value of β = 0.26 to β = 0.2, and then to β = 0.093 and β = 0.053. Also we set γ = 7.5 −1
We chose these values because the R t are similar to the one in the NYS dataset. After carrying out the simulation we compute the log-odds= log(I/S) and add noise N (0, 0.05) We show the results applying the framework to this data in Fig. 9 and the real and estimated values in Table 4 . Initially, the framework estimates R t = 1.7 (1.95, 2.0) 95% CI. This value is below the R 0 = 2 of the model, which is reasonable given that the framework should estimate the mean value of R t in a segment, which is always below R 0 (26). The second discontinuity is found at t = 39 and last one at t = 59, two days before the real change.
This analysis shows that the framework can detect the number of discontinuities correctly, and it also estimates the value of R t accurately and the times of the breaks with no more than 2 With the selected coefficients, we carry out a linear fit, and we iteratively remove the coefficients with p-value * < 0.01. The dots indicate the parameters that were selected, either for being the first of a chunk or being significant. Third row: BIC as a function of α. Since we performed these steps for all the α values, we compute the BIC for each of them. Then, we pick the one with the lowest BIC value. 