In mrc-ide/reestimate_covidIFR_analysis: Reestimation of IFRs using Serology Data

knitr::opts_chunk$set(echo = FALSE, warning = FALSE, 
                      message = FALSE, eval = TRUE, 
                      results = 'hide', fig.align = 'center', dpi = 200)
knitr::opts_knit$set(root.dir = here::here())

library(COVIDCurve)
library(tidyverse)
source("R/my_themes.R")

quick_sim_sero_diagnostics <- function(modout, spec, sens) {
  maxmachain <- drjacoby::plot_par(modout$mcmcout, modout$inputs$IFRmodel$maxMa, display = FALSE)
  spechain <- drjacoby::plot_par(modout$mcmcout, "spec", display = FALSE)
  senschain <- drjacoby::plot_par(modout$mcmcout, "sens", display = FALSE)
  modchain <- drjacoby::plot_par(modout$mcmcout, "mod", display = FALSE)
  sodchain <- drjacoby::plot_par(modout$mcmcout, "sod", display = FALSE)
  seroratechain <- drjacoby::plot_par(modout$mcmcout, "sero_con_rate", display = FALSE)

  maxmachain <- maxmachain[[1]][["trace"]] + 
    geom_hline(yintercept = 0.2, linetype = "dashed", size = 1.25) +
    theme(legend.position = "none")
  spechain <- spechain[[1]][["trace"]] + geom_hline(yintercept = spec, linetype = "dashed", size = 1.25) +
    labs(caption = paste0("Prior was Beta(",
                          paste(modout$inputs$IFRmodel$paramdf[modout$inputs$IFRmodel$paramdf$name == "spec", c("dsc1", "dsc2")], collapse = ","),
                          ")")) +
    theme(legend.position = "none")
  senschain <- senschain[[1]][["trace"]] + geom_hline(yintercept = sens, linetype = "dashed", size = 1.25) +
    labs(caption = paste0("Prior was Beta(",
                          paste(modout$inputs$IFRmodel$paramdf[modout$inputs$IFRmodel$paramdf$name == "sens", c("dsc1", "dsc2")], collapse = ","),
                          ")")) +
    theme(legend.position = "none")
  modchain <- modchain[[1]][["trace"]] + 
    labs(caption = paste0("Prior was Norm+(",
                          paste(modout$inputs$IFRmodel$paramdf[modout$inputs$IFRmodel$paramdf$name == "mod", c("dsc1", "dsc2")], collapse = ","),
                          ")")) +
    theme(legend.position = "none")
  sodchain <- sodchain[[1]][["trace"]] +
    labs(caption = paste0("Prior was Beta(",
                          paste(modout$inputs$IFRmodel$paramdf[modout$inputs$IFRmodel$paramdf$name == "sod", c("dsc1", "dsc2")], collapse = ","),
                          ")")) +
    theme(legend.position = "none")
  seroratechain <- seroratechain[[1]][["trace"]] +
    labs(caption = paste0("Prior was Norm+(",
                          paste(modout$inputs$IFRmodel$paramdf[modout$inputs$IFRmodel$paramdf$name == "sero_con_rate", c("dsc1", "dsc2")], collapse = ","),
                          ")")) +
    theme(legend.position = "none")
  # get legend
  legend_bt <- cowplot::get_legend(maxmachain + theme(legend.position = "bottom",
                                                      legend.title = element_blank()))

  if (modout$inputs$account_seroreversion) {
    # weibull shape
    rev_shapechain <- drjacoby::plot_par(modout$mcmcout, "sero_rev_shape", display = FALSE)
    rev_shapechain <- rev_shapechain[[1]][["trace"]] +
      labs(caption = paste0("Prior was Norm+(",  paste(modout$inputs$IFRmodel$paramdf[modout$inputs$IFRmodel$paramdf$name == "sero_rev_shape", c("dsc1", "dsc2")], collapse = ","),
                            ")")) +
      theme(legend.position = "none")
    # weibull scale
    rev_scalechain <- drjacoby::plot_par(modout$mcmcout, "sero_rev_scale", display = FALSE)
    rev_scalechain <- rev_scalechain[[1]][["trace"]] +
      labs(caption = paste0("Prior was Norm+(",  paste(modout$inputs$IFRmodel$paramdf[modout$inputs$IFRmodel$paramdf$name == "sero_rev_scale", c("dsc1", "dsc2")], collapse = ","),
                            ")")) +
      theme(legend.position = "none")


    topp <- cowplot::plot_grid(spechain, senschain, modchain, sodchain, maxmachain, seroratechain,
                               rev_shapechain, rev_scalechain,
                               ncol = 2, nrow = 4)
  } else {
    # out
    topp <- cowplot::plot_grid(spechain, senschain, modchain, sodchain, maxmachain, seroratechain,
                               ncol = 2, nrow = 3)

  }

  cowplot::plot_grid(title, topp, legend_bt, nrow = 3, rel_heights = c(0.1, 1, 0.1))

}

