{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Global acute malnutrition by commune\n",
    "\n",
    "*MUAC screening — Artibonite, 2024*\n",
    "\n",
    "Cassion · data-analysis.cassion.dev\n",
    "\n",
    "[![Open In Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/thecassion/cassion-learning-platform/blob/main/apps/data-analysis/public/datasets/examples/muac-screening-artibonite-2024/gam-by-commune.python.en.ipynb)"
   ],
   "id": "cell-000"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## What this produces\n",
    "\n",
    "GAM and SAM prevalence by commune, with 95% confidence intervals, from a\n",
    "community mass screening register. The reference figures are in the dataset's\n",
    "quality notes: overall GAM near 8.6% and SAM near 2.2%, ranging from roughly 5%\n",
    "to 15% by commune. If your numbers land far outside that, you have a bug rather\n",
    "than a finding.\n",
    "\n",
    "Every dataset on this platform is synthetic. Nothing here describes a real\n",
    "child, and these figures must never be cited as real prevalence.\n",
    "\n",
    "## Setup\n",
    "\n",
    "The download makes this runnable in Colab, where there is no local file. Running\n",
    "it locally against your own copy is the same code with a different path."
   ],
   "id": "cell-001"
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "\n",
    "URL = (\n",
    "    \"https://data-analysis.cassion.dev/datasets/files/\"\n",
    "    \"muac-screening-artibonite-2024.v1.csv\"\n",
    ")\n",
    "\n",
    "muac = pd.read_csv(\n",
    "    URL,\n",
    "    dtype={\"child_id\": \"string\", \"commune\": \"string\", \"sex\": \"string\"},\n",
    "    na_values={\"muac_mm\": [\"-99\"]},\n",
    ")\n",
    "\n",
    "muac[\"screening_date\"] = pd.to_datetime(muac[\"screening_date\"], format=\"%Y-%m-%d\")\n",
    "\n",
    "print(muac.shape)\n",
    "muac.head()"
   ],
   "id": "cell-002"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Two arguments there are doing the work. `na_values` scoped to `muac_mm` stops\n",
    "the missing-value sentinel `-99` entering the mean — without it the average MUAC\n",
    "comes out several millimetres low, and nothing about the result looks wrong.\n",
    "`dtype` keeps identifiers as text, which matters the moment you join to another\n",
    "file.\n",
    "\n",
    "## Clean what would corrupt the indicator\n",
    "\n",
    "Two defects in this register change the answer. Seven records were left in\n",
    "centimetres and never converted; the plausible millimetre and centimetre ranges\n",
    "do not overlap, so the correction is unambiguous. Two communes recorded oedema\n",
    "as `Y`/`N` in the first quarter rather than `true`/`false`."
   ],
   "id": "cell-003"
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "unit_error = muac[\"muac_mm\"].notna() & (muac[\"muac_mm\"] < 40)\n",
    "muac.loc[unit_error, \"muac_mm\"] = muac.loc[unit_error, \"muac_mm\"] * 10\n",
    "\n",
    "implausible = muac[\"muac_mm\"].notna() & (\n",
    "    (muac[\"muac_mm\"] < 80) | (muac[\"muac_mm\"] > 220)\n",
    ")\n",
    "muac.loc[implausible, \"muac_mm\"] = np.nan\n",
    "\n",
    "oedema_map = {\n",
    "    \"true\": True, \"TRUE\": True, \"Y\": True, \"y\": True, \"yes\": True,\n",
    "    \"false\": False, \"FALSE\": False, \"N\": False, \"n\": False, \"no\": False,\n",
    "}\n",
    "muac[\"oedema\"] = muac[\"oedema\"].astype(\"string\").str.strip().map(oedema_map)\n",
    "\n",
    "print(f\"unit errors corrected : {int(unit_error.sum())}\")\n",
    "print(f\"implausible to missing: {int(implausible.sum())}\")\n",
    "print(f\"oedema still missing  : {int(muac['oedema'].isna().sum())}\")"
   ],
   "id": "cell-004"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Define the indicator before computing it\n",
    "\n",
    "- **Numerator** — children with MUAC below 125 mm, or with bilateral pitting\n",
    "  oedema.\n",
    "- **Denominator** — children with a valid MUAC measurement **or** a recorded\n",
    "  oedema assessment.\n",
    "- **Disaggregation** — commune.\n",
    "\n",
    "Oedema is severe acute malnutrition regardless of the measurement, so filtering\n",
    "on `muac_mm` alone understates the caseload precisely among the most severe\n",
    "cases."
   ],
   "id": "cell-005"
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "SAM_MM, GAM_MM = 115, 125\n",
    "\n",
    "muac[\"has_assessment\"] = muac[\"muac_mm\"].notna() | muac[\"oedema\"].notna()\n",
    "\n",
    "muac[\"sam\"] = np.where(\n",
    "    ~muac[\"has_assessment\"],\n",
    "    np.nan,\n",
    "    ((muac[\"muac_mm\"] < SAM_MM) | (muac[\"oedema\"] == True)).astype(float),\n",
    ")\n",
    "muac[\"gam\"] = np.where(\n",
    "    ~muac[\"has_assessment\"],\n",
    "    np.nan,\n",
    "    ((muac[\"muac_mm\"] < GAM_MM) | (muac[\"oedema\"] == True)).astype(float),\n",
    ")"
   ],
   "id": "cell-006"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The indicator table\n",
    "\n",
    "`screened` and `denominator` are separate columns on purpose: they differ by the\n",
    "rows with no assessment at all, and a reader who sees both can judge how much of\n",
    "the register the rate rests on."
   ],
   "id": "cell-007"
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "table = (\n",
    "    muac.groupby(\"commune\", dropna=False)\n",
    "    .agg(\n",
    "        screened=(\"child_id\", \"size\"),\n",
    "        denominator=(\"has_assessment\", \"sum\"),\n",
    "        sam_cases=(\"sam\", \"sum\"),\n",
    "        gam_cases=(\"gam\", \"sum\"),\n",
    "    )\n",
    "    .reset_index()\n",
    ")\n",
    "\n",
    "table[\"gam_rate\"] = table[\"gam_cases\"] / table[\"denominator\"]\n",
    "table[\"sam_rate\"] = table[\"sam_cases\"] / table[\"denominator\"]\n",
    "\n",
    "table = table.sort_values(\"gam_rate\", ascending=False)\n",
    "table.round(4)"
   ],
   "id": "cell-008"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Say how uncertain you are\n",
    "\n",
    "A rate from 340 children is not the same claim as a rate from 40. The\n",
    "Clopper-Pearson interval is exact for a proportion and does not misbehave at\n",
    "small counts, which matters for the smaller communes."
   ],
   "id": "cell-009"
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.stats import beta\n",
    "\n",
    "def clopper_pearson(successes, n, confidence=0.95):\n",
    "    if n == 0:\n",
    "        return (np.nan, np.nan)\n",
    "    alpha = 1 - confidence\n",
    "    lower = beta.ppf(alpha / 2, successes, n - successes + 1) if successes > 0 else 0.0\n",
    "    upper = beta.ppf(1 - alpha / 2, successes + 1, n - successes) if successes < n else 1.0\n",
    "    return (lower, upper)\n",
    "\n",
    "bounds = table.apply(\n",
    "    lambda r: clopper_pearson(r[\"gam_cases\"], r[\"denominator\"]), axis=1\n",
    ")\n",
    "table[\"gam_low\"] = [b[0] for b in bounds]\n",
    "table[\"gam_high\"] = [b[1] for b in bounds]\n",
    "\n",
    "table[[\"commune\", \"denominator\", \"gam_rate\", \"gam_low\", \"gam_high\"]].round(4)"
   ],
   "id": "cell-010"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plot the ranking with its uncertainty"
   ],
   "id": "cell-011"
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "EMERGENCY_THRESHOLD = 0.15\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(8, 5))\n",
    "order = table.sort_values(\"gam_rate\")\n",
    "\n",
    "ax.errorbar(\n",
    "    order[\"gam_rate\"] * 100,\n",
    "    range(len(order)),\n",
    "    xerr=[\n",
    "        (order[\"gam_rate\"] - order[\"gam_low\"]) * 100,\n",
    "        (order[\"gam_high\"] - order[\"gam_rate\"]) * 100,\n",
    "    ],\n",
    "    fmt=\"o\",\n",
    "    color=\"#2F5D50\",\n",
    "    ecolor=\"#9AA8A3\",\n",
    "    capsize=3,\n",
    ")\n",
    "\n",
    "ax.axvline(EMERGENCY_THRESHOLD * 100, color=\"#B5533C\", linestyle=\"--\", linewidth=1)\n",
    "ax.text(\n",
    "    EMERGENCY_THRESHOLD * 100 + 0.3, 0.2,\n",
    "    \"emergency threshold (15%)\", color=\"#B5533C\", fontsize=9,\n",
    ")\n",
    "\n",
    "ax.set_yticks(range(len(order)))\n",
    "ax.set_yticklabels(order[\"commune\"])\n",
    "ax.set_xlabel(\"GAM prevalence by MUAC (%)\")\n",
    "ax.set_title(\"Global acute malnutrition by commune, with 95% CI\")\n",
    "ax.spines[[\"top\", \"right\"]].set_visible(False)\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ],
   "id": "cell-012"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Reading the result\n",
    "\n",
    "Several intervals overlap. That means the ordering between those communes is not\n",
    "supported by the data, and a decision that allocates one additional CMAM site by\n",
    "rank alone is reading precision the screening does not have.\n",
    "\n",
    "Note also what the interval does *not* cover. This is a census of the children\n",
    "who came to be screened, not a probability sample, so the interval describes\n",
    "sampling variation only. Whether those children resemble the ones who did not\n",
    "come is usually the larger source of error, and it belongs in the limitations\n",
    "section of any report built on this."
   ],
   "id": "cell-013"
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "name": "python",
   "file_extension": ".py"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}