assess to request using discrete distributions for generation time or discretize internally

[avallecam]

This is regarding the epiparameter-epichains interoperability. The current feature could lead to getting unwanted aggregate() outputs. I think {cfr} does internal discretisation when epiparameter is used as input.

Find below a reprex of the current behaviour with a default continuous and discretized epiparameter:

# default: continuous epiparameter ---------------------------------------

epiparameter_default <- epiparameter::epiparameter_db(
  disease = "mers",
  epi_name = "serial",
  single_epiparameter = TRUE
)
#> Using Cowling B, Park M, Fang V, Wu P, Leung G, Wu J (2015). "Preliminary
#> epidemiological assessment of MERS-CoV outbreak in South Korea, May to
#> June 2016." _Eurosurveillance_.
#> doi:10.2807/1560-7917.es2015.20.25.21163
#> <https://doi.org/10.2807/1560-7917.es2015.20.25.21163>.. 
#> To retrieve the citation use the 'get_citation' function

mers_chains <- epichains::simulate_chains(
  n_chains = 1000,
  statistic = "size",
  offspring_dist = rnbinom,
  mu = 0.6,
  size = 0.02,
  generation_time = as.function(epiparameter_default, func_type = "generate")
)

head(aggregate(mers_chains,by = "time"))
#>       time cases
#> 1 0.000000  1000
#> 2 6.196621     1
#> 3 6.804927     1
#> 4 6.872387     1
#> 5 7.041034     1
#> 6 7.160494     1


# discrete epiparameter --------------------------------------------------

epiparameter_discrete <- epiparameter::discretise(epiparameter_default)

mers_chains <- epichains::simulate_chains(
  n_chains = 1000,
  statistic = "size",
  offspring_dist = rnbinom,
  mu = 0.6,
  size = 0.02,
  generation_time = as.function(epiparameter_discrete, func_type = "generate")
)

head(aggregate(mers_chains,by = "time"))
#>   time cases
#> 1    0  1000
#> 2    6     1
#> 3    7     6
#> 4    8    21
#> 5    9    40
#> 6   10    56

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

[jamesmbaazam]

Thanks, Andree.

Quite interesting looking at the results from the discretisation. I've not looked into how the discretise method works in epiparameter but it seems to be rounding up the times to the next time point? That will introduce some kind of censoring and affect downstream analyses. An option would maybe to give the user control over whether they want to round up or down or only retain the whole part of the time, so that 6.9 is not different from 6.1.

@sbfnk Thoughts here?

[sbfnk]

Which part of this is unexpected? You're passing different generation times and get different outputs (in particular: you pass a continuous distribution and get continuous time vs. you pass a discrete distribution and get discrete time).

[avallecam]

Which part of this is unexpected?

My motivation to create this issue was that I expect discrete outputs when using aggregate(mers_chains,by = "time"). Though I found that this depends on the generation time distribution (continuous or discrete).

Then:

• Should we always have discrete outputs independent of the input generation time? or
• Should we suggest that the input always be discrete? or
• keep it open as it currently is

[sbfnk]

Should we always have discrete outputs independent of the input generation time? or

I don't think so - sometimes you might want to have the exact infection times. Your use case of plotting the time series is a specific one but I think there might be others.

Should we suggest that the input always be discrete?

I don't think we need to - "true" generation times are continuous and so is the underlying model.

keep it open as it currently is

I think so. Just noting that in your example you could add a processing step before aggregating where you convert the continuous infection times into discrete days (which is statistically the better way than discretising the generation time). Something like

mers_chains$day <- floor(mers_chains$time)
aggregate(mers_chains,by = "day")

[avallecam]

Just noting that in your example you could add a processing step before aggregating where you convert the continuous infection times into discrete days (which is statistically the better way than discretising the generation time). Something like

The suggestion looks appropriate. May this lead to making aggregate() open to other columns?

# default: continuous epiparameter ---------------------------------------

epiparameter_default <- epiparameter::epiparameter_db(
  disease = "mers",
  epi_name = "serial",
  single_epiparameter = TRUE
)
#> Using Cowling B, Park M, Fang V, Wu P, Leung G, Wu J (2015). "Preliminary
#> epidemiological assessment of MERS-CoV outbreak in South Korea, May to
#> June 2016." _Eurosurveillance_.
#> doi:10.2807/1560-7917.es2015.20.25.21163
#> <https://doi.org/10.2807/1560-7917.es2015.20.25.21163>.. 
#> To retrieve the citation use the 'get_citation' function

mers_chains <- epichains::simulate_chains(
  n_chains = 1000,
  statistic = "size",
  offspring_dist = rnbinom,
  mu = 0.6,
  size = 0.02,
  generation_time = as.function(epiparameter_default, func_type = "generate")
)

head(aggregate(mers_chains,by = "time"))
#>       time cases
#> 1 0.000000  1000
#> 2 6.265476     1
#> 3 6.523928     1
#> 4 6.533494     1
#> 5 6.700116     1
#> 6 6.919337     1


# alternative -------------------------------------------------------------

class(mers_chains)
#> [1] "epichains"  "data.frame"

mers_chains$day <- floor(mers_chains$time)

aggregate(mers_chains, by = "day")
#> Error in match.arg(by): 'arg' should be one of "time", "generation"

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

[sbfnk]

Ah yes, maybe -- @jamesmbaazam what do you think? Can we defer to aggegregate.data.frame()? I for one was confused (as I had forgotten) that this has reduced functionality compared to the base R version.

As a workaround for now you could do

mers_chains$time <- floor(mers_chains$time)

[jamesmbaazam]

Ah yes, maybe -- @jamesmbaazam what do you think? Can we defer to aggegregate.data.frame()? I for one was confused (as I had forgotten) that this has reduced functionality compared to the base R version.

As a workaround for now you could do

mers_chains$time <- floor(mers_chains$time)

An <epichains> object with extra columns is still the same class. The class is only invalidated when you remove the protected columns, noting that the time column is optional and only appears when you specify the generation_time argument. I created an aggregate.epichains method as syntactic sugar for the purpose of creating a time series based on onset times or generation of infection (outlined in the design vignette), so your suggestion to modify the time column as in mers_chains$time <- floor(mers_chains$time) in order to use the aggregate.epichains method makes sense to me. This is more of an opinionated workflow but I'm happy to consider modifications.

@sbfnk I'm quite reluctant to defer the method to data.frame because the method will have a confusing scope (sometimes creates a time series, if the right argument is specified, and sometimes creates a user defined aggregation).

My question is, in terms of deferring to aggregate.data.frame, should we still retain the opinionated scope of creating a time series? How will we ensure that given that the user can supply any column now?

[sbfnk]

That makes sense - if you prefer being opinionated here I don't feel strongly.