search

epiverse-trace / epidemics

A library of published compartmental epidemic models, and classes to represent demographic structure, non-pharmaceutical interventions, and vaccination regimes, to compose epidemic scenarios.
https://epiverse-trace.github.io/epidemics/
Other
8 stars 2 forks source link

Vector inputs for ODE models #176

Closed pratikunterwegs closed 5 months ago

pratikunterwegs commented 6 months ago

This is a substantial user-facing and internal update to {epidemics}. Please treat any review as a full-package review.

Context

Changes in this PR

Planned changes not in this PR

github-actions[bot] commented 5 months ago

This pull request:

(Note that results may be inacurrate if you branched from an outdated version of the target branch.)

github-actions[bot] commented 5 months ago

This pull request:

(Note that results may be inacurrate if you branched from an outdated version of the target branch.)

github-actions[bot] commented 5 months ago

This pull request:

(Note that results may be inacurrate if you branched from an outdated version of the target branch.)

github-actions[bot] commented 5 months ago

This pull request:

(Note that results may be inacurrate if you branched from an outdated version of the target branch.)

pratikunterwegs commented 5 months ago

Thanks @jamesmbaazam for the feedback - happy to hear more on some of the points you've raised.

Other comments

  • In .check_prepare_args_default(), lots of subsetting are done in a pattern that could potentially be turned into helper functions. Can subsetting operations like mod_args[["vaccination"]][["time_begin"]] be combined into a subsetting function that take an "intervention" and "element" argument? This would reduce the hardcoding of the list elements and make the code more readable and maintainable.

  • I see get_parameter() has been removed but I see the replaced subsetting operations are identical to what get_parameter() would do. Why was it removed instead of being updated to do the new list subsetting? Could you convert all the x$<some element> into a getter function like get_parameter() previously? This would make the code more readable.

