Analyze Child Malnutrition in Nigeria

Clean real DHS data and apply Bayesian probability to Nigerian nutrition data.

Introduction

30 Second Summary

When governments decide where to send food aid, they rely on data. Specifically, they need numbers that summarize thousands of health measurements into clear regional patterns.

In this project, you will take real Nigerian health survey data, clean it from scratch, compute core statistics by hand, and use Bayes' theorem to answer the question: if a region has high child stunting, what is the probability it is in Northern Nigeria?

What You'll Build

A Jupyter Notebook that processes raw Demographic and Health Survey nutrition data through a complete analysis pipeline, from data cleaning to Bayesian inference.

By the end of this project, you'll have:

  • A cleaned and merged dataset combining anthropometry and nutrition indicator data for Nigeria into an analysis-ready format.
  • Custom Python functions computing mean, median, mode, variance, standard deviation, skewness, and kurtosis without using built-in statistical shortcuts.
  • A complete Bayes' theorem calculation with joint, marginal, and conditional probability tables answering a real public health question about regional malnutrition.
  • Secret Mission: Analyze the wealth-stunting relationship across quintiles and formally test for statistical independence.

Are there any prerequisites?

You should have basic Python familiarity (variables, loops, functions). You also need Anaconda and VS Code installed on your machine.

Anaconda ships with Jupyter Notebook, pandas, NumPy, and matplotlib pre-installed, so no additional package installs are needed.

Before We Start

Before We Start...

Set Up Your Project and Load the Data

Before you can analyze malnutrition data, you need a workspace to hold your files and a tool to run your code interactively. Jupyter Notebook lets you write Python in small cells and see results instantly, which is perfect for exploratory data analysis.

In this step, you will create a project folder, launch Jupyter, download real nutrition survey data from Nigeria, and load it into pandas DataFrames for inspection.

In this step, get ready to:
  • Create a project folder and launch Jupyter Notebook.
  • Download two CSV files from the Humanitarian Data Exchange.
  • Import your libraries, load the data, and inspect it.
Create your project folder and launch Jupyter Notebook

You need a dedicated folder to keep your notebook and data files together. You will create this folder on your Desktop so it is easy to find.

  • Press Cmd+Space (macOS) or the Windows key (Windows) to open your search bar.
  • Type Terminal (macOS) or Anaconda Prompt (Windows) and press Enter to open it.

Why Anaconda Prompt on Windows?

Anaconda Prompt is a terminal that already has Python, Jupyter, and all Anaconda packages on its PATH. On macOS, the default Terminal works because Anaconda adds itself to your shell during installation.

  • Move to your Desktop by running this command:
cd ~/Desktop
  • Create your project folder and move into it by running these commands:
mkdir nigeria-nutrition-analysis
cd nigeria-nutrition-analysis
  • Confirm the folder exists by running:
ls

The folder will be empty for now. That is expected.

  • Launch Jupyter Notebook by running:
jupyter notebook

Your default browser opens to the Jupyter dashboard, showing files in the current directory.

  • Click the New Notebook button in the top-right area of the dashboard.
  • Select Python 3 as the kernel if prompted.
  • Click the notebook name at the top of the page (it will say Untitled) to open the rename dialog.
  • Type nigeria_nutrition_analysis and press Enter to rename it.

Your notebook is now saved as nigeria_nutrition_analysis.ipynb inside your project folder.

Browser did not open?

If the browser does not open automatically, look in the terminal output for a URL starting with http://localhost:8888. Copy and paste that full URL into your browser manually.

If you see a "command not found" error, make sure you opened Anaconda Prompt (Windows) or that Anaconda is on your PATH (macOS). Help me fix my Jupyter launch error.

Download the nutrition data from HDX

The data you will analyze comes from Nigeria's Demographic and Health Survey, hosted on the Humanitarian Data Exchange (HDX). HDX is a free, open platform for sharing humanitarian datasets. No account is required to download.

  • Open this link in your browser: https://data.humdata.org/dataset/dhs-data-for-nigeria
  • Scroll down to the list of data files (the Data and Resources section).
  • Find the file named Anthropometry Data for Nigeria (CSV, 21.7K).
  • Click the Download button next to it to save the file to your computer.
  • Find the file named Select Nutrition Indicators Data for Nigeria (CSV, 10.4K).
  • Click Download next to it as well.
  • Move both downloaded files into your nigeria-nutrition-analysis folder on the Desktop.
  • Rename the anthropometry file to anthropometry_nigeria.csv and the nutrition indicators file to nutrition_indicators_nigeria.csv.

Why rename the files?

HDX download filenames can be long or include version numbers. Renaming to short, consistent names means your Python code can reference them reliably without worrying about filename changes.

Don't see the files on the HDX page?

Make sure you are on the correct dataset page for "DHS Data for Nigeria". The page lists multiple CSV resources grouped by topic. Scroll through the full list to locate the Anthropometry and Nutrition Indicators files.

Help me find the right files on the HDX page.

Import libraries and load the data

With your notebook open and CSV files in place, you can now bring in the Python libraries and load the data. You will use pandas for data manipulation, NumPy for numerical computation, and matplotlib for visualization.

  • In the first cell of your notebook, import your libraries by typing this code and pressing Shift+Enter to run it:
# Cell 1: Imports and Setup
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt

What does this code do?

Each import statement loads a library and gives it a short alias. pd is the standard alias for pandas, np for NumPy, and plt for matplotlib's pyplot module.

These three libraries come pre-installed with Anaconda, so no additional installation is needed.

  • In a new cell, load both CSV files and print basic info by typing this code and pressing Shift+Enter:
# Cell 2: Load Data
anthropometry = pd.read_csv('anthropometry_nigeria.csv')
nutrition = pd.read_csv('nutrition_indicators_nigeria.csv')

print("Anthropometry shape:", anthropometry.shape)
print("Nutrition shape:", nutrition.shape)
print("\nAnthropometry columns:", anthropometry.columns.tolist())
print("\nFirst few rows:")
anthropometry.head()

What does this code do?

pd.read_csv() reads each CSV file into a DataFrame, which is pandas' table-like data structure.

.shape tells you how many rows and columns each dataset has. .columns.tolist() lists all column names so you know what data is available. .head() displays the first 5 rows as a preview.

You should see output showing the shape of both DataFrames and a table of the first few rows of the anthropometry data.

  • In a new cell, inspect the data structure in more detail by running:
# Cell 3: Inspect Data Structure
print(anthropometry.info())
print("\nUnique indicators:", anthropometry['Indicator'].nunique())
print("\nSample indicators:")
print(anthropometry['Indicator'].unique()[:10])
print("\nCharacteristic categories:")
print(anthropometry['CharacteristicCategory'].unique())
print("\nSurvey years:", anthropometry['SurveyYear'].unique())

What does this code do?

.info() shows the data types and non-null counts for every column, helping you spot missing data early.

.nunique() counts how many distinct indicators exist in the dataset. .unique() shows the actual category values so you can understand how the data is organized (by region, wealth quintile, urban/rural, etc.).

You should see a list of characteristic categories like "Total", "Residence", "Region", and "Wealth quintile". These represent the different ways the DHS data is broken down.

Seeing a FileNotFoundError?

This means Python cannot find the CSV files. Check that both anthropometry_nigeria.csv and nutrition_indicators_nigeria.csv are in the same folder where you launched Jupyter.

If you moved the files after launching Jupyter, they need to be in the nigeria-nutrition-analysis folder on your Desktop. Also confirm that the filenames match exactly, including the underscore and lowercase letters.

Help me fix my FileNotFoundError when loading CSVs in Jupyter.

✔️ Awesome, I've got everything!

Great. Double check that you have saved your notebook by pressing Cmd+S (macOS) or Ctrl+S (Windows).

ⓧ I'd like to double check the full code

Here is the complete code for all three cells so far:

# Cell 1: Imports and Setup
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt

# Cell 2: Load Data
anthropometry = pd.read_csv('anthropometry_nigeria.csv')
nutrition = pd.read_csv('nutrition_indicators_nigeria.csv')

print("Anthropometry shape:", anthropometry.shape)
print("Nutrition shape:", nutrition.shape)
print("\nAnthropometry columns:", anthropometry.columns.tolist())
print("\nFirst few rows:")
anthropometry.head()

# Cell 3: Inspect Data Structure
print(anthropometry.info())
print("\nUnique indicators:", anthropometry['Indicator'].nunique())
print("\nSample indicators:")
print(anthropometry['Indicator'].unique()[:10])
print("\nCharacteristic categories:")
print(anthropometry['CharacteristicCategory'].unique())
print("\nSurvey years:", anthropometry['SurveyYear'].unique())

Your workspace is set up and you have real Nigerian health survey data loaded and ready to explore. Next up, you will clean and reshape this data to focus on the malnutrition indicators that matter for your analysis.

Clean and Reshape the Data

You have two raw DataFrames loaded and you know what columns they contain. The problem is that this data is messy. It contains multiple survey years, several indicator types you do not need, potential missing values, and a "long" row-per-indicator format that makes comparison difficult.

