Geospatial Variability in Excess Death Rates during the COVID-19 Pandemic in Mexico: Examining Socio Demographic, Climate and Population Health Characteristics.

OBJECTIVES: This study examined how socio-demographic, climate and population health characteristics shaped the geospatial variability in excess mortality patterns during the COVID-19 pandemic in Mexico.
METHODS: We used Serfling regression models to estimate all-cause excess mortality rates for all 32 Mexican states. The association between socio-demographic, climate, health indicators and excess mortality rates were determined using multiple linear regression analyses. Functional data analysis characterized clusters of states with distinct excess mortality growth rate curves.
RESULTS: The overall all-cause excess deaths rate during the COVID-19 pandemic in Mexico until April 10, 2021 was estimated at 39.66 per 10 000 population. The lowest excess death rates were observed in southeastern states including Chiapas (12.72) and Oaxaca (13.42), whereas Mexico City had the highest rate (106.17), followed by Tlaxcala (51.99). We found a positive association of excess mortality rates with aging index, marginalization index, and average household size (P < 0.001) in the final adjusted model (Model R2=77%). We identified four distinct clusters with qualitatively similar excess mortality curves.
CONCLUSION: Central states exhibited the highest excess mortality rates, whereas the distribution of aging index, marginalization index, and average household size explained the variability in excess mortality rates across Mexico.

Introduction

Monitoring all-cause excess mortality, above an expected level of total deaths, as a pandemic unfolds is one of the key ways to evaluate its mortality impact (Weinberger et al., 2020). All-cause excess mortality estimates include deaths that are directly or indirectly attributed to the pandemic (CDC, Serfling, 1963). Besides direct deaths due to COVID-19, deaths indirectly attributed to the COVID-19 pandemic include those related to denied or delayed care for acute emergencies (Maringe et al., 2020, Schirmer et al., 2020) or other chronic conditions (Douglas et al., 2020), the disruption of routine health care services in an overburdened health care system (Roberton et al., 2020), unaddressed mental health concerns including suicide and self-harm (Kawohl and Nordt, 2020, Sahoo et al., 2020), and drug overdoses (Currie et al., 2021). Detailed analyses of excess mortality can inform prevention and mitigation strategies by determining where the mortality impact of the pandemic has been most significant.

Mexico is one of the countries in Latin America that is bearing the brunt of the COVID-19 pandemic with the fourth-highest number of COVID-19 deaths in the world, after the USA, Brazil, and India, as of late June 2021 (Statistica, 2021). In fact, Mexico has reported a total of 2 487 747 (1.38% of global cases) confirmed cases of COVID-19, including 231 847 deaths (5.96% of global deaths), as of June 25, 2021 (WHO, 2021). A previous study reported a high all-cause excess death rate of 26.10 per 10 000 population in Mexico in 2020, with COVID-19 deaths accounting for only 38.64% of the estimated excess deaths (Dahal et al., 2021). Additionally, Mexico was identified as one of the countries with highest excess deaths in terms of absolute numbers, excess deaths per 100 000 population, and excess deaths as percent of annual deaths in recent research (Karlinsky and Kobak, 2021). While the relatively low proportion of COVID-19 deaths, out of all excess deaths, could be the result of low testing rates, misclassification of COVID-19 deaths, and delays in reporting COVID-19 deaths (Gutierrez et al., 2020), a substantial number of deaths during the pandemic could be due to the indirect causes (CDC, 2021) and need to be examined in more depth.

The distribution of indirect causes of deaths depends on several factors such as sociodemographic characteristics, population health and the selection, timing and intensity of any public health interventions, in addition to the efficiency and reach of the health and social care system (Kontis et al., 2020). In Mexico, pandemic control measures have varied widely (Knaul et al., 2021). Therefore, a more detailed understanding of the mortality burden of the pandemic can be obtained by quantifying spatial heterogeneity in excess deaths at a state level and by examining the influence of underlying sociodemographic, economic, and health system related factors, and also climate factors. In this study, we pose the question whether in a country such as Mexico, with very high COVID-19 mortality, potential spatial variability in the excess deaths can be explained by underlying sociodemographic, climate and population health indicators. To answer that research question, we first estimated the all-cause excess deaths during the COVID-19 pandemic in Mexico comprising 31 states and Mexico City. Next, we evaluated the potential associations between different sociodemographic factors, climate, and excess mortality patterns at the state level in Mexico. Furthermore, we also conducted a cluster analysis to characterize the shapes of the excess mortality curves into different groups that describe the potential geospatial variability in excess mortality. Analyses such as these are critically important for understanding excess mortality and for guiding intervention strategies.

