--- title: "corncrake(): Correcting Surveillance Counts for Under-Ascertainment" subtitle: "From an observed count to an estimated true total, with uncertainty" author: "Dr Nicolas Smoll, SCPHU, Sunshine Coast Hospital and Health Service" date: "`r Sys.Date()`" output: html_document: toc: true toc_depth: 3 toc_float: true theme: flatly pdf_document: toc: true toc_depth: 3 number_sections: true latex_engine: xelatex vignette: > %\VignetteIndexEntry{corncrake(): Correcting Surveillance Counts for Under-Ascertainment} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r setup, include=FALSE} knitr::opts_chunk$set(collapse = TRUE, comment = "#>", warning = FALSE, message = FALSE) library(mudnester) ``` ## What `corncrake()` does Corncrakes are famously detected far more often by ear than by eye — the overwhelming majority of field records are calls, not sightings. Surveyors have long applied call-based correction factors to convert a count of detections into an estimate of the true, largely-unheard population behind it. `corncrake()` does the same thing to a surveillance count: every notification system under-ascertains true disease burden to some degree, and that degree is rarely constant — it typically varies by age group, by region, and over time as testing behaviour, clinical thresholds, or case definitions change. `corncrake()` takes an observed count, together with a factor that may vary by stratum and by time, and returns an estimate of the true total sitting behind it — with uncertainty bounds wherever they can be derived. It is designed to run directly on `roost()` output (`count_col` defaults to `"n"`, and `time_col` auto-detects the aggregation column from a `roost_tbl`), but works on any tidy data frame with a count column, per the ecosystem's "no function requires another to have run first" design. --- ## Three ways to get a factor `corncrake()` supports three sources for the ascertainment factor, controlled by `method`: 1. **`"user_supplied"`** (default): you supply `factor_table`, a lookup table of factors that can vary by `group_by` stratum and by a `date_start`/`date_end` validity window. This is the right choice when a factor comes from a published estimate, an external evaluation study (e.g. a capture-recapture study), or expert judgement. 2. **`"ratio_estimate"`**: the factor is derived internally as `secondary_count / count` at each stratum/time point, from a second, more-complete data stream. This is the classic surveillance "multiplier method" — for example, dividing all positive laboratory tests by notified cases to estimate under-notification. 3. **`"severity_anchor"`**: the factor is derived by comparing an *observed* severity ratio already in your data (e.g. a case-fatality or case-hospitalisation rate) against a `reference_rate` representing the believed-true rate from a well-ascertained source. This is the case-fatality/infection-fatality-rate anchor inversion described below, and it's the method to reach for early in a novel outbreak, before any seroprevalence survey exists. --- ## Synthetic data ```{r data} set.seed(42) n <- 400 diag <- data.frame( onset_date = as.Date("2024-01-01") + sample(0:59, n, replace = TRUE), age = sample(0:90, n, replace = TRUE), stringsAsFactors = FALSE ) diag <- preening(diag, age_col = "age", scheme = "flucan_sentinel") cases_monthly <- roost( diag, date_col = "onset_date", time_unit = "month", group_cols = "age_group" ) cases_monthly ``` --- ## User-supplied factors Suppose an evaluation study estimated ascertainment separately for two broad age bands, with a slightly higher (and less certain) multiplier for younger ages, tightening over time as testing improved: ```{r factor-table} factors <- data.frame( age_group = rep(c("0-4", "5-15", "16-49", "50-64", "65+"), each = 2), date_start = rep(as.Date(c("2024-01-01", "2024-02-01")), 5), date_end = rep(as.Date(c("2024-01-31", "2024-02-29")), 5), factor = c(3.2, 2.8, 2.6, 2.3, 1.8, 1.6, 1.5, 1.4, 1.3, 1.2), factor_lower = c(2.4, 2.1, 2.0, 1.8, 1.4, 1.3, 1.2, 1.1, 1.1, 1.0), factor_upper = c(4.2, 3.6, 3.3, 2.9, 2.3, 2.0, 1.9, 1.7, 1.6, 1.5), source = "Illustrative multiplier, SCPHU surveillance evaluation 2025" ) knitr::kable(factors) ``` `group_by` and `time_col` tell `corncrake()` how to match rows in `cases_monthly` against rows in `factors`. Here `time_col` is left `NULL` and auto-detects as `"month"` from `cases_monthly`'s `roost_meta`: ```{r apply-corncrake} cases_corrected <- corncrake( cases_monthly, factor_table = factors, group_by = "age_group" ) cases_corrected[, c("age_group", "month", "n", "ascertainment_factor", "corrected_count", "corrected_count_lower", "corrected_count_upper")] ``` `cases_corrected` is still a `roost_tbl` — `corncrake()` preserves the input class, so it drops straight into any downstream code (or a future `bowerbird::roost_plot()`) written against `roost()` output. --- ## Uncertainty: two different questions `ci_method` controls what `corrected_count_lower`/`corrected_count_upper` actually represent, because "uncertainty" can mean two different things here: - **`"table"`** (default): the factor itself is uncertain (as in `factors` above); the observed count is treated as fixed. Bounds come from `factor_lower`/`factor_upper`. - **`"propagate"`**: the factor is treated as fixed; the observed count is treated as a realisation of a Poisson process with its own sampling uncertainty (exact Poisson confidence interval), which is then propagated through the point factor. - **`"none"`**: point estimate only. ```{r ci-methods} corncrake(cases_monthly, factor_table = factors, group_by = "age_group", ci_method = "propagate")[ , c("age_group", "month", "n", "corrected_count", "corrected_count_lower", "corrected_count_upper") ] ``` Note the `factor_table` doesn't need bounds at all for `ci_method = "propagate"` to work — the uncertainty here comes entirely from the count, not the factor. --- ## Deriving a factor instead: the ratio (multiplier) method If a more-complete secondary data stream is available — say, all positive laboratory results, independent of whether a notification was ever made — `corncrake()` can derive the factor directly rather than requiring you to supply one: ```{r ratio-estimate} lab_positive_monthly <- cases_monthly lab_positive_monthly$n <- round(cases_monthly$n * runif(nrow(cases_monthly), 1.3, 2.5)) cases_corrected2 <- corncrake( cases_monthly, method = "ratio_estimate", group_by = "age_group", secondary_data = lab_positive_monthly, secondary_count_col = "n" ) cases_corrected2[, c("age_group", "month", "n", "ascertainment_factor", "corrected_count")] ``` `secondary_data` must share the same `time_col` and `group_by` column names as the primary data. The derived factor is a point estimate only (`ascertainment_factor_lower`/`_upper` are `NA`); use `ci_method = "propagate"` if you want count-based bounds alongside a ratio-estimate factor. --- ## Severity anchor: the CFR/IFR-anchor inversion Both methods above need a second *count* data stream. Early in a novel outbreak — the scenario WHO pandemic-preparedness planning calls "Disease X", where the pathogen is real but its identity, and therefore any tailored surveillance stream, doesn't yet exist — that second count stream usually isn't available yet. What often *is* available is an externally published severity estimate from a reference jurisdiction or a global body, together with your own locally observed severity ratio. This is the case-fatality/infection-fatality-rate anchor inversion described in Smoll et al.'s Queensland COVID-19 under-ascertainment analysis. The logic: if surveillance ascertained every true infection, the observed case-fatality rate (deaths ÷ notified cases) would equal the true infection-fatality rate. Under-ascertainment inflates the observed rate above the true one — by exactly the ascertainment factor: $$\text{UAF} = \frac{\text{CFR}_{\text{obs}}}{\text{IFR}_{\text{ref}}} = \frac{\text{deaths}_{\text{obs}} / \text{cases}_{\text{obs}}}{\text{IFR}_{\text{ref}}}$$ The same identity holds for any other severity outcome — a case-hospitalisation rate against a reference infection-hospitalisation rate works identically. `corncrake()` implements this generally as `method = "severity_anchor"`, comparing `severity_count_col / count_col` in your data against a `reference_rate`. ### A worked example, reproducing the paper's own numbers Suppose a jurisdiction early in a Disease X outbreak observes 100 registered deaths against 5,000 notified cases, and a reference IFR of 1.0% (with a plausible range of 0.5%–2.0%) is available from a high-ascertainment reference jurisdiction: ```{r severity-anchor-worked} disease_x <- data.frame( month = as.Date("2024-01-01"), n_cases = 5000L, n_deaths = 100L ) disease_x_corrected <- corncrake( disease_x, count_col = "n_cases", method = "severity_anchor", time_col = "month", severity_count_col = "n_deaths", reference_rate = 0.01, reference_rate_lower = 0.005, reference_rate_upper = 0.020, reference_source = "WHO Disease X planning scenario, IFR 1.0% (0.5-2.0%)" ) disease_x_corrected[, c("n_cases", "n_deaths", "ascertainment_factor", "ascertainment_factor_lower", "ascertainment_factor_upper", "corrected_count")] ``` The observed CFR here is `100/5000 = 2.0%`, double the 1.0% reference IFR, so `UAF = 2.0` — implying the true infection burden was twice the notified case count, exactly matching the paper's own worked example. ### The bounds invert — and `corncrake()` handles that for you This is the detail worth being deliberate about. Because `UAF` is *divided* by `reference_rate`, it's a **decreasing** function of it: a *higher* reference rate implies *less* under-ascertainment, not more. That means the usual intuition — "lower bound in, lower bound out" — is backwards here: - `ascertainment_factor_lower` is computed from `reference_rate_upper` - `ascertainment_factor_upper` is computed from `reference_rate_lower` In the example above, the 0.5%–2.0% reference range produces a UAF range of `[1.0, 4.0]` — and the *lower* UAF bound (`1.0`) comes from the *higher* reference rate (2.0%), not the lower one. Getting this backwards by hand is an easy mistake to make (the source paper calls it out as a dedicated remark), which is exactly why `corncrake()` encodes it once rather than leaving it as an instruction to re-derive on every use. ### A stratified or time-varying reference rate `reference_rate` doesn't have to be a single scalar. Supply a `factor_table`-shaped data frame with a `rate` column (optionally `rate_lower`/`rate_upper`) for a reference rate that varies by `group_by` stratum or by time window, using exactly the same `date_start`/`date_end` validity-window mechanism as `method = "user_supplied"`'s `factor_table`: ```{r severity-anchor-table} disease_x_stratified <- data.frame( month = as.Date(c("2024-01-01", "2024-01-01")), age_group = c("0-17", "18+"), n_cases = c(1000L, 4000L), n_deaths = c(1L, 99L) ) reference_rates <- data.frame( age_group = c("0-17", "18+"), date_start = as.Date(NA), # open-ended: one reference rate per age group, all time date_end = as.Date(NA), rate = c(0.001, 0.02), source = "Illustrative age-stratified reference IFR" ) corncrake( disease_x_stratified, count_col = "n_cases", method = "severity_anchor", group_by = "age_group", time_col = "month", severity_count_col = "n_deaths", reference_rate = reference_rates )[, c("age_group", "n_cases", "n_deaths", "ascertainment_factor")] ``` `severity_count_col` and `count_col` need to sit in the same table at the same stratification/time — see `vignette("flyway")` for building exactly that shape from a linked cohort's onset and fatality dates in one step. --- ## Rates alongside corrected counts If a population denominator is available, `denominator_col` adds `corrected_rate` (+ bounds) directly: ```{r rates} cases_with_pop <- cases_monthly cases_with_pop$pop <- ifelse(cases_with_pop$age_group == "0-4", 8000, ifelse(cases_with_pop$age_group == "5-15", 15000, ifelse(cases_with_pop$age_group == "16-49", 45000, ifelse(cases_with_pop$age_group == "50-64", 20000, 18000)))) corncrake( cases_with_pop, factor_table = factors, group_by = "age_group", denominator_col = "pop", rate_multiplier = 100000 )[, c("age_group", "month", "corrected_count", "corrected_rate")] ``` --- ## A note on missing factors Not every stratum/time combination in your data needs to be covered by `factor_table` — but if one isn't, `corncrake()` needs to know what to do about it. By default (`on_missing = "warn_na"`) it leaves the row `NA` and issues a single warning naming how many rows were affected; `on_missing = "error"` stops outright, which is useful when you want to be certain your factor table has full coverage before proceeding. ```{r missing-example, error=TRUE} sparse_factors <- factors[factors$age_group != "0-4", ] corncrake(cases_monthly, factor_table = sparse_factors, group_by = "age_group", on_missing = "error") ```