corncrake() doesCorncrakes 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.
corncrake() supports three sources for the ascertainment
factor, controlled by method:
"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."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."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.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
#> # A tibble: 10 × 3
#> age_group month n
#> <ord> <date> <int>
#> 1 0-4 2024-01-01 7
#> 2 0-4 2024-02-01 13
#> 3 5-15 2024-01-01 20
#> 4 5-15 2024-02-01 28
#> 5 16-49 2024-01-01 87
#> 6 16-49 2024-02-01 80
#> 7 50-64 2024-01-01 37
#> 8 50-64 2024-02-01 28
#> 9 65+ 2024-01-01 48
#> 10 65+ 2024-02-01 52
#>
#> -- roost_meta --------------------------------------
#> time_unit : month
#> date_range : 2024-01-01 to 2024-02-29
#> group_cols : age_group
#> hemisphere : southern
#> n_rows_in : 400
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:
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)
| age_group | date_start | date_end | factor | factor_lower | factor_upper | source |
|---|---|---|---|---|---|---|
| 0-4 | 2024-01-01 | 2024-01-31 | 3.2 | 2.4 | 4.2 | Illustrative multiplier, SCPHU surveillance evaluation 2025 |
| 0-4 | 2024-02-01 | 2024-02-29 | 2.8 | 2.1 | 3.6 | Illustrative multiplier, SCPHU surveillance evaluation 2025 |
| 5-15 | 2024-01-01 | 2024-01-31 | 2.6 | 2.0 | 3.3 | Illustrative multiplier, SCPHU surveillance evaluation 2025 |
| 5-15 | 2024-02-01 | 2024-02-29 | 2.3 | 1.8 | 2.9 | Illustrative multiplier, SCPHU surveillance evaluation 2025 |
| 16-49 | 2024-01-01 | 2024-01-31 | 1.8 | 1.4 | 2.3 | Illustrative multiplier, SCPHU surveillance evaluation 2025 |
| 16-49 | 2024-02-01 | 2024-02-29 | 1.6 | 1.3 | 2.0 | Illustrative multiplier, SCPHU surveillance evaluation 2025 |
| 50-64 | 2024-01-01 | 2024-01-31 | 1.5 | 1.2 | 1.9 | Illustrative multiplier, SCPHU surveillance evaluation 2025 |
| 50-64 | 2024-02-01 | 2024-02-29 | 1.4 | 1.1 | 1.7 | Illustrative multiplier, SCPHU surveillance evaluation 2025 |
| 65+ | 2024-01-01 | 2024-01-31 | 1.3 | 1.1 | 1.6 | Illustrative multiplier, SCPHU surveillance evaluation 2025 |
| 65+ | 2024-02-01 | 2024-02-29 | 1.2 | 1.0 | 1.5 | Illustrative multiplier, SCPHU surveillance evaluation 2025 |
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:
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")]
#> # A tibble: 10 × 7
#> age_group month n ascertainment_factor corrected_count
#> <ord> <date> <int> <dbl> <dbl>
#> 1 0-4 2024-01-01 7 3.2 22.4
#> 2 0-4 2024-02-01 13 2.8 36.4
#> 3 5-15 2024-01-01 20 2.6 52
#> 4 5-15 2024-02-01 28 2.3 64.4
#> 5 16-49 2024-01-01 87 1.8 157.
#> 6 16-49 2024-02-01 80 1.6 128
#> 7 50-64 2024-01-01 37 1.5 55.5
#> 8 50-64 2024-02-01 28 1.4 39.2
#> 9 65+ 2024-01-01 48 1.3 62.4
#> 10 65+ 2024-02-01 52 1.2 62.4
#> # ℹ 2 more variables: corrected_count_lower <dbl>, corrected_count_upper <dbl>
#>
#> -- roost_meta --------------------------------------
#> time_unit : month
#> date_range : 2024-01-01 to 2024-02-29
#> group_cols : age_group
#> hemisphere : southern
#> n_rows_in : 400
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.
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.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")
]
#> # A tibble: 10 × 6
#> age_group month n corrected_count corrected_count_lower
#> <ord> <date> <int> <dbl> <dbl>
#> 1 0-4 2024-01-01 7 22.4 9.01
#> 2 0-4 2024-02-01 13 36.4 19.4
#> 3 5-15 2024-01-01 20 52 31.8
#> 4 5-15 2024-02-01 28 64.4 42.8
#> 5 16-49 2024-01-01 87 157. 125.
#> 6 16-49 2024-02-01 80 128 101.
#> 7 50-64 2024-01-01 37 55.5 39.1
#> 8 50-64 2024-02-01 28 39.2 26.0
#> 9 65+ 2024-01-01 48 62.4 46.0
#> 10 65+ 2024-02-01 52 62.4 46.6
#> # ℹ 1 more variable: corrected_count_upper <dbl>
#>
#> -- roost_meta --------------------------------------
#> time_unit : month
#> date_range : 2024-01-01 to 2024-02-29
#> group_cols : age_group
#> hemisphere : southern
#> n_rows_in : 400
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.
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:
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")]
#> # A tibble: 10 × 5
#> age_group month n ascertainment_factor corrected_count
#> <ord> <date> <int> <dbl> <dbl>
#> 1 0-4 2024-01-01 7 1.86 13
#> 2 0-4 2024-02-01 13 1.31 17
#> 3 5-15 2024-01-01 20 1.95 39
#> 4 5-15 2024-02-01 28 2.18 61
#> 5 16-49 2024-01-01 87 1.59 138
#> 6 16-49 2024-02-01 80 2.28 182
#> 7 50-64 2024-01-01 37 1.81 67
#> 8 50-64 2024-02-01 28 1.96 55
#> 9 65+ 2024-01-01 48 1.48 71
#> 10 65+ 2024-02-01 52 1.54 80
#>
#> -- roost_meta --------------------------------------
#> time_unit : month
#> date_range : 2024-01-01 to 2024-02-29
#> group_cols : age_group
#> hemisphere : southern
#> n_rows_in : 400
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.
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.
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:
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")]
#> n_cases n_deaths ascertainment_factor ascertainment_factor_lower
#> 1 5000 100 2 1
#> ascertainment_factor_upper corrected_count
#> 1 4 10000
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.
corncrake() handles that for
youThis 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_upperascertainment_factor_upper is computed from
reference_rate_lowerIn 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.
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:
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")]
#> age_group n_cases n_deaths ascertainment_factor
#> 1 0-17 1000 1 1.0000
#> 2 18+ 4000 99 1.2375
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.
If a population denominator is available,
denominator_col adds corrected_rate (+ bounds)
directly:
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 tibble: 10 × 4
#> age_group month corrected_count corrected_rate
#> <ord> <date> <dbl> <dbl>
#> 1 0-4 2024-01-01 22.4 280
#> 2 0-4 2024-02-01 36.4 455
#> 3 5-15 2024-01-01 52 347.
#> 4 5-15 2024-02-01 64.4 429.
#> 5 16-49 2024-01-01 157. 348
#> 6 16-49 2024-02-01 128 284.
#> 7 50-64 2024-01-01 55.5 278.
#> 8 50-64 2024-02-01 39.2 196
#> 9 65+ 2024-01-01 62.4 347.
#> 10 65+ 2024-02-01 62.4 347.
#>
#> -- roost_meta --------------------------------------
#> time_unit : month
#> date_range : 2024-01-01 to 2024-02-29
#> group_cols : age_group
#> hemisphere : southern
#> n_rows_in : 400
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.
sparse_factors <- factors[factors$age_group != "0-4", ]
corncrake(cases_monthly, factor_table = sparse_factors, group_by = "age_group",
on_missing = "error")
#> Error:
#> ! (*)> mudnester::corncrake() — 2 row(s) had no matching ascertainment factor for their stratum/time and were left NA.