Transmission dynamics and forecasts of the COVID-19 pandemic in Mexico, March 20-November 11, 2020.

The ongoing coronavirus pandemic reached Mexico in late February 2020. Since then Mexico has observed a sustained elevation in the number of COVID-19 deaths. Mexicos delayed response to the COVID-19 pandemic until late March 2020 hastened the spread of the virus in the following months. However, the government followed a phased reopening of the country in June 2020 despite sustained virus transmission. In order to analyze the dynamics of the COVID-19 pandemic in Mexico, we systematically generate and compare the 30-day ahead forecasts of national mortality trends using various growth models in near real-time and compare forecasting performance with those derived using the COVID-19 model developed by the Institute for Health Metrics and Evaluation. We also estimate and compare reproduction numbers for SARS-CoV-2 based on methods that rely on both the genomic data as well as case incidence data to gauge the transmission potential of the virus. Moreover, we perform a spatial analysis of the COVID-19 epidemic in Mexico by analyzing the shapes of COVID-19 growth rate curves at the state level, using techniques from functional data analysis. The early estimate of reproduction number indicates sustained disease transmission in the country with R~1.3. However, the estimate of R as of September 27, 2020 is ~0.91 indicating a slowing down of the epidemic. The spatial analysis divides the Mexican states into four groups or clusters based on the growth rate curves, each with its distinct epidemic trajectory. Moreover, the sequential forecasts from the GLM and Richards model also indicate a sustained downward trend in the number of deaths for Mexico and Mexico City compared to the sub-epidemic and IHME model that outperformed the others and predict a more stable trajectory of COVID-19 deaths for the last three forecast periods.


Model calibration and forecasting approach 236
We conducted 30-day ahead short-term forecasts utilizing thirteen data sets spanned over a 237 period of four months (July 4-October 9, 2020) ( Table 1). Each forecast was fitted to the daily 238 death counts from the IHME smoothed data estimates reported between March 20-September 27, 239 2020 for (i) Mexico and (ii) Mexico City. The first calibration period relies on data from March 240 20-July 4, 2020. Sequentially models are recalibrated each week with the most up-to-date data, 241 meaning the length of the calibration period increases by one week up to August 2, 2020. 242 However, owing to irregular publishing of data estimates by the IHME, after August 2, 2020 the 243 length of calibration period increased by 2 weeks, followed by a one week increase from August 244 17-September 27, 2020 as the data estimates were again published every week. 245 . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

(which was not certified by peer review)
The copyright holder for this preprint this version posted January 13, 2021. We compare the cumulative mortality forecasting results obtained from our models to the 251 cumulative smoothed death estimates reported by IHME reference scenario as of November 11, 252 2020. Weekly forecasts obtained from our models were compared to the forecasts obtained from 253 the three IHME modeling scenarios for the same calibration and forecasting periods.  corresponds to the set of parameters of the GLM model. For the GLM and sub-epidemic wave 263 model, we provide initial best guesses of the parameter estimates. However, for the Richards 264 growth model, we initialize the parameters for the nonlinear least squares' method [53] over a 265 wide range of feasible parameters from a uniform distribution using Latin hypercube sampling 266 . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

(which was not certified by peer review)
The copyright holder for this preprint this version posted January 13, 2021. ;https://doi.org/10.1101https://doi.org/10. /2021 [54] in order to test the uniqueness of the best fit model. Moreover, the initial conditions are set 267 at the first data point for each of the three models [55]. Uncertainty bounds around the best-fit 268 solution are generated using parametric bootstrap approach with replacement of data assuming a 269 Poisson error structure for the GLM and sub-epidemic model, whereas a negative binomial error 270 structure was used to generate the uncertainty bounds of the Richards growth model; where we 271 assume the mean to be three times the variance based on the noise in the data. Detailed 272 description of this method is provided in ref [55]. 273 274 Each of the M best-fit parameter sets are used to construct the 95% confidence intervals for each 275 parameter by refitting the models to each of the M = 300 datasets generated by the bootstrap 276 approach during the calibration phase. Further, each M best fit model solution is used to generate 277 m= 30 additional simulations with Poisson error structure for GLM and sub-epidemic model and 278 negative binomial error structure for Richards model extended through a 30-day forecasting 279 period. We construct the 95% prediction intervals with these 9000 (M × m) curves for the 280 forecasting period. Detailed description of the methods of parameter estimation can be found in 281 references [55-57] 282 283

