key: cord-0709210-32c2cbtd
authors: Middya, Asif Iqbal; Roy, Sarbani
title: Geographically varying relationships of COVID-19 mortality with different factors in India
date: 2021-04-12
journal: Sci Rep
DOI: 10.1038/s41598-021-86987-5
sha: 73a318faa9c8418025122a702b4253df68495ed4
doc_id: 709210
cord_uid: 32c2cbtd

COVID-19 is a global crisis where India is going to be one of the most heavily affected countries. The variability in the distribution of COVID-19-related health outcomes might be related to many underlying variables, including demographic, socioeconomic, or environmental pollution related factors. The global and local models can be utilized to explore such relations. In this study, ordinary least square (global) and geographically weighted regression (local) methods are employed to explore the geographical relationships between COVID-19 deaths and different driving factors. It is also investigated whether geographical heterogeneity exists in the relationships. More specifically, in this paper, the geographical pattern of COVID-19 deaths and its relationships with different potential driving factors in India are investigated and analysed. Here, better knowledge and insights into geographical targeting of intervention against the COVID-19 pandemic can be generated by investigating the heterogeneity of spatial relationships. The results show that the local method (geographically weighted regression) generates better performance ([Formula: see text] ) with smaller Akaike Information Criterion (AICc [Formula: see text] ) as compared to the global method (ordinary least square). The GWR method also comes up with lower spatial autocorrelation (Moran’s [Formula: see text] and [Formula: see text] ) in the residuals. It is found that more than 86% of local [Formula: see text] values are larger than 0.60 and almost 68% of [Formula: see text] values are within the range 0.80–0.97. Moreover, some interesting local variations in the relationships are also found.

India, no comprehensive study is performed at the local level to investigate geographical relationships between COVID-19 deaths and associated potential factors. To bridge the gap, a local method (GWR) is employed to explore the geographical distribution and associated potential socio-economic, demographic, and environmental factors for COVID-19 deaths. Note that, the GWR model helps us to identify whether there is geographical heterogeneity present in the relationships. Moreover, a comparison between local (OLS) and global (GWR) models are also performed. This paper offers further knowledge and insight into geographical targeting of intervention and control strategies against the COVID-19 epidemic. In summary, the key objectives of this study are (i) to explore the potential socio-economic, demographic, and environmental driving factors for COVID-19 deaths in India; (ii) to investigate geographically varying relationships of COVID-19 deaths with the driving factors by employing local (GWR) model. (iii) comparing the results of the local (GWR) model with the global (OLS) model to validate its suitability.

Data description. The geographical variabilities of COVID-19 deaths are modeled based on the districtlevel data across India. Note that, the COVID-19 mortality data are acquired for more than 400 districts in India. The geographical distributions of COVID-19 deaths are shown in Fig. 1 . The largest number of COVID-19 deaths are observed in the districts of the state Maharashtra. A total of 9 among 28 states contains at least one district that reports more than 1000 COVID-19 deaths. Table 1 summarizes all the raw datasets, their descriptions, the sources including the links from where these data can be found, and potential factors (independent variables) that are extracted from the raw datasets.

