Comparing a facility to its own past — Data Quality Assessment — Data Analysis

Consistency is a comparison, and you get to choose what against

The consistency dimension asks whether the numbers hang together. In practice that means comparing each reported value against something, and there are two candidates.

Against its peers. Facility A reported 40 doses; the district median is 25. Tempting and usually wrong — facilities differ enormously in catchment size, and a large facility will be flagged every month for being large.

Against its own past. Facility A reported 40 doses this month, against its own median of 25 for the year. That is a comparison worth making, because a facility is a reasonable control for itself.

Use the second. Every check in this lesson is a facility against its own history.

Ratio to the median

The simplest useful statistic, and it survives the fact that these series are short.

import pandas as pd

reported = vax[vax["reported"]].copy()

reported["own_median"] = reported.groupby(["facility_id", "antigen"])[
    "doses_administered"
].transform("median")
reported["ratio"] = reported["doses_administered"] / reported["own_median"]

print(reported.nlargest(5, "ratio")[
    ["facility_id", "antigen", "period", "doses_administered", "own_median", "ratio"]
])

library(dplyr)

reported <- vax |>
  filter(report_submitted) |>
  mutate(own_median = median(doses_administered),
         ratio = doses_administered / own_median,
         .by = c(facility_id, antigen))

reported |> slice_max(ratio, n = 5)

Median, not mean. The mean is dragged by the very outlier you are looking for, so a value can inflate its own comparison baseline and hide. This is the same reason SMART plausibility checks use robust statistics.

Where a series is long enough, prefer a rolling median so a genuine seasonal level shift does not flag every month after it:

reported = reported.sort_values(["facility_id", "antigen", "period"])
reported["rolling_median"] = (
    reported.groupby(["facility_id", "antigen"])["doses_administered"]
    .transform(lambda s: s.rolling(6, min_periods=3, center=True).median())
)

reported <- reported |>
  arrange(facility_id, antigen, period) |>
  mutate(rolling_median = zoo::rollapply(doses_administered, 6, median,
                                         partial = TRUE, align = "center"),
         .by = c(facility_id, antigen))

The small-count problem

Run the ratio check on this extract and the top of the list is this:

Facility	Antigen	Month	Value	Own median	Ratio
FAC013	mcv2	December	5	3	1.67
FAC037	opv3	September	13	8	1.62
FAC030	mcv2	October	8	5	1.60

The largest outlier in the entire district is a facility that gave five second doses of measles vaccine instead of three. Two doses. Nobody is going to investigate that, and if your DQA report leads with it, nobody will read the rest.

This is the defining failure of ratio-based flagging, and it is guaranteed rather than unlucky: a ratio has a small denominator at the bottom of the range, so the most extreme ratios always come from the smallest counts. The check is doing exactly what it was asked to do and the question was wrong.

Put an absolute floor beside the ratio

RATIO_HIGH, RATIO_LOW, MIN_ABSOLUTE = 2.0, 0.5, 10

reported["difference"] = reported["doses_administered"] - reported["own_median"]
reported["flag"] = (
    ((reported["ratio"] > RATIO_HIGH) | (reported["ratio"] < RATIO_LOW))
    & (reported["difference"].abs() >= MIN_ABSOLUTE)
)
print(f"{reported['flag'].sum()} flags from {len(reported)} facility-antigen-months")

reported <- reported |>
  mutate(
    difference = doses_administered - own_median,
    flag = (ratio > 2 | ratio < 0.5) & abs(difference) >= 10
  )

Applied here, nothing at all exceeds twice its own median, and eleven facility-antigen-months fall below half of theirs. The largest genuine movement is FAC005’s penta3 in October — 13 doses against a median of 43, a real drop of thirty.

That is the finding, and it is one line long. A check that produces one investigable item from 2,094 observations is working; a check that produces two hundred is a check nobody will run twice.

The direction is not symmetric

A drop and a spike mean different things, and a DQA that treats them the same misses most of what is actually happening.

A drop is usually real or reporting. Stock-out, staff absence, a strike, a facility that started reporting to a different system. Ask about the month, not about the number.
A spike is usually a campaign, an outreach round, or catch-up after a stock-out — all legitimate — or double counting, which is not. Check whether the months either side dip correspondingly, which is the signature of a period boundary problem rather than of extra activity.

reported["previous"] = reported.groupby(["facility_id", "antigen"])[
    "doses_administered"].shift(1)
reported["next"] = reported.groupby(["facility_id", "antigen"])[
    "doses_administered"].shift(-1)

spike_and_dip = reported[
    (reported["ratio"] > 1.5)
    & (reported["previous"] < reported["own_median"] * 0.7)
]

reported <- reported |>
  mutate(previous = lag(doses_administered),
         next_value = lead(doses_administered),
         .by = c(facility_id, antigen))

A high month preceded by a low one is very often one month’s work recorded in the next month’s return. That is a timeliness finding wearing a consistency costume, and the fix is a deadline, not a recount.

Series that are too stable

The opposite check, and the one people forget. Real service delivery is noisy. A facility reporting a nearly identical figure every month is worth a look.

stability = (
    reported.groupby(["facility_id", "antigen"])["doses_administered"]
    .agg(["mean", "std", "size"])
)
stability["cv"] = stability["std"] / stability["mean"]
print(stability[stability["size"] >= 8].nsmallest(5, "cv"))

reported |>
  summarise(mean = mean(doses_administered),
            cv = sd(doses_administered) / mean(doses_administered),
            n = n(),
            .by = c(facility_id, antigen)) |>
  filter(n >= 8) |>
  slice_min(cv, n = 5)

The most stable series in this extract has a coefficient of variation of about 0.06 — six percent month to month. That is low, and it is not evidence of anything. A large facility with a stable catchment and a regular clinic day can easily produce that.

State the limit plainly, because this check is the one most easily misused: a low coefficient of variation puts a facility on the visit list. It does not go in the report as a finding, and it certainly does not go in a sentence containing the word “fabricated”.

Rank facilities, do not flag rows

The output of this lesson is not a list of flagged months. It is a ranking of reporting units, which is what the sampling lesson needed.

risk = (
    reported.groupby("facility_id")
    .agg(flags=("flag", "sum"), months=("flag", "size"))
    .assign(flag_rate=lambda d: d["flags"] / d["months"])
    .sort_values("flag_rate", ascending=False)
)
print(risk.head(6))

risk <- reported |>
  summarise(flags = sum(flag), months = n(), .by = facility_id) |>
  mutate(flag_rate = flags / months) |>
  arrange(desc(flag_rate))

A facility with four flags across twelve months is a candidate for a visit. A single flagged month in a facility that is otherwise steady is a question for a phone call. The unit of action is the facility, because the vehicle goes to a facility.

What comes next

Trend checks compare a number to its own past. The next lesson looks inside the number itself — the last digit, the heaping on multiples of five, and the honest limits of what those can tell you about how a measurement was produced.