Stratified and partially-pooled CFR
Why stratify, and why pool
The CFR often varies by group – age, treatment centre, calendar period – and a per-group estimate is frequently required. Two limiting cases bracket the options. Fitting each group independently (no pooling) gives high-variance estimates where a group’s sample is small. Fitting one CFR for all groups (complete pooling) ignores between-group variation.
Partial pooling lies between the two: each group’s CFR is
drawn from a common distribution, so groups with more data are estimated
close to independently while groups with less data are shrunk toward the
shared mean. cfrnow is an epidist model, so
prob takes a brms formula and pooling is a
one-line change: prob ~ (1 | group).
A six-site line list
We simulate six treatment centres: two large (n = 450, 250) and four small (n = 15–20), with true CFRs spread around a common mean. The large sites pin the shared mean; the small ones carry little information on their own, so shrinkage acts mainly on them. Site labels are synthetic.
library(cfrnow)
set.seed(20260714)
sites <- data.frame(
site = c("A", "B", "C", "D", "E", "F"),
n = c(450, 250, 20, 15, 15, 18),
cfr = c(0.45, 0.50, 0.30, 0.65, 0.35, 0.60)
)
death_delay <- LogNormal(mean = 11, sd = 6.5)
recovery_delay <- LogNormal(mean = 16, sd = 8)
ll <- do.call(rbind, Map(function(s, n, cfr) {
cbind(
simulate_linelist(
n = n, cfr = cfr, delay = death_delay,
recovery = recovery_delay, onset_days = 45
),
site = s
)
}, sites$site, sites$n, sites$cfr))
d <- prepare_cfr_data(ll, obs_time = as.Date("2026-02-20"), covariates = "site")
c(cases = d$n_cases, deaths = d$n_deaths)
#> cases deaths
#> 768 331Three fits: complete, no, and partial pooling
The only thing that changes across the three fits is the CFR formula; the delay and recovery delay are co-estimated in each. (Two short chains keep the build time short; use more for inference.)
otd <- LogNormal(meanlog = Normal(2.4, 0.2), sdlog = Normal(0.5, 0.15))
otr <- LogNormal(meanlog = Normal(2.7, 0.2), sdlog = Normal(0.5, 0.15))
args <- list(
delay = otd, recovery_delay = otr, prob_prior = Beta(1, 1),
backend = "cmdstanr", chains = 2, iter = 1000, refresh = 0, seed = 1
)
# complete pooling: a single CFR
fit_pool <- do.call(fit_cfr, c(list(d, formula = brms::bf(mu ~ 1, prob ~ 1)), args))
# no pooling: an independent CFR per site (one estimated logit-CFR per site)
fit_sep <- do.call(fit_cfr, c(list(d, formula = brms::bf(mu ~ 1, prob ~ 0 + site)), args))
# partial pooling: per-site CFR shrunk to a common mean
fit_partial <- do.call(fit_cfr, c(list(d, formula = brms::bf(mu ~ 1, prob ~ (1 | site))), args))Reading the estimates
For the complete-pool fit, summary() reports a single
prob row:
summary(fit_pool)
#> quantity mean q2.5 q50 q97.5 rhat ess_bulk
#> 1 prob 0.5068873 0.4670587 0.5068254 0.5471019 1.0024634 1057.3273
#> 2 delay_mean 11.8285807 10.8445526 11.7952104 13.0539195 0.9985004 708.5507
#> 3 delay_sd 8.2678154 6.9779992 8.1925073 9.9966112 1.0007103 758.6395
#> 4 recovery_mean 15.3811710 14.4376727 15.3502430 16.4113405 0.9985705 867.9546
#> 5 recovery_sd 7.7321604 6.6970264 7.7229427 8.9205239 1.0006240 755.1232For a grouped fit it reports one prob[<site>] row
per site:
summary(fit_sep)
#> quantity mean q2.5 q50 q97.5 rhat
#> 1 prob[A] 0.4782209 0.42521542 0.4777632 0.5330519 1.0034819
#> 2 prob[B] 0.5647843 0.49379949 0.5646614 0.6334943 1.0002918
#> 3 prob[C] 0.1924828 0.05191085 0.1808920 0.4033201 1.0011974
#> 4 prob[D] 0.5173554 0.26551994 0.5189106 0.7469972 1.0018583
#> 5 prob[E] 0.5926362 0.28755794 0.6033133 0.8542252 1.0018477
#> 6 prob[F] 0.7721063 0.51701374 0.7841490 0.9508909 1.0101150
#> 7 delay_mean 11.8757707 10.91193100 11.8627544 12.9790122 0.9994797
#> 8 delay_sd 8.3073686 6.97887826 8.2663047 9.9396996 1.0028251
#> 9 recovery_mean 15.3532639 14.40193424 15.3430244 16.5476570 1.0017358
#> 10 recovery_sd 7.7271904 6.61247749 7.6712733 9.0576431 1.0019149
#> ess_bulk
#> 1 1786.055
#> 2 1472.503
#> 3 1272.856
#> 4 1797.341
#> 5 1532.157
#> 6 1572.089
#> 7 1245.251
#> 8 1315.848
#> 9 1507.885
#> 10 1477.967Collect those per-site CFR rows from the two grouped fits for the figure below; the complete-pool fit supplies the shared mean.
cfr_rows <- function(fit, approach) {
s <- summary(fit)
r <- s[grepl("^prob\\[", s$quantity), ] # the prob[<site>] rows
data.frame(
approach = approach,
site = sub("^prob\\[(.*)\\]$", "\\1", r$quantity),
median = r$q50, lower = r$q2.5, upper = r$q97.5
)
}
est <- rbind(cfr_rows(fit_sep, "no pool"), cfr_rows(fit_partial, "partial pool"))
est <- merge(est, sites[, c("site", "n", "cfr")], by = "site")
names(est)[names(est) == "cfr"] <- "truth"
s_pool <- summary(fit_pool)
pool_mean <- s_pool$q50[s_pool$quantity == "prob"]Shrinkage
Each site shows its no-pool and partial-pool estimate against the truth (×) and the complete-pool mean (dotted). For the large sites (A, B) the two estimates nearly coincide. For the small sites (C–F) the no-pool interval is wide and the estimate noisy; partial pooling pulls it toward the mean and narrows the interval sharply – lower variance at the cost of some bias.
library(ggplot2)
est$approach <- factor(est$approach, levels = c("no pool", "partial pool"))
est$site <- factor(est$site, levels = rev(sites$site))
ggplot(est, aes(median, site, colour = approach)) +
geom_vline(xintercept = pool_mean, linetype = 3) +
geom_point(aes(x = truth), shape = 4, size = 2.4, colour = "grey30") +
geom_pointrange(aes(xmin = lower, xmax = upper), position = position_dodge(width = 0.55)) +
scale_y_discrete(labels = function(s) paste0(s, " (n=", sites$n[match(s, sites$site)], ")")) +
labs(
x = "CFR: median & 95% CrI (× = truth, dotted = complete-pool mean)",
y = NULL, colour = NULL
) +
theme_minimal()
Varying the delay by group
A pooled onset-to-death delay that appears irregular is often a mixture of several populations. If sites resolve on different timescales, pool the delay by site as well, so each retains its own onset-to-death distribution rather than contributing to a single multimodal one:
