As of today, there are almost 14 million of confirmed cases of coronavirus disease 2019 in the world, and a total of 591,228 reported deaths (Worldometer) . Mexico has 317,635 COVID-19 confirmed cases with 36,906 deaths. The numbers are growing rapidly in the country, along with plans for a gradual reopening which started in June. In addition to this, the testing rate in Mexico is the lowest among the OECD countries, implying a likely large sub-reporting of cases. The Ministry of Health of Mexico has declared that the data available are inconsistent and, therefore, difficult to use to define the reopening of activities. Given this, a better understanding of the current epidemiological situation at the state level, and the impact of government-imposed mitigation measures is necessary.
Here, we use the semi-mechanistic Bayesian hierarchical model of COVID-19 epidemiological dynamics found in Flaxman et al. (2020) to assess the effect of different control interventions in Mexico. We estimate the number of infections and the time-varying reproduction number (the expected number of secondary cases caused by each infected individual) as a function of human mobility. Figure  1 shows the cases and deaths reported in Mexico. The distribution of deaths among states is highly heterogeneous, with Mexico City and State of Mexico accounting for almost 40% of the deaths reported to date.
Mexico's government determined three phases of contingency for the COVID-19 pandemic. The first two cases were confirmed on February 27, from travellers returning from Italy. Later, on March 17, the suspension of school activities began nationally. A week later, on March 23, The Comité Nacional para la Seguridad en Salud announced the start of the second phase of the epidemic. Almost a month later, on April 21, phase three was declared. Distancing measures included suspension of all non essential activities of public, private and social sectors. The measures were initially planned to last until April 30, but were later extended depending on the local situation of every municipality of the country.
In June, the government announced a monitoring system to regulate the use of public space in accordance with the risk of COVID-19 contagion, with the aim of moving towards a new normality. They defined the traffic light of state epidemiological risk. This traffic light is made up of four colors: red, orange, yellow and green. Red is the color with the highest epidemiological risk, where only essential economic activities are allowed. The traffic light moves from color to color until reaching green, when all activities are allowed. In the last update, 17 entities were in orange, with high risk, and 15 in red, with maximum risk. It remains unclear the extent to which this traffic light has been effective. This week, the government did not present an epidemiological traffic light because state data is not consistent and the reopening has caused outbreaks. Given this, a better understanding of the current situation would be highly beneficial to guide policy decisions.
The outline of the paper is as follows. In Section 2 we explain briefly the model and some important modelling assumptions. In Section 3 we describe our data sources. In Section 4 we show the application of the model results. Finally, in Section 5 we give some concluding remarks and future research directions.
The main idea of Flaxman et al. (2020) is to model the expected number of deaths in a given country on a given day as a function of the number of infections occurring in previous days, linking the deaths to the infected cases by using two quantities. The first is the infection fatality ratio (IFR), which is the probability of death given infection. The second is the distribution of times from infection to death. The number of deaths every day is the sum of the past infections weighted by their probability of death, where the probability of death depends on the number of days since infection and the state IFR.
For the state-level IFR, we made adjustments accounting for the substantial heterogeneity we expect to observe in health outcomes across states, due to variation in healthcare quality and capacity (as in Mellan et al., 2020) . We consider the IFR of Mexico City as the IFR estimated from the Chinese epidemic Verity et al., 2020) , and the IFR of Chiapas and Oaxaca, which are the most marginalized states in Mexico, more similar to those that would be expected in a Lower Middle Income Country (see Walker et al., 2020) . Then, we weight according to age groups. Finally, we interpolated the rest of the states using five levels of the State Marginalization Index. The marginalization index is an official standard indicator that allows states to be differentiated according to the deprivations suffered by the population, such as education, housing, monetary income, or access to health services. We also allow the IFR for every state to have additional noise, multiplying it by a N (1, 0.1) distribution.
The distribution of times from infection to death is the sum of two independent random times: the incubation period and the time between onset of symptoms and death. As in Flaxman et al. (2020) , we assume the infection-toonset distribution is a Gamma(5.1, 0.86), and the onset-to-death distribution is a Gamma(18.8, 0.45).
Now, in Mexico, there are two evident and significant problems that must be considered with death data. One of them is the death-delayed reporting and the other is the under-reporting. The death-delayed reporting refers to COVID-19 deaths that will take time to report (see Castañeda and Garrido) . Death records must be validated by two information systems, which obtain information from medical units and health jurisdictions. Many of the COVID-19 deaths that were tested are awaiting the proper certification and classification. The assessment and confirmation processes take days. That is why the official figures of deaths from COVID-19, published by the Ministry of Health, have an important lag. Figure 2 illustrates the cumulative proportion of deaths reported as the delay increases, with different colours corresponding to different dates of death. Very few deaths (around 10%) are reported within one day, and the vast majority of deaths (around 90%) are reported after a month. To account for this, we are adjusting the deaths reported in the following way. We assume the delay has the same distribution in all the country, but this can be corrected when more data becomes available. Let p t be the proportion of cases reported on day t, then
where the parameters η 1 , ..., η k are the empirical proportion of cases reported, and α ∼ Gamma(100, 1). Notice that, since η 1 + ...
Afterwards, we are multiplying the case numbers of each day by the estimated proportion of cases reported. We do not report the last two days, since underreporting in those two days is very high. Also, to avoid biasing the model, we only include observed deaths from the day after a state has accumulated 10 deaths.
The under-reporting refers to COVID-19 deaths that will not be reported (or it will take more than a year). Estimating the extent of the under-reporting remains very difficult, especially in areas outside of Mexico City where less studies have been conducted. Estimates of excess deaths can provide information about the burden of mortality potentially related to the COVID-19 pandemic, including deaths that are directly or indirectly attributed to COVID-19. However, in Mexico, there are no official published data on deaths during 2020. The publication of death statistics is annual and the last record is from 2018.
In Mexico City, there are several death certificates stating that the confirmed or probable cause of death was COVID-19, three times higher than the COVID-19 deaths reported by the government (see Adam et al.) . Another study of the death certificates shows that the excess mortality derived from the COVID-19 health crisis in Mexico City is four times greater than that reported (Romero and Despeghel) . This excess includes both people who died from COVID-19 and also those who died from other causes derived from the health crisis. According to the Registro Nacional de Población (Renapo), the number of deaths from COVID-19 in the country doubles what is reported every day (Navarro and López) . However, the government rejects that these death figures are official.
We adopted the extension of Mellan et al. (2020) to reflect the uncertainty about under-reported deaths. Meaning that we address the effect of underreporting in the data by setting a ψ m ∼ Beta(θ, ρ) prior distribution to death under-reporting in each state m, where both of the hyperparameters of the beta density are fixed in 80, in order to get a mode of 50% (meaning that, on average, half of the COVID-19 deaths are being reported). Figure 3 shows the plot of this prior distribution.
With all of this considerations, if we let D t,m be the daily deaths attributed to COVID-19 at day t and state m, we model D t,m using a positive real-valued function d t,m = E[ψ m p t D t,m ] that represents the expected number of deaths, taking into account the designated state-under-reported rate ψ m and the deathdelayed reporting adjustment p t .
The rest of the mathematical model follows the original manuscript of Flaxman et al. (2020) . Daily deaths are assumed to follow a negative binomial distribution with mean d t,m and variance d t,m + dt,m φ , where φ follows a positive half normal distribution, φ ∼ N + (0, 5), and d t,m is given by the following discrete sum:
IFR m is the IFR of the state m, c τ,m is the number of new infections on day τ in state m, and π is the discretized distribution of times from infection to death.
To model the number of infections over time, we need to specify a serial interval distribution g with density g(τ ) (the time between when a person gets infected and when they subsequently infect another other people), which we set as in Flaxman et al. (2020) , g ∼ Gamma(6.5, 0.62). With this, the number of infections c t,m on a given day t, and a state m, is then given by:
m is a population adjustment to account for susceptible individuals, g m is the discretized serial interval distribution, and R t,m is the state-specific time-varying reproduction number which was parametrized as a function of Google mobility data as explained below.
Denote by X k,t,m to the k-th Google mobility indicator at time t for state m, smoothed by a seven day moving average (meaning that each point is the mean of the day and the previous 6 days). To allow flexibility, we define dummy variables after June, were relaxations started. Denote this variables by Z k,t,m , then Z k,t,m = 1 if the day t is after June 1 or Z k,t,m = 0 otherwise, for k = 1, ...4 and for any state m. The state-specific time-varying reproduction number R t,m is modelled by:
The basic reproduction number, R 0,m , has a normal prior for every state m, R 0,m ∼ N (3.28, k), where k ∼ N + (0, 0.5). This premise is tested in Mellan et al. (2020) through a sensitivity analysis. The conclusion is that, for this type of model, the results vary slightly depending on it. However, our major concern must be two highly influential assumptions, each state IFR and the under-reporting distribution.
As input of deaths and reported cases, our model uses Open data from the Dirección General de Epidemiología, and the Historical bases on cases associated with COVID-19 nationwide. For population counts we are using data provided by Víctor M. García as the Consejo Nacional de Población population projections for 2020 are not disaggregated for the age groups over 65 years old. Regarding intervention data, the values taken into account are the dates in which interventions were effectively applied, even though they were encouraged at earlier dates.
The attack rate (AR) is the proportion of individuals that has been infected to date. The country's AR is estimated in 18% (95% Confidence Interval (CI): 16%-19%) on July 7. Of this infections, 7,235,600 (95% CI: 6,281,900-8,062,000) have occurred in the previous 14 days (which is the official definition of active cases in Mexico). Table 1 shows the estimated IFR, state population, reported deaths and deaths per million population, estimated number of infections in thousands, estimated number of infections in the previous 14 days in thousands, and estimated AR for all the states. Our computed IFRs for each of the states range from 0.4% to 1.1%, reflecting substantial differences between states in their demographic structure and healthcare provision. The percentage of people that have been infected with SARS-CoV-2 ranges from 47% (95% Cl: 44%-50%) in Baja California to 5% (95% Cl: 4%-6%) in Michoacán. Figure 4 shows the estimates of infections, deaths and reproduction number for each state. The reproduction number is an important parameter for assessing whether additional interventions are required. At the start of the epidemic, the basic reproduction number meant that an infected individual would infect three others on average. If the reproduction number is above one, a sustained outbreak is likely. The aim of control interventions is typically to reduce the reproduction number below one. Following non-pharmaceutical interventions such as school closures and decreases in population mobility, our results show substantial reductions in the estimated value of R t,m in each state m. However, we can notice growth from June to this date, which has translated to a linear growth of the infections. In some states, the number of infections goes down and then comes back up again, making what many call a second wave.
We applied the semi-mechanistic Bayesian hierarchical model found in Flaxman et al. (2020) as an alternative model to estimate the number of infections and the reproduction number in Mexico. In this context, we have incomplete and biased data. We expect the model to provide more accurate estimates, since the estimation is made with deaths observed over time, on the basis that the death data is more reliable than the data of infected people. Weekly updates can be checked at https://covid19mex.netlify.app/.
The results presented here suggest an ongoing epidemic in which substantial reductions in the average reproduction number have been achieved through nonpharmaceutical interventions. However, our results also show that so far the changes in mobility have not been enough to reduce the reproduction number below one. Therefore we predict continued growth of the epidemic across Mexico and increases in the associated number of cases and deaths unless further actions are taken.
Our results also reveal extensive heterogeneity in predicted attack rates between states, suggesting that the epidemic is at a far more advanced stage in some states compared to others. Despite this, in most states, a small proportion of individuals has been infected to date, indicating that herd immunity is not close yet. The exception is perhaps Baja California, where almost half of the people have already been infected.
It is now established that SARS-CoV-2 pre-existing immune reactivity exists to some degree in the general population (see, e.g., Sette and Crotty, 2020) . Pre-existing cross-reactivity against COVID-19 in a fraction of the human population is relevant in our modelling, that considers at every time a population adjustment for susceptible individuals. This can be introduced naturally in the model and is worthy of further exploration when there is more information. Also, there is still a lot to be done regarding the adjustment in death-delayed reporting, we will be able to set a distribution to death-delay for each state when more data is available.
To conclude, let us highlight that this model uses mobility to predict the rate of transmission, neglecting the potential effect of additional behavioural changes or interventions such as increased mask wearing, changes in age specific movement, testing and tracing. Figure 4 : For each state, the first plot has the daily number of infections, brown bars are reported cases, blue bands are predicted infections, dark blue 50% CI, and light blue 95% CI. The second plot shows daily number of deaths, brown bars are reported deaths, pink bars are estimated deaths adjusted for the reporting lag and the under-reporting, blue bands are predicted deaths, and CI as in the first plot. The third plot is the time-varying reproduction number, dark green 50% CI, light green 95% CI. The dotted line with a square is the date of school closures. The dotted line with the circle marks the start of sanitary emergency. Lastly, the line with the rhombus is the start of phase 3.