xjtulyc Global Health Data Access and Analysis
This skill provides tools for accessing and analyzing global health data from major open sources: the WHO Global Health Observatory API, IHME Global Burden of Disease study data, and Our World in Data health datasets. It covers computing burden-of-disease metrics (DALYs, YLLs, YLDs), age-standardized rates, health inequality measures, and publication-quality choropleth maps.
Prerequisites
pip install requests pandas numpy scipy matplotlib geopandas shapely seaborn statsmodels
For GBD data downloads you may need an IHME account. Set API tokens via environment variables (no hardcoded credentials):
export WHO_API_KEY="<paste-your-who-api-key>" # currently optional for public endpoints
export IHME_API_TOKEN="<paste-your-ihme-token>" # for authenticated GBD API endpoints
Core Functions
1. WHO Global Health Observatory (GHO) API
import os
import time
import requests
import pandas as pd
import numpy as np
from typing import Optional
# WHO GHO OData API v3 base URL (public, no API key required for most indicators)
WHO_GHO_BASE = "https://ghoapi.azureedge.net/api"
def get_who_indicator(
indicator_code: str,
countries: Optional[list[str]] = None,
years: Optional[list[int]] = None,
dimension_filter: Optional[dict] = None,
page_size: int = 1000,
verbose: bool = True,
) -> pd.DataFrame:
"""
Retrieve data for a WHO GHO indicator via the OData REST API.
Parameters
----------
indicator_code : str
WHO indicator code, e.g. 'MDG_0000000007' (under-5 mortality rate),
'NCD_BMI_30A' (obesity prevalence), 'WHOSIS_000001' (life expectancy).
Browse codes at: https://ghoapi.azureedge.net/api/Indicator
countries : list[str] | None
List of ISO 3166-1 alpha-3 country codes, e.g. ['NGA', 'KEN', 'ETH'].
If None, retrieves all countries.
years : list[int] | None
List of years to filter, e.g. [2000, 2005, 2010, 2015, 2019].
If None, retrieves all available years.
dimension_filter : dict | None
Extra OData filter key-value pairs, e.g. {'Dim1': 'MLE'} for males.
page_size : int
Number of records per API page request.
verbose : bool
Whether to print progress messages.
Returns
-------
pd.DataFrame
Columns: indicator_code, country_code, country_name, year, value,
low, high, sex, dim1, dim2, comments.
"""
filters = []
if countries:
country_filter = " or ".join([f"SpatialDim eq '{c}'" for c in countries])
filters.append(f"({country_filter})")
if years:
year_filter = " or ".join([f"TimeDim eq {y}" for y in years])
filters.append(f"({year_filter})")
if dimension_filter:
for key, val in dimension_filter.items():
filters.append(f"{key} eq '{val}'")
odata_filter = " and ".join(filters) if filters else None
records = []
skip = 0
session = requests.Session()
headers = {}
# Use API key if provided (currently optional for most WHO endpoints)
api_key = os.environ.get("WHO_API_KEY")
if api_key:
headers["Ocp-Apim-Subscription-Key"] = api_key
while True:
params = {"$top": page_size, "$skip": skip}
if odata_filter:
params["$filter"] = odata_filter
url = f"{WHO_GHO_BASE}/{indicator_code}"
try:
resp = session.get(url, params=params, headers=headers, timeout=30)
resp.raise_for_status()
except requests.HTTPError as e:
print(f"HTTP error for indicator {indicator_code}: {e}")
break
except requests.RequestException as e:
print(f"Request error: {e}")
time.sleep(5)
continue
data = resp.json()
batch = data.get("value", [])
if not batch:
break
for item in batch:
records.append({
"indicator_code": indicator_code,
"country_code": item.get("SpatialDim", ""),
"country_name": item.get("SpatialDimType", ""),
"year": item.get("TimeDim"),
"value": item.get("NumericValue"),
"low": item.get("Low"),
"high": item.get("High"),
"sex": item.get("Dim1", ""),
"dim2": item.get("Dim2", ""),
"comments": item.get("Comments", ""),
})
skip += page_size
if verbose:
print(f" Fetched {len(records)} records so far...")
if len(batch) < page_size:
break
time.sleep(0.3)
df = pd.DataFrame(records)
df["year"] = pd.to_numeric(df["year"], errors="coerce")
df["value"] = pd.to_numeric(df["value"], errors="coerce")
df["low"] = pd.to_numeric(df["low"], errors="coerce")
df["high"] = pd.to_numeric(df["high"], errors="coerce")
if verbose:
print(f"Total records: {len(df)} for indicator '{indicator_code}'")
return df
def list_who_indicators(search_term: str = "") -> pd.DataFrame:
"""
List available WHO GHO indicators, optionally filtered by search term.
Returns DataFrame with columns: indicator_code, indicator_name.
"""
url = f"{WHO_GHO_BASE}/Indicator"
params = {}
if search_term:
params["$filter"] = f"contains(IndicatorName, '{search_term}')"
resp = requests.get(url, params=params, timeout=30)
resp.raise_for_status()
items = resp.json().get("value", [])
return pd.DataFrame([{
"indicator_code": it["IndicatorCode"],
"indicator_name": it["IndicatorName"],
} for it in items])
2. IHME Global Burden of Disease Data
def download_gbd_data(
cause: str,
location: str | list[str],
metric: str = "Rate",
year: int | list[int] = 2019,
measure: str = "DALYs",
sex: str = "Both",
age_group: str = "All Ages",
fallback_csv_path: Optional[str] = None,
) -> pd.DataFrame:
"""
Download IHME GBD data via API or fall back to a pre-downloaded CSV.
The IHME GBD API requires registration at healthdata.org.
For unauthenticated use, download CSVs from https://ghdx.healthdata.org/gbd-results
and pass the path via `fallback_csv_path`.
Parameters
----------
cause : str
GBD cause of death/disability name (e.g., 'Diabetes mellitus',
'Lower respiratory infections', 'Ischemic heart disease').
location : str | list[str]
Country name(s) or GBD super-region name.
metric : str
'Rate' (per 100k), 'Number', or 'Percent'.
year : int | list[int]
Year(s) to retrieve.
measure : str
'DALYs', 'Deaths', 'YLLs', 'YLDs', 'Prevalence', 'Incidence'.
sex : str
'Both', 'Male', 'Female'.
age_group : str
GBD age group (e.g., 'All Ages', '<5 years', '70+ years').
fallback_csv_path : str | None
Path to a locally downloaded GBD CSV as fallback.
Returns
-------
pd.DataFrame
Standardized DataFrame with columns: location, year, cause, measure,
metric, sex, age_group, value, lower_ci, upper_ci.
"""
ihme_token = os.environ.get("IHME_API_TOKEN")
if ihme_token:
headers = {"Authorization": f"Bearer {ihme_token}"}
base_url = "https://api.healthdata.org/gbd/v1/results"
locations = [location] if isinstance(location, str) else location
years = [year] if isinstance(year, int) else year
params = {
"cause_name": cause,
"location_name": ",".join(locations),
"metric_name": metric,
"year_id": ",".join(str(y) for y in years),
"measure_name": measure,
"sex_name": sex,
"age_name": age_group,
"format": "json",
}
try:
resp = requests.get(base_url, params=params,
headers=headers, timeout=30)
resp.raise_for_status()
data = resp.json()
records = data.get("results", [])
if records:
df = pd.DataFrame(records)
df = df.rename(columns={
"location_name": "location",
"year_id": "year",
"cause_name": "cause",
"measure_name": "measure",
"metric_name": "metric",
"sex_name": "sex",
"age_name": "age_group",
"val": "value",
"lower": "lower_ci",
"upper": "upper_ci",
})
return df[["location", "year", "cause", "measure",
"metric", "sex", "age_group",
"value", "lower_ci", "upper_ci"]]
except requests.RequestException as e:
print(f"IHME API request failed: {e}. Falling back to CSV.")
if fallback_csv_path:
df = pd.read_csv(fallback_csv_path)
# Standardize column names from IHME CSV export format
col_map = {
"location_name": "location",
"year_id": "year",
"cause_name": "cause",
"measure_name": "measure",
"metric_name": "metric",
"sex_name": "sex",
"age_name": "age_group",
"val": "value",
"lower": "lower_ci",
"upper": "upper_ci",
}
df = df.rename(columns={k: v for k, v in col_map.items() if k in df.columns})
# Apply filters
locations = [location] if isinstance(location, str) else location
years = [year] if isinstance(year, int) else year
if "location" in df.columns:
df = df[df["location"].isin(locations)]
if "year" in df.columns:
df = df[df["year"].isin(years)]
if "cause" in df.columns and cause:
df = df[df["cause"].str.contains(cause, case=False, na=False)]
if "measure" in df.columns:
df = df[df["measure"] == measure]
return df
raise RuntimeError(
"No IHME_API_TOKEN set and no fallback_csv_path provided. "
"Download GBD data from https://ghdx.healthdata.org/gbd-results "
"and pass the CSV path as fallback_csv_path."
)
3. DALYs, YLLs, and YLDs
def compute_dalys(
yll_series: pd.Series,
yld_series: pd.Series,
index_cols: Optional[pd.DataFrame] = None,
) -> pd.DataFrame:
"""
Compute DALYs as the sum of YLLs and YLDs.
DALYs (Disability-Adjusted Life Years) = YLLs + YLDs.
YLLs (Years of Life Lost) = premature mortality burden.
YLDs (Years Lived with Disability) = morbidity burden.
Parameters
----------
yll_series : pd.Series
YLL values (same index as yld_series).
yld_series : pd.Series
YLD values.
index_cols : pd.DataFrame | None
Optional metadata DataFrame to attach (same index).
Returns
-------
pd.DataFrame
DataFrame with yll, yld, daly columns plus metadata if provided.
"""
result = pd.DataFrame({
"yll": yll_series.values,
"yld": yld_series.values,
"daly": yll_series.values + yld_series.values,
})
if index_cols is not None:
result = pd.concat([index_cols.reset_index(drop=True), result], axis=1)
return result
4. Age-Standardized Rates
# WHO World Standard Population 2000-2025 age weights (18 age groups)
WHO_WORLD_STANDARD_POPULATION = pd.DataFrame({
"age_group": [
"0-4", "5-9", "10-14", "15-19", "20-24", "25-29",
"30-34", "35-39", "40-44", "45-49", "50-54", "55-59",
"60-64", "65-69", "70-74", "75-79", "80-84", "85+",
],
"weight": [
0.0886, 0.0869, 0.0860, 0.0847, 0.0822, 0.0793,
0.0761, 0.0715, 0.0659, 0.0604, 0.0537, 0.0455,
0.0372, 0.0296, 0.0221, 0.0152, 0.0091, 0.0063,
],
})
def age_standardize(
crude_rates: pd.Series,
age_weights: pd.Series,
reference_pop: Optional[pd.Series] = None,
per_unit: float = 100_000,
) -> dict:
"""
Compute directly age-standardized rate using the WHO world standard population.
Parameters
----------
crude_rates : pd.Series
Age-specific crude rates (cases per person, NOT per 100k).
Index should match age_weights index.
age_weights : pd.Series
Proportional weight for each age group (must sum to 1.0).
Use WHO_WORLD_STANDARD_POPULATION['weight'] as default.
reference_pop : pd.Series | None
If provided, use direct standardization with this reference population
(absolute counts); age_weights is ignored.
per_unit : float
Multiplier for output (default 100,000 for rates per 100k).
Returns
-------
dict
Keys: age_standardized_rate, variance (approximate), ci_95_lower, ci_95_upper.
"""
rates = np.asarray(crude_rates, dtype=float)
weights = np.asarray(age_weights, dtype=float)
if reference_pop is not None:
ref = np.asarray(reference_pop, dtype=float)
weights = ref / ref.sum()
# Normalize weights just in case
weights = weights / weights.sum()
asr = np.sum(weights * rates) * per_unit
# Approximate variance (Chiang 1961 method)
variance = np.sum((weights**2) * rates * (1 - rates + 1e-12)) * per_unit**2
se = np.sqrt(variance)
ci_lower = max(0.0, asr - 1.96 * se)
ci_upper = asr + 1.96 * se
return {
"age_standardized_rate": asr,
"standard_error": se,
"ci_95_lower": ci_lower,
"ci_95_upper": ci_upper,
"variance": variance,
}
5. Health Inequality Measures
def compute_concentration_index(
health_outcome: pd.Series,
wealth_rank: pd.Series,
weights: Optional[pd.Series] = None,
) -> dict:
"""
Compute the concentration index (CI) for measuring socioeconomic health inequality.
The concentration index ranges from -1 (all disease burden in the poorest) to
+1 (all disease burden in the richest). CI = 0 means perfect equality.
Parameters
----------
health_outcome : pd.Series
Health variable (e.g., prevalence, mortality rate per individual or group).
wealth_rank : pd.Series
Socioeconomic rank variable (higher = wealthier). Same index.
weights : pd.Series | None
Population weights for grouped data.
Returns
-------
dict
Keys: concentration_index, standard_error, t_statistic, p_value,
erreygers_ci (normalized version for bounded outcomes).
"""
df = pd.DataFrame({
"y": health_outcome.values,
"rank": wealth_rank.values,
}).dropna()
if weights is not None:
df["w"] = weights.reindex(df.index).fillna(1.0).values
else:
df["w"] = 1.0
total_w = df["w"].sum()
df["w_norm"] = df["w"] / total_w
# Fractional rank (Erreygers & van Ourti method)
df = df.sort_values("rank")
df["cum_w"] = df["w_norm"].cumsum()
df["frac_rank"] = df["cum_w"].shift(1).fillna(0) + df["w_norm"] / 2
y_mean = (df["y"] * df["w_norm"]).sum()
# Covariance form of CI
cov_xy = (df["w_norm"] * (df["y"] - y_mean) * (df["frac_rank"] - 0.5)).sum()
ci = 2.0 * cov_xy / (y_mean + 1e-12)
# Standard error via OLS regression (Kakwani et al.)
from statsmodels.api import OLS, add_constant
X_reg = add_constant(df["frac_rank"])
y_reg = 2.0 * df["y"] / y_mean
try:
model = OLS(y_reg, X_reg, weights=df["w_norm"]).fit(
cov_type="HC3"
)
se = model.bse["frac_rank"] / 2.0
t_stat = ci / (se + 1e-12)
p_val = float(2 * (1 - stats_module.t.cdf(abs(t_stat), df=len(df) - 2)))
except Exception:
se = np.nan
t_stat = np.nan
p_val = np.nan
# Erreygers normalized CI (for bounded variables like prevalence 0-1)
y_max = df["y"].max()
y_min = df["y"].min()
erreygers_ci = 4 * y_mean / (y_max - y_min + 1e-12) * ci if y_max > y_min else np.nan
return {
"concentration_index": ci,
"standard_error": se,
"t_statistic": t_stat,
"p_value": p_val,
"erreygers_ci": erreygers_ci,
"mean_health": y_mean,
"n": len(df),
}
import scipy.stats as stats_module
def compute_inequality_gaps(
df: pd.DataFrame,
value_col: str,
group_col: str,
reference_group: str,
) -> pd.DataFrame:
"""
Compute absolute and relative health gaps relative to a reference group.
Parameters
----------
df : pd.DataFrame
Must contain `value_col` (health metric) and `group_col` (socioeconomic group).
value_col : str
Column with health values (e.g., mortality rate).
group_col : str
Column with group labels (e.g., income quintile).
reference_group : str
Group to use as reference (usually the most advantaged).
Returns
-------
pd.DataFrame
Group-level table with absolute_gap and relative_gap columns.
"""
group_means = df.groupby(group_col)[value_col].mean().reset_index()
ref_value = group_means.loc[
group_means[group_col] == reference_group, value_col
].values[0]
group_means["absolute_gap"] = group_means[value_col] - ref_value
group_means["relative_gap"] = group_means[value_col] / (ref_value + 1e-12)
return group_means
6. Choropleth Map Visualization
def plot_global_health_map(
df: pd.DataFrame,
value_col: str,
year: int,
title: str,
country_col: str = "country_code",
year_col: str = "year",
cmap: str = "YlOrRd",
legend_label: str = "Value",
vmin: Optional[float] = None,
vmax: Optional[float] = None,
output_file: Optional[str] = None,
missing_color: str = "#cccccc",
) -> None:
"""
Plot a global choropleth map of a health indicator using geopandas.
Parameters
----------
df : pd.DataFrame
Health data with country ISO-3 codes and values.
value_col : str
Column to map.
year : int
Year to filter.
title : str
Map title.
country_col : str
Column with ISO 3166-1 alpha-3 country codes.
year_col : str
Column with year values.
cmap : str
Matplotlib colormap.
legend_label : str
Label for the colorbar.
vmin, vmax : float | None
Color scale min/max. Auto-scaled if None.
output_file : str | None
Save figure to this path if provided.
missing_color : str
Fill color for countries with no data.
"""
import geopandas as gpd
import matplotlib.pyplot as plt
import matplotlib.colors as mcolors
# Load Natural Earth world shapefile (bundled with geopandas)
world = gpd.read_file(
gpd.datasets.get_path("naturalearth_lowres")
)
world = world.rename(columns={"iso_a3": "iso3"})
# Filter data to the requested year
df_year = df[df[year_col] == year][[country_col, value_col]].copy()
df_year = df_year.drop_duplicates(subset=[country_col])
merged = world.merge(
df_year,
left_on="iso3",
right_on=country_col,
how="left",
)
fig, ax = plt.subplots(1, 1, figsize=(18, 10))
# Plot countries with missing data
missing = merged[merged[value_col].isna()]
missing.plot(ax=ax, color=missing_color, edgecolor="white", linewidth=0.3)
# Plot countries with data
present = merged[~merged[value_col].isna()]
present.plot(
column=value_col,
ax=ax,
cmap=cmap,
edgecolor="white",
linewidth=0.3,
vmin=vmin,
vmax=vmax,
legend=True,
legend_kwds={
"label": legend_label,
"orientation": "horizontal",
"fraction": 0.03,
"pad": 0.04,
"shrink": 0.6,
},
)
ax.set_title(f"{title}\n(Year: {year})", fontsize=16, fontweight="bold")
ax.axis("off")
# Coverage note
n_covered = len(present)
n_total = len(world)
ax.text(
0.01, 0.02,
f"Countries with data: {n_covered}/{n_total} | "
f"Missing: {n_total - n_covered}",
transform=ax.transAxes,
fontsize=9,
color="grey",
)
plt.tight_layout()
if output_file:
plt.savefig(output_file, dpi=150, bbox_inches="tight")
print(f"Map saved to {output_file}")
plt.show()
Example 1: Under-5 Mortality Trends in Sub-Saharan Africa (2000–2019)
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
# ------------------------------------------------------------------
# WHO GHO indicator: MDG_0000000007 = Under-5 mortality rate (per 1000 live births)
# ------------------------------------------------------------------
# Sub-Saharan African country ISO-3 codes (representative sample)
SSA_COUNTRIES = [
"NGA", "ETH", "COD", "TZA", "KEN", "UGA", "MOZ", "GHA",
"ZMB", "ZWE", "SEN", "MLI", "NER", "BFA", "CMR", "AGO",
"MWI", "RWA", "MDG", "SOM",
]
YEARS = list(range(2000, 2020))
print("Fetching under-5 mortality data from WHO GHO...")
u5mr_df = get_who_indicator(
indicator_code="MDG_0000000007",
countries=SSA_COUNTRIES,
years=YEARS,
verbose=True,
)
# Filter to both-sex estimates
u5mr_df = u5mr_df[u5mr_df["sex"].isin(["", "BTSX", None])].copy()
u5mr_df = u5mr_df.dropna(subset=["value", "year"])
u5mr_df["year"] = u5mr_df["year"].astype(int)
print(f"\nData shape: {u5mr_df.shape}")
print(u5mr_df.groupby("country_code")["value"].describe().round(1))
# ------------------------------------------------------------------
# Time series plot
# ------------------------------------------------------------------
fig, axes = plt.subplots(1, 2, figsize=(16, 7))
# Individual country trends
pivot = u5mr_df.pivot_table(
index="year", columns="country_code", values="value", aggfunc="mean"
)
for country in pivot.columns:
axes[0].plot(pivot.index, pivot[country], alpha=0.4, linewidth=1.0, color="steelblue")
# Regional mean
regional_mean = u5mr_df.groupby("year")["value"].mean()
regional_median = u5mr_df.groupby("year")["value"].median()
axes[0].plot(regional_mean.index, regional_mean.values,
color="navy", linewidth=3, label="Regional mean")
axes[0].plot(regional_median.index, regional_median.values,
color="darkred", linewidth=2, linestyle="--", label="Regional median")
axes[0].set_xlabel("Year")
axes[0].set_ylabel("Deaths per 1,000 live births")
axes[0].set_title("Under-5 Mortality Rate\nSub-Saharan Africa (2000–2019)")
axes[0].legend()
axes[0].grid(alpha=0.3)
# Percent reduction per country
start_val = u5mr_df[u5mr_df["year"] == 2000].set_index("country_code")["value"]
end_val = u5mr_df[u5mr_df["year"] == 2019].set_index("country_code")["value"]
common_idx = start_val.index.intersection(end_val.index)
pct_change = ((end_val[common_idx] - start_val[common_idx]) / start_val[common_idx] * 100)
pct_change = pct_change.sort_values()
colors = ["#d73027" if x > -30 else "#4dac26" for x in pct_change.values]
axes[1].barh(range(len(pct_change)), pct_change.values, color=colors, alpha=0.8)
axes[1].set_yticks(range(len(pct_change)))
axes[1].set_yticklabels(pct_change.index, fontsize=8)
axes[1].axvline(0, color="black", linewidth=0.8)
axes[1].set_xlabel("% Change in U5MR (2000–2019)")
axes[1].set_title("Percent Reduction in Under-5 Mortality\n(Green = >30% reduction)")
axes[1].grid(axis="x", alpha=0.3)
plt.tight_layout()
plt.savefig("ssa_u5mr_trends.png", dpi=150, bbox_inches="tight")
plt.show()
# ------------------------------------------------------------------
# Choropleth: U5MR in 2019
# ------------------------------------------------------------------
plot_global_health_map(
df=u5mr_df,
value_col="value",
year=2019,
title="Under-5 Mortality Rate",
country_col="country_code",
cmap="YlOrRd",
legend_label="Deaths per 1,000 live births",
output_file="u5mr_2019_map.png",
)
Example 2: DALYs from NCDs by WHO Region with Age-Standardization
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
# ------------------------------------------------------------------
# Simulate GBD-style NCD DALY data by WHO region
# (Replace with download_gbd_data() call when credentials are available)
# ------------------------------------------------------------------
who_regions = {
"AFR": "African Region",
"AMR": "Region of the Americas",
"SEAR": "South-East Asia Region",
"EUR": "European Region",
"EMR": "Eastern Mediterranean Region",
"WPR": "Western Pacific Region",
}
ncd_causes = [
"Cardiovascular diseases",
"Diabetes mellitus",
"Chronic respiratory diseases",
"Neoplasms",
]
age_groups = WHO_WORLD_STANDARD_POPULATION["age_group"].tolist()
age_weights = WHO_WORLD_STANDARD_POPULATION["weight"].values
np.random.seed(2024)
records = []
for region_code, region_name in who_regions.items():
for cause in ncd_causes:
for year in [2000, 2010, 2019]:
for age_group in age_groups:
base_yll = np.random.exponential(200)
base_yld = np.random.exponential(150)
records.append({
"region_code": region_code,
"region_name": region_name,
"cause": cause,
"year": year,
"age_group": age_group,
"yll": base_yll,
"yld": base_yld,
"population": np.random.randint(50_000, 5_000_000),
})
gbd_df = pd.DataFrame(records)
# ------------------------------------------------------------------
# Compute DALYs
# ------------------------------------------------------------------
daly_df = compute_dalys(
yll_series=gbd_df["yll"],
yld_series=gbd_df["yld"],
index_cols=gbd_df[["region_code", "region_name", "cause", "year",
"age_group", "population"]],
)
daly_df["daly_rate"] = daly_df["daly"] / daly_df["population"]
# ------------------------------------------------------------------
# Age-standardize DALYs by region, cause, and year
# ------------------------------------------------------------------
asr_records = []
for (region_code, cause, year), group in daly_df.groupby(
["region_code", "cause", "year"]
):
group = group.set_index("age_group")
# Align age groups
age_order = WHO_WORLD_STANDARD_POPULATION["age_group"].tolist()
crude_rates = group.reindex(age_order)["daly_rate"].fillna(0)
weights = pd.Series(age_weights, index=age_order)
asr_result = age_standardize(
crude_rates=crude_rates,
age_weights=weights,
per_unit=100_000,
)
asr_records.append({
"region_code": region_code,
"region_name": who_regions[region_code],
"cause": cause,
"year": year,
**asr_result,
})
asr_df = pd.DataFrame(asr_records)
print("Age-standardized DALY rates by region and cause (2019):")
print(
asr_df[asr_df["year"] == 2019]
.pivot_table(index="region_name", columns="cause",
values="age_standardized_rate")
.round(1)
.to_string()
)
# ------------------------------------------------------------------
# Inequality analysis: concentration index across regions
# ------------------------------------------------------------------
# Use GDP per capita proxy as wealth rank
gdp_rank = {
"AFR": 1, "SEAR": 2, "EMR": 3, "AMR": 4, "WPR": 5, "EUR": 6
}
asr_df_2019 = asr_df[asr_df["year"] == 2019].copy()
asr_df_2019["wealth_rank"] = asr_df_2019["region_code"].map(gdp_rank)
for cause in ncd_causes:
sub = asr_df_2019[asr_df_2019["cause"] == cause]
if len(sub) >= 3:
ci_result = compute_concentration_index(
health_outcome=sub["age_standardized_rate"],
wealth_rank=sub["wealth_rank"],
)
direction = "pro-rich" if ci_result["concentration_index"] > 0 else "pro-poor"
print(f"\n{cause}: CI={ci_result['concentration_index']:.3f} ({direction})")
# ------------------------------------------------------------------
# Visualization: grouped bar chart by WHO region for 2019
# ------------------------------------------------------------------
fig, axes = plt.subplots(1, 2, figsize=(16, 7))
pivot_2019 = asr_df_2019.pivot_table(
index="region_code",
columns="cause",
values="age_standardized_rate",
)
pivot_2019.plot(
kind="bar",
ax=axes[0],
colormap="tab10",
edgecolor="white",
alpha=0.85,
)
axes[0].set_xlabel("WHO Region")
axes[0].set_ylabel("Age-Standardized DALY Rate (per 100k)")
axes[0].set_title("NCD Burden by WHO Region (2019)\nAge-Standardized DALY Rates")
axes[0].legend(title="Cause", bbox_to_anchor=(1.0, 1.0), fontsize=8)
axes[0].tick_params(axis="x", rotation=30)
axes[0].grid(axis="y", alpha=0.3)
# Trend lines: total NCD DALYs over time
trend = asr_df.groupby(["region_code", "year"])["age_standardized_rate"].sum().reset_index()
for region in trend["region_code"].unique():
sub = trend[trend["region_code"] == region]
axes[1].plot(
sub["year"], sub["age_standardized_rate"],
marker="o", linewidth=2, label=region,
)
axes[1].set_xlabel("Year")
axes[1].set_ylabel("Total Age-Standardized DALY Rate (per 100k)")
axes[1].set_title("Total NCD DALY Trend by WHO Region\n(2000–2019)")
axes[1].legend(title="Region", fontsize=9)
axes[1].grid(alpha=0.3)
plt.suptitle("Global NCD Burden Analysis — IHME GBD + Age-Standardization", fontsize=13)
plt.tight_layout()
plt.savefig("ncd_daly_burden.png", dpi=150, bbox_inches="tight")
plt.show()
Notes and Best Practices
- WHO GHO API: Most indicators are publicly accessible without authentication. Pass
WHO_API_KEYenv var for rate-limit exemptions. Browse indicator codes athttps://ghoapi.azureedge.net/api/Indicator. - GBD data: For large-scale GBD downloads, use the IHME GHD Results Tool at
https://ghdx.healthdata.org/gbd-results. Save the CSV and pass viafallback_csv_path. - Our World in Data: OWID health datasets can be fetched directly as CSVs, e.g.,
pd.read_csv("https://ourworldindata.org/grapher/child-mortality.csv?tab=table"). - Age standardization: Always report the reference population used (WHO World Standard 2000-2025 vs. Segi world population). Results are not comparable across different standards.
- Concentration index: Interpret with caution for non-continuous wealth measures. Wagstaff's normalized CI is preferred for binary outcomes.
- Mapping:
geopandas.datasets.get_path('naturalearth_lowres')is deprecated in newer geopandas. Usegeodatasets.get_path('naturalearth.land')or download from Natural Earth directly for production code. - Causality: Cross-country observational health data cannot establish causal effects. Use DAG-based confounder control or difference-in-differences designs for causal inference.
- DALY weights: Disability weights for YLD calculations should come from the GBD study directly; do not use outdated WHO 1990 weights.