reporting delays bias CFR

delays_cfr.R

``` r
# packages ----------------------------------------------------------------

library(tidyverse)
library(incidence2) 

# input -------------------------------------------------------------------

onset <- tibble(
  date = seq(ymd(20200101),ymd(20200104),1),
  cases = c(1,2,3,1),
  group = "onset"
)

report <- tibble(
  date = seq(ymd(20200103),ymd(20200106),1),
  cases = c(1,2,3,1),
  group = "report"
)

death <- tibble(
  date = seq(ymd(20200105),ymd(20200108),1),
  cases = c(1,2,3,1),
  group = "death"
)

# cfr ---------------------------------------------------------------------

bind_rows(
  onset,
  report,
  death
) %>% 
  pivot_wider(
    names_from = group,
    values_from = cases
  ) %>% 
  replace_na(
    list(
      report = 0,
      death = 0)) %>%  
  mutate(
    cum_report = cumsum(report),
    cum_death = cumsum(death)
  ) %>% 
  mutate(cum_cfr = cum_death/cum_report) %>% 
  select(date,starts_with("cum_"))
#> # A tibble: 8 × 4
#>   date       cum_report cum_death cum_cfr
#>   <date>          <dbl>     <dbl>   <dbl>
#> 1 2020-01-01          0         0 NaN    
#> 2 2020-01-02          0         0 NaN    
#> 3 2020-01-03          1         0   0    
#> 4 2020-01-04          3         0   0    
#> 5 2020-01-05          6         1   0.167
#> 6 2020-01-06          7         3   0.429
#> 7 2020-01-07          7         6   0.857
#> 8 2020-01-08          7         7   1

# confidence intervals ---------------------------------------------------

cumcfr <- bind_rows(
  onset,
  report,
  death
) %>% 
  pivot_wider(
    names_from = group,
    values_from = cases
  ) %>% 
  replace_na(
    list(
      Report = 0,
      Death = 0)) %>%  
  mutate(
    cum_report = cumsum(Report),
    cum_death = cumsum(Death)
  ) %>% 
  mutate(cum_cfr = cum_death/cum_report) %>% 
  select(date,starts_with("cum_"))


cumcfr_ic <- cumcfr %>% 
  filter(cum_death>0) %>% 
  mutate(raw=pmap(.l = select(., x=cum_death,n=cum_report),
                  .f = binom.test)) %>%
  mutate(raw=map(.x = raw,.f = broom::tidy)) %>% 
  unnest(raw) %>% 
  select(date,starts_with("cum_"),starts_with("conf"))

cumcfr %>% 
  left_join(cumcfr_ic) %>% 
  ggplot(aes(date,cum_cfr)) +
  geom_point() +
  geom_line() +
  geom_hline(aes(yintercept = 1), lty = 2) +
  geom_ribbon(
    aes(ymin = conf.low, ymax = conf.high),
    alpha = 0.2,
    fill = "blue"
  ) +
  theme_bw() +
  labs(x = NULL, 
       y = NULL,
       title = "Naive CFR")

# ggsave(
#   "figure/delays_cfr-trend.png",
#   width = 4,
#   height = 1.7,
#   dpi = "retina"
# )

# plot --------------------------------------------------------------------

# useful
# https://stats.stackexchange.com/a/352351/127936

bind_rows(
  onset,
  report,
  death
) %>% 
  uncount(cases) %>% 
  group_by(group) %>% 
  mutate(id = seq(1,n(),1)) %>% 
  ungroup() %>% 
  pivot_wider(
    names_from = "group",
    values_from = "date"
  ) %>% 
  incidence(
    date_index = c(
      onset = "onset",
      report = "report",
      death = "death"
    )
  ) %>% 
  mutate(
    count_variable = fct_relevel(
      count_variable,"onset","report"
    )
  ) %>% 
  plot(
    i_epiet, 
    show_cases = TRUE, 
    angle = 45, 
    n_breaks = 10
  )
```

![](https://i.imgur.com/hrJpDW5.png)<!-- -->

<sup>Created on 2024-04-16 with [reprex v2.1.0](https://reprex.tidyverse.org)</sup>