The pipeline: from incomplete counts to published numbers

Surveillance counts arrive late. A case from week W can be reported in week W, W+1, W+2 or later. So the newest weeks look lower than they will be. A trend, an intensity level or an alert measured on them has the same dip.

This vignette takes one synthetic series through csalert. It completes the weeks that are still arriving, measures how far to trust that completion, and ends with numbers that you can publish. Three formats carry the run:

csfmt_reporting_triangle_v3, the input. It holds one cell per reference week and reporting DATE, so it records what was reported when. A delay is a number of days.
csfmt_ensemble_v3, the working format. It holds $data, and one matrix of Monte-Carlo $draws per measure, with rows = weeks and columns = draws. Each stage adds draw columns, so the uncertainty travels forward.
a quantile collapse of those draws, the output. It can be healed into cstidy::csfmt_rts_data_v3 for the usual plots and tables.

One rule: every stage takes the ensemble and returns it

Every analysis stage takes a csfmt_ensemble_v3 and returns a csfmt_ensemble_v3. ens_collapse() is the last stage: it reduces the draws to quantiles, and a collapsed table is output, never input.

So the pipeline is one chain through the ensemble:

csfmt_reporting_triangle_v3
  -> nowcast_delay_ecdf_v1  OR  nowcast_passthrough_to_ensemble_v1
  -> [ens_add_rate]              # when a denominator exists
  -> [short_term_trend]          # growth rate + P(increasing)
  -> [mem_thresholds_v1]         # MEM intensity
  -> [signal_detection_hlm]      # per-draw exceedance
  -> ens_collapse(heal = TRUE)   # ensemble -> quantiles -> csfmt_rts_data_v3

The stages in brackets are optional, and their order is free. Each adds draw columns and passes the ensemble on. A collapsed table holds quantiles, not draws, so short_term_trend() on it cannot get back a per-draw trend or a P(increasing). The draws are the uncertainty, and the collapse is where you spend them.

The run below has nine stages:

Nowcast: complete the recent weeks that are still being reported.
Validation: replay the method on what was known in the past.
Reporting completion: read the reporting delay from the triangle.
Rate: a nowcast numerator over a nowcast denominator, per draw.
Short-term trend: the recent slope of that rate, per draw.
MEM intensity: seasonal intensity thresholds, per draw.
Signal detection: historical-limits exceedance, per draw.
Collapse: quantiles, healed into cstidy::csfmt_rts_data_v3.
Quality control: check the input, and compare two runs.

Stages 4 to 8 are the ensemble chain above. Stages 2 and 3 take the triangle, because they ask how the data arrives, not what the completed counts say. A question about how the data arrives reads the triangle, and a question about what the data says reads the ensemble.

For a new indicator, run stage 3 first. It supports the choice of max_delay_days, which the other stages then use.

library(data.table)
#> 
#> Attaching package: 'data.table'
#> The following object is masked from 'package:base':
#> 
#>     %notin%
library(csalert)
#> csalert 2026.10.23
#> https://niphr.github.io/csalert/

A synthetic reporting triangle

The generator below simulates a laboratory indicator from ISO week 2019-35 to 2024-18. That span holds five winter seasons, which stages 6 and 7 need. Reporting speeds up at the start of ISO year 2024, and stage 3 finds that change from the triangle alone.

The reporting axis is a calendar date, so a delay is a number of days. The generator draws a whole-week delay, then a day inside that reporting week. So a report lands from day 0, the Monday of the reference week, to day 34.

It emits a numerator and a denominator: positive tests, and tests taken. It draws the positives as a binomial given the tests, so the numerator is a subset of the denominator and their ratio is a proportion. Two independent Poissons would not give one. Each test has one reporting delay, so a result and its test land in the same triangle cell.

The generator calls set.seed() inside itself. With the seed outside, a second call runs on a different random stream, and that Monte-Carlo noise then reads as signal.

sim_reports <- function(first_ref = "2019-35",
                        last_ref = "2024-18",
                        delay_slow = c(0.45, 0.30, 0.15, 0.07, 0.03),
                        delay_fast = c(0.70, 0.20, 0.06, 0.03, 0.01),
                        seed = 1L) {
  set.seed(seed)                                    # pinned inside, see above
  weeks <- cstime::dates_by_isoyearweek$isoyearweek
  monday <- as.Date(cstime::dates_by_isoyearweek$mon)   # each week's own Monday
  i0 <- match(first_ref, weeks)
  i1 <- match(last_ref, weeks)
  rows <- lapply(i0:i1, function(ri) {
    ref <- weeks[ri]
    s <- cos(2 * pi * (cstime::isoyearweek_to_isoweek_n(ref) - 5) / 52)  # winter peak
    n_tests <- rpois(1, 400 + 120 * s)
    positive <- rbinom(n_tests, 1, stats::plogis(-2.2 + 1.1 * s))        # subset of tests
    p <- if (cstime::isoyearweek_to_isoyear_n(ref) <= 2023) delay_slow else delay_fast
    data.table(
      isoyearweek_reference = ref,
      # a whole-week delay, then a day inside that reporting week: delay day
      # 7 * weeks_late + day_of_week, so 0 to 34
      reporting_date = monday[ri] + 7L * sample(0:4, n_tests, TRUE, p) +
        sample(0:6, n_tests, TRUE),
      positive = positive
    )
  })
  # "today" is the last day of the last reference week: drop what has not been
  # reported yet, so the most recent weeks are still incomplete -- the problem a
  # nowcast solves
  rbindlist(rows)[reporting_date <= monday[i1] + 6L]
}

reports <- sim_reports()
triangle_long <- reports[, .(numerator = sum(positive), denominator = .N),
                         by = .(isoyearweek_reference, reporting_date)]
triangle_long[, `:=`(indicator_tag = "lab_flu", location_code = "nation",
                     age = "total", sex = "total")]
head(triangle_long, 4)
#>    isoyearweek_reference reporting_date numerator denominator indicator_tag
#>                   <char>         <Date>     <int>       <int>        <char>
#> 1:               2019-35     2019-08-30         0          24       lab_flu
#> 2:               2019-35     2019-08-26         1          16       lab_flu
#> 3:               2019-35     2019-08-29         0          21       lab_flu
#> 4:               2019-35     2019-08-27         0          16       lab_flu
#>    location_code    age    sex
#>           <char> <char> <char>
#> 1:        nation  total  total
#> 2:        nation  total  total
#> 3:        nation  total  total
#> 4:        nation  total  total
c(cells = nrow(triangle_long),
  numerator_never_exceeds_denominator =
    all(triangle_long$numerator <= triangle_long$denominator))
#>                               cells numerator_never_exceeds_denominator 
#>                                8079                                   1

Wrap it as a csfmt_reporting_triangle_v3. The as-of boundary comes from the data: it is the newest reporting date, not the system clock. The constructor stops on a reporting column that is not a Date, and on a missing reporting date.

tri <- csfmt_reporting_triangle_v3(
  triangle_long,
  id_cols       = c("indicator_tag", "location_code", "age", "sex"),
  reference_col = "isoyearweek_reference",
  reporting_col = "reporting_date",
  value_col     = "numerator"
)
c(as_of = attr(tri, "as_of"),
  newest_reporting_date = max(triangle_long$reporting_date))
#>                 as_of newest_reporting_date 
#>          "2024-05-05"          "2024-05-05"

max_delay_days is the delay horizon of stages 1 to 3, in DAYS. Stage 3 shows how to choose it. For this series it is 35: delay days 0 to 34, the five weeks that start on the Monday of the reference week.

max_delay_days <- 35L

1. Nowcast

