Add functions to simulate datasets to the seroprevalence_data module.…

[ntorresd]

This PR partly closes #57

• fix: add exception to function plot_foi() to plot a FOI trend with different length along with the data for the case when their sizes don't coincide

• clean test_visualisation

• feature: add three functions to simulate datasets. get_sim_counts() generates a list with simulated counts by age following a binomial distribution. generate_sim_data() uses the counts generated by get_sim_counts() to create a dataframe with the necessary structure to use other functions of the package. group_sim_data() serves to group the previously generated dataset by age group; right now it groups the data by periods of 5 years.

• add test_simulate_data to test the data simulation functions in the seroprevalence_data module

[ben18785]

Thanks @ntorresd -- looking better. We still need much more extensive tests of the functions though as these functions are really key to lots of other functionality (including how we test that the Stan models are themselves working):

• The get_sim_probability function needs value-based-testing. I.e. do the probabilities it outputs match theoretically derived ones?
• Same for get_sim_data.
• There are a variety of internal functions that also need checking, which I think are currently lacking unit tests.

[ben18785]

Minor thing but can we move the group_by to the next line for consistency when using the pipe?

[ben18785]

I know this function is stochastic but we can still test it works by varying the FOIs and checking that in limits it behaves as we'd like it to.

[ntorresd]

Right now I'm relaying on model implementation to test this functionality. My idea is to make use of suitable models for different FOI trends. It is to be expected for a constant foi_simto be correctly approximated by the "constant" model for all times, so what I do is to make sure that this lies in the confidence interval obtained by implementing this model:

# Define constant FoI for simulations
case_label <- "constant_foi_"
foi_model <- "constant"
foi_sim <- rep(0.02, tsur - birth_year_min)
max_lambda <- 0.035

# Generate simulated data and run the "constant" model
sim_data <- generate_sim_data(foi = foi_sim,
                              sample_size_by_age = sample_size_by_age,
                              tsur = tsur,
                              birth_year_min = birth_year_min,
                              survey_label = 'foi_sim',
                              seed = seed)
sim_seromodel <- run_seromodel(sim_data, foi_model = foi_model, n_iters = n_iters)

# Check consistency between sim_foi and the fitted foi
foi <- rstan::extract(sim_seromodel$seromodel_fit, "foi", inc_warmup = FALSE)[[1]]
foi_lower <- apply(foi, 2, function(x) quantile(x, 0.05))
foi_upper <- apply(foi, 2, function(x) quantile(x, 0.95))
expect_true(all((foi_sim >= foi_lower) & (foi_sim <= foi_upper)))

I implemented this for a constant FOI and for the following smooth-decreasing FOI (red lines in the images):

case_label <- "smth_dec_foi_" #Smooth-decendent FoI
foi_model = "tv_normal"
foi_max = 0.2
stretch = 0.15
x <- 1:(tsur - birth_year_min)
foi_sim <- (-foi_max * (atan(stretch * (x - 25))) / (0.5 * pi) + foi_max) / 2

I think it's worth to open a separate issue to discuss this in more detail, in case that you're not convinced by my approach.

[ben18785]

We should be able to test whether data we obtain from simulating from these models is as expected without resorting to solving the inverse problem (as this leaves us liable to a number of issues, e.g. the FOIs aren't identifiable given the data). We have analytical results for all of the models which should allow us to test these directly. E.g. for a constant FOI model then we will know approximately the probability that someone aged x is seropositive; if we use a large sample size, we can check that the simulated proportion is near to that value.