R/dictionary.R

In yairgoldy/sl4hm: Simulated likelihood for modelling household epidemics

#' For each family type and duration calculates the distribution on the number of positives
#' @param seed seed for the simulation
#' @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_cc proabilty that a child will infect a child
#' @param beta_ac proabilty that a child will infect an adult
#' @param beta_ca proabilty that an adult will infect a child
#' @param beta_aa proabilty that an adult will infect an adult
#' @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 ssn the number of households to simuluate
#'
#' @return a tibble with prabaility to see number of positive adults and number of postive children as a function of the duration and the index case
#' @export
#' @import magrittr
#' @import data.table
#' @importFrom tibble as_tibble
#' @importFrom dplyr mutate select group_by summarise ungroup arrange inner_join rename
#' @importFrom tidyr spread
#' @importFrom purrr pmap
#' @seealso `load_P` `create_family_structure_data`
family_dictionary <- function(seed = -1, family_structure_data,beta_cc,beta_ac,beta_ca,beta_aa,
                              P= NULL,T_max=45,ssn=1000){
  beta <- matrix(c(beta_cc,beta_ca,beta_ac,beta_aa),2,2,byrow =T)

  family_type_dat <- dplyr::rename(family_structure_data,true_day = day) %>%
    dplyr::mutate(seed = 1000*seed+1:dplyr::n(), family_prob = 1)

  res <- dplyr::mutate(family_type_dat,
                       sims = purrr::pmap(list(adults,childs,first.is.adult,seed),
                                          function(adultsi,childsi,first.is.adulti,seedi){
                                            simulate_household(seed = seedi,adultsi,childsi,first.is.adulti,
                                                               beta_cc,beta_ac,beta_ca,beta_aa,P,T_max,ssn)}))

  res <- dplyr::mutate(res,
                       probs =  purrr::pmap(list(adults,childs,first.is.adult,sims),
                                            function(adultsi,childsi,first.is.adulti,simi){
                                              prob_table(adultsi,childsi,first.is.adulti,simi,T_max,ssn)})) %>%
    dplyr::select(-sims,-seed,-family_prob)

  res <- tidyr::unnest(res,cols = "probs") %>%
    tidyr::gather( day, prob, -c(1:6)) %>%
    dplyr::mutate(prob  = ifelse(prob ==0,1e-4,prob),
           logprob = log(prob),
           day = as.numeric(day))

  return(res)
}



prob_table <- function(adults,childs,first.is.adult,sim,T_max,ssn){

  if(length(sim)== 0) return(NULL)
  if(first.is.adult){
    prob_dat <- tidyr::expand_grid(pos.adults= 0:(adults-1),pos.childs = 0:childs)
  } else {
    prob_dat <- tidyr::expand_grid(pos.adults= 0:adults,pos.childs = 0:(childs-1))
  }


  for(i in 10:T_max){
    infi <- dplyr::filter(sim, day ==i)
    day <- as.character(i)
    prob_dat <- dplyr::mutate(prob_dat,!!day := purrr::map2_dbl(pos.adults,pos.childs,
                                                                function(a,c){sum(a==infi$pos.adults & c== infi$pos.childs)/ssn}))

  }

  return(prob_dat)

}