#...................... 
# read in data
#......................
modout <- readRDS(params$path)

# get sim details
simnm <- sub(".RDS", "", basename(params$path))
simnm <- sub("_SeroRev", "", simnm)
simnm <- sub("_NoSeroRev", "", simnm)

if (modout$inputs$account_seroreversion) {
  sim_param_map <- readRDS("data/param_map/SimCurves_serorev/simfit_param_map.RDS") %>% 
    dplyr::filter(sim == simnm)
} else {
  sim_param_map <- readRDS("data/param_map/SimCurves_noserorev/simfit_param_map.RDS") %>%
    dplyr::filter(sim == simnm)
}

truecurve <- sim_param_map %>% 
  dplyr::pull("curve")
# pull out of list
truecurve <- tibble::tibble(time = 1:length(truecurve[[1]]),
                            infxns = truecurve[[1]])

trueseroprev <- sim_param_map %>% 
  dplyr::pull("simdat")
trueseroprev <- trueseroprev[[1]]$StrataAgg_Seroprev


# update title information
type <- switch(sim_param_map$nm,
               "expgrowth" = {"Exp. Growth"},
               "intervene" = {"Outbreak Control"},
               "secondwave" = {"Second Wave"})
sens <- sim_param_map$sens
spec <- sim_param_map$spec

r paste0("Simulation run: ", type, "Sens: ", sens, "; Spec: ", spec)

Fit Diagnostics

Log Likelihood/Prior Plot & Metropolis Coupling Acceptance Plot

#......................
# like and mc accept
#......................
likeplot <- modout$mcmcout$output %>%
  dplyr::filter(rung == "rung50") %>%
  dplyr::filter(stage == "sampling") %>%
  dplyr::select(c("chain", "iteration", "loglikelihood", "logprior")) %>%
  tidyr::gather(., key = "like", value = "val", 3:4) %>%
  ggplot() +
  geom_line(aes(x = iteration, y = val, color = chain), size = 0.25, alpha = 0.8) +
  scale_color_viridis_d() +
  ylab("") + xlab("Iteration") +
  facet_wrap(.~like, scales = "free_y") +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
        legend.position = "bottom",
        legend.title = element_blank())
mcaccplot <- drjacoby::plot_mc_acceptance(modout$mcmcout)

cowplot::plot_grid(likeplot, mcaccplot, ncol = 1, rel_heights = c(1, 0.4))

Gelman-Rubin Convergence Diagnostic

gdiag <- COVIDCurve::get_gelman_rubin_diagnostic(modout)
tibble::tibble(gelman_mpsrf = gdiag$mpsrf) %>% 
  dplyr::mutate(
    gelman_mpsrf = round(gelman_mpsrf, 2),
    gelman_mpsrf = kableExtra::cell_spec(gelman_mpsrf, "latex", 
                                         color = ifelse(gelman_mpsrf >= 1.1, "red", "green"),
                                         bold = T,
                                         font_size = 16)) %>% 
  knitr::kable(., format = "simple")

\pagebreak

Serological Mixing/Prior Diagnostics

quick_sim_sero_diagnostics(modout, 
                           spec = sim_param_map$spec, 
                           sens = sim_param_map$sens)

\pagebreak

IFRs & Incidence Curve

# truth fatality data from simulation workers
fatalitydata <- tibble::tibble(Strata = c("ma1", "ma2", "ma3", "ma4", "ma5"),
                               IFR = c(1e-3, 1e-3, 0.05, 0.1, 0.2),
                               Rho = 1) %>% 
  dplyr::rename(param = Strata)