Performance metrics 284
We utilized the following four performance metrics to assess the quality of our model fit and the 285 30-day ahead short term forecasts: the mean absolute error (MAE) [58], the mean squared error 286 (MSE) [59], the coverage of the 95% prediction intervals [59], and the mean interval score (MIS) 287 [59] for each of the three models: GLM, Richards model and the sub-epidemic model. For 288 calibration performance, we compare the model fit to the observed smoothed death data 289 . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

(which was not certified by peer review)
The copyright holder for this preprint this version posted January 13, 2021. ;https://doi.org/10.1101https://doi.org/10. /2021 estimates fitted to the model, whereas for the performance of forecasts, we compare our forecasts 290 with the smoothed death data reported on November 11, 2020 for the time period of the forecast. 291 The mean squared error (MSE) and the mean absolute error (MAE) assess the average deviations 292 of the model fit to the observed death data. The mean absolute error (MAE) is given by: 293 The mean squared error (MSE) is given by: 294 is the time series of reported smoothed death estimates, ‫ݐ‬ is the time stamp and Θ is the 296 set of model parameters. For the calibration period, n equals the number of data points used for 297 calibration, and for the forecasting period, n = 30 for the 30-day ahead short-term forecast. 298 299 Moreover, in order to assess the model uncertainty and performance of prediction interval, we 300 use the 95% PI and MIS. The prediction coverage is defined as the proportion of observations 301 that fall within 95% prediction interval and is calculated as: 302 where Y t are the smoothed death data estimates, L t and U t are the lower and upper bounds of the 303 95% prediction intervals, respectively, n is the length of the period, and I is an indicator variable 304 that equals 1 if value of Y t is in the specified interval and 0 otherwise 305 . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) The copyright holder for this preprint this version posted January 13, 2021. ;https://doi.org/10.1101https://doi.org/10. /2021

306
The mean interval score addresses the width of the prediction interval as well as the coverage. 307 The mean interval score (MIS) is given by: 308 In order to analyze the time series data for Mexico from March 20-December 5, 2020 for three 317 modes of mobility; driving, walking and public transport we utilize the R code developed by 318 Kieran Healy [60]. We analyze the mobility trends to look for any pattern in sync with the curve 319 of COVID-19 death counts. The time series for mobility requests is decomposed into trends, 320 weekly and remainder components. The trend is a locally weighted regression fitted to the data 321 and remainder is any residual left over on any given day after the underlying trend and normal 322 daily fluctuations have been accounted for. 323 324

Reproduction number 325
. CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) The copyright holder for this preprint this version posted January 13, 2021. is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) The copyright holder for this preprint this version posted January 13, 2021. ;https://doi.org/10.1101https://doi.org/10. /2021 interval of SARS-CoV-2 is modeled assuming gamma distribution with a mean of 5.2 days and a 348

367
. CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

(which was not certified by peer review)
The copyright holder for this preprint this version posted January 13, 2021.  and we assume that the 380 generation interval equals the SI. The SI was assumed to follow a gamma distribution with a 381 mean of 5.2 days and a standard deviation of 1.72 days [65]. Analytical estimates of R t were 382 obtained within a Bayesian framework using EpiEstim R package in R language [73]. R t was 383 estimated at 7-day intervals. We reported the median and 95% credible interval (CrI). 384 385 Estimating the reproduction number, R, from the genomic analysis. 386 In order to estimate the reproduction number for the SARS-CoV-2 between February 27-May 387 29, 2020, from the genomic data, 111 SARS-CoV-2 genomes sampled from infected patients 388 from Mexico and their sampling times were obtained from GISAID repository [46]. Short 389 sequences and sequences with significant number of gaps and non-identified nucleotides were 390 . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

(which was not certified by peer review)
The copyright holder for this preprint this version posted January 13, 2021. . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