Datasets. Three raw datasets are mainly utilized to investigate geographically varying relationships of COVID-19 deaths with different environmental, demographic, and socio-economic factors. The first dataset includes district wise COVID-19 death counts in India. The cumulative number of COVID-19-related deaths for each district is collected up to February 24, 2021, from the COVID19INDIA website (https:// www. covid 19ind ia. org/). COVID19INDIA is a crowdsourced initiative to document the COVID-19 data from the states and union territories of India. In this study, the district-level COVID-19 death count is considered as the dependent variable. The second dataset pertaining to environmental pollution includes the daily concentration of different air pollutants (e.g. PM 2.5 , SO 2 , NO 2 , etc.). The concentration of air pollutants (from January 2016 to January 2020) for a total of 130 monitoring stations are obtained from the Central Pollution Control Board (CPCB 15 ), INDIA. The third dataset contains socio-economic and demographic data that may have an association with COVID-19 mortality. The district-level socio-economic and demographic data are obtained from the last census in India that was conducted in 2011. www.nature.com/scientificreports/ Additionally, the district-level data of each district needs to be linked with the GPS coordinate of the centroid of that district. The dataset containing GPS coordinates of the districts of India are collected from Kaggle (https:// www. kaggle. com/ sirpu nch/ indian-census-data-with-geosp atial-index ing).

Data preparation. From the raw datasets, a total of eleven potential demographic, socioeconomic, and environmental pollution related factors (see Table 1 ) are selected to explain the district-level geographical variation of COVID-19 mortality. The district-level demographic and socioeconomic factors that are selected in this study are: population; households with at least 9 persons; growth rate; sex ratio; persons with age 50 years or more; households having TV, computer (or laptop), mobile phones and car; number of persons having higher education; the percentage of the urban population. On the other hand, the environmental pollution related variables that are selected are as follows: PM 2.5 exposure; NO 2 exposure; SO 2 exposure.

The district-level long-term exposure to three air pollutants namely PM2.5, NO 2 , and SO 2 are calculated from the raw data of 130 pollution monitoring stations. The mean concentration of each of the above-mentioned air pollutants of all the 130 monitoring stations is computed for the period 2016-2020. For each pollutant, the computed values are spatially aggregated by averaging the values of all monitoring stations of a district. If a district doesn't contain any monitoring stations, then it's exposure to that pollutant is computed using Nearest Neighbour interpolation (NNI).

A multicollinearity verification is performed via the Variance Inflation Factor (VIF) to remove unnecessary redundancy among the explanatory variables. VIF can be expressed as follows [Eq. (1)]:

where, R 2 k denotes the coefficient of determination that is computed by regressing the k th variable on remaining explanatory variables. The mathematical expression for R 2 k is given in Eq. (2) . Here, SSE k and SST k denote the sum of squares of total variation and sum of squares of errors respectively. Firstly, regression analysis is conducted among all the 11 explanatory variables to compute the VIFs that are shown in Table 2 . It is observed that the variable HH_Abv_8_P has high Variance Inflation Factor (VIF = 12.4). Now, if VIFs are larger than 10, it indicates that there is multicollinearity 16 . Eventually, the variable HH_Abv_8_P is removed from the set of explanatory variables. After that, the regression is again performed on the remaining 10 variables, with the (1) www.nature.com/scientificreports/ VIFs given in Table 3 . Now, it is observed that no VIF exceeds 10 eventually this set of 10 variables can be used for model building.

Modeling spatial relationship. In this paper, the OLS (Ordinary Least Square) and GWR (Geographically Weighted Regression) models are utilized to determine the geographical relationship of COVID-19 mortality with potential risk factors. The OLS method generally attempts to understand the global relationships between the dependent and independent variables. In this case, the regression and its parameters are unchanged over the geographic space. Mathematically, Eq. (3) represents a global regression model as follows:

where, Y i denotes the dependent or response variable; X ik is the ith observation of kth independent variable; η k the global regression coefficient for kth independent variable; η 0 represent the intercept parameter; and δ 0 denotes the error term.

GWR technique extends the global regression [Eq. (3)] by enabling local parameter estimation 13 . It allows regression coefficients to be a function of geographical location. In other words, the regression coefficients are quantified independently in different geographical locations. A GWR model [Eq. (4)] can be represented as follows:

where, Y i , X ki , and δ i denote the dependent (or response) variable, kth independent (or predictor) variable, and error at location i respectively; (µ i , ν i ) denotes coordinates of location i; ξ k (µ i , ν i ) represent local coefficient for kth predictor at location i. Note that, GWR model allows regression parameters to vary continuously across the geographic space. For each location i, a set of regression parameters is estimated. The estimation of parameters can be performed as follows: www.nature.com/scientificreports/ where, X denotes a matrix containing the values of independent variables and a column of all 1s; Y represents a vector of values of the dependent variable; ξ (µ, ν) is a vector of local regression parameters; W(µ, ν) is a diagonal matrix whose diagonal elements represent the geographical weighting of the observations for regression location. The weights in W(µ, ν) assigns greater weights to the observations that are closer to the regression point than the observations that are farther away. In this work, the weights are computed using a Gaussian kernel function which is defined as follows:

where, B represents the bandwidth and D j i denotes the distance between the regression point i and the location of observation j. Note that, the bandwidth can be defined either by a fixed number of closest neighbors (known as adaptive bandwidth) or by a fixed distance (known as fixed bandwidth). Golden Section search 17 is utilized to find the optimum size of the bandwidth for GWR.

Performance metrics. The performance of the models are assessed by three metrics namely R 2 , adjusted R 2 , and AICc. Here, AICc is a corrected version of the Akaike Information Criterion (AIC). AICc can be defined as follows 13 : where, N denotes the sample size, S is the hat matrix, tr(S) denotes the trace of S, and σ represents the estimated standard deviation of the error term. AICc denotes model's accuracy and lower AICc indicates better model quality. It is usually used to find the best-fit model. The value of R 2 represents the ability of a model to explain the variance in the dependent variable and therefore a larger R 2 signifies the better performance of the model. It is computed from the estimated and the actual values of the dependent variable. Moreover, Moran's I index is computed to investigate the spatial autocorrelation of the model residuals. Mathematically, it is defined as follows:

where, N denotes total number of observations, y i and y j are variable values at location i and j respectively,ȳ represents the mean value, and w ij denotes a weight between location i and j. The value of Moran's I index can vary between − 1 (perfect dispersion) to + 1 (a perfect positive autocorrelation). Note that, a zero value indicates perfect spatial randomness.

Here, a step-wise GWR model selection using AICc is presented that can be utilized to investigate geographically varying relationships of COVID-19 mortality with different driving factors. The following are the steps to build an appropriate GWR model 18 .

• Step 1 Suppose there are n explanatory variables (in our case n = 10 ). For each of the explanatory variables, fit a separate GWR model by regressing that variable against the COVID19_Death variable. Compute AICc for each of the n = 10 models. Find the model that generates the lowest AICc and permanently include the corresponding explanatory variable in subsequent model building. • Step 2 Subsequently select a variable from the remaining ( n − 1 ) variables, build a model with the permanently included variables along with the newly selected variable. Find the explanatory variable that produces the lowest AICc and permanently include it in subsequent model building. Set n = n − 1.

Step 2 until it is observed that there is no reduction in AICc.

The above-mentioned steps are carried out using MGWR 2.2 software 19 . When calibrating the GWR, an adaptive bisquare spatial kernel is applied. Moreover, in order to select an optimal bandwidth, the Golden Section search 17 is employed. Figure 2 shows the changes in AICc during the step-by-step selection of explanatory variables for model building. It is observed that after the inclusion of a total of five variables, the AICc values start increasing when further new variables are included. Note that, both a global (OLS) model and a local (GWR) model are calibrated with these five explanatory variables.

In this section, firstly the performance of the global model (OLS) and local model (GWR) are discussed. Next, the geographically varying relationships of COVID-19 mortality with different factors are presented. Table 4 .

The variables Tot_population, HH_With_TCMC, Age_ABV_50, and PM 2.5 returns significant t values of 2.91, 12.114, − 1.225 and − 2.485 respectively. Moreover, the Moran's I of the residuals of the global OLS model are also analysed. It is found that there is significant spatial autocorrelation (Moran's I = 0.348 and p < 0.01). The Table 5 .

The performance of OLS and GWR model in terms of R 2 , Adj R 2 , and AICc are also provided in Table 6 . Moreover, Fig. 3a As shown in Fig. 5a , the GWR model produces local intercept that can vary within the range -0.78 to 1.13 with a mean of 0.013. In Fig. 5b , the regions with deep red color (mainly the state of West Bengal) denote those areas where the variable HH_With_TCMC has a strong positive relationship with COVID-19 death. The variable Higher_Edu is a strong predictor (See Fig. 5c ) for COVID-19 death in some parts of western India (mainly the state of Gujarat), southern India (mainly the states of Tamil Nadu and Kerala), and Eastern India (mainly the state of West Bengal). On the other hand, in the southern and the south-western part of India, a positive relationship between population and COVID-19 death is found (see Fig. 5d ). However, in some regions of central and western India (the states of Madhya Pradesh and Gujarat), a strong negative relationship between population and COVID-19 death is also observed. Fig. 5e shows that mainly in the western part of India there is a strong positive relationship between PM 2.5 and COVID-19 death, whereas in the other parts of India there is no such www.nature.com/scientificreports/ strong relationship. The explanatory variable Age_Abv_50 shows a positive relationship in central, eastern, and northern parts of India (see Fig. 5f ). Moreover, Table 7 represents the district-level results of the local model (GWR) for some of the districts that are severely affected by COVID-19 disease. The local R 2 values revealed district-level variability in GWR model performance. Specifically, the local R 2 values could be helpful here to see where geographically weighted regression predicts well and where it predicts poorly. It is observed that the GWR model yields high local R 2 value for most of the heavily affected districts. For instance, very high local R 2 values are found for the following districts: Pune ( R 2 = 0.982 ), Thane ( R 2 = 0.986 ), Lucknow ( R 2 = 0.984 ), Chittor ( R 2 = 0.973 ), Nasik ( R 2 = 0.986 ), Solapur ( R 2 = 0.974 ), Kolhapur ( R 2 = 0.972 ), Sangli ( R 2 = 0.973 ), Satara ( R 2 = 0.978 ), Latur ( R 2 = 0.967 ), Mumbai ( R 2 = 0.981 ), Kolkata ( R 2 = 0.987 ), Chennai ( R 2 = 0.978 ), Jalgaon ( R 2 = 0.966 ), Nanded ( R 2 = 0.960 ). On the other hand, moderate local R 2 values are found for Dharwad ( R 2 = 0.948 ), Nagpur ( R 2 = 0.955 ), Srikakulam ( R 2 = 0.920 ), Ludhiana ( R 2 = 0.921 ), Guntur ( R 2 = 0.95 ), Kurnool ( R 2 = 0.856 ), Coimbatore ( R 2 = 0.934 ), West Godavari ( R 2 = 0.95 ), Anantapur ( R 2 = 0.886 ), Bhopal ( R 2 = 0.916 ), and Krishna ( R 2 = 0.937 ). The lowest R 2 values are observed for the following districts: Hassan ( R 2 = 0.683 ), Indore ( R 2 = 0.798 ), and Jaipur ( R 2 = 0.797 ). Note that, for most of the highly COVID-19-affected districts, the variables PM 2.5 and HH_With_TCMC are usually exhibited positive relationships in regression modeling. On the other hand, the variable Higher_Edu usually exhibits negative relationships for most of the highly affected districts. 

