Spatial Survival Model for COVID-19 in México

A spatial survival analysis was performed to identify some of the factors that influence the survival of patients with COVID-19 in the states of Guerrero, México, and Chihuahua. The data that we analyzed correspond to the period from 28 February 2020 to 24 November 2021. A Cox proportional hazards frailty model and a Cox proportional hazards model were fitted. For both models, the estimation of the parameters was carried out using the Bayesian approach. According to the DIC, WAIC, and LPML criteria, the spatial model was better. The analysis showed that the spatial effect influences the survival times of patients with COVID-19. The spatial survival analysis also revealed that age, gender, and the presence of comorbidities, which vary between states, and the development of pneumonia increase the risk of death from COVID-19.


Introduction
The coronavirus disease COVID-19, which is caused by the SARS-CoV-2 virus (Severe Acute Respiratory Syndrome Coronavirus type 2), was first reported in Wuhan City, China, on 31 December 2019, and on 11 March 2020, COVID-19 was declared a pandemic by the World Health Organization (WHO) [1].The first person with COVID-19 in México was detected on February 27 of the same year [2].
México is a country with a high prevalence of comorbidities, such as hypertension, obesity, and diabetes.Given that it has been shown that the existence of comorbidities associated with SARS-CoV-2 infection increases the risk of mortality [3], México faced a severe crisis during the COVID-19 pandemic.
The risk factors of SARS-CoV-2 have been investigated using the survival analysis methodology [4].In México, previous studies have investigated the risk factors associated with deaths due to COVID-19 using survival analysis.It was reported that male gender; advanced age; the presence of comorbidities such as chronic kidney disease, diabetes, and arterial hypertension; the development of pneumonia; hospitalization; admission to an intensive care unit (ICU); intubation; geographic location; and the health sector where the patient was treated are associated with lower survival in patients affected by SARS-CoV-2 [5][6][7].Furthermore, [8] found that poorer population groups have lower COVID-19 survival rates.
Survival data are analyzed using proportional hazards regressions [9,10].Considering that data can come from different locations, regions, or zones and that the number of cases varies between locations, it is appropriate to include location information in a survival model.The variability of data in space can be included in a survival model by means of a random effect or term [11].
Statistical models that take geographic location into account are increasingly used in survival analysis [12].This is because they (i) incorporate spatial variation in survival times [13]; (ii) implement the Bayesian approach for parameter estimation and the new tools of Geographic Information Systems [13][14][15]; (iii) surpass the classical models by reducing the bias in the estimates and incorporating geographic information of the data [16,17]; (iv) introduce information about the spatial locations of the data, which plays an important role in predicting survival as it serves as a proxy for unmeasured regional characteristics such as access to medical care [13,18]; and (v) provide relevant information for public health decision making, such as the extent and direction of disease spread or the locations of disease hotspots, allowing control measures to be more effective and facilitating the adequate allocation of health resources [19].
Some papers using survival analysis in the health field that take the spatial information of data into account include Bayesian spatial survival modeling for dengue fever in Indonesia [20]; modeling the time to detection of urban tuberculosis in Portugal [21]; modeling of spatially correlated survival data for people with different types of cancer [22]; modeling of spatial variation in leukemia survival data [23]; parametric normal transformation models for spatially correlated, right-censored survival data [24]; Bayesian Weibull and Cox semi-parametric spatial models to describe a data set on dengue hospitalization [25]; an estimate of recovery times of COVID-19 patients in India [26]; a hierarchical conditional autoregressive model for colorectal cancer survival data [27]; a semi-parametric model with the cure fraction for multivariate time-to-event data [28]; a proportional hazards spatial frailty model to determine risk factors associated with under-five mortality in Kenya [29]; a hierarchical Bayesian model to jointly model longitudinal and survival data considering the effects of spatial and temporal frailty for AIDS data [30]; and a proportional hazards models with marginal cure rates to identify risk factors for tooth loss and predict the remaining useful life of a patient's teeth [31].
Previous studies conducted in México did not consider the spatial effect on the survival of patients with COVID-19.Therefore, this paper aims to investigate the risk factors associated with COVID-19 deaths in three Mexican states, considering the spatial variability of the data in the statistical model.

