INTRODUCTION

The SARS-CoV-2 coronavirus, which caused the COVID-19 disease and was originated in China in 2019, has dramatically impacted the world. It caused many infections and deaths.

During the first five months after the World Health Organization (WHO) declared the pandemic, Mexico became the fourth location with the most deaths caused by this coronavirus, at least during the first wave of infections. Mexico reported 79, 791 deaths up to 27 August 2020, and 592, 495 infected cases as of the same date (^{Worldometer, 2021}). Mexico, like other countries, adopted non-pharmaceutical interventions (e.g., social-distancing measures) to reduce the spread of the COVID-19 pandemic and its effects; that is why it is crucial to identify factors that increase the spread of the disease.

For the reasons mentioned above, it is unsurprising that since the start of the pandemic, many researchers have turned their attention to exploring the different risk factors for both the spread of the disease and its mortality rate.

In the absence of information regarding the impact of mobility on the spread of the disease, the risk of contagion remains only as a perception for individuals, causing them to avoid using some public transport modes or to ignore the recommendations by governments (^{Chen et al., 2022}). ^{Hörcher, Singh and Graham (2022)} affirm that mixed messages have been disseminated to travelers on the risk of infections and safety in public transport, so they discussed five potential approaches to social distancing in public transport and concluded that multiple demand management measures would have to be implemented simultaneously.

Motivated by the lack of research in terms of mobility and the spread of COVID-19, this paper explores the association of five modes of transportation as well as other poverty indicators with the incidence rate of COVID-19 in the Metropolitan Area of the Valley of Mexico (MAVM) within the period that includes the first wave of the outbreak. The purposes of this work are two. The first one is to identify and quantify the effect of public transportation modes on the incidence rate of COVID-19 infections in a hypothetic scenario where the trips using public transport were not affected. The second purpose is to identify and quantify the comorbidities that are more related to the death rate of COVID-19 in a zone of the Metropolitan Area of the Valley of Mexico during the first wave of the coronavirus pandemic.

In this paper two hypotheses were tested: if the number of arrivals to each zone, using any transit mode considered, is significantly associated with the incident rate of cases, and if the comorbidities are significantly associated with the rate of death caused by COVID-19.

This study is made on a borough level considering 26 municipalities of the State of Mexico and the 16 boroughs of Mexico City. A multi-linear regression model was used to analyze and quantify the effect of the control variables and data from the Health Ministry. The number of reported cases of infection and death due to COVID-19 was considered since the first case in Mexico City up to August 27, 2020.

It is important to emphasize that due to the shortage of mobility data corresponding to the study period, a hypothetical scenario was considered in which no containment measures would affect the number of trips made in the study area and then visualize the results that this situation would lead us to. Although the results obtained in this work should be interpreted under that consideration, this study could be relevant to gain experience and identify the factors that are most important in the COVID-19 incidence and have greater control in the future.

This article is organized as follows: first, the literature review is presented, where some studies that were carried out using data from different areas of Mexico were incorporated. Next, the methodology used in this research is described, indicating the sources of the data used and redefining the variables into consideration. After that, the two models proposed to investigate the relationship between variables through a functional form were introduced and the three assumptions of a linear regression in the two models studied were reviewed, namely, normality, homoscedasticity and non-autocorrelation of the errors. Then, the results obtained from the models are presented. After that, a discussion of the results is done, and some conclusions are made. Finally, this work is complemented with Appendix A, where local spatial clusters are shown using a bivariate local indicator spatial association analysis (LISA).

LITERATURE REVIEW