Methods

Data

We obtained weekly all-cause death counts updated on May 25, 2021 for Mexico at the state level and for Mexico City, based on epidemiological weeks from January 2020 until April 10, 2021 and for the preceding 5 years (2015–2019) to establish a baseline mortality level (Government of Mexico). We accessed publicly available weekly mortality data from the National Institute of Statistics and Geography (INEGI) for the years from 2015 to 2018, and data from National Population Registry (RENAPO) for the years 2020 and 2021(Government of Mexico). For the year 2019, we chose either INEGI or RENAPO as the data source, based on the value of weekly mortality of the last week of 2018 and the first week of 2019 for each state. We obtained the national and state-level population size estimates from the National Population Council (CONAPO) of Mexico (CONAPO, 2015-2030). Mortality data was not available for the state of Tlaxcala for the last six weeks of the study period. For this reason, this state was excluded from our regression and functional cluster analyses.

For each state, including Mexico City, we obtained data on seven variables: population density (2020), aging index (2020), average household size (2020), marginalization index (2020), rate of new cases of depression per 100 000 population (2019), public spending on health as percent of GDP (2019), and climate zone. Data on population density, aging index, average household size, and rate of new case of depression were obtained from INEGI (INEGI), data on public spending on health was obtained from the subsystem of health accounts at the federal and state level (SICUENTAS) (General Directorate of Health Information, 2021), and the data on the marginalization index was available from CONAPO (CONAPO, 2020). To model climate variation, we used the Köppen-Geiger classification system (Méndez-Arriaga, 2020) which divides Mexican states into three climatic groups: A, B, and C as follows:

Group A: Tropical/megathermal climates: warm humid climate; warm sub-humid climate; warm, semi-warm humid climate; and semi-warm sub-humid climate.

Group B: Dry (desert and semi-arid) climate: dry, warm and semi-dry climate; dry, temperate and semi-dry semi-cold climate; dry, warm dry climate; dry, temperate dry climate; dry, temperate dry winter rains; dry, warm very dry climate; dry, temperate and very dry semi-cold climate.

Group C: Temperate/mesothermal climates: temperate, semi-warm humid climate; temperate, semi-warm sub-humid climate; temperate, humid climate; temperate, sub-humid climate; temperate, semi-cold humid climate; temperate, semi-cold sub-humid climate. According to this grouping group A, B, and C include 11, 14, and 7 states respectively. States in Mexico based on this climate categorization are presented in Supplemental Figure 2.

Summary statistics of these variables from 31 states and Mexico City are provided in Table 1 .

Table 1

Descriptive statistics for six continuous variables included in the multiple regression analysis of excess mortality in Mexico (n=32).

Variable	Minimum	Mean (SD)	Median (IQR)	Maximum
Population density (2020) (habitants per km2)	10.80	309.68 (1078.69)	67.15 (127.20)	6163.30
Aging index (2020)	28.70	46.41 (10.44)	45.45 (7.40)	90.20
Average household size (2020)	3.6	3.91 (0.18)	3.90 (0.20)	4.4
Marginalization index (2020)	11.32	18.89 (2.73)	19.43 (3.23)	23.01
Rate of new case of depression per 100,000 population (2019)	22.06	114.88 (76.05)	92.14 (67.37)	348.17
Public spending on health as a percent of GDP (2019)	0.94	3.12 (1.02)	2.96 (1.32)	5.81

Pandemic period and excess deaths

