---
title: "Computing prevalence with WHO growth standards"
subtitle: "SMART nutrition survey, 2024"
author: "Cassion · data-analysis.cassion.dev"
format:
  html:
    toc: true
    code-fold: false
jupyter: python3
---

## What this produces

Weight-for-height z-scores against the WHO 2006 growth standards, computed from
the LMS reference rather than read from a package — because the three decisions
that move the answer are all outside the package call.

Z-scores are deliberately not shipped in this dataset. Computing them is the
exercise.

Every dataset on this platform is synthetic. No real child is described, and
these results must not be cited as a real nutrition situation.

## The reference table

R has the official WHO `anthro` package. Python's options are thinner, so this
example reads the WHO 2006 weight-for-length and weight-for-height LMS tables
directly. They ship with the platform, and the values are the same ones the WHO
package uses — the R version of this example produces identical prevalence to the
first decimal.

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

BASE = "https://data-analysis.cassion.dev/datasets/"

smart = pd.read_csv(BASE + "files/smart-nutrition-survey-2024.v1.csv",
                    dtype={"child_id": "string", "sex": "string"})
reference = pd.read_csv(BASE + "reference/who-2006-weight-for-lenhei.csv")

print(f"survey: {len(smart)} children")
print(f"reference: {len(reference)} rows")
reference.head()
```

`lorh` is the standard the row belongs to: `L` for weight-for-length, used for
children under two, and `H` for weight-for-height above that. `l`, `m` and `s`
are the LMS parameters at that length or height.

## Range-check before anything else

Fourteen records hold impossible measurements — weights out by a factor of ten in
both directions, and heights entered in metres. These are data entry errors, not
statistical outliers, and they must go before any flagging rule, because a single
114 kg child moves the survey mean that the SMART flag is measured against.

```{python}
impossible = (
    (smart["weight_kg"] < 2) | (smart["weight_kg"] > 30)
    | (smart["height_cm"] < 45) | (smart["height_cm"] > 130)
)
smart.loc[impossible, ["child_id", "team", "weight_kg", "height_cm"]]
```

```{python}
plausible = smart[~impossible.fillna(False)].copy()
print(f"kept {len(plausible)} of {len(smart)}")
print(f"missing age: {int(plausible['age_months'].isna().sum())}, "
      f"missing weight: {int(plausible['weight_kg'].isna().sum())}")
```

## Measurement position, which is a rule about age

This is the step people skip. The WHO standard is not "use whatever position was
measured" — it is **length below 24 months, height at 24 months and above**.
Where the recorded position differs from the rule, the measurement is converted
by about 0.7 cm, which is the systematic difference between recumbent length and
standing height for the same child.

```{python}
plausible["standard"] = np.where(plausible["age_months"] < 24, "L", "H")

should_be_length = (plausible["standard"] == "L") & ~plausible["measured_lying"]
should_be_height = (plausible["standard"] == "H") & plausible["measured_lying"]

print(f"measured standing, should be length: {int(should_be_length.sum())}")
print(f"measured lying, should be height   : {int(should_be_height.sum())}")

plausible["lenhei"] = plausible["height_cm"]
plausible.loc[should_be_length, "lenhei"] = plausible.loc[should_be_length, "height_cm"] + 0.7
plausible.loc[should_be_height, "lenhei"] = plausible.loc[should_be_height, "height_cm"] - 0.7
```

Skip this and you bias every z-score for the younger half of the sample, in a
direction that depends on how the teams happened to work.

## The LMS calculation

```{python}
plausible["lenhei_key"] = (plausible["lenhei"] * 10).round() / 10

lms = reference.set_index(["sex", "lorh", "lenhei"])[["l", "m", "s"]]
scored = plausible.join(lms, on=["sex", "standard", "lenhei_key"])

# The LMS transform: z = ((weight/M)^L - 1) / (L * S)
raw_z = ((scored["weight_kg"] / scored["m"]) ** scored["l"] - 1) / (
    scored["l"] * scored["s"]
)
```

Beyond ±3 the LMS curve is extrapolated, so WHO replaces it with a linear
extension anchored on the distance between the 2nd and 3rd standard deviations.
Without this, extreme children get z-scores that are too extreme, and a SAM
prevalence built on them is overstated.

```{python}
def sd_at(row, n):
    return row["m"] * (1 + row["l"] * row["s"] * n) ** (1 / row["l"])

def who_adjust(row, z):
    if pd.isna(z):
        return z
    if z > 3:
        sd3, sd2 = sd_at(row, 3), sd_at(row, 2)
        return 3 + (row["weight_kg"] - sd3) / (sd3 - sd2)
    if z < -3:
        sd3, sd2 = sd_at(row, -3), sd_at(row, -2)
        return -3 + (row["weight_kg"] - sd3) / (sd2 - sd3)
    return z