Taking these two together as they're related - I felt get_parameter() wasn't offering much over inbuilt accessors. We have a preference for list-based S3 classes, so the [[ accessor will likely remain robust. get_parameter() was also exported, which meant adding input checks, which I felt was making it slow down internal functions when it was used in them (as in .check_prepare_args()). Stripping those checks away for an internal-use version, say .get_parameter() did not seem to offer much more over [[. In .check_prepare_args_*() for example, get_parameter() use looks like get_parameter(mod_args[["some element"]], "some member").

  • Can the details section of all the check_args_*() functions be combined into a single template? They seem to be quite similar and repeated and might be easier to maintain if turned into one.

Do you mean the overall documentation for these functions named .check_prepare_args_*()? They could be combined. The Return section would be quite long, since it would repeat elements a few times, once per model in case of slight differences.

  • Can the R/dummy_elements.R script be renamed to something more descriptive? The no_time_dependence() function could be moved to the time_dependence.R script and the no_population_change() function moved somewhere with related content. This will make it easier to find the functions when needed.

I can return the dummy-elements functions to files that hold the classes and do away with this one.

  • output_to_df() should probably be renamed to output_to_dt() since it returns a data.table. The return value should also be changed from data.frame to data.table to reflect the code.

Thanks - I'll change the documented return type. I can rename this fn to .output_to_table(), as I'm not sure that we'll stick with <data.table>.

  • epidemic_size():

    • The return value needs to be updated to say it returns a data.table.

Thanks, got it.

  • Could you enhance it to return the epidemic size by a selection of age groups?

This is a nice feature in theory, but in practice it may be good to restrict it to the current implementation. Demography groups in the population do not have to be named, but might be depending on the whether the demography vector or contact matrix have names. It might be easier for users to filter a resulting data.frame of demography group, time, and size, on the names of the demography group as shown in the column. I can add an issue for epidemic_size() to return a table-like rather than a vector.

  • Could you return the "stage" as part of the output? I can easily imagine that being useful in cases where a user wants to plot the epidemic size at various stages of the outbreak.
  • I would suggest renaming "stage" to timepoint or something similar and filter by time instead of stage. We often think of epidemic timelines in terms of dates and not stages. Using stages forces you to think in terms of percentages of the whole timeline. Additionally, you could make "stage" an optional argument alongside the new "time" argument.

Yes to returning the stage as well as renaming it - collecting these feature requests in #190.

The reason epidemic_size() uses a percentage term is to quickly access e.g. the last size, the size at 50% time, etc., in a way that's agnostic to the actual time the epidemic runs. I can shift to a timepoint-based implementation (with optional time, returning final size by default) if that's more useful/more understandable. It would help to consider the expected output when a simulation does not run up to the requested time - should it error, or return a size of 0, or return NA, or return the last size available? Note that the last size might not be the 'final size' as an epidemic might not be complete by the simulation end time.

  • Is the term "demography_vector" chosen for any reason? I think in the contexts it's used here, it mostly refers to the population size per age group. Could you rename it to some like "group_pop_size" or something similar since demography is too broad a term?

I think this was adopted from {finalsize}, which adopted it from earlier implementations. I think it's not too difficult to understand so unless there's a good reason to change it, I'll keep it to conform with {finalsize}.

  • Can population_change in model_diphtheria() be structured in a manner to take a time series of population changes? For example, I might want to model a population change that occurs at a certain time point and then another change that occurs at a later time point. This would allow for more flexibility in the model.

Yes; here's an example showing both increases and decreases in certain groups, t = 70, and t = 100:

# the `values` list must always have as many elements as the `time` vector
# and each element of `values` must be of the same length as `demography_vector`
pop_change <- list(
  time = c(70, 100)
  values = list(
    change_1 = c(1e4, 1e5, 2e5),
    change_2 = c(-9e3, 1e3, -1e5)
  )
)

I considered adding this to the diphtheria vignette but this functionality might need more thinking through before a more detailed example is added.

  • Can model_vacamole() be renamed to model_covid19() or something descriptive so as to align with the other models named after the disease they represent? I can understand if you're trying to use the original name to conserve the credit (or maybe I'm not privy to the reuse agreement) but credit can be given in the function documentation and README as is currently done.

It could indeed be. One reason I would not do that is that the philosophy of {epidemics} is now more towards collecting published models that can target types of diseases rather than specific diseases, e.g., the default model was recently used for a range of 11 diseases. Vacamole might be suitable for flu etc. as well, and the name refers to the two-dose leaky vaccination scenario with a vaccination-mediated infection pathway it implements, so I think it would be good to retain it.

  • I would suggest doing the input checks outside of the model_*() functions to make the code more focused on the model rather than being overpowered with input checks.

Do you mean that in an analysis script users should run input checking on their arguments first, and then pass the checked arguments to model functions? I am not against that - it would certainly make the model functions lighter and more readable.

My feeling has been that users would appreciate self-contained functions more, where they only really have to create a <population> and pass it to a function which will do the rest. Passing a checked argument list would probably require adding default_args_model_*() functions, which users would have to modify to pass their own values. This is the approach I've taken in {noromod} for @bolthikaru. Happy to hear inputs on the benefits on this approach. Also cc-ing @TimTaylor in case interested.

jamesmbaazam commented 5 months ago

Thanks @jamesmbaazam for the feedback - happy to hear more on some of the points you've raised.

Other comments

  • In .check_prepare_args_default(), lots of subsetting are done in a pattern that could potentially be turned into helper functions. Can subsetting operations like mod_args[["vaccination"]][["time_begin"]] be combined into a subsetting function that take an "intervention" and "element" argument? This would reduce the hardcoding of the list elements and make the code more readable and maintainable.
  • I see get_parameter() has been removed but I see the replaced subsetting operations are identical to what get_parameter() would do. Why was it removed instead of being updated to do the new list subsetting? Could you convert all the x$<some element> into a getter function like get_parameter() previously? This would make the code more readable.

Taking these two together as they're related - I felt get_parameter() wasn't offering much over inbuilt accessors. We have a preference for list-based S3 classes, so the [[ accessor will likely remain robust. get_parameter() was also exported, which meant adding input checks, which I felt was making it slow down internal functions when it was used in them (as in .check_prepare_args()). Stripping those checks away for an internal-use version, say .get_parameter() did not seem to offer much more over [[. In .check_prepare_args_*() for example, get_parameter() use looks like get_parameter(mod_args[["some element"]], "some member").

Thanks for the explanation. It makes sense to me

  • Can the details section of all the check_args_*() functions be combined into a single template? They seem to be quite similar and repeated and might be easier to maintain if turned into one.

Do you mean the overall documentation for these functions named .check_prepare_args_*()? They could be combined. The Return section would be quite long, since it would repeat elements a few times, once per model in case of slight differences.

Yes, I meant the overall documentation. I suggested that after noticing the patterns, but it might not make sense to combine them.

  • Can the R/dummy_elements.R script be renamed to something more descriptive? The no_time_dependence() function could be moved to the time_dependence.R script and the no_population_change() function moved somewhere with related content. This will make it easier to find the functions when needed.

I can return the dummy-elements functions to files that hold the classes and do away with this one.

Yes, it would make it easier to navigate the code base.

  • output_to_df() should probably be renamed to output_to_dt() since it returns a data.table. The return value should also be changed from data.frame to data.table to reflect the code.

Thanks - I'll change the documented return type. I can rename this fn to .output_to_table(), as I'm not sure that we'll stick with <data.table>.

Alright.

  • epidemic_size():
  • The return value needs to be updated to say it returns a data.table.

Thanks, got it.

  • Could you enhance it to return the epidemic size by a selection of age groups?

This is a nice feature in theory, but in practice it may be good to restrict it to the current implementation. Demography groups in the population do not have to be named, but might be depending on the whether the demography vector or contact matrix have names. It might be easier for users to filter a resulting data.frame of demography group, time, and size, on the names of the demography group as shown in the column. I can add an issue for epidemic_size() to return a table-like rather than a vector.

  • Could you return the "stage" as part of the output? I can easily imagine that being useful in cases where a user wants to plot the epidemic size at various stages of the outbreak.
  • I would suggest renaming "stage" to timepoint or something similar and filter by time instead of stage. We often think of epidemic timelines in terms of dates and not stages. Using stages forces you to think in terms of percentages of the whole timeline. Additionally, you could make "stage" an optional argument alongside the new "time" argument.

Yes to returning the stage as well as renaming it - collecting these feature requests in #190.

Great

It would help to consider the expected output when a simulation does not run up to the requested time - should it error, or return a size of 0, or return NA, or return the last size available?

It is often a difficult situation but sticking to a simple and easy-to-debug solution for now might do for a first pass. It can be improved with user feedback. I'm inclined to suggest returning the last size with a warning and some text on how to interpret it.

Note that the last size might not be the 'final size' as an epidemic might not be complete by the simulation end time.

Wouldn't this be the case for stage = 1 too?

  • Is the term "demography_vector" chosen for any reason? I think in the contexts it's used here, it mostly refers to the population size per age group. Could you rename it to some like "group_pop_size" or something similar since demography is too broad a term?

I think this was adopted from {finalsize}, which adopted it from earlier implementations. I think it's not too difficult to understand so unless there's a good reason to change it, I'll keep it to conform with {finalsize}.

  • Can population_change in model_diphtheria() be structured in a manner to take a time series of population changes? For example, I might want to model a population change that occurs at a certain time point and then another change that occurs at a later time point. This would allow for more flexibility in the model.

Yes; here's an example showing both increases and decreases in certain groups, t = 70, and t = 100:

# the `values` list must always have as many elements as the `time` vector
# and each element of `values` must be of the same length as `demography_vector`
pop_change <- list(
  time = c(70, 100)
  values = list(
    change_1 = c(1e4, 1e5, 2e5),
    change_2 = c(-9e3, 1e3, -1e5)
  )
)

Ah, nice!

I considered adding this to the diphtheria vignette but this functionality might need more thinking through before a more detailed example is added.

  • Can model_vacamole() be renamed to model_covid19() or something descriptive so as to align with the other models named after the disease they represent? I can understand if you're trying to use the original name to conserve the credit (or maybe I'm not privy to the reuse agreement) but credit can be given in the function documentation and README as is currently done.

It could indeed be. One reason I would not do that is that the philosophy of {epidemics} is now more towards collecting published models that can target types of diseases rather than specific diseases, e.g., the default model was recently used for a range of 11 diseases. Vacamole might be suitable for flu etc. as well, and the name refers to the two-dose leaky vaccination scenario with a vaccination-mediated infection pathway it implements, so I think it would be good to retain it.

That makes sense for the generic models but there are other models named after diseases here as well. No?

  • I would suggest doing the input checks outside of the model_*() functions to make the code more focused on the model rather than being overpowered with input checks.

Do you mean that in an analysis script users should run input checking on their arguments first, and then pass the checked arguments to model functions? I am not against that - it would certainly make the model functions lighter and more readable.

Oh no, I see I didn't phrase that well. I meant that, for readability, you should consider bundling up the checks into helper functions as is done in our other packages like simulist. But this might be a matter of style.

Your suggestion here is also good but could potentially lead to downstream issues if you don't put in measures to ensure the user checks their inputs first.

My feeling has been that users would appreciate self-contained functions more, where they only really have to create a <population> and pass it to a function which will do the rest. Passing a checked argument list would probably require adding default_args_model_*() functions, which users would have to modify to pass their own values. This is the approach I've taken in {noromod} for @bolthikaru. Happy to hear inputs on the benefits on this approach. Also cc-ing @TimTaylor in case interested.

pratikunterwegs commented 5 months ago

Note that the last size might not be the 'final size' as an epidemic might not be complete by the simulation end time.

Wouldn't this be the case for stage = 1 too?

Yes, it's simply the last size in the time series. Some user care is required in determining whether it is really the 'final size', after looking at the trajectory.

That makes sense for the generic models but there are other models named after diseases here as well. No?

There are - we don't have good naming options for these. We have not made much progress in implementing more models, so general rules for model names have not really been developed.

pratikunterwegs commented 5 months ago

Do you mean that in an analysis script users should run input checking on their arguments first, and then pass the checked arguments to model functions? I am not against that - it would certainly make the model functions lighter and more readable.

Oh no, I see I didn't phrase that well. I meant that, for readability, you should consider bundling up the checks into helper functions as is done in our other packages like simulist. But this might be a matter of style.

This could indeed be done. In {finalsize} some of the discussion were around how much to show this checking and preparation, and the idea was to strike a balance to prevent the codebase spreading over too many files. To some extent this is the role of the .check_prepare_args_*() functions.

pratikunterwegs commented 5 months ago

Can the R/dummy_elements.R script be renamed to something more descriptive? The no_time_dependence() function could be moved to the time_dependence.R script and the no_population_change() function moved somewhere with related content. This will make it easier to find the functions when needed.

I can return the dummy-elements functions to files that hold the classes and do away with this one.

These elements don't have a better place to go, as they are not classes.

github-actions[bot] commented 5 months ago

This pull request:

(Note that results may be inacurrate if you branched from an outdated version of the target branch.)

github-actions[bot] commented 5 months ago

This pull request:

(Note that results may be inacurrate if you branched from an outdated version of the target branch.)

pratikunterwegs commented 5 months ago

Thanks both @jamesmbaazam and @adamkucharski for your reviews. @CarmenTamayo, you mentioned you might have some feedback as well, would you like to add it? Otherwise I think this PR is ready to merge.

CarmenTamayo commented 5 months ago

Thanks both @jamesmbaazam and @adamkucharski for your reviews. @CarmenTamayo, you mentioned you might have some feedback as well, would you like to add it? Otherwise I think this PR is ready to merge.

Hi @pratikunterwegs I provided some feedback as well, specifically for the new vignette for modelling parameter uncertainty, can you see my comments? I did this last week -> I hadn't submitted the comments but they're available now- apologies for this

pratikunterwegs commented 5 months ago

Sorry Pratik, I had pending comments but hadn't actually submitted them

Thanks! I remembered you'd asked me questions about some parts of this PR so I knew there must be comments somewhere!

github-actions[bot] commented 5 months ago

This pull request:

(Note that results may be inacurrate if you branched from an outdated version of the target branch.)

github-actions[bot] commented 5 months ago

This pull request:

(Note that results may be inacurrate if you branched from an outdated version of the target branch.)

pratikunterwegs commented 5 months ago

Thanks all for your reviews - merging this now so we can make smaller PRs in future for related issues.