---
title: "Water quality against chlorination"
subtitle: "WASH household survey, 2024"
author: "Cassion · data-analysis.cassion.dev"
format:
  html:
    toc: true
    code-fold: false
jupyter: python3
---

## The question

Does free residual chlorine at or above 0.2 mg/L predict a lower E. coli risk
class? And — the part that decides whether the answer means anything — how much
of the sample can actually answer it?

Every dataset on this platform is synthetic. No real household is described.

## Setup

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

URL = (
    "https://data-analysis.cassion.dev/datasets/files/"
    "wash-household-survey-2024.v1.csv"
)

wash = pd.read_csv(URL, dtype={"household_id": "string"})
wash["district"] = (
    wash["district"].str.strip().str.lower().str.replace(" ", "-", regex=False)
)
len(wash)
```

## Start with the denominator, not the finding

Water quality is tested on a subset. Reporting a quality result against the full
sample overstates how much testing was done, and it is the first thing a
reviewer will check.

```{python}
tested = pd.DataFrame({
    "households": [
        len(wash),
        int(wash["free_residual_chlorine_mgl"].notna().sum()),
        int(wash["ecoli_cfu_100ml"].notna().sum()),
        int(
            (
                wash["free_residual_chlorine_mgl"].notna()
                & wash["ecoli_cfu_100ml"].notna()
            ).sum()
        ),
    ]
}, index=["surveyed", "chlorine tested", "E. coli tested", "both tested"])
tested["share of sample"] = (tested["households"] / len(wash)).round(3)
tested
```

**Only the last row can answer the question.** Everything below is computed on
that subset, and the report has to say so — a chlorination finding presented
against 2,403 households when 320 were tested for both is a misstatement of the
evidence, whatever the finding is.

## The missingness is not random

About 530 households report treating their water and have no chlorine
measurement. That is a cross-field inconsistency, not a missing-at-random value,
and it matters: dropping those rows removes households that treat their water
more often than average, which biases the comparison in a predictable direction.

```{python}
treated_untested = wash["water_treated_at_home"] & wash[
    "free_residual_chlorine_mgl"
].isna()
print(f"treat water but no chlorine reading: {int(treated_untested.sum())}")

pd.crosstab(
    wash["water_treated_at_home"],
    wash["free_residual_chlorine_mgl"].notna().map(
        {True: "chlorine tested", False: "not tested"}
    ),
    normalize="index",
).round(3)
```

Households that treat their water are *more* likely to be tested here, so the
tested subset over-represents treatment. Any effect estimated on it is an effect
among the tested, not among the population.

## The comparison

WHO's guideline for free residual chlorine at the point of delivery is
0.2 mg/L. E. coli is reported in the standard risk classes rather than as a raw
count, because the count is over-dispersed and the classes are what a programme
acts on.

```{python}
both = wash.dropna(subset=["free_residual_chlorine_mgl", "ecoli_cfu_100ml"]).copy()

both["adequate_chlorine"] = both["free_residual_chlorine_mgl"] >= 0.2
both["risk_class"] = pd.cut(
    both["ecoli_cfu_100ml"],
    bins=[-1, 0, 10, 100, np.inf],
    labels=["0 (conforms)", "1-10 (low)", "11-100 (intermediate)", ">100 (high)"],
)

table = pd.crosstab(
    both["adequate_chlorine"],
    both["risk_class"],
    normalize="index",
).round(3)
table.index = ["below 0.2 mg/L", "at or above 0.2 mg/L"]
table
```

```{python}
counts = pd.crosstab(both["adequate_chlorine"], both["risk_class"])
counts.index = ["below 0.2 mg/L", "at or above 0.2 mg/L"]
counts
```

The pattern is strong: households at or above the guideline are far more likely
to show no detectable E. coli, and none of them reach the high-risk class.

## Is it more than sampling noise?

```{python}
from scipy.stats import chi2_contingency

conforms = pd.crosstab(
    both["adequate_chlorine"], both["ecoli_cfu_100ml"] == 0
)
chi2, p, dof, expected = chi2_contingency(conforms)

rate_low = conforms.loc[False, True] / conforms.loc[False].sum()
rate_high = conforms.loc[True, True] / conforms.loc[True].sum()

print(f"conforming, chlorine below 0.2 : {rate_low:.1%}")
print(f"conforming, chlorine at/above  : {rate_high:.1%}")
print(f"difference                     : {rate_high - rate_low:+.1%}")
print(f"chi-square p                   : {p:.2e}")
```

Two cautions on that p-value. It says the association is unlikely to be chance;
it does not say chlorination *caused* it, because households that chlorinate
differ from those that do not in ways this survey does not record. And with a
few hundred observations a p-value this small mostly reflects how large the
difference is — the effect size is the number to report, not the p.

## Where the untested households are

```{python}
by_district = wash.groupby("district").agg(
    households=("household_id", "size"),
    chlorine_tested=("free_residual_chlorine_mgl", lambda s: s.notna().mean()),
    ecoli_tested=("ecoli_cfu_100ml", lambda s: s.notna().mean()),
).round(3)
by_district
```

If testing coverage differs by district, a district comparison of water quality
is partly a comparison of which district got tested. Check this before ranking
anything.

## What to report

The finding, the denominator it rests on, and the fact that the tested subset
over-represents households that treat their water. Then the recommendation the
programme can actually act on — which here is to test more households, not to
conclude anything firm about chlorination from 320 of them.