Estimand. For each reference week, the count reported within the horizon, at delay days 0 to max_delay_days - 1. A report at delay max_delay_days or later is past the horizon. reporting_triangle_matrix() keeps it apart, in late. The observed count, original, still includes it. The nowcast does not estimate it, and the backtest does not score it. For a settled week the estimand is observed. For the most recent weeks it is not observed, and the nowcast is a predictive distribution of it. Stage 3 shows how much arrives on the last delay day or later, which is how you check that the horizon is wide enough.

nowcast_delay_ecdf_v1() completes each incomplete week from a pool of settled reference weeks. Take a reference week of age d days, so delay days 0 to d are observed. For each pool week s, T_s is its settled total and O_s is its count by delay day d. The draws scale the count of the week by the quantiles of the completion ratio of the pool:

draw[ref] = observed_so_far[ref] * quantile(T_s / O_s)

The quantiles sit at n_sim evenly spaced probabilities from 0 to 1, in random order. So the 5 to 95 band is the observed count times the 5% and 95% pool ratios, to within the spacing of those probabilities. Nothing parametric is added, and no regression is fitted. A regression on 35 daily delay columns would have more predictors than settled training rows.

The engine also pools a daily delay ECDF, p(d), the share of the counts of a week that has arrived by delay day d. p(d) cancels out of every draw. It only decides whether delay day d is completed at all. So the draws are not built from observed_so_far / p(d), the closed-form maximum likelihood estimate under n[ref, d] ~ Poisson(lambda[ref] * p[d]).

delay_window, default 26 WEEKS and the one quantity here that is not days, limits the pool to the settled weeks of about that span. The completion ratios can then follow a reporting pattern that changes, as this one does.

denominator_col nowcasts a second measure on the same draw axis. Stage 4 needs it: a rate is coherent only if both its parts are completed in the same Monte-Carlo draw.

set.seed(2)
ens <- nowcast_delay_ecdf_v1(tri, max_delay_days = max_delay_days, n_sim = 500,
                             denominator_col = "denominator")
ens
#> <csfmt_ensemble_v3> 245 rows | 1 series | draws: numerator_nowcasted, denominator_nowcasted

Collapse the draws to quantiles. original is the count reported so far. A settled week sits exactly on it, and a recent week is completed above it.

q <- ens_collapse(ens, probs = c(0.05, 0.5, 0.95))
tail(q[, .(isoyearweek, original,
           lo  = numerator_nowcasted_q05x0,
           med = numerator_nowcasted_q50x0,
           hi  = numerator_nowcasted_q95x0)], 8)
#>    isoyearweek original       lo      med       hi
#>         <char>    <num>    <num>    <num>    <num>
#> 1:     2024-11       87 87.00000 87.00000 87.00000
#> 2:     2024-12       94 94.00000 94.00000 94.00000
#> 3:     2024-13       80 80.00000 80.00000 80.00000
#> 4:     2024-14       66 66.00000 66.00000 66.00000
#> 5:     2024-15       53 53.00000 53.79078 58.10309
#> 6:     2024-16       52 52.61077 54.53215 63.12222
#> 7:     2024-17       44 47.18014 51.90651 62.87138
#> 8:     2024-18       29 38.81167 46.83598 83.81000

The settled weeks have no band. Only the incomplete weeks at the right-hand edge carry one:

show <- tail(q, 30)
xs   <- seq_len(nrow(show))
plot(xs, show$numerator_nowcasted_q50x0, type = "n",
     ylim = range(show$numerator_nowcasted_q05x0, show$numerator_nowcasted_q95x0),
     xlab = "week (last 30)", ylab = "weekly total (nowcast)")
polygon(c(xs, rev(xs)),
        c(show$numerator_nowcasted_q05x0, rev(show$numerator_nowcasted_q95x0)),
        col = grDevices::adjustcolor("steelblue", 0.25), border = NA)
lines(xs, show$numerator_nowcasted_q50x0, lwd = 2, col = "steelblue4")

Keep the reference weeks and their Mondays. Every delay and every age below is a number of days between two dates, so the later stages need the Mondays, not the week labels.

weeks      <- cstime::dates_by_isoyearweek$isoyearweek
ref_weeks  <- q$isoyearweek
ref_monday <- as.Date(cstime::dates_by_isoyearweek$mon[match(ref_weeks, weeks)])
age        <- as.integer(attr(tri, "as_of") - ref_monday)   # age in DAYS
c(weeks = length(ref_weeks), oldest_age_days = max(age), newest_age_days = min(age))
#>           weeks oldest_age_days newest_age_days 
#>             245            1714               6

2. Validation: replay the method on the past

This stage and the next take the triangle, because they ask how the data arrives. That is a property of the reporting, and a collapse of the delay axis into weekly totals removes it.

The triangle records the day on which each count arrived. So you can rebuild what was known on any past date and replay the engine on it, with no second dated extract.

That rebuild is exact only for an append-only reporting system, in which a count keeps its reporting date for ever. Real systems also correct counts, delete records and reclassify cases, and none of those leave a trace in the current triangle. A case reclassified last month looks as though it always had its new class. Where such revisions matter, the replay understates how much the published numbers moved. Only a dated archive of extracts can measure that.

nowcast_censor() goes back in time. It takes a Date, and returns the triangle without the cells reported after that date, with its as-of boundary moved back:

week_end <- ref_monday + 6L                       # the Sunday that ends each week
past <- nowcast_censor(tri, as_of = week_end[length(week_end) - 8L])
# format() first: c() on a Date and an integer coerces the integer to a Date
c(now  = format(attr(tri, "as_of")), then     = format(attr(past, "as_of")),
  rows_now = format(nrow(tri)),      rows_then = format(nrow(past)))
#>          now         then     rows_now    rows_then 
#> "2024-05-05" "2024-03-10"       "8079"       "7836"

nowcast_truth() gives the target: a data.table with reference and truth, the total of each reference week over delay days 0 to max_delay_days - 1. It keeps only the weeks old enough for that total to be settled, so the newest weeks have no truth yet.

truth <- nowcast_truth(tri, max_delay_days = max_delay_days)
tail(truth, 3)
#>    reference truth
#>       <char> <num>
#> 1:   2024-12    94
#> 2:   2024-13    80
#> 3:   2024-14    66
c(reference_weeks = length(ref_weeks), settled = nrow(truth))
#> reference_weeks         settled 
#>             245             241

A method is any function that takes a triangle and returns a csfmt_ensemble_v3, with its own arguments fixed. That one-argument contract lets the same code replay and compare engines with different arguments.

as_of_weeks MUST be a Date vector, and any other class is an error. Replay as of the Sunday that ends each week, which is what “as of week W” meant when delays were counted in weeks.

method_ecdf <- function(x) {
  nowcast_delay_ecdf_v1(x, max_delay_days = max_delay_days, n_sim = 500)
}
as_of_dates <- tail(week_end, 30)

nowcast_backtest() runs the replay and returns the quantiles: one row per reference, as_of, horizon and quantile_level, with the predicted value. horizon is the number of whole weeks from the reference week to the as-of date. Horizon 0 is the current week, which is the least observed.

bt <- nowcast_backtest(
  tri, method_ecdf,
  max_delay_days = max_delay_days,
  as_of_weeks    = as_of_dates,
  horizons       = 0:3,
  probs          = c(0.05, 0.25, 0.5, 0.75, 0.95),
  seed           = 1
)
head(bt, 4)
#>    reference      as_of horizon quantile_level predicted
#>       <char>     <Date>   <int>          <num>     <num>
#> 1:   2023-38 2023-10-15       3           0.05  19.00000
#> 2:   2023-39 2023-10-15       2           0.05  14.00000
#> 3:   2023-40 2023-10-15       1           0.05  17.62526
#> 4:   2023-41 2023-10-15       0           0.05  15.80000
nrow(bt)
#> [1] 600

