---
title: "Coverage and dropout rate by facility"
subtitle: "Routine vaccination coverage, 2024 · R"
author: "Cassion · data-analysis.cassion.dev"
format:
  html:
    toc: true
    code-fold: false
---

## What this produces

The same dropout calculation and flagging logic as the Python example, written
with tidyverse verbs.

Every dataset on this platform is synthetic. Coverage here describes no real
district.

## Setup

```{r}
#| message: false
library(readr)
library(dplyr)
library(tidyr)

URL <- paste0(
  "https://data-analysis.cassion.dev/datasets/files/",
  "vaccination-coverage-2024.v1.csv"
)

epi <- read_csv(URL, col_types = cols(
  facility_id = col_character(),
  period      = col_date(),
  .default    = col_guess()
)) |>
  mutate(month = format(period, "%Y-%m"))

glimpse(epi)
```

## A non-report is not zero children vaccinated

```{r}
epi |>
  filter(!report_submitted) |>
  summarise(
    rows            = n(),
    rows_with_zero  = sum(doses_administered == 0),
    facility_months = n_distinct(paste(facility_id, month))
  )
```

Every non-reporting row carries a zero. Sum without filtering and you have
asserted that no child in that catchment was vaccinated that month.

## The denominator is annual, the reporting is monthly

`target_population` is the annual target repeated on every row. Monthly coverage
divides by that figure over twelve; dividing by the full annual figure
understates coverage twelvefold.

```{r}
MONTHS_IN_YEAR <- 12

epi |>
  filter(antigen == "penta3") |>
  group_by(month) |>
  summarise(
    doses         = sum(doses_administered),
    annual_target = sum(target_population),
    completeness  = round(mean(report_submitted), 3),
    .groups = "drop"
  ) |>
  mutate(
    wrong_annual_denominator    = round(doses / annual_target, 3),
    correct_monthly_denominator = round(doses / (annual_target / MONTHS_IN_YEAR), 3)
  )
```

## Dropout between penta1 and penta3

Dropout does not depend on a population estimate at all — both terms come from
the same register, which is why it is a better programme signal than coverage.

```{r}
series <- epi |>
  filter(report_submitted, antigen %in% c("penta1", "penta3")) |>
  group_by(facility_id, antigen) |>
  summarise(doses = sum(doses_administered), .groups = "drop") |>
  pivot_wider(names_from = antigen, values_from = doses) |>
  filter(!is.na(penta1), !is.na(penta3)) |>
  mutate(dropout = (penta1 - penta3) / penta1) |>
  arrange(desc(dropout))

cat("median dropout:", round(100 * median(series$dropout), 1), "%\n")
head(series, 10)
```

A facility reporting near-zero dropout deserves as much suspicion as one
reporting a very high figure. Near-zero usually means penta3 was reconstructed
from penta1 rather than counted.

## The over-reporting flag, applied at the right level

```{r}
by_month <- epi |>
  filter(report_submitted, antigen %in% c("penta1", "penta3")) |>
  select(facility_id, month, antigen, doses_administered) |>
  pivot_wider(names_from = antigen, values_from = doses_administered) |>
  filter(!is.na(penta1), !is.na(penta3)) |>
  mutate(impossible = penta3 > penta1)

by_month |>
  group_by(facility_id) |>
  summarise(ever_flagged = any(impossible), .groups = "drop") |>
  summarise(
    flagged    = sum(ever_flagged),
    facilities = n()
  )
```

Thirty-seven of thirty-eight. A check that flags almost the whole district
identifies nothing: the two doses go to different children in different months
and the counts have no reason to move together, so ordinary noise crosses the
line constantly.

```{r}
persistent <- series |> filter(dropout < 0)
persistent
```

```{r}
by_month |>
  filter(facility_id %in% persistent$facility_id) |>
  group_by(facility_id) |>
  summarise(
    months_flagged  = sum(impossible),
    months_reported = n(),
    .groups = "drop"
  ) |>
  arrange(desc(months_flagged))
```

Six facilities on the annual total. That is a list a supervisor can act on.

## Coverage by antigen

```{r}
epi |>
  filter(report_submitted) |>
  group_by(antigen) |>
  summarise(
    coverage = round(
      sum(doses_administered) / (sum(target_population) / MONTHS_IN_YEAR), 3
    ),
    .groups = "drop"
  ) |>
  arrange(desc(coverage))
```

The denominator counts only facilities that reported, keeping numerator and
denominator on the same set of facilities. That is the subject of the
completeness-adjustment example.

## What to report

Coverage with its denominator stated, dropout as the signal that does not depend
on a population estimate, and the over-reporting list with the aggregation level
it was computed at.