(which was not certified by peer review)
The copyright holder for this preprint this version posted January 13, 2021.
To identify the clusters by comparing the curves, we used a simple metric. For any two 417 rate curves, h i and h j , we compute the norm ||h i −h j ||, where the double bars denote the L 2 418 This process is depicted in S17 Fig To perform clustering of 32 curves into smaller groups, we apply the 423 dendrogram function in Matlab using the "ward" linkage as done in the previous work [80]. The 424 number of clusters is decided empirically based on the display of overall clustering results. After 425 clustering the states into different groups, we derived average curve for each cluster after using a 426 time wrapping algorithm as done in previous work [80,81]. 427 428 Twitter data analysis. 429 To observe any relationship between the COVID-19 cases by dates of symptom onset and the 430 frequency of tweets indicating stay-at-home orders we used a public dataset of 698 Million 431 tweets of COVID-19 chatter [47]. The frequency of tweets indicating stay-at-home order is used 432 to gauge the compliance of people with the orders of staying at home to avoid spread of the virus 433 by maintaining social distance. Tweets indicate the magnitude of the people being pro-lockdown 434 and show how these numbers have dwindled over the course of the pandemic. To get to the 435 . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

(which was not certified by peer review)
The copyright holder for this preprint this version posted January 13, 2021. ; https://doi.org/10.1101/2021.01.11.21249561 doi: medRxiv preprint plotted data, we removed all retweets and tweets not in the Spanish language. We also filtered by 436 the following hashtags: #quedateencasa, and #trabajardesdecasa, which are two of the most used 437 hashtags when users refer to the COVID-19 pandemic and their engagement with health 438 measures. Lastly, we limited the tweets to the ones that originated from Mexico, via its 2-code 439 country code: MX. A set of 521,359 unique tweets were gathered from March 12 to November 440 11, 2020. We then overlay the curve of tweets over the epi-curve in Mexico to observe any 441 relation between the shape of the epidemic case trajectory and the frequency of tweets during the 442 time period, March 12-November 11, 2020. We also estimate the correlation coefficient between 443 the cases and frequency of tweets. The driving and walking trend subsequently increased in June with the reopening of the non-455 essential services. Fig 1 (upper panel) shows that reopening coincides with the highest levels of 456 daily deaths. These remain at a high level for just over two months (June and July). Then, from 457 . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

(which was not certified by peer review)
The copyright holder for this preprint this version posted January 13, 2021. ; https://doi.org/10.1101/2021.01.11.21249561 doi: medRxiv preprint mid-August, the number of deaths begins to fall, reaching a reduction of nearly 50% by mid-458 October. But at the end of October a new growth begins. In the subsequent sections, we first present the results for the short-term forecasting, followed by 468 the estimation of the reproduction numbers from the early phase of the COVID-19 epidemic in 469 Mexico and the instantaneous reproduction number. Then we present the results of the spatial 470 analysis followed by the results of the twitter data analysis. 471 472

Model calibration and forecasting 473
Here we compare the calibration and 30-day ahead forecasting performance of our three models: 474 the GLM, Richards growth model and the sub-epidemic wave model between March 20-475 September 27, 2020 for (i) Mexico and (ii) Mexico City. We also compare the results of our 476 cumulative forecasts with the IHME forecasting estimates (for the three scenarios) obtained from 477 the IHME data. 478 479 . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

(which was not certified by peer review)
The copyright holder for this preprint this version posted January 13, 2021. ; https://doi.org/10.1101/2021.01.11.21249561 doi: medRxiv preprint Calibration performance. Across the thirteen sequential model calibration phases for Mexico 480 over a period of seven months (March-September), based on the calibration phases as provided 481 in Table S1 and (Table S1) CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

