The hierarchy that changes underneath you — Routine Data and DHIS2 — Data Analysis

One tree, and everything hangs off it

Every value in the system is attached to an org unit, and org units form a single tree: country, region, district, facility. The level a unit sits at is a property of the tree, not a column on the unit.

Level 1  Country
Level 2  Region
Level 3  District
Level 4  Facility        <- our 38 units

Aggregation is the tree walked upward. Ask for a district figure and the system sums every facility below it; ask for a national one and it sums every district. This is the same “aggregate to the grain” discipline the joining course taught, with the grain supplied by a hierarchy somebody else maintains.

import pandas as pd

vax = pd.read_csv("vaccination-coverage-2024.v1.csv", parse_dates=["period"])
print(f"{vax['facility_id'].nunique()} facilities in the extract")

library(dplyr)
n_distinct(vax$facility_id)

Thirty-eight, all at level 4. The extract carries no parent, which is the usual case and the reason most district totals in this sector are computed by joining to a facility list rather than by the system.

Groups are not levels

The extract has a facility_type column, and it is worth being clear about what that is.

by_type = (
    vax.groupby("facility_type")["facility_id"].nunique().sort_values(ascending=False)
)
print(by_type)

vax |> summarise(facilities = n_distinct(facility_id), .by = facility_type)

Facility type	Facilities
Health post	20
Health centre	16
District hospital	2

That is an org unit group, not a level. The distinction matters in three practical ways:

A unit belongs to exactly one parent and any number of groups. Health post is a group; it does not tell you where the facility is.
Groups can overlap. A facility can be in “health post”, “hard to reach” and “supported by partner X” at once, and aggregating by two overlapping groups double-counts.
Group sets are the safe version. A group set is a set of mutually exclusive groups — facility type is usually one — and aggregating by a group set is safe in the way aggregating by an arbitrary group is not.

Ask whether the thing you are grouping by is a group set. If it is not, check that the groups partition before you sum.

The property nobody warns you about

Here is the one that costs people a quarter.

DHIS2 stores the current hierarchy, not a historical one. A data value is attached to a facility. The facility is attached to a parent now. There is no record, in the ordinary data model, of which parent it had in March.

So when a facility is reassigned from District A to District B — a boundary change, a reorganisation, a correction — all of its historical data moves with it. District A’s total for last January, already reported and already in a donor’s spreadsheet, changes the next time anybody asks for it.

Three symptoms, all of which look like data quality problems and none of which is:

A district total that differs from the same query run six months ago, with no data entry in between.
A published annual figure that cannot be reproduced.
Two reports of the same period disagreeing, both correct on the day they were run.

What to do about it

You cannot stop the hierarchy moving. You can stop it moving your numbers silently.

Pin the hierarchy with the extract. Pull the org unit list at the same time as the data and store it beside the values.

org_units = pd.read_csv("outputs/extract/2026-07/org-units.csv")
    # ea/facility id, name, parent_id, level, extracted_on

pinned = vax.merge(org_units[["facility_id", "district_id"]],
                   on="facility_id", how="left", validate="many_to_one")
assert pinned["district_id"].notna().all(), "facility with no district in the pinned tree"

org_units <- readr::read_csv(here::here("outputs", "extract", "2026-07", "org-units.csv"))

pinned <- vax |>
  left_join(select(org_units, facility_id, district_id), by = "facility_id",
            relationship = "many-to-one")

stopifnot(!any(is.na(pinned$district_id)))

The join is many_to_one and the assertion is not decorative. A facility present in the data and absent from the pinned tree means the tree was pulled at a different time from the data, which is exactly the situation this is meant to prevent.

Aggregate from the pinned tree, not from the live one. Then a figure published in March is reproducible in September, because the tree it used is on disk.

Record the extraction date on both. The next lesson but one makes that a property of the extract script rather than a discipline.

Reassignment, reconstructed

Where the reassignment date is known — and somebody always knows it, even if the system does not — the honest handling is to carry it explicitly.

reassignments = pd.DataFrame([
    {"facility_id": "FAC017", "from_district": "D01", "to_district": "D02",
     "effective": "2024-04-01", "reason": "boundary revision"},
])

def district_at(facility, period, log=reassignments):
    moves = log[(log["facility_id"] == facility)
                & (pd.to_datetime(log["effective"]) <= period)]
    return moves["to_district"].iloc[-1] if len(moves) else None

reassignments <- tibble::tribble(
  ~facility_id, ~from_district, ~to_district, ~effective,     ~reason,
  "FAC017",     "D01",          "D02",        as.Date("2024-04-01"), "boundary revision"
)

That is a slowly changing dimension, and treating it as one is the correct answer. It is also more work than most teams will do, so the fallback is honest labelling:

District totals are computed on the org unit hierarchy as at 2026-07-28. One facility was reassigned from D01 to D02 in April 2024; its 2024 data appears under D02 throughout, including for periods before the reassignment.

That sentence is the deliverable. It costs nothing and it turns an irreproducible number into a reproducible one.

Aggregating up without double-counting

Two failure modes when you do the aggregation yourself.

A facility appearing twice in the tree. Impossible in DHIS2 — a unit has one parent — but entirely possible in the CSV of facilities somebody emailed you, where a facility that moved appears under both districts.

duplicated = org_units["facility_id"].duplicated(keep=False)
assert not duplicated.any(), org_units.loc[duplicated, ["facility_id", "district_id"]]

stopifnot(!any(duplicated(org_units$facility_id)))

Summing across levels. A district total plus its facilities is the district counted twice. This happens when an extract mixes levels — ask for level 4 only, and check.

print(org_units["level"].value_counts())

org_units |> count(level)

Every extract should be at exactly one level. If it is not, the first thing your script does is filter to one, and say which.

Coverage is a district-level indicator

A point from the survey course arriving here with a mechanism. Facility catchment populations overlap, and people cross boundaries to be vaccinated, so a facility-level coverage figure is a numerator from one population over a denominator from another.

by_facility = (
    vax[vax["report_submitted"] == True]
    .query("antigen == 'penta3'")
    .groupby("facility_id")
    .agg(doses=("doses_administered", "sum"), target=("target_population", "sum"))
)
by_facility["coverage"] = by_facility["doses"] / by_facility["target"]
print(by_facility["coverage"].describe()[["min", "max"]].round(3))

vax |>
  filter(report_submitted, antigen == "penta3") |>
  summarise(coverage = sum(doses_administered) / sum(target_population),
            .by = facility_id) |>
  summarise(min = min(coverage), max = max(coverage))

The spread across facilities is far wider than any real difference in service delivery, because the boundary crossing is noise at facility level and cancels at district level. Aggregate to the level where the crossing happens inside the unit, and report facility figures as workload rather than coverage.

What comes next

You can now place a value in space. The next unit places it in time — period types, which period a late entry belongs to, and the completeness registration that decides whether the value was ever expected at all.