For both the national data and the data for each state, we separately estimated the baseline mortality level by fitting cyclical Serfling regression models to all-cause deaths in the non-COVID-19 period, after excluding data from March 2020 to April 2021 by employing established methodology (Chowell et al., 2014, Chowell et al., 2012, Dahal et al., 2018a, Serfling, 1963, Viboud et al., 2013). Details on the model equation that was used can be found in (Dahal et al., 2021). After establishing a weekly baseline and the corresponding 95% CI at the national level, we defined the periods of COVID-19 pandemic as the weeks in 2020 and 2021 where the observed all-cause mortality rate at the national level in Mexico exceeded the upper 95% confidence limit of the national baseline mortality rate. The excess mortality rate was estimated at the state level and for Mexico City for the same defined period of the COVID-19 pandemic. Excess all-cause mortality rate was estimated as the difference between the observed and model adjusted baseline mortality rates for each week constituting the pandemic period. The overall pandemic excess mortality in 2020 and 2021 was calculated by summing the excess death rates across the pandemic weeks in the given year (Chowell et al., 2014, Dahal et al., 2021, Dahal et al., 2018a). Negative excess mortality estimates were replaced by zeros in our analyses to account for underreporting due to reporting delays (Aron and Muellbauer, 2020, CDC, 2021).

Multiple regression analysis

After estimating the total excess mortality rate for each state, we explored the association between the total excess mortality rate and the predictor variables. Because the population density and rate of new case of depression distributions were skewed, we transformed these variables to log base 10. Since we identified Mexico City as a potential influential point, we performed sensitivity analysis by comparing the results of different models, including and excluding Mexico City. Since there was no significant change in the statistical inference of the parameters, we included Mexico City in the model, and parameters were estimated using ordinary least squares (OLS) method.

Cluster analysis

We followed the analytic methods described in (Srivastava and Chowell, 2020) to pre-process the weekly cumulative all-cause excess deaths for 30 states and Mexico City (excluding Tlaxcala, refer to study setting for details). Then, we analyzed the shapes of the excess all-cause death rate curves to compare, cluster, and summarize growth rates.

We employed the following steps to smooth and normalize the weekly all-cause excess death data:

Smoothing: Cumulative excess deaths curves were smoothed using the smooth function in Matlab, which uses a moving average filter over a 10-week span.

Time differencing: If denotes the given cumulative number of excess deaths for state i on week t, then per week growth at time t is given by .

Re-scaling: We rescaled each curve by dividing each by the total excess number of deaths for a given state i, which is equivalent to computing where and K is the number of weeks in the period.

Smoothing: We then smoothed the normalized curves over a 5-weeks span, using the smooth function in Matlab.

To identify the clusters by comparing the curves, we used a simple metric. For any two rate curves, h and h, we compute the norm ||hi −hj||, where the double bars denote the L norm of the difference function, i.e., ||hi −hj|| = which is approximated by , where K is the number of weeks in the period.

To perform clustering of thirty-one curves into smaller groups, we applied the dendrogram function in Matlab using the “Ward's” linkage as explained in (Srivastava and Chowell, 2020). The Ward's linkage minimizes the total within-cluster variance and tends to produce more compact clusters. It is also less sensitive to outliers than other linkages. The number of clusters was decided empirically by inspecting the overall clustering results. After clustering the states into different groups, we derived the average curve for each cluster using a time wrapping algorithm (Srivastava and Chowell, 2020, Srivastava and Klassen, 2016).

Results

From March 1, 2020 to April 10, 21 (total of 58 weeks), the observed death rate was greater than the upper 95% confidence interval of the baseline starting from week of April 12-18, 2020 until the week of April 4-10, 2021 (total of 52 weeks) (Figure 1 ). For this period starting from April 12, 2020, the first peak of weekly excess mortality rate occurred in the week of July 12-18, 2020, with the excess death rate of 1.04 per 10 000 population, then declined slightly for a few weeks and then increased again from the week of December 20-26, 2020, and reached a peak with an all-cause excess death rate of 1.99 per 10 000 population in the week of January 17-23, 2021. The excess death rate remained below 0.5 from the week of February 28, 2021 until the end of the study period.

Figure 1

Mortality rate per 10 000 population, Mexico, January 2015–March 2021. The black curve is the observed weekly death rate. The grey curve is the predicted baseline death rate. Square dotted curves indicate the upper and lower 95% confidence intervals of the baseline death rate. The long-dashed line indicates the COVID-19 pandemic period.

