The v4 nowcasts return a
csfmt_ensemble_v4. The rate, trend, HLM and MEM steps take
one and return it with their results added. A full ensemble holds draws
in memory. A compressed ensemble holds a table with a point value and 9
percentiles for each base measure, and you store that table.
This vignette takes one simulated series from weekly counts to a stored table, and back again.
The series has 320 weeks of positive tests out of all tests. Each
week has 100 draws. Here they are simulated, but in a real run they come
from a nowcast such as nowcast_ecdf_v4().
library(data.table)
library(csalert)
w <- cstime::dates_by_isoyearweek$isoyearweek
iyw <- w[match("2018-01", w) + 0:319]
season_week <- cstime::isoyearweek_to_seasonweek_n(iyw)
share <- 0.05 + 0.15 * exp(-(season_week - 25)^2 / 40)
set.seed(1)
tested <- matrix(rpois(320 * 100, 2000), nrow = 320)
positive <- matrix(rbinom(length(tested), tested, share), nrow = 320)
ens <- csfmt_ensemble_v4(
data.table(location_code = "nation", isoyearweek = iyw),
id_cols = "location_code",
draws = list(positive = positive, tested = tested)
)
ens
#> <csfmt_ensemble_v4> full | 320 rows | 1 series | draws: positive, testedEach draw matrix has one row per week and one column per draw. Column
j is the same draw in every measure.
ens <- ens_add_rate_v4(
ens,
numerator = "positive",
denominator = "tested",
per = 100,
name = "positive_pr100"
)
ens <- short_term_trend_v4(
ens,
measure = "positive",
denominator = "tested",
name = "positive_pr100",
seed = 4
)
ens <- signal_detection_hlm_v4(ens, measure = "positive_pr100")
ens$data[
316:320,
.(isoyearweek, positive_pr100_trend_pinc, positive_pr100_hlm_prob_high)
]
#> isoyearweek positive_pr100_trend_pinc positive_pr100_hlm_prob_high
#> <char> <num> <num>
#> 1: 2024-03 0.9258002 0.00
#> 2: 2024-04 0.9228678 0.00
#> 3: 2024-05 0.9110528 0.00
#> 4: 2024-06 0.8697581 0.08
#> 5: 2024-07 0.7609514 0.19The trend models the counts of positive in a
quasi-Poisson regression, with the log of tested as the
offset. positive_pr100_trend_pinc is P(increasing). The HLM
step compares each draw of the rate with a limit from the 5 earlier
years. positive_pr100_hlm_prob_high is the share of the
draws above that limit.
stored <- ens_compress_v4(ens)
dim(stored)
#> [1] 320 49
max(nchar(names(stored), type = "bytes"))
#> [1] 29The table holds the columns of $data. It also holds a
median and 9 percentiles for each of the 3 base measures,
positive, tested and
positive_pr100. It holds no draws, and each column is a
plain vector. ens_compress_v4() stops with an error when a
name is longer than 63 bytes, the limit of a PostgreSQL name. So you can
write the table to PostgreSQL.
back <- as_ensemble_v4(stored, id_cols = "location_code", n_sim = 100)
back
#> <csfmt_ensemble_v4> compressed | 320 rows | 1 series | draws: none
back <- short_term_trend_v4(
back,
measure = "positive",
denominator = "tested",
name = "positive_pr100",
seed = 4
)
pinc <- data.table(
isoyearweek = iyw,
full = ens$data$positive_pr100_trend_pinc,
rebuilt = back$data$positive_pr100_trend_pinc
)
pinc[316:320]
#> isoyearweek full rebuilt
#> <char> <num> <num>
#> 1: 2024-03 0.9258002 0.9315201
#> 2: 2024-04 0.9228678 0.9098614
#> 3: 2024-05 0.9110528 0.8931547
#> 4: 2024-06 0.8697581 0.8187807
#> 5: 2024-07 0.7609514 0.6879770
max(abs(pinc$full - pinc$rebuilt), na.rm = TRUE)
#> [1] 0.1053404Each step rebuilds the draws of the base measures that it reads, from
the stored percentiles, with seed 4. n_sim = 100 gives the
same number of draws as the full ensemble. The rebuild is an
approximation. On this series, the two P(increasing) values differ by at
most 0.11 over the 318 weeks that have a trend.
Each base measure gets its own u, so the rebuilt draws of
positive and tested are independent. The parts
of a full ensemble from nowcast_ecdf_v4() are independent
too. A rate added again from its parts is their ratio. Here
positive depends on tested, so that ratio is
too wide. Its median interquartile range over the 320 weeks is 1.09
times that of the full rate. A step on the stored rate rebuilds it from
its own percentiles, at 1.01 times. The chunk below adds the rate again
from its parts.
back <- ens_add_rate_v4(
back,
numerator = "positive",
denominator = "tested",
per = 100,
name = "positive_pr100"
)
back <- signal_detection_hlm_v4(back, measure = "positive_pr100")
back
#> <csfmt_ensemble_v4> compressed | 320 rows | 1 series | draws: positive_pr100_trend_beta1, positive_pr100_trend_se, positive_pr100, positive_pr100_hlm