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/extra_plotting_functions.R")
source("R/my_themes.R")

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

# pieces for below
studyidchar <- title <- toupper(stringr::str_split(basename(params$path), "_age|_carehomes", simplify = T)[[1]])

# update title
if (grepl("_SeroRev", basename(params$path))) { # catch serorevsion
  title <- paste(title, "with the Potential for Seroreversion Included")
  studylvl <- "SeroRev"
} else if (grepl("_carehomes", basename(params$path))) { # catch care homes
  title <- paste(title, "with Care Homes Excluded")
  studylvl <- "CareHome"
} else {
  title <- paste(title, "with Standard Input")
  studylvl <- "standard"
}

r paste(title)

Fit Diagnostics

Log Likelihood & Prior Plot

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())

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

Metropolis Coupling Acceptance Plot

mcaccplot <- drjacoby::plot_mc_acceptance(modout$mcmcout)
mcacclogplot1 <- drjacoby::plot_rung_loglike(modout$mcmcout, x_axis_type = 2, y_axis_type = 2)
mcacclogplot2 <- drjacoby::plot_rung_loglike(modout$mcmcout, x_axis_type = 2, y_axis_type = 3)

cowplot::plot_grid(mcaccplot, mcacclogplot1, mcacclogplot2, nrow = 3)

\pagebreak

Serological Diagnostics

quick_sero_diagnostics(modout)

\pagebreak

IFRs & Incidence Curve

# common name for knit and joins
dictkey <- modout$inputs$IFRmodel$IFRdictkey %>% 
  dplyr::rename(param = Strata,
                strata = ageband)
#...................... 
# get ifrs 
#......................
ifrs <- COVIDCurve::get_cred_intervals(IFRmodel_inf = modout, 
                                       whichrung = "rung50",
                                       what = "IFRparams", by_chain = FALSE)
ifrs <- dplyr::left_join(dictkey, ifrs) %>% 
  dplyr::mutate(strata = forcats::fct_reorder(strata, as.numeric(stringr::str_extract(param, "[0-9]+"))))

#...................... 
# get crude ifrs 
#......................
dscdat <- readRDS("results/descriptive_results/descriptive_results_datamap.RDS") 
crudedat <- dscdat %>% 
  dplyr::filter(study_id == studyidchar & breakdown == "ageband")  %>%
  dplyr::select(c("study_id", "plotdat")) %>%
  tidyr::unnest(cols = "plotdat") %>% 
  dplyr::rename(strata = ageband)

crudeplotdat <- crudedat %>%
  dplyr::filter(seromidpt == obsday) %>%
  dplyr::mutate(infxns = popn * seroprev,
                crudeIFR =  cumdeaths/infxns,
                crudeIFR = ifelse(crudeIFR > 1, NA, crudeIFR)) %>%  # protect against small denom
  dplyr::select(c("strata", "crudeIFR")) %>%
  dplyr::group_by(strata) %>%
  dplyr::mutate(serodayfct = factor(1:dplyr::n()),
                strata = factor(strata, levels = dictkey$strata)) %>% # add factors in for plot order
  dplyr::ungroup(.)


#...................... 
# get incidence curve
#......................
infxncurve <- COVIDCurve::draw_posterior_infxn_cubic_splines(IFRmodel_inf = modout,
                                                             whichrung = "rung50",
                                                             dwnsmpl = 1e2,
                                                             by_chain = FALSE, 
                                                             by_strata = TRUE)
#...................... 
# make ifr and incidence plots
#......................
plot1 <- ggplot() +
  geom_pointrange(data = ifrs, aes(x = strata, ymin = LCI, ymax = UCI, y = median), 
                  color = "#969696", size = 1.2) +
  geom_point(data = crudeplotdat, aes(x = strata, y = crudeIFR, shape = serodayfct), 
             color = "#000000", size = 3, alpha = 0.75, show.legend = F) +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 45, vjust = 0.90, hjust= 1, face = "bold"),
        legend.position = "right") +
  xlab("") + ylab("Median (95% CIs)") +
  labs(caption = "Grey and black points represent model fits and crude IFR estimates, respectively")
plot2 <- infxncurve$plotObj 
# main plot
main_plotObj <- cowplot::plot_grid(plot1, plot2, ncol = 1, nrow = 2)

# quick out
jpgsnapshot(paste0("figures/", studyidchar, "-", studylvl, "_IFRplot.jpg"), 
            plot = main_plotObj)

plot(main_plotObj)

Strata IFR

\pagebreak

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

Overall IFR

# overall
COVIDCurve::get_overall_IFR_cred_intervals(IFRmodel_inf = modout, 
                                         whichrung = "rung50", 
                                         whichstandard = "arpop",
                                         by_chain = FALSE) %>% 
  dplyr::mutate_if(is.numeric, round, 4) %>% 
  knitr::kable(., format = "simple")