Study Area
México is a country located in North America; it has 32 Federal Entities, also known as states.Each state is made up of a certain number of municipalities; there are a total of 2475 municipalities throughout the country.The state with the most municipalities is Oaxaca with 570, followed by Puebla with 217, Veracruz with 212, and Jalisco and the State of México with 125 municipalities each.The states with the fewest municipalities are Baja California with five and Baja California Sur with five (Figure 1).In the year 2020, México's total population was 126,014,024 people [32].The states with the highest populations were the State of México with 16,992,418 people, México City with 9,209,944 people, Veracruz with 8,062,579 people, Jalisco with 8,348,151 people, and Puebla with 6,583,278 people, while the states with the lowest populations were Colima with 731,391 people and Baja California Sur with 798,447 people (Figure 1).
According to the 2020 Population and Housing Census conducted by the National Institute of Statistics and Geography (INEGI), the population density in México was 64 inhabitants per square kilometer (hab/km 2 ) [32].The INEGI categorizes the 32 Federal Entities as high, medium, or low population density (Figure 1).
The states with high population densities are Colima, Aguascalientes, Guanajuato, Querétaro, Hidalgo, the State of México, Morelos, México City, Puebla, and Tlaxcala; the states with medium population densities are Baja California, Sinaloa, Nuevo León, Jalisco, Michoacán de Ocampo, Guerrero, Veracruz de Ignacio de la Llave (Veracruz), Tabasco, Chiapas, and Yucatán; and the states with low population densities are Baja California Sur, Campeche, Chihuahua, Coahuila, Durango, Nayarit, Oaxaca, Quintana Roo, San Luis Potosí, Sonora, Tamaulipas, and Zacatecas (Figure 1).Sur, Campeche, Chihuahua, Coahuila, Durango, Nayarit, Oaxaca, Quintana Roo, San Luis Potosí, Sonora, Tamaulipas, and Zacatecas (Figure 1).To analyze the survival times of people infected with SARS-CoV-2 in space, we worked with the State of México (red polygon in Figure 2) with high population density, 125 municipalities, and a population of 16,992,418; Guerrero (orange polygon in Figure 2) with medium population density, 81 municipalities, and a population of 3,540,685; and Chihuahua (yellow polygon in Figure 2), with low population density, 67 municipalities, and a population of 3,741,869.
Statistical analyses were performed in R software version 4.1.2[33] (The R Foundation for Statistical Computing, Vienna, Austria).To analyze the survival times of people infected with SARS-CoV-2 in space, we worked with the State of México (red polygon in Figure 2) with high population density, 125 municipalities, and a population of 16,992,418; Guerrero (orange polygon in Figure 2) with medium population density, 81 municipalities, and a population of 3,540,685; and Chihuahua (yellow polygon in Figure 2), with low population density, 67 municipalities, and a population of 3,741,869.
Statistical analyses were performed in R software version 4.1.2[33] (The R Foundation for Statistical Computing, Vienna, Austria).

Database
This study used the open-access database of suspected cases of COVID-19 published by the Ministry of Health of México through the Epidemiological Surveillance System for Respiratory Viral Diseases [34].In México, as of 24 November 2021, a cumulative total of 3,872,263 cases had been confirmed.However, the open-access database contains 2,028,000 records of confirmed COVID-19 patients reported by the laboratory of the National Network of Epidemiological Surveillance Laboratories and private laboratories endorsed by the Institute of Epidemiological Diagnosis and Reference [33].The analyzed data correspond to the period from 28 February 2020, to 24 November 2021.

Database
This study used the open-access database of suspected cases of COVID-19 published by the Ministry of Health of México through the Epidemiological Surveillance System for Respiratory Viral Diseases [34].In México, as of 24 November 2021, a cumulative total of 3,872,263 cases had been confirmed.However, the open-access database contains 2,028,000 records of confirmed COVID-19 patients reported by the laboratory of the National Network of Epidemiological Surveillance Laboratories and private laboratories endorsed by the Institute of Epidemiological Diagnosis and Reference [33].The analyzed data correspond to the period from 28 February 2020, to 24 November 2021.
The variables considered in this study were age, gender, pneumonia, diabetes, chronic obstructive pulmonary disease (COPD), cardiovascular disease (CVD), obesity, asthma, chronic kidney disease (CKD), and hypertension (Table 1).The patients' sociodemographic variables were not considered because the Ministry of Health, which collected the information, did not include them and there was no way to obtain them for the patients considered in this study.The variables considered in this study were age, gender, pneumonia, diabetes, chronic obstructive pulmonary disease (COPD), cardiovascular disease (CVD), obesity, asthma, chronic kidney disease (CKD), and hypertension (Table 1).The patients' sociodemographic variables were not considered because the Ministry of Health, which collected the information, did not include them and there was no way to obtain them for the patients considered in this study.
The response variable was survival time, which was defined as the time between the date of symptom onset and patient death.The data were censored on 24 November 2021, for people who were alive at the end of the study period.

