{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Computing prevalence with WHO growth standards\n",
    "\n",
    "*SMART nutrition survey, 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/smart-nutrition-survey-2024/whz-prevalence.python.en.ipynb)"
   ],
   "id": "cell-000"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## What this produces\n",
    "\n",
    "Weight-for-height z-scores against the WHO 2006 growth standards, computed from\n",
    "the LMS reference rather than read from a package — because the three decisions\n",
    "that move the answer are all outside the package call.\n",
    "\n",
    "Z-scores are deliberately not shipped in this dataset. Computing them is the\n",
    "exercise.\n",
    "\n",
    "Every dataset on this platform is synthetic. No real child is described, and\n",
    "these results must not be cited as a real nutrition situation.\n",
    "\n",
    "## The reference table\n",
    "\n",
    "R has the official WHO `anthro` package. Python's options are thinner, so this\n",
    "example reads the WHO 2006 weight-for-length and weight-for-height LMS tables\n",
    "directly. They ship with the platform, and the values are the same ones the WHO\n",
    "package uses — the R version of this example produces identical prevalence to the\n",
    "first decimal."
   ],
   "id": "cell-001"
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "\n",
    "BASE = \"https://data-analysis.cassion.dev/datasets/\"\n",
    "\n",
    "smart = pd.read_csv(BASE + \"files/smart-nutrition-survey-2024.v1.csv\",\n",
    "                    dtype={\"child_id\": \"string\", \"sex\": \"string\"})\n",
    "reference = pd.read_csv(BASE + \"reference/who-2006-weight-for-lenhei.csv\")\n",
    "\n",
    "print(f\"survey: {len(smart)} children\")\n",
    "print(f\"reference: {len(reference)} rows\")\n",
    "reference.head()"
   ],
   "id": "cell-002"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`lorh` is the standard the row belongs to: `L` for weight-for-length, used for\n",
    "children under two, and `H` for weight-for-height above that. `l`, `m` and `s`\n",
    "are the LMS parameters at that length or height.\n",
    "\n",
    "## Range-check before anything else\n",
    "\n",
    "Fourteen records hold impossible measurements — weights out by a factor of ten in\n",
    "both directions, and heights entered in metres. These are data entry errors, not\n",
    "statistical outliers, and they must go before any flagging rule, because a single\n",
    "114 kg child moves the survey mean that the SMART flag is measured against."
   ],
   "id": "cell-003"
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "impossible = (\n",
    "    (smart[\"weight_kg\"] < 2) | (smart[\"weight_kg\"] > 30)\n",
    "    | (smart[\"height_cm\"] < 45) | (smart[\"height_cm\"] > 130)\n",
    ")\n",
    "smart.loc[impossible, [\"child_id\", \"team\", \"weight_kg\", \"height_cm\"]]"
   ],
   "id": "cell-004"
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plausible = smart[~impossible.fillna(False)].copy()\n",
    "print(f\"kept {len(plausible)} of {len(smart)}\")\n",
    "print(f\"missing age: {int(plausible['age_months'].isna().sum())}, \"\n",
    "      f\"missing weight: {int(plausible['weight_kg'].isna().sum())}\")"
   ],
   "id": "cell-005"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Measurement position, which is a rule about age\n",
    "\n",
    "This is the step people skip. The WHO standard is not \"use whatever position was\n",
    "measured\" — it is **length below 24 months, height at 24 months and above**.\n",
    "Where the recorded position differs from the rule, the measurement is converted\n",
    "by about 0.7 cm, which is the systematic difference between recumbent length and\n",
    "standing height for the same child."
   ],
   "id": "cell-006"
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plausible[\"standard\"] = np.where(plausible[\"age_months\"] < 24, \"L\", \"H\")\n",
    "\n",
    "should_be_length = (plausible[\"standard\"] == \"L\") & ~plausible[\"measured_lying\"]\n",
    "should_be_height = (plausible[\"standard\"] == \"H\") & plausible[\"measured_lying\"]\n",
    "\n",
    "print(f\"measured standing, should be length: {int(should_be_length.sum())}\")\n",
    "print(f\"measured lying, should be height   : {int(should_be_height.sum())}\")\n",
    "\n",
    "plausible[\"lenhei\"] = plausible[\"height_cm\"]\n",
    "plausible.loc[should_be_length, \"lenhei\"] = plausible.loc[should_be_length, \"height_cm\"] + 0.7\n",
    "plausible.loc[should_be_height, \"lenhei\"] = plausible.loc[should_be_height, \"height_cm\"] - 0.7"
   ],
   "id": "cell-007"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Skip this and you bias every z-score for the younger half of the sample, in a\n",
    "direction that depends on how the teams happened to work.\n",
    "\n",
    "## The LMS calculation"
   ],
   "id": "cell-008"
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plausible[\"lenhei_key\"] = (plausible[\"lenhei\"] * 10).round() / 10\n",
    "\n",
    "lms = reference.set_index([\"sex\", \"lorh\", \"lenhei\"])[[\"l\", \"m\", \"s\"]]\n",
    "scored = plausible.join(lms, on=[\"sex\", \"standard\", \"lenhei_key\"])\n",
    "\n",
    "# The LMS transform: z = ((weight/M)^L - 1) / (L * S)\n",
    "raw_z = ((scored[\"weight_kg\"] / scored[\"m\"]) ** scored[\"l\"] - 1) / (\n",
    "    scored[\"l\"] * scored[\"s\"]\n",
    ")"
   ],
   "id": "cell-009"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Beyond ±3 the LMS curve is extrapolated, so WHO replaces it with a linear\n",
    "extension anchored on the distance between the 2nd and 3rd standard deviations.\n",
    "Without this, extreme children get z-scores that are too extreme, and a SAM\n",
    "prevalence built on them is overstated."
   ],
   "id": "cell-010"
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def sd_at(row, n):\n",
    "    return row[\"m\"] * (1 + row[\"l\"] * row[\"s\"] * n) ** (1 / row[\"l\"])\n",
    "\n",
    "def who_adjust(row, z):\n",
    "    if pd.isna(z):\n",
    "        return z\n",
    "    if z > 3:\n",
    "        sd3, sd2 = sd_at(row, 3), sd_at(row, 2)\n",
    "        return 3 + (row[\"weight_kg\"] - sd3) / (sd3 - sd2)\n",
    "    if z < -3:\n",
    "        sd3, sd2 = sd_at(row, -3), sd_at(row, -2)\n",
    "        return -3 + (row[\"weight_kg\"] - sd3) / (sd2 - sd3)\n",
    "    return z\n",
    "\n",
    "scored[\"whz\"] = [\n",
    "    who_adjust(row, z) for row, z in zip(scored.to_dict(\"records\"), raw_z)\n",
    "]\n",
    "scored[\"whz\"].describe().round(3)"
   ],
   "id": "cell-011"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Two flagging rules\n",
    "\n",
    "WHO flags are fixed bounds: a weight-for-height z-score outside -5 to +5 is\n",
    "biologically implausible. SMART flags are relative — more than 3 SD from the\n",
    "*survey* mean. They exclude different children, and a plausibility report states\n",
    "which was used."
   ],
   "id": "cell-012"
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "scored[\"who_flag\"] = scored[\"whz\"].abs() > 5\n",
    "\n",
    "mean_z, sd_z = scored[\"whz\"].mean(), scored[\"whz\"].std()\n",
    "scored[\"smart_flag\"] = (scored[\"whz\"] - mean_z).abs() > 3 * sd_z\n",
    "\n",
    "print(f\"WHO flagged  : {int(scored['who_flag'].sum())}\")\n",
    "print(f\"SMART flagged: {int(scored['smart_flag'].sum())}\")"
   ],
   "id": "cell-013"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The SMART rule is relative to a mean the flagged observations themselves\n",
    "influence, which is why the range check has to come first.\n",
    "\n",
    "## Prevalence\n",
    "\n",
    "**Oedema overrides anthropometry.** A child with bilateral pitting oedema is\n",
    "severely acutely malnourished whatever their weight-for-height, so the SAM\n",
    "numerator is not simply the count below -3 z-scores."
   ],
   "id": "cell-014"
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "analysable = scored[~scored[\"who_flag\"] & scored[\"whz\"].notna()]\n",
    "\n",
    "gam = ((analysable[\"whz\"] < -2) | analysable[\"oedema\"]).mean()\n",
    "sam = ((analysable[\"whz\"] < -3) | analysable[\"oedema\"]).mean()\n",
    "\n",
    "print(f\"analysable children: {len(analysable)}\")\n",
    "print(f\"GAM: {gam:.1%}\")\n",
    "print(f\"SAM: {sam:.1%}\")\n",
    "print(f\"mean z: {analysable['whz'].mean():.2f}   SD: {analysable['whz'].std():.2f}\")"
   ],
   "id": "cell-015"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Global acute malnutrition near 14.9%, severe near 3.9%. That sits just under the\n",
    "15% WHO emergency threshold — which is exactly where the decisions above stop\n",
    "being academic, because a skipped position adjustment or a different flagging\n",
    "rule moves the figure across the line.\n",
    "\n",
    "The standard deviation of the z-score is a quality signal in its own right. SMART\n",
    "expects it between about 0.8 and 1.2; above that suggests measurement error\n",
    "inflating the spread. This survey sits at the top of that range, for a reason the\n",
    "next section identifies.\n",
    "\n",
    "## The team effect, which is not a nutrition finding"
   ],
   "id": "cell-016"
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "by_team = analysable.groupby(\"team\").agg(\n",
    "    children=(\"whz\", \"size\"),\n",
    "    mean_z=(\"whz\", \"mean\"),\n",
    ")\n",
    "by_team[\"gam\"] = analysable.groupby(\"team\").apply(\n",
    "    lambda g: ((g[\"whz\"] < -2) | g[\"oedema\"]).mean(), include_groups=False\n",
    ")\n",
    "by_team.round(3)"
   ],
   "id": "cell-017"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Team 3 reports GAM near 22% against 10 to 16% for the others, with a mean z-score\n",
    "of -1.10 against -0.43 to -0.69. A real difference of that size between randomly\n",
    "assigned clusters would be extraordinary. It is a measurement artefact, and\n",
    "reporting it as a geographic finding would send resources to the wrong clusters.\n",
    "\n",
    "## The design effect\n",
    "\n",
    "This is a cluster sample, so children within a cluster resemble each other and\n",
    "the effective sample size is smaller than the count of children."
   ],
   "id": "cell-018"
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "analysable = analysable.assign(case=(analysable[\"whz\"] < -2) | analysable[\"oedema\"])\n",
    "\n",
    "clusters = analysable.groupby(\"cluster\").agg(m=(\"case\", \"size\"), y=(\"case\", \"sum\"))\n",
    "k, M = len(clusters), clusters[\"m\"].sum()\n",
    "p_bar = clusters[\"y\"].sum() / M\n",
    "\n",
    "# Ultimate-cluster variance of a ratio estimator: deviation of each cluster's\n",
    "# case count from what the overall rate predicts for its size.\n",
    "var_cluster = (k / ((k - 1) * M**2)) * ((clusters[\"y\"] - p_bar * clusters[\"m\"]) ** 2).sum()\n",
    "var_srs = p_bar * (1 - p_bar) / M\n",
    "\n",
    "deff = var_cluster / var_srs\n",
    "icc = (deff - 1) / (clusters[\"m\"].mean() - 1)\n",
    "\n",
    "print(f\"clusters: {k}   children: {M}   mean cluster size: {clusters['m'].mean():.1f}\")\n",
    "print(f\"design effect: {deff:.2f}   ICC: {icc:.3f}\")\n",
    "print(f\"effective sample size: {M / deff:.0f} of {M}\")"
   ],
   "id": "cell-019"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A design effect near 2.3 and an ICC around 0.05 are ordinary for a nutrition\n",
    "cluster survey. A DEFF below 1 or above about 4 usually means the calculation is\n",
    "wrong rather than the survey unusual — check your own arithmetic before you\n",
    "report it."
   ],
   "id": "cell-020"
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "se_cluster, se_srs = np.sqrt(var_cluster), np.sqrt(var_srs)\n",
    "\n",
    "print(f\"GAM {p_bar:.1%}  (95% CI {p_bar - 1.96*se_cluster:.1%} - {p_bar + 1.96*se_cluster:.1%})  with clustering\")\n",
    "print(f\"          (95% CI {p_bar - 1.96*se_srs:.1%} - {p_bar + 1.96*se_srs:.1%})  ignoring clustering\")"
   ],
   "id": "cell-021"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The correct interval is wider, and its upper bound crosses 15%. The point estimate\n",
    "sits below the WHO emergency threshold; the interval does not rule out being above\n",
    "it. That is the sentence the report needs, not a bare \"14.9%, below the threshold\".\n",
    "\n",
    "## What to report\n",
    "\n",
    "Prevalence with its interval and the design effect used, the flagging rule named,\n",
    "the exclusions counted, and the team comparison — because a survey where one team\n",
    "differs from the others by half a z-score has a measurement problem that outranks\n",
    "every prevalence figure in it."
   ],
   "id": "cell-022"
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "name": "python",
   "file_extension": ".py"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}