In this step, you will filter, clean, and reshape the data into a single wide-format DataFrame where each row is a characteristic (like a region or wealth quintile) and each column is a nutrition indicator. This is the standard workflow for preparing real-world survey data for analysis.

In this step, get ready to:
  • Filter to the most recent survey year and key nutrition indicators.
  • Handle missing values, standardize labels, and remove duplicates.
  • Reshape from long format to wide format for analysis.
Filter to the most recent survey year and key indicators

The HDX dataset contains data from multiple survey years spanning 1990 to 2024. You only want the most recent year for your analysis. The dataset also contains dozens of indicator types, but you only need three: stunting, wasting, and underweight.

What are stunting, wasting, and underweight?

These are the three standard measures of child malnutrition used by the WHO and DHS surveys. Stunting measures chronic malnutrition (low height-for-age). Wasting measures acute malnutrition (low weight-for-height). Underweight combines both (low weight-for-age).

In the DHS indicator system, these are identified by codes: CN_NUTS_C_HA2 for stunting, CN_NUTS_C_WH2 for wasting, and CN_NUTS_C_WA2 for underweight.

  • Create a new cell in your notebook (Cell 4) and add the following filtering code:
# Cell 4: Filter and Clean
latest_year = anthropometry['SurveyYear'].max()
anthro_recent = anthropometry[anthropometry['SurveyYear'] == latest_year].copy()
nutr_recent = nutrition[nutrition['SurveyYear'] == latest_year].copy()

stunting_ids = anthro_recent['IndicatorId'].str.contains('CN_NUTS_C_HA2', na=False)
wasting_ids = anthro_recent['IndicatorId'].str.contains('CN_NUTS_C_WH2', na=False)
underweight_ids = anthro_recent['IndicatorId'].str.contains('CN_NUTS_C_WA2', na=False)

nutrition_data = anthro_recent[stunting_ids | wasting_ids | underweight_ids].copy()
print(f"Filtered to {len(nutrition_data)} rows for stunting/wasting/underweight")

What does this code do?

  • anthropometry['SurveyYear'].max() finds the most recent survey year in the dataset automatically.
  • .str.contains('CN_NUTS_C_HA2', na=False) creates a boolean mask that is True for rows where the IndicatorId contains the stunting code.
  • The | operator combines the three masks with OR logic, keeping any row that matches stunting, wasting, or underweight.
  • .copy() creates independent DataFrames so edits in later cells do not modify the originals.
  • Run the cell by pressing Shift+Enter.
  • You should see a print statement showing how many rows remain after filtering (something like "Filtered to X rows for stunting/wasting/underweight").

Seeing 0 rows after filtering?

If nutrition_data is empty, the indicator IDs in your data might not match the filter strings exactly. Run anthro_recent['IndicatorId'].unique() in a new cell to inspect what indicator codes your dataset actually contains.

help me debug why my indicator filter returned zero rows

Handle missing values and standardize labels

Real-world data often has missing or non-numeric entries hiding in columns that should be numbers. The Value column might contain blank cells or text placeholders that would break your calculations. You also need consistent text labels for grouping later.

  • Create a new cell (Cell 5) and add this cleaning code:
# Cell 5: Handle Missing Values and Type Conversion
nutrition_data['Value'] = pd.to_numeric(nutrition_data['Value'], errors='coerce')
print(f"Missing values in Value column: {nutrition_data['Value'].isna().sum()}")
nutrition_data = nutrition_data.dropna(subset=['Value'])

nutrition_data['CharacteristicLabel'] = nutrition_data['CharacteristicLabel'].str.strip()
nutrition_data['CharacteristicCategory'] = nutrition_data['CharacteristicCategory'].str.strip()

duplicates = nutrition_data.duplicated(subset=['IndicatorId', 'CharacteristicCategory', 'CharacteristicLabel'])
print(f"Duplicate rows: {duplicates.sum()}")
nutrition_data = nutrition_data.drop_duplicates(subset=['IndicatorId', 'CharacteristicCategory', 'CharacteristicLabel'])

print(f"\nCleaned dataset shape: {nutrition_data.shape}")
print(f"Characteristic categories: {nutrition_data['CharacteristicCategory'].unique()}")

What does this code do?

  • pd.to_numeric(..., errors='coerce') converts the Value column to numbers. Any text or blank entries become NaN instead of crashing your script.
  • .dropna(subset=['Value']) removes any rows where the value could not be converted to a number.
  • .str.strip() removes leading and trailing whitespace from text labels. Without this, "North West" and "North West " would be treated as different categories.
  • .drop_duplicates() removes rows that have the same indicator, category, and label combination. This prevents double-counting during analysis.
  • Run the cell with Shift+Enter.
  • You should see the count of missing values, the count of duplicates, the cleaned dataset shape, and a list of characteristic categories (like "Total", "Residence", "Region", "Wealth quintile", etc.).

Getting a SettingWithCopyWarning?

If you see a warning about setting values on a copy of a slice, make sure you used .copy() when creating nutrition_data in Cell 4. The .copy() method creates an independent DataFrame that you can safely modify.

help me fix this SettingWithCopyWarning in my pandas code

Reshape from long to wide format

Right now, each row represents one indicator for one characteristic. A single region like "North East" has three separate rows for stunting, wasting, and underweight. For analysis, you want one row per characteristic with columns for each indicator. This is called pivoting from long to wide format.

  • Create a new cell (Cell 6) and add this reshaping code:
# Cell 6: Reshape Data - Create Indicator Column
def label_indicator(indicator_id):
    if 'HA2' in str(indicator_id):
        return 'Stunting'
    elif 'WH2' in str(indicator_id):
        return 'Wasting'
    elif 'WA2' in str(indicator_id):
        return 'Underweight'
    return 'Other'

nutrition_data['NutritionIndicator'] = nutrition_data['IndicatorId'].apply(label_indicator)

nutrition_wide = nutrition_data.pivot_table(
    index=['CharacteristicCategory', 'CharacteristicLabel'],
    columns='NutritionIndicator',
    values='Value',
    aggfunc='first'
).reset_index()

print("Reshaped data:")
print(nutrition_wide.head(10))
print(f"\nFinal shape: {nutrition_wide.shape}")
print(nutrition_wide.describe())

What does this code do?

  • label_indicator() is a helper function that maps the long DHS indicator codes to readable names like "Stunting", "Wasting", and "Underweight".
  • .apply(label_indicator) runs this function on every row to create a new NutritionIndicator column.
  • .pivot_table() reshapes the data so each unique combination of category and label becomes one row, with separate columns for Stunting, Wasting, and Underweight values.
  • .describe() gives you a quick summary (count, mean, std, min, max) to confirm the data looks reasonable.
  • Run the cell with Shift+Enter.
  • You should see a table where each row has a CharacteristicCategory, a CharacteristicLabel, and numeric columns for Stunting, Wasting, and Underweight.
  • The .describe() output should show reasonable percentage values (roughly between 0 and 60) for each indicator.

Seeing all NaN values in your pivot table?

If the Stunting, Wasting, or Underweight columns are entirely NaN, the NutritionIndicator column might not have been created correctly. Run nutrition_data['NutritionIndicator'].unique() to check that it contains "Stunting", "Wasting", and "Underweight" (not "Other").

help me debug why my pivot table has NaN values

✔️ Awesome, I've got everything!

Great. Make sure you have saved your notebook with all six cells before moving on.

ⓧ I'd like to double check the full code

Here is the full notebook up to this point. Compare your cells against this:

# Cell 1: Imports and Setup
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt

# Cell 2: Load Data
anthropometry = pd.read_csv('anthropometry_nigeria.csv')
nutrition = pd.read_csv('nutrition_indicators_nigeria.csv')

print("Anthropometry shape:", anthropometry.shape)
print("Nutrition shape:", nutrition.shape)
print("\nAnthropometry columns:", anthropometry.columns.tolist())
print("\nFirst few rows:")
anthropometry.head()

# Cell 3: Inspect Data Structure
print(anthropometry.info())
print("\nUnique indicators:", anthropometry['Indicator'].nunique())
print("\nSample indicators:")
print(anthropometry['Indicator'].unique()[:10])
print("\nCharacteristic categories:")
print(anthropometry['CharacteristicCategory'].unique())
print("\nSurvey years:", anthropometry['SurveyYear'].unique())

# Cell 4: Filter and Clean
latest_year = anthropometry['SurveyYear'].max()
anthro_recent = anthropometry[anthropometry['SurveyYear'] == latest_year].copy()
nutr_recent = nutrition[nutrition['SurveyYear'] == latest_year].copy()

stunting_ids = anthro_recent['IndicatorId'].str.contains('CN_NUTS_C_HA2', na=False)
wasting_ids = anthro_recent['IndicatorId'].str.contains('CN_NUTS_C_WH2', na=False)
underweight_ids = anthro_recent['IndicatorId'].str.contains('CN_NUTS_C_WA2', na=False)

nutrition_data = anthro_recent[stunting_ids | wasting_ids | underweight_ids].copy()
print(f"Filtered to {len(nutrition_data)} rows for stunting/wasting/underweight")

# Cell 5: Handle Missing Values and Type Conversion
nutrition_data['Value'] = pd.to_numeric(nutrition_data['Value'], errors='coerce')
print(f"Missing values in Value column: {nutrition_data['Value'].isna().sum()}")
nutrition_data = nutrition_data.dropna(subset=['Value'])

nutrition_data['CharacteristicLabel'] = nutrition_data['CharacteristicLabel'].str.strip()
nutrition_data['CharacteristicCategory'] = nutrition_data['CharacteristicCategory'].str.strip()

duplicates = nutrition_data.duplicated(subset=['IndicatorId', 'CharacteristicCategory', 'CharacteristicLabel'])
print(f"Duplicate rows: {duplicates.sum()}")
nutrition_data = nutrition_data.drop_duplicates(subset=['IndicatorId', 'CharacteristicCategory', 'CharacteristicLabel'])

print(f"\nCleaned dataset shape: {nutrition_data.shape}")
print(f"Characteristic categories: {nutrition_data['CharacteristicCategory'].unique()}")

# Cell 6: Reshape Data - Create Indicator Column
def label_indicator(indicator_id):
    if 'HA2' in str(indicator_id):
        return 'Stunting'
    elif 'WH2' in str(indicator_id):
        return 'Wasting'
    elif 'WA2' in str(indicator_id):
        return 'Underweight'
    return 'Other'

nutrition_data['NutritionIndicator'] = nutrition_data['IndicatorId'].apply(label_indicator)

nutrition_wide = nutrition_data.pivot_table(
    index=['CharacteristicCategory', 'CharacteristicLabel'],
    columns='NutritionIndicator',
    values='Value',
    aggfunc='first'
).reset_index()

print("Reshaped data:")
print(nutrition_wide.head(10))
print(f"\nFinal shape: {nutrition_wide.shape}")
print(nutrition_wide.describe())

Your data is now clean, filtered, and in a format that makes analysis straightforward. Next up, you will write your own statistical functions from scratch to compute descriptive statistics on this data.

Compute Descriptive Statistics from Scratch

Your data is clean and reshaped. Now the real question is: what story do these numbers tell about child malnutrition in Nigeria?

In this step, you will write every core statistical function from scratch using only NumPy. Building these by hand forces you to understand the math behind each measure, which is exactly what employers expect in data science interviews.

In this step, get ready to:
  • Build custom functions for mean, median, mode, variance, standard deviation, range, and IQR.
  • Compute skewness, kurtosis, and percentiles to understand distribution shape.
  • Compare mean vs median and explain which better summarizes skewed nutrition data.
Write custom functions for central tendency and spread

Descriptive statistics break into two families: measures of center (where the data clusters) and measures of spread (how far apart the values are). You will write functions for both.

  • Create a new cell in your notebook (Cell 7). Start by extracting the stunting values into a NumPy array, then define your central tendency functions by pasting this code:
# Cell 7: Descriptive Statistics from Scratch
# Get stunting values for analysis
stunting_values = nutrition_wide['Stunting'].dropna().values

# Mean
def compute_mean(data):
    return np.sum(data) / len(data)