Spatial Autocorrelation
To study the survival times in the three states of México, we first determined, using the Moran index, the degree of spatial association that existed in the confirmed cases of COVID-19.Moran's index (I G ) is defined as where Y i and Y j are the values of the variable Y in localities i and j, respectively; Ȳ is the average of the variable Y in the m localities under study; and w ij indicates the proximity between units i and j.When locations i and j are neighbors, w ij = 1.Otherwise, w ij = 0.
For i = j, w ij = 1.The hypotheses set for I G are given as follows: In this work, the localities are the municipalities of each of the three states.Variable Y is the number of confirmed cases of COVID-19.Two municipalities are neighbors if they share at least one point on the border, which can be a vertex and/or at least one border (Queen).Thus, the distance function is given by an indicator function, where w ij = 1 if the municipalities share a vertex and/or at least one border.Also, w ij = 1 if i = j.Otherwise, w ij = 0 and i, j = 1, . . ., m.With the w ij values, we constructed the weight matrix (W) of dimension m × m, assuming that ∑ m i=1 ∑ m j=1 = m.For the State of México, m = 125; for Guerrero, m = 81; and for Chihuahua, m = 67.The calculation of the Moran index was performed with the moran.mcfunction of the spded package in R statistical software version 4.1.2(The R Foundation for Statistical Computing, Vienna, Austria) [33].

Statistical Model
Very often, time-to-event data are grouped into strata (clusters), such as clinical sites, geographic regions, and so on.t ij is the time to death or censoring for the j-th subject at location s i ; i = 1, . . ., m; j = 1, . . ., n i ; and the total number of subjects under study is n = ∑ m i=1 n i .x ij is a vector of covariates of dimension p, and β = β 1 , . . ., β p ′ is a vector of regression coefficients.The proportional hazard (PH) model (Model ( 2)) is shown below: where h 0 is the baseline hazard.
The proportional hazard frailty model (Model ( 3)) is shown below: where v i is an unobserved frailty associated with s i and is designed to capture differences among the strata.In the spatial survival analysis, subjects are in m distinct regions (s 1 , . . ., s m ).
The PH frailty model has a survival function: with a density function: where S 0 (•) is the baseline survival function and f 0 (•) is the baseline density function, which are assumed to be unique for all individuals.For this study, The likelihood of a PH frailty model is given by

Prior Distributions
The parameters of the PH model have the following prior distributions: where TBP L refers to the transformed Bernstein polynomial [35].For a fixed positive integer (L), the prior TBP L (α, S θ (•)) is defined as where w L = (w 1 , . . . ,w L ) ′ is a vector of positive weights, I(• | a, b) is a beta cumulative distribution function with parameters (a, b), and {S θ (•) : θ ∈ Θ} is a parametric family of survival functions with support on positive reals.Furthermore, a log-logistic distribution for S θ (t) is assumed; that is, where measures the amount of spatial variation across locations, and the (i, j) element of R is modeled as which is a correlation function controlling the spatial dependence of v(s).ϕ > 0 is a range parameter controlling the spatial decay over distance, k ∈ (0, 2] is a shape parameter, and ∥ s i − s j ∥ is the distance between s i and s j .The prior GRF τ 2 , ϕ is defined as where r ij is the (i, j) element of R −1 [35].
The spatial dependence of COVID-19 survival times is captured in the covariance structure (R) of the Gaussian random field v, which is assumed to be a stationary process.