(which was not certified by peer review)
The copyright holder for this preprint this version posted January 13, 2021. ; https://doi.org/10.1101/2021.01.11.21249561 doi: medRxiv preprint models in four calibration phases (3/20-7/04, 3/20-7/11,3/20-7/17, 3/20-8/02). The Richards 503 model has much higher estimates for the MIS compared to the other two models. The 95% PI 504 across all thirteen calibration phases lies between 91.6-99.6% for the sub-epidemic model, 505 followed by the Richards model (85.9-100%) and the GLM model (53.2-100%) ( Table S2, Table  518 S4. For Mexico, the sub-epidemic model consistently outperforms the GLM and Richards 519 growth model for ten out of the thirteen forecasting phases in terms of RMSE and MAE, eight 520 forecasting phases in terms of MIS and nine forecasting phases in terms of the 95% PI coverage. 521 This is followed by the GLM and then the Richards growth model. CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

(which was not certified by peer review)
The copyright holder for this preprint this version posted January 13, 2021.

Comparison of daily death forecasts 540
The thirteen sequentially generated daily death forecasts from GLM and Richards growth model, 541 for Mexico and Mexico City indicate towards a sustained decline in the number of deaths (S1 542 Fig, S2 Fig, S3 Fig and S4 Fig). However, the IHME model (actual smoothed death data 543 estimates) shows a decline in the number of deaths for the first six forecasts periods followed by 544 a stable epidemic trajectory for the last seven forecasts, for Mexico City and Mexico. Unlike the 545 GLM and Richards models, the sub-epidemic model is able to reproduce the observed 546 stabilization of daily deaths observed after the first six forecast periods for Mexico and the last 547 three forecasts for Mexico City (S5 Fig, S6 Fig, S7 Fig and S8 Fig)  548 . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

Comparison of cumulative mortality forecasts 550
The total number of COVID-19 deaths is an important quantity to measure the progression of an 551 epidemic. Here we present and compare the results of our 30-day ahead cumulative forecasts 552 generated using GLM, Richards and sub-epidemic growth model against the actual reported 553 smoothed death data estimates using the estimates of deaths obtained from the IHME death data 554 reported as of November 11, 2020. Figs 6 and 7 present a comparative assessment of the 555 estimated cumulative death counts derived from our three models along with their comparison 556 with three IHME modeling scenarios; current projection, universal masks and mandates easing; 557 for our thirteen sequentially generated forecasts. projections (IHME C.P), IHME universal masks (IHME U.M) and IHME mandates easing 568 (IHME M.E) to predict the cumulative COVID-19 deaths for the Mexico City in the thirteen 569 sequential forecasts. The blue circles represent the mean deaths and the magenta vertical line 570 . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

(which was not certified by peer review)
The copyright holder for this preprint this version posted January 13, 2021. ; https://doi.org/10.1101/2021.01.11.21249561 doi: medRxiv preprint indicates the 95% PI around the mean death count. The horizontal dashed line represents the 571 actual death count reported by that date in the November 11, 2020 IHME estimates files. 572 573 Mexico. The 30 day ahead forecast results for the thirteen sequentially generated forecasts for 574 Mexico utilizing GLM, Richards model, sub-epidemic growth model and IHME model are 575 presented in S9 Fig, S10 Fig, S11 Fig and S12 Fig. The  CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

(which was not certified by peer review)
The copyright holder for this preprint this version posted January 13, 2021. is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

(which was not certified by peer review)
The copyright holder for this preprint this version posted January 13, 2021. ; https://doi.org/10. 1101/2021 In summary, the Richards growth model consistently under-estimates the actual death count 593 compared to the GLM, sub-epidemic and three IHME model scenarios. The GLM model also 594 provides lower estimates of the mean death counts compared to the sub-epidemic and three 595 IHME model scenarios, but higher mean death estimates than the Richards model. The 95% PI 596 for the Richards model is substantially wider than the other two models. Moreover, the three 597 IHME model scenarios predict approximately similar cumulative death counts across the thirteen 598 generated forecasts, hence indicating that they do not differ substantially. 599 600 Mexico City. The 30 day ahead forecast results for thirteen sequentially generated forecasts for 601 the Mexico City utilizing GLM, Richards model, sub-epidemic growth model and IHME model 602 are presented in S13 Fig, S14 Fig, S15 Fig and S16 Fig. The cumulative death comparison is 603 given in Fig 7 and Table 3. For the fourth and fifth generated forecast all models under-predict 604 the true death counts (11,326, 11,769 deaths respectively). For the first and second generated 605 forecast the sub-epidemic model and the IHME scenarios closely approximate the actual death 606 count (~10,081, ~10,496 deaths). For the third and sixth generated forecast GLM and Richards 607 model underestimate the actual death count (~10,859, ~12,615 deaths respectively) whereas the 608 sub-epidemic model over predicts the death count for the third forecast and under-predicts the 609 death count for the sixth forecast. The three IHME model scenarios seem to predict the actual 610 death counts closely. From the seventh-thirteenth generated forecasts, all models under-predict 611 the actual death counts. 612 613 Table 3: Cumulative mortality estimates for each forecasting period of the COVID-19 pandemic 614 in Mexico City (2020). 615 . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