nowcast_score_v1() scores each forecast of a replay against its settled truth from nowcast_truth(). It returns one row per horizon:

n: the number of scored forecasts in the row.
coverage_50, coverage_90: the measured share of settled truths inside the nominal 50% and 90% central intervals on this replay. These are sample quantities, not a property of the method.
median_signed, median_abs, q05, q95, p_gt_25, p_gt_50: the revision of the published median against the settled truth, as a fraction of the truth. median_signed is the median signed relative revision over the replayed forecasts. It is a median, not a mean, so it is not the bias in the usual expected-error sense. p_gt_25 and p_gt_50 are empirical exceedance proportions: the share of replayed forecasts whose absolute relative revision was above 0.25 and 0.50. They are not probabilities of anything.
wis, wis_log, coverage_95, width_95_rel_median: the weighted interval score and a 95% interval summary. The help page of nowcast_score_v1() defines them.

To compare methods, replay each method on the same reference and as-of weeks with the same seed, and score each replay. The comparison is then paired by forecast. seed makes the run of each method reproducible. It does not give common random numbers across methods, which would need the methods to use their random numbers in the same way. Here they cannot: nowcast_passthrough_to_ensemble_v1() draws no random numbers at all. It does no completion and republishes the counts reported so far, which makes the numbers easy to read:

methods <- list(
  delay_ecdf  = method_ecdf,
  passthrough = function(x) {
    nowcast_passthrough_to_ensemble_v1(x, max_delay_days = max_delay_days)
  }
)
ev <- rbindlist(lapply(names(methods), function(m) {
  bt_m <- nowcast_backtest(
    tri, methods[[m]],
    max_delay_days = max_delay_days,
    as_of_weeks    = as_of_dates,
    horizons       = 0:3,
    seed           = 1
  )
  cbind(method = m, nowcast_score_v1(bt_m, truth))
}))
ev[, .(method, horizon, n, coverage_50, coverage_90, median_signed, median_abs)]
#>         method horizon     n coverage_50 coverage_90 median_signed median_abs
#>         <char>   <int> <int>       <num>       <num>         <num>      <num>
#> 1:  delay_ecdf       3    29       0.828       1.000        0.0089     0.0154
#> 2:  delay_ecdf       2    28       0.536       0.857        0.0328     0.0612
#> 3:  delay_ecdf       1    27       0.444       0.741        0.1135     0.1252
#> 4:  delay_ecdf       0    26       0.346       0.731        0.3889     0.4376
#> 5: passthrough       3    29       0.241       0.241       -0.0147     0.0147
#> 6: passthrough       2    28       0.036       0.036       -0.0630     0.0630
#> 7: passthrough       1    27       0.000       0.000       -0.1515     0.1515
#> 8: passthrough       0    26       0.000       0.000       -0.3543     0.3543

Read that table as a measurement on this sample only. Each row rests on 26 to 29 scored weeks of one synthetic series, far too few to characterise either engine. A coverage estimate from 26 scored weeks has a standard error of about 0.06 at 0.9 and 0.10 at 0.5. That is before any allowance for dependence between overlapping windows.

The comparison does show the shape of the problem. For the passthrough, construction guarantees one thing: each forecast is no greater than the settled truth. A count still arriving cannot exceed its own total, so every single revision is zero or negative.

That alone does not force a negative median at a given horizon. If more than half the weeks were already complete, the median would be exactly zero. Nor does it order the horizons. The negative median_signed at all four horizons, and horizon 0 as the worst, are findings on this sample, not results of the construction.

Its coverage_50 and coverage_90 are equal, because a single draw gives no interval, only a point. The point covers the truth only when the count reported so far already equals it. At horizon 3 that happens on the weeks where nothing arrived on delay days 28 to 34.

The engine completes too much at horizon 0

The delay_ecdf engine covers 0.731 of the settled truths at horizon 0 against a nominal 0.90, and its median_signed is positive, at 0.389. That is not a random dip. The replay window ends at the as-of date, so it spans the reporting speed-up at the start of ISO 2024. delay_window trains on the settled weeks of about the 26 weeks before. So the engine learns the delay curve of the slow regime, and applies it to the counts of the fast one. It scales up counts that were already nearly complete.

A replay on as-of weeks inside the slow regime shows the contrast:

slow_dates <- tail(week_end[cstime::isoyearweek_to_isoyear_n(ref_weeks) == 2023], 30)
c(from = slow_dates[1], to = slow_dates[length(slow_dates)])
#>         from           to 
#> "2023-06-11" "2023-12-31"
bt_slow <- nowcast_backtest(tri, method_ecdf, max_delay_days = max_delay_days,
                            as_of_weeks = slow_dates, horizons = 0:3, seed = 1)
nowcast_score_v1(bt_slow, truth)[
  , .(horizon, n, coverage_50, coverage_90, median_signed, median_abs)]
#>    horizon     n coverage_50 coverage_90 median_signed median_abs
#>      <int> <int>       <num>       <num>         <num>      <num>
#> 1:       3    30       0.367       0.733        0.0141     0.0394
#> 2:       2    30       0.467       0.767       -0.0154     0.0453
#> 3:       1    30       0.500       0.800       -0.0027     0.0451
#> 4:       0    30       0.467       0.800       -0.0072     0.1896

Horizon 0 recovers, and the sign of the revision flips. Read that as a demonstration of the mechanism, not as a measured effect size. The two windows also differ in their weeks and in their part of the season. 26 to 30 scored weeks cannot separate those from the regime change. The general lesson is the reason delay_window exists: a nowcast engine is only as current as the reporting it was trained on. Stage 3 shows when that reporting changed.

3. Reporting completion: how fast does the data arrive?

What `pct_delayD` counts

The delay axis is DAYS. pct_delay0 is the share of the cases of a reference week reported on its own Monday. In general:

pct_delayD is the pooled share of the cases of a reference week reported by the end of the day reference Monday + D.

The index in the name is the delay day, from 0. There are exactly max_delay_days columns, and the last is pct_delay<max_delay_days - 1>. Each is the delay ECDF read at one day, with no interpolation. mean_delay is on the same axis, in DAYS.

Five of the 35 columns fall on the end of an ISO week, and a weekly reader wants those: pct_delay6, pct_delay13, pct_delay20, pct_delay27 and pct_delay34. pct_delay6 is the share in by the end of the reference week itself. The other four are the ends of the four weeks after it.

A triangle with a known answer settles it. Each reference week below reports exactly 50 on delay day 0, 30 on delay day 7 and 20 on delay day 14. The data stops at the newest reporting date, as real data does:

pin_monday <- as.Date("2023-01-02") + 7 * rep(0:29, each = 3)
pin <- data.table(
  isoyearweek_reference = format(pin_monday, "%G-%V"),
  reporting_date = pin_monday + rep(c(0, 7, 14), 30),
  numerator = rep(c(50, 30, 20), 30),
  indicator = "pinned", location = "nation", age = "total", sex = "total"
)
pin <- pin[reporting_date <= as.Date("2023-01-02") + 7 * 29]
tri_pin <- csfmt_reporting_triangle_v3(
  pin, id_cols = c("indicator", "location", "age", "sex")
)
reporting_completion_v1(tri_pin, max_delay_days = 15)[
  , .(n_settled, mean_delay, complete_by_md, pct_delay0, pct_delay7, pct_delay14)]
#>    n_settled mean_delay complete_by_md pct_delay0 pct_delay7 pct_delay14
#>        <int>      <num>          <num>      <num>      <num>       <num>
#> 1:        28        4.9              1         50         80         100

pct_delay0 is 50, not 80: it is the share reported on the Monday of the reference week, delay day 0. pct_delay7 is 80, the cumulative share through delay day 7. mean_delay is 0*0.50 + 7*0.30 + 14*0.20 = 4.9 days.

