{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Computing prevalence with WHO growth standards\n",
    "\n",
    "*SMART nutrition survey, 2024 · R*\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.r.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, and GAM and SAM\n",
    "prevalence with a design effect for the cluster sample.\n",
    "\n",
    "Z-scores are deliberately **not** shipped in this dataset. Computing them is the\n",
    "exercise, and reading a precomputed column teaches nothing about the three\n",
    "decisions that move the answer: the range check, the measurement position, and\n",
    "which flagging rule you apply.\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",
    "## Use the official package, not your own LMS code\n",
    "\n",
    "The WHO standards are an LMS table, and reimplementing the interpolation is a\n",
    "well-known source of quiet error. `anthro` is maintained by WHO and applies the\n",
    "length/height adjustment itself."
   ],
   "id": "cell-001"
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#| message: false\n",
    "# install.packages(\"anthro\")\n",
    "library(readr)\n",
    "library(dplyr)\n",
    "library(anthro)\n",
    "\n",
    "URL <- paste0(\n",
    "  \"https://data-analysis.cassion.dev/datasets/files/\",\n",
    "  \"smart-nutrition-survey-2024.v1.csv\"\n",
    ")\n",
    "\n",
    "smart <- read_csv(URL, col_types = cols(\n",
    "  child_id = col_character(),\n",
    "  .default = col_guess()\n",
    "))\n",
    "\n",
    "glimpse(smart)"
   ],
   "id": "cell-002"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Range-check before flagging\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 not outliers to be\n",
    "flagged statistically; they are data entry errors, and they must be excluded\n",
    "before any flagging rule is applied, because a single 114 kg child moves the\n",
    "survey mean that the SMART flag is computed against."
   ],
   "id": "cell-003"
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "impossible <- smart |>\n",
    "  filter(weight_kg < 2 | weight_kg > 30 | height_cm < 45 | height_cm > 130)\n",
    "\n",
    "impossible |> select(child_id, team, weight_kg, height_cm)"
   ],
   "id": "cell-004"
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plausible <- smart |>\n",
    "  filter(\n",
    "    is.na(weight_kg) | is.na(height_cm) |\n",
    "      !(weight_kg < 2 | weight_kg > 30 | height_cm < 45 | height_cm > 130)\n",
    "  )\n",
    "\n",
    "cat(\"kept:\", nrow(plausible), \"of\", nrow(smart), \"\\n\")\n",
    "cat(\"missing age:\", sum(is.na(plausible$age_months)),\n",
    "    \" missing weight:\", sum(is.na(plausible$weight_kg)), \"\\n\")"
   ],
   "id": "cell-005"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Measurement position is not optional\n",
    "\n",
    "Children under two years are measured lying — recumbent length — and older\n",
    "children standing. Length reads about 0.7 cm greater than height for the same\n",
    "child, so mixing the two without adjustment biases every z-score for the younger\n",
    "half of the sample.\n",
    "\n",
    "Do **not** adjust the column by hand. Pass the position to `anthro_zscores` via\n",
    "`measure` and let it apply the WHO rule, which depends on age as well as\n",
    "position."
   ],
   "id": "cell-006"
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "z <- anthro_zscores(\n",
    "  sex             = ifelse(plausible$sex == \"m\", 1, 2),\n",
    "  age             = plausible$age_months,\n",
    "  is_age_in_month = TRUE,\n",
    "  weight          = plausible$weight_kg,\n",
    "  lenhei          = plausible$height_cm,\n",
    "  measure         = ifelse(plausible$measured_lying, \"l\", \"h\")\n",
    ")\n",
    "\n",
    "scored <- plausible |>\n",
    "  mutate(whz = z$zwfl, who_flag = z$fwfl)\n",
    "\n",
    "table(scored$who_flag, useNA = \"ifany\")"
   ],
   "id": "cell-007"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Two flagging rules, two slightly different surveys\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-008"
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mean_z <- mean(scored$whz, na.rm = TRUE)\n",
    "sd_z   <- sd(scored$whz, na.rm = TRUE)\n",
    "\n",
    "scored <- scored |>\n",
    "  mutate(smart_flag = abs(whz - mean_z) > 3 * sd_z)\n",
    "\n",
    "scored |>\n",
    "  summarise(\n",
    "    who_flagged   = sum(who_flag == 1, na.rm = TRUE),\n",
    "    smart_flagged = sum(smart_flag, na.rm = TRUE),\n",
    "    both          = sum(who_flag == 1 & smart_flag, na.rm = TRUE)\n",
    "  )"
   ],
   "id": "cell-009"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The SMART rule is relative to a mean that 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-010"
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "analysable <- scored |> filter(who_flag == 0, !is.na(whz))\n",
    "\n",
    "prevalence <- analysable |>\n",
    "  summarise(\n",
    "    children = n(),\n",
    "    gam = mean(whz < -2 | oedema),\n",
    "    sam = mean(whz < -3 | oedema),\n",
    "    mean_z = mean(whz),\n",
    "    sd_z   = sd(whz)\n",
    "  )\n",
    "\n",
    "prevalence |> mutate(across(c(gam, sam), ~ round(100 * .x, 1)),\n",
    "                     across(c(mean_z, sd_z), ~ round(.x, 2)))"
   ],
   "id": "cell-011"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Global acute malnutrition near 14.9% and severe near 3.9%. That sits just under\n",
    "the 15% WHO emergency threshold — which is exactly the position where the\n",
    "analytical choices above stop being academic, because a different flagging rule\n",
    "or a skipped position adjustment moves the figure across the line.\n",
    "\n",
    "Note the standard deviation of the z-score. SMART expects it between about 0.8\n",
    "and 1.2; a value above that suggests measurement error inflating the spread, and\n",
    "this survey sits at the top of the acceptable range for a reason the next section\n",
    "identifies.\n",
    "\n",
    "## The team effect, which is not a nutrition finding"
   ],
   "id": "cell-012"
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "analysable |>\n",
    "  group_by(team) |>\n",
    "  summarise(\n",
    "    children = n(),\n",
    "    mean_z   = round(mean(whz), 2),\n",
    "    gam      = round(100 * mean(whz < -2 | oedema), 1),\n",
    "    .groups = \"drop\"\n",
    "  )"
   ],
   "id": "cell-013"
  },
  {
   "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.09 against -0.43 to -0.69. A real difference in nutrition status between\n",
    "randomly assigned clusters of that size would be extraordinary. This is a\n",
    "measurement artefact — a team measuring height long or weight light — 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. Children within a cluster resemble each other, so the\n",
    "effective sample size is smaller than the count of children and a confidence\n",
    "interval computed as though the sample were simple random understates itself."
   ],
   "id": "cell-014"
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "clusters <- analysable |>\n",
    "  mutate(case = whz < -2 | oedema) |>\n",
    "  group_by(cluster) |>\n",
    "  summarise(m = n(), y = sum(case), .groups = \"drop\")\n",
    "\n",
    "k     <- nrow(clusters)\n",
    "M     <- sum(clusters$m)\n",
    "p_bar <- sum(clusters$y) / M\n",
    "\n",
    "# Ultimate-cluster variance of a ratio estimator. Note it uses the deviation of\n",
    "# each cluster's case count from what the overall rate predicts for its size —\n",
    "# not the variance of the cluster rates, which ignores that clusters differ in\n",
    "# size.\n",
    "var_cluster <- (k / ((k - 1) * M^2)) * sum((clusters$y - p_bar * clusters$m)^2)\n",
    "var_srs     <- p_bar * (1 - p_bar) / M\n",
    "\n",
    "deff <- var_cluster / var_srs\n",
    "icc  <- (deff - 1) / (mean(clusters$m) - 1)\n",
    "\n",
    "cat(sprintf(\"clusters: %d   children: %d   mean cluster size: %.1f\\n\",\n",
    "            k, M, mean(clusters$m)))\n",
    "cat(sprintf(\"design effect: %.2f   ICC: %.3f\\n\", deff, icc))\n",
    "cat(sprintf(\"effective sample size: %.0f of %d\\n\", M / deff, M))"
   ],
   "id": "cell-015"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A design effect near 2.3 and an ICC around 0.045 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 — this is a check worth running on your own\n",
    "arithmetic before you report it."
   ],
   "id": "cell-016"
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "se_cluster <- sqrt(var_cluster)\n",
    "se_srs     <- sqrt(var_srs)\n",
    "\n",
    "cat(sprintf(\"GAM %.1f%%  (95%% CI %.1f - %.1f)  accounting for clustering\\n\",\n",
    "            100 * p_bar, 100 * (p_bar - 1.96 * se_cluster), 100 * (p_bar + 1.96 * se_cluster)))\n",
    "cat(sprintf(\"           (95%% CI %.1f - %.1f)  if clustering is ignored\\n\",\n",
    "            100 * (p_bar - 1.96 * se_srs), 100 * (p_bar + 1.96 * se_srs)))"
   ],
   "id": "cell-017"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Two things to read there. The correct interval is wider, and a report that ignores\n",
    "clustering claims precision the design cannot deliver — the most common way a\n",
    "survey overstates what it knows.\n",
    "\n",
    "And the upper bound crosses 15%. The point estimate sits below the WHO emergency\n",
    "threshold; the interval does not rule out being above it. That is the sentence\n",
    "the report needs, not a bare \"14.9%, below the emergency 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 the report."
   ],
   "id": "cell-018"
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "R",
   "language": "R",
   "name": "ir"
  },
  "language_info": {
   "name": "R",
   "file_extension": ".r"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
