Using SPI with POLIS data • spi

This article shows how to prepare a POLIS AFP line list for the simple SPI (spi_simple()) and the model-based SPI (spi_index()). It is for users with POLIS access under a GPEI data-sharing agreement who clean their data with polished. The examples require your own POLIS data and are not run when this article is built. Preparing data for SPI covers data preparation in general.

1. Load the cleaned data

polished downloads and cleans the AFP line list, the population and the district boundaries.

cleaned <- polished::load_polished(country = "AFG")

afp <- cleaned$afp
population <- cleaned$pop_adm2
districts <- readRDS("data/gpei_adm2_shape.rds")

2. What the tables contain

Table	Rows	Columns used here
`afp`	One per AFP case	`adm2_guid`, the paralysis onset date (`paralysis_onset_date`), the onset year (`year_onset`) and the classification (`classification_all`)
`population`	One per district and year	`adm2_guid`, `year` and the under-15 population (`u15_pop`)
`districts`	One per district	`adm2_guid` and the geometry

districts is the GPEI district shapefile used when running polished. All three tables use the POLIS district GUID (adm2_guid), so they link to each other without matching names.

Cases are assigned to a month and a year by their date of paralysis onset.

3. Decide how to handle pending classifications

Cases still pending classification are mostly recent cases. The expected history can be biased if recent years contain many unresolved cases but previous years do not.

Decide whether pending cases count as NPAFP before aggregating, and apply the same rule to every year. Define the rule once and use it for both the simple and the model-based SPI:

# NPAFP with pending cases; use "NPAFP" alone to leave them out
npafp_classes <- c("NPAFP", "PENDING", "LAB PENDING")

afp_npafp <- afp |>
  dplyr::filter(classification_all %in% npafp_classes)

Leaving pending cases out lowers the most recent counts relative to earlier years. Counting them includes some cases that may later receive another classification.

4. NPAFP counts for the simple SPI

Count the cases by district and year of onset:

npafp_year <- afp_npafp |>
  dplyr::count(adm2_guid, year = year_onset, name = "npafp_cases")

spi_results <- spi_simple(
  npafp_year,
  population = population,
  id_col = "adm2_guid",
  pop_col = "u15_pop"
)

A count from the line list has no row for a district-year without a case. Passing the population table means those district-years are included with zero cases.

polished can also produce the counts directly. polished::calc_polio_indicators() returns, at district level, one row per district and year with the district GUID in guid and the number of NPAFP cases in npafp_cases:

indicators <- polished::calc_polio_indicators(
  cases = afp,
  population = population,
  levels = "adm2"
)

npafp_year <- indicators$adm2 |>
  dplyr::select(adm2_guid = guid, year, npafp_cases)

By default calc_polio_indicators() counts pending cases as NPAFP (include_pending = TRUE). Set include_pending to match the rule chosen in section 3. Its district table lists only district-years with AFP activity, so pass the population table to spi_simple() as above.

The simple SPI calculates expected counts from each district’s own earlier NPAFP rate. If a district has no earlier year or no earlier NPAFP case, history_check records the reason and SPI is NA. Include these rows in the surveillance review, using the NPAFP rate and the available history.

5. Monthly counts for the model-based SPI

The model needs one row per district and month, including months with no case. Count the cases by district and month of onset:

months <- seq(as.Date("2016-01-01"), as.Date("2025-12-01"), by = "month")

npafp_month <- afp_npafp |>
  dplyr::mutate(
    month = as.Date(format(paralysis_onset_date, "%Y-%m-01"))
  ) |>
  dplyr::count(adm2_guid, month, name = "count")

Then build the monthly grid. Start from all districts in the boundary file, not only districts that appear in the case data. This keeps districts that reported zero NPAFP cases during the entire period.

cases <- tibble::tibble(adm2_guid = unique(districts$adm2_guid)) |>
  dplyr::cross_join(tibble::tibble(month = months)) |>
  dplyr::left_join(npafp_month, by = c("adm2_guid", "month")) |>
  dplyr::mutate(count = dplyr::coalesce(count, 0L))

6. Check the inputs and calculate the model-based SPI

spi_check_inputs() reports any case or population GUID that is not in the district layer, along with missing months and population problems:

spi_check_inputs(
  cases = cases,
  population = population,
  shapefile = districts,
  id_col = "adm2_guid",
  pop_col = "u15_pop"
)

Then fit the model and calculate the SPI. This needs INLA and takes several minutes; see Calculating the model-based SPI.

spi_year <- spi_index(
  cases = cases,
  population = population,
  adjacency = districts,
  first_assessment = 2019,
  id_col = "adm2_guid",
  pop_col = "u15_pop"
)