Coming from an older script? These columns used to be counted in WEEKS. On the old weekly axis, pct_delay1 meant “in by the end of the week after the reference week”, which is now pct_delay13. The names kept their shape and changed their unit. So old code reads a column that exists and gets a different quantity. Check every pct_delayD that you use.

n_settled is 28, not 30. A reference week counts only once as_of - reference Monday >= max_delay_days - 1 days, which excludes the two newest of the 30. A second filter then drops a counted week with a total of zero, and n_settled counts what passes both. No week here is empty, so the age rule alone sets the number. The same holds on the working triangle:

completion <- reporting_completion_v1(tri, max_delay_days = max_delay_days)
c(age_eligible = sum(age >= max_delay_days - 1L), reported = completion$n_settled)
#> age_eligible     reported 
#>          241          241

The two agree only because this series has no week with a zero count. On an indicator with quiet weeks, such as a rare pathogen or a small stratum, n_settled is the smaller. A period with fewer than three weeks left gets no row at all.

A worked reference week

Take reference week 2023-07. cstime gives its dates, and those of the two weeks after it:

cstime::dates_by_isoyearweek[
  isoyearweek %in% c("2023-07", "2023-08", "2023-09"),
  .(isoyearweek, isoyear, mon, thu, sun)
]
#>    isoyearweek isoyear        mon        thu        sun
#>         <char>   <int>     <Date>     <Date>     <Date>
#> 1:     2023-07    2023 2023-02-13 2023-02-16 2023-02-19
#> 2:     2023-08    2023 2023-02-20 2023-02-23 2023-02-26
#> 3:     2023-09    2023 2023-02-27 2023-03-02 2023-03-05

For a case with reference week 2023-07:

statistic	covers reports up to	which is
`pct_delay0`	the Monday of the reference week	Monday 2023-02-13
`pct_delay6`	end of ISO week 2023-07	Sunday 2023-02-19
`pct_delay13`	end of ISO week 2023-08	Sunday 2023-02-26
`pct_delay20`	end of ISO week 2023-09	Sunday 2023-03-05

So pct_delay6 covers the seven days from Monday 2023-02-13. pct_delay13 covers the 14 days from the same Monday, not the seven days of week 2023-08 alone. The columns are cumulative.

The weekday of the extract

The weekday of the extract changes the delay arithmetic. A delay is the number of days from the Monday of the reference week to the reporting date. So the seven days of one reporting week land on seven different delay days:

days <- seq(as.Date("2023-02-13"), as.Date("2023-02-19"), by = "day")
data.table(date = days, weekday = weekdays(days),
           delay_day_for_2023_07 = as.integer(days - as.Date("2023-02-13")))
#>          date   weekday delay_day_for_2023_07
#>        <Date>    <char>                 <int>
#> 1: 2023-02-13    Monday                     0
#> 2: 2023-02-14   Tuesday                     1
#> 3: 2023-02-15 Wednesday                     2
#> 4: 2023-02-16  Thursday                     3
#> 5: 2023-02-17    Friday                     4
#> 6: 2023-02-18  Saturday                     5
#> 7: 2023-02-19    Sunday                     6

On the old weekly axis, all seven of those days shared one label and one delay.

Nothing in the pipeline reads the system clock. The as-of boundary is the newest reporting date in the data: as_of is max(reporting_date). So a Monday extract and a Sunday extract of the same week have different as-of dates. The age in days of every reference week moves with them.

nowcast_censor() shows exactly what a Monday extract would have held. It drops every cell reported after the given date:

monday_as_of <- attr(tri, "as_of") - 6L         # the Monday of the as-of week
tri_mon <- nowcast_censor(tri, as_of = monday_as_of)

rbind(
  cbind(extract = "Sunday (week complete)",
        reporting_completion_v1(tri, max_delay_days = max_delay_days)[
          , .(n_settled, mean_delay, pct_delay6, pct_delay13, pct_delay20, pct_delay27)]),
  cbind(extract = "Monday (six days earlier)",
        reporting_completion_v1(tri_mon, max_delay_days = max_delay_days)[
          , .(n_settled, mean_delay, pct_delay6, pct_delay13, pct_delay20, pct_delay27)])
)
#>                      extract n_settled mean_delay pct_delay6 pct_delay13
#>                       <char>     <int>      <num>      <num>       <num>
#> 1:    Sunday (week complete)       241       9.15       47.9        76.8
#> 2: Monday (six days earlier)       240       9.16       47.8        76.7
#>    pct_delay20 pct_delay27
#>          <num>       <num>
#> 1:        90.6          97
#> 2:        90.6          97

The pooled curve barely moves. n_settled moves, and that is the age test at work. A week counts once as_of - reference Monday >= max_delay_days - 1, so an as_of six days earlier moves the settled boundary back with it.

On this triangle, one reference week is at risk of a total that is too low: the newest settled one. Its age is exactly max_delay_days - 1, so its last delay cell falls on the extract day and is still filling. Every newer week is not settled. Every older week has its last delay day before the extract day. A report at delay max_delay_days or later goes to late, and nowcast_truth() does not count late. So a late report does not move the settled total of an older week, also on real data. Count the weeks that moved:

truth_sun <- nowcast_truth(tri, max_delay_days)
truth_mon <- nowcast_truth(tri_mon, max_delay_days)
cmp <- merge(truth_sun, truth_mon, by = "reference",
             suffixes = c("_sunday", "_monday"))
c(settled_sunday = nrow(truth_sun), settled_monday = nrow(truth_mon),
  in_both = nrow(cmp), totals_that_moved = nrow(cmp[truth_sunday != truth_monday]))
#>    settled_sunday    settled_monday           in_both totals_that_moved 
#>               241               240               240                 0

The reason bounds the whole effect on this triangle. The only cell at risk is the delay-34 numerator of the newest settled week, which holds what arrived on the extract day itself:

newest_settled <- ref_weeks[max(which(age >= max_delay_days - 1L))]
triangle_long[isoyearweek_reference == newest_settled &
              reporting_date == attr(tri, "as_of"),
              .(isoyearweek_reference, reporting_date, numerator, denominator)]
#> Empty data.table (0 rows and 4 cols): isoyearweek_reference,reporting_date,numerator,denominator

Nothing arrived in that cell, so the settled total of this week cannot move. The weekday effect is bounded by the count in the last delay cell. Widen max_delay_days past the day where reporting really stops, and that cell empties. A daily axis makes the cell one day wide instead of one week wide, so there is less count in it to lose.

The effect is small here because of scale, and that does not carry over. The cell at risk is the last delay day of one week, pooled with hundreds of settled weeks. A shorter series, a shorter horizon or a heavier reporting tail makes the same mechanism material. A period with the minimum of three weeks gives that one cell a third of the weight.

In general, a mid-week extract undercounts the newest settled week by the part of its last delay cell that has not arrived yet. nowcast_truth() scores a backtest against that total. A report past the horizon is not in it, so a late report does not change the score. Run the backtest on an end-of-day extract, and record which day it was.

“As of today”, for the weeks on screen

With as_of = 2024-05-05 and max_delay_days = 35, the reference weeks fall into three groups:

data.table(
  isoyearweek = ref_weeks,
  days_observed = age + 1L,
  status = fifelse(age >= max_delay_days - 1L, "settled",
           fifelse(age >= 7L, "still filling", "current week"))
)[, .N, keyby = status]
#> Key: <status>
#>           status     N
#>           <char> <int>
#> 1:  current week     1
#> 2:       settled   241
#> 3: still filling     3