All-cause excess death rates for the national level, Mexico City, and 31 states of Mexico are presented in Table 1. The map displaying state-level estimates is depicted in Figure 2 . While the total excess death rate in Mexico was at 39.66 per 10 000 population, equivalent to a total of ∼508 289 excess deaths, the excess mortality rate in Mexico City was the highest and estimated at 106.17 per 10 000 population (∼95 690 total number of excess deaths). Among 31 states, Tlaxcala (51.99), Morelos (45.90), Puebla (45.12), and Mexico (44.43) were among the states with the highest excess mortality rates. The states with the lowest death rates included Chiapas (12.72), Oaxaca (13.42), Quintana Roo (19.41), and Yucatan (21.11) (Table 2 ). Only one state, Chiapas, had no excess deaths in 2021. COVID-19 accounted for only 42.16% of total excess deaths at the national level ranging from 20.97% in Chiapas to 76.05% in Quintana Roo.

Figure 2

Map depicting excess death rate per 10 000 population by state in Mexico.

Table 2

Estimates for all cause excess mortality rate by state per 10 000 population during COVID-19 pandemic in Mexico, March 1, 2020-April 10, 2021.

State/Region	All cause excess death rate per 10 000 population (March 1, 2020-April 10, 2021)	All cause excess death rate per 10 000 population in 2020 (Includes weeks starting from March 1, 2020, to December 27, 2020)	All cause excess death rate per 10 000 population in 2021 (Includes week starting on January 3, 2021, to April 4, 2021)	Total number of all cause excess deaths	COVID-19 deaths (Percentage of all cause excess deaths)
National	39.66	27.25	12.41	508 288.78	214 298 (42.16)
Aguascalientes	29.37	21.61	7.76	4228.25	2300 (54.39)
Baja California Sur	26.13	16.36	9.77	2118.75	1326 (62.58)
Baja California	38.55	31.78	6.77	14 051.54	8048 (57.27)
Campeche	22.25	21.02	1.22	2228.05	1183 (53.09)
Chihuahua	27.76	25.37	2.39	10 561.33	6467 (61.23)
Chiapas	12.72	12.72	0	7287.73	1528 (20.97)
Mexico City	106.17	62.93	43.24	95 689.73	32 166 (33.61)
Coahuila	38.33	30.88	7.45	12 369.31	6189 (50.03)
Colima	28.32	19.23	9.09	2234.78	1158 (51.82)
Durango	27.48	22.94	4.53	5142.23	2372 (46.13)
Guerrero	29.96	21.51	8.45	10 966.69	4231 (38.58)
Guanajuato	36.54	20.44	16.09	22 842.82	10 568 (46.26)
Hidalgo	36.39	22.86	13.52	11 277.96	5995 (53.16)
Jalisco	32.92	18.66	14.26	27 799.37	11 759 (42.30)
Mexico	44.43	28.67	15.76	77 705.35	33 571 (43.20)
Michoacan	31.12	16.69	14.43	15 063.24	5492 (36.46)
Morelos	45.90	24.16	21.75	9428.76	3052 (32.37)
Nayarit	23.64	15.08	8.55	3060.99	1766 (57.69)
Nuevo Leon	36.80	24.98	11.82	20 735.12	9305 (44.87)
Oaxaca	13.42	10.55	2.87	5565.51	3494 (62.78)
Puebla	45.12	36.50	8.62	29 849.45	11 142 (37.33)
Queretaro	35.86	18.79	17.06	8241.92	4063 (49.29)
Quintana Roo	19.41	16.95	2.46	3354.42	2551 (76.05)
San Luis Potosi	32.26	22.29	9.97	9266.51	5190 (56.01)
Sinaloa	31.53	24.97	6.56	9969.06	5927 (59.45)
Sonora	32.22	27.11	5.11	9924.27	6485 (65.34)
Tabasco	22.49	20.34	2.15	5791.32	3994 (68.96)
Tamaulipas	30.33	25.38	4.95	11 086.40	4800 (43.29)
Tlaxcala	51.99	34.93	17.06	7200.90	2367 (32.87)
Veracruz	22.35	17.46	4.89	19 111.79	9524 (49.83)
Yucatan	21.11	16.49	4.62	4781.42	3563 (74.52)
Zacatecas	43.22	29.81	13.41	7217.52	2721 (37.69)

