Exploring subnational human development in India

Setup

# Load libraries import numpy as np import pandas as pd pd.set_option('display.max_columns', None) import matplotlib.pyplot as plt plt.style.use('ggplot') from matplotlib.ticker import FormatStrFormatter from matplotlib_scalebar.scalebar import ScaleBar from pylab import rcParams import matplotlib as mpl mpl.rcParams['figure.dpi'] = 72 import seaborn as sns #sns.set_style("darkgrid") #sns.set_context(context="paper", font_scale=1.5, rc=None) #sns.set(font="serif") import plotly.express as px import plotly.graph_objects as go import geopandas as gpd import libpysal from libpysal import weights from libpysal.weights import Queen import esda from esda.moran import Moran, Moran_Local import splot from splot.esda import moran_scatterplot, plot_moran, lisa_cluster, plot_local_autocorrelation from splot.libpysal import plot_spatial_weights from giddy.directional import Rose import statsmodels.api as sm import statsmodels.formula.api as smf from spreg import OLS from spreg import MoranRes from spreg import ML_Lag from spreg import ML_Error from mgwr.gwr import GWR, MGWR from mgwr.sel_bw import Sel_BW from mgwr.utils import shift_colormap, truncate_colormap import warnings warnings.filterwarnings('ignore')

/opt/conda/lib/python3.9/site-packages/esda/getisord.py:636: SyntaxWarning: "is" with a literal. Did you mean "=="?
  if __name__ is "__main__":
/opt/conda/lib/python3.9/site-packages/spglm/utils.py:367: SyntaxWarning: "is not" with a literal. Did you mean "!="?
  if resetlist is not ():

Import data

df = pd.read_csv("https://github.com/quarcs-lab/data-open/raw/master/shdi/SHDIold.csv").sort_values(by=['year', 'SubContinent'])

df.info()

<class 'pandas.core.frame.DataFrame'>
Int64Index: 54668 entries, 2087 to 45190
Data columns (total 27 columns):
 #   Column           Non-Null Count  Dtype  
---  ------           --------------  -----  
 0   iso_code         54668 non-null  object 
 1   country          54668 non-null  object 
 2   year             54668 non-null  int64  
 3   GDLcode          54668 non-null  object 
 4   continent2       54668 non-null  object 
 5   level            54668 non-null  object 
 6   region           54668 non-null  object 
 7   continent        54668 non-null  object 
 8   CountryName_std  54668 non-null  object 
 9   shdi             54668 non-null  float64
 10  healthindex      54668 non-null  float64
 11  incindex         54668 non-null  float64
 12  edindex          54668 non-null  float64
 13  esch             54668 non-null  float64
 14  msch             54668 non-null  float64
 15  lifexp           54668 non-null  float64
 16  gnic             54668 non-null  float64
 17  lgnic            54668 non-null  float64
 18  pop              54164 non-null  float64
 19  POLY_IDcountry   54588 non-null  float64
 20  iso3             54668 non-null  object 
 21  status           54668 non-null  object 
 22  CountryName      54668 non-null  object 
 23  SubContinent     54668 non-null  object 
 24  Continent        54668 non-null  object 
 25  Continent_std    54668 non-null  object 
 26  _merge           54668 non-null  object 
dtypes: float64(11), int64(1), object(15)
memory usage: 11.7+ MB

df

iso_codeobject

USA2.9%

IND2%

184 others95.1%

countryobject

2087

AUS

Australia

2106

AUS

Australia

2165

AUS

Australia

2171

AUS

Australia

2198

AUS

Australia

2205

AUS

Australia

2260

AUS

Australia

2291

AUS

Australia

2320

AUS

Australia

35159

NZL

New Zealand

map = gpd.read_file("https://gist.github.com/cmg777/81b6e2b1ea5865c09daf3b2bbf35bbf8/raw/5e0a0370f453421bc503594cfae1e466c89e449f/GDL%2520Shapefiles%2520V4.geojson")

map.plot('shdi');

map.info()

<class 'geopandas.geodataframe.GeoDataFrame'>
RangeIndex: 1745 entries, 0 to 1744
Data columns (total 7 columns):
 #   Column    Non-Null Count  Dtype   
---  ------    --------------  -----   
 0   GDLcode   1745 non-null   object  
 1   constant  1745 non-null   object  
 2   iso_code  1742 non-null   object  
 3   country   1742 non-null   object  
 4   region    1742 non-null   object  
 5   shdi      1742 non-null   float64 
 6   geometry  1745 non-null   geometry