The weeks that are not settled are the weeks that the nowcast of stage 1 completed, and the weeks that the completion table leaves out. The completed weeks are those whose nowcast rises above the count reported so far:

q[numerator_nowcasted_q95x0 > original, isoyearweek]
#> [1] "2024-15" "2024-16" "2024-17" "2024-18"

A change in reporting speed: `period = "year"`

period = "year" returns one row per ISO year that qualifies, and it shows a reporting system that speeds up or slows down. One pooled curve cannot: it averages the regimes and describes neither.

“Qualifies” matters. A period needs at least three settled weeks with a total above zero. A period below that gets no row, and no warning. The year in which your data starts or ends is the usual loss. A missing year means too few weeks, never zero delay.

reporting_completion_v1(tri, max_delay_days = max_delay_days, period = "year")[
  , .(period, n_settled, mean_delay, pct_delay6, pct_delay13, pct_delay20, pct_delay27)]
#>    period n_settled mean_delay pct_delay6 pct_delay13 pct_delay20 pct_delay27
#>    <char>     <int>      <num>      <num>       <num>       <num>       <num>
#> 1:   2019        18       9.38       47.3        73.8        91.0        96.9
#> 2:   2020        53       9.61       44.7        74.8        89.7        96.8
#> 3:   2021        52       9.27       45.9        76.5        90.7        97.1
#> 4:   2022        52       9.67       44.3        74.5        89.6        96.8
#> 5:   2023        52       9.65       44.6        75.0        89.1        96.3
#> 6:   2024        14       6.06       70.7        90.2        96.4        98.9

That is the change built into sim_reports(), found from the triangle. Read pct_delay6, the share in by the end of the reference week itself. It runs 47.3, 44.7, 45.9, 44.3 and 44.6 across ISO 2019 to 2023, then steps to 70.7 in 2024. mean_delay steps the other way, from about 9.5 days to 6.06. A flat run and then a step is what a real regime change looks like. A single year out of line with its neighbours is usually noise.

The pooled pct_delay6 of 47.9 sits just above every one of the slow years. The slow weeks outnumber the fast ones 227 to 14. So the pooled curve is the old regime with a little of the new one in it, not a mean of the two. It will drift year by year as 2024 gains weeks.

The slices are by ISO year, not calendar year. The ISO year of a week is the calendar year of its Thursday. That rule decides the year of a week at the turn of the year:

cstime::dates_by_isoyearweek[
  isoyearweek %in% c("2022-52", "2023-01"), .(isoyearweek, isoyear, mon, thu, sun)
]
#>    isoyearweek isoyear        mon        thu        sun
#>         <char>   <int>     <Date>     <Date>     <Date>
#> 1:     2022-52    2022 2022-12-26 2022-12-29 2023-01-01
#> 2:     2023-01    2023 2023-01-02 2023-01-05 2023-01-08

ISO week 2022-52 has its Thursday on 2022-12-29, so it belongs to ISO year 2022, although its Sunday, 2023-01-01, is in calendar 2023. ISO week 2023-01 has its Thursday on 2023-01-05, and belongs to 2023.

period = "month" slices finer, and shows when a change happened. It uses the same Thursday rule to give a week that spans two months to one of them:

tail(reporting_completion_v1(tri, max_delay_days = max_delay_days, period = "month")[
  , .(period, n_settled, mean_delay, pct_delay6, pct_delay13)], 6)
#>     period n_settled mean_delay pct_delay6 pct_delay13
#>     <char>     <int>      <num>      <num>       <num>
#> 1: 2023-10         4      10.36       39.1        73.9
#> 2: 2023-11         5      10.10       39.4        71.6
#> 3: 2023-12         4       9.77       43.5        75.4
#> 4: 2024-01         4       5.68       72.9        92.4
#> 5: 2024-02         5       6.40       69.3        88.7
#> 6: 2024-03         4       5.90       70.1        90.7

The step falls between 2023-12 and 2024-01. Note the small n_settled per month, four or five weeks, so one month alone is noisy. Read the sequence, not one row.

Every number here depends on `max_delay_days`

reporting_completion_v1() adds late to the last delay column, max_delay_days - 1, before it computes the shares. So its denominator is the total reported by as_of, late reports included. The nowcast estimand leaves late out, so this denominator is larger whenever a late report exists. Two results follow, whatever the real reporting tail looks like:

complete_by_md is the last cumulative share of that total, so it is 1.
The last column, pct_delay<max_delay_days - 1>, is that share as a percentage, so it is 100.

Neither shows the tail. The column before the last one does: 100 - pct_delay<max_delay_days - 2> is the share reported at delay day max_delay_days - 1 or later. mean_delay counts each of those reports at delay day max_delay_days - 1, so it is a lower bound. The sweep below runs one to eight weeks of horizon, in days, and reads the last column by name each time:

sens <- rbindlist(lapply(7L * (1:8), function(md) {
  r <- reporting_completion_v1(tri, max_delay_days = md)
  data.table(max_delay_days = md, n_settled = r$n_settled,
             mean_delay = r$mean_delay, complete_by_md = r$complete_by_md,
             last_col = paste0("pct_delay", md - 1L),
             last_pct = r[[paste0("pct_delay", md - 1L)]],
             pct_delay6 = r$pct_delay6, pct_delay13 = r$pct_delay13)
}), fill = TRUE)
sens
#>    max_delay_days n_settled mean_delay complete_by_md    last_col last_pct
#>             <int>     <int>      <num>          <num>      <char>    <num>
#> 1:              7       245       4.54              1  pct_delay6      100
#> 2:             14       244       7.31              1 pct_delay13      100
#> 3:             21       243       8.52              1 pct_delay20      100
#> 4:             28       242       9.01              1 pct_delay27      100
#> 5:             35       241       9.15              1 pct_delay34      100
#> 6:             42       240       9.16              1 pct_delay41      100
#> 7:             49       239       9.18              1 pct_delay48      100
#> 8:             56       238       9.20              1 pct_delay55      100
#>    pct_delay6 pct_delay13
#>         <num>       <num>
#> 1:      100.0          NA
#> 2:       48.2       100.0
#> 3:       48.2        76.9
#> 4:       48.0        76.8
#> 5:       47.9        76.8
#> 6:       47.8        76.7
#> 7:       47.7        76.6
#> 8:       47.6        76.6
c(complete_by_md_always_1 = all(sens$complete_by_md == 1),
  last_pct_always_100     = all(sens$last_pct == 100))
#> complete_by_md_always_1     last_pct_always_100 
#>                    TRUE                    TRUE

The same holds in every period, where it is easiest to mistake for a finding:

per <- rbindlist(lapply(c("all", "year", "month"), function(p) {
  r <- reporting_completion_v1(tri, max_delay_days = max_delay_days, period = p)
  data.table(period_arg = p, rows = nrow(r),
             complete_by_md_all_1 = all(r$complete_by_md == 1),
             pct_delay34_all_100 = all(r$pct_delay34 == 100))
}))
per
#>    period_arg  rows complete_by_md_all_1 pct_delay34_all_100
#>        <char> <int>               <lgcl>              <lgcl>
#> 1:        all     1                 TRUE                TRUE
#> 2:       year     6                 TRUE                TRUE
#> 3:      month    55                 TRUE                TRUE

The diagnostic that works is the max_delay_days sweep itself. Read the sens table above from top to bottom. mean_delay climbs from 4.54 days at a 7-day horizon to 9.15 at 35 days, then moves only to 9.20 at 56. That flat end is what supports max_delay_days <- 35L here. It is a sensitivity analysis, not a proof. pct_delay6 is 100 at a 7-day horizon, where it is the last column. From 14 days on it moves only from 48.2 to 47.6, because the total at every horizon holds every report.

