3 Simulating biased published estimates
We simulate an outbreak exactly as in vignette("epidist"), a stochastic epidemic with a lognormal delay between primary (symptom onset) and secondary (case notification) events, both interval censored to a day.
as_epidist_linelist_data() converts the dates to times in days from the first primary event.
We keep the simulated times as well, to work out what each study would have measured.
as_epidist_marginal_model() adds the censoring windows and the bounds of each delay.
Figure 3.1 shows the event windows of the simulated cases.
Code
set.seed(101)
meanlog <- 1.8
sdlog <- 0.5
growth_rate <- 0.2
outbreak_start_date <- as.Date("2024-01-01")
outbreak <- simulate_gillespie(r = growth_rate, seed = 101)
obs <- simulate_secondary(
outbreak,
dist = rlnorm, meanlog = meanlog, sdlog = sdlog
)
obs_cens <- simulate_dates(
obs,
outbreak_start_date = outbreak_start_date, keep_times = TRUE
)
linelist <- as_epidist_linelist_data(obs_cens)
#> ℹ No observation time column provided, using 2024-04-30 as the observation date (the maximum of the secondary event upper bound).
delays <- as_epidist_marginal_model(linelist)
#> ! Setting 8341 relative observation times (`relative_obs_time`) greater than 74
#> (2x the maximum delay) to Inf.
#> ℹ This improves model efficiency by reducing the number of unique observation
#> times in the data.
#> ℹ The impact on model accuracy should be negligible because these relative
#> observation times are high enough to cause very limited right truncation.
#> ℹ The original relative observation times are available in
#> `orig_relative_obs_time`.
#> ℹ Raise `obs_time_threshold` to avoid this behaviour.
head(delays)
#> # A tibble: 6 × 22
#> ptime_lwr ptime_upr stime_lwr stime_upr obs_time case ptime delay stime
#> <dbl> <dbl> <dbl> <dbl> <dbl> <int> <dbl> <dbl> <dbl>
#> 1 0 1 5 6 120 1 0.0488 5.59 5.64
#> 2 0 1 9 10 120 2 0.0658 9.76 9.82
#> 3 0 1 8 9 120 3 0.219 7.90 8.12
#> 4 0 1 10 11 120 4 0.250 10.4 10.6
#> 5 0 1 4 5 120 5 0.301 3.98 4.28
#> 6 0 1 9 10 120 6 0.314 9.36 9.68
#> # ℹ 13 more variables: pdate_lwr <date>, pdate_upr <date>, sdate_lwr <date>,
#> # sdate_upr <date>, obs_date <date>, pwindow <dbl>, swindow <dbl>,
#> # relative_obs_time <dbl>, orig_relative_obs_time <dbl>, delay_lwr <dbl>,
#> # delay_upr <dbl>, n <dbl>, delay_min <dbl>
Figure 3.1: The primary and secondary event windows of 200 evenly spaced cases from the simulated outbreak, ordered by primary event time.
Now imagine several studies each analysing this outbreak at a different point in time, reporting integer date differences with no adjustment for right truncation, that is cens_adjusted = 0 and trunc_adjusted = FALSE.
We mimic this by taking cases whose delay is less than a study specific cutoff, since a case with a longer delay would not yet have been observed by that study.
Code
naive_snapshot <- function(cutoff) {
seen <- delays |>
filter(.data$delay_upr <= cutoff) |>
pull(delay_lwr)
return(tibble(cutoff = cutoff, naive_mean = mean(seen), n = length(seen)))
}
true_mean <- exp(meanlog + sdlog^2 / 2)
bias_illustration <- map(c(4, 6, 8, 12, 16, 25), naive_snapshot) |>
list_rbind()Code
Figure 3.2: The naive mean delay from a snapshot taken with a given cutoff (points), compared to the true mean delay (dashed line). Studies analysing the outbreak earlier see a more heavily right truncated, and so more biased, sample of delays.
Real reviews collect studies that differ in how they measured delays and in what they chose to report.
We build ten such studies from the same line list with simulate_study(), which works out which cases a study would have seen and what it would have measured from them.
Each row of the design below is one study.
Two of them report a mean and standard deviation with the covariance between the two, one adjusted for both biases and one adjusted for neither.
Click to expand for the study designs and the code to build the table
Code
set.seed(2)
# A study that stopped at a calendar date needs the growth rate of primary
# events over its collection window. The simulated outbreak grows more slowly
# than its rate parameter once it starts to deplete susceptibles, so we take
# the realised rate rather than the one we simulated with.
primary_counts <- as_tibble(delays) |>
filter(.data$ptime_lwr < 20) |>
count(day = .data$ptime_lwr)
accrual_growth_rate <- coef(
glm(n ~ day, family = poisson, data = primary_counts)
)[["day"]]
study_designs <- tribble(
~study, ~report, ~probs, ~cens_adjusted,
"naive cohort", "moments", NA, 0,
"naive IQR", "quantiles", c(0.25, 0.5, 0.75), 0,
"calendar stop", "quantiles", c(0.2, 0.5, 0.8), 0,
"uniform window", "moments", NA, 2,
"midpoint", "quantiles", c(0.3, 0.6, 0.9), 3,
"adjusted (MVN)", "multivariate", NA, 1,
"naive (MVN)", "multivariate", NA, 0,
"delays over 2d", "moments", NA, 0,
"mean and se", "mean_se", NA, 0,
"midpoint window", "moments", NA, 4
) |>
mutate(
trunc_adjusted = c(
FALSE, FALSE, FALSE, FALSE, FALSE, TRUE, FALSE, FALSE, FALSE, FALSE
),
relative_obs_time = c(12, 16, 20, 25, 30, Inf, 14, 18, 22, 26),
trunc_design = c(
"cohort", "cohort", "accrual", "cohort", "cohort", "cohort",
"cohort", "cohort", "cohort", "cohort"
),
delay_min = c(0, 0, 0, 0, 0, 0, 0, 2, 0, 0),
growth_rate = c(0, 0, accrual_growth_rate, 0, 0, 0, 0, 0, 0, 0),
n = c(180, 55, 240, 95, 40, 300, 150, 130, 25, 110)
)
studies_table <- study_designs |>
transmute(
Study = study,
Reports = case_when(
report == "moments" ~ "mean, sd",
report == "mean_se" ~ "mean with a standard error",
report == "multivariate" ~ "mean, sd with their covariance",
.default = paste0("quantiles at p = ", map_chr(probs, toString))
),
Censoring = case_when(
cens_adjusted == 0 ~ "0: date differences",
cens_adjusted == 1 ~ "1: fully adjusted",
cens_adjusted == 2 ~ "2: uniform window",
cens_adjusted == 3 ~ "3: midpoint",
.default = "4: midpoint primary, uniform secondary"
),
Truncation = case_when(
trunc_adjusted ~ "adjusted",
trunc_design == "accrual" ~ "calendar stop",
.default = "cohort cutoff"
),
`Obs time` = relative_obs_time,
`Min delay` = delay_min,
N = n
)Code
| Study | Reports | Censoring | Truncation | Obs time | Min delay | N |
|---|---|---|---|---|---|---|
| naive cohort | mean, sd | 0: date differences | cohort cutoff | 12 | 0 | 180 |
| naive IQR | quantiles at p = 0.25, 0.5, 0.75 | 0: date differences | cohort cutoff | 16 | 0 | 55 |
| calendar stop | quantiles at p = 0.2, 0.5, 0.8 | 0: date differences | calendar stop | 20 | 0 | 240 |
| uniform window | mean, sd | 2: uniform window | cohort cutoff | 25 | 0 | 95 |
| midpoint | quantiles at p = 0.3, 0.6, 0.9 | 3: midpoint | cohort cutoff | 30 | 0 | 40 |
| adjusted (MVN) | mean, sd with their covariance | 1: fully adjusted | adjusted | Inf | 0 | 300 |
| naive (MVN) | mean, sd with their covariance | 0: date differences | cohort cutoff | 14 | 0 | 150 |
| delays over 2d | mean, sd | 0: date differences | cohort cutoff | 18 | 2 | 130 |
| mean and se | mean with a standard error | 0: date differences | cohort cutoff | 22 | 0 | 25 |
| midpoint window | mean, sd | 4: midpoint primary, uniform secondary | cohort cutoff | 26 | 0 | 110 |
simulate_study() takes the line list and the design of one study, so purrr::pmap() builds them one at a time from the rows of the design.
as_epidist_estimates_data() takes the list and combines them.
The same pattern works for a review, with one function that converts a row of the review to the summaries and metadata of a study.
Code
biased_estimates <- study_designs |>
pmap(simulate_study, data = linelist, max_delay = 60) |>
as_epidist_estimates_data()
#> ! "naive IQR" (row 3), "calendar stop" (row 6), and "midpoint" (row 11) report
#> quantiles less than ten days above their smallest counted delay. Quantiles
#> rounded to whole days that close to the origin carry a bias that a larger
#> sample does not shrink.
#> ℹ See the Checks section of `?as_epidist_estimates_data`.
biased_estimates
#> # A tibble: 22 × 17
#> study type value se n p pwindow swindow relative_obs_time
#> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 naive cohort mean 6.27 NA 180 NA 1 1 12
#> 2 naive cohort sd 2.50 NA 180 NA 1 1 12
#> 3 naive IQR quan… 4 NA 55 0.25 1 1 16
#> 4 naive IQR quan… 6 NA 55 0.5 1 1 16
#> 5 naive IQR quan… 9 NA 55 0.75 1 1 16
#> 6 calendar stop quan… 3 NA 240 0.2 1 1 20
#> 7 calendar stop quan… 5 NA 240 0.5 1 1 20
#> 8 calendar stop quan… 7 NA 240 0.8 1 1 20
#> 9 uniform wind… mean 6.99 NA 95 NA 1 1 25
#> 10 uniform wind… sd 3.79 NA 95 NA 1 1 25
#> # ℹ 12 more rows
#> # ℹ 8 more variables: trunc_adjusted <lgl>, trunc_design <chr>,
#> # cens_adjusted <int>, delay_min <dbl>, growth_rate <dbl>,
#> # growth_rate_sd <dbl>, max_delay <dbl>, mvn_id <chr>max_delay sets how far the grid behind the implied summaries extends for a study that adjusted for right truncation, since nothing else bounds its delays.
The default suits most delays, so raise it only for a delay with a long tail.
We set it here, and simulate_study() passes it on to as_epidist_estimates_data(), where ?as_epidist_estimates_data.data.frame gives the default.
