---
title: "Global acute malnutrition by commune"
subtitle: "MUAC screening — Artibonite, 2024"
author: "Cassion · data-analysis.cassion.dev"
format:
  html:
    toc: true
    code-fold: false
jupyter: python3
---

## What this produces

GAM and SAM prevalence by commune, with 95% confidence intervals, from a
community mass screening register. The reference figures are in the dataset's
quality notes: overall GAM near 8.6% and SAM near 2.2%, ranging from roughly 5%
to 15% by commune. If your numbers land far outside that, you have a bug rather
than a finding.

Every dataset on this platform is synthetic. Nothing here describes a real
child, and these figures must never be cited as real prevalence.

## Setup

The download makes this runnable in Colab, where there is no local file. Running
it locally against your own copy is the same code with a different path.

```{python}
import pandas as pd
import numpy as np

URL = (
    "https://data-analysis.cassion.dev/datasets/files/"
    "muac-screening-artibonite-2024.v1.csv"
)

muac = pd.read_csv(
    URL,
    dtype={"child_id": "string", "commune": "string", "sex": "string"},
    na_values={"muac_mm": ["-99"]},
)

muac["screening_date"] = pd.to_datetime(muac["screening_date"], format="%Y-%m-%d")

print(muac.shape)
muac.head()
```

Two arguments there are doing the work. `na_values` scoped to `muac_mm` stops
the missing-value sentinel `-99` entering the mean — without it the average MUAC
comes out several millimetres low, and nothing about the result looks wrong.
`dtype` keeps identifiers as text, which matters the moment you join to another
file.

## Clean what would corrupt the indicator

Two defects in this register change the answer. Seven records were left in
centimetres and never converted; the plausible millimetre and centimetre ranges
do not overlap, so the correction is unambiguous. Two communes recorded oedema
as `Y`/`N` in the first quarter rather than `true`/`false`.

```{python}
unit_error = muac["muac_mm"].notna() & (muac["muac_mm"] < 40)
muac.loc[unit_error, "muac_mm"] = muac.loc[unit_error, "muac_mm"] * 10

implausible = muac["muac_mm"].notna() & (
    (muac["muac_mm"] < 80) | (muac["muac_mm"] > 220)
)
muac.loc[implausible, "muac_mm"] = np.nan

oedema_map = {
    "true": True, "TRUE": True, "Y": True, "y": True, "yes": True,
    "false": False, "FALSE": False, "N": False, "n": False, "no": False,
}
muac["oedema"] = muac["oedema"].astype("string").str.strip().map(oedema_map)

print(f"unit errors corrected : {int(unit_error.sum())}")
print(f"implausible to missing: {int(implausible.sum())}")
print(f"oedema still missing  : {int(muac['oedema'].isna().sum())}")
```

## Define the indicator before computing it

- **Numerator** — children with MUAC below 125 mm, or with bilateral pitting
  oedema.
- **Denominator** — children with a valid MUAC measurement **or** a recorded
  oedema assessment.
- **Disaggregation** — commune.

Oedema is severe acute malnutrition regardless of the measurement, so filtering
on `muac_mm` alone understates the caseload precisely among the most severe
cases.

```{python}
SAM_MM, GAM_MM = 115, 125

muac["has_assessment"] = muac["muac_mm"].notna() | muac["oedema"].notna()

muac["sam"] = np.where(
    ~muac["has_assessment"],
    np.nan,
    ((muac["muac_mm"] < SAM_MM) | (muac["oedema"] == True)).astype(float),
)
muac["gam"] = np.where(
    ~muac["has_assessment"],
    np.nan,
    ((muac["muac_mm"] < GAM_MM) | (muac["oedema"] == True)).astype(float),
)
```

## The indicator table

`screened` and `denominator` are separate columns on purpose: they differ by the
rows with no assessment at all, and a reader who sees both can judge how much of
the register the rate rests on.

```{python}
table = (
    muac.groupby("commune", dropna=False)
    .agg(
        screened=("child_id", "size"),
        denominator=("has_assessment", "sum"),
        sam_cases=("sam", "sum"),
        gam_cases=("gam", "sum"),
    )
    .reset_index()
)

table["gam_rate"] = table["gam_cases"] / table["denominator"]
table["sam_rate"] = table["sam_cases"] / table["denominator"]

table = table.sort_values("gam_rate", ascending=False)
table.round(4)
```

## Say how uncertain you are

A rate from 340 children is not the same claim as a rate from 40. The
Clopper-Pearson interval is exact for a proportion and does not misbehave at
small counts, which matters for the smaller communes.

```{python}
from scipy.stats import beta

def clopper_pearson(successes, n, confidence=0.95):
    if n == 0:
        return (np.nan, np.nan)
    alpha = 1 - confidence
    lower = beta.ppf(alpha / 2, successes, n - successes + 1) if successes > 0 else 0.0
    upper = beta.ppf(1 - alpha / 2, successes + 1, n - successes) if successes < n else 1.0
    return (lower, upper)

bounds = table.apply(
    lambda r: clopper_pearson(r["gam_cases"], r["denominator"]), axis=1
)
table["gam_low"] = [b[0] for b in bounds]
table["gam_high"] = [b[1] for b in bounds]

table[["commune", "denominator", "gam_rate", "gam_low", "gam_high"]].round(4)
```

## Plot the ranking with its uncertainty

```{python}
import matplotlib.pyplot as plt

EMERGENCY_THRESHOLD = 0.15

fig, ax = plt.subplots(figsize=(8, 5))
order = table.sort_values("gam_rate")

ax.errorbar(
    order["gam_rate"] * 100,
    range(len(order)),
    xerr=[
        (order["gam_rate"] - order["gam_low"]) * 100,
        (order["gam_high"] - order["gam_rate"]) * 100,
    ],
    fmt="o",
    color="#2F5D50",
    ecolor="#9AA8A3",
    capsize=3,
)

ax.axvline(EMERGENCY_THRESHOLD * 100, color="#B5533C", linestyle="--", linewidth=1)
ax.text(
    EMERGENCY_THRESHOLD * 100 + 0.3, 0.2,
    "emergency threshold (15%)", color="#B5533C", fontsize=9,
)

ax.set_yticks(range(len(order)))
ax.set_yticklabels(order["commune"])
ax.set_xlabel("GAM prevalence by MUAC (%)")
ax.set_title("Global acute malnutrition by commune, with 95% CI")
ax.spines[["top", "right"]].set_visible(False)

plt.tight_layout()
plt.show()
```

## Reading the result

Several intervals overlap. That means the ordering between those communes is not
supported by the data, and a decision that allocates one additional CMAM site by
rank alone is reading precision the screening does not have.

Note also what the interval does *not* cover. This is a census of the children
who came to be screened, not a probability sample, so the interval describes
sampling variation only. Whether those children resemble the ones who did not
come is usually the larger source of error, and it belongs in the limitations
section of any report built on this.