A mean_delay that kept climbing would be clear evidence that the last column still caps a real tail. A flat end is weaker evidence the other way. A sparse tail and a change in the mix of settled weeks both flatten the curve.

Two cautions on the sweep:

The short horizons are biased. At a 7-day horizon every report from day 6 on counts at day 6. So pct_delay6 is 100 by construction, and mean_delay is too low.
Not all the movement past 35 days comes from the tail. n_settled falls down the rows, because a longer horizon settles fewer weeks. The weeks it drops are the newest, which report fast on this series. That lowers the pooled pct_delay6 a little, on composition alone. Compare rows with equal n_settled, or read period = "year", before you call a small drift a tail.

sim_reports() emits no delay beyond day 34. So on this triangle, a horizon of 35 days gives a late of 0 for every week, and the flat end is exact. We know that only because we can read the generator. On a real series you cannot check it. Widen the horizon until mean_delay stops moving. The share reported on the last delay day or later is then 100 - pct_delay<max_delay_days - 2>. Only a still wider horizon shows how that share spreads over the later days.

4. Rate: a nowcast numerator over a nowcast denominator

From here on, every stage takes the ensemble and returns the ensemble.

Estimand. The percentage of tests that were positive, per reference week. It is computed per draw, so it carries the uncertainty of both the numerator and the denominator. A division of the collapsed medians would lose that, and give a ratio that no draw produced.

Column j is the same Monte-Carlo draw in every measure, so the division is element by element:

ens <- ens_add_rate(ens,
                    numerator   = "numerator_nowcasted",
                    denominator = "denominator_nowcasted",
                    per         = 100)
names(ens$draws)
#> [1] "numerator_nowcasted"                               
#> [2] "denominator_nowcasted"                             
#> [3] "numerator_nowcasted_vs_denominator_nowcasted_pr100"

The naming grammar builds the name of the new measure, so build it the same way and do not type it out:

rate <- csfmt_var("numerator_nowcasted", denom = "denominator_nowcasted", per = 100)
rate
#> [1] "numerator_nowcasted_vs_denominator_nowcasted_pr100"

qr <- ens_collapse(ens, probs = c(0.05, 0.5, 0.95))
tail(qr[, .(isoyearweek,
            positives = numerator_nowcasted_q50x0,
            tests     = denominator_nowcasted_q50x0,
            pct_lo    = round(get(csfmt_var(rate, q = 0.05)), 2),
            pct       = round(get(csfmt_var(rate, q = 0.50)), 2),
            pct_hi    = round(get(csfmt_var(rate, q = 0.95)), 2))], 6)
#>    isoyearweek positives    tests pct_lo   pct pct_hi
#>         <char>     <num>    <num>  <num> <num>  <num>
#> 1:     2024-13  80.00000 482.0000  16.60 16.60  16.60
#> 2:     2024-14  66.00000 440.0000  15.00 15.00  15.00
#> 3:     2024-15  53.79078 433.1894  12.05 12.41  13.28
#> 4:     2024-16  54.53215 445.1033  11.23 12.39  14.23
#> 5:     2024-17  51.90651 392.3792  10.26 12.88  16.19
#> 6:     2024-18  46.83598 427.9788   6.04 10.42  20.19

Two guards apply. A zero denominator gives NA, not a false 0% that would read as a real drop. The numerator is a subset of the denominator, so the rate is capped at per, and a draw above the cap gives a warning.

No warning appeared above, so no draw was above the cap. But the engine nowcasts the numerator and the denominator independently, and nothing keeps them coherent. On a series where the two are close, expect that warning.

5. Short-term trend on the rate

Estimand. The OLS slope of the completed series over the last trend_isoyearweeks weeks, and that slope as a percentage of the current level. It is descriptive: a summary of the six numbers in the window, not an estimate of a latent growth parameter. Per draw, it carries the uncertainty that the nowcast put into those six numbers.

Fit the trend to the nowcast, not to the reported counts. That is the point of the pipeline. The reported counts turn down at the right-hand edge only because the reports have not arrived. A trend fitted to them shows a fall that comes from the reporting.

The method returns a slope, a growth rate and the share of draws with a positive slope. It does not return an increasing or not-increasing label: choose your own cut-off on that share. The cstidy::csfmt_rts_data_v1 method of short_term_trend() does return a status factor, from a different estimator, as the end of this section explains.

The ensemble method fits a closed-form OLS line to the last trend_isoyearweeks nowcast values, separately in every draw column. It adds two draw matrices and one column in $data:

ens_trend <- short_term_trend(
  ens, measure = rate, trend_isoyearweeks = 6,
  denominator = "denominator_nowcasted", per = 100
)
setdiff(names(ens_trend$draws), names(ens$draws))
#> [1] "numerator_nowcasted_vs_denominator_nowcasted_pr100_trend_beta1"
#> [2] "numerator_nowcasted_vs_denominator_nowcasted_pr100_trend_gr"
grep("_trend_", names(ens_trend$data), value = TRUE)
#> [1] "numerator_nowcasted_vs_denominator_nowcasted_pr100_trend_increasing_pr"

..._trend_beta1: the OLS slope, in percentage points per week, per draw.
..._trend_gr: the growth rate 100 * slope / level, per draw, in percent of the current level per week.
..._trend_increasing_pr: the share of draws with a positive slope. It is already a summary over the draws, so it goes in $data, and the collapse leaves it as it is. The method adds it by reference, so ens$data has it too.

Collapse as usual:

gr <- csfmt_var(rate, role = "trend", suffix = "_gr")
qt <- ens_collapse(ens_trend, probs = c(0.05, 0.5, 0.95))
trend <- qt[, .(isoyearweek,
                gr    = round(get(csfmt_var(gr, q = 0.50)), 2),
                gr_lo = round(get(csfmt_var(gr, q = 0.05)), 2),
                gr_hi = round(get(csfmt_var(gr, q = 0.95)), 2),
                p_increasing = get(csfmt_var(rate, role = "trend",
                                             suffix = "_increasing_pr")))]
tail(trend, 10)
#>     isoyearweek     gr  gr_lo gr_hi p_increasing
#>          <char>  <num>  <num> <num>        <num>
#>  1:     2024-09  -1.95  -7.86  4.07        0.262
#>  2:     2024-10  -5.17 -11.05  1.34        0.082
#>  3:     2024-11  -7.37 -13.76 -1.13        0.036
#>  4:     2024-12  -8.59 -14.58 -2.67        0.016
#>  5:     2024-13 -10.82 -18.05 -5.02        0.002
#>  6:     2024-14  -7.47 -13.44 -1.24        0.042
#>  7:     2024-15 -12.46 -20.04 -4.56        0.014
#>  8:     2024-16 -12.61 -21.26 -4.84        0.018
#>  9:     2024-17 -10.13 -20.66 -1.18        0.034
#> 10:     2024-18  -8.83 -32.10  5.18        0.152

The last four rows are the nowcast weeks, and every one of their point estimates is negative. Three of the four bands lie entirely below zero. The newest band is the widest, from -32.10 to 5.18, with p_increasing at 0.152. The series is past its winter peak, on the way into spring.

That result says the completed rate fell over each six-week window, given this nowcast. It does not say that the fall will continue. It depends on the nowcast being right about the four weeks that are still filling. Stage 2 measured those weeks as the weakest point of the engine.

Every week gets an interval, the settled weeks too

Two things are uncertain, and both reach the draws.

The first is the level. A week still filling has a nowcast, so its numerator and denominator differ from draw to draw. A settled week does not: every draw of its observed total is the same number.

The second is the line fitted through those levels. The line is never certain, not even when the six weekly rates are known exactly. Six known points do not fix the slope of the process that made them. So short_term_trend() always adds the sampling error of the slope, se * reference, to each draw before it forms the growth rate.