Table 3 shows the results of fitting a taxonomy of multiple regression models of excess mortality rate at the state level in Mexico. To select a final model from among models 4, 6, 7, and 8, we performed a multiple partial F-test. In the multiple partial F-test, we failed to find a significant contribution of adding population density, depression rate, and public expenditure on health on predicting excess mortality rate after accounting for the contribution of aging index, marginalization index, and average household size (F-value3,24=0.39, P-value=0.7631).

Table 3

Results of fitting a taxonomy of multiple regression models of excess mortality rate at the state level in Mexico (n=31).

	Parameter estimate (se)
	Model 1	Model 2	Model 3	Model 4	Model 5	Model 6	Model 7	Model 8
intercept	-19.21* (9.03)	-18.18 (18.09)	-47.03⁎⁎ (13.61)	-250.80⁎⁎⁎(50.67)	-200.36* (74.15)	-215.08⁎⁎ (63.11)	-230.65⁎⁎ (69.68)	-229.07⁎⁎ (72.97)
Aging index	1.12⁎⁎⁎ (0.19)		1.01⁎⁎⁎ (0.18)	1.12⁎⁎⁎(0.15)	1.07⁎⁎⁎ (0.16)	0.99⁎⁎⁎ (0.20)	0.94⁎⁎⁎ (0.22)	0.94⁎⁎⁎ (0.22)
Marginalization index		2.72⁎⁎ (0.95)	1.76* (0.68)	3.30⁎⁎⁎(0.66)	2.79⁎⁎ (0.87)	2.98⁎⁎⁎ (0.74)	2.94⁎⁎⁎ (0.76)	2.87⁎⁎ (1.02)
Average household size				43.47⁎⁎⁎(10.56)	34.57* (14.31)	35.64* (13.39)	37.95* (14.18)	37.95* (14.47)
Tropical/megathermal climate					-5.60 (5.84)
Dry climate					-3.56 (5.00)
log10popdensity						3.79 (3.98)	4.38 (4.17)	4.49 (4.41)
log10depression_rate							4.20 (7.47)	4.24 (7.63)
Public expenditure on health as a percent of GDP								-0.22 (2.29)
Root MSE	10.97	14.41	10.04	8.01	8.18	8.02	8.13	8.30
R2	0.55	0.22	0.63	0.77	0.78	0.78	0.78	0.79
Model F-test	35.18⁎⁎⁎	8.18⁎⁎	24.31⁎⁎⁎	31.09⁎⁎⁎	18.09⁎⁎⁎	23.46⁎⁎⁎	18.34⁎⁎⁎	14.68⁎⁎⁎
(df1, df2)	(1, 29)	(1,29)	(2, 28)	(3, 27)	(5, 25)	(4, 26)	(5, 25)	(6, 24)

P<0.05

P <0.01

P<0.001

We also tested for spatial autocorrelation using Moran I statistics. The result indicated the presence of spatial autocorrelation for the dependent variable (p-value =0.001). However, we failed to find statistically significant Moran I for the residuals for the Model with aging index, marginalization index, and average household size (p-value=0.4133). Hence, we fitted a spatial lag model with three predictors. The lag parameter (Rho) from the spatial lag model was not statistically significant (Rho=0.209, p-value=0.240). In addition, the value of AIC for the lag model (223.22) was slightly higher than that of the OLS model for Model 4 (222.71). Therefore, we chose Model 4 (Table 3) as our final model. Our final model was able to explain 77% of the observed variance in the excess mortality rate (Coefficient of determination (R2)=0.77).

As shown in Table 4 , we found a positive association of excess mortality rate with aging index, marginalization index, and average household size in the adjusted model at 0.05 level of significance.

Table 4

Results for the Final Regression Model 4 of excess mortality rate at the state level in Mexico (n=31).

	Parameter estimate	Standard error	P-value	95% Confidence Limits	Standardized estimate
intercept	-250.80	50.67	<.0001	-354.77, -146.83	0
Ageing index	1.12	0.15	<.0001	0.82, 1.42	0.74
Marginalization index	3.30	0.66	<.0001	1.95, 4.66	0.57
Average household size	43.47	10.56	0.0003	21.80, 65.13	0.48

