epiverse-trace / finalsize Goto Github PK

Implement code for final_size with susceptibility groups. This is the R-only equivalent of the Rcpp code in #39.

The idea here is to provide finalsize with a function exported to users, and with other internal functions and documentation, that allow the calculation of the final size of an epidemic in a population with heterogeneous mixing and also heterogeneous susceptibility to infection.

This is to be achieved by a translation of the Cpp code here https://gitlab.com/epidemics-r/code_snippets/-/blob/feature/newton_solver/include/finalsize.hpp#L20 into R.

This issue includes the following sub-issues: #44 #45 #46 #47 #48

It would be very useful to test that errors / warnings / messages (currently, I think only errors) do happen when they are supposed to, with the correct message. See this for an example.

Originally posted by @thibautjombart in #19 (comment)

No bleeding of infections from one age group to the other if the contact rate is set to 0

Originally posted by @Bisaloo in #32 (review)

Iterative solver loop from L. 42 onwards can probably be vectorised which would be quicker if so.

Originally posted by @TimTaylor in #42 (comment)

There are use cases for using both R0 and transmissibiltiy, so it would be good to have functions to convert between these.

In fluevidencesynthesis it was done: https://github.com/MJomaba/flu-evidence-synthesis/blob/c5cbd127ff717f3637c5ed298dd3428e0fa6b05d/src/rapi.cc#L740-L787

This approach should only need minor changes for our needs.

An implementation of final_size called final_size_cpp that uses the Newton solver from @BlackEdder here: https://gitlab.com/epidemics-r/code_snippets/-/blob/feature/newton_solver/include/finalsize.hpp with an example here: https://gitlab.com/epidemics-r/code_snippets/-/blob/feature/newton_solver/tests/catch_final_size.cpp.

The initial proportions of each age group infected in final_size are hard coded to a numeric. Suggest instead a function argument that allows this proportion to be specified by the user. Could be a single value for all age groups or one per age group.
This issue addresses @thibautjombart in #19 (comment)

Not really an issue - should we start a discussion board?

Just to mention that some related (probably simpler) functionality to what I think is being implemented here is available at https://github.com/sbfnk/epimixr/blob/master/R/epidemic_age_dist.r - haven't looked yet in detail how the different methods are related but either way there might be useful code/inspiration there.

There are related functions there for estimating reproduction numbers and projecting immunity from seroprevalence/vaccination data. These could one day be brought under the same hat and/or linked up.

Replace double rbind statements with:

  # replicate contact matrix
  contact_matrix_ <- kronecker(
    X = matrix(1, nrow = n_susc_groups, ncol = n_susc_groups),
    Y = contact_matrix
  )

to replicate contact matrix in epi_spread

Originally posted by @TimTaylor in #42 (comment)

Integrate code from @BlackEdder at https://gitlab.com/epidemics-r/code_snippets/-/blob/feature/newton_solver/include/finalsize.hpp into package to allow calculation of final size with multiple susceptibility groups.

Compare speeds for final_size and final_size_cpp.

Final size is currently noted as pi in the code, which conflicts with R native pi value.

Mostly likely fine as I doubt you need $\pi$ the mathematical constant, but by calling this object pi you're overwriting the default constant loaded by R

Originally posted by @joshwlambert in #42 (comment)

Add unit tests for the functions f1 and f2 which are copied from final_size.R

Adding functionality to handle populations with non uniform distribution of susceptibilities in current finalsize function final_size in R.

This will add a p_susceptibility (or similar) argument to final_size.

End point is a function whose functionality is in line with the proposed final_size_cpp function whose underlying Cpp code already implements heterogenous susceptibility distributions across age groups.

Writing vignettes for final_size use, which is expected to be the single function in finalsize

Make pkgdown website for finalsize; must include vignettes from #9

Designing an output format for finalsize that is shared by outputs from the ODE-based simulations in the upcoming epidemics package. Implementation likely to be as S3 and may be moved to another package.

final_size returns final epidemic sizes of 0.0 for any prop_initial_infected (> 0.0) when r0 < 1.0. This is implausible as some individuals are already initially infected. The algorithm implemented here (most likely) does not work for r0 < 1.0, and final_size should stop when r0 < 1.0 is provided.

Add basic checks for function argument assumptions in final_size.

Here for instance, there seems to be an assumption that contact_matrix is a matrix, but this is not formally checked.

Originally posted by @thibautjombart in #19 (comment)

stopifnot allows multiple expressions seperated by commas. That would probably be a little more efficient

Originally posted by @TimTaylor in #42 (comment)

Remove Mac specific .DS_Store file.

Allow users to specify a single prop_suscep for all age groups.

Add tests for final_size_cpp, the Rcpp implementation of final_size.
Check:

1. Basic functionality with synthetic data,
2. Functionality with polymod data,
3. Increasing r0 leads to larger epidemics,
4. Passing incorrect length vector arguments throws specific errors.

Structure the finalsize Cpp Newton solver code as a header in an include directory so that it can be called from other Rcpp projects similar to the BH (Boost Headers) package (which provides the Boost libraries as headers).

Function solve_final_size_iterative must also take adapt_step and step_rate arguments passed from wrapping function.

Originally posted by @pratikunterwegs in #42 (comment)

Main effort to implement a Cpp Newton solver for finalsize is catalogued under a separate project here https://github.com/orgs/epiverse-trace/projects/2

Iterative solver function fn_f looks like it can also be vectorised

Originally posted by @TimTaylor in #42 (comment)

Check that final_size returns a vector with values in the range prop_initial_infected and 1.0 for r0 > 1.0.