dtypes: float64(1), geometry(1), object(5)
memory usage: 95.6+ KB

map.head(9)

GDLcodeobject

AFGr10111.1%

AFGr10211.1%

7 others77.8%

constantobject

World100%

0

AFGr101

World

1

AFGr102

World

2

AFGr103

World

3

AFGr104

World

4

AFGr105

World

5

AFGr106

World

6

AFGr107

World

7

AFGr108

World

8

AGOr201

World

fig, ax = plt.subplots(figsize=(12,8)) map.plot(column="shdi", scheme='BoxPlot', cmap='coolwarm', legend=True, ax=ax, legend_kwds={'bbox_to_anchor':(1.00, 0.70)}) plt.title('Spatial distribution of the subnational HDI: \nWhich regions are below/above the median?') plt.tight_layout() ax.axis("off") #plt.savefig('myMap.png',dpi=150, bbox_inches='tight') plt.show()

map['shdi'].describe().round(2)

Prepare data

# Remove national observations and focus on India df1 = df[(df.level != 'National') & (df.country == 'India')] df1

iso_codeobject

IND100%

countryobject

India100%

21064

IND

India

21119

IND

India

21167

IND

India

21170

IND

India

21216

IND

India

21300

IND

India

21314

IND

India

21331

IND

India

21346

IND

India

21355

IND

India

df1.columns

df1["year"]= df["year"].astype(str)

df1990 = df1[df1['year'] == '1990']

df2019 = df1[df1['year'] == '2019']

df2019['Continent'].unique()

df2019['SubContinent'].unique()

df2019['country'].unique()

Descriptive statistics

df1.loc[:, df1.columns != 'year'].describe().round(2)

shdifloat64

healthindexfloat64

count

1080

mean

0.6

0.74

std

0.09

0.07

min

0.38

0.52

25%

0.53

0.69

50%

0.6

0.75

75%

0.67

0.8

max

0.78

0.88

df1990.describe().round(2)

shdifloat64

healthindexfloat64

count

36

mean

0.51

0.67

std

0.09

0.08

min

0.38

0.52

25%

0.44

0.61

50%

0.49

0.66

75%

0.56

0.74

max

0.71

0.82

df2019.describe().round(2)

shdifloat64

healthindexfloat64

count

36

mean

0.68

0.8

std

0.05

0.04

min

0.57

0.7

25%

0.65

0.77

50%

0.68

0.8

75%

0.72

0.82

max

0.78

0.88

Regional differences

Strip plot

px.strip(df1, x = ['shdi'], y = 'country', range_x= [0.35, 0.80], color = 'continent2', hover_name= 'region', hover_data = ['country'], animation_frame= 'year', labels=dict(continent2 = "Continent", shdi ="Subnational human development index") )

Box plot

px.box(df1, x = 'shdi', y = 'country', range_x= [0.35, 0.80], color = 'continent2', hover_name= 'region', hover_data = ['country'], animation_frame= 'year', labels=dict(continent2 = "Continent", shdi ="Subnational human development index") )

Histogram