#...................... 
# get ifrs 
#......................
ifrs <- COVIDCurve::get_cred_intervals(IFRmodel_inf = modout, 
                                       whichrung = "rung50",
                                       what = "IFRparams", by_chain = FALSE) %>% 
  dplyr::left_join(., fatalitydata, by = "param")

ifrs %>% 
  ggplot() +
  geom_pointrange(aes(x = param, ymin = LCI, ymax = UCI, y = median), 
                  color = "#969696", size = 1.2) +
  geom_hline(aes(yintercept = IFR),
             color = "#3182bd", size = 3, alpha = 0.75, show.legend = F) +
  facet_wrap(~param, scales = "free") + 
  theme_bw() +
  theme(axis.text.x = element_blank(),
        axis.ticks.x = element_blank(),
        legend.position = "right") +
  xlab("") + ylab("Median (95% CIs)") +
  labs(caption = "Grey points delineate model fits with Median and 95% Credible Intervals. Blue horizontal line is true simulated value.")

\pagebreak

#.......................
# infection curve and serology
#.......................
infxncurve <- COVIDCurve::draw_posterior_infxn_cubic_splines(IFRmodel_inf = modout,
                                                             whichrung = "rung50",
                                                             dwnsmpl = 1e2,
                                                             by_chain = TRUE,
                                                             by_strata = FALSE)$curvedata
pA <- infxncurve %>%
  dplyr::select(c("chain", "sim", "time", "totinfxns")) %>%
  ggplot() +
  geom_line(aes(x = time, y = totinfxns, group = sim, color = chain), alpha = 0.8, size = 0.8) +
  geom_line(data = truecurve, aes(x = time, y = infxns), color = "#bdbdbd", linetype = "dashed", size = 1.1) +
  scale_color_viridis_d("Chain") +
  labs(caption = "Viridis is inferred, Grey Line is Simulated Infection Curve") +
  xlab("Time") + ylab("Num. Infxns") +
  xyaxis_plot_theme

#...................... 
# seroprev plot
#......................
seropnts <- COVIDCurve::draw_posterior_sero_curves(IFRmodel_inf = modout, 
                                                   whichrung = "rung50",
                                                   dwnsmpl = 1e2,
                                                   by_chain = F)
serocurvedat <- seropnts %>% 
  dplyr::select(c("sim", "ObsDay", dplyr::starts_with("RG_pd_"),
                  dplyr::starts_with("crude_pd_"))) %>% 
  tidyr::pivot_longer(., cols = -c("sim", "ObsDay"),
                      names_to = "seroprev_strata_lvl", values_to = "seroprev") %>% 
  dplyr::mutate(seroprevlvl = ifelse(stringr::str_detect(seroprev_strata_lvl, "RG_"), "RGCorr", "Crude"),
                param = stringr::str_extract(seroprev_strata_lvl, "ma[0-9]+")) %>% 
  dplyr::select(-c("seroprev_strata_lvl")) %>% 
  tidyr::pivot_wider(., names_from = "seroprevlvl", values_from = "seroprev")

# plot
labcaption <- paste(strwrap("Grey boxes indicate model inputs (serosuvery dates and seroprevalences). N.B. If seroreversion is considered, the true prevalence (green) will also be greater than the observed/inferred prevalences.", width = 180), collapse = "\n") 