scored["whz"] = [
    who_adjust(row, z) for row, z in zip(scored.to_dict("records"), raw_z)
]
scored["whz"].describe().round(3)
```

## Two flagging rules

WHO flags are fixed bounds: a weight-for-height z-score outside -5 to +5 is
biologically implausible. SMART flags are relative — more than 3 SD from the
*survey* mean. They exclude different children, and a plausibility report states
which was used.

```{python}
scored["who_flag"] = scored["whz"].abs() > 5

mean_z, sd_z = scored["whz"].mean(), scored["whz"].std()
scored["smart_flag"] = (scored["whz"] - mean_z).abs() > 3 * sd_z

print(f"WHO flagged  : {int(scored['who_flag'].sum())}")
print(f"SMART flagged: {int(scored['smart_flag'].sum())}")
```

The SMART rule is relative to a mean the flagged observations themselves
influence, which is why the range check has to come first.

## Prevalence

**Oedema overrides anthropometry.** A child with bilateral pitting oedema is
severely acutely malnourished whatever their weight-for-height, so the SAM
numerator is not simply the count below -3 z-scores.

```{python}
analysable = scored[~scored["who_flag"] & scored["whz"].notna()]

gam = ((analysable["whz"] < -2) | analysable["oedema"]).mean()
sam = ((analysable["whz"] < -3) | analysable["oedema"]).mean()

print(f"analysable children: {len(analysable)}")
print(f"GAM: {gam:.1%}")
print(f"SAM: {sam:.1%}")
print(f"mean z: {analysable['whz'].mean():.2f}   SD: {analysable['whz'].std():.2f}")
```

Global acute malnutrition near 14.9%, severe near 3.9%. That sits just under the
15% WHO emergency threshold — which is exactly where the decisions above stop
being academic, because a skipped position adjustment or a different flagging
rule moves the figure across the line.

The standard deviation of the z-score is a quality signal in its own right. SMART
expects it between about 0.8 and 1.2; above that suggests measurement error
inflating the spread. This survey sits at the top of that range, for a reason the
next section identifies.

## The team effect, which is not a nutrition finding

```{python}
by_team = analysable.groupby("team").agg(
    children=("whz", "size"),
    mean_z=("whz", "mean"),
)
by_team["gam"] = analysable.groupby("team").apply(
    lambda g: ((g["whz"] < -2) | g["oedema"]).mean(), include_groups=False
)
by_team.round(3)
```

Team 3 reports GAM near 22% against 10 to 16% for the others, with a mean z-score
of -1.10 against -0.43 to -0.69. A real difference of that size between randomly
assigned clusters would be extraordinary. It is a measurement artefact, and
reporting it as a geographic finding would send resources to the wrong clusters.

## The design effect

This is a cluster sample, so children within a cluster resemble each other and
the effective sample size is smaller than the count of children.

```{python}
analysable = analysable.assign(case=(analysable["whz"] < -2) | analysable["oedema"])

clusters = analysable.groupby("cluster").agg(m=("case", "size"), y=("case", "sum"))
k, M = len(clusters), clusters["m"].sum()
p_bar = clusters["y"].sum() / M

# Ultimate-cluster variance of a ratio estimator: deviation of each cluster's
# case count from what the overall rate predicts for its size.
var_cluster = (k / ((k - 1) * M**2)) * ((clusters["y"] - p_bar * clusters["m"]) ** 2).sum()
var_srs = p_bar * (1 - p_bar) / M

deff = var_cluster / var_srs
icc = (deff - 1) / (clusters["m"].mean() - 1)

print(f"clusters: {k}   children: {M}   mean cluster size: {clusters['m'].mean():.1f}")
print(f"design effect: {deff:.2f}   ICC: {icc:.3f}")
print(f"effective sample size: {M / deff:.0f} of {M}")
```

A design effect near 2.3 and an ICC around 0.05 are ordinary for a nutrition
cluster survey. A DEFF below 1 or above about 4 usually means the calculation is
wrong rather than the survey unusual — check your own arithmetic before you
report it.

```{python}
se_cluster, se_srs = np.sqrt(var_cluster), np.sqrt(var_srs)

print(f"GAM {p_bar:.1%}  (95% CI {p_bar - 1.96*se_cluster:.1%} - {p_bar + 1.96*se_cluster:.1%})  with clustering")
print(f"          (95% CI {p_bar - 1.96*se_srs:.1%} - {p_bar + 1.96*se_srs:.1%})  ignoring clustering")
```

The correct interval is wider, and its upper bound crosses 15%. The point estimate
sits below the WHO emergency threshold; the interval does not rule out being above
it. That is the sentence the report needs, not a bare "14.9%, below the threshold".

## What to report

Prevalence with its interval and the design effect used, the flagging rule named,
the exclusions counted, and the team comparison — because a survey where one team
differs from the others by half a z-score has a measurement problem that outranks
every prevalence figure in it.