px.histogram(df1, x = 'shdi', range_x= [0.35, 0.80], nbins= 15, color = 'country', hover_name= 'region', hover_data= ['country'], marginal='rug', # box is also a good option animation_frame = 'year', labels=dict(continent2 = "Continent", shdi ="Subnational human development index") )

ECDF

px.ecdf(df1, x="shdi", range_x= [0.30, 0.80], color="country", hover_name= 'region', hover_data= ['country'], #markers=True, #lines=False, marginal="rug", # histogram animation_frame = 'year', labels=dict(continent2 = "Continent", shdi ="Subnational human development index") )

Time series

Continents

fig = px.line( df1, x="year", y="shdi", log_y= True, color="region", hover_name="country", hover_data= ['country'], labels=dict(shdi="HDI", continent2 ="Continent", year = ""), facet_col="country", facet_col_wrap = 2, facet_row_spacing = 0.01, height= 500 ) fig.update_layout(showlegend=False)

Correlations

px.scatter(df1, y = "shdi", x = "incindex", range_y= [0.2, 0.9], range_x= [0.2, 0.9], hover_name = "region", hover_data= ['country'], color = "country", #size = "pop", size_max = 60, trendline= 'ols', #trendline_scope= "overall", #marginal_x= 'box', #marginal_y= 'box', animation_frame= 'year', labels=dict(country = "Country", continent2 = "Continent", lgnic ="Log of gross national income per capita", lifexp="Life Expectancy (years)") )

Maps

gdf1990 = map.merge(df1990, on = 'GDLcode', how = 'left') gdf2019 = map.merge(df2019, on = 'GDLcode', how = 'left')

gdf2019.columns

fig, ax = plt.subplots(figsize=(12,8)) gdf1990.plot(column="shdi_y", scheme='BoxPlot', cmap='coolwarm', legend=True, ax=ax, edgecolor='black', legend_kwds={'bbox_to_anchor':(0.80, 0.95)}) plt.title('Spatial distribution of human development in 1990: \nWhich regions are below/above the median?') plt.tight_layout() ax.axis("off") #plt.savefig('myMap.png',dpi=300, bbox_inches='tight') plt.show()

fig, ax = plt.subplots(figsize=(12,8)) gdf2019.plot(column="shdi_y", scheme='BoxPlot', cmap='coolwarm', legend=True, ax=ax, edgecolor='black', legend_kwds={'bbox_to_anchor':(0.80, 0.95)}) plt.title('Spatial distribution of human development in 1990: \nWhich regions are below/above the median?') plt.tight_layout() ax.axis("off") #plt.savefig('myMap.png',dpi=300, bbox_inches='tight') plt.show()

Spatial dependence

gdf2019India = gdf2019.query("country_x == 'India'").reset_index(drop=True)

#W = weights.KNN.from_dataframe(gdf2019India, k=4) W = weights.KNN.from_dataframe(gdf2019India, k=6) #W = weights.KNN.from_dataframe(gdf2019India, k=8) #W = weights.Queen.from_dataframe(gdf2019India) #W = weights.Rook.from_dataframe(gdf2019India) #>>>> Compute inverse distance squared based on maximum nearest neighbor distance between the n observations #maxMin = weights.min_threshold_dist_from_shapefile("nepal.shp") #W = weights.DistanceBand.from_dataframe(nepal, threshold = maxMin, binary=False, alpha=-2) # >>> Row-standardize W W.transform = 'r'

plot_spatial_weights(W, gdf2019India);

gdf2019India['Wshdi_y'] = weights.lag_spatial(W, gdf2019India['shdi_y'])

globalMoran = Moran(gdf2019India['shdi_y'], W) MoranI = globalMoran.I MoranI = "{:.2f}".format(MoranI) localMoran = Moran_Local(gdf2019India['shdi_y'], W, permutations = 999, seed=12345)

fig, axes = plt.subplots(nrows=1, ncols=2, figsize=(12,5)) moran_scatterplot(localMoran, p=0.10, aspect_equal=False, zstandard = False, ax=axes[0]) lisa_cluster(localMoran, gdf2019India, p=0.10, legend_kwds={'bbox_to_anchor':(0.02, 0.90)}, ax=axes[1]) #cx.add_basemap(axes[1], crs=gdf.crs.to_string(), source=cx.providers.Stamen.TonerLite, attribution=False) #cx.add_basemap(axes[1], crs=gdf.crs.to_string(), source=cx.providers.Stamen.TonerLabels, attribution=False) axes[0].set_xlabel('Indicator') axes[0].set_ylabel('Spatial lag of indicator') axes[0].set_title(f"(a) Moran scatterplot (Moran's I = {MoranI})") axes[1].set_title("(b) Spatial clusters and outliers (p<0.10)") plt.show()

gdf2019India.explore( column= 'shdi_y', tooltip=['shdi_y','country_y', 'region_y'], k = 3, scheme='FisherJenks', cmap='coolwarm', legend=True, tiles='CartoDB dark_matter', style_kwds =dict(color="gray", weight=0.4), legend_kwds=dict(colorbar=False) )

px.scatter( gdf2019India, x='shdi_y', y='Wshdi_y', hover_name="region_y", hover_data = ['shdi_y', 'Wshdi_y', 'country_y'], trendline="ols", marginal_x="box", marginal_y="box")

Exploring subnational human development in India

Setup

Import data

Prepare data

Descriptive statistics

Regional differences

Strip plot

Box plot

Histogram

ECDF

Time series

Continents

Correlations

Maps

Spatial dependence

Spatial inequality