The result of our clustering analyses is displayed in a dendrogram plot (Supplemental Figure 1). Specifically, we identified the following four prominent clusters based on the shapes of excess growth rate curves at state level:

Cluster 1: Baja California, Coahuila, Guanajuato, Hidalgo, Jalisco, Mexico, Mexico City, Michoacan, Morelos, Nayarit, Nuevo Leon, San Luis Potosi

Cluster 2: Aguascalientes, Chihuahua, Durango, Queretaro, Zacatecas

Cluster 3: Baja California Sur, Colima, Guerrero, Oaxaca, Puebla, Quintana Roo, Sinaloa, Sonora, Tabasco, Tamaulipas, Veracruz, Yucatan

Cluster 4: Campeche and Chiapas

Figure 3 shows the average growth rate curves and one standard deviation band around it. The growth patterns in each cluster are very distinct. For cluster 1, we see two different peaks in growth rate, first small peak in July 2020 and the second big peak in January 2021. For cluster 2, there is a rapid increase in growth rate since July 2021 that peaks on around December 2020. Unlike cluster 2, in cluster 3, the first big peak in July is followed by a small peak in January. Finally, in cluster 4, the growth rate rapidly increases from April to July followed by a rapid fall and a small rise in January 2021. Overall, the first peak in most of the states occurred in around July, 2020 and the second peak occurred in around January, 2021.

Figure 3

Average growth rate in each cluster, the dotted blue lines are the one standard deviation band around the average growth rate.

Discussion

In this study we investigated the excess mortality patterns during the COVID-19 pandemic at the national and subnational level in Mexico from March 1, 2020, to April 10, 2021. We estimated an excess all-cause mortality rate of 39.66 per 10 000 population at the national level (a total of ∼508 289 excess deaths), indicating a devastating mortality impact of the COVID-19 pandemic in Mexico. Mexico City alone accounted for about 19% of total excess deaths in Mexico, with an excess mortality rate of 106.17 per 10 000 population. We found that the excess mortality rate has continuously declined after the second COVID-19 peak during the week of January 17-23, 2021.

Interestingly, we found that the states with the highest excess death rate (i.e., Mexico City, Tlaxcala, Morelos, Puebla, Mexico) were in the central states in Mexico, while the lowest excess death rates were observed in the southern states (i.e., Chiapas, Oaxaca, Quintana Roo, Yucatan, Campeche). In Mexico, the majority of the indigenous population are located in the southern states. According to intercensal Survey of 2015, 75% of the country's indigenous population lived in 8 states (highest in Oaxaca (14.42%), followed by Chiapas (14.19%), Veracruz (9.16%), Mexico (9.13), Puebla (9.10%), Yucatan (8.75%), Guerrero (5.67%), and Hidalgo (5.04%) (National Institute of Indigenous Peoples, July 9, 2017). Similarly, states with the highest proportion of native population as a proportion of state population in 2015 were Yucatan (50.2%), Oaxaca (43.7%), Chiapas (32.7%), Quintana Roo (32.5%), and Campeche (22.2%) (National Institute of Indigenous Peoples, July 9, 2017.) Compared to non-native groups, the indigenous populations across continents have suffered significant health disparities and a greater burden of diseases, including higher infant mortality, and lower life expectancy (Curtice and Choo, 2020). During the pandemic, the indigenous populations have remained at higher risk of infection and death (CIDRAP, Power et al., 2020). A previous study demonstrated higher excess deaths in U.S. states with higher concentration of Native Americans during the 1918 influenza pandemic (Dahal et al., 2018b). Similarly, studies from New Zealand, Norway, and Alaska have also found that the indigenous populations in those regions were disproportionately affected by the 1918 influenza pandemic compared to the non-indigenous population (Mamelund, 2003, Mamelund et al., 2013, Rice, 2018). In contrast and very intriguingly, in this study we found a lower excess mortality rate during the pandemic in Mexican states with a higher proportion of the native indigenous population. This is an unexpected finding that warrants further inquiry and examination as it may provide great insight to factors that may potentially buffer against the impact of the pandemic and other adverse health events, including natural disasters, if the findings can be replicated in other studies.