The slope error cannot be turned off. p_increasing is the share of draws with a positive slope. Without the slope error, a settled week has the same slope in every draw. The share is then exactly 0 or 1 for every settled week, and a threshold on it does nothing.

settled_weeks <- truth$reference
tr_settled <- trend[isoyearweek %in% settled_weeks & !is.na(gr_lo)]
p <- tr_settled$p_increasing
data.table(
  settled_rows_with_a_full_window = nrow(tr_settled),
  any_degenerate_interval = any(tr_settled$gr_lo == tr_settled$gr_hi),
  p_increasing_strictly_inside_0_1 = sum(p > 0 & p < 1),
  p_increasing_at_a_boundary = sum(p %in% c(0, 1))
)
#>    settled_rows_with_a_full_window any_degenerate_interval
#>                              <int>                  <lgcl>
#> 1:                             236                   FALSE
#>    p_increasing_strictly_inside_0_1 p_increasing_at_a_boundary
#>                               <int>                      <int>
#> 1:                              236                          0

Every settled week has a real interval, and each of the 236 has a p_increasing strictly between 0 and 1. A week can still get exactly 0 or 1. Its window then has a slope so many standard errors from zero that every draw falls on the same side. That is a probability of 1 at the Monte-Carlo resolution that the number of draws buys, not a sign. Raise n_sim to tell 0.999 from 1.

What that interval is, and what it is not

The slope error answers a model-based question. Treat the six weekly values as noisy observations around a latent straight line, and report the uncertainty of its slope. A trend arrow asks that question. It commits you to three assumptions:

the weekly mean is linear across the window,
the deviations are independent across weeks,
the deviations are homoskedastic, with the distribution that error_reference names.

None of the three holds for this series. sim_reports() draws tests from a Poisson and positives from a binomial, both with a sinusoidal mean. So the mean in the window is curved, and the variance of the rate changes with the level and with the number of tests. Read the interval as the sensitivity of the slope, not as a validated confidence interval, and say which one you mean when you publish it.

A practical trap comes on top of the assumptions. Under error_reference = "auto", the identity and quasi-Poisson families use a t on trend_isoyearweeks - 2 degrees of freedom. That is 4 at the width of 6 used here. At the default width of 3 it is 1, a Cauchy: the growth-rate quantiles then have no finite variance and very heavy tails. Widen the window, or pass error_reference = "normal". A Wald interval uses the normal, and so does a pipeline built on stats::confint().

Width 2 is an error for those two families. Both estimate a dispersion from width - 2 residual degrees of freedom, so at width 2 there are none and se has no value. The quasi-Poisson dispersion is floored at 1, so its se is never smaller than the Poisson one. The identity residual variance has a floor only when you pass denominator and per, as the call above does. The floor is the Poisson variance of the rate, mean(rate) * per / mean(denominator). For a count, pass a denominator of 1 and per = 1, which gives the mean of the window. Without a denominator there is no floor, and values on an exact line get a p_increasing of exactly 0 or 1. family = "binomial" fixes the dispersion at 1, and accepts width 2.

The ensemble method is not the same estimator as short_term_trend() on a cstidy::csfmt_rts_data_v1. That one fits a quasi-Poisson log-link model over a moving window, and returns a status factor, increasing or notincreasing, after Benedetti (2019); see ?short_term_trend. The ensemble method defaults to an OLS slope, because it is a fixed linear filter that runs down 500 draw columns at once. family = "quasipoisson" gives it the log link of the table method, at the same cost, because it runs the iteration as a filter too. The csfmt_rts_data_v1 table method is for counts without delayed reporting. For counts with delayed reporting, run short_term_trend() on the ensemble, as above, before ens_collapse().

Carry the trend ensemble forward:

ens <- ens_trend

6. MEM intensity thresholds

Estimand. Which of five seasonal intensity levels the rate of this week falls in. The thresholds come from earlier seasons of the same series, by the Moving Epidemic Method. The classification is per draw, so the answer is a distribution over the levels, not one label.

This is the first stage that needs history, and a lot of it. mem_thresholds_v1() needs complete earlier seasons to fit the thresholds of a season. min_seasons (default 2) is a hard floor. prefer_seasons (default 5) is the depth below which a fit is provisional. A season needs min_weeks_per_season (default 30) weeks to count as training at all. That is why this series starts in 2019, not last year.

ens <- mem_thresholds_v1(ens, measure = rate)
#> mem_thresholds_v1: 3 season(s) fit on < 5 training seasons (provisional); see mem_n_seasons.
tail(ens$data[, .(isoyearweek, mem_n_seasons, mem_preepidemic,
                  mem_medium, mem_high, mem_veryhigh)], 3)
#>    isoyearweek mem_n_seasons mem_preepidemic mem_medium mem_high mem_veryhigh
#>         <char>         <int>           <num>      <num>    <num>        <num>
#> 1:     2024-16             4        15.82896   24.07276 26.55251     27.72845
#> 2:     2024-17             4        15.82896   24.07276 26.55251     27.72845
#> 3:     2024-18             4        15.82896   24.07276 26.55251     27.72845

Read the message. It says that some seasons were fit on fewer than prefer_seasons training seasons, and mem_n_seasons records how many each one got. Thresholds fit on two seasons are worth much less than thresholds fit on five, and this is where you see which you have.

After the collapse, the classification is a probability per level:

qm <- ens_collapse(ens, probs = 0.5)
pcols <- grep("_status_prob_", names(qm), value = TRUE)
intensity <- qm[, c("isoyearweek", pcols), with = FALSE]
setnames(intensity, pcols, sub(".*_status_prob_", "", pcols))
tail(intensity, 6)
#>    isoyearweek preepidemic   low medium  high veryhigh
#>         <char>       <num> <num>  <num> <num>    <num>
#> 1:     2024-13       0.000 1.000  0.000 0.000        0
#> 2:     2024-14       1.000 0.000  0.000 0.000        0
#> 3:     2024-15       1.000 0.000  0.000 0.000        0
#> 4:     2024-16       1.000 0.000  0.000 0.000        0
#> 5:     2024-17       0.900 0.100  0.000 0.000        0
#> 6:     2024-18       0.826 0.164  0.008 0.002        0

Some weeks belong to a season with no thresholds: the early seasons, which have no earlier seasons to learn from. Every draw of those weeks is NA, so they get no level:

status <- ens$draws[[csfmt_var(rate, role = "status")]]
c(weeks = nrow(status),
  weeks_with_no_threshold = sum(apply(status, 1, function(r) all(is.na(r)))))
#>                   weeks weeks_with_no_threshold 
#>                     245                     100

A real influenza series would also use exclude_seasons. It keeps unusual seasons, such as a pandemic year or a season with a data gap, out of the training data. An excluded season still gets thresholds; only the data that they are fit on changes.

7. Signal detection: historical limits

Estimand. Whether the rate of this week is above what the same week looked like in earlier years: the historical-limits method. The baseline is the mean and standard deviation of the same week, plus and minus one week, across baseline_isoyears earlier years. The limit is the 99.5th percentile of a normal with that mean and standard deviation. The method compares every draw with the limit, so the output is an exceedance probability, not a yes or no.

The default baseline_isoyears = 5 needs five full years of history before the first week gets a limit. This series is shorter, so it uses three, which is the kind of trade every new indicator faces:

ens <- signal_detection_hlm(ens, measure = rate, baseline_isoyears = 3)
qh <- ens_collapse(ens, probs = 0.5)
hcols <- grep("_hlmstatus_prob_", names(qh), value = TRUE)
signal <- qh[, c("isoyearweek", "hlm_threshold", hcols), with = FALSE]
setnames(signal, hcols, sub(".*_hlmstatus_prob_", "p_", hcols))
tail(signal, 6)
#>    isoyearweek hlm_threshold p_null p_high
#>         <char>         <num>  <num>  <num>
#> 1:     2024-13      23.88968  1.000  0.000
#> 2:     2024-14      19.26117  1.000  0.000
#> 3:     2024-15      18.20797  1.000  0.000
#> 4:     2024-16      16.16278  1.000  0.000
#> 5:     2024-17      16.16115  0.946  0.054
#> 6:     2024-18      15.66598  0.818  0.182

p_high is the share of draws above the limit. It describes the nowcast uncertainty: how likely the completed rate is above the historical limit. It is not a p-value, a posterior probability or a false-alarm rate. The baseline comes from the point history with no uncertainty attached, so hlm_threshold is treated as known.

The first weeks have too little history for a baseline, and every draw of them is NA:

hlm <- ens$draws[[csfmt_var(rate, role = "hlmstatus")]]
c(weeks = nrow(hlm),
  weeks_with_no_baseline = sum(apply(hlm, 1, function(r) all(is.na(r)))))
#>                  weeks weeks_with_no_baseline 
#>                    245                    157

8. Collapse: the end of the chain

ens_collapse() reduces every draw matrix to quantile columns, and every status matrix to a probability per level. With heal = TRUE it passes the result to cstidy::set_csfmt_rts_data_v3(), which adds the standard calendar columns and returns a csfmt_rts_data_v3.

final <- ens_collapse(ens, probs = c(0.05, 0.5, 0.95), heal = TRUE)
class(final)
#> [1] "csfmt_rts_data_v3" "data.table"        "data.frame"
c(rows = nrow(final), cols = ncol(final))
#> rows cols 
#>  245   52
tail(final[, .(isoyearweek, isoyear, isoweek, season, seasonweek)], 3)
#>    isoyearweek isoyear isoweek    season seasonweek
#>         <char>   <int>   <int>    <char>      <num>
#> 1:     2024-16    2024      16 2023/2024         34
#> 2:     2024-17    2024      17 2023/2024         35
#> 3:     2024-18    2024      18 2023/2024         36

This is where the pipeline ends. A csfmt_rts_data_v3 is a format for plots, tables and reports, not for analysis, and there is no way back. The draws are gone, so no per-draw trend, level probability or exceedance probability can come after this point. Every stage that needs draws runs before the collapse, which is why stages 4 to 7 come before it.

csdb ships table validators for csfmt_rts_data_v1 and v2 only, so it cannot validate a healed v3 table.

9. Quality control

The quality-control functions check a run and change nothing. They return a verdict, and the caller decides what to do with it.

qc_surveillance_data_v1() checks the input of one indicator before the nowcast. It checks that the input has rows and a reference column, and that its newest week is not older than expected:

qc_surveillance_data_v1(triangle_long)
#> $ok
#> [1] TRUE
#> 
#> $reasons
#> character(0)
qc_surveillance_data_v1(triangle_long, expect_latest = "2024-19")
#> $ok
#> [1] FALSE
#> 
#> $reasons
#> [1] "latest reference 2024-18 < expected 2024-19 (feed not updated)"

qc_week_over_week_v1() compares this week’s run with last week’s run. Here last week’s run is the same pipeline on the triangle as of seven days earlier, which nowcast_censor() rebuilds. Each run nowcasts the count and adds signal detection, which gives it a status column:

run <- function(x) {
  set.seed(3)
  e <- nowcast_delay_ecdf_v1(x, max_delay_days = max_delay_days, n_sim = 500)
  e <- signal_detection_hlm(e, measure = "numerator_nowcasted", baseline_isoyears = 3)
  ens_collapse(e, probs = 0.5)
}
this_week <- run(tri)
last_week <- run(nowcast_censor(tri, as_of = attr(tri, "as_of") - 7L))

qc <- qc_week_over_week_v1(this_week, last_week, max_delay_weeks = 5)
qc$integrity
#> Empty data.table (0 rows and 5 cols): indicator_tag,isoyearweek,column,prv,cur
qc$signal
#>    indicator_tag isoyearweek                              column  from    to
#>           <char>      <char>                              <char> <num> <num>
#> 1:       lab_flu     2024-18 numerator_nowcasted_hlmstatus_q50x0    NA     1
#>    change
#>    <char>
#> 1:    new

$integrity lists the settled weeks whose published median changed between the two runs. It is empty here, because no report in this series arrives after its week settles. $signal lists the status changes in the newer weeks, and a new week counts as a change.

Both tables come from compare_results(), which holds every value column of every week in both runs. Here are the nowcast medians of the newer weeks, which move from one run to the next:

changes <- compare_results(this_week, last_week)
changes[column == "numerator_nowcasted_q50x0" & isoyearweek >= "2024-14",
        .(isoyearweek, cur = round(cur, 1), prv = round(prv, 1))]
#>    isoyearweek   cur   prv
#>         <char> <num> <num>
#> 1:     2024-14  66.0  67.0
#> 2:     2024-15  53.8  56.5
#> 3:     2024-16  54.5  62.7
#> 4:     2024-17  51.9  53.3
#> 5:     2024-18  46.8    NA

Both functions need identity columns with csfmt schema names, such as indicator_tag and location_code, as this triangle has. ?compare_results explains the error you get otherwise.

The naming grammar

The pipeline builds its column names from named parts. So code can find a column by its parts, and does not need the name written out:

<measure>[_vs_<denom>][_<role>][_<q-coord> | _prob_<level>][_pr<per>][<suffix>]

csfmt_var("numerator", role = "nowcasted", q = 0.5)
#> [1] "numerator_nowcasted_q50x0"
csfmt_parse("numerator_nowcasted_q50x0")
#> $measure
#> [1] "numerator"
#> 
#> $role
#> [1] "nowcasted"
#> 
#> $q
#> [1] 0.5

csfmt_parse() is not the inverse of csfmt_var(). It removes parts from the right, and matches a role against a fixed list. So it cannot tell which of two role words was the role, and it gets the denominator of the rate name wrong:

nm <- csfmt_var("numerator_nowcasted", denom = "denominator_nowcasted", per = 100)
nm
#> [1] "numerator_nowcasted_vs_denominator_nowcasted_pr100"
csfmt_parse(nm)$denom
#> [1] "denominator"

The _nowcasted of the denominator was read as the role, so "denominator" comes back instead of "denominator_nowcasted". The parse is reliable for a name with one role word, such as numerator_nowcasted_q50x0. Check the result when a measure or a denominator itself ends in a role word.

q_label() and q_value() map a probability to its column label and back. They are not an exact inverse pair either. q_label() holds one decimal-percent digit and exactly two integer-percent digits. So a finer probability is rounded, and p = 1 gives a three-digit label that q_value() cannot read.

q_value(q_label(0.0125))   # rounded to the label grid
#> [1] 0.012
q_label(1)
#> [1] "q100x0"
q_value(q_label(1))        # NA: three integer digits do not parse
#> [1] NA

Where next

Run the chain on several series at once. Every stage works over the time_series_id axis, so a triangle with many locations or age groups runs unchanged. The trend masks the windows that would cross two series, and MEM fits each series on its own.
nowcast_passthrough_to_ensemble_v1() replaces the engine at stage 1 for an indicator that SHOULD NOT be nowcast. Every later stage is the same, so an indicator can skip the nowcast and keep the rest of the pipeline.
reporting_completion_trend_v1() puts the year and month slices of stage 3 into one table with a scope column.
The cstidy::csfmt_rts_data_v1 methods of short_term_trend() and signal_detection_hlm() are for counts without delayed reporting. No vignette runs them, and their help pages hold the only examples. vignette("csalert", package = "csalert") explains when to use them.

Richard Aubrey White