{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Accept, exclude, or repeat — analysis notebook\n",
    "\n",
    "*SMART survey analysis and plausibility report · project deliverable*\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/projects/smart-survey-analysis/notebooks/whz-prevalence.python.en.ipynb)"
   ],
   "id": "cell-000"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The decision this produces\n",
    "\n",
    "Whether to accept the survey, exclude a measurement team's data, or repeat the\n",
    "fieldwork — before the estimate informs a scale-up.\n",
    "\n",
    "The estimate came back close to the threshold that would trigger one, and one\n",
    "team's measurements looked different from the others. That is the whole\n",
    "question: if the team is right the scale-up is justified, and if the team is\n",
    "measuring badly it is not.\n",
    "\n",
    "Audience: the nutrition cluster and the national nutrition coordination body.\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",
    "## Setup"
   ],
   "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",
    "EMERGENCY = 0.15\n",
    "PLAUSIBLE_WEIGHT = (2, 30)\n",
    "PLAUSIBLE_HEIGHT = (45, 130)\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\"{len(smart)} children, {smart['cluster'].nunique()} clusters, \"\n",
    "      f\"{smart['team'].nunique()} teams\")"
   ],
   "id": "cell-002"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Z-scores, computed rather than read"
   ],
   "id": "cell-003"
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# `between()` returns False for NaN, so a bare `~between(...)` counts a missing\n",
    "# measurement as an impossible one. Those are different failures and only one of\n",
    "# them is a data-entry error.\n",
    "impossible = (\n",
    "    (smart[\"weight_kg\"].notna() & ~smart[\"weight_kg\"].between(*PLAUSIBLE_WEIGHT))\n",
    "    | (smart[\"height_cm\"].notna() & ~smart[\"height_cm\"].between(*PLAUSIBLE_HEIGHT))\n",
    ")\n",
    "plausible = smart[~impossible].copy()\n",
    "\n",
    "print(f\"impossible measurements excluded: {int(impossible.sum())}\")\n",
    "print(f\"missing weight (kept, scored as NA): {int(smart['weight_kg'].isna().sum())}\")\n",
    "print(f\"missing age (kept, scored as NA)   : {int(smart['age_months'].isna().sum())}\")\n",
    "\n",
    "# WHO uses length below 24 months and height at 24 months and above; where the\n",
    "# recorded position differs, the measurement is converted by 0.7 cm.\n",
    "plausible[\"standard\"] = np.where(plausible[\"age_months\"] < 24, \"L\", \"H\")\n",
    "to_length = (plausible[\"standard\"] == \"L\") & ~plausible[\"measured_lying\"]\n",
    "to_height = (plausible[\"standard\"] == \"H\") & plausible[\"measured_lying\"]\n",
    "\n",
    "plausible[\"lenhei\"] = plausible[\"height_cm\"]\n",
    "plausible.loc[to_length, \"lenhei\"] += 0.7\n",
    "plausible.loc[to_height, \"lenhei\"] -= 0.7\n",
    "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",
    "raw = ((scored[\"weight_kg\"] / scored[\"m\"]) ** scored[\"l\"] - 1) / (scored[\"l\"] * scored[\"s\"])\n",
    "\n",
    "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\"] = [who_adjust(r, z) for r, z in zip(scored.to_dict(\"records\"), raw)]\n",
    "analysable = scored[(scored[\"whz\"].abs() <= 5) & scored[\"whz\"].notna()]\n",
    "len(analysable)"
   ],
   "id": "cell-004"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The estimate, and where it sits against the threshold"
   ],
   "id": "cell-005"
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def prevalence(df):\n",
    "    return pd.Series({\n",
    "        \"children\": len(df),\n",
    "        \"GAM\": ((df[\"whz\"] < -2) | df[\"oedema\"]).mean(),\n",
    "        \"SAM\": ((df[\"whz\"] < -3) | df[\"oedema\"]).mean(),\n",
    "        \"mean z\": df[\"whz\"].mean(),\n",
    "        \"SD z\": df[\"whz\"].std(),\n",
    "    })\n",
    "\n",
    "all_teams = prevalence(analysable)\n",
    "all_teams.round(3)"
   ],
   "id": "cell-006"
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "clusters = analysable.assign(\n",
    "    case=(analysable[\"whz\"] < -2) | analysable[\"oedema\"]\n",
    ").groupby(\"cluster\").agg(m=(\"case\", \"size\"), y=(\"case\", \"sum\"))\n",
    "\n",
    "def with_interval(cl):\n",
    "    k, M = len(cl), cl[\"m\"].sum()\n",
    "    p = cl[\"y\"].sum() / M\n",
    "    var_cluster = (k / ((k - 1) * M**2)) * ((cl[\"y\"] - p * cl[\"m\"]) ** 2).sum()\n",
    "    deff = var_cluster / (p * (1 - p) / M)\n",
    "    se = np.sqrt(var_cluster)\n",
    "    return p, p - 1.96 * se, p + 1.96 * se, deff\n",
    "\n",
    "p, low, high, deff = with_interval(clusters)\n",
    "print(f\"GAM {p:.1%}  (95% CI {low:.1%} - {high:.1%})   design effect {deff:.2f}\")\n",
    "print(f\"emergency threshold {EMERGENCY:.0%}: point estimate \"\n",
    "      f\"{'above' if p > EMERGENCY else 'below'}, interval \"\n",
    "      f\"{'crosses it' if high > EMERGENCY > low else 'does not cross it'}\")"
   ],
   "id": "cell-007"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The point estimate sits just below the threshold and the interval crosses it. On\n",
    "that alone the survey cannot say whether the scale-up criterion is met — which is\n",
    "why the team question has to be settled before the number is used at all.\n",
    "\n",
    "## The team question"
   ],
   "id": "cell-008"
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "by_team = analysable.groupby(\"team\").apply(prevalence, include_groups=False).round(3)\n",
    "by_team"
   ],
   "id": "cell-009"
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "spread = by_team[\"mean z\"].max() - by_team[\"mean z\"].min()\n",
    "suspect = by_team[\"mean z\"].idxmin()\n",
    "\n",
    "print(f\"spread in mean z across teams: {spread:.2f}\")\n",
    "print(f\"lowest-measuring team: {suspect}\")\n",
    "print(\"SMART treats a between-team spread above ~0.3 z as a supervision problem.\")"
   ],
   "id": "cell-010"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Clusters were assigned to teams independently of nutrition status, so a real\n",
    "difference of half a z-score between teams would be extraordinary. Two further\n",
    "checks say which explanation fits."
   ],
   "id": "cell-011"
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "scored[\"last_digit\"] = ((scored[\"height_cm\"] * 10).round() % 10).astype(\"Int64\")\n",
    "digits = pd.crosstab(scored[\"team\"], scored[\"last_digit\"], normalize=\"index\") * 100\n",
    "rounded = digits[[0, 5]].sum(axis=1).round(1)\n",
    "rounded.to_frame(\"% of heights on .0 or .5\").assign(expected=20.0)"
   ],
   "id": "cell-012"
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sizes = analysable.groupby(\"team\")[\"cluster\"].nunique()\n",
    "sizes.to_frame(\"clusters measured\")"
   ],
   "id": "cell-013"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Digit preference names a *different* team as its problem. So the survey has two\n",
    "independent measurement faults, in two teams, and only one of them moves the\n",
    "prevalence estimate.\n",
    "\n",
    "## The three options, costed"
   ],
   "id": "cell-014"
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "options = {}\n",
    "\n",
    "options[\"accept as measured\"] = with_interval(clusters)\n",
    "\n",
    "without = analysable[analysable[\"team\"] != suspect]\n",
    "cl_without = without.assign(\n",
    "    case=(without[\"whz\"] < -2) | without[\"oedema\"]\n",
    ").groupby(\"cluster\").agg(m=(\"case\", \"size\"), y=(\"case\", \"sum\"))\n",
    "options[f\"exclude team {suspect}\"] = with_interval(cl_without)\n",
    "\n",
    "summary = pd.DataFrame(\n",
    "    {k: {\"GAM\": v[0], \"CI low\": v[1], \"CI high\": v[2], \"design effect\": v[3]}\n",
    "     for k, v in options.items()}\n",
    ").T.round(3)\n",
    "summary[\"children\"] = [len(analysable), len(without)]\n",
    "summary[\"crosses threshold\"] = summary[\"CI high\"] > EMERGENCY\n",
    "summary"
   ],
   "id": "cell-015"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The recommendation"
   ],
   "id": "cell-016"
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(f\"as measured      : GAM {options['accept as measured'][0]:.1%}, \"\n",
    "      f\"n={len(analysable)}\")\n",
    "print(f\"excluding team {suspect} : GAM {options[f'exclude team {suspect}'][0]:.1%}, \"\n",
    "      f\"n={len(without)}\")\n",
    "print(f\"\\ndifference: \"\n",
    "      f\"{options[f'exclude team {suspect}'][0] - options['accept as measured'][0]:+.1%}\")"
   ],
   "id": "cell-017"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Do not accept as measured.** The estimate is inflated by a team whose mean\n",
    "z-score sits half a standard deviation below the others across independently\n",
    "assigned clusters, and the scale-up decision turns on exactly that difference.\n",
    "\n",
    "**Do not silently exclude the team either.** Excluding a quarter of clusters\n",
    "changes the sampling frame, and a prevalence computed on the remainder is no\n",
    "longer the survey that was designed. It is reported here as a sensitivity\n",
    "analysis, not as the answer.\n",
    "\n",
    "**Recommend re-measuring that team's clusters** if the teams are still in the\n",
    "field, and publishing with the team comparison in the body of the report rather\n",
    "than an annex if they are not. A reader shown only the district figure cannot\n",
    "know that a quarter of the sample was measured by someone whose results match\n",
    "nobody else's.\n",
    "\n",
    "## What this does not establish\n",
    "\n",
    "Which team is right. The comparison shows one team differs; it does not show the\n",
    "others are correct. If a standardisation exercise was run before fieldwork, its\n",
    "results settle this question and this notebook does not have them.\n",
    "\n",
    "Age heaping on whole years affects which growth standard applies, because the\n",
    "length/height rule switches at exactly 24 months. That is a second, smaller\n",
    "source of error running through every team, and it is in the plausibility report\n",
    "rather than here."
   ],
   "id": "cell-018"
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "name": "python",
   "file_extension": ".py"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}