R/simulation.R

In yairgoldy/sl4hm: Simulated likelihood for modelling household epidemics

if(getRversion() >= "2.15.1")
  utils::globalVariables(c("scenario", "data","simseed","index","type"))

#' Simulate nsim simulations of an epidemic on a dataset of households
#' @param nsim number of simulation to run
#' @param family_structure_data dataset of families including the number of adults, number of children, and the index case (see `create_family_structure_data`)
#' @param beta_aa proabilty that an adult will infect an adult (must be appear in the dictionary)
#' @param gamma_val relative susceptibility of children to adults (must be appear in the dictionary)
#' @param delta_val relative infectivity of children to adults (must be appear in the dictionary)
#' @param P a vector of generation time distribution (default is generated by `load_P``)
#' @param T_max the length of the epidemic within the family
#' @param future_plan the plan for using the future package (default is `sequential`)
#' @param dictionary_path path to the folder that contain the dictionaries
#' @return a tibble with the maximum likelihood for each simulated dataset
#' @export
#' @import magrittr
#' @importFrom tibble tibble
#' @importFrom purrr safely map
#' @importFrom dplyr mutate select group_by summarise ungroup arrange inner_join sample_frac
#' @importFrom tidyr  spread  gather nest unnest
#' @importFrom future plan
#' @importFrom furrr future_map
#' @importFrom stats  pgamma rbinom
#' @seealso `load_P` `create_family_structure_data`
run_simulation <- function(nsim = 1, family_structure_data,
                           beta_aa,gamma_val,delta_val,
                           P =NULL,
                           T_max = 45,
                           future_plan = "sequential",
                           dictionary_path = NULL){
  future::plan(future_plan)

  results <- tibble::tibble(scenario = 1:nsim,true_beta_aa= beta_aa,
                            true_gamma_val =gamma_val, true_delta_val = delta_val, simseed= scenario)

 sim_fn <- purrr::safely(calc_sim_lik)
  results <- results %>%
    dplyr::group_by(scenario) %>%
    tidyr::nest() %>%
    dplyr::ungroup() %>%
    ##Randomise the order of scenarios - helps share the load across cores
    dplyr::sample_frac(size = 1, replace = FALSE) %>%
    dplyr::mutate(nll =
                    furrr::future_map(data,
                                      ~sim_fn(
                                        seed = .$simseed,
                                        beta_aa = .$true_beta_aa,
                                        gamma_val = .$true_gamma_val,
                                        delta_val = .$true_delta_val,
                                        family_structure_data = family_structure_data,
                                        P = P,
                                        T_max = T_max,
                                        dictionary_path = dictionary_path),.progress = TRUE
                    )) %>%
    tidyr::unnest(cols = "data") %>%
    dplyr::select(-simseed)


  sres <- dplyr::mutate(results,nll = purrr::map(nll,function(nlli){return(nlli$result)})) %>%
    tidyr::unnest(nll) %>%
    dplyr::select(-scenario)
  future::plan("sequential")
  return(sres)
}

calc_sim_lik <- function(seed = -1, family_structure_data,P,beta_aa,gamma_val,delta_val,T_max = 45,dictionary_path){

  beta_cc <- gamma_val*delta_val*beta_aa
  beta_ca <- gamma_val*beta_aa
  beta_ac <- delta_val*beta_aa

  dat <- create_simulated_data(seed = seed , family_structure_data,beta_cc,beta_ca,beta_ac,beta_aa,P,T_max,ssn =1)
  res <- calc_likelihood(dat,dictionary_path)

  res <-  dplyr::filter(res,nll <= min(nll))
  return(res)
}