In our analyses of data from Mexico, we found that COVID-19 specific deaths accounted for only 42.16% of total excess death at the national level, lowest in Chiapas (20.97%) and highest in Quintana Roo (76.05%). At the state level in Mexico, the timing and the rigor of implementation of public policies to contain the virus varied widely (Knaul et al., 2021). For example, some of the states, such as Veracruz, Yucatan, Nuevo Leon, and Tamaulipas, established policies to promote social distancing before the federal government enacted those policies (Knaul et al., 2021). While states such as Chiapas, Tabasco, San Luis Potosi, and Zacatecas underperformed in implementation of public policy measures (Knaul et al., 2021), some other states, with a relatively low excess mortality rate, such as Baja California Sur and Nayarit, implemented public information campaigns and international travel restrictions for longer periods, despite the potential adverse impact on tourism, which is a major economic activity (Knaul et al., 2021).

We found a positive association between the aging index and excess mortality in the adjusted model confirming previous studies linking older age and COVID-19. The aging index is defined as the number of older adults (60 years of age and older) for every 100 children and youth (0 to 14 years of age) (Instituto Nacional de Estadística y Geografía (INEGI)), and it increases as the population ages. Older age is a significant predictor of COVID-19 mortality as well as mortality from other causes (Bello-Chavolla et al., 2020, Ho et al., 2020, CDC, 2021). Previous studies that analyzed excess mortality patterns during the first wave of COVID-19 in 21 industrialized countries have shown that those aged 65 years and above comprised 94% of all excess deaths, indicating a very high risk of death among older aged population specifically due to COVID-19 (Kontis et al., 2020).

Similarly, our finding of a positive association between the marginalization index and excess mortality supports previous research and reviews underscoring the close link between social disadvantage and COVID mortality (Saini et al., 2021) and the overall increased burden of the pandemic in marginalized populations. The marginalization index that we used is an indicator of the inequities in quality of housing, access to basic public services like electricity and drinking water, schooling, proportion of poorly paid population and other sociodemographic and population health characteristics (CONAPO 2020, World Bank. World Development Report 2009). There may be several explanations for our findings. For example, public health measures such as social distancing and sheltering in place to combat the COVID-19 pandemic resulted in a disproportionate burden to vulnerable and marginalized populations (Anderson et al., 2020, Benfer and Wiley, March 19 2020, Kantamneni, 2020). Marginalized groups are also more likely to be infected by the coronavirus due to the context of their living arrangements, which may limit the ability to self-isolate and socially distance. Similarly, it is well demonstrated that marginalized populations tend to have a higher prevalence of chronic conditions such as obesity, hypertension and diabetes, which are all strong risk factors associated with poor prognostic outcomes among those infected with COVID-19 (Anderson et al., 2020). To complicate matters further, marginalized populations are also often at greater risk of dying due to other indirect causes such as limited access to already-stressed health care systems, poor mental health outcomes, food insecurity, lower health literacy and lower consumption of health services, abuse, and violence, among other social ills (Anderson et al., 2020, Benfer and Wiley, March 19 2020, Evans et al., 2021).

Interestingly, Chiapas, a southern Mexico state with a high marginalization index (CONAPO, 2020), a higher concentration of indigenous population (32.7% indigenous population as of 2015) (National Institute of Indigenous Peoples, July 9, 2017), and lower average performance in implementation of public policies to combat COVID-19 (Knaul et al., 2021), had the lowest excess mortality rate among the states examined in these analyses. While there is no clear explanation for these findings, the varying climate across the states examined may be a contributing factor. For example, the southern states in Mexico have a weak seasonality and a tropical climate throughout the year (Burton, 2013). According to previous studies conducted in Mexico (Méndez-Arriaga, 2020), tropical climate delayed the local transmission of SARS-CoV-2 at regional level. As such, the temperate climate regions like Tlaxcala and Jalisco may have been more vulnerable for local transmission than the tropical climate regions such as Chiapas and Veracruz (Méndez-Arriaga, 2020). However, it is interesting that in our analysis we did not find a significant difference in excess mortality rates across states that fall in three distinct climate groups in Mexico. These findings should be replicated in other settings that comprise multiple climate regions to determine the impact of seasonality in the transmission spread and impact of excess mortality patterns. Additionally, more research is needed to elucidate the factors associated with lower all-cause excess death rate in the relatively marginalized southern states as observed in this study. Such insight may provide mitigation strategies for other regions with higher impact.

