{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Global acute malnutrition by commune\n",
    "\n",
    "*MUAC screening — Artibonite, 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/muac-screening-artibonite-2024/gam-by-commune.r.en.ipynb)"
   ],
   "id": "cell-000"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## What this produces\n",
    "\n",
    "The same analysis as the Python example, in R. Python and R are peers on this\n",
    "platform, and a technique that exists in only one of them is incomplete — you\n",
    "will inherit whichever your predecessor used.\n",
    "\n",
    "GAM and SAM prevalence by commune with 95% confidence intervals. The reference\n",
    "figures are in the dataset's quality notes: overall GAM near 8.6% and SAM near\n",
    "2.2%, ranging from roughly 5% to 15% by commune. If your numbers land far\n",
    "outside that, you have a bug rather 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"
   ],
   "id": "cell-001"
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#| message: false\n",
    "library(readr)\n",
    "library(dplyr)\n",
    "library(ggplot2)\n",
    "\n",
    "URL <- paste0(\n",
    "  \"https://data-analysis.cassion.dev/datasets/files/\",\n",
    "  \"muac-screening-artibonite-2024.v1.csv\"\n",
    ")\n",
    "\n",
    "muac <- read_csv(\n",
    "  URL,\n",
    "  col_types = cols(\n",
    "    child_id       = col_character(),\n",
    "    commune        = col_character(),\n",
    "    screening_date = col_date(format = \"%Y-%m-%d\"),\n",
    "    age_months     = col_integer(),\n",
    "    sex            = col_character(),\n",
    "    muac_mm        = col_integer(),\n",
    "    oedema         = col_character(),\n",
    "    outcome        = col_character()\n",
    "  ),\n",
    "  na = c(\"\", \"NA\", \"-99\")\n",
    ")\n",
    "\n",
    "dim(muac)\n",
    "glimpse(muac)"
   ],
   "id": "cell-002"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note the difference from pandas. The `na` argument of `read_csv` applies across\n",
    "the whole file rather than to one column, so if another column ever held a\n",
    "legitimate `-99` you would scope the treatment afterwards. Declaring `col_types`\n",
    "explicitly is what stops `child_id` being read as a number and losing a leading\n",
    "zero.\n",
    "\n",
    "## Clean what would corrupt the indicator\n",
    "\n",
    "Seven records were left in centimetres and never converted. The plausible\n",
    "millimetre and centimetre ranges do not overlap, so the correction is\n",
    "unambiguous. Two communes recorded oedema as `Y`/`N` in the first quarter."
   ],
   "id": "cell-003"
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n_unit_error <- sum(!is.na(muac$muac_mm) & muac$muac_mm < 40)\n",
    "\n",
    "muac <- muac |>\n",
    "  mutate(\n",
    "    muac_mm = if_else(!is.na(muac_mm) & muac_mm < 40, muac_mm * 10L, muac_mm),\n",
    "    muac_mm = if_else(\n",
    "      !is.na(muac_mm) & (muac_mm < 80 | muac_mm > 220),\n",
    "      NA_integer_,\n",
    "      muac_mm\n",
    "    ),\n",
    "    oedema = case_when(\n",
    "      tolower(trimws(oedema)) %in% c(\"true\", \"y\", \"yes\") ~ TRUE,\n",
    "      tolower(trimws(oedema)) %in% c(\"false\", \"n\", \"no\") ~ FALSE,\n",
    "      TRUE ~ NA\n",
    "    )\n",
    "  )\n",
    "\n",
    "cat(\"unit errors corrected :\", n_unit_error, \"\\n\")\n",
    "cat(\"oedema still missing  :\", sum(is.na(muac$oedema)), \"\\n\")"
   ],
   "id": "cell-004"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Define the indicator before computing it\n",
    "\n",
    "- **Numerator** — MUAC below 125 mm, or bilateral pitting oedema.\n",
    "- **Denominator** — a valid MUAC measurement **or** a recorded 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.\n",
    "\n",
    "`oedema %in% TRUE` rather than a bare `oedema` test: a missing oedema assessment\n",
    "is not an absent oedema, and `NA | FALSE` is `NA` in R, which would silently\n",
    "drop the row from the numerator."
   ],
   "id": "cell-005"
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "SAM_MM <- 115\n",
    "GAM_MM <- 125\n",
    "\n",
    "muac <- muac |>\n",
    "  mutate(\n",
    "    has_assessment = !is.na(muac_mm) | !is.na(oedema),\n",
    "    sam = if_else(has_assessment, (muac_mm < SAM_MM) | oedema %in% TRUE, NA),\n",
    "    gam = if_else(has_assessment, (muac_mm < GAM_MM) | oedema %in% TRUE, NA)\n",
    "  )"
   ],
   "id": "cell-006"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The indicator table, with exact intervals\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.\n",
    "\n",
    "`binom.test` gives the Clopper-Pearson interval, which is exact for a proportion\n",
    "and behaves at the small counts the less populous communes produce."
   ],
   "id": "cell-007"
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ci <- function(cases, n) {\n",
    "  if (n == 0 || is.na(cases)) return(c(NA_real_, NA_real_))\n",
    "  binom.test(cases, n)$conf.int\n",
    "}\n",
    "\n",
    "table <- muac |>\n",
    "  group_by(commune) |>\n",
    "  summarise(\n",
    "    screened    = n(),\n",
    "    denominator = sum(has_assessment),\n",
    "    sam_cases   = sum(sam, na.rm = TRUE),\n",
    "    gam_cases   = sum(gam, na.rm = TRUE),\n",
    "    .groups = \"drop\"\n",
    "  ) |>\n",
    "  mutate(\n",
    "    gam_rate = gam_cases / denominator,\n",
    "    sam_rate = sam_cases / denominator\n",
    "  ) |>\n",
    "  rowwise() |>\n",
    "  mutate(\n",
    "    gam_low  = ci(gam_cases, denominator)[1],\n",
    "    gam_high = ci(gam_cases, denominator)[2]\n",
    "  ) |>\n",
    "  ungroup() |>\n",
    "  arrange(desc(gam_rate))\n",
    "\n",
    "table"
   ],
   "id": "cell-008"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plot the ranking with its uncertainty"
   ],
   "id": "cell-009"
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#| fig-width: 8\n",
    "#| fig-height: 5\n",
    "EMERGENCY_THRESHOLD <- 0.15\n",
    "\n",
    "ggplot(table, aes(x = gam_rate * 100, y = reorder(commune, gam_rate))) +\n",
    "  geom_vline(\n",
    "    xintercept = EMERGENCY_THRESHOLD * 100,\n",
    "    colour = \"#B5533C\", linetype = \"dashed\"\n",
    "  ) +\n",
    "  geom_errorbarh(\n",
    "    aes(xmin = gam_low * 100, xmax = gam_high * 100),\n",
    "    height = 0, colour = \"#9AA8A3\"\n",
    "  ) +\n",
    "  geom_point(colour = \"#2F5D50\", size = 2) +\n",
    "  annotate(\n",
    "    \"text\", x = EMERGENCY_THRESHOLD * 100 + 0.3, y = 1,\n",
    "    label = \"emergency threshold (15%)\", colour = \"#B5533C\",\n",
    "    hjust = 0, size = 3\n",
    "  ) +\n",
    "  labs(\n",
    "    x = \"GAM prevalence by MUAC (%)\",\n",
    "    y = NULL,\n",
    "    title = \"Global acute malnutrition by commune, with 95% CI\"\n",
    "  ) +\n",
    "  theme_minimal(base_size = 11) +\n",
    "  theme(panel.grid.minor = element_blank())"
   ],
   "id": "cell-010"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Reading the result\n",
    "\n",
    "Several intervals overlap, so the ordering between those communes is not\n",
    "supported by the data. A decision that allocates one additional CMAM site by\n",
    "rank alone is reading precision the screening does not have.\n",
    "\n",
    "The interval also describes sampling variation only. This is a census of the\n",
    "children who came to be screened rather than a probability sample, so whether\n",
    "those children resemble the ones who did not come is a separate — and usually\n",
    "larger — source of error, belonging in the limitations section of any report\n",
    "built on this."
   ],
   "id": "cell-011"
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "R",
   "language": "R",
   "name": "ir"
  },
  "language_info": {
   "name": "R",
   "file_extension": ".r"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}