aggregate `<epichains>` object by time per chain

[avallecam]

Hi @jamesmbaazam, I found that it is possible to aggregate() multiple chains within one simulation. However, Is it possible to have time aggregates for each chain in that same simulation?

From the example at https://epiverse-trace.github.io/epichains/reference/aggregate.epichains.html, I generated two alternatives, but I am wondering if epichains already has an optimal way to do this without using any data wrangling (neither tidy nor base).

My goal is to get a similar output as in the tutorial episode introduction figure: https://epiverse-trace.github.io/tutorials-middle/superspreading-simulate.html and minimize as much as possible the data wrangling steps in the tutorial

Reprex:

library(epichains)
library(tidyverse)

set.seed(32)

chains <- simulate_chains(
  n_chains = 10,
  statistic = "size",
  offspring_dist = rpois,
  stat_threshold = 10,
  generation_time = function(n) rep(3, n),
  lambda = 2
)

# aggregate all chains
aggregate(chains, by = "time")
#>   time cases
#> 1    0    10
#> 2    3    25
#> 3    6    61
#> 4    9    33
# chains %>% filter(chain == 1) # filter keeps the <epichain> class object

# alternative 1: aggregate per chain
alt_1 <- expand_grid(
  chain_seq = 1:10,
  epichains_input = list(chains)
) %>% 
  rowwise() %>% 
  mutate(epichains_filter = list(filter(epichains_input, chain == chain_seq))) %>% 
  mutate(epichains_aggregate = list(aggregate(epichains_filter, by = "time"))) %>% 
  mutate(epichains_cumsum = list(mutate(epichains_aggregate, cases_cumsum = cumsum(cases)))) %>% 
  unnest(epichains_cumsum)

# alternative 2: aggregate per chain
alt_2 <- chains %>% 
  as_tibble() %>%
  mutate(day = ceiling(time)) %>%
  count(chain, day, name = "cases") %>%
  group_by(chain) %>%
  mutate(cases_cumsum = cumsum(cases)) %>%
  ungroup()

# compare
waldo::compare(
  alt_1 %>% dplyr::select(chain = chain_seq, time:cases_cumsum),
  alt_2 %>% dplyr::rename(time = day)
)
#> ✔ No differences

# output of alternative 1
alt_1 %>% dplyr::select(chain = chain_seq, time:cases_cumsum)
#> # A tibble: 34 × 4
#>    chain  time cases cases_cumsum
#>    <int> <dbl> <int>        <int>
#>  1     1     0     1            1
#>  2     1     3     2            3
#>  3     1     6     3            6
#>  4     1     9     6           12
#>  5     2     0     1            1
#>  6     2     3     2            3
#>  7     2     6     5            8
#>  8     2     9    14           22
#>  9     3     0     1            1
#> 10     3     3     3            4
#> # ℹ 24 more rows

^{Created on 2025-04-04 with reprex v2.1.1}

Desired output, that was not able to render in the same reprex:

# goal: visualize grouped chain trajectories
alt_1 %>%
ggplot(aes(
  x = time,
  y = cases_cumsum,
  group = chain_seq
)) +
geom_line() +
labs(
  x = "Day",
  y = "Cumulative cases"
)

[jamesmbaazam]

Hi @avallecam It is currently not possible to do that internally with aggregate() but I could add an argument like by_chain = TRUE/FALSE so that offspring are aggregated by time or generation per chain. Do let me know if that's what you mean.

[avallecam]

Hi @avallecam It is currently not possible to do that internally with aggregate() but I could add an argument like by_chain = TRUE/FALSE so that offspring are aggregated by time or generation per chain. Do let me know if that's what you mean.

thanks @jamesmbaazam, yes, that would be a good fit.

In order to provide an argument as flexible as possible, does doing an aggregate(x = chains, by = c("chain", "time")) sounds appropriate? (edit: but I can adapt to the easiest solution)

pasting output as a reference of the columns available:

chains
#> `<epichains>` object
#> 
#> < epichains head (from first known infector) >
#> 
#>    chain infector infectee generation time
#> 11     1        1        2          2    3
#> 12     1        1        3          2    3
#> 13     2        1        2          2    3
#> 14     2        1        3          2    3
#> 15     3        1        2          2    3
#> 16     3        1        3          2    3
#> 
#> 
#> Number of chains: 10
#> Number of infectors (known): 9
#> Number of generations: 4
#> Use `as.data.frame(<object_name>)` to view the full output in the console.

[jamesmbaazam]

My suggested option is straightforward to implement with very little code as its just an ifelse ladder. Moreover, aggregate() doesn't have many arguments for it to be an issue.

I'm just trying to think through your suggestion of by = c("chain", "time"/generation):

• Does it mean you want by to be either c(chain, time) or c(chain, generation) if aggregating by chain and the second element? That would require more typing on the user side and is prone to mistakes. It'll require several checks to ensure they arguments are well specified and will expand the codebase with no significant benefit.

• It might cause users to think they can specify it in any order or want to try out other orders like c(time, generation).

If you don't feel strongly about it, I can implement the first and we can see if it causes any problems.

[avallecam]

Thanks @jamesmbaazam I agree with your assessment. It is misleading to have c(time, generation). Agree with the simple ifelse you propose.

[jamesmbaazam]

I'll go ahead and implement it.

[avallecam]

hi @jamesmbaazam thanks for adding this to short-term. I'm wondering if this feature request could also include the plot() function, or better to fill in a separate issue? I just noticed that the output of aggregate could be nicely printed by plot() from the last part of the getting started vignette.