In order to better understand how different driving factors influence the overall fatalities caused by COVID-19, the geographical distribution of COVID-19-related deaths are investigated. The highest number of COVID-19-related deaths are found primarily in the western part of India (Pune, Thane, Mumbai, Nagpur, Nashik, Raigad, Jalgaon, Kolhapur, Sangli, Satara, Solapur, Ahmedabad, Surat). On the other hand, the number of COVID-19-related deaths is relatively low in the northern and eastern parts of India. This study identified considerable geographical variability of COVID-19 deaths and their heterogeneous relationship at the local level with the driving factors in India. More specifically, the utilization of the GWR method successfully found the geographically varying relationship of COVID-19 mortality with various potential socio-economic, demographic, and environmental pollution related factors. This study reveals five important local factors are significantly related with district-level COVID-19 deaths as follows: (i) population (ii) PM 2.5 level (iii) households having TV, computer (or laptop), mobile phones and car (iv) persons with age 50 years or more (v) number of persons having higher education. Furthermore, this study also validates the effectiveness of local parameter estimation by comparing the global OLS method with the local GDR method. To the best of our knowledge, this is the first study that explores geographically varying relationships of COVID-19 deaths with various potential driving factors in India. Rigorous analyses are performed to demonstrate the shortcomings of global technique (OLS) as compared to the local technique (GWR) in terms of several performance metrics. The OLS model only explains 71.9% of the variance of district-level COVID-19 deaths. It is found that the predictive efficiency and model accuracy are further enhanced by implementing the GWR method. The GWR model explains 97% of the variance of district-level COVID-19 deaths. Moreover, Moran's I index verifies that no significant spatial autocorrelation is present in the residuals of the GWR model. Note that, a key advantage of such a local method is its capability to visualize the geographically varying heterogeneous relationships between the dependent and the independent variables. In other words, it enables us for a better understanding of relationships based on geographical contexts and study area's known features.

The findings of this study reveal that there are strong positive relationships of COVID-19 deaths with the explanatory variables PM2.5 and Tot_population across the regions of the COVID-19 death hotspots in the western part of India. The positive association of COVID-19 deaths with long term exposure of PM 2.5 is consistent with the previous works 20, 21 . Note that, long-term PM 2.5 exposure is substantially associated with some of the comorbidities (e.g. chronic lung disease, cardiovascular disease, etc.) that may lead to COVID-19 deaths 22, 23 . Similarly, a positive association between COVID-19-related deaths and Tot_population is also observed in other studies 6, 24 . However, the reverse association is found for these two variables ( PM 2.5 and Tot_population ) in the other parts of India. The explanatory variable HH_With_TCMC is found to be an important factor that may be a measure of the number of households with the upper class and rich people. A strong positive relationship is observed between HH_With_TCMC and COVID-19 death in the hotspots of eastern and western parts of India (Kolkata, North 24 Parganas, Pune, Thane, Surat, Nagpur, etc.). Note that, in those hotspots, the value of HH_With_TCMC is substantially high. An interesting observation reveals that a strong negative relationship exists between COVID-19 death and Higher_Edu in the eastern, central, and southern parts of India. It is expected that the higher educated people are well aware of the symptoms and the complications of COVID-19 that may lead to the fewer number of fatalities in those regions. Now, in some regions of the south-eastern part of India, the number of COVID-19 deaths is also seen to be high.