Posterior Distributions
For Equation (2), the likelihood function is Therefore, the posterior distribution is For Equation (3), the likelihood function is Therefore, the posterior distribution is A Bayesian fitting of the proportional hazard frailty model was obtained using R software and several libraries, such as spBayesSurv [36].The function survregbayes set the following hyperparameters as defaults: To analyze the survival times of people with COVID-19, the proportional hazard frailty model was used (Model (3)).Right censoring was used, and the model's covariates were age, gender, pneumonia, diabetes, chronic obstructive pulmonary disease, cardiovascular disease, obesity, asthma, chronic kidney disease, and hypertension.For Guerrero, the analysis was carried out with the 37,278 registered cases of COVID-19; for Chihuahua, we worked with the 45,954 observed cases of COVID-19; and for the State of México, a sample of 100,000 observations out of the 176,268 available observations were used.The initial values of the hyperparameters of the prior distributions (Equations ( 5) and ( 6)) were set at a 0 = 4 and b 0 = 4 for Chihuahua and the State of México and at a 0 = 1 and b 0 = 1 for Guerrero.For the three states, the hyperparameters were fixed at L = 15, a τ = 5, b τ = 5, a ϕ = 6, and b ϕ = 3.For the estimation of the parameters of Model 3, Markov Chain Monte Carlo (MCMC) algorithms were used, which consisted of two chains of 25,000 iterations, with a thinning of 10 and a burning of 5000 samples.
The MCMC procedure is carried out with the Gelman and Rubin convergence diagnostic ( R) [37,38], which is implemented in the library coda [39].Values of R close to one give evidence of the convergence of the MCMC [37].Another diagnostic is the MCMC traceplot (plot of the values in the simulated chains vs. the iteration).Good mixing of the chains indicates convergence of the MCMC.Cox-Snell residuals [35] can be used to verify the proportional hazards assumption [40].However, testing this assumption numerically or graphically is complicated since the proportional hazards hypothesis only approximates the correct model for one covariate and any formal test based on a sufficiently large sample will reject the null hypothesis of proportionality [41].In any case, Cox-Snell residuals should be considered to determine the model fit [42].
A classical survival analysis was also performed using the Cox model (Model ( 2)).For the three states, a 0 = 6, b 0 = 6, L = 15, a τ = 5, b τ = 5, a ϕ = 6, and b ϕ = 3 were considered.However, only the results of the best model are reported, which for the three states, was the proportional hazard frailty model.We used three criteria to compare the fitted models: the deviance information criteria (DIC) [43], the Watanabe-Akaike information criterion (WAIC) [44], and the log pseudo marginal likelihood (LPLM) [45].Generally, smaller DIC and WAIC values show good model fitting, while larger LPML values indicate a better predictive performance of a model.

State of México
As of 24 November 2021, 176,268 confirmed cases of COVID-19 were registered in the State of México.The municipalities with more than 10,000 cases were located in the east and center of the state: Ecatepec de Morelos (22,613 cases), Nezahualcóyotl (16,992 cases), Toluca (13,647 cases), Naucalpan de Juárez (12,450), and Tlalneplantla de Baz (10,605 cases).The municipalities with less than 20 cases were located in the southeast of the state, including Zacazonapan (18 cases), Ixtapan del Oro (14 cases), and Otzoloapan (12 cases) (Figure 3).
According to the Moran index, in the State of México, confirmed COVID-19 cases were more correlated in the municipalities that are neighbors than in the municipalities that are not (I G = 0.363, p − value < 0.05).Therefore, it was convenient to use the proportional hazard frailty model to model COVID-19 survival times.
The estimations of the posterior means and medians of the parameters, as well as their 95% credible intervals (CrIs), and the Gelman and Rubin convergence diagnostic values for the spatial model (proportional hazard frailty model) of the State of México are presented in Table 2.Because the values of R are close to 1, there is convergence in the model chains [37].The chains show good mixing in Figure A1 (Appendix A).However, there is no indication that the Cox assumption is violated since no large deviations from the Cox-Snell residuals are observed in Figure A2 (Appendix A).
Healthcare 2024, 12, x FOR PEER REVIEW 8 of 21 cases), Toluca (13,647 cases), Naucalpan de Juárez (12,450), and Tlalneplantla de Baz (10,605 cases).The municipalities with less than 20 cases were located in the southeast of the state, including Zacazonapan (18 cases), Ixtapan del Oro (14 cases), and Otzoloapan (12 cases) (Figure 3).  2. Because the values of  are close to 1, there is convergence in the model chains [37].The chains show good mixing in Figure A1 (Appendix A).However, there is no indication that the Cox assumption is violated since no large deviations from the Cox-Snell residuals are observed in Figure A2 (Appendix A).
The estimates of TBP L and GRF were also significant.The variance of the spatial term v, which was given by τ2 = 1.540, was significant, with a 95% CrI of 0.936-2.668(Table 2).Therefore, the inclusion of the random effect was relevant, which means that the risk of death from COVID-19 was not homogeneous across the 125 municipalities of the State of México.
Figure 4 shows the average values of posterior sampling for frailties (v i ).We can identify clusters of municipalities.The inhabitants of the eastern municipalities had the highest risk of mortality adjusted for the effects of covariates.Figure 4 shows the average values of posterior sampling for frailties ( ).We can identify clusters of municipalities.The inhabitants of the eastern municipalities had the highest risk of mortality adjusted for the effects of covariates.