# by chain
COVIDCurve::get_overall_IFR_cred_intervals(IFRmodel_inf = modout, 
                                         whichrung = "rung50", 
                                          whichstandard = "arpop",
                                         by_chain = TRUE) %>% 
  dplyr::mutate_if(is.numeric, round, 4) %>% 
  knitr::kable(., format = "simple")

#...................... 
# global ifr iteration plot
#......................
# make weighted demog
demogwi <- modout$inputs$IFRmodel$demog %>%
  dplyr::mutate(wi = popN/sum(popN))

modout$mcmcout$output %>%
  dplyr::filter(stage == "sampling" & rung == "rung50") %>%
  tidyr::pivot_longer(., cols = modout$inputs$IFRmodel$IFRparams, # if chain isn't included in vector, grepl won't do anything
                      names_to = "Strata", values_to = "est") %>%
  dplyr::left_join(., demogwi, by = "Strata") %>%
  dplyr::mutate(est = est*wi) %>%
  dplyr::group_by(iteration, chain, rung) %>% # need to make sure we capture only the strata levels
  dplyr::summarise(est = sum(est)) %>% 
  ggplot() +
  geom_line(aes(x = iteration, y = est, color = chain), size = 0.25, alpha = 0.8) +
  scale_color_viridis_d() +
  ylab("Global IFR") + xlab("Iteration") + 
  theme_bw() + 
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
        legend.position = "none")

\pagebreak

Seroprevalence Observed vs. Inferred True Prev

# population denom 
demog <- modout$inputs$IFRmodel$demog %>% 
  dplyr::rename(param = Strata)
demog <- demog %>% 
  dplyr::left_join(., dictkey)

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") %>% 
  dplyr::left_join(., dictkey) %>% 
  dplyr::left_join(., demog) %>% 
  dplyr::mutate(strata = forcats::fct_reorder(strata, as.numeric(stringr::str_extract(param, "[0-9]+"))))

# crude data
SeroPrevObs <- crudedat %>% 
  dplyr::select(c("obsdaymin", "obsdaymax", "strata", "seroprev")) %>% 
  dplyr::filter(!duplicated(.)) %>% 
  dplyr::left_join(., demog) %>% 
  dplyr::mutate(obsmidday = (obsdaymin + obsdaymax)/2,
                strata = forcats::fct_reorder(strata, as.numeric(stringr::str_extract(param, "[0-9]+"))))

serocurvedatPlotObj <- serocurvedat %>% 
  ggplot() +
  geom_line(aes(x = ObsDay, y = Crude, group = sim, color = "Crude"), 
            alpha = 0.5) +
    geom_line(aes(x = ObsDay, y = RGCorr,  group = sim, color = "RGCorr"),
            alpha = 0.5) +
  geom_rect(data = SeroPrevObs, aes(xmin = obsdaymin, xmax = obsdaymax, ymin = -Inf, ymax = Inf),
            color = "#d9d9d9", fill = "#d9d9d9", alpha = 0.8) +
  geom_point(data = SeroPrevObs, aes(x = obsmidday, y = seroprev, group = strata),
             color = "#000000", size = 1.2) +
  facet_wrap(.~strata) +
  scale_color_manual("Seroprev. \n Adjustment", 
                     values = c("Crude" = "#FFD301", "RGCorr" = "#246BCF"),
                     labels = c("Inferred 'Truth'", "Inferred 'Observed' - \n Rogan-Gladen Corrected")) +
  labs(caption = "Grey box is observed seroprevalence across study period") +
  xyaxis_plot_theme + 
  theme(axis.text.x = element_text(angle = 45, vjust = 0.90, hjust= 1, face = "bold"))

jpgsnapshot(paste0("figures/", studyidchar, "-", studylvl, "_SeroPrevPDplot.jpg"), 
            plot = serocurvedatPlotObj)

plot(serocurvedatPlotObj)

\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::left_join(., y = dictkey) %>% 
  dplyr::mutate(strata = forcats::fct_reorder(strata, as.numeric(stringr::str_extract(param, "[0-9]+"))))

#......................
# get deaths posterior pred check
#......................
datclean <-  deathrecast %>% 
  dplyr::left_join(., modout$inputs$IFRmodel$IFRdictkey, by = "Strata") %>% 
  dplyr::rename(strata = ageband)
datclean <- datclean %>% 
  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() +
  labs(caption = "Grey Lines are Draws from Posterior, Blue Line is Observed Data (from JHU CSSE COVID-19 Data)") +
  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)) 

jpgsnapshot(paste0("figures/", studyidchar, "-", studylvl, "_PPCplot.jpg"), 
            plot = ppc_plot)

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.