R/outbreak_summary.R

In dachuwu/DTQbp: A package for modeling case detection, contact tracing and quarantine with branching processes

#' Summarize scenario analysis
#' @author Ta-Chou Ng
#'
#' @param res a list containing the case_data data.tables generated by scenario_sim
#'
#' @return data.table of outbreak summary statistics
#' @export
#'
#' @importFrom data.table rbindlist
#' @importFrom igraph graph_from_data_frame bfs
#'
outbreak_summary <- function(dat){

  calc.nb <- function(xx){

    mu <- mean(xx)
    nllik <- function(dsp){
          -sum(dnbinom(xx, mu = mu, size = dsp, log = T))}

    res <- optim(c(1), nllik, method = "L-BFGS-B", lower = 0.01, upper = 30)
    k <- res$par

    if(F){
      pdat1 <- function(x) ecdat(xx)(x+.5)-ecdat(xx)(x-.5)
      pdat2 <- function(x)  dpois(x, mu) #dnbinom(x, mu = mu, size = k)
      xsq <- 0:15
      d <- data.frame(x = xsq, y1 = pdat1(xsq), y2 = pdat2(xsq))
      d$y1 <- d$y1/sum(d$y1)
      d$y2 <- d$y2/sum(d$y2)
      ggplot(d)+
        geom_bar(aes(x=x, y=y1), stat="identity", width = .4)+
        geom_bar(aes(x=x+.4, y=y2), stat="identity", width = .4, fill="blue3")

    }


    return(list(mu = mu, k=k))
  }

  calc.p1 <- function(xx){
    return(mean(xx>1))
  }

  total_gen <- max(dat$generation)

  sum_cluster_index <- outbreak_summary_cluster(x0 = dat$infector, x1 = dat$caseid,
                                                lasym = dat$asym,
                                                init.src=dat$caseid[dat$infector==0])

  asympid <- dat[asym & generation>=0 & generation<=max(0,total_gen-3), caseid]
  asympid <- sample(asympid, min(100,length(asympid)))
  if(length(asympid)>0){
    sum_cluster_asymp <- outbreak_summary_cluster(x0 = dat$infector, x1 = dat$caseid,
                                                  lasym = dat$asym, init.src = asympid)
  }else{
    sum_cluster_asymp <- list()
  }

  sympid <- dat[!asym & generation>=0 & generation<=max(0,total_gen-3), caseid]
  sympid <- sample(sympid, min(100,length(sympid)))
  if(length(sympid)>0){
    sum_cluster_symp <- outbreak_summary_cluster(x0 = dat$infector, x1 = dat$caseid,
                                                  lasym = dat$asym, init.src = sympid)
  }else{
    sum_cluster_symp <- list()
  }


  dat <- dat[!is.na(R0)&!is.na(Re),]

  d_re <- calc.nb(xx = dat$Re)
  d_r0 <- calc.nb(xx = dat$R0)

  p1_re <- calc.p1(dat$Re)
  p1_r0 <- calc.p1(dat$R0)

  sum_outbreak <- data.frame(mu_Re = d_re$mu, k_Re = d_re$k ,
                    mu_R0 = d_r0$mu, k_R0 = d_r0$k ,
                    pgt1_re = p1_re,
                    pgt1_r0 = p1_r0,
                    psrc_ngen_gt0=mean(sum_cluster_index$ngen>0),
                    psrc_ngen_gt1=mean(sum_cluster_index$ngen>1),
                    psrc_ngen_gt2=mean(sum_cluster_index$ngen>2),
                    pasymp_ngen_gt0=mean(sum_cluster_asymp$ngen>0),
                    pasymp_ngen_gt1=mean(sum_cluster_asymp$ngen>1),
                    pasymp_ngen_gt2=mean(sum_cluster_asymp$ngen>2),
                    psymp_ngen_gt0=mean(sum_cluster_symp$ngen>0),
                    psymp_ngen_gt1=mean(sum_cluster_symp$ngen>1),
                    psymp_ngen_gt2=mean(sum_cluster_symp$ngen>2)
                      )

  return(list(sum_outbreak = sum_outbreak, sum_cluster_index = sum_cluster_index))


}