In those regions, significant positive relationships are found between COVID − 19 deaths and Tot_population , whereas significant negative relationships are observed for the variable Higher_Edu.

This research work inherits certain shortcomings that need to be resolved in future research. For instance, there may have high possibilities of under-reporting in COVID-19 death counts that may introduce bias in the study 25 . Moreover, due to data unavailability, we were not able to include some significant district-level driving factors in our study, such as health care system quality, number of hospital beds, household income, and poverty data. Despite the above-mentioned shortcomings, this is the first study that explores geographically varying relationships of COVID-19 mortality with different socioeconomic, demographic, and environmental pollution related factors in India. This research work also highlights the significance of the geographically weighted regression in the geographical analysis of the health outcome of COVID-19 disease.

COVID-19 pandemic is one of the most serious global public health catastrophe of the century. In this work, the geographically varying relationships between COVID-19 deaths and different potential driving factors are assessed across India. The geographical distribution of reported COVID-19 death cases is found to be heterogeneous over India. This heterogeneity in distribution is related to many underlying factors, including demographic, socioeconomic, and environmental pollution related variations between different parts of India. The GWR model makes it possible for the regression coefficients to differ across the geospace, creating geographical patterns about the strength of the relationship. The geographical heterogeneity and non-stationary of the relationships between COVID-19 deaths and the driving factors are demonstrated by mapping the local parameter estimates. The local parameter estimates reflect the quality of local model fitting and the nature of the association. The local method (GWR) yields better performance with smaller AICc as compared to the global method (OLS).

It should be noted that the impacts of other influencing factors (e.g. Meteorological factors) are not included in this work. This might be the direction for future studies. Moreover, in this study, currently we do not consider time evolution of variables, it is because for the dependent and the independent variables we may require more time series data for the effective temporal modelling. However, we plan to consider the time evolution of the variables for the future studies when more time series data will be available. 

WHO Coronavirus Disease (COVID-19) Dashboard

Clinical course and risk factors for mortality of adult inpatients with covid-19 in Wuhan,China: A retrospective cohort study

Hospitalization rates and characteristics of patients hospitalized with laboratory-confirmed coronavirus disease 2019-covid-net, 14 states

The clinical and chest ct features associated with severe and critical covid-19 pneumonia

The socioeconomic determinants of covid-19: A spatial analysis of german county level data. medRxiv

Examining the association between socio-demographic composition and covid-19 fatalities in the european region using spatial regression approach

Significance of geographical factors to the covid-19 outbreak in india

A spatial analysis of the covid-19 period prevalence in us counties through june 28, 2020: Where geography matters?

Ordinary Least-Squares Regression 224-228 (L. Moutinho and GD Hutcheson, The SAGE dictionary of quantitative management research

Geographically weighted regression: A method for exploring spatial nonstationarity

Geographically weighted regression

Geographically Weighted Regression: The Analysis of Spatially Varying Relationships

Local indicators of spatial association-lisa

Central Pollution Control Board, Ministry of Environment, Forest and Climate Change

Applied logistic regression analysis

An Extension of Geographically Weighted Regression with Flexible Bandwidths

A python implementation of multiscale geographically weighted regression for investigating process spatial heterogeneity and scale

Exposure to air pollution and covid-19 mortality in the United States

The relationship between air pollution and covid-19-related deaths: An application to three french cities

Particulate matter air pollution and cardiovascular disease: An update to the scientific statement from the american heart association

Comorbidity and its impact on patients with covid-19. SN Compreh

Influence of socio-ecological factors on covid-19 risk: a cross-sectional study based on 178 countries/regions worldwide

Is India missing covid-19 deaths?

S.R. proposed the research topic, provided conceptual and technical guidance. A.I.M. designed the research plan, wrote the manuscript, collected the data, performed the statistical analysis. A.I.M. and S.R. both involved in the revision of the manuscript and interpretation of the results. Both the authors read and approved the final manuscript.

The authors declare no competing interests.

Correspondence and requests for materials should be addressed to S.R.Reprints and permissions information is available at www.nature.com/reprints.Publisher's note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.