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

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

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:

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.


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.
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.


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:

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.


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:

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.

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:

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.


Rates alongside corrected counts

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

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.

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.