State of Guerrero
As of 24 November 2021, 37,278 confirmed cases of COVID-19 were registered in Guerrero.The municipalities with more than 1000 cases were located in the south, center, and north of the state: Acapulco de Juárez (14,639 cases), Chilpancingo de los Bravo (6773 cases), Zihuatanejo de Azueta (1984 cases), Iguala de la Independencia (1664 cases), and Taxco de Alarcon (1007 cases).The municipalities with the fewest cases were located in the east and north of the state: General Canuto A. Neri (nine cases), Atlamajalcingo del Monte (eight cases), Tlacoapa (six cases), Pedro Ascencio Alquisiras (four cases), and Iliatenco (four cases) (Figure 5).
According to the Moran index, in the state of Guerrero, the confirmed COVID-19 cases were more correlated in the municipalities that are neighbors than in the municipalities that are not (I G = 0.135, p − value < 0.05).
Estimates of the means and posterior medians of the parameters, their corresponding 95% CrIs, and the values of R (which indicate convergence) are presented in Table 3. Figure A3 (Appendix A) reveals that the MCMC chains of the estimated parameters show good mixing.On the other hand, there is no indication that the Cox assumption is violated since no large deviations from the Cox-Snell residuals are observed in Figure A4  (Appendix A). and north of the state: Acapulco de Juárez (14,639 cases), Chilpancingo de los Bravo (6773 cases), Zihuatanejo de Azueta (1984 cases), Iguala de la Independencia (1664 cases), and Taxco de Alarcon (1007 cases).The municipalities with the fewest cases were located in the east and north of the state: General Canuto A. Neri (nine cases), Atlamajalcingo del Monte (eight cases), Tlacoapa (six cases), Pedro Ascencio Alquisiras (four cases), and Iliatenco (four cases) (Figure 5).According to the Moran index, in the state of Guerrero, the confirmed COVID-19 cases were more correlated in the municipalities that are neighbors than in the municipalities that are not ( 0.135,  value 0.05).Estimates of the means and posterior medians of the parameters, their corresponding 95% CrIs, and the values of  (which indicate convergence) are presented in Table 3. Figure A3 (Appendix A) reveals that the MCMC chains of the estimated parameters show good mixing.On the other hand, there is no indication that the Cox assumption is violated since no large deviations from the Cox-Snell residuals are observed in Figure A4 (Appendix A).
Asthma was not significant, while male gender and age were factors that increased the risk of death in patients with COVID-19.
Based on Table 3, the comorbidities associated with an increased risk of death from COVID- 19    Asthma was not significant, while male gender and age were factors that increased the risk of death in patients with COVID-19.
The estimates of TBP L and GRF were also significant.The posterior mean of the variance ( τ2 = 1.753) of the spatial PH model was significant, with a 95% CrI of 1.015-3.255(Table 3).Therefore, the random effects had an influence.This means that the risk of COVID-19 was not homogeneous across the 81 municipalities of the state of Guerrero.
Figure 6 shows the quantiles of the average values of the posterior samples of the spatial term v.We can identify clusters of municipalities.People with COVID-19 in the southern municipalities of Guerrero and, in general, those located on the periphery of the state presented the highest mortality risk adjusted for covariate effects.
The estimates of  and  were also significant.The posterior mean of the variance (̂ 1.753) of the spatial PH model was significant, with a 95% CrI of 1.015-3.255(Table 3).Therefore, the random effects had an influence.This means that the risk of COVID-19 was not homogeneous across the 81 municipalities of the state of Guerrero.
Figure 6 shows the quantiles of the average values of the posterior samples of the spatial term .We can identify clusters of municipalities.People with COVID-19 in the southern municipalities of Guerrero and, in general, those located on the periphery of the state presented the highest mortality risk adjusted for covariate effects.