(which was not certified by peer review)
The copyright holder for this preprint this version posted January 13, 2021. . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

(which was not certified by peer review)
The copyright holder for this preprint this version posted January 13, 2021.  Estimate of reproduction number, from case incidence data. The reproduction number 638 from the case incidence data (February 27-May 29, 2020) using GGM was estimated at 639 . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

(which was not certified by peer review)
The copyright holder for this preprint this version posted January 13, 2021. ; https://doi.org/10.1101/2021.01.11.21249561 doi: medRxiv preprint ܴ ௧ 1.1(95% CI:1.1,1.1), in accordance with the estimate of ܴ ௧ obtained from the genomic data 640 analysis. The growth rate parameter, r, was estimated at 1.2 (95% CI: 1.1, 1.4) and the 641 deceleration of growth parameter, p, was estimated at 0.7 (95% CI: 0.68,0.71), indicating early 642 sub-exponential growth dynamics of the epidemic (Fig 9). reproduction number for Mexico and the red shaded area represents the 95% credible interval 662 . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

(which was not certified by peer review)
The copyright holder for this preprint this version posted January 13, 2021. . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

(which was not certified by peer review)
The copyright holder for this preprint this version posted January 13, 2021. ; https://doi.org/10.1101/2021.01.11.21249561 doi: medRxiv preprint clearly visible. For cluster 1, the rate rises rapidly from April to July and then shows small 686 fluctuations. For cluster 2, there is rapid increase in growth rate from April to July followed by a 687 rapid decline. Chihuahua in cluster 3 shows a slow growth rate until September followed by a 688 rapid rise until mid-September which then declines rapidly. For cluster 4, the rate rises slowly, 689 from April to September, and then shows a rapid rise (Fig S20). 690 691 From the colormap (Fig 12) we can see that the cases were concentrated from the beginning in CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

(which was not certified by peer review)
The copyright holder for this preprint this version posted January 13, 2021. ; https://doi.org/10.1101/2021.01.11.21249561 doi: medRxiv preprint hashtag (stay-at-home order hashtag) has been gradually declining as the number of cases have 709 continued to increase or remain at a steady pace, showing the frustration of public on the 710 relaxation of lock downs. Mostly the non-government public health experts are calling for more 711 lockdowns or social distancing measures but are not being heard by the authorities. It could also 712 imply that the population is not following the government's stay at home orders and hence we 713 continue to observe the cases. S22 Fig shows that  The results of our GLM model fit to the death data for all the thirteen calibration phases and 721 GGM fit to the case incidence data indicate sub-exponential growth dynamics of COVID-19 722 epidemic in the Mexico and Mexico City with the parameter p estimated between (p~0.6-0.8). 723 Moreover, the early estimates of R indicate towards sustained disease transmission in the country 724 with ܴ ௧ estimated at 0.9 as of September 27, 2020. As the virus transmission continues in 725 Mexico, the twitter analysis implies the relaxation of lockdowns, with not much decline in the 726 mobility patterns observed over the last few weeks as evident from the Apple's mobility trends. 727 Moreover, the systematic comparison of our models across thirteen sequential forecasts deem 728 sub-epidemic model as the most appropriate model for mortality forecasting. The sub-epidemic 729 model is also able to reproduce the stabilization in the trajectory of mortality forecasts as 730 observed by the most recent IHME forecasts. 731 . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