Among the different available investigations on the identification of factors that impact COVID-19 incidence and deaths through the first wave of the coronavirus pandemic, the work of ^{Sannigrahi et al. (2020)} can be found, where spatial regression models (geographically weighted regression (GWR), spatial error model (SEM), spatial lag model (SLM), ordinary least square (OLS), partial least squares (PLS) and principal components regression (PCR) are used. There, the association between income, poverty, and total population, is analyzed with the number of people infected and the number of deaths due to COVID-19. There are 31 European regions considered as a study area with data up to 29 April 2020.

In ^{Ehlert (2021)} the association of socioeconomic, demographic, and health-related variables at the regional level with COVID-19-related cases and deaths in Germany through the first wave of the coronavirus pandemic is explored. A comparison of spatial auto-regressive SAR, SEM, and spatial auto-regressive combined (SAC) models is done there. Also, there are some explanatory variables considered like income, age, employment, population, gender, heart conditions, such as chronic obstructive pulmonary disease (COPD), and the number of people in healthcare.

In ^{Stojkoski et al. (2020)}, 31 factors are included, such as employment, age, income, population, hospital beds, and international influence on social interactions to explain the outcome of the first wave of the coronavirus pandemic.

In other studies, like the one by ^{Irawan et al. (2022)}, the relationship between different variables of travel behavior and individual characteristics is studied. The role of the behavior of humans in protecting themselves from contracting COVID-19 disease in Indonesia is particularly analyzed. A structural model is proposed, which was adjusted using SEM. It is shown that the more people perceive the severity of the disease, the more they try to reduce their activities outside the home. Table 1 summarizes the models and some variables considered in each study mentioned above.

MATERIALS AND METHODS

Mobility in the area considered

This work was focused on the conurbation around Mexico City, known as the Metropolitan Area of the Valley of Mexico (MAVM), described by the National Population Council (^{CONAPO, 2015}). It comprises the 16 boroughs of Mexico City, 59 municipalities of the State of Mexico, and one municipality of the state of Hidalgo. The MAVM is the most populated area in the country, it is home to more than 18 million inhabitants, and has the highest concentration of resources and national institutions. According to the National Institute of Statistic and Geography (INEGI), in 2020, around 21.8 million people lived in the MAVM in 7,866 km². It has many offices and an extensive transport network, with more than 20 different modes of transport.

Although the federal government made available a database regarding the COVID-19 disease during the first year of the pandemic, the database was incomplete or had some inconsistencies. Due to the unavailability of data for some explicative variables considered in our model, 34 municipalities were not considered in this study. Figure 1 (a) shows the MAVM; the remaining 42 boroughs considered can be seen in (b), 16 from Mexico City (in white) and 26 from the State of Mexico in (light orange).

Source: Author´s own elaboration.

Figure 1 Area of study 

The origin-destination movements in Mexico City are usually obtained every ten years from surveys; the last was carried out in 2017 (^{INEGI, 2017}). According to this survey, in 2017, about 52.7 million trips were carried out in the MAVM on a weekday. During the first months of the pandemic, there was great interest by the Mexico City government in estimating the reduction in trips made in the city using some modes of transportation; however, it was not clear what the reference value was, and there was no information about all modes of transportation or the surrounding municipalities.

This paper considered the information obtained from the 2017 origin-destination survey carried out by the INEGI. Also, it is assumed that no significant changes were made until 2020. So, in this analysis, a hypothetical scenario is contemplated where the containment measures did not affect the number of trips. The number of weekly trips for each transport mode that originated in the MAVM and ended at each of the municipalities in the study area is used.

In the MAVM more than twenty different transport modes operate; however, using the selection method forward and backward stepwise, the final model includes collective/micro, bus, walk, taxi, and subway transit modes. Here, collective/micro refers to small public transport vehicles such as microbuses and other vehicles with a maximum capacity of 18 people.

Dependent variables

The first case of COVID-19 in Mexico was reported on February 27, 2020; since then, the Health Ministry daily reported the number of infected people. Likewise, the reported cases were tracked at the municipality level nationwide. In this work, each municipality of the MAVM was considered as the observation unit; however, as mentioned before, because of the absence of data for some variables included in the model, 34 municipalities in the area were eliminated from the study, leading to 42 municipalities for the analysis. The study period for the number of people infected with COVID-19 comprehends from the first day of contagion in Mexico until August 22, 2020. The reported deaths due to the infection are considered from the first death in MAVM to August 22, 2020.

One of the proposed regression models contemplates the incidence rate as a dependent variable; it is the number of people infected with COVID-19 in the 42 municipalities mentioned above during the described period. This incidence rate is constructed by dividing the total number of people infected with COVID-19 (from February 27 to August 22, 2020) by the total population of each municipality, provided by the ^{INEGI (2020b)}. Similarly, the other proposed regression model contemplates the rate of death caused by the disease, given by the number of deaths reported divided by the total population of each municipality. Figure 2 shows the spatial distribution of the incidence and death rates in the study area. In both cases, the high levels of infected and dead people are focused on the north of Mexico City and the south of the State of Mexico.

Source: Author´s own elaboration.

Figure 2 COVID-19 incidence and death rates classified in quartiles 

THE MATHEMATICAL MODEL

In the regression analysis, the relationship between variables through a functional form is investigated. This relationship is approximated by a model of the form:

where Y is the dependent variable (output), X₀, X₁, …,X_n are the independent variables (explanatory), and ε is a random error representing the discrepancy in the approximation. Next, an example of a functional form is shown:

where β₀, β₁, …, β_n are the coefficients/parameters to be estimated from data. The estimation allows evaluation of the importance of each explanatory variable by knowing its effects on the dependent variable when changes are present. The most used method to estimate these parameters in regression analysis is the ordinary least squares (OLS) method, which, under specific hypotheses, produces estimates with desirable properties. Other superior methods are used when any of those hypotheses are not met. An important issue in regression analysis is to specify the functional form that will relate the response variable to the set of explanatory variables. Functional forms can be linear in the parameters, which, in turn, can be classified as linear or non-linear in the variables.

In this work, simple functional forms for the models were proposed. That allowed us to interpret the estimated coefficients without significant difficulties. The method used to estimate the parameters in each proposed model is ordinary least squares.

Analysis of linear regression assumptions

The three assumptions of a linear regression in the two models studied were reviewed, namely, normality, homoscedasticity, and non-autocorrelation of the errors (^{Rao, 2009}). First, it was verified that model (2) does not present multicollinearity problems obtaining a coefficient of determination of 0.9532. The F-value is 108.6 with p-value lower than 2.2e-16. The hypotheses under which the linear regression model is governed using OLS were also verified. Through the Breusch-Pagan test (^{Kennedy, 2003}), no evidence of heteroscedasticity was found at a confidence level of 99.9% and the value of the Breusch-Pagan statistic equal to BP=5.2889. The normality of the errors was verified using the Shapiro-Wilk test (^{Kennedy, 2003}) at a confidence level of 99.9% and with the value W=0.97681 for the Shapiro-Wilk statistic. Finally, no evidence of error correlation was found using the Durbin-Watson test at a confidence level of 99.9% and the Durbin Watson statistic equal to DW = 1.7187 (Chaterjee and Hadi, 2012). The results were visually corroborated by Figure 3.

Source: Author´s own elaboration.

Figure 3 Graphic diagnosis for model (2) 

For model (4), R² = 0.8544 and the F-value =42.26 with p-value < 4.39e-14 were obtained. Through the Breusch-Pagan test, no evidence of heteroscedasticity was found at a confidence level of 99.9% and the value of the Breusch-Pagan statistic equal to BP=6.1759. The normality of the errors was verified using the Shapiro-Wilk test at a confidence level of 99.9% and with the value W=0.98452 for the Shapiro-Wilk statistic. Finally, no evidence of error correlation was found using the Durbin-Watson test at a confidence level of 99.9% and the Durbin Watson statistic equal to DW = 2.0388. The results were visually corroborated by Figure 4.

Source: Author´s own elaboration.

Figure 4 Graphic diagnosis for model (4) 

RESULTS

With the methodology described above, the results were interpreted and the effect of the significant explanatory variables in the incidence of COVID-19 was inferred.

Table 4 shows the coefficients estimated by OLS. There can be seen for the three considered modes of transport in the model (2): the In(t_cm) has a positive coefficient, indicating that the greater the number of trips made using collective/micro, the greater the increase in the incidence rates. More precisely, its coefficient βtcm=0.14560, indicates that for every 1% of increase in t_cm, the incidence rate, I_r, increases by 0.14560%. In(t_b) also has a positive sign and its coefficient, βtb=0.03485, indicates that for every 1% of increase in the number of trips made using a bus t_b, the incidence rate, I_r, increases 0.03485%. In addition, In(t_w), had a negative coefficient, which means that the number of trips done by walking negatively affects the incidence rates of COVID-19, that is, its coefficient, βtw=-0.22369, indicates that for every 1% increase in t_w, the incidence rate decreases by 0.22369%. The transport mode collective/micro has the highest effect on the I_r for the transport mode-related variables since it has the biggest coefficient. The variable p_m was significant with a positive coefficient, indicating that the higher the percentage of the male population in the municipality, the higher level of incidence rate was observed. Its coefficient, βpm=0.09683, indicates that for each unit increase in pm the incidence rate, I_r, increases by 9.683%. The variable In(p_d) was significant with a positive coefficient βpd=0.04170, which indicates that the higher the population density, the higher the incidence rate of COVID-19, and for every 1% increase in p_d, the incidence rate, I_r, increases 0.04170%.

Regarding the three comorbidity variables considered, all significantly positively affect the incidence rate I_r. The coefficient βhp=0.51038 indicates that, for every 1% increase in h_p , the incidence rate increases by 0.51038%. The variable In(h_a) has a coefficient βha=0.26514, which indicates that for every 1% increase in h_a, the variable I_r increases 0.26514%. For the variable In(h_t), its coefficient βht=0.19527, indicates that for every 1% increase in h_t, there is an increase of 0.19527% in I_r. It can be observed that the effect of the prevalence of pneumonia on the incidence rate of COVID-19 is greater than the effect of the prevalence of hypertension, and this in turn has a greater effect than the prevalence of smoking.

As well as model (2), model (4) does not present multicollinearity problems. In addition, homoscedasticity, non-autocorrelation, and normality were investigated and verified using the same tests used for the model (2). This allows us to make inference and quantify the effects of significant independent variables on the dependent variable. With model (4), it is proven that the comorbidity variables, In(h_p) and In(h_d) are significantly associated with the COVID-19 death rate.

Table 5 shows the estimated coefficients by OLS for model (4); also in this case, all the variables are significant. The coefficient of ln(h_p), δhp=0.28116, indicates that for every 1% increase in pneumonia prevalence rate, h_p, the death rate, D_r, increases by 0.28116%. While the coefficient of the variable h_d, δhd=1034.71408, indicates that for each unit increase in diabetes prevalence rate, h_d, the rate of deaths due to COVID-19 disease increases by 103471.08%. Regarding the variable ln(p_dw), its coefficient δpdw=-1.26979, indicates that for a 1% increase in p_dw, the variable D_r decreases 1.26979%. The population density variable with a coefficient δpd=0.05472, indicating that for every 1% of increase in population density, there is an increase of 0.05472% in the rate of deaths due to COVID-19 disease. Also, it can be seen that p_w is significant and its coefficient δpw=0.01460 indicates that for every unit of increase in the percentage of the population that earns less than twice the minimum wage, the rate of death due to COVID-19 disease increases by 1.460%.

DISCUSSION AND CONCLUSIONS

This paper examined both associations of transport variables with the incidence rate of COVID-19 and the comorbidity variables with the COVID-19 death rate. The two studies were carried out during the first wave of infections. The hypotheses were tested using mathematical models with log-transformed variables. It was found that the number of trips to the municipalities using some modes of transport is significantly associated with COVID-19 incidence rates, and some comorbidities are significantly associated with COVID-19 deaths. The signs of the significant variable coefficients in both models resulted as expected.

It is important to highlight that the magnitude of the transport variable coefficients is consistent with reality. For example, the fact that the variable related to the collective/micro transport mode will result in a higher coefficient than the coefficient of the variable t_b could be explained by saying that the public transport vehicles corresponding to the collective/micro mode have no control over their capacity because people can travel very close to each other and stand for long periods, so the social distancing is almost impossible in them; contrary to the bus mode of transport, where all passengers are assigned to a seat. The problem with the bus mode of transportation is that they have no natural ventilation, making it easy for the pathogen to spread. We also find a negative coefficient for the transport variable related to the trips made by walking t_w, which is also consistent with reality since it allows social distancing and natural ventilation, so this mode of transport contributes to the reduction in the incidence rate of COVID-19. Comparing the effects of the transport variables quantitatively, it turns out that the impact of the variable related to the trips made by walking is more significant than the other variables related to collective/micro and bus transport modes.

Regarding the comorbidity variables, to our knowledge, it is a novel finding that not only the comorbidities are related to mortality due to COVID-19, but they are also risk factors for contracting the disease. This may be because people with these types of diseases go to hospitals more frequently; therefore, they have greater exposure to people infected with COVID-19. It has been shown that more men than women are dying (of COVID-19), potentially due to sex-based immunological or gendered differences, such as patterns and prevalence of smoking (^{Wenham et al., 2020}), here our estimation on the coefficient of the variable p_m suggests that men are more likely to be infected than women.

^{Ejaz et al. (2020)} showed that comorbidities lead the COVID-19 patient into a vicious infectious cycle of life and are substantially associated with significant morbidity and mortality; thus, our findings on the association of comorbidities with the increase in the rate of COVID deaths are consistent with the studies carried out. Also, in ^{Gasmi et al. (2021)}, it is shown that comorbidities contribute to the prognosis of acute disease and the increased risk of severe symptoms, where it is observed that 70% of the patients in intensive care have comorbidities. In our study, it is shown that the prevalence variables of diabetes and pneumonia are associated with the increase in the rate of COVID-19 deaths, the rest of the comorbidities do not appear in the model since they present multicollinearity.

With respect to the home occupancy rate variable, it is consistent with reality since the lower the saturation in dwelling homes, the less probability of contagion is expected, and therefore, it could imply less mortality.

Two machine learning models were also fitted to justify the models presented in this paper. The models based on ordinary least squares (OLS) were compared to neuronal network regression and principal component regression (PCA) models. The coefficient of determination R² was used as a performance measure. It can be consulted at ^{Hernández and Wences (2023)} for a more detailed explanation of the comparison analysis used here.

Concerning model (2), for each machine learning algorithm, 10, 30, 50, 100, 300, 500 and 1000 iterations were performed. In Figure 5, the box plot for the performance measure is shown.

Source: Author´s own elaboration.

Figure 5 Box Diagram for the R² in the models of incidence rate 

It can be observed that the models are equivalent with respect to the performance measure since the whiskers in the box plot intersect each other. Figure 6 presents the graphs of the times, as a function of the number of iterations, for each model. The graphs for principal component regression (PCA) and OLS regressions overlap, and the computational cost is considerably cheaper for PCA and OLS than for neuronal network regression.

Source: Author´s own elaboration.

Figure 6 Computational cost for each model in incidence rate 

For model (4), it is also proven that the two machine learning models, neuronal network regression and principal component regression, are equivalent to model (4), see Figure 7. Figure 8 shows the computational cost for the three models. The results are similar to the previous case.

Source: Author´s own elaboration.

Figure 7 Box Diagram for the R² in the models of death rate 

Source: Author´s own elaboration.

Figure 8 Computational cost for each model in death rate 

Let us recall that this study was carried out considering a hypothetical scenario where the containment measures did not affect the number of trips made in the study area due to the difficulty of obtaining updated data. Although this could be a limitation in our research, the importance of this analysis is to assess the effects that would have occurred if mobility in the MAVM would not had been restricted. Also, let us notice that, according to ^{Prieto et al. (2022)}, the mobility trends were similar for all the transit modes considered in this work. So, it may be expected that if current data were available, the study could be repeated, and the results may be proportional to the observed data.

Despite the limitations of this work, the models presented had a good fit with high determination coefficients. They fulfilled the multilinear regression hypotheses, which allowed us to perform statistical inference since the ordinary least squares estimators are of minimum variance in the whole class of unbiased estimators. So, through models, which could be considered simple, the change in the incidence and death rates of COVID-19 in the study area has been explained. Future work that could be carried out with available data would be to consider variables related to environmental factors in a model and explain the indirect impact on incidence and death rates.

This study is complemented by Appendix A, where local spatial associations between the COVID-19 incidence/mortality rate variables are shown. Each variable was found to be significant in the proposed models by computing the bivariate local Moran’s I value. It is worth mentioning that, in this work, we use the I Moran test to determine the spatial autocorrelation, generating the matrix of spatial weights using the Queen criterion. Geographically speaking, two municipalities are neighbors if they share a point on their borders; however, the test returned a p-value beyond 0.01. This means that the hypothesis of non-autocorrelation of errors is not rejected, so we did not find enough evidence of spatial autocorrelation. Other criteria could be sought for future work to build the weight matrix and investigate spatial autocorrelation. For example, in ^{Naveen and Gurtoo (2022)}, passenger clusters are considered based on their travel dynamics to guide those responsible for transport planning to provide a more efficient service, in the presence of COVID-19, rather than control access. In ^{Iacus et al. (2021)}, they define the concept of Mobility Functional Areas as the geographic zones highly interconnected according to the analysis of mobile positioning data and affirm that they can be used to inform local transportation, spatial planning, and health and economic policies.