pB <- serocurvedat %>% 
  dplyr::left_join(., trueseroprev) %>% 
  ggplot() +
  geom_line(aes(x = ObsDay, y = Crude, color = "Crude"), alpha = 0.5) + # yellow
  geom_line(aes(x = ObsDay, y = RGCorr, color = "RGCorr"), alpha = 0.5) + # blue
  geom_line(aes(x = ObsDay, y = TruePrev, color = "TruePrev")) + # green
  geom_line(aes(x = ObsDay, y = ObsPrev, color = "ObsPrev")) +
  geom_rect(aes(xmin = 135, xmax = 145, ymin = 0, ymax = Inf), 
            alpha = 0.05, color = "#bdbdbd", fill = "#bdbdbd") + 
  geom_rect(aes(xmin = 155, xmax = 165, ymin = 0, ymax = Inf), 
            alpha = 0.05, color = "#bdbdbd", fill = "#bdbdbd") + 
  scale_color_manual("", 
                     values = c("Crude" = "#FFD301", "RGCorr" = "#246BCF",
                                "TruePrev" = "#34A853", "ObsPrev" = "#EA4335"),
                      labels = c("Inferred 'Truth'", "Simulated Obs. Prev.",
                                "Inferred 'Observed' - \n Rogan-Gladen Corrected",
                                "Simulated True Prev.")) + # labels follow alphabetical
  facet_wrap(~Strata) + 
  xyaxis_plot_theme + 
  theme(axis.text.x = element_text(angle = 45, vjust = 0.90, hjust= 1, face = "bold")) +
  labs(caption = labcaption) +
  ylab("Seroprevalence") + xlab("Time")

#...................... 
# full out
#......................
cowplot::plot_grid(pA, pB, ncol = 1)

\pagebreak

Death Posterior Predictive Check

#......................
# get deaths posterior pred check
#......................
# recast deaths
proplist <- split(modout$inputs$IFRmodel$data$prop_deaths, 1:nrow(modout$inputs$IFRmodel$data$prop_deaths))

deathrecast <- lapply(proplist,
                      function(x){
                        tibble::tibble(
                          Strata = x$Strata,
                          Deaths = x$PropDeaths * modout$inputs$IFRmodel$data$obs_deaths$Deaths)}) %>% 
  dplyr::bind_rows(.) %>% 
  dplyr::group_by(Strata) %>% 
  dplyr::mutate(ObsDay = 1:nrow(modout$inputs$IFRmodel$data$obs_deaths)) %>% 
  dplyr::ungroup(.)


#...................... 
# get posterior draws
#......................
postdat <- COVIDCurve::posterior_check_infxns_to_death(IFRmodel_inf = modout,
                                                       whichrung = "rung50",
                                                       dwnsmpl = 1e2,
                                                       by_chain = FALSE)
postdat_long <- postdat %>%
  dplyr::select(c("sim", "time", dplyr::starts_with("deaths"))) %>%
  tidyr::gather(., key = "param", value = "deaths", 3:ncol(.)) %>%
  dplyr::mutate(param = gsub("deaths_", "", param)) %>% 
  dplyr::mutate(strata = forcats::fct_reorder(param, as.numeric(stringr::str_extract(param, "[0-9]+"))))

#......................
# get deaths posterior pred check
#......................
datclean <-  deathrecast %>% 
  dplyr::mutate(strata = forcats::fct_reorder(Strata, as.numeric(stringr::str_extract(Strata, "[0-9]+")))) 


#...................... 
# make plot
#......................
ppc_plot <- ggplot() +
  geom_line(data = postdat_long, aes(x= time, y = deaths, group = sim), 
            size = 1.2, color = "#bdbdbd") +
  geom_line(data = datclean, aes(x=ObsDay, y = Deaths), 
            color = "#3182bd") +
  facet_wrap(.~strata) +
  theme_bw() +
  xlab("Time") + ylab("Deaths") + 
  theme(plot.title = element_text(hjust = 0.5, vjust = 0.5),
        plot.subtitle = element_text(hjust = 0.5, vjust = 0.5)) 

plot(ppc_plot)

\pagebreak

Additional Parameters

# knots ifrs 
knots <- COVIDCurve::get_cred_intervals(IFRmodel_inf = modout, 
                                        whichrung = "rung50",
                                        what = "Knotparams", by_chain = FALSE)

ypos <- COVIDCurve::get_cred_intervals(IFRmodel_inf = modout, 
                                       whichrung = "rung50",
                                       what = "Infxnparams", by_chain = FALSE)

nes <-  COVIDCurve::get_cred_intervals(IFRmodel_inf = modout, 
                                       whichrung = "rung50",
                                       what = "Noiseparams", by_chain = FALSE)

addparams <- dplyr::bind_rows(knots, ypos, nes)

addparams %>% 
  dplyr::mutate_if(is.numeric, round, 4) %>% 
  knitr::kable(., format = "simple")

mrc-ide/reestimate_covidIFR_analysis documentation built on Nov. 19, 2021, 12:07 p.m.