We also found a positive association between average household size and excess death rate in the adjusted model. Although the links between average family size and excess death rate at the state level have not been reported elsewhere, the average family size could interact with other social determinants of health such as poverty, food insecurity, and lack of access to health care. However, since in our model we have controlled for marginalization index and aging index, household size could just be picking up a higher degree of exposure to COVID-19. Further studies are needed to understand the potential mechanism underlying this association and to more specifically consider family size as a potential population-level indicator of communities at risk for increased impact.

Our classification of excess deaths growth rate curves at the state level reflects four distinct categories of Mexican states. In all of the clusters the first peak of the excess deaths growth rate curve occurred in around July, 2020 which happened after the phased reopening of non-essential services in June, 2020 in Mexico. The reopening of the country coincided with an increase in both driving and walking trends, and the highest levels of COVID-19 deaths that remained at a high level during June and July 2020 (Tariq et al., 2021). The visual analysis of the growth rate curve indicates that the coastal tropical southeastern states were most affected during the first few months of the pandemic compared to other states. However, these states exhibited the lower overall excess death rate, which could indicate the effect of temperature and other environmental factors. Moreover, information on the growth rate curves can be utilized at the state level to guide the implementation of medical and public health measures. Besides, learning from the public health measures implemented in states of one cluster (for example, cluster 4) can be helpful to the other states (for example, states in cluster 1).

We observed that most of the states in clusters 1 and 3 were adjacent to each other compared with the states in clusters 2 and 4. Clusters 1 and 3 are the clusters with double peaks and are similar compared to clusters 2 and 4. These similarities could be due to shorter distances or greater population connectivity. We note that a previous study on the 1918 influenza pandemic found that nearby areas have similar excess mortality patterns. For example, Northern counties in Arizona, USA, had higher excess deaths compared to the Southern counties (Dahal et al., 2018b). It is worth noting that the overall wave pattern of mortality in Mexico (Figure 1) is comparable with the wave pattern displayed by neighboring US states including Texas, Arizona, and California (CDC, 2021). This wave pattern is also consistent with cluster 1, which comprises 12 states, including Baja California, and which borders with California and Coahuila bordering Texas. Cluster 1 also includes the central states such as Morelos, Mexico, and Mexico City, which share high air traffic connectivity with the US. Therefore, the overall wave pattern in Mexico could have been dominated by those states that are highly connected to the United States. Likewise, the variations in the multiple wave patterns across other states could also reflect how Mexico is integrated into a more global epidemiologic system.

To our knowledge, this is the first study that assesses the growth rate curves of excess deaths at the state level in Mexico. In our study, the estimates of excess death rate, as well as the proportion of COVID-19 attributed deaths, could be underestimated due to factors such as low testing rates in Mexico, misclassification of COVID-19 deaths, and delay in reporting COVID-19 deaths. In our analyses, we replaced negative excess mortality by zeros to account for underreporting of deaths which may not adjust all potential periods of negative excess mortality. Since our analyses are based on state-level data, ecological fallacy should be considered when interpreting the results at the individual level in a particular state. Additionally, other potential confounders that were not measured may explain the patterns of excess mortality across states. These limitations should be taken into account when interpreting our findings.

Conclusion

Our estimate of all-cause excess death rate in Mexico was 39.66 per 10 000, with central states exhibiting higher rates and southern states exhibiting lower rates. Our study highlights that several population measures including the aging index, marginalization index, and average household size were significantly associated with the all-cause excess mortality rates across Mexican states during the COVID-19 pandemic. Our excess mortality estimates can help tailor state specific medical and public health interventions to prevent excess mortality in vulnerable areas but targeting specific regions and socio-economic indicators. We also recommend further studies that investigate the lower excess death rate in southern states, and studies that explore the role of environmental factors, particularly the social determinants of health, in spatial variation in excess death rate in Mexico and other regions heavily impacted by COVID-19.