732
The sub-exponential growth pattern of the COVID-19 epidemic in Mexico can be attributed to a 733 myriad of factors including non-homogenous mixing, spatial structure, population mobility, 734 behavior changes and interventions [83]. These results are also consistent with the sub-735 exponential growth patterns of COVID-19 outbreaks observed in Mexico [84] and Chile [85]. 736 Along with the sub-exponential growth dynamics, the reproduction number estimated from the 737 genomic sequence analysis and the case incidence data (ܴ In general, Mexico has seen a sustained SARS-CoV-2 transmission and an increase or a 750 sustained steady number of cases despite the social distancing interventions including the stay-at-751 home orders that were declared on April 12, 2020 and eased around June 2020. As our twitter 752 data analysis also shows, the number of cases by onset dates was negatively correlated to the 753 stay-at-home orders, indicating population's frustration towards the stay-at-home orders and 754 . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

(which was not certified by peer review)
The copyright holder for this preprint this version posted January 13, 2021. ;https://doi.org/10.1101https://doi.org/10. /2021 reopening the country resulting in pandemic fatigue. Hence, people might have stopped 755 following the government's preventive orders to stay at home. Mexico has been one of the 756 countries where the stay-at-home orders have been least respected, with an average reduction in 757 mobility in Mexico reported to be ~35. 4% compared to 71% in Brazil,and 86% in Argentina and 758 Colombia [90]. The intensity of these orders have affected the population disproportionately, 759 with some proportion of the population showing aggression towards quarantine and stay-at-home 760 orders [32]. We can also appreciate the variable spatial-temporal dynamics of the COVID-19 761 epidemic in Mexico. Our classification of epidemic pattern at the state level in Mexico shows 762 distinct variation of growth rates across states. For instance, cluster 1 including Baja California, 763 Colima and Mexico City has stable growth at a higher rate and cluster 4 including 764 Aguascalientes, Durango, Queretaro, and Zacatecas shows a rising pattern in the growth rate 765 Richards) are less reactive. This model can also be utilized as a potential forecasting tool for 776 other cities in Mexico and compared with other models. Moreover, shorter forecasts (5,10 days) 777 . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

(which was not certified by peer review)
The copyright holder for this preprint this version posted January 13, 2021. ;https://doi.org/10.1101https://doi.org/10. /2021 could be also be conducted with the sub-epidemic model using the consecutive calibration 778 phases to further reduce the error metrics [48]. 779 780 Overall, the sequential forecasts based on the daily smoothed death estimates for Mexico from 781 the two models (GLM, Richards) suggest a decline in over-all deaths (S1 Fig and S2 Fig)  Similarly, for Mexico City, the sequential forecasts obtained from the GLM and Richards model 799 fit to the daily death data estimates indicate a decline in the overall deaths (S3 Fig, S4 Fig). The 800 . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

(which was not certified by peer review)
The copyright holder for this preprint this version posted January 13, 2021. ;https://doi.org/10.1101https://doi.org/10. /2021 IHME and sub-epidemic models on the other hand indicate a stabilized trajectory of deaths for 801 the last three forecast periods (S7 Fig and S8 Fig)  the three IHME modeling scenarios. Whereas the GLM slightly under predicts the cumulative 813 death counts (Fig 6, Fig 7) The three models (GLM, Richards, sub-epidemic model) used in this study generally provide 819 good fits to the mortality curve, based on the residuals, with the Richards model unable to 820 capture the early sub-exponential dynamics of the death curve. Moreover, these 821 phenomenological models are particularly valuable for providing rapid predictions of the 822 epidemic trajectory in complex scenarios that can be used for real-time preparedness since these 823 . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

(which was not certified by peer review)
The copyright holder for this preprint this version posted January 13, 2021. ; https://doi.org/10.1101/2021.01.11.21249561 doi: medRxiv preprint models do not require specific disease transmission processes to account for the interventions. 824 Since these models do not explicitly account for behavioral changes, the results should be 825 interpreted with caution. Importantly, since the death curves employed in this study are reported 826 according to the date of reporting, they are likely influenced by variation in the testing rates and 827 related factors including the case fatality rates. Further, delays in reporting of deaths due to the 828 magnitude of the epidemic could also influence our predictions. Moreover, using the reporting 829 date is not ideal, due to the time difference between the death date and the reporting date of 830 death, which at a given moment can give a false impression of the ongoing circumstances. 831 832 This paper is not exempt from limitations. First, the IHME (current projections, mandated mask, 833 and worst-case scenario) model utilized has been revised multiple times over the course of the 834 pandemic and differs substantially in methodology, assumptions, range of predictions and 835 quantities estimated. Second, the IHME has been irregular in publishing the downloadable 836 estimates online for some periods. Third, we model the death estimates by date of reporting 837 rather than by the date of death. Lastly, the unpredictable social component of the epidemic on 838 ground is also a limiting factor for the study as we do not know the ground truth death pattern 839 when the forecasts are generated. 840

841
In conclusion, the reproduction number has been fluctuating around ~1.0 from end of July-end of 842 September 2020, indicating sustained virus transmission in the region, as the country has seen 843 much lower mobility reduction and mixed reactions towards the stay-at-home orders. Moreover, 844 the spatial analysis indicates that states like Mexico, Michoacán, Morelos, Nuevo Leon, Baja 845 California need much strict public health measures to contain the epidemic. The GLM and sub-846 . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

(which was not certified by peer review)
The copyright holder for this preprint this version posted January 13, 2021. ; https://doi.org/10.1101/2021.01.11.21249561 doi: medRxiv preprint epidemic model applied to mortality data in Mexico provide reasonable estimates for short-term 847 projections in near real-time. While the GLM and Richards models predict that the COVID-19 848 outbreak in Mexico and Mexico City may be on a sustained decline, the sub-epidemic model and 849 IHME model predict a stabilization of daily deaths. However, our forecasts need to be 850 References 870 . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

(which was not certified by peer review)
The copyright holder for this preprint this version posted January 13, 2021. ; https://doi.org/10.1101/2021.01.11.21249561 doi: medRxiv preprint CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) The copyright holder for this preprint this version posted January 13, 2021. ; https://doi.org/10.1101/2021.01.11.21249561 doi: medRxiv preprint CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) The copyright holder for this preprint this version posted January 13, 2021. CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) The copyright holder for this preprint this version posted January 13, 2021. ; https://doi.org/10.1101/2021.01.11.21249561 doi: medRxiv preprint 1142 1143 Fig 6: Systematic comparison of the six models (GLM, Richards, sub-epidemic model, IHME 1144 current projections (IHME C.P), IHME universal masks (IHME U.M) and IHME mandates 1145 easing (IHME M.E)) to predict the cumulative COVID-19 deaths for Mexico in the thirteen 1146 sequential forecasts. The blue circles represent the mean deaths and the magenta vertical line 1147 indicates the 95% PI around the mean death count. The horizontal dashed line represents the 1148 actual death count reported by that date in the November 11, 2020 IHME estimates files. CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) The copyright holder for this preprint this version posted January 13, 2021. ; https://doi.org/10.1101/2021.01.11.21249561 doi: medRxiv preprint projections (IHME C.P), IHME universal masks (IHME U.M) and IHME mandates easing 1159 (IHME M.E))to predict the cumulative COVID-19 deaths for the Mexico City in the thirteen 1160 sequential forecasts. The blue circles represent the mean deaths and the magenta vertical line 1161 indicates the 95% PI around the mean death count. The horizontal dashed line represents the 1162 actual death count reported by that date in the November 11, 2020 IHME estimates files. CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) The copyright holder for this preprint this version posted January 13, 2021. ; https://doi.org/10.1101/2021.01.11.21249561 doi: medRxiv preprint CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) The copyright holder for this preprint this version posted January 13, 2021. ; https://doi.org/10.1101/2021.01.11.21249561 doi: medRxiv preprint CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) The copyright holder for this preprint this version posted January 13, 2021. CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) The copyright holder for this preprint this version posted January 13, 2021. ; https://doi.org/10.1101/2021.01.11.21249561 doi: medRxiv preprint