# Median
def compute_median(data):
    sorted_data = np.sort(data)
    n = len(sorted_data)
    if n % 2 == 0:
        return (sorted_data[n//2 - 1] + sorted_data[n//2]) / 2
    else:
        return sorted_data[n//2]

# Mode
def compute_mode(data):
    values, counts = np.unique(np.round(data, 1), return_counts=True)
    max_count = np.max(counts)
    modes = values[counts == max_count]
    return modes

What does this code do?

stunting_values extracts all non-null stunting percentages as a raw NumPy array. This is the data you will feed into every function.

  • compute_mean sums all values and divides by the count. The simplest measure of center.
  • compute_median sorts the data and picks the middle value. If the count is even, it averages the two middle values.
  • compute_mode rounds values to one decimal, counts occurrences, and returns the most frequent value(s).
  • In the same cell (Cell 7), add the spread functions and the print statements that apply everything by pasting this code directly below:
# Variance (population and sample)
def compute_variance(data, sample=True):
    mean = compute_mean(data)
    squared_diffs = (data - mean) ** 2
    if sample:
        return np.sum(squared_diffs) / (len(data) - 1)
    else:
        return np.sum(squared_diffs) / len(data)

# Standard Deviation
def compute_std(data, sample=True):
    return np.sqrt(compute_variance(data, sample))

# Range
def compute_range(data):
    return np.max(data) - np.min(data)

# IQR
def compute_iqr(data):
    sorted_data = np.sort(data)
    q1 = np.percentile(sorted_data, 25)
    q3 = np.percentile(sorted_data, 75)
    return q3 - q1

# Apply all functions
print("=" * 50)
print("DESCRIPTIVE STATISTICS: Stunting Prevalence (%)")
print("=" * 50)
print(f"Mean:                {compute_mean(stunting_values):.2f}")
print(f"Median:              {compute_median(stunting_values):.2f}")
print(f"Mode:                {compute_mode(stunting_values)}")
print(f"Sample Variance:     {compute_variance(stunting_values):.2f}")
print(f"Sample Std Dev:      {compute_std(stunting_values):.2f}")
print(f"Range:               {compute_range(stunting_values):.2f}")
print(f"IQR:                 {compute_iqr(stunting_values):.2f}")

What do the spread functions measure?

  • compute_variance measures how far values deviate from the mean on average (squared). The sample=True parameter divides by n-1 (Bessel's correction) to avoid underestimating spread from a sample.
  • compute_std is the square root of variance. It brings units back to the same scale as the original data (percentage points).
  • compute_range is the gap between the highest and lowest values. Simple but sensitive to outliers.
  • compute_iqr is the range of the middle 50% of values (75th percentile minus 25th). It ignores extremes entirely.
  • Run Cell 7 by pressing **Shift+Enter**.
  • You should see a formatted table showing all seven statistics with numeric values. Compare the mean and standard deviation to what .describe() showed at the end of Cell 6. They should match.

Seeing a NameError or KeyError?

  • If you see NameError: name 'nutrition_wide' is not defined, make sure you ran all previous cells (Cells 1 through 6) first. Jupyter needs each cell executed in order.
  • If you see KeyError: 'Stunting', check that your pivot in Cell 6 produced a column named Stunting. Run nutrition_wide.columns.tolist() in a scratch cell to verify.

Help me debug my descriptive statistics cell

Why build these from scratch instead of using pandas?

Pandas has .mean(), .median(), and .std() built in. But calling a function you do not understand is just pushing buttons.

Writing the formulas yourself means you can explain what sample variance actually measures in an interview. You also catch bugs faster when you know what the function should produce.

✔️ Awesome, I've got everything!

Great. Make sure you have saved your notebook before moving on.

ⓧ I'd like to double check the full code

Here is the complete Cell 7:

# Cell 7: Descriptive Statistics from Scratch
# Get stunting values for analysis
stunting_values = nutrition_wide['Stunting'].dropna().values

# Mean
def compute_mean(data):
    return np.sum(data) / len(data)

# Median
def compute_median(data):
    sorted_data = np.sort(data)
    n = len(sorted_data)
    if n % 2 == 0:
        return (sorted_data[n//2 - 1] + sorted_data[n//2]) / 2
    else:
        return sorted_data[n//2]

# Mode
def compute_mode(data):
    values, counts = np.unique(np.round(data, 1), return_counts=True)
    max_count = np.max(counts)
    modes = values[counts == max_count]
    return modes

# Variance (population and sample)
def compute_variance(data, sample=True):
    mean = compute_mean(data)
    squared_diffs = (data - mean) ** 2
    if sample:
        return np.sum(squared_diffs) / (len(data) - 1)
    else:
        return np.sum(squared_diffs) / len(data)

# Standard Deviation
def compute_std(data, sample=True):
    return np.sqrt(compute_variance(data, sample))

# Range
def compute_range(data):
    return np.max(data) - np.min(data)

# IQR
def compute_iqr(data):
    sorted_data = np.sort(data)
    q1 = np.percentile(sorted_data, 25)
    q3 = np.percentile(sorted_data, 75)
    return q3 - q1

# Apply all functions
print("=" * 50)
print("DESCRIPTIVE STATISTICS: Stunting Prevalence (%)")
print("=" * 50)
print(f"Mean:                {compute_mean(stunting_values):.2f}")
print(f"Median:              {compute_median(stunting_values):.2f}")
print(f"Mode:                {compute_mode(stunting_values)}")
print(f"Sample Variance:     {compute_variance(stunting_values):.2f}")
print(f"Sample Std Dev:      {compute_std(stunting_values):.2f}")
print(f"Range:               {compute_range(stunting_values):.2f}")
print(f"IQR:                 {compute_iqr(stunting_values):.2f}")
Compute skewness, kurtosis, and percentiles

Center and spread tell you where data sits and how wide it stretches. But they do not tell you whether the data is lopsided or has extreme tails. That is the job of skewness and kurtosis.

  • Create a new cell (Cell 8). Define the skewness and kurtosis functions, then compute percentiles by pasting this code:
# Cell 8: Skewness and Kurtosis
def compute_skewness(data):
    n = len(data)
    mean = compute_mean(data)
    std = compute_std(data, sample=True)
    skew = (n / ((n-1) * (n-2))) * np.sum(((data - mean) / std) ** 3)
    return skew

def compute_kurtosis(data):
    n = len(data)
    mean = compute_mean(data)
    std = compute_std(data, sample=True)
    kurt = ((n * (n+1)) / ((n-1) * (n-2) * (n-3))) * np.sum(((data - mean) / std) ** 4)
    kurt -= (3 * (n-1)**2) / ((n-2) * (n-3))
    return kurt

print(f"\nSkewness:            {compute_skewness(stunting_values):.3f}")
print(f"Kurtosis (excess):   {compute_kurtosis(stunting_values):.3f}")

# Percentiles
for p in [25, 50, 75, 90]:
    val = np.percentile(stunting_values, p)
    print(f"{p}th percentile:     {val:.2f}%")

What do skewness and kurtosis tell you?

  • Skewness measures asymmetry. A positive value means the tail stretches to the right (a few very high stunting regions pull the distribution upward). A negative value means the tail stretches left.
  • Kurtosis (excess) measures tail heaviness. A positive value means more extreme values than a normal distribution would predict. Nigeria's North-South disparity creates exactly this pattern.
  • Percentiles split the data into ordered slices. The 25th percentile means 25% of breakdowns have stunting below that value. The 90th percentile flags the most extreme regions.
  • Run Cell 8 by pressing **Shift+Enter**.
  • You should see a positive skewness value (confirming the distribution is right-skewed due to high northern stunting rates) and four percentile values printed below it.

Getting a ZeroDivisionError?

The kurtosis formula divides by (n-3). If your dataset has 3 or fewer values, this will fail. Check that Cell 7 ran successfully and that stunting_values has more than 3 elements by running len(stunting_values) in a scratch cell.

Help me fix my skewness or kurtosis calculation

✔️ Awesome, I've got everything!

Great. Save your notebook before continuing.

ⓧ I'd like to double check the full code

Here is the complete Cell 8:

# Cell 8: Skewness and Kurtosis
def compute_skewness(data):
    n = len(data)
    mean = compute_mean(data)
    std = compute_std(data, sample=True)
    skew = (n / ((n-1) * (n-2))) * np.sum(((data - mean) / std) ** 3)
    return skew

def compute_kurtosis(data):
    n = len(data)
    mean = compute_mean(data)
    std = compute_std(data, sample=True)
    kurt = ((n * (n+1)) / ((n-1) * (n-2) * (n-3))) * np.sum(((data - mean) / std) ** 4)
    kurt -= (3 * (n-1)**2) / ((n-2) * (n-3))
    return kurt

print(f"\nSkewness:            {compute_skewness(stunting_values):.3f}")
print(f"Kurtosis (excess):   {compute_kurtosis(stunting_values):.3f}")

# Percentiles
for p in [25, 50, 75, 90]:
    val = np.percentile(stunting_values, p)
    print(f"{p}th percentile:     {val:.2f}%")
Compare mean vs median for skewed nutrition data

You just confirmed the stunting distribution is positively skewed. When data is skewed, the mean gets dragged toward the extreme tail. This makes it a misleading summary for policy decisions.

  • Create a new cell (Cell 9). Document the mean vs median comparison with a concrete numerical example by pasting this code:
# Cell 9: Mean vs Median Comparison
print("\n" + "=" * 50)
print("WHY MEDIAN BEATS MEAN FOR SKEWED DATA")
print("=" * 50)
mean_val = compute_mean(stunting_values)
median_val = compute_median(stunting_values)
print(f"Mean stunting:   {mean_val:.2f}%")
print(f"Median stunting: {median_val:.2f}%")
print(f"Difference:      {mean_val - median_val:.2f} percentage points")
print(f"\nThe mean is pulled {'upward' if mean_val > median_val else 'downward'} by extreme values.")
print("For policy decisions, the median better represents the 'typical' region.")

Why does this matter for real policy?

Nigeria's North West has stunting rates around 50%, while South East and South West sit around 15-20%. The mean gets pulled upward by these extreme northern values.

If a policymaker uses the mean to describe the "typical" region, they overstate the problem in the South and understate it in the North. The median gives a more honest picture of where most regions actually sit.

  • Run Cell 9 by pressing **Shift+Enter**.
  • You should see the mean and median printed side by side, a positive difference (confirming the mean is higher than the median), and the interpretation that the mean is pulled upward by extreme values.
  • Save your notebook.

Mean and median are the same value?

If your mean equals your median exactly, your data may not have loaded correctly in earlier cells. Run print(stunting_values) to check that you have a range of different values. If you only see one or two unique numbers, go back and verify Cell 4's indicator filtering worked correctly.

Help me understand why my mean and median are identical

✔️ Awesome, I've got everything!

Great. Your notebook now has all core descriptive statistics computed from scratch.

ⓧ I'd like to double check the full code

Here are all three cells added in this step (Cells 7, 8, and 9):

# Cell 7: Descriptive Statistics from Scratch
# Get stunting values for analysis
stunting_values = nutrition_wide['Stunting'].dropna().values

# Mean
def compute_mean(data):
    return np.sum(data) / len(data)

# Median
def compute_median(data):
    sorted_data = np.sort(data)
    n = len(sorted_data)
    if n % 2 == 0:
        return (sorted_data[n//2 - 1] + sorted_data[n//2]) / 2
    else:
        return sorted_data[n//2]

# Mode
def compute_mode(data):
    values, counts = np.unique(np.round(data, 1), return_counts=True)
    max_count = np.max(counts)
    modes = values[counts == max_count]
    return modes

# Variance (population and sample)
def compute_variance(data, sample=True):
    mean = compute_mean(data)
    squared_diffs = (data - mean) ** 2
    if sample:
        return np.sum(squared_diffs) / (len(data) - 1)
    else:
        return np.sum(squared_diffs) / len(data)

# Standard Deviation
def compute_std(data, sample=True):
    return np.sqrt(compute_variance(data, sample))

# Range
def compute_range(data):
    return np.max(data) - np.min(data)

# IQR
def compute_iqr(data):
    sorted_data = np.sort(data)
    q1 = np.percentile(sorted_data, 25)
    q3 = np.percentile(sorted_data, 75)
    return q3 - q1

# Apply all functions
print("=" * 50)
print("DESCRIPTIVE STATISTICS: Stunting Prevalence (%)")
print("=" * 50)
print(f"Mean:                {compute_mean(stunting_values):.2f}")
print(f"Median:              {compute_median(stunting_values):.2f}")
print(f"Mode:                {compute_mode(stunting_values)}")
print(f"Sample Variance:     {compute_variance(stunting_values):.2f}")
print(f"Sample Std Dev:      {compute_std(stunting_values):.2f}")
print(f"Range:               {compute_range(stunting_values):.2f}")
print(f"IQR:                 {compute_iqr(stunting_values):.2f}")
# Cell 8: Skewness and Kurtosis
def compute_skewness(data):
    n = len(data)
    mean = compute_mean(data)
    std = compute_std(data, sample=True)
    skew = (n / ((n-1) * (n-2))) * np.sum(((data - mean) / std) ** 3)
    return skew

def compute_kurtosis(data):
    n = len(data)
    mean = compute_mean(data)
    std = compute_std(data, sample=True)
    kurt = ((n * (n+1)) / ((n-1) * (n-2) * (n-3))) * np.sum(((data - mean) / std) ** 4)
    kurt -= (3 * (n-1)**2) / ((n-2) * (n-3))
    return kurt

print(f"\nSkewness:            {compute_skewness(stunting_values):.3f}")
print(f"Kurtosis (excess):   {compute_kurtosis(stunting_values):.3f}")

# Percentiles
for p in [25, 50, 75, 90]:
    val = np.percentile(stunting_values, p)
    print(f"{p}th percentile:     {val:.2f}%")
# Cell 9: Mean vs Median Comparison
print("\n" + "=" * 50)
print("WHY MEDIAN BEATS MEAN FOR SKEWED DATA")
print("=" * 50)
mean_val = compute_mean(stunting_values)
median_val = compute_median(stunting_values)
print(f"Mean stunting:   {mean_val:.2f}%")
print(f"Median stunting: {median_val:.2f}%")
print(f"Difference:      {mean_val - median_val:.2f} percentage points")
print(f"\nThe mean is pulled {'upward' if mean_val > median_val else 'downward'} by extreme values.")
print("For policy decisions, the median better represents the 'typical' region.")

You now have every core descriptive statistic computed from scratch and an understanding of why distribution shape matters for interpretation. Next up, you will turn these numbers into visualizations that make the North-South disparity impossible to miss.

Visualize Distribution Shape

You have a full set of descriptive statistics for Nigeria's stunting data. Numbers like skewness and IQR are powerful, but they become much easier to interpret when you can see the shape of the distribution.

In this step, you will create visualizations that make abstract statistics tangible. A histogram reveals how values cluster together. A box plot highlights outliers. An empirical CDF shows what fraction of observations fall below any given threshold.

In this step, get ready to:
  • Create a three-panel distribution figure with histogram, box plot, and empirical CDF.
  • Build a grouped bar chart comparing stunting, wasting, and underweight across Nigeria's six geopolitical zones.
  • Identify which regions exceed the WHO 30% "very high" stunting threshold.
Create the distribution shape figure

A single figure with three panels gives you three complementary views of the same data. The histogram shows how stunting values are distributed. The box plot highlights the median, IQR, and any outliers. The empirical CDF shows the cumulative probability at every value.

  • Create a new cell in your notebook (Cell 10) and add the following visualization code:
# Cell 10: Visualization - Distribution Shape
fig, axes = plt.subplots(1, 3, figsize=(15, 5))

# Histogram with mean and median lines
axes[0].hist(stunting_values, bins=10, edgecolor='black', alpha=0.7, color='steelblue')
axes[0].axvline(mean_val, color='red', linestyle='--', linewidth=2, label=f'Mean: {mean_val:.1f}%')
axes[0].axvline(median_val, color='green', linestyle='-', linewidth=2, label=f'Median: {median_val:.1f}%')
axes[0].set_xlabel('Stunting Prevalence (%)')
axes[0].set_ylabel('Frequency')
axes[0].set_title('Distribution of Stunting Prevalence')
axes[0].legend()

# Box plot showing IQR and outliers
axes[1].boxplot(stunting_values, vert=True)
axes[1].set_ylabel('Stunting Prevalence (%)')
axes[1].set_title('Box Plot: Stunting Prevalence')

# Empirical CDF with WHO threshold
sorted_vals = np.sort(stunting_values)
cdf = np.arange(1, len(sorted_vals) + 1) / len(sorted_vals)
axes[2].plot(sorted_vals, cdf, marker='o', markersize=4, color='steelblue')
axes[2].axvline(30, color='red', linestyle='--', label='WHO "Very High" threshold (30%)')
axes[2].set_xlabel('Stunting Prevalence (%)')
axes[2].set_ylabel('Cumulative Probability')
axes[2].set_title('Empirical CDF')
axes[2].legend()

plt.tight_layout()
plt.savefig('distribution_shape.png', dpi=150, bbox_inches='tight')
plt.show()

What does this code do?

This creates a single figure with three side-by-side panels using plt.subplots(1, 3).

  • Panel 1 (histogram): bars show how many stunting values fall into each bin. The red dashed line marks the mean, and the green solid line marks the median. If the red line sits to the right of the green line, the distribution has positive skew.
  • Panel 2 (box plot): the box spans the IQR (25th to 75th percentile), the line inside the box is the median, and dots outside the whiskers are outliers.
  • Panel 3 (empirical CDF): each point shows the fraction of values at or below a given stunting percentage. The red vertical line at 30% marks the WHO threshold for "very high" public health significance.
  • Run the cell by pressing **Shift+Enter**.
  • You should see a wide figure with three charts. In the histogram, notice whether the red mean line sits to the right of the green median line. This confirms the positive skew you calculated in the previous step.
  • The file distribution_shape.png is now saved in your project folder.

Not seeing the figure?

  • Make sure you ran all previous cells first. This cell relies on stunting_values, mean_val, and median_val from earlier cells.
  • If you see a blank output, check that plt.show() is at the very end of the cell.
  • If the histogram looks empty, verify that stunting_values has data by running len(stunting_values) in a new cell. You should see a number greater than 0.

Why isn't my matplotlib figure appearing in Jupyter Notebook?

Compare malnutrition indicators across regions

Nigeria has six geopolitical zones. A grouped bar chart places all three malnutrition indicators side by side for each zone, making the North-South disparity immediately visible.

  • Create a new cell (Cell 11) and add this regional comparison chart:
# Cell 11: Regional Bar Chart
region_data = nutrition_wide[nutrition_wide['CharacteristicCategory'] == 'Region'].copy()
if not region_data.empty:
    # Extract regional values for each indicator
    regions = region_data['CharacteristicLabel'].values
    stunting_by_region = region_data['Stunting'].values
    wasting_by_region = region_data['Wasting'].values
    underweight_by_region = region_data['Underweight'].values

    # Set up grouped bar positions
    x = np.arange(len(regions))
    width = 0.25

    fig, ax = plt.subplots(figsize=(12, 6))
    bars1 = ax.bar(x - width, stunting_by_region, width, label='Stunting', color='#e74c3c')
    bars2 = ax.bar(x, wasting_by_region, width, label='Wasting', color='#f39c12')
    bars3 = ax.bar(x + width, underweight_by_region, width, label='Underweight', color='#3498db')

    ax.set_xlabel('Geopolitical Zone')
    ax.set_ylabel('Prevalence (%)')
    ax.set_title('Child Malnutrition by Region in Nigeria')
    ax.set_xticks(x)
    ax.set_xticklabels(regions, rotation=45, ha='right')
    ax.legend()
    # WHO reference line for "very high" stunting
    ax.axhline(30, color='gray', linestyle=':', alpha=0.5, label='WHO Very High (30%)')

    plt.tight_layout()
    plt.savefig('regional_comparison.png', dpi=150, bbox_inches='tight')
    plt.show()

What does this code do?

  • First, it filters nutrition_wide to rows where CharacteristicCategory == 'Region', giving you one row per geopolitical zone.
  • The grouped bar chart places three bars (Stunting in red, Wasting in orange, Underweight in blue) next to each other for every zone. The width = 0.25 controls bar thickness so they fit side by side.
  • The horizontal gray dotted line at 30% marks the WHO threshold for "very high" stunting prevalence. Any zone whose red bar reaches above this line has a critical stunting problem.
  • Run the cell by pressing **Shift+Enter**.
  • You should see a bar chart with six groups along the x-axis (one per geopolitical zone). The Northern zones (North East, North West) should have visibly taller red bars than the Southern zones.
  • The file regional_comparison.png is now saved in your project folder.

Seeing an empty chart or no bars?

  • If the chart is blank, the if not region_data.empty: check may be failing. Run region_data in a new cell to see if it has rows. If it is empty, check that your CharacteristicCategory column contains the exact string 'Region' (case-sensitive).
  • If x-axis labels overlap or appear cut off, verify that plt.tight_layout() is included before plt.show().

My regional bar chart is empty or has missing bars. Help me debug.

Identify regions above the WHO threshold

The WHO classifies stunting prevalence above 30% as "very high" public health significance. Both of your visualizations now mark this threshold.

  • Look at your empirical CDF (the rightmost panel in your distribution figure). The red vertical line sits at 30%. Any data point to the right of that line represents a characteristic breakdown with "very high" stunting.
  • Now look at your regional bar chart. The horizontal gray dotted line marks the same 30% threshold. Count how many zones have their red (stunting) bar above this line.
  • Compare the Northern zones (North Central, North East, North West) against the Southern zones (South East, South South, South West). The gap between them is the North-South disparity that drives the positive skew you calculated earlier.

Why does this disparity cause skewness?

When you pool all stunting values together, most fall in a moderate range. But the Northern zones pull the tail of the distribution far to the right.

This is exactly why the mean (pulled toward the high Northern values) sits above the median (the central value unaffected by extremes). For a policymaker, reporting the median is more honest because it represents the typical region rather than being inflated by outliers.

✔️ Awesome, I've got everything!

Great. Make sure both distribution_shape.png and regional_comparison.png are saved in your project folder, then move on.

ⓧ I'd like to double check the full code

Here is the complete code for Cells 10 and 11. Compare it to your notebook:

# Cell 10: Visualization - Distribution Shape
fig, axes = plt.subplots(1, 3, figsize=(15, 5))

axes[0].hist(stunting_values, bins=10, edgecolor='black', alpha=0.7, color='steelblue')
axes[0].axvline(mean_val, color='red', linestyle='--', linewidth=2, label=f'Mean: {mean_val:.1f}%')
axes[0].axvline(median_val, color='green', linestyle='-', linewidth=2, label=f'Median: {median_val:.1f}%')
axes[0].set_xlabel('Stunting Prevalence (%)')
axes[0].set_ylabel('Frequency')
axes[0].set_title('Distribution of Stunting Prevalence')
axes[0].legend()

axes[1].boxplot(stunting_values, vert=True)
axes[1].set_ylabel('Stunting Prevalence (%)')
axes[1].set_title('Box Plot: Stunting Prevalence')

sorted_vals = np.sort(stunting_values)
cdf = np.arange(1, len(sorted_vals) + 1) / len(sorted_vals)
axes[2].plot(sorted_vals, cdf, marker='o', markersize=4, color='steelblue')
axes[2].axvline(30, color='red', linestyle='--', label='WHO "Very High" threshold (30%)')
axes[2].set_xlabel('Stunting Prevalence (%)')
axes[2].set_ylabel('Cumulative Probability')
axes[2].set_title('Empirical CDF')
axes[2].legend()

plt.tight_layout()
plt.savefig('distribution_shape.png', dpi=150, bbox_inches='tight')
plt.show()

# Cell 11: Regional Bar Chart
region_data = nutrition_wide[nutrition_wide['CharacteristicCategory'] == 'Region'].copy()
if not region_data.empty:
    regions = region_data['CharacteristicLabel'].values
    stunting_by_region = region_data['Stunting'].values
    wasting_by_region = region_data['Wasting'].values
    underweight_by_region = region_data['Underweight'].values

    x = np.arange(len(regions))
    width = 0.25

    fig, ax = plt.subplots(figsize=(12, 6))
    bars1 = ax.bar(x - width, stunting_by_region, width, label='Stunting', color='#e74c3c')
    bars2 = ax.bar(x, wasting_by_region, width, label='Wasting', color='#f39c12')
    bars3 = ax.bar(x + width, underweight_by_region, width, label='Underweight', color='#3498db')

    ax.set_xlabel('Geopolitical Zone')
    ax.set_ylabel('Prevalence (%)')
    ax.set_title('Child Malnutrition by Region in Nigeria')
    ax.set_xticks(x)
    ax.set_xticklabels(regions, rotation=45, ha='right')
    ax.legend()
    ax.axhline(30, color='gray', linestyle=':', alpha=0.5, label='WHO Very High (30%)')

    plt.tight_layout()
    plt.savefig('regional_comparison.png', dpi=150, bbox_inches='tight')
    plt.show()

Your visualizations now tell the story that your numbers hinted at. The North-South disparity is no longer an abstract skewness value. It is a visible gap in the bar chart. Next up, you will use probability theory to quantify exactly how much knowing a region's stunting status updates your belief about whether that region is in Northern Nigeria.

Apply Probability and Bayes' Theorem

Your visualizations made the North-South malnutrition gap impossible to ignore. Now you can quantify that relationship using probability.

In this step, you will define events from your data, compute marginal and conditional probabilities, test for statistical independence, and apply Bayes' theorem to answer a real question: given that a region has high stunting, what is the probability it is in Northern Nigeria?

In this step, get ready to:
  • Construct a probability table defining events from your regional stunting data.
  • Test whether region and stunting are statistically independent.
  • Apply Bayes' theorem by hand to compute a posterior probability.
Construct a probability table from your regional data

Probability gives you a formal way to express how likely events are and how they relate to each other. You will define two events directly from your cleaned data and compute their probabilities.

  • In your notebook, create a new cell below Cell 11 (the regional bar chart). Add the following code to define your events and compute probabilities:
# Cell 12: Probability and Bayes' Theorem
print("=" * 50)
print("PROBABILITY AND BAYES' THEOREM")
print("=" * 50)

# Define events using regional data
# A = Region is in Northern Nigeria (North Central, North East, North West)
# B = Stunting prevalence is above the national median

northern_zones = ['North Central', 'North East', 'North West']
southern_zones = ['South East', 'South South', 'South West']

all_regions = region_data.copy()
national_median_stunting = compute_median(all_regions['Stunting'].values)
print(f"\nNational median stunting: {national_median_stunting:.1f}%")

# Classify regions
all_regions['IsNorth'] = all_regions['CharacteristicLabel'].isin(northern_zones)
all_regions['HighStunting'] = all_regions['Stunting'] > national_median_stunting

total_regions = len(all_regions)
n_north = all_regions['IsNorth'].sum()
n_high_stunting = all_regions['HighStunting'].sum()
n_north_and_high = ((all_regions['IsNorth']) & (all_regions['HighStunting'])).sum()

# Marginal probabilities
P_A = n_north / total_regions
P_B = n_high_stunting / total_regions
P_A_and_B = n_north_and_high / total_regions

# Conditional probability
P_B_given_A = n_north_and_high / n_north if n_north > 0 else 0

print(f"\nP(North) = P(A) = {n_north}/{total_regions} = {P_A:.3f}")
print(f"P(High Stunting) = P(B) = {n_high_stunting}/{total_regions} = {P_B:.3f}")
print(f"P(North AND High Stunting) = P(A and B) = {n_north_and_high}/{total_regions} = {P_A_and_B:.3f}")
print(f"P(High Stunting | North) = P(B|A) = {n_north_and_high}/{n_north} = {P_B_given_A:.3f}")

What does this code do?

You are defining two events from your real data. Event A is "region is in Northern Nigeria" (the three northern geopolitical zones combined). Event B is "stunting prevalence is above the national median."

  • P(A) is the marginal probability that a randomly chosen region is in the North.
  • P(B) is the marginal probability that a region has above-median stunting.
  • P(A and B) is the joint probability that a region is both Northern AND has high stunting.
  • P(B|A) is the conditional probability of high stunting given the region is Northern.
  • Run the cell by pressing Shift+Enter.
  • You should see the national median stunting value printed, followed by four probability values (P(A), P(B), P(A and B), P(B|A)) each expressed as fractions and decimals.

Verify probabilities are valid

All probability values should be between 0 and 1. If any value is negative or greater than 1, something went wrong with your data filtering in earlier cells.

Seeing a KeyError or empty output?

  • Make sure Cell 11 ran successfully. The variable region_data must exist before this cell can use it.
  • Check that your CharacteristicCategory column contains the value Region (case-sensitive). If the filter returned zero rows, region_data will be empty.

Help me debug why my region_data DataFrame is empty or causing a KeyError.

Test whether region and stunting are independent

Two events are independent if knowing one gives you no information about the other. The formal test is whether P(A and B) equals P(A) * P(B). If they differ, the events are statistically associated.

  • Create a new cell and add the independence test by running this code:
# Cell 13: Independence Test
print("\n" + "-" * 50)
print("INDEPENDENCE TEST")
print("-" * 50)
P_A_times_P_B = P_A * P_B
print(f"P(A) * P(B) = {P_A:.3f} * {P_B:.3f} = {P_A_times_P_B:.3f}")
print(f"P(A and B) = {P_A_and_B:.3f}")
print(f"\nIf independent, P(A and B) should equal P(A)*P(B)")
if abs(P_A_and_B - P_A_times_P_B) > 0.01:
    print("RESULT: Events are NOT independent. Region and stunting are associated.")
else:
    print("RESULT: Events appear approximately independent.")

What does this code do?

If events A and B were independent (knowing one tells you nothing about the other), then P(A and B) would exactly equal P(A) * P(B). This code compares the two values.

A difference greater than 0.01 means the events are associated. In plain language: being in Northern Nigeria is statistically linked to having high stunting. This is not a coincidence in the data.

  • Run the cell. You should see the computed P(A) * P(B) value, the actual P(A and B) value, and a result statement confirming the events are NOT independent.

Why does independence matter here?

If region and stunting were independent, targeting interventions by geography would be pointless. The fact they are associated confirms that geographic targeting (directing nutrition programs to the North) is statistically justified.

Seeing 'approximately independent' when you expect association?

  • Double check that national_median_stunting was computed from your regional data in Cell 12. If it came from a different subset, the threshold may not split regions meaningfully.
  • Confirm that your northern_zones list matches the exact labels in your data (case and spacing matter).

Help me debug why my independence test is showing the wrong result.

Apply Bayes' theorem to compute the posterior probability

Bayes' theorem lets you flip a conditional probability. You already know P(B|A), the probability of high stunting given a Northern region. Now you will compute P(A|B), the probability a region is Northern given it has high stunting.

  • Create a new cell and apply Bayes' theorem by running this code:
# Cell 14: Bayes' Theorem
print("\n" + "-" * 50)
print("BAYES' THEOREM: P(North | High Stunting)")
print("-" * 50)
# P(A|B) = P(B|A) * P(A) / P(B)
P_A_given_B = (P_B_given_A * P_A) / P_B if P_B > 0 else 0

print(f"\nP(North | High Stunting) = P(B|A) * P(A) / P(B)")
print(f"                        = {P_B_given_A:.3f} * {P_A:.3f} / {P_B:.3f}")
print(f"                        = {P_A_given_B:.3f}")
print(f"\nInterpretation:")
print(f"  Prior probability of North: P(A) = {P_A:.3f}")
print(f"  Posterior probability:      P(A|B) = {P_A_given_B:.3f}")
print(f"  Update: Knowing a region has high stunting {'increases' if P_A_given_B > P_A else 'decreases'} ")
print(f"  our belief it is Northern by {abs(P_A_given_B - P_A)*100:.1f} percentage points.")
print(f"\n  This means: If you learn a region has above-median stunting,")
print(f"  the probability it is in Northern Nigeria jumps from {P_A*100:.0f}% to {P_A_given_B*100:.0f}%.")

What does this code do?

Bayes' theorem computes P(A|B) = P(B|A) * P(A) / P(B). In plain language: the probability a region is Northern, given you know it has high stunting.

  • The prior P(A) is your belief before seeing any stunting data. It is simply the fraction of regions that are Northern.
  • The posterior P(A|B) is your updated belief after learning that a region has above-median stunting.
  • The difference between prior and posterior tells you how much the stunting evidence shifts your belief about geographic location.
  • Run the cell. You should see the full Bayes' theorem calculation printed step by step, followed by an interpretation comparing the prior to the posterior probability.
  • Save your notebook by pressing Cmd+S (macOS) or Ctrl+S (Windows).

Why does the posterior jump matter for policy?

If the posterior is much higher than the prior, it means stunting data is highly informative about geography. A public health official who learns a region has high stunting can be significantly more confident it is in the North. This justifies directing nutrition interventions toward Northern Nigeria specifically.

Getting a ZeroDivisionError?

  • This happens if P_B is zero, meaning no regions have stunting above the median. Make sure you ran Cell 12 first so that P_B is defined.
  • If P_B prints as 0.000, check that your HighStunting column uses a strict greater-than comparison (not greater-than-or-equal). The median region itself should not count as "above median."

Help me fix my Bayes' theorem calculation that is failing with a division error.

✔️ Awesome, I've got everything!

Great. Double check you have saved your notebook.

ⓧ I'd like to double check the full code

Here is the complete notebook code up to this point. Compare your cells against this reference:

# Cell 1: Imports and Setup
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt

# Cell 2: Load Data
anthropometry = pd.read_csv('anthropometry_nigeria.csv')
nutrition = pd.read_csv('nutrition_indicators_nigeria.csv')

print("Anthropometry shape:", anthropometry.shape)
print("Nutrition shape:", nutrition.shape)
print("\nAnthropometry columns:", anthropometry.columns.tolist())
print("\nFirst few rows:")
anthropometry.head()

# Cell 3: Inspect Data Structure
print(anthropometry.info())
print("\nUnique indicators:", anthropometry['Indicator'].nunique())
print("\nSample indicators:")
print(anthropometry['Indicator'].unique()[:10])
print("\nCharacteristic categories:")
print(anthropometry['CharacteristicCategory'].unique())
print("\nSurvey years:", anthropometry['SurveyYear'].unique())

# Cell 4: Filter and Clean
latest_year = anthropometry['SurveyYear'].max()
anthro_recent = anthropometry[anthropometry['SurveyYear'] == latest_year].copy()
nutr_recent = nutrition[nutrition['SurveyYear'] == latest_year].copy()

stunting_ids = anthro_recent['IndicatorId'].str.contains('CN_NUTS_C_HA2', na=False)
wasting_ids = anthro_recent['IndicatorId'].str.contains('CN_NUTS_C_WH2', na=False)
underweight_ids = anthro_recent['IndicatorId'].str.contains('CN_NUTS_C_WA2', na=False)

nutrition_data = anthro_recent[stunting_ids | wasting_ids | underweight_ids].copy()
print(f"Filtered to {len(nutrition_data)} rows for stunting/wasting/underweight")

# Cell 5: Handle Missing Values and Type Conversion
nutrition_data['Value'] = pd.to_numeric(nutrition_data['Value'], errors='coerce')
print(f"Missing values in Value column: {nutrition_data['Value'].isna().sum()}")
nutrition_data = nutrition_data.dropna(subset=['Value'])

nutrition_data['CharacteristicLabel'] = nutrition_data['CharacteristicLabel'].str.strip()
nutrition_data['CharacteristicCategory'] = nutrition_data['CharacteristicCategory'].str.strip()

duplicates = nutrition_data.duplicated(subset=['IndicatorId', 'CharacteristicCategory', 'CharacteristicLabel'])
print(f"Duplicate rows: {duplicates.sum()}")
nutrition_data = nutrition_data.drop_duplicates(subset=['IndicatorId', 'CharacteristicCategory', 'CharacteristicLabel'])

print(f"\nCleaned dataset shape: {nutrition_data.shape}")
print(f"Characteristic categories: {nutrition_data['CharacteristicCategory'].unique()}")

# Cell 6: Reshape Data - Create Indicator Column
def label_indicator(indicator_id):
    if 'HA2' in str(indicator_id):
        return 'Stunting'
    elif 'WH2' in str(indicator_id):
        return 'Wasting'
    elif 'WA2' in str(indicator_id):
        return 'Underweight'
    return 'Other'

nutrition_data['NutritionIndicator'] = nutrition_data['IndicatorId'].apply(label_indicator)

nutrition_wide = nutrition_data.pivot_table(
    index=['CharacteristicCategory', 'CharacteristicLabel'],
    columns='NutritionIndicator',
    values='Value',
    aggfunc='first'
).reset_index()

print("Reshaped data:")
print(nutrition_wide.head(10))
print(f"\nFinal shape: {nutrition_wide.shape}")
print(nutrition_wide.describe())

# Cell 7: Descriptive Statistics from Scratch
stunting_values = nutrition_wide['Stunting'].dropna().values

def compute_mean(data):
    return np.sum(data) / len(data)

def compute_median(data):
    sorted_data = np.sort(data)
    n = len(sorted_data)
    if n % 2 == 0:
        return (sorted_data[n//2 - 1] + sorted_data[n//2]) / 2
    else:
        return sorted_data[n//2]

def compute_mode(data):
    values, counts = np.unique(np.round(data, 1), return_counts=True)
    max_count = np.max(counts)
    modes = values[counts == max_count]
    return modes

def compute_variance(data, sample=True):
    mean = compute_mean(data)
    squared_diffs = (data - mean) ** 2
    if sample:
        return np.sum(squared_diffs) / (len(data) - 1)
    else:
        return np.sum(squared_diffs) / len(data)

def compute_std(data, sample=True):
    return np.sqrt(compute_variance(data, sample))

def compute_range(data):
    return np.max(data) - np.min(data)

def compute_iqr(data):
    sorted_data = np.sort(data)
    q1 = np.percentile(sorted_data, 25)
    q3 = np.percentile(sorted_data, 75)
    return q3 - q1

print("=" * 50)
print("DESCRIPTIVE STATISTICS: Stunting Prevalence (%)")
print("=" * 50)
print(f"Mean:                {compute_mean(stunting_values):.2f}")
print(f"Median:              {compute_median(stunting_values):.2f}")
print(f"Mode:                {compute_mode(stunting_values)}")
print(f"Sample Variance:     {compute_variance(stunting_values):.2f}")
print(f"Sample Std Dev:      {compute_std(stunting_values):.2f}")
print(f"Range:               {compute_range(stunting_values):.2f}")
print(f"IQR:                 {compute_iqr(stunting_values):.2f}")

# Cell 8: Skewness and Kurtosis
def compute_skewness(data):
    n = len(data)
    mean = compute_mean(data)
    std = compute_std(data, sample=True)
    skew = (n / ((n-1) * (n-2))) * np.sum(((data - mean) / std) ** 3)
    return skew

def compute_kurtosis(data):
    n = len(data)
    mean = compute_mean(data)
    std = compute_std(data, sample=True)
    kurt = ((n * (n+1)) / ((n-1) * (n-2) * (n-3))) * np.sum(((data - mean) / std) ** 4)
    kurt -= (3 * (n-1)**2) / ((n-2) * (n-3))
    return kurt

print(f"\nSkewness:            {compute_skewness(stunting_values):.3f}")
print(f"Kurtosis (excess):   {compute_kurtosis(stunting_values):.3f}")

for p in [25, 50, 75, 90]:
    val = np.percentile(stunting_values, p)
    print(f"{p}th percentile:     {val:.2f}%")

# Cell 9: Mean vs Median Comparison
print("\n" + "=" * 50)
print("WHY MEDIAN BEATS MEAN FOR SKEWED DATA")
print("=" * 50)
mean_val = compute_mean(stunting_values)
median_val = compute_median(stunting_values)
print(f"Mean stunting:   {mean_val:.2f}%")
print(f"Median stunting: {median_val:.2f}%")
print(f"Difference:      {mean_val - median_val:.2f} percentage points")
print(f"\nThe mean is pulled {'upward' if mean_val > median_val else 'downward'} by extreme values.")
print("For policy decisions, the median better represents the 'typical' region.")

# Cell 10: Visualization - Distribution Shape
fig, axes = plt.subplots(1, 3, figsize=(15, 5))

axes[0].hist(stunting_values, bins=10, edgecolor='black', alpha=0.7, color='steelblue')
axes[0].axvline(mean_val, color='red', linestyle='--', linewidth=2, label=f'Mean: {mean_val:.1f}%')
axes[0].axvline(median_val, color='green', linestyle='-', linewidth=2, label=f'Median: {median_val:.1f}%')
axes[0].set_xlabel('Stunting Prevalence (%)')
axes[0].set_ylabel('Frequency')
axes[0].set_title('Distribution of Stunting Prevalence')
axes[0].legend()

axes[1].boxplot(stunting_values, vert=True)
axes[1].set_ylabel('Stunting Prevalence (%)')
axes[1].set_title('Box Plot: Stunting Prevalence')

sorted_vals = np.sort(stunting_values)
cdf = np.arange(1, len(sorted_vals) + 1) / len(sorted_vals)
axes[2].plot(sorted_vals, cdf, marker='o', markersize=4, color='steelblue')
axes[2].axvline(30, color='red', linestyle='--', label='WHO "Very High" threshold (30%)')
axes[2].set_xlabel('Stunting Prevalence (%)')
axes[2].set_ylabel('Cumulative Probability')
axes[2].set_title('Empirical CDF')
axes[2].legend()

plt.tight_layout()
plt.savefig('distribution_shape.png', dpi=150, bbox_inches='tight')
plt.show()

# Cell 11: Regional Bar Chart
region_data = nutrition_wide[nutrition_wide['CharacteristicCategory'] == 'Region'].copy()
if not region_data.empty:
    regions = region_data['CharacteristicLabel'].values
    stunting_by_region = region_data['Stunting'].values
    wasting_by_region = region_data['Wasting'].values
    underweight_by_region = region_data['Underweight'].values

    x = np.arange(len(regions))
    width = 0.25

    fig, ax = plt.subplots(figsize=(12, 6))
    bars1 = ax.bar(x - width, stunting_by_region, width, label='Stunting', color='#e74c3c')
    bars2 = ax.bar(x, wasting_by_region, width, label='Wasting', color='#f39c12')
    bars3 = ax.bar(x + width, underweight_by_region, width, label='Underweight', color='#3498db')

    ax.set_xlabel('Geopolitical Zone')
    ax.set_ylabel('Prevalence (%)')
    ax.set_title('Child Malnutrition by Region in Nigeria')
    ax.set_xticks(x)
    ax.set_xticklabels(regions, rotation=45, ha='right')
    ax.legend()
    ax.axhline(30, color='gray', linestyle=':', alpha=0.5, label='WHO Very High (30%)')

    plt.tight_layout()
    plt.savefig('regional_comparison.png', dpi=150, bbox_inches='tight')
    plt.show()

# Cell 12: Probability and Bayes' Theorem
print("=" * 50)
print("PROBABILITY AND BAYES' THEOREM")
print("=" * 50)

northern_zones = ['North Central', 'North East', 'North West']
southern_zones = ['South East', 'South South', 'South West']

all_regions = region_data.copy()
national_median_stunting = compute_median(all_regions['Stunting'].values)
print(f"\nNational median stunting: {national_median_stunting:.1f}%")

all_regions['IsNorth'] = all_regions['CharacteristicLabel'].isin(northern_zones)
all_regions['HighStunting'] = all_regions['Stunting'] > national_median_stunting

total_regions = len(all_regions)
n_north = all_regions['IsNorth'].sum()
n_high_stunting = all_regions['HighStunting'].sum()
n_north_and_high = ((all_regions['IsNorth']) & (all_regions['HighStunting'])).sum()

P_A = n_north / total_regions
P_B = n_high_stunting / total_regions
P_A_and_B = n_north_and_high / total_regions
P_B_given_A = n_north_and_high / n_north if n_north > 0 else 0

print(f"\nP(North) = P(A) = {n_north}/{total_regions} = {P_A:.3f}")
print(f"P(High Stunting) = P(B) = {n_high_stunting}/{total_regions} = {P_B:.3f}")
print(f"P(North AND High Stunting) = P(A and B) = {n_north_and_high}/{total_regions} = {P_A_and_B:.3f}")
print(f"P(High Stunting | North) = P(B|A) = {n_north_and_high}/{n_north} = {P_B_given_A:.3f}")

# Cell 13: Independence Test
print("\n" + "-" * 50)
print("INDEPENDENCE TEST")
print("-" * 50)
P_A_times_P_B = P_A * P_B
print(f"P(A) * P(B) = {P_A:.3f} * {P_B:.3f} = {P_A_times_P_B:.3f}")
print(f"P(A and B) = {P_A_and_B:.3f}")
print(f"\nIf independent, P(A and B) should equal P(A)*P(B)")
if abs(P_A_and_B - P_A_times_P_B) > 0.01:
    print("RESULT: Events are NOT independent. Region and stunting are associated.")
else:
    print("RESULT: Events appear approximately independent.")

# Cell 14: Bayes' Theorem
print("\n" + "-" * 50)
print("BAYES' THEOREM: P(North | High Stunting)")
print("-" * 50)
P_A_given_B = (P_B_given_A * P_A) / P_B if P_B > 0 else 0

print(f"\nP(North | High Stunting) = P(B|A) * P(A) / P(B)")
print(f"                        = {P_B_given_A:.3f} * {P_A:.3f} / {P_B:.3f}")
print(f"                        = {P_A_given_B:.3f}")
print(f"\nInterpretation:")
print(f"  Prior probability of North: P(A) = {P_A:.3f}")
print(f"  Posterior probability:      P(A|B) = {P_A_given_B:.3f}")
print(f"  Update: Knowing a region has high stunting {'increases' if P_A_given_B > P_A else 'decreases'} ")
print(f"  our belief it is Northern by {abs(P_A_given_B - P_A)*100:.1f} percentage points.")
print(f"\n  This means: If you learn a region has above-median stunting,")
print(f"  the probability it is in Northern Nigeria jumps from {P_A*100:.0f}% to {P_A_given_B*100:.0f}%.")

You have now used real nutrition data to compute probabilities from scratch, test for statistical independence, and apply Bayes' theorem to answer a genuine public health question. Next up, you will extend this analysis in the Secret Mission by examining how wealth quintiles relate to stunting.

Secret mission

Wealth Quintile Analysis and Independence Test

In the main project, you tested whether region and stunting are independent using just two categories (North vs South). But what about wealth? Nigeria's DHS data breaks households into five wealth quintiles. Can you extend your Bayesian analysis across all five categories, formally test for independence, and compute posterior probabilities that reveal how strongly poverty predicts malnutrition?

Clean Up Your Resources

Clean Up Your Resources

This project runs entirely locally, so there are no ongoing costs. Decide whether to keep your files for future reference or delete them to free up disk space.

Resources you used:

  • Jupyter notebook file (nigeria_nutrition_analysis.ipynb).
  • Downloaded CSV files (anthropometry_nigeria.csv and nutrition_indicators_nigeria.csv).
  • Generated visualization files (distribution_shape.png, regional_comparison.png, wealth_stunting.png).

Keep everything running

No action needed. Choose this if you want to revisit or extend your analysis later.

  • Your notebook and data files are stored locally on your machine.
  • There are no ongoing costs or cloud resources to worry about.
  • You can reopen nigeria_nutrition_analysis.ipynb any time to review your work or add new analysis cells.

Pause - I'll come back to this later

Shut down your Jupyter Notebook server to free up system memory. Your files remain untouched on disk.

  • In the terminal where Jupyter is running, press Ctrl+C to stop the notebook server.
  • Confirm the shutdown when prompted.
  • Your notebook, CSV files, and PNG images will still be in your project folder when you return.

Delete - I don't want to use this again

Remove all project files from your machine. This frees disk space and removes the downloaded data.

  • Stop the Jupyter Notebook server by pressing Ctrl+C in the terminal.
  • Navigate to the folder containing your project files and delete all of them by running:
rm nigeria_nutrition_analysis.ipynb anthropometry_nigeria.csv nutrition_indicators_nigeria.csv distribution_shape.png regional_comparison.png wealth_stunting.png

On Windows, use this command instead:

del nigeria_nutrition_analysis.ipynb, anthropometry_nigeria.csv, nutrition_indicators_nigeria.csv, distribution_shape.png, regional_comparison.png, wealth_stunting.png
  • Confirm the files are gone by running ls (macOS/Linux) or dir (Windows). You should see none of the project files listed.

Nice Work!

Nice Work!

You just built a complete data analysis notebook that takes raw Nigerian health survey data all the way from messy CSVs to Bayesian inference. That is the exact workflow expected in data science and public health analytics roles.

You've learned how to:

  • Clean and reshape real-world DHS nutrition data from the Humanitarian Data Exchange, handling missing values, standardizing categories, filtering survey years, merging multiple files, and pivoting from long to wide format.
  • Compute all core descriptive statistics from scratch using only NumPy: mean, median, mode, variance, standard deviation, range, IQR, skewness, kurtosis, and percentiles.
  • Apply Bayes' theorem by hand to answer a genuine public health question, computing joint, marginal, and conditional probabilities to determine how knowing a region has high stunting updates your belief about whether it is in Northern Nigeria.
  • Secret Mission: Analyzed the wealth-stunting relationship across 5 quintiles and formally tested for statistical independence.

Ready to quiz yourself?