State of Chihuahua
As of 24 November 2021, 45,954 confirmed COVID-19 cases were registered in the state of Chihuahua.The municipalities with more than 1000 cases were located in the central and eastern parts of the state: Chihuahua (12,582 cases), Delicias (1842 cases), Hidalgo del Parral (1750 cases), and Cuauhtémoc (1687 cases).The municipalities with less than five cases were located in the west of the state: Matachí (three cases), Huejotitlán (two cases), and Maguarichi (one case) (Figure 7).
Although the confirmed cases of COVID-19 in the state of Chihuahua presented a random pattern (I G = −0.0423,p − value = 0.72), we used the proportional hazard frailty model to study the COVID-19 survival times.
Table 4 presents the means and medians of the estimated parameters with their respective 95% CrIs and R values.The Gelman and Rubin statistic ( R) has values close to 1, which implies that there is convergence in the Markov chains.Figure A5 (Appendix A) shows a good mix of the chains of the parameters of the PH frailty model.On the other hand, there is no indication that the Cox assumption is violated since no large deviations from the Cox-Snell residuals are observed in Figure A6 (Appendix A).
As in the State of México and Guerrero, male gender and age were associated with a higher risk of death.The variables that were not significant were COPD, CVD, and asthma.
The variance of v ( τ2 = 2.175) was significant, with a 95% CrI of 1.221-3.706 (Table 4).The random effect had an influence, which means that the risk of COVID-19 was not homogeneous across the 67 municipalities of the state of Chihuahua.In general, the municipalities with the most cases of COVID-19 presented larger values of v.This was the case for the municipalities in the north and southeast of the state of Chihuahua, which had the highest risk of mortality adjusted for the effects of covariates (Figure 8).
As of 24 November 2021, 45,954 confirmed COVID-19 cases were registered in the state of Chihuahua.The municipalities with more than 1000 cases were located in the central and eastern parts of the state: Chihuahua (12,582 cases), Delicias (1842 cases), Hidalgo del Parral (1750 cases), and Cuauhtémoc (1687 cases).The municipalities with less than five cases were located in the west of the state: Matachí (three cases), Huejotitlán (two cases), and Maguarichi (one case) (Figure 7).  4 presents the means and medians of the estimated parameters with their respective 95% CrIs and  values.The Gelman and Rubin statistic ( ) has values close to 1, which implies that there is convergence in the Markov chains.Figure A5 (Appendix A) shows a good mix of the chains of the parameters of the PH frailty model.On the other hand, there is no indication that the Cox assumption is violated since no large deviations from the Cox-Snell residuals are observed in Figure A6 (Appendix A).
As in the State of México and Guerrero, male gender and age were associated with a higher risk of death.The variables that were not significant were COPD, CVD, and asthma.Finally, according to the DIC and WAIC, models that incorporate the spatial effect have better fits than models that do not (Table 5).
For the LPML, the three spatial models also have better predictive fits than the nonspatial models.Therefore, the relevance of using the proportional hazard frailty models to model COVID-19 survival times is verified.Finally, according to the DIC and WAIC, models that incorporate the spatial effect have better fits than models that do not (Table 5).
For the LPML, the three spatial models also have better predictive fits than the nonspatial models.Therefore, the relevance of using the proportional hazard frailty models to model COVID-19 survival times is verified.

Discussion
In this study, we present a different approach to modeling the spatial dependence of the survival times of patients with COVID-19 in three Mexican states.Spatial heterogeneity, sociodemographic variables, and individual characteristics, such as comorbidities, were found to be associated with the severity of COVID-19, as reported in several investigations.
According to the DIC and WAIC, models incorporating a spatial term have a better predictive fit than Cox models without a frailty term.This is because classical survival models cannot account for the spatial correlation of the data and, therefore, should not be con-

Discussion
In this study, we present a different approach to modeling the spatial dependence of the survival times of patients with COVID-19 in three Mexican states.Spatial heterogeneity, sociodemographic variables, and individual characteristics, such as comorbidities, were found to be associated with the severity of COVID-19, as reported in several investigations.
According to the DIC and WAIC, models incorporating a spatial term have a better predictive fit than Cox models without a frailty term.This is because classical survival models cannot account for the spatial correlation of the data and, therefore, should not be considered as the standard model in research when geographic information is available [12,16].These results are in agreement with those reported by Thamrin et al. [20], who compared a Bayesian spatial survival model with a nonspatial one to analyze the factors influencing the survival of dengue patients and found that the spatial model was more appropriate.Daniel et al. [29] found that a spatial Cox proportional hazards model performed better compared with a nonspatial model in identifying risk factors for under-five mortality.In a study by Mahanta et al. [26] related to recovery times in patients with COVID-19, a spatial survival model presented a better fit than a model without a frailty term.
According to the obtained results, men have a higher risk of dying than women, age is associated with a higher risk, and certain comorbidities can increase this risk, which is consistent with several studies reported in the literature [5][6][7][8]47].
Diabetes, hypertension, obesity, and CKD were the only statistically significant comorbidities in the three spatial survival models.The impact of comorbidities on the risk of death from COVID-19 is well known; however, the Mexican population has a high prevalence of metabolic diseases, which makes it vulnerable to developing complications from COVID-19.Mexico is second worldwide in the prevalence of obesity; according to the ENSANUT 2018 survey, 75.2% of the Mexican population over 20 years of age is overweight or obese.In addition, the prevalence of diabetes in Mexicans over 20 years of age is 10.3%, and the prevalence of hypertension is 18.4% in patients over 20 years of age [48].For Chihuahua and the State of Mexico, CKD was associated with a higher probability of dying.According to other studies on comorbidities recorded in Mexican data sets, CKD posed the highest risk of severe COVID-19 in Mexico [49].In contrast, in the state of Guerrero, the highest risk was posed by diabetes.
It is worth mentioning that for the State of México, people who had asthma had a lower risk of dying compared with people who did not have asthma.Previous studies have reported this "protective" effect [50,51].Some authors report that asthma may protect against the fatal outcomes of COVID-19 due to several possible mechanisms, such as the use of inhaled corticosteroids, chronic inflammation, reduced viral exposure due to protection, and/or mucus hypersecretion [51,52].
People with pneumonia in the state of Guerrero had a higher risk of dying compared with people in the states of Chihuahua and México (HR = 16.776,95% CrI: 15.150-18.411).A probable cause of this high risk of death is that Guerrero is one of the states with the worst health services in México, both in the perception of its population and due to a lack of supplies and medicines, and the lack of medical coverage is particularly severe in rural areas, which could have led to poor care for people who developed severe pneumonia.In fact, living in the southern region of the country is related to the severity of COVID-19.These disparities based on geographic location and ethnicity are closely linked to socioeconomic inequality: the southern region has higher rates of poverty and a concentration of indigenous people, which shows how different forms of inequality intersect [53].
There is evidence of geographic disparities in COVID-19 survival times in the three states that were analyzed, which may be influenced by socioeconomic factors, demographic factors, and health-related lifestyle factors.The impacts of these comorbidities vary by geographic location, and some are more important predictors of the risk of COVID-19 death in some states than in others.This demonstrates the importance of using global and local models, in this case, at the state level, to investigate the determinants of geographic disparities in health outcomes and health services utilization [54,55].
To our knowledge, this is the first study in México to use the spatial survival analysis approach to study risk factors associated with COVID-19 deaths that exhibit spatial variability across these three states.
The limitations of this work are as follows: First, the data set did not include clinical variables related to the evolution of the diseases due to COVID-19, which could have been useful for adjusting the model.Secondly, the variables that allowed us to identify the presence of comorbidities were self-reported; therefore, misclassification bias is very likely.Thirdly, the sociodemographic variables of the patients, such as age and income deciles, or access to basic services, such as drinking water and drainage, among others, were not recorded in the database.In future work, we intend to consider information related to socioeconomic aspects at the municipal level, such as the level of social progress, the level of urbanization, the level of poverty, and climatic conditions.In addition, this study will be extended to more states in the country.

Conclusions
In this paper, the survival times of patients with COVID-19 in three states of the Mexican Republic were studied, taking into account the geographical locations of the cases, comorbidities, age, gender, and the development of pneumonia.Patient survival times differed by geographic location.Therefore, patient location is an important factor for COVID-19 survival times.The proportional hazard frailty model performed better than the Cox model.According to the results, the COVID-19 survival times of male patients were shorter compared with women; age also influenced the survival times of patients with COVID-19.Obesity, diabetes, CKD, and hypertension were comorbidities that increased the risk of death from COVID-19 in all three states.In the State of México, patients with asthma had a lower risk of dying than patients without asthma.

Figure 1 .
Figure 1.States of México with number of municipalities (bars and left axis) and total population (lines and right axis).

Figure 1 .
Figure 1.States of México with number of municipalities (bars and left axis) and total population (lines and right axis).

Figure 2 .
Figure 2. States of México selected for this study.

Figure 2 .
Figure 2. States of México selected for this study.

Figure 3 .
Figure 3. Confirmed cases of COVID-19 in the municipalities of the State of México from 28 February 2020 to 24 November 2021.According to the Moran index, in the State of México, confirmed COVID-19 cases were more correlated in the municipalities that are neighbors than in the municipalities that are not ( 0.363,  value 0.05).Therefore, it was convenient to use the proportional hazard frailty model to model COVID-19 survival times.The estimations of the posterior means and medians of the parameters, as well as their 95% credible intervals (CrIs), and the Gelman and Rubin convergence diagnostic values for the spatial model (proportional hazard frailty model) of the State of México are presented in Table2.Because the values of  are close to 1, there is convergence in the model chains[37].The chains show good mixing in FigureA1 (Appendix A).However, there is no indication that the Cox assumption is violated since no large deviations from the Cox-Snell residuals are observed in FigureA2 (Appendix A).Male gender and age were factors that increased the risk of death from COVID-19, while the COPD variable was not significant.According to Table2and taking into account the hazard ratios (HRs) of the regression coefficients, CKD (HR = 1.588, 95% CrI: 1.482-1.707),diabetes (HR = 1.178, 95% CrI: 1.135-1.217),obesity (HR = 1.137, 95% CrI: 1.088-1.1852),and hypertension (HR = 1.065, 95% CrI: 1.026-1.102)were the comorbidities identified as risk factors.People who developed pneumonia had an increased risk of death (HR = 5.365, 95% CrI: 5.181-5.578).The estimated risk of death was 20% lower for people with asthma (HR = 0.796, 95% CrI: 0.677-0.922).The estimates of  and  were also significant.The variance of the spatial term , which was given by ̂ 1.540, was significant, with a 95% CrI of 0.936-2.668

Figure 3 .
Figure 3. Confirmed cases of COVID-19 in the municipalities of the State of México from 28 February 2020 to 24 November 2021.

Figure 4 .
Figure 4. Quantiles of posterior samples for frailties in the State of México.

Figure 4 .
Figure 4. Quantiles of posterior samples for frailties in the State of México.

Figure 5 .
Figure 5. Confirmed cases of COVID-19 in the municipalities of Guerrero from 28 February 2020 to 24 November 2021.

Figure 5 .
Figure 5. Confirmed cases of COVID-19 in the municipalities of Guerrero from 28 February 2020 to 24 November 2021.

Figure 6 .
Figure 6.Quantiles of posterior samples for frailties in Guerrero.

Figure 6 .
Figure 6.Quantiles of posterior samples for frailties in Guerrero.

Figure 7 .
Figure 7. Confirmed cases of COVID-19 in the municipalities of Chihuahua from 28 February 2020 to 24 November 2021.Although the confirmed cases of COVID-19 in the state of Chihuahua presented a random pattern ( 0.0423,  value 0.72), we used the proportional hazard frailty model to study the COVID-19 survival times.Table4presents the means and medians of the estimated parameters with their respective 95% CrIs and  values.The Gelman and Rubin statistic ( ) has values close to 1, which implies that there is convergence in the Markov chains.FigureA5 (Appendix A) shows a good mix of the chains of the parameters of the PH frailty model.On the other hand, there is no indication that the Cox assumption is violated since no large deviations from the Cox-Snell residuals are observed in FigureA6 (Appendix A).As in the State of México and Guerrero, male gender and age were associated with a higher risk of death.The variables that were not significant were COPD, CVD, and asthma.

Figure 7 .
Figure 7. Confirmed cases of COVID-19 in the municipalities of Chihuahua from 28 February 2020 to 24 November 2021.

Figure 8 .
Figure 8. Quantiles of posterior samples for frailties in Chihuahua.

Figure 8 .
Figure 8. Quantiles of posterior samples for frailties in Chihuahua.

Figure A1 .
Figure A1.Traceplots of the parameters in the proportional hazard spatial model for the State of México.: Transformed Bernstein polynomial parameter (Equation (10));  : spatial variation across locations of GRF; : range parameter.

Figure A2 .
Figure A2.Cox-Snell residuals for the PH frailty model of the State of México.

Figure A2 . 21 Figure A3 .
Figure A2.Cox-Snell residuals for the PH frailty model of the State of México.Healthcare 2024, 12, x FOR PEER REVIEW 18 of 21

Figure A3 .
Figure A3.Traceplots of the parameters in the proportional hazard spatial model for the state of Guerrero.α: Transformed Bernstein polynomial parameter (Equation (10)); τ 2 : spatial variation across locations of GRF; φ: range parameter.

Figure A3 .
Figure A3.Traceplots of the parameters in the proportional hazard spatial model for the state of Guerrero.: Transformed Bernstein polynomial parameter (Equation (10));  : spatial variation across locations of GRF; : range parameter.

Figure A4 .
Figure A4.Cox-Snell residuals for the PH frailty model of Guerrero.

Figure A5 .
Figure A5.Traceplots of the parameters in the proportional hazard spatial model for the state of Chihuahua.α: Transformed Bernstein polynomial parameter (Equation (10)); τ 2 : spatial variation across locations of GRF; φ: range parameter.

Figure A5 .
Figure A5.Traceplots of the parameters in the proportional hazard spatial model for the state of Chihuahua. : Transformed Bernstein polynomial parameter (Equation (10));  : spatial variation across locations of GRF; : range parameter.

Figure A6 .
Figure A6.Cox-Snell residuals for the PH frailty model of Chihuahua.

Figure A6 .
Figure A6.Cox-Snell residuals for the PH frailty model of Chihuahua.

Table 1 .
Variables considered in this study.

Table 1 .
Variables considered in this study.

Table 2 .
Estimated parameters of the spatial PH model for the State of México.

Table 3 .
Estimated parameters of the spatial PH model for the state of Guerrero.
CrI: Credible interval, R: Gelman and Rubin's convergence diagnostic, CrI-Upper: Upper limit of the credible interval, ** 95% credible interval does not include zero.

Table 4 .
Estimated parameters